Variance Analysis of Linear SIMO Models with Spatially Correlated Noise

0.7 MB·English

by  Niklas Everitt

VarianceAnalysisofLinearSIMOModelswithSpatially
(cid:63)
CorrelatedNoise
NiklasEveritta,GiulioBottegala,CristianR.Rojasa,H˚akanHjalmarssona
aACCESSLinnaeusCenter,SchoolofElectricalEngineering,KTH-RoyalInstituteofTechnology,Sweden
5
1
0 Abstract
2
Substantialimprovementinaccuracyofidentifiedlineartime-invariantsingle-inputmulti-output(SIMO)dynamicalmodelsis
n
possiblewhenthedisturbancesaffectingtheoutputmeasurementsarespatiallycorrelated.Usinganorthogonalrepresentation
a
J forthemodulescomposingtheSIMOstructure,inthispaperweshowthatthevarianceofaparameterestimateofamodule
is dependent on the model structure of the other modules, and the correlation structure of the disturbances. In addition, we
3
quantifythevariance-errorfortheparameterestimatesforfinitemodelorders,wheretheeffectofnoisecorrelationstructure,
1
model structure and signal spectra are visible. From these results, we derive the noise correlation structure under which the

mentioned model parameterization gives the lowest variance, when one module is identified using less parameters than the
]
Y othermodules.
S
. Keywords: Systemidentification,Asymptoticvariance,LinearSIMOmodels,Least-squares.
s
c
[

1 1 Introduction available information to the identification process, we
v ask the fundamental questions: will, and how much, an
7
added sensor improve the accuracy of an estimate of a
7 Recently,systemidentificationindynamicnetworkshas
certaintargettransferfunctioninthenetwork?Weshall
0 gained popularity, see e.g., Van den Hof et al. (2013);
attempt to give some answers by focusing on a special
3 Dankers et al. (2013b,a, 2014); Ali et al. (2011); Mat-
case of dynamic networks, namely single-input multi-
0 erassietal.(2011);Torresetal.(2014);HaberandVer-
1. haegen (2014); Gunes et al. (2014); Chiuso and Pil- output(SIMO)systems,andinawidercontext,dynamic
networks.
0 lonetto (2012). In this framework, signals are modeled
5 asnodesinagraphandedgesmodeltransferfunctions.
1 To estimate a transfer function in the network, a large SIMO models are interesting in themselves. They find
: numberofmethodshavebeenproposed.Somehavebeen applications in various disciplines, such as signal pro-
v
showntogiveconsistentestimates,providedthatacer- cessing and speech enhancement (Benesty et al., 2005),
i
X tainsubsetofsignalsisincludedintheidentificationpro- (Doclo and Moonen, 2002), communications (Bertaux
r cess(Dankersetal.,2013b,a).Inthesemethods,theuser etal.,1999),(Schmidt,1986),(Trudnowskietal.,1998),
a has the freedom to include additional signals. However, biomedical sciences (McCombie et al., 2005) and struc-
little is known on how these signals should be chosen, turalengineering(Ulusoyetal.,2011).Someoftheseap-
and how large the potential is for variance reduction. plications are concerned with spatio-temporal models,
Fromatheoreticalpointofviewthereareonlyafewre- inthesensethatthemeasuredoutputcanbestrictlyre-
sultsregardingspecificnetworkstructurese.g.,H¨agget latedtothelocationatwhichthesensorisplaced(Sto-
al.(2011);Ramazietal.(2014);Wahlbergetal.(2009). ica et al., 1994), (Viberg et al., 1997), (Viberg and Ot-
Togetabetterunderstandingofthepotentialofadding tersten, 1991). In these cases, it is reasonable to expect
that measurements collected at locations close to each
otherareaffectedbydisturbancesofthesamenature.In
(cid:63)
This work was partially supported by the Swedish Re- otherwords,noiseontheoutputscanbecorrelated;un-
searchCouncilundercontract621-2009-4017,andbytheEu-
derstanding how this noise correlation affects the accu-
ropeanResearchCouncilundertheadvancedgrantLEARN,
racyoftheestimatedmodelisakeyissueindata-driven
contract267381.
modelingofSIMOsystems.
Emailaddresses: [email protected](NiklasEveritt),
[email protected](GiulioBottegal),[email protected](Cristian
R.Rojas),[email protected](H˚akanHjalmarsson). For SIMO systems, our aim in this contribution is to
PreprintsubmittedtoAutomatica 14January2015
quantify the model error induced by stochastic distur- 1.1 Contributionofthispaper
bances and noise in prediction error identification. We
consider the situation where the true system can accu-
Asamotivation,letusfirstintroducethefollowingtwo
rately be described by the model, i.e., the true system
outputexample.Considerthemodel:
lies within the set of models used, and thus the bias
(systematic)erroriszero.Then,themodelerrormainly
y (t)=θ u(t 1)+e (t),
consistsofthevariance-error,whichiscausedbydistur- 1 1,1 − 1
y (t)=θ u(t 2)+e (t),
bances and noise when the model is estimated using a 2 2,2 − 2
finitenumberofinput-outputsamples.Inparticular,we
shall quantify the variance error in terms of the noise where the input u(t) is white noise and ek,k = 1,2 is
covariancematrix,inputspectrumandmodelstructure. measurement noise. We consider two different types of
These quantities are also crucial in answering the ques- measurementnoise(uncorrelatedwiththeinput).Inthe
tionswehaveposedabove,namely,whenandhowmuch, firstcase,thenoiseisperfectlycorrelated,letusforsim-
adding a sensor pays off in terms of accuracy of the es- plicity assume that e1(t) = e2(t). For the second case,
timatedtargettransferfunction. e1(t)ande2(t)areindependent.Itturnsoutthatinthe
first case we can perfectly recover the parameters θ
1,1
andθ ,while,inthesecondcasewedonotimprovethe
2,2
There are expressions for the model error, in terms of accuracy of the estimate of θ by also using the mea-
1,1
the asymptotic (in sample size) (co-) variance of the surementy (t).Thereasonforthisdifferenceisthat,in
2
estimated parameters, for a variety of identification thefirstcase,wecanconstructthenoisefreeequation
methodsformulti-outputsystems(Ljung,1999),(Ljung
and Caines, 1980). Even though these expressions cor-
y (t) y (t)=θ u(t 1) θ u(t 2)
respondtotheCram´er-Raolowerbound,theyaretypi- 1 − 2 1,1 − − 2,2 −
callyratheropaque,inthatitisdifficulttodiscernhow
and we can perfectly recover θ and θ , while in the
themodelaccuracyisinfluencedbytheaforementioned 1,1 2,2
second case neither y (t) nor e (t) contain information
quantities.Thereisawell-knownexpressionforthevari- 2 2
aboute (t).
ance of an estimated frequency response function that 1
lends itself to the kind of analysis we wish to facilitate
Alsothemodelstructureplaysanimportantroleforthe
(Ljung, 1985; Yuan and Ljung, 1984; Zhu, 1989). This
benefit of the second sensor. To this end, we consider
formula is given in its SIMO version by (see e.g., Zhu
a third case, where again e (t) = e (t). This time, the
(2001)) 1 2
modelstructureisslightlydifferent:
n
CovGˆ(ejω)≈ NΦu(ω)−1Φv(ω), (1) y1(t)=θ1,1u(t−1)+e1(t),
y (t)=θ u(t 1)+e (t).
2 2,1 2
−
whereΦ istheinputspectrumandΦ isthespectrum
u v Inthiscase,wecanconstructthenoisefreeequation
of the noise affecting the outputs. Notice that the ex-
pression is valid for large number of samples N and
y (t) y (t)=(θ θ )u(t 1).
large model order n. For finite model order there are 1 2 1,1 2,2
− − −
mainly results for SISO models (Hjalmarsson and Nin-
ness, 2006; Ninness and Hjalmarsson, 2004, 2005a,b) The fundamental difference is that now only the differ-
andrecently,multi-input-single-output(MISO)models. ence (θ1,1 θ2,1) can be recovered exactly, but not the
−
For MISO models, the concept of connectedness (Gev- parametersθ1,1andθ2,1themselves.Theycanbeidenti-
ersetal.,2006)givesconditionsonwhenoneinputcan fiedfromy1(t)andy2(t)separately,aslongasy1 andy2
help reduce the variance of an identified transfer func- aremeasured.AsimilarconsiderationismadeinLjung
tion. These results were refined in M˚artensson (2007). et al. (2011), where SIMO cascade systems are consid-
For white inputs, it was recently shown in Ramazi et ered.
al. (2014) that an increment in the model order of one
transferfunctionleadstoanincrementinthevarianceof This paper will generalize these observations in the fol-
another estimated transfer function only up to a point, lowingcontributions:
after which the variance levels off. It was also quan-
tified how correlation between the inputs may reduce Weprovidenovelexpressionsforthevariance-error
•
the accuracy. The results presented here are similar of an estimated frequency response function of a
in nature to those in Ramazi et al. (2014), while they SIMO linear model in orthonormal basis form, or
regard another special type of multi-variable models, equivalently,infixeddenominatorform(Ninnesset
namelymulti-inputsingle-output(MISO)models.Note al.,1999).Theexpressionrevealshowthenoisecor-
that variance expressions for the special case of SIMO relationstructure,modelordersandinputvariance
cascade structures are found in Wahlberg et al. (2009), affectthevariance-erroroftheestimatedfrequency
Everittetal.(2013). responsefunction.
2
(cid:104) (cid:105)
• Fthoerfarenqoune-nwchyitsepiencptruutmsptehcetrbumen,ewfiteoshfotwhewchoerrreelain- wfohresroemne1o≤rt.h.o.n≤ormnmalabnadsisΓif(uqn)c=tionBs1(q),(.q.).,nBmni.(qO)r-,
tionstructureisfocused. {Bk }k=1
thonormalwithrespecttothescalarproductdefinedfor
• Winphuetnaonndeomnoedouflteh,ie.eo.u,ttphuetrse,laistiiodnesnhtiipfiebdetuwsienegnltehses c1om(cid:82)pπlfex(efiuωn)gcti(oeniωs)fd(ωz.),Lge(tz)us:iCntr→oduCc1e×mtheasve(cid:104)cft,ogr(cid:105)n:o=-
parameters, we derive the noise correlation struc- 2π π ∗
tatio−n
tureunderwhichthementionedmodelparameter-
izationgivesthelowesttotalvariance.    
y (t) e (t)
1 1
   
Thepaperisorganizedasfollows:inSection2wedefine y2(t) e2(t)
the SIMO model structure under study and provide an y(t):= .. , e(t):= .. .
 .   . 
expression for the covariance matrix of the parameter    
estimates. Section 3 contains the main results, namely y (t) e (t)
m m
a novel variance expression for LTI SIMO orthonormal
basisfunctionmodels.TheconnectionwithMISOmod-
The noise sequence e(t) is zero mean and temporally
elsisexploredinSection4.InSection5,themainresults { }
white,butmaybecorrelatedinthespatialdomain:
areappliedtoanon–whiteinputspectrum.Insection6
we derive the correlation structure that gives the mini-
E[e(t)]=0
mumtotalvariance,whenoneblockhaslessparameters
than the other blocks. Numerical experiments illustrat- E(cid:2)e(t)e(s)T(cid:3)=δt sΛ (4)
−
ingtheapplicationofthederivedresultsarepresentedin
Section7.AfinaldiscussionendsthepaperinSection8. for some positive definite matrix covariance matrix Λ,
and where E[] is the expected value operator. We ex-
·
pressΛintermsofitsCholeskyfactorization
2 ProblemStatement
Λ=Λ ΛT , (5)
CH CH
whereΛ Rm m islowertriangular,i.e.,
e (t) CH ×
1 y (t) ∈
1
G1(q) Σ  
γ 0 ... 0
. 11
.  
u(t) Gk.(q) ek(t)Σ yk(t) ΛCH =γ2...1 γ.2.2. ...... 00  (6)
... γm1 γm2 ... γmm
e (t)
m y (t)
Gm(q) Σ m forsome{γij}.AlsonoticethatsinceΛ>0,
Fig.1.BlockschemeofthelinearSIMOsystem. Λ−1 =Λ−CHTΛ−CH1 . (7)
We summarize the assumptions on input, noise and
Weconsiderlineartime-invariantdynamicsystemswith
modelasfollows:
one input and m outputs (see Fig. 1). The model is de-
scribedasfollows:
Assumption1 The input u(t) is zero mean station-
{ }
ary white noise with finite moments of all orders, and
     
y1(t) G1(q) e1(t) varianceσ2 >0.Thenoise e(t) iszeromeanandtem-
      porally white, i.e, (4) holds{with}Λ > 0, a positive def-
y2(t) G2(q) e2(t) (cid:2) (cid:3)
 ... = ... u(t)+ ... , (2) isnoimteemρa>tr0ix..TIhteisdaatsasuismgeednethraattedEin|eo(tp)e|n4+lρoop<, th∞atfiosr,
     
the input u(t) is independent of the noise e(t) . The
ym(t) Gm(q) em(t) trueinput{-outp}utbehaviorofthedatagenera{tings}ystem
can be captured by our model, i.e the true system can be
whereq denotestheforwardshiftoperator,i.e.,qu(t)= described by (2) and (3) for some parameters θo Rni,
i ∈
u(t+1)andtheGi(q)arecausalstablerationaltransfer i = 1,...,m, where n1 ... nm. The orthonormal
functions.TheGi aremodeledas basisfunctions k(q) a≤reassu≤medstable. (cid:50)
{B }
G (q,θ )=Γ (q)θ , θ Rni, i=1,...,m, (3) Remark2 Generalizationofthemodelstructure:
i i i i i
∈
3
Theassumptionthatthemoduleshavethesameor- (1999); S¨oderstr¨om and Stoica (1989)). Λ is assumed
•
thogonal parameterization in (3) is less restrictive known,however,thisassumptionisnotrestrictivesince
than it might appear at first glance, and made for Λcanbeestimatedfromdata.Theestimateofθisgiven
clarity and ease of presentation. A model consist- by
ingofnon-orthonormalbasisfunctionscanbetrans-
formedintothisformatbyalineartransformation, (cid:32) (cid:33) 1
whichcanbecomputedbytheGram-Schmidtproce- θˆ = (cid:88)N ϕ(t)Λ 1ϕT(t) − (cid:88)N ϕ(t)Λ 1y(t). (9)
N − −
dure (Trefethen and Bau, 1997). Notice that fixed
t=1 t=1
denominatormodels,withthesamedenominatorin
alltransferfunctions,alsofitthisframework(Nin-
Inserting (8)in(9)gives
ness et al., 1999). However, it is essential for our
analysisthatthemodulessharepoles.
(cid:32) (cid:33) 1
It is not necessary to restrict the input to be (cid:88)N − (cid:88)N
• white noise. A colored input introduces a weight- θˆN = θ+ ϕ(t)Λ−1ϕT(t) ϕ(t)Λ−1e(t).
ing by its spectrum Φu, which means that the t=1 t=1
basis functions need to be orthogonal with re-
spect to the inner product f,g := fΦ ,g . Under Assumption 1, the noise sequence is zero mean,
If Φu(z) = σ2R(z)R∗(z), w(cid:104)here(cid:105)RΦu(z) is(cid:104)a muoni(cid:105)c hence θˆN is unbiased. It can be noted that this is the
stable minimum phase spectral factor; the trans- sameestimateastheoneobtainedbythepredictioner-
formation Γ˜(q) = σ 1R(q) 1Γ (q) is a procedure ror method and, if the noise is Gaussian, by the max-
i − − i
imum likelihood method (Ljung, 1999). It also follows
that gives a set of orthogonal basis functions in the
thattheasymptoticcovariancematrixoftheparameter
weighted space. If we use this parameterization, all
estimatesisgivenby
the main results of the paper carry over naturally.
However,ingeneral,thenewparameterizationdoes
not contain the same set of models as the original AsCovθˆN =(cid:0)E(cid:2)ϕ(t)Λ−1ϕT(t)(cid:3)(cid:1)−.1 (10)
parameterization.AnotherwayistousetheGram-
Schmidt method, which maintains the same model
Here AsCov θˆ is the asymptotic covariance matrix
set and the main results are still valid. If we would N
of the parameter estimates, in the sense that the
like to keep the original parameterization, we may,
asymptotic covariance matrix of a stochastic sequence
in some cases, proceed as in Section 5 where the
f ,f C1 q isdefinedas1
inputsignalisgeneratedbyanAR-spectrum. { N}∞N=1 N ∈ ×
Notethattheassumptionthatn ... n isnot
1 m
• restrictive as it only represents a≤n orde≤ring of the AsCovfN := limN E[(fN E[fN])∗(fN E[fN])].
modulesinthesystem. N→∞ · − −
In the problem we consider, using Parseval’s formula
2.1 Weightedleast-squaresestimate and(7),theasymptoticcovariancematrix,(10),canbe
writtenas2
Banydinthtreonducimngtθra=nsf(cid:104)eθr1Tf,u.n.c.t,iθomTn(cid:105)mTa∈triRxn, n := (cid:80)mi=1ni AsCovθˆ =(cid:20) 1 (cid:90) πΨ(ejω)Ψ (ejω)dω(cid:21)−1
N ∗
× 2π
π
Γ 0 0  = Ψ,Ψ −−1, (11)
1 (cid:104) (cid:105)
 
Ψ˜(q)=0 ... 0 , where
 
0 0 Γ
m 1
Ψ(q)= σΨ˜(q)Λ−CHT. (12)
wecanwritethemodel(2)asalinearregressionmodel
y(t)=ϕT(t)θ+e(t), (8) Note that Ψ(q) is block upper triangular since Ψ˜(q) is
blockdiagonalandΛ−CHT isuppertriangular.
where
1 Thisdefinitionisslightlynon-standardinthatthesecond
ϕT(t)=Ψ˜(q)Tu(t).
term is usually conjugated. For the standard definition, in
general,allresultshavetobetransposed,however,allresults
An unbiased and consistent estimate of the parameter inthispaperaresymmetric.
vector θ can be obtained from weighted least-squares, 2 Non-singularityof(cid:104)Ψ,Ψ(cid:105)usuallyrequiresparameteriden-
with optimal weighting matrix Λ 1 (see, e.g., Ljung tifiabilityandpersistenceofexcitation(Ljung,1999).
−
4
2.2 Theintroductoryexample correlation of the noise processes, i.e., they are in-
dependentofthevalueofβ.However,notethatthe
crosscorrelationbetweenθˆ andθˆ in(16):
With formal assumptions in place, we now consider the 1,1 2,1
introductory example in greater detail. Consider the
(cid:112)
model Var(θˆ1,1 1 β2θˆ2,1)
− −
(cid:34) (cid:35)T (cid:34) (cid:35)
1 1 1
yy1((tt))==θθ1,1qq−11uu((tt))++θe1(tq), 2u(t)+e (t) ((1143)) = (cid:112)1 β2 σ2Λ (cid:112)1 β2
2 2,1 − 2,2 − 2 − − − −
1
= β2. (17)
whichusesthedelaysq 1 andq 2 asorthonormalbasis σ2
− −
functions. With θ = [θ θ θ ]T; the corresponding
1,1 2,1 2,2
regressionmatrixis Thiscrosscorrelationwillinduceacrosscorrelation
inthetransferfunctionestimatesaswell.
(cid:34) (cid:35) (2) As seen in (16), the variance of θˆ1,2 strongly de-
u(t 1) 0 0
ϕ(t)T = − . pends on β. In particular, when β tends to 0, one
0 u(t 1) u(t 2) is ideally able to estimate θˆ1,2 perfectly. Note that
− − in the limit case β = 0 one has e (t) = e (t), so
1 2
that (2) can be rearranged to obtain the noise-free
Thenoisevectorisgeneratedby
equation
(cid:34) (cid:35) (cid:34) (cid:35)(cid:34) (cid:35)
e1(t) =Lw(t)= (cid:112) 1 0 w1(t) , (15) y1(t)−y2(t)=(θ1,1−θ2,1)u(t−1)+θ1,2u(t−2),
e (t) 1 β2 β w (t)
2 2
− whichshowsthatbothθ andthedifferenceθ
1,2 1,1
−
θ can be estimated perfectly. This can of course
wherew (t)andw (t)areuncorrelatedwhiteprocesses 2,1
1 2 alsobeseenfrom(16),cf.(17).
with unit variance. The parameter β [0, 1] tunes the
∈
correlationbetweene (t)ande (t).Whenβ =0,thetwo
1 2 The example shows that correlated measurements can
processes are perfectly correlated (i.e., identical); con-
be favorable for estimating for estimating certain pa-
versely, when β = 1, they are completely uncorrelated.
(cid:2) (cid:3) rameters, but not necessarily for all. The main focus of
Note that, for every β [0, 1], one has E e (t)2 =
(cid:2) (cid:3) ∈ 1 this article is to generalize these observations to arbi-
E e (t)2 =1.Infact,thecovariancematrixofe(t)be-
2 trarybasisfunctions,numberofsystemsandnumberof
comes
estimatedparameters.Additionally,theresultsareused
(cid:34) (cid:112) (cid:35) to derive the optimal correlation structure of the noise.
1 1 β2 Butfirst,weneedsometechnicalpreliminaries.
Λ=LLT = − .
(cid:112)
1 β2 1
−
2.3 Reviewofthegeometryoftheasymptoticvariance
Then,whencomputing (11)inthisspecificcasegives
Thefollowinglemmaisinstrumentalinderivingthere-
 (cid:112)  sultsthatfollow.
1 1 β2 0
AsCovθˆN = σ12 (cid:112)1−β2 1− 0 . (16) (L2e0m11m))aL3et(JLe:mRmnaIIC.91inqHbjeadlmiffaerrsesnotniaabnledwMi˚athrtreensspseocnt
0 0 β2 stopaθn,naedndbyΨth∈e rLon2w×s→mo;f Ψle×taSnΨd be thre su,brspacneboefaLn1×orm-
Wenotethat: thonormal basis for . Suppo{sBekSt}hka=t1J (θ≤o) Cn q is
Ψ (cid:48) ×
the gradient of J witSh respect to θ and J (θo)∈= Ψ(z )L
AsCov (cid:104)θˆ1,1 θˆ2,1(cid:105)T = σ12Λ, forsomez0 ∈CandL∈Cm×q.Then (cid:48) o
1
AsCovθˆ2,2 = σ2β2. AsCovJ(θˆ )=L (cid:88)r (z ) (z )L. (18)
N ∗ BkS o ∗BkS o
k=1
Theaboveexpressionsrevealstwointerestingfacts:
2.4 Non-estimablepartofthenoise
(1) The (scalar) variances of θˆ and θˆ , namely the
1,1 2,1
estimatesofparametersofthetwomodulesrelated As seen in the example of Section 2.2, strong noise cor-
to the same time lag, are not affected by possible relation may be helpful in the estimation. In fact, the
5
variance error will depend on the non-estimable part of Similartotheabove,fori m,wealsodefine
≤
thenoise,i.e.,thepartthatcannotbelinearlyestimated
fromothernoisesources.Tobemorespecific,definethe (cid:104) (cid:105)T
signal vector e (t) to include the noise sources from ei:m(t):= ei(t) ... em(t) ,
j i
module 1 to mo\dule j, with the one from module i ex-
cluded,i.e., and for j < i we define eˆi:mj(t) as the linear minimum
variance estimate of e (t)|based on the other signals
i:m
e (t):= ej i(t),and
j i \
\  (cid:104) (cid:105)T
 (cid:104) e1(t),...,ej(t) (cid:105)T j <i, Λi:m|j :=Cov[ei:m(t)−eˆi:m|j(t)].
e (t),...,e (t) j =i,
1 i 1
(cid:104)e (t),...,e (t),e −(t),...,e (t)(cid:105)T j >i. As a small example of why this formulation is useful,
1 i 1 i+1 j considerthecovariancematrixbelow,wherethereiscor-
−
relationbetweenanypair(e (t),e (t)):
i j
Now, the linear minimum variance estimate of e given
i
e (t),isgivenby     
j i 1 0.6 0.9 1 0 0 1 0.6 0.9
\
    
eˆ (t):=(cid:37)Te (t). (19) Λ=0.6 1 0.54=0.6 0.8 0 0 0.8 0 .
ij ij j i     
| \
0.9 0.54 1 0.9 0 0.44 0 0 0.44
Introducethenotation
(cid:124) (cid:123)(cid:122) (cid:125)
ΛCH
λ :=Var[e (t) eˆ (t)], (20)
i|j i − i|j FromtheCholeskyfactorizationaboveweseethat,since
γ is zero, Lemma 5 gives that e (t) is orthogonal to
with the convention that λ := λ . The vector (cid:37) in 32 3
i0 i ij e (t) given e (t), i.e., there is no information about
(19)isgivenby | 2 2 3
\
e (t) in e (t) if we already know e (t). This is not ap-
3 2 1
parentfromΛwhereeveryentryisnon-zero.Ifweknow
(cid:2) (cid:3) 1 (cid:2) (cid:3)
(cid:37)ij = Covej i(t) − E ej i(t)ei(t) . e (t) a considerable part of e (t) and e (t) can be esti-
\ \ 1 2 3
mated.Withoutknowinge (t),λ =λ =λ =1,while
1 1 2 3
Wecall ifweknowe (t),λ =0.64andλ =0.19.
1 21 31
| |
e (t) eˆ (t)
i − i|j 3 Mainresults
thenon-estimablepartofe (t)givene (t).
i j i
\ In this section, we present novel expressions for the
variance-errorofanestimatedfrequencyresponsefunc-
Definition4 When eˆ (t) does not depend on e (t),
ij k tion. The expression reveals how the noise correlation
where 1 k j, k = i|, we say that e (t) is orthogonal
≤ ≤ (cid:54) i structure, model orders and input variance affect the
toe (t)conditionallytoe (t).
k j i variance-erroroftheestimatedfrequencyresponsefunc-
\
tion. We will analyze the effect of the correlation struc-
The variance of the non-estimable part of the noise is
ture of the noise on the transfer function estimates. To
closely related to the Cholesky factor of the covariance
thisend,collectallmtransferfunctionsinto
matrixΛ.Wehavethefollowinglemma.
(cid:104) (cid:105)
Lemma5 Let e(t) Rm have zero mean and covari- G:= G1 G2 ... Gm .
∈
ance matrix Λ > 0. Let Λ be the lower triangular
CH
Cholesky factor of Λ, i.e., Λ satisfies (5), with γ For convenience, we will simplify notation according to
CH ik
asitsentriesasdefinedby (6).Thenforj <i, { } thefollowingdefinition:
i Definition6 The asymptotic covariance of Gˆ(ejω0) :=
(cid:88)
λij = γi2k. G(ejω0,θˆN)forthefixedfrequencyω0 isdenotedby
|
k=j+1
AsCovGˆ.
Furthermore,γ =0isequivalenttothat e (t)isorthog-
ij i
onaltoej(t)conditionallytoej i(t). In particular, the variance of Gˆi(ejω0) := Gi(ejω0,θˆiN)
\ forthefixedfrequencyω isdenotedby
0
PROOF. SeeAppendixA. (cid:4) AsVarGˆ .
i
6
Defineχ astheindexofthefirstsystemthatcontainsthe
k
boafssiyssftuenmcstitohnatBdko(enjωo0t)c.oNnotatiicnetthheatbaχskis−fu1niscttihoen.number AsCovθˆ˜1|B1(eiω)|2
.
.
Theorem7 Let Assumption 1 hold. Suppose that the .
parameters θi Rni, i = 1,...,m, are estimated using
∈
weightedleast-squares(9).Lettheentriesofθbearranged AsCovθˆ˜ (eiω)2
asfollows: k|Bk |
θ¯=[θ ... θ θ ... θ ... θ ... ..
1,1 m,1 1,2 m,2 1,n1 .
... θ θ ... θ ... θ ]T.(21)
m,n1 2,n1+1 m,n1+1 m,nm
AsCovθˆ˜ (eiω)2
andthecorrespondingweightedleast-squaresestimatebe nm|Bnm |
denotedbyθˆ¯.Then,thecovarianceofθˆ¯is
AsCovθˆ¯= 1 diag(Λ ,Λ ,...,
σ2 1:m χ2:m|χ2−1 Fig. 2. A graphical representation of AsCovGˆ where each
...,Λ ). (22) term of the sum in (25) is represented by a layer. A basis
χnm:m|χnm−1 function only affects the covariance between modules that
Inparticular,thecovarianceoftheparametersrelatedto alsocontainthatbasisfunction.Thus,thefirstbasisfunction
basisfunctionnumberk isgivenby affects the complete covariance matrix while the last basis
functionn onlyaffectsmodulesχ ,...,m.
m nm
AsCovθˆ¯ = 1 Λ , (23) Remark9 The orthogonal basis functions correspond
k σ2 χk:m|χk−1 to a decomposition of the output signals into orthogonal
components and the problem in a sense becomes decou-
where
pled.Asanexample,considerthesystemdescribedby
θˆ¯k =(cid:104)θˆχk,k ... θˆm,k(cid:105)T , y1(t)=θ1,1 1(q)u(t)+e1(t),
B
y (t)=θ (q)u(t)+e (t),
2 2,1 1 2
andwhere,forχ i m, B
k ≤ ≤ y3(t)=θ3,1 1(q)u(t)+θ3,2 2(q)u(t)+e3(t), (26)
B B
λ
AsVarθˆi,k = i|σχk2−1. (24) pSaurpapmoseetetrh,a(t2w4e)abreecoimnteesrestedinestimatingθ3,2.Forthis
Italsoholdsthat
λ
AsCovGˆ =(cid:88)nm (cid:34)0χ0k−1AsCo0vθˆ¯ (cid:35)|Bk(ejω0)|2, (25) To understand thAesVmaerchθˆ3a,n2i=smsσ32b|2ehind this express(i2o7n),
k=1 k
let u (t) = (q)u(t), and u (t) = (q)u(t) so that
1 1 2 2
B B
the system can be visualized as in Figure 3, i.e., we can
1wherχe AsC1ovmθˆ¯aktriisxgwivitehnablyl e(2n3tr)ieasndeq0uχakl−t1oizserao.χkFo−r consideru1 andu2 asseparateinputs.
k
× −
zχekro=m1a,tr0iχceks−1ofisdiamneenmsipotnysmcoamtrpiax.tibInle(t2o5t)h,e0ddiaegnoonteasl Ftioirnstawboeuotbθserv,eatnhdatthietitseromnlyθy3uthactocnotrnitbauitnisngintfoorymais-
3,2 3,1 1 3
blocks. a nuisance from the perspective of estimating θ . This
3,2
term vanishes when u = 0 and we will not be able to
1
achievebetteraccuracythantheoptimalestimateofθ
PROOF. SeeAppendixB. (cid:4) for this idealized case. So let us study this setting fir3s,t2.
Straightforward application of the least-squares method,
Remark8 Thecovarianceofθˆ¯ ,whichcontainthepa- using u and y , gives an estimate of θ with variance
k 2 3 3,2
rameters related to basis function k, is determined by λ/σ2, which is larger than (27) when e depends on e
3 1
which other models share the basis function ; cf. (17) and e . However, in this idealized case, y = e and
Bk 2 1 1
of the introductory example. The asymptotic covariance y = e , and these signals can thus be used to estimate
2 2
ofGˆcanbeunderstoodasasumofthecontributionsfrom e3. This estimate can then be subtracted from y3 before
eachofthen basisfunctions.Thecovariancecontribu- theleast-squaresmethodisapplied.Theremainingnoise
m
tion from a basis function k is weighted by k(ejω0)2 iny3 willhavevarianceλ32,ife3 isoptimallyestimated
andonlyaffectsthecovarianBcebetweensystem|Bsthatcon|- (see(19)–(20)),andhence|theleast-squaresestimatewill
tainthatbasisfunction,asvisualizedinFigure2. nowhavevarianceλ /σ2,i.e.,thesameas (27).
32
|
7
Tounderstandwhyitispossibletoachievethesameac-
curacy as this idealized case when u is non-zero, we e1(t)
1
need to observe that our new inputs u (t) and u (t) are y (t)
1 2 1
orthogonal(uncorrelated)3.Returningtothecasewhen θ1,1 Σ
onlytheoutputy isusedforestimatingθ ,thisimplies
3 3,2
thatwepaynopriceforincludingthetermθ3,1u1 inour e2(t)
model,andthenestimatingθ andθ jointly,i.e.,the u (t) y (t)
3,1 3,2 1 2
variance of θˆ will still be λ/σ24. The question now is θ2,1 Σ
3,2
ifwecanusey andy asbeforetoestimatee .Perhaps
1 2 3 e (t) eˆ (t)
surprisingly, we can use the same estimate as when u1 3 − 3|2
wnoaws,zeinroa.dTdhiteiorneatdoetrhemparyevoibojuecstotphtaimt tahliessetismtimataetoefweil,l u1(t) θ˜3,1 Σ y˜3(t)
3
contain a term which is a multiple of u . However, due
1
to the orthogonality between u and u , this term will
1 2
only affect the estimate of θ (which we anyway were u (t)
3,1 2
θ
not interested in, in this example), and the accuracy of 3,2
the estimate of θ will be λ /σ2, i.e. (27). Figure 4
3,2 32
illustrates the setting with y˜ d|enoting y subtracted by
3 3 Fig. 4. The SIMO system of Remark 9, described by (26),
theoptimalestimateofe .Inthefigure,thenewparam-
3 butwithy˜ denotingy subtractedbytheoptimalestimate
eter θ˜3,1 reflects that the relation between u1 and y˜3 is ofe3 andθ˜33,1 reflectsth3attherelationbetweenu1 andy˜3 is
differentfromθ3,1asdiscussedabove.Akeyinsightfrom differentfromθ3,1.
thisdiscussionisthatfortheestimateofaparameterin
the path from input i to output j, it is only outputs that Corollary10 Let the same assumptions as in Theo-
are not affected by input i that can be used to estimate rem7hold.Then,foranyfrequencyω ,itholdsthat
0
the noise in output j; when this particular parameter is
estimated,usingoutputsinfluencedbyinputiwillintro- (cid:88)ni
duce a bias, since the noise estimate will then contain a AsVarGˆ = (ejω0)2AsVarθˆ , (28)
i k i,k
term that is not orthogonal to this input. In (24), this |B |
k=1
manifestsitselfinthatthenumeratorisλ ,onlythe
χ 1firstsystemsdonotcontainu . i|χk−1 where
k i
−
λ
AsVarθˆi,k = i|σχk2−1, (29)
e (t)
1 y (t)
1
θ1,1 Σ andλ isdefinedin (20).
ij
|
e (t)
u (t) 2 y (t)
1 θ2,1 Σ 2 PROOF. Follows from Theorem 7, since (28) is a di-
agonalelementof (25). (cid:4)
e (t)
3 y (t)
3
θ3,1 Σ FromCorollary10,wecantellwhenincreasingthemodel
orderofG willincreasetheasymptoticvarianceofGˆ .
j i
u2(t) Corollary11 Under the same conditions as in Theo-
θ
3,2 rem7,ifweincreasethenumberofestimatedparameters
ofG fromn ton +1,theasymptoticvarianceofG will
j j j i
Fig.3.TheSIMOsystemofRemark9,describedby (26). increase,ifandonlyifallthefollowingconditionshold:
(1) n <n ,
j i
(2) e (t) is not orthogonal to e (t) conditioned on
We now turn our attention to the variance of the indi- i j
e (t),
vidualtransferfunctionestimates. j i
(3) |B\nj+1(ejω0)|2 (cid:54)=0.
3 Thissinceu(t)iswhiteandB andB areorthonormal.
1 2
4 Withu andu correlated,thevariancewillbehigher,see
Section41forafu2rtherdiscussionofthistopic. PROOF. SeeAppendixC. (cid:4)
8
Remark12 Corollary 11 explicitly tells when an in- Ifthemodelordersaredistinct,thecorrespondinggraph
crease in the model order of Gj from nj to nj +1 will is given in Figure 5, where one can see that AsVarGˆ4
increase the variance of Gi. If nj ≥ ni, there will be no dependsony2(andony4ofcourse),butdependsneither
increase in the variance of Gi, no matter how many ad- on y nor y since γ =γ =0, AsVarGˆ depends on
3 1 43 41 3
ditionalparametersweintroducetothemodelGj,which y ,butnotony sinceγ =0andAsVarGˆ dependson
was also seen the introductory example in Section 2.2. 1 2 32 2
y ,whilethevarianceofGˆ isgivenby(30).Ifn =n ,
Naturally, if e (t) is orthogonal to e (t) conditioned on 1 1 2 4
i j thefirstconditionofCorollary11isnotsatisfiedandwe
e (t), eˆ (t) does not depend on e (t) and there is no
inj\cireaseiin|jvarianceofGˆ ,cf.Remajrk9. have to cut the edge between Gˆ4 and Gˆ2. Similarly, if
i n = n , we have to cut the edge between Gˆ and Gˆ ,
1 2 2 1
and if additionally n = n = n , we have to cut the
1 2 3
3.1 AgraphicalrepresentationofCorollary11 edgebetweenGˆ andGˆ .
3 1
FollowingthenotationinBayesianNetworks(Koskiand
Noble, 2012), Conditions 1) and 2) in Corollary 11 can Gˆ Gˆ Gˆ Gˆ
1 2 3 4
be interpreted graphically. Each module is represented
by a vertex in a weighted directed acyclic graph . Let
G
theverticesbeorderedbymodelorder,i.e.,letthefirst
vertexcorrespondtoGˆ .UnderAssumption1,withthe Fig. 5. Graphical representation of Conditions 1) and 2) in
1
additionalassumptionthatmoduleiisthefirstmodule Corollary11fortheCholeskyfactorgivenin(31).
with order n , let there be an edge, denoted by j i,
i
→
fromvertexjtoi,j <i,ife (t)isnotorthogonaltoe (t)
i j
conditioned on e (t). Notice that this is equivalent to
j i 4 ConnectionbetweenMISOandSIMO
γ = 0. Let the\weight of the edge be γ and define
ij ij
(cid:54)
the parents of vertex i to be all nodes with a link to
vertexi,i.e.,pa (i):= j :j i .Then,(29),together There is a strong connection between the results pre-
withLemma5,sGhowsth{atonl→you}tputscorrespondingto sentedhereandthoseregardingMISOsystemspresented
parentsofnodeiaffecttheasymptoticvariance.Indeed, in Ramazi et al. (2014). We briefly restate the problem
avertexwithoutparentshasvariance formulationandsomeresultsfromRamazietal.(2014)
toshowtheconnection.TheMISOdatageneratingsys-
AsVarGˆ = λi (cid:88)ni (ejω0)2, (30) tmemspiastiinallsyomcoerrseelnasteedthinepduutaslaonfdthoeneSIoMutOpucta,sae.MWISitOh
i σ2 |Bk |
k=1 systemisdescribedby
whichcorrespondsto(28)with  
u (q)
1
 
λ =...=λ =λ . (cid:104) (cid:105)u2(q)
i|0 i|i−1 i y(t)= G1(q) G2(q) ... Gm(q)  .. +e(t).
 . 
Thus,AsVarGˆ isindependentofthemodelorderofthe  
i u (q)
othermodules. m
Theinputsequence u(t) iszeromeanandtemporally
As an example, consider four systems with the lower { }
white,butmaybecorrelatedinthespatialdomain,
Choleskyfactorofthecovarianceofe(t)givenby:
  E[u(t)]=0
1 0.4 0.2 0 E(cid:2)u(t)u(s)T(cid:3)=δ Σ ,
  t s u
0.4 1.16 0.08 0.3 −
Λ= 
  forsomepositivedefinitematrixcovariancematrixΣ =
0.2 0.08 1.04 0  u
  Σ ΣT ,whereΣ isthelowertriangularCholesky
CH CH CH
0 0.3 0 1.09 factorofΣ.Thenoisee(t)iszeromeanandhasvariance
  T λ. The asymptotic covariance of the estimated parame-
1 0 0 0 1 0 0 0
terscanbeexpressedusing (11)with
  
0.4 1 0 00.4 1 0 0
=   (31) Ψ =ΨMISO :=Ψ˜Σ . (32)
0.2 0 1 00.2 0 1 0 CH
  
0 0.3 0 1 0 0.3 0 1 We make the convention that (cid:80)k2 x = 0 whenever
(cid:124) (cid:123)(cid:122) (cid:125) k=k1 k
k >k .
ΛCH 1 2
9
Theorem13(Theorem4inRamazietal.(2014)) MISO systems. Let the first χ systems contain basis
k
UnderAssumption1,butwithn n ..., n ,for functionk,so
1 2 m
≥ ≥ ≥
anyfrequencyω itholdsthat
0
(cid:88)m λ (cid:88)nj AsVarθˆ¯kMISO =λΣ1−:χ1k
AsVarGˆ = (ejω0)2, (33)
i σ2 |Bk | where Σ denotes the covariance matrix of the first
j=i i|j k=nj+1+1 χ input1s:.χHkence
k
ewshteimreabnlme+pa1rt:=of0uia(nt)dgσivi2|ejniusjt\hi(et)v.ariance of the non- AsCovGˆ =λ(cid:88)n1 (cid:34)Σ1−0:χ1k 0 0 (cid:35)|Bk(ejω0)|2,
Corollary14(Corollary6inRamazietal.(2014)) k=1 m−χk
Under Assumption 1, but with n1 n2 ..., nm. and
≥ ≥ ≥
Suppose that the order of block j is increased from n to
j
nanj+ce1.oTfGhˆeniftahnedreoinslyanifianllctrheaesfeolilnowthinegacsoynmdpittiootnicshvaorldi-: AsCovθˆ¯MISO =λ diag(Σ1−:χ11,Σ1−:χ12,...,Σ1−:χ1nm).(36)
i
Note that, while the correlation between the noise
(1) n <n ,
j i sourcesisbeneficial,thecorrelationintheinputisdetri-
(2) u (t) is not orthogonal to u (t) conditioned on
i j mentalfortheestimationaccuracy.Intuitively,ifweuse
u (t),
j i thesameinputtoparallelsystems,andonlyobservethe
(3) |Bn\j+1(ejω0)|2 (cid:54)=0. sum of the outputs, there is no way to determine the
contribution from the individual systems. On the other
Remark15 The similarities between Corollary 10 and hand, as observed in the example in Section 2.2, if the
Theorem13,andbetweenCorollary11andCorollary14 noise is correlated, we can construct equations with re-
arestriking.Inbothcasesitisthenon-estimablepartof ducednoiseandimprovetheaccuracyofourestimates.
theinputandnoise,respectively,alongwiththeestimated
basis functions that are the key determinants for the re-
This difference may also be understood from the struc-
sulting accuracy. Just as in Corollary 10, Theorem 13
ture of Ψ, which through (11) determines the variance
canbeexpressedwithrespecttothebasisfunctions:
properties of any estimate. Consider a single SISO sys-
tem G as the basic case. For the SIMO structure con-
AsVarGˆ =(cid:88)ni AsVarθˆ (ejω0)2. (34) sidered1in this paper, as noted before, ΨSIMO of (12) is
i i,k|Bk | block upper triangular with m columns (the number of
k=1 outputs), while ΨMISO is block lower triangular with as
However,now many columns as inputs. ΨMISO is block lower triangu-
lar since Ψ˜ is block diagonal and Σ is lower triangu-
CH
AsVarθˆi,k = σ2λ (35) lcaorrriensp(3o2n)d.sAtdodeinxgteanndionugtpΨuStIMyOj twoitthheoSnIeMcoOlusmtrnuc(taunrde
i|χk nj rows):
whereσ2 isdeterminedbythecorrelationstructureofthe
inputs ui|il(t) to the systems Gi(q,θi) that do share basis ΨSIMO =(cid:34)ΨSIMO (cid:63)(cid:35), (37)
function k(q) (i = 1,...,χk). Note that in the SIMO e 0 (cid:63)
B
casewehad
λ wherethezerocomesfromthatΨSIMO alsoisblockup-
AsVarθˆi,k = σi|χ2k per triangular. Since ΨMISO is bleock lower triangular,
addinganinputu totheMISOstructureextendsΨMISO
j
wthheerneoλisie|χksoiusrdceetserem(itn)eadffbeycttihnegcsoyrsrteelamtisonGs(tqru,θct)urtehaotf withnj rows(andonecolumn):
i i i
(cid:34) (cid:35)
do not share basis function k(q) (i = 1,...,χk). Note ΨMISO 0
that (33) found in Ramazi eBt al. (2014) does not draw ΨMISO = , (38)
e
theconnectiontothevarianceoftheparameters.Thisis (cid:63) (cid:63)
madeexplicitinthealternateexpressions (35)and (34).
where (cid:63) denotes the added column and added row re-
Thecorrelationbetweenparametersrelatedtothesame spectively.AdditionofcolumnstoΨ decreasesthevari-
basis functions is not explored in Ramazi et al. (2014). anceofG ,whileadditionofrowsincreasesthevariance.
1
Infact,itispossibletofollowthesamelineofreasoning First, a short motivation of this will be given. Second,
leadingtoTheorem7andarriveatthecounter-partfor wewilldiscusstheimplicationforthevarianceanalysis.
10