A frequency approach to topological identification and graphical modeling

A frequency approac h to top ologica l iden tification and grap hical mo deli ng Giacomo Inno c en ti ∗ June 22, 2018 Abstract This w orks explores and illustrates recent results dev eloped by the author in field of dynamical netw ork analysis [12, 21, 13, 20, 22, 11, 19, 10]. The considered approac h is blind, i.e., n o a priori assumptions on the in terconnected systems are a v ailable. Moreov er, the p erspective is that of a simple “observe r” who can perform no kind of tes t on the netw ork in order to study the related resp onse, that is no action or forcing input aimed to reveal particular resp onses of th e system can b e p erformed. In such a scenario a frequen cy based meth o d of inv estiga tion is developed t o obtain u seful insights on the netw ork. The information thus derived can b e fruitfully exploited to build acyclic graphical models, which can b e seen as extension of Ba yesian Netw orks or Mark o v Chains. Moreov er, it is shown that the top ology of p olytree linear netw orks can b e exactly identi fied v ia the same mathematical tools. In this resp ect, it is w orth observing th at imp ortant real systems, suc h as all t he transp ortation netw orks, fit this class. 1 In tro duction 1.1 Problem bac kground In the la st decade the a nalysis and the study of netw orks has b ecome ubiquitous. In particular, in the study of co mplex systems the comprehensio n o f the internal connections among the elementary units, which define the underlying structure of the pro cess, turns out to play a k ey role to fully understand its principal mechanisms. While net works of dynamical systems have been deeply studied and a nalysed in sy stems theory (see [27, 1 6, 17] and references within for a s ample of recent studies), a great need still exists for metho ds to infer their structure and their top olo gical characteris tics [2 3]. Indeed, in most engineering scenarios the net work is given or it is the very ob jective o f design. On the other hand, dy na mical net works represe nt a powerful to ol to model large sca le systems, dividing the problem into a set of reduced complexity subsystems to b e designed. In this scenar io inferr ing a suitable in ternal topo logy is of fundamen tal imp ortance. In other situations the structure of the netw ork ca n not b e a prior i kno wn, but its agent s ca n be designed to iden tify what their neig hbo urs a r e, t o int eract with them followin g a distr ibuted strateg y and, more generally , to establish any sor t o f consensus [26]. In this case, typical examples are given by senso r netw orks or co op era tiv e con trollers [25]. On the other hand, there are als o interesting s ituations where the agents ar e not directly manipulable, the link structure is actually unknown and it is imp ortant to rec onstruct it along with the under ly ing dynamics. In this p ersp ective, Graphical Mo dels are a mar riage b etw een pro babilit y theory and graph theo ry . They provide a natural to ol for dealing with t w o problems that o ccur throughout applied mathematics and engi- neering – uncer tain ty and co mplexit y – and in particular they ar e playing an increas ingly imp ortant role in ∗ The author is with Dipartiment o di Ingegneria dell’Informazione, Università di Siena, Via Roma 56, 53100, Siena and with Dipartimen to di Sistemi e Informatica, Università di Firenze, Via di S. Marta 3, 50139, Fir enze, e-mail: innocenti@dii. unisi.it 1 the design and analysis of Machine Learning algor ithms and Neural Net w orks. F undamen tal to the idea o f a graphical mo del is the no tion of mo dularity , i.e., a complex systems turn out b y combining simpler parts. Probability theory provides the g lue wher eb y the par ts are combined, ensuring that the system a s a whole is consistent, and providing wa ys to in terface mo dels to data . Man y o f the classica l multiv ariate probabilistic systems studied in fields such as Statistics, Sy s tems Engineering, Information theory , Pattern Recognition and Statistical Mechanics a re sp ecial cases of the genera l graphica l mo del formalism; ex amples include mix- ture mo dels, factor analysis, hidden Marko v models, K alman filters and Is ing models. The graphical mo del framework pr ovides a way to view all of these systems as instances of a common underlying formalism. Relev ant exa mples ab out the imp ortance of the co nnection str ucture in defining the prop erties of a system can b e o bserved in different fields, s uc h as Economics [18], Biology [9, 8], Ecology [33] and Neura l Sciences [2]. W e commonly find this kind of problems in Biolog ical Neural Netw orks [2], Bio chemical Meta bo lic P athw ays [8], Ecology [33] a nd a lso Financial Markets (esp e cially when tra de is happ ening at a high frequency) [24, 32]. In all these ca ses very little is known abo ut the net w ork top ology and the dynamics among the elementary agents. Hence, even inferring so me relation pr op erties betw een tw o differen t no des constitutes a precious information [5]. T o the best knowledge of the author, there a r e only few theore tica l o r methodo lo gical results ab out the reconstruction of a n unknown top olog y related to a set of net work ed dy na mical systems when the agents a r e not manipulable. Recently , this c hallenging pro blem has s ta rted drawing the atten tion of many res earchers, esp ecially in the ph ysics communi ty . Interesting and novel r esults a re given in [3 1], [23] and [5 ], even though they ar e not sy s tematic and do not provide theo retical g uarantees ab out the cor rectness of the identi fied links and the related dynamics. In particular, in [31], [23] and [5] b oth theoretical and numerical pro cedures to detect the active links inside a dynamica l netw o rks are developed for sp ecific cla sses of mo dels. In [23] the attent ion is fo cused on sparse netw orks a ffected by noise and the aim is obtaining a final top olo gy with a reduced num ber of links. No guarantee of an exac t identification is provided. In [31] a more general class of net works is addres sed, but the pr op o sed method to identify their topolo gy is ba sed on the assumption that the input o f every single no de can b e manipulated and that it is p ossible to p erfor m many exp eriments to detect the adjacent links. Unfortunately , such a n hypo thesis is no t feas ible in most real- w orld situations and it is not pra ctical for la rge netw orks. In[5 ] a generalization of Gr a nger causa lit y is introduced to quantify the amount of shared information among times series. Ev en though the ab ov e mentioned works can be settled in the problem of reconstructing the in ternal con- nections of a dynamical netw ork, these approaches a re no t yet sys tematic and they do not provide theore tica l guarantees a bo ut the corr ectness of the identified links and, p os s ibly , o f the co rresp onding dynamics. In this pap er, the res ults developed in the rece nt years by the author in the ab ov e field are pres en ted and summarized. A sho rt introduction o f the r e fer ence scenar io a nd the main a ssumption and idea s ar e rep or ted hereafter. 1.2 Blind frequency approac h In the following a metho d to infer a simplified co nnected gra ph, based on standard to ols fro m identification and filtering theory [14, 1], is presented. The prop osed appr o ach is blind, in the sense that no a priori knowledge ab out the or iginal netw o rk top ology is given. Moreover, every observed pro cess is as so ciated with a distinct ag en t and only their outputs can b e measured. F urthermore, no action can b e per formed to alter the net work b ehaviour in order to test its resp onse to externa l perturba tio ns. Then, we a ssume to mea s ure the s calar outputs { X j } j =1 ,...,N of a net work system compos ed of man y interacting agents without knowing what the internal interconnections are and without any p ossibility to a lter the functioning via a so rt of testing input. The main idea is to estimate each signal b y using the others as the inputs o f a dynamical system to b e iden tified. Indeed, the corresp onding mo dels provide a description o f the dependencies among the pro cess es and of the in ternal dynamics of the whole system, a s well. Moreov er, for each agent, every pro c e s s acting as input in its dynamical mo del can b e asso ciated to a link in the reco nstructed net w ork top o logy . Then, 2 the o ptimal s olution in terms of the minimization of the mode ling err or would b e given by an a pr iori la rge n umber of dynamical mo dels, namely one for every po ssible link. How ev er, this a pproach is not useful b ecause of b oth the a prior i la rge num b er of mo del to b e identified and the lack of prop er distinctions a mong the links. The pro po s ed metho d, instead, aims to recons tr uct a spa rse link structure, that is, to select only a limited n umber of chosen inputs. Suc h a purpo se is motiv ated by the need to provide a final mo del able to provide immediate informatio n ab out the net w ork s tructure by detecting the “most impor tant ” links. This a ppr oach is suggested by the observ ation that in a lar ge netw ork one could ex p ect that most links could be neglected without a ffecting in a significant way the quality of the identification. Consequently , sub optimal stra tegies not only need to b e taken into account, but can a ls o b e mor e informative to detect the int ernal mechanisms of in terdep endence among the signals. Indeed, one may obta in a better unders ta nding o f the net w ork top ology b y forcing a prop er selection o f a reduced num ber o f links. T o this a im, we show ho w a simple sub optimal strateg y can be employ ed to provide a graphical mo del of the netw ork with a r educed co mplexit y , where the links asso ciated to each agent could a lso b e used to single out a selection of inputs for its dyna mical mo del. The main goa l is to obtain a simplified mo del of the netw ork, which can b e used to r eveal some insights about its internal structure. T her efore, g iven a s et of N signals { X j } j =1 ,...,N , the problem of iden tifying, for each of them, the m j signals providing the “b est” estimate is defined. This optimization problem is com binatorially in tractable and it can b e appro ached only for small v alues of m j . In this pap er, how e ver, it is shown that solv ing the problem in the simplest case, where m j = 1 for a n y j , provides information a bo ut the connection structure of the pro cesses. Such a result also suggests to ex plo it the corresp onding topo logy to derive for each j a sub optimal set of m j pro cesses, p ossibly with m j ≥ 1 , in o rder to estimate via mu ltiv ariate linea r filtering the dyna mics of the rela ted X j . W e also examine some genera l prop erties of the netw ork structure, which ca n b e obtained in this spec ia l case, and we inv es tigate how it can b e employed a s a useful to o l to un veil the netw o rk s tr ucture ga ining insigh ts into its top olog y . Such an infor mation not o nly provides a dee p er understanding of the in ternal dynamics, but it also repr esents a useful to ol in the realization of a suitable mo del of the whole system. A class of dyna mical net works that can be co mpletely and exactly disclose d via the prop osed method is also present. Such a class turns out quite co mmon and, for instance, it contains distribution netw orks. The a pplication of the ab ov e approach to clustering problems will also be presented. 1.3 Notation Let R and C respe c tively b e the real and complex num bers sets. Then, for any z ∈ C such tha t z = α + j β , with α, β ∈ R , the complex conjugate o f z is defined as z ∗ = α − j β . Let X ( t ) and Y ( t ) b e time-discrete zero-mean wide-sense stationary pro cess es. W e denote by E [ · ] the mean op erator. It follows that E [ X ( t )] = E [ Y ( t )] = 0 . Since the t wo pr o cesses are wide-sense sta tionary it is p ossible to define the cros s-cov ariance and the cov ariance functions, R X Y ( τ ) := E [ X ( t ) Y ( t + τ )] and R X ( τ ) := R X X ( τ ) , with no a m biguity since there is no dep e ndence on t . Giv en a time-discre te sequence a ( t ) , we also define Z [ a ( · )]( z ) := + ∞ X t = −∞ a ( t ) z − t as the bilatera l Zeta-transfor m of the signal a ( t ) for a n y z ∈ C such that the ab ov e sum exists. Moreover, we denote by Φ X Y ( z ) := Z [ R X Y ( · )]( z ) the cro ss-p ow er sp ectral density b etw een X ( t ) and Y ( t ) and by Φ X ( z ) := Φ X X ( z ) the power spectra l densit y of X ( t ) assuming that the bilateral Zeta- transform con verges in a neig h bo rho o d of the unit circle | z | = 1 . Finally , let φ X ( ω ) = Φ X ( e iω ) and also let {Z [ a ( · )]( z ) } C := + ∞ X t =0 a ( t ) z − t 3 be the causa l truncation of a filter or a signal. 2 Problem set up Pro ofs a nd mathematical details are omitted for the sake of simplicity . See [12, 21, 1 0, 11] for a further reading. 2.1 Mathematical to ols In this section some basic notions of graph theory , which a re functiona l to the following developmen ts ar e recalled as well as the main mathematical to ols. F or an extensive ov erview on graph theory see [7]. First w e provide the standard definition of a graph. Definition 1. An indir e cte d (or u ndir e ct e d) gr aph G is a p air ( V , A ) c omprising a finite set V of vertexes (or no des) to gether with a set A of e dges (or ar cs), which ar e unor der e d s u bsets of two distinct elements of V . A dir e cte d gr aph differs fr om the indir e ct typ e, b e c ause the e dges ar e or der e d c ouple of n o des, the first one b eing denote d as the dep arting no de of the ar c and t he s e c ond the arriving vertex. The or der of the c ouple define the dir e ct ion of the e dge. In a graph it can also b e defined the notion o f “pa th” linking tw o no des as an ordered sequence of contiguous edges having them at the extremities. A path formed by coherently directed e dg es is said a directed pa th. The following definition introduces the notion of “tre e”. Definition 2. A gr aph T = ( V , A ) is a tree if for any p air of distinct no des N i , N j ∈ V ther e is a unique p ath linking them. The mor e standard definition of a tree as an acycl ic and c onne cte d gr aph is equiv alent to Definition 2 a s shown in [7]. How ev er, in the development of our results we will extensively employ the prop erty of uniqueness of a path linking t wo no des. It is also impo rtant to recall the notion of co nnected graph, that consists in the existence of a path b etw een ev ery couple of its no des. Lemma 3. Given a c onne cte d gr aph G = ( V , A G ) , ther e ex ist s a tr e e T = ( V , A T ) such that A T ⊆ A G . The tr e e T is also said to b e a spanning tree for G . The use of a weighting function on the edges of a graph a llows one to define a weigh t for any of its “subgraphs” as the sum of the v alues asso cia ted to each of the corres ponding arcs. Definition 4. Given a c onne cte d gr aph G = ( V , A ) and a weighting function w : A → ℜ , we define Minimum Sp anning T r e e (MST) any sp anning tr e e T of G such that w ( T ) ≤ w ( T ′ ) (1) for any other sp anning tr e e T ′ of G . In this work directed g raphs are a lso a ddressed a nd s o the co ncepts of “ro oted tree” a nd “p olytr ee” are also introduced hereafter. Definition 5. A ro oted tree is a dir e cte d gr aph T = ( V , A ) such that its e dges defines a tr e e and they satisfy the c ondition that only a unique no de N r ∈ V , r eferr e d to as the ro o t , has ex clusively dep arting e dges. A po lytree ext ends the r o ote d tr e e c onc ept in the s en s e that the pr op erty of only dep arting e dges is satisfie d by a finite subset R of the no de set. 4 The presence of ro ot no des denotes a par tial order ing amo ng them in terms of couples of no des linked b y directed paths. F or instance, if tw o vertexes are connected by a directed arc, the departing o ne is s a id a p ar ent of the other, that in turn is deno ted as a child of the fir st no de. In the following the set of a ll the pa rents of the j -th no de will b e addressed as Π( j ) . The notions of a nce s tor and descendant extend and generalize the a bove situation. The previous definitions and lemmas are functional to introd uce a rigorous model to address no is y linea r dynamical systems interconnected to form a tree-like top olo gy . Definition 6. Consider the triple T = ( T , X , H ) wher e • T = ( V , A ) is a p olytr e e c ounting n no des { N j } j =1 ,...,n and k r o ots among them; • X = { X i } i =1 ,...,n is a set of n zer o-me an wide-sense st ationary and er go dic discr ete sto cha stic pr o c esses; • H = { H j i ( z ) } j,i =1 ,...,n is a set of n time-discr ete S IS O tra nsfer functions H j i ( z ) r epr esente d in the domain of the Z -tr ansform and functional ly linking the ou t put of the i -t h no de with the input of t he j -th one; it is also assume d that H j i ( z ) = 0 if and only if ther e is no c onne ctions b et we en the i -th and the j -th system. F or j = 1 , ..., n , define the fol lowing sto chastic pr o c esses   j := X j if N j ∈ R  j := X j − P i ∈ Π( j ) H j i ( z ) X i if N j is child of N i for e ach i ∈ Π( j ) . (2) W e say that T is a Acy clic Line ar Network (ALN) if the fol lowing c ondition is satisfie d E [  j  i ] = 0 when i 6 = j . (3) It is worth underlining that a num b er of r eal world systems can b e represented as ALNs despite the p eculiar link str ucture o f such a mo del. F or instance, a spe c ia l class of ALNs is represented by the tra nsp o rtation systems. A spe c ia l sub class, in pa rticular, consis ts in the p ow er distr ibution netw o rk []. 2.2 ALN p ersp ectiv e Let us as sume that the obs erved system can b e seen as a set o f interconnected pro ce s ses. Then, let us represent such a sy s tem as a set of N agents interacting together according to an unknown structure. Mor eov e r, assume that it is pos sible to obtain scala r mea surements from these a gents in the form of N scala r time series. F urther more, let it b e p oss ible to remove any deterministic co mponent from these time series, in order to obtain N stochastic pro cesse s { X i } i =1 ,...,N which are wide sense stationary a nd with zero mea n [30]. W e in tend to derive g raphical mo del in term of a ALN repr e s en ting the stronges t int erconnections amo ng the pro cesses in o rder to obtain a simplified r epresentation of their unknown connectivity structure. This can be used as s tarting po int for dev eloping more complex to polo gical descriptions, but the links of ea ch a gent, representing its strongest dep e ndences , can a lso be exploited to design a suitable dynamical mo del for it. The appro ach we follow to single out the net work connectivit y is given in terms o f filtering theory [14]. Given tw o s to chastic pr o cesses X i and X j , w e say tha t X j do es no t contain more information that X i if it is p ossible to fully determine X j from the knowledge of X i b y the application of a suitable oper ator. This observ ation pr ovides us with a tool to quanti fy the “re la tedness” of t w o sto chastic pro cess e s in the following wa y: they ar e “related”, if it is possible to provide a “go o d” estimate of one (acco rding to a sp ecific criterion) b y knowing the other, or vicev ersa. A sp ecial situation o ccurs when the op erator mapping X i in to X j is inv e rtible. In this paper, we develop an approa ch intended to b e used in situations where the op erator, providing the be s t estimate of one pro cess given the other, needs to be derived from the data itself (for exa mple b y using a n identification or a matc h filter tech nique). In these ca ses, it is typical to ado pt a 5 parametric a pproach by choo sing a fixed clas s of p ossible mo dels. In this pap er we limit o urselves to the c la ss of linear time-inv ariant op erato rs mapping stationary pro cesses into stationa ry ones and we consider a least square criterion fo r the optimization. Indeed, in such a case the solution of the b est estimate is well-kno wn and it is given b y the Wiener filter. W e examine tw o scenario s: non-causa l a nd causal Wiener filtering. In the non-causal s cenario the cons ider ed trans fo rmations ar e alw ays invertible. This prop erty will lead to the definition of a (pseudo)- metr ic on the set o f stationa ry pro c e sses. In the causal scenario, instead, this pr op erty is los t, but the genera l a pproach to determine an acyclic co nnectivit y structure can s till be applied. A ccording to the ab ov e formulation, we would lik e to develop a suitable graphical mo del des cribing the connectivity structure of the netw ork. In particular, we intend to link each no de to the most “related” ones, where the “ relatedness” cr iterion co mes from SISO Wiener filtering along with the a cyclic top olog y requirement. W e highlight that a similar and more general r esult can b e obtained by mo deling every single pro cess X j as the o utput of a MISO (Multiple Input Sing le Output) linear system generalizing (1 3), X j ( z ) = e j ( z ) + X i ∈ I j W j i ( z ) X i ( z ) , (4) whose set of inputs I ( j ) a re the pro cesses { X τ } τ ∈ I ( j ) providing, a ccording to a chosen cr iterion, the b est trade-off be t w een their num ber and the lar gest reduction of the cost function (5). How ever, such an appro ach relies on several com binatorial sea rches in or der to p oint out the most appropr ia te input set I ( j ) . Therefore, we rather aim to develop a subo ptimal so lution, exploiting the weigh t matrices defined in the previous section, to derive gr aphs, whose top ologies could b e p ossibly used to design the linear MISO mo dels o f each pro ce s s. T o this purp ose, let us represent the pro cess X j ∈ Θ as the j -th node of the set V . Moreover, let us weigh t the edg es connecting the N elemen ts of V acco rding to the matrix W V ∈ R N × N . Mor eov er, co nsider a set A ( W ) ⊂ V × V of (un)directed arcs co nnecting no des of V a nd denote by G = G ( V , A ( W )) the cor resp onding (un)directed g raph. Hence, if the no des, which are closer in the g raph to the j -th one, a re chosen so tha t they represent the pro cess es having the strongest relations with X j , then the input set I ( j ) o f a pos s ible MIMO mo del of X j can b e co ns tructed starting from its neighbours. 2.3 Preliminary results First, we need to introduce a few preliminar y results whic h will be exploited to define a mathematical to ol for the quantitativ e characterization of the connections among the systems of the netw ork. The main idea developed in this section is to define a distance among time series in o r der to establish a useful concept of “closeness” . Sp ecifically , we will a ssume a pro cess X i as the input of a linear sys tem W j i ( Z ) whose output will be used to e s timate X j . The t w o pro cesses will be considered “clos e” if it is p ossible to find a linea r transfer function W j i which provides a “g o o d” estimate of X j . In this r espe c t, o ur approach w ill rely o n standar d least squares techniques, na mely Wiener filter ing. The distance function will b e defined as the v aria nce of the s uitably filtered err or signa l in or der to guar ant ee the prop erties of a metric. In particular, we distinguish t wo ca ses, which we refer to as the “causal” and the “non-caus al” scena r io. W e ar e treating these tw o cases separately , becaus e they will lead to tw o differ e nt kinds of graphs. 2.3.1 Non-causal scenario First, for the sake of completeness, we provide a for m ulation of the Wiener filter that will b e useful in the developmen t o f our metho d. Prop osition 7 (Non-causa l Wiener filter) . Given a fr e quency weighting tr ansfer function Q ( z ) analytic on the unit cir cle | z | = 1 , a zer o-me an wide-sense stationary sc alar pr o c ess Y and a zer o-me an wide-sense stationary ve ctor pr o c ess Z , the Wiener filter mo deling Y by Z is the line ar stable filter ˆ W ( z ) minimizing the 6 varianc e of the err or ε Q := Q ( z )[ Y − W ( z ) Z ] , that is ˆ W ( z ) = arg min W ( z ) E [ ε 2 Q ] . (5) Its expr ession is given by ˆ W ( z ) = Φ Y Z ( z )Φ − 1 Z ( z ) (6) and it do es not dep end up on Q ( z ) . Mor e over, the minimize d c ost is e qual to min E  ε 2 Q  = 1 2 π Z π − π | Q ( e j ω ) | 2  Φ Y ( ω ) − Φ ∗ Y Z ( ω )Φ − 1 Z ( ω )Φ Y Z ( ω )  dω . Given t wo scalar sto chastic pro ces s es X i , X j , let W j i ( z ) b e the Wiener filter mo deling X j from X i . As a consequence of Pro po sition 7 we ha ve that, for any weigh ting transfer function Q j ( z ) E h [ Q j ( z )( X j − W j i ( z ) X i )] 2 i = Z π − π | Q j ( e j ω ) | 2  Φ X j ( ω ) − | Φ X j X i ( ω ) | 2 Φ X i ( ω )  dω . (7) By choosing Q j ( z ) equal to the in v erse of the spectr al factor of Φ X j ( z ) , that is the stable and causally in vertible filter F j ( z ) , such that [29 ]: Φ X j ( z ) = F − 1 j ( z )( F − 1 j ( z )) ∗ , (8) we obtain min E [ ε 2 F j ] = Z π − π  1 − | φ X j X i ( ω ) | 2 φ X i ( ω ) φ X j ( ω )  dω . (9) This sp ecial choice o f Q j ( z ) makes the cost dep end explicitly on the coherence function of the t w o pro cesses [14]: C X i X j ( ω ) := | φ X j X i ( ω ) | 2 φ X i ( ω ) φ X j ( ω ) , (10) which is no n- negative and sy mmetric with resp ect to ω . It is a lso well-known that the Cro s s-Sp e ctral Densit y satisfies the Sc h w artz Inequalit y . Hence, the coherence function is limited b etw een 0 a nd 1 . The ch oice Q j ( z ) = F j ( z ) can be now understo o d as motiv ated by the nec e s sit y to achiev e a dimensionless cost function not dep ending on the p ow er o f the signals . Let us consider, now, a set Θ of discrete zero mean and wide sense stationary sto chastic pro cesses. F or the sake of simplicity , herea fter we will neglec t the s ignal argument t or z , when it can’t be misunder sto o d. The cost obtained by the minimization of the error ε F j using the Wiener filter as b efore a llows us to define on Θ the function d ( X i , X j ) :=  1 2 π Z π − π  1 − C X i X j ( ω )  dω  1 / 2 , ∀ X i , X j ∈ Θ . (11) Prop osition 8. The function d ( · , · ) as define d in (11) is a pseudo-metric. Remark 9. Observe that (11) assumes values b etwe en 0 and 1 . The first c ase happ ens only when C X i X j ( ω ) = 1 for any ω ; the se c ond situation, inste ad, is r elate d to the sc enario wher e C X i X j ( ω ) is always e qual to 0 , i.e., when the cr oss-sp e ctr al density b etwe en X i and X j never assumes values differ ent fr om zer o. F o r a g iv en set Θ counting N pro ce sses we also in tro duce the sy mmetric weigh t matr ix D Θ ∈ R N × N such that D Θ ( j, i ) := d ( X j , X i ) , ∀ X j , X i ∈ Θ . (12) 7 It is worth observing that D Θ ( j, i ) ex pr esses how m uc h g o o d is the linear mo del of X j ( X i ) consisting of the only input pro cess X i ( X j ). Ther efore, the infor mation recorded in the j -th raw (column) allows o ne to so rt the elements of Θ a ccording to their capacity of describing X j in term of the linear SISO (Single Input Single Output) mo del X j = e j + W j i ( z ) X i , i = 1 , . . . , N . (13) In particula r , it holds that D Θ ( j, j ) = 0 ∀ j = 1 , . . . , N . 2.3.2 Causal scenario Given tw o sto chastic pro cess es X i , X j and a transfer function W j i ( z ) , let us consider a g ain the qua dr atic cost ε Q := Q ( z )[ X j − W j i ( z ) X i ] . Analog ously to the pr evious case it is p ossible to derive a causal linear filter minimizing (5). The causal filter pr oviding such a minimization is referred to as the caus al Wiener filter [14]. Prop osition 10 (Caus al Wiener filter) . The Causal Wiener filter m o deling X j by X i is the c ausal stable line ar filter ˆ W C j i ( z ) minimizing the filter e d quantity ε Q j := Q j ( z )[ X j − W j i ( z ) X i ] . It do es not dep end up on Q j ( z ) and its expr ession is given by ˆ W C j i ( z ) = F − 1 j ( z )  F j ( z ) Φ X i X j ( z ) Φ X i ( z )  C , (14) wher e F j ( z ) is the inverse of the sp e ctr al factor of Φ X j ( z ) . Since the weigh ting function Q j ( z ) do es not affect the Wiener filter, but only the energ y of the filtered error, we ca n choose again Q j ( z ) equal to F j ( z ) , the in verse of the sp ectral factor of Φ X j ( z ) . This ch oice simply op erates a normalization of the err or spectrum at every frequency . Th us, similar ly to what we hav e done in the non-causal framework, w e can define the function d C ( X j , X i ) :=  E   F j ( z )( X j − ˆ W C ij ( z ) X i )  2  1 2 (15) to repr esent the modeling erro r o n the pro cess X j due to the a pplication of the causal Wiener filter to X i . T o highlight the prop erties of the pr esent a pproach, it is imp ortant to notice a difference b etw een d C and the function d defined by (11). Indeed, it is straig h tforward to observe that d C is not symmetric and, therefore, that it can’t b e a metric, pr oviding us with a less g e ner al to ol to handle the pro cesses. Ho w ever, d C still repres en ts a mo deling error and in this re s pect it well suits the or iginal problem of quantifying the “relatedness” b et w een tw o pro cesses. As in the previo us case, we intro duce for a given s et Θ counting N pro cesses the w eight matrix D C Θ ∈ R N × N such that D C Θ ( j, i ) := d C ( X j , X i ) , ∀ X j , X i ∈ Θ . (16) It is w orth observing that D C Θ ( j, i ) s till repres e n ts the capacity of X i in describing X j according to the SISO mo del (13), but only via c a usal filters. 3 Graphical mo del s 3.1 ALN form ulation A ccording to the ab ov e int ro duced scenar io, let us depict ea ch sto chastic pro cess X j ( t ) as the sup erp osition of linear dynamic transformations of the other pro cesses ’ o utputs: X j ( z ) = e j ( z ) + N X j =1 ,j 6 = i W j i ( z ) X i ( z ) , (17) 8 where W j i ( z ) is a suitable tr a nsfer function and e j is the mo del err or. In this framework, it can b e considered in teresting to find the set of { W j i ( z ) } i,j =1 ...N ,i 6 = j which allows us to bes t des crib e the time s eries a ccording to the least squares criterio n min X j E  ε 2 j ( t )  , (18) where ε j ( z ) = Q j ( z ) e j ( z ) a nd the Q j ( z ) a re dynamic weigh t functions, i.e. transfer functions. The best description provided by (17) according to (18) c ould b e explo ited in or der to infer dynamic r e lations b et ween the r andom pro ces ses. The problem with such an approa ch is due to the complexity of the final mo del, since it may b e given by an a priori la rge num b er of transfer functions, na mely N ( N − 1 ) . Hence, it is quite natural to develop a subo ptimal strateg y to reduce its complexity . In pa r ticular, this b oils down to the pr o blem of choosing for every j = 1 , . . . , N a set I ( j ) ⊂ { 1 , . . . , j − 1 , j + 1 , . . . , n } , in o rder to satisfy a cer tain tr ade-off betw een the n um ber of elemen ts in eac h I ( j ) and the g lobal cost (18) asso ciated to the mo deling er rors e j ( z ) = X j ( z ) − X τ ∈ I ( j ) W j τ ( z ) X τ ( z ) . (19) The simplified linear mo del (19) repr esents a sub optimal strategy with resp ect to considering all the p ossible contributions pr esent in (17). Nonetheless, if no a prior i assumptions are introduced, the choice of the most suited sets I ( j ) may turn o ut a complex pro ce dure, since se veral combinatorial sea r ches are needed. Therefore, in the following we dev elop a simplified approa c h inspired by graphical mo dels, such as B ay e s ian Net w orks (BN) a nd Markov Random Fields (MRF). These mo dels consist o f interconnected no des representing ra ndom v aria bles , which are emplo yed to provide a sta tistica l description of the system. In particular, in b oth cases the infer e nce pro blem asso cia ted to each no de ca n b e solved by the knowledge of a set o f conditionally r elated random v a riables, which are depicted as the adjac en t no des (MRF) or as a prop er extension of this set (BN). Similarly , we intend to desig n a pro cedure to build fr o m the org inal set of pro ces s es a gra ph, s uch that its top o logy has a direct interpretation in terms of the sets I ( j ) . 3.2 Undirected graph Undirected ar cs a re well suited to repr e s en t mutual relations hips, where one s ub ject affects an other and viceversa. In this p ersp ective, b etw een tw o no des j and i ther e exists only o ne edge and, then, the matrix weigh ting the arcs results symmetric . W e c hoo se W V equal to D Θ as in (12), so that the w eight of each edge is a direct expressio n of the mo deling capability b etw een the connected no des accor ding to the (13) pa radigm. Then, a ssume that, acco rding to the ab ov e W V , the i -th node is the best neighbor of the j -th, while the k -th is the best co nnectio n for the i -th. Observe that the w eight D Θ ( i, k ) is the minimum among all the v alues { D Θ ( i, r ) } , r = 1 , . . . , i − 1 , i + 1 , . . . , N , while D Θ ( i, j ) ma y be ranked in any o ther pos ition. Howev er, even though the j -th no de may not be the most suited o ne to represent by itself t he i -th, the fact, that the i -th is the b est choice to describ e with a s ing le pro cess the j -th, is meaning ful. In particular, we can read this information as the fact that the j -th no de sha res with the i -th a certain amount of infor mation, whic h no other pr o c e ss can mo del. T his idea derive from a common obser v ation in Infor matics: if t w o different signa ls are able to repr esent the 90% of another pro cess, but they do describ e exactly the same 90%, then there is no improvem ent in using them bo th to mo del it. Co n versely , it would b e mo r e useful to find out a signa l describing, for instance, only a 20% of the ob jective pro ce s s, if it can ex pla in so me miss ing information from the fir st chosen s ig nal. In such a scenar io, w e will take b oth pro ces ses X j and X k as inputs for a p oss ible MISO mo del (4) of X i . Observe that the a bove reasoning implicitly a ssumes tha t, from the mo deling point of view, the amount of sha r ed information b etw een tw o pro cesse s decreases as the length o f the (weigh ted) path linking them increases. Figuratively , this prop erty can b e reg arded as though the presence of int ermediate nodes in the path weakens the exchanged information. In this p ersp ective, then, the influence of distant no des r eaches the o b jectiv e pro ce s s as filtered by neighbor s of this latter. Due to this fact, w e wan t to des ign a connected 9 graph, since without any a priori infor mation ab out the signals nature, the be s t pra ctice is to assume that every pro cess is able to affect all the others to a certain extent. Con versely , the choice of enforcing a unique path b etw e e n every c o uple of no des, that is of providing an acyclic graphical mo del, turns out r easonable in order to k eep the complexity low, acco rding to the sub optimal solution persp ective. It is worth observing that a connected graph, such that for every couple of no des it exists a unique path linking them and such that every no de is directly attached to the be s t neighbour accor ding to the w eight matrix, is a Minim um Spanning T r ee [7]. An alternative definition is the connected acyclic gra ph with the least total w eight . Efficient algorithms exist for its computation once the w eight matrix is given. 3.3 Directed graph Directed arcs a re best suited to deal with binar y rela tionships which differ a ccording to the chosen order of the ob jects. Therefor e, since we ar e also co ncerned with ca usal dependencies among the pro cesses , we are going to develop a directed gra ph, who se purp ose is the same as the pr evious one, i.e. providing a re a sonable description of the connectivity structure o f the entire net work. As in the previo us case, moreov er, the re s ulting graphical mo del can b e p ossible exploited in or der to single out a sub optimal solution for the inputs of the MISO mo del (4) of each pr o c e s s. T o this a im, recalling the previous a pproach, it is worth o bserving that D C Θ provides distinguished infor mation for every couple o f no de, according to the choice of the input. How ever, in order to derive a directed g raph, let us c ho o se first W V such that W V ( j, i ) = min  D C Θ ( j, i ) , D C Θ ( i, j )  , ( 20) so that the w eight o f each edge describes the best ca pabilit y of one of the linked no des in des cribing the other acco rding to the (13) paradig m, also taking in to account the differences due to the ca usality pro per t y . Moreov er, to keep tra ck of the b est input in ea ch couple, let us define the matrix H ( W V , D C Θ ) such that H ( W V , D C Θ )( j, i ) =  +1 if W V ( j, i ) = D C Θ ( j, i ) − 1 if W V ( j, i ) = D C Θ ( i, j ) . (21) Similarly to the previous c a se, let us apply the Minimum Spanning T ree a lgorithm to the weigh t matrix W V as defined in (2 0). Then, choo se for eac h edge ( j, i ) the dir ection from the i -th node to the j -th one, if H ( W V , D C Θ )( j, i ) = +1 , and the opp osite other wise. The resulting gr aph is a po lytree with minim um cost [15]. Analogously to the previous ca se, we use the arcs to r epresent a input/output relationships. In particular, observe that the dir ection o f each edge defines whic h no de would act as the input when mo deling the pro cesses via MISO linear representations. Hence, in terms o f the mo del (4), a sub optimal though natura l choice for each I ( j ) turns out the set of the parents no des of the j - th pro c e ss. Observe that according to the p olytree graph so me pro cesses may miss any input, due to the caus a l prop erty of the considered rela tionships. The above pro cedures, developed to desig n simple graphical mo dels of a s e t of pro cesses, a re summar ized in T able 1 . It is w orth highlig h ting that, once such a model has b een der ived, the a nalysis of the connections of its j -th no de can b e p ossibly exploited in order to design a reaso n able input set I ( j ) of a representation (4). 4 T op ol ogical iden tification See [12] for a comprehens ive reading ab out all the (omitted) mathematical details of this section. 4.1 Wiener filtering approac h In this section we rigoro usly inv es tigate the relia bilit y a nd corre c tness of the prop ose d g raphical mo deling pro cedure when it is emplo yed to r econstruct the top ology of a g iv en netw ork. T o this aim, co nsidering the 10 Undirected graph: 1. represent each pr o cess as a no de of the gra ph, obtaining the no des set V 2. compute the matrix D Θ as in (12) 3. c hoo se the weigh t ma trix of the g raph W V equal to D Θ 4. compute the Minim um Spanning T r e e algorithm on the weight s of W V 5. c hoo se as graphica l mo del the Minim um Spanning T ree. Directed graph: 1. represent each pr o cess as a no de of the gra ph, obtaining the no des set V 2. compute the matrix D C Θ as in (16) 3. c hoo se the weigh t ma trix of the g raph W V as in (20) 4. compute the Minim um Spanning T r e e algorithm on the weight s of W V 5. for each arc cho ose its direction a ccording to (21) 6. c hoo se as graphica l mo del the resulting p olytree. T a ble 1: Pro ce dur es to design simple graphical mo dels. ab ov e theory , we take into acco unt the class o f the ALN mo dels a s the unknown netw o rk whose top olog y has to b e identified. Mor eov er, we pro po s e a compariso n b etw een our appro ach and the results of the application of the standard MISO Wiener theor y , whic h is briefly summar ized hereafter. Definition 11 (MISO Wiener filter) . L et us c onsider the pr o c esses Y and { X i } i =1 ,...,n and assume that they ar e wide sense stationary and with zer o m e an. Then, c onsider the MISO line ar filter n X i =1 ˆ W i ( z ) X i that describ es the pr o c ess Y by { X i } i =1 ,...,n and that minimizes the mo del ling err or min W i i =1 ,...,n E  ε 2  = min W i i =1 ,...,n 1 2 π Z + π − π φ ε ( ω ) dω , (22) ε ( z ) = Y ( z ) − n X i =1 W i ( z ) X i ( z ) . (23) Such a filter is r eferr e d to as the MISO Wiener filter. Let us now cons ider the class of Acyclic Linear Netw orks (ALN), that will be us e d to highlight the correctness of both our graphical method a nd the MISO Wiener filter in iden tifying a given topo logy . Recalling the previous definition, we ca n s ee a netw ork in s uc h a class a s the interconnection of m dynamical linear SISO sys tems, whose output sig nals { X j } j =1 ,...,m are linked via relationships describ ed as:    X 1 ( z ) = e 1 ( z ) + P i : X i ∈ P ( X 1 ) G 1 i ( z ) X i ( z ) , . . . X n ( z ) = e n ( z ) + P i : X i ∈ P ( X n ) G ni ( z ) X i ( z ) , (24) where P ( X j ) denotes the set of the inputs of the j -th system. A ccording to the conv ent ion previously ado pted, let us desc ribe the links str ucture of the netw ork b y means of a graph, where the j -th no de represents the j -th sy stem and its o utput X j and wher e the edge b etw e en the j -th and the i -th no des is present, if and only if the signal X j is an input of X i , or viceversa. Ther efore, r ecalling the ALN definition, we ca n c a st the concepts of parent and child directly in ter ms of input and output of a system. 11 Definition 12. The set P ( X j ) of the signals, which ar e inputs of the j -th one, is re ferr e d t o as the set of the p ar ents of X j . The set C ( X j ) of the signals which have X j as a p ar ent is r eferr e d to as the set of the childr en of X j . A signal that has no p ar ent is said a r o ot, while a signal without any childr en is said a le af. Observe that in the ab ov e netw ork (24) the disturbances { e j } j =1 ,...,n are a set of jointly uncorr elated, z e r o- mean and wide-sense stationar y pro cesses, as defined in the ALN definition. Moreover, define the transfer function set G = { G j i ( z ) } i,j =1 ,...,n so that G j i ( z ) 6 = 0 if and o nly if X i ∈ P ( X j ) . An analogo us rea soning can be pe rformed for the path concept. Definition 13. Given two signals X j and X k , if it exists a (unique) se quenc e T j k = { τ i } i =1 ,...,m , with τ 1 = j and τ m = k , r epr esenting m > 1 signals X τ i , such that H j k ( z ) = Q m − 1 i =1 G τ i τ i +1 ( z ) 6 = 0 , then T j k is said a (dir e ct e d) p ath of length l = m − 1 b et we en X j and X k and t he r elate d H j k ( z ) is the total t r ansfer function asso ciate d to it. In the following it will b e useful to define H ii ( z ) = 1 for an y v alue of i . Because of the p olytree s tructure, developing each X j b y iteratively substituting (24) in the expres sion of its comp onents, we have that the no de X j itself never a ppea rs in the right hand side o f the expression. In fact, the ab ov e dev elopmen t leads to descr ibe e very signal via disturbances only: X j = n X i =1 H j i ( z ) e i . (25) Without loss of generality , deno te Y = X n and let us fo r m ulate the MISO Wiener filter ing problem (23)- (22) rega rding the mo delling of Y via the pro cess es { X j } j =1 ,...,n − 1 . T o this aim, let us define the following quantit ies E ( z ) = [ e 1 ( z ) . . . e n ( z )] T , W ( z ) = [ W 1 ( z ) . . . W n − 1 ( z )] T , W n ( z ) = 0 , K ( z ) = [ H 1 n ( z ) . . . H nn ( z )] T = = [ K 1 ( z ) . . . K n ( z )] T , H nn ( z ) = 1 H ( z ) =    H 11 ( z ) . . . H 1 n − 1 ( z ) . . . . . . H n 1 ( z ) . . . H nn − 1 ( z )    . Therefore, in the minimization problem (22) the term (23) can be rewritten as ε ( z ) = K T ( z ) E ( z ) − W T ( z ) H ( z ) E ( z ) = E T ( z )  K ( z ) − H ( z ) T W ( z )  . Observe that E T ( z )  K ( z ) − H ( z ) T W ( z )  is a linear combination o f the signals e j , which ar e uncorrelated. Hence, the problem bo ils down to min W E  ε 2  = min W Z + π − π n X j =1      K j ( ω ) − n − 1 X i =1 H ij ( ω ) W i ( ω )      2 φ e j ( ω ) dω , (26) which is formally equiv alent to solve the linear lea st square problem min W   K ( ω ) − H ( ω ) T W ( ω )   2 , (27) 12 being k·k the euclidean no rm. Intuit ively , then, the s ignals e j act as a “basis” for the pro cess es X j and Y , being the incor relation pro p erty forma lly analogo us to orthogo nalit y . A so particula r for mulation will lead to futher developing in the following sections. This latter formulation of the problem, suitably co m bined with concepts b orrowed from B ay e s ian Netw ork theory [1 5], can b e fruitfully explo ited to inv estigate the res ults obtained via direc t application of the MISO Wiener filtering to the top olog ical identification problem. Definition 14. The Markov Blanket of t he no de Y of a p olytr e e network is the ensemble of al l the p ar ents and childr en of Y and al l the p ar ents of no des which ar e childr en of Y . Prop osition 15. The MISO Wiener filter b etwe en Y and al l the other pr o c esses { X j } j =1 ,...,n − 1 may have non zer o c omp onent s, i.e. differ ent form the z er o t r ansfer function, only with r esp e ct to the signals b elonging to its Markov Blanket. Then, since the parents o f Y have distance equal to o ne from the parents of its children, the following result holds. Corollary 1 6. The u ndir e cte d version of t he top olo gy of a line ar p olytr e e network define d as in (24) c an b e exactly single d out via mult iple applic ations of the MISO Wiener fi lter only by using also t he c orr esp onding distanc e matrix (12) to pur ge the p ar ents of the child r en of the c onside r e d no de. 4.2 The lo cal coherence-based approac h In the following we inv estigate the results obtained by the a pplication to ALNs of the gra phical mo deling pro cedure introduced in the pr e v ious sections. Mo r eov er, a co mparison with the Wiener filtering approach will b e pr ovided in or der to highlight the adv a n tages of our metho d. T o this aim, let us co nsider a netw ork of interconnect linear systems of the form (24), a lso assuming that their link structure is describ ed by a po lytree top ology . Consider the non-ca usal dista nce matrix (12) and fo r each p ossible tuple ( i, j, k ) assume that it exists an ar bitrary small interv al I s uch that φ X i X j ( ω ) φ e k ( ω ) 6 = 0 ∀ ω ∈ I . (28) Then, the follo wing result holds. Prop osition 17. If the i -t h and j -t h systems ar e actu al ly linke d in t he r e al network, i.e., if X i is an input of X j or vic eversa, t hen t he indir e ct e dge c onne cting the i -th and the j -th no des b elongs to the MST c ompute d via the non-c aus al distanc e matrix (12) . Corollary 18 (Main result) . If c ondition (28) is satisfie d, t hen the MST c ompute d on D Θ pr ovides the indir e ct top olo gy of the original p olytr e e network. It is interesting to observe that for p olytree linear netw o r k the prop osed technique, based o n the com- putation of the MST with resp e ct to the weigh ts matrix D Θ , provide the exact identification of the indirect version of the orig ina l to p olo gy , while the application of the Wiener filter do es no t. 5 Compressiv e se nsing See [22] for details. 5.1 Preliminary considerations Recent ly , in [23] and [4] in teresting equiv a lences b etw een the identification of a dynamical net w ork a nd a l 0 sparsification problem a re highlighted, suggesting the difficulty of the re c o nstruction pro cedure [3]. 13 The main idea is to cast the problem a s the estimation of a “ sparse Wiener filter”. Recalling the problem formulation in the pr e v ious sections , g iv en a set of N sto chastic pro cess es X = { x 1 , ..., x n } , w e consider ea c h x j as the output of a n unknown dynamical system, the input of which is given by at most m j sto chastic pro cesses { x α j, 1 , ..., x α j,m j } selected fr om X \ { x j } . The choice of { x α j, 1 , ..., x α j,m j } is r ealized a c c ording to a criterion that tak es in to acc o un t the mean sq ua re of the modeling error . The par a meters m j can b e a priori defined, if we intend to imp ose a certain degr ee of spa rsity on the net work, o r a strategy fo r self-tuning ca n b e in tro duced p enalizing the int ro duction of an additional link if it do es not provide a significant reduction of the cost. F or a n y p ossible choice o f { x α j, 1 , ..., x α j,m j } , the computation of the Wiener Filter leads to the definition of a mo deling error whic h is a natural mea sure o f how the time s eries { x α j, 1 , ..., x α j,m j } ca n “ des crib e” the output x j in ter ms of predictive/smoo thing capability . Once this step has b een p erformed, each system is represented by a no de of graph a nd, then, the arcs linking a ny x α j,m k to x j are int ro duced for each no de x j . A t the e nd of this pro cedur e a gra ph, r epresenting the netw ork top ology , has b een obtained. W e will show that this way of casting the pr oblem has strong similar ities with l 0 -minimization problems, which hav e been a v ery a ctiv e topic of research in Signal Pro c essing during the la st few years. Indeed, a l 0 - minimization problem amounts to finding the “spars e st” solutio n o f a set of linear equatio ns [4]. Unfortunately , with no additiona l ass umptions on the solution, the pr oblem is combinatorially intractable [4]. This has prop elled the s tudy of rela xed problems in volving, for example, the minimization of the no rm- 1 which is a conv ex pro blem a nd it is known to provide solutions with at least a cer tain order of s pa rsity [23]. W e exploit the similarity of the tw o problems bo rrowing some algorithmic to ols which hav e b een developed in the are a of Signal Pro cessing adapting them to o ur needs. W e start intro ducing a pre- Hilb ert space for wide-s ense sto chastic pro cesses , wher e the inner pro duct defines the notion of p erp endicularity b etw een tw o sto chastic pro cesses. This allows us to sea mlessly adopt well-known greedy algorithms for a (sub optimal) s olution of the problem, such as Ma tc hing Pursuit (MP) or Ortho g onal Least Squar e s (OLS), which are based on iterated pro jections. 5.2 Pre-Hilb ert space Hereafter a metho d to build up a pre-Hilb ert space for stationar y ra ndom pro cesses is pr esent ed. The mathematical details are omitted for the sake of s implicit y (see [22]) Definition 19. L et e ( t ) := ( e 1 ( t ) , .., e N ( t )) T a N -dimensional time-discr ete, zer o-me an, wide-sense station- ary r andom pr o c ess define d for t ∈ Z , su ch that, for any i, j ∈ { 1 , ..., N } , the p ower sp e ctr al density Φ e i e j ( z ) exists on the unit cir cle | z | = 1 of the c omplex plane and is r e al-r ational. F ormal ly, we write Φ e i e j ( z ) = A ( z ) B ( z ) for i, j = 1 , ..., N with A ( z ) , B ( z ) r e al c o efficient p olynomials su ch that B ( z ) 6 = 0 for any z ∈ C , | z | = 1 . W e say that e is a ve ctor of r ational ly r elate d r andom pr o c esses. Definition 20. Define the sets 0 F := { W ( z ) | W ( z ) is a r e al-r ational sc alar function of z ∈ C define d for | z | = 1 } 0 F m × n := { W ( z ) | W ( z ) ∈ C m × n and any of its entries is in 0 F } . Prop osition 21. The ensemble ( 0 F e, + , · , R ) is a ve ctor sp ac e. Definition 22. W e define a sc alar op er ation < · , · > on 0 F e in the fol lowing way < x 1 , x 2 > := R x 1 x 2 (0) . Prop osition 23. The set 0 F e , along with the op er ation < · , · > is a pr e-Hilb ert sp ac e (with the te chnic al assumption that x 1 and x 2 ar e the same pr o c esses if x 1 ∼ x 2 ). 14 Definition 24. F or any x ∈ 0 F e we denote the norm induc e d by the inner pr o duct as k x k := < x, x > . W e provide an ad- ho c v ersion of the Wiener Filter (gua ranteeing that the filter will b e real ra tional) with an int erpretation in terms of the Hilb ert pro jection theo rem. Indeed, g iven signals y , x 1 , ..., x n ∈ 0 F e , the Wiener Filter estimating y from x := ( x 1 , ..., x n ) can b e in terpreted as the op era tor that determines the pro jection of y onto the subspace 0 F x . Prop osition 25. L et e b e a ve ctor of r ationally r elate d pr o c esses. L et y and x 1 , ..., x n b e pr o c esses in the sp ac e 0 F e . Define x := ( x 1 , ..., x n ) T and c onsider the pr oblem inf W ∈ 0 F 1 × n k y − W ( z ) x k 2 . (29) If Φ x ( ω ) > 0 , for al l ω ∈ [ − π , π ] , the solution exists, is unique and has the form W ( z ) = Φ y x ( z )Φ xx ( z ) − 1 . Mor e over, for any W ′ ( z ) ∈ 0 F 1 × n x , it holds that < y − W ( z ) x, W ′ ( z ) x > = 0 . (30) Let e b e a v ector o f r ationally rela ted pro cesses. Consider a s e t X := { x 1 , ..., x n } ⊂ 0 F e of n rationally related pr o cesses with zero mea n and known (cros s)-power sp ectral densities Φ x i x j ( z ) . F or any given pro cess x j , with j ∈ { 1 , ..., n } , fix m j ∈ N suc h that 0 < m j < n − 1 . W e intend to identify the m j pro cesses x α j,k , α j,k 6 = j , with k = 1 , ..., m j , which, filtered by suitable rational transfer functions W j,α j,k ( z ) , provide the bes t estimate of X j in the s e nse of the mean squar e s. The mathematica l formulation of the problem is the following: min α j, 1 ,...,α j,m j 6 = j W j,α j,k ( z ) ∈ 0 F      x j − m j X k =1 W j,α j,k ( z ) x α j,k      2 , (31) where every W j,α j,k ( z ) , with k = 1 , ..., m j , is a p ossibly non-c a usal trans fer function. Given any set { α j,k } m j k =1 , the ab ove problem is immediately solved by a multiple input Wiener filter. It is the determination of the parameters α j,k that makes the problem com binatorial. The signals x α j,k ∈ X repr esent the m j signals with the “stronges t affinity” with x j according to a least square criterion. Moreov er, a s previously o bserved, the result of this optimization problem lends itself to a natural interpretation in terms o f Graph Theo ry . Recalling the approach o f the previous sectio ns , as s ume that the pro cesses x j , j = 1 , ..., n , are repres en ted r esp e c tively by the no des N j , j = 1 , ..., n , in a graph. F or j = 1 , ..., n , draw an or ie nted edge from each of the no des N α j,k , k = 1 , ..., m j , to the no de N j . F o llowing this pro cedure, the o btained edge s et ca n be readily visualized and provides a qualitative description of the in ternal s tructure of the whole system in terms of the mo st r elev ant linear relations. If the goa l is to mo del a complex sys tem as a netw ork of suitable agents, the para meters m j ’s, defining the maximum num ber of ent ering e dg es for the no des N j ’s, can be used to a djust the deg ree of complexity of the gr a ph. Indeed, when m j = n − 1 for any j = 1 , ..., n , we exp ect to obtain a co mplete gra ph, that is all p ossible links will b e present. Conv ersely , a different ch oice for the pa rameters m j ’s will lea d to a reduction of co mplexity of the identified net work mo del. 5.3 Links with compressiv e sensing In this section we highlight the connections b etw een the problem of reconstructing a topolo g y and the com- pressive sensing problem. Such a co nnection is p ossible b e c ause of the pre-Hilber t structure constructed 15 befo r e. Indeed, the concept of inner pro duct defines a notion of “pro jection” amo ng stochastic pro cesses and makes it p os s ible to seamlessly imp or t to ols developed for the compressive sensing pr oblem in order to tackle that of iden tifying a top ology . In the recent few years sparsity problems have a ttra cted the attention of r esearchers in the are a of Signal Pro cessing. The rea s on is mainly due to the p ossibility of repres en ting a signal using only few elements (w ords) of a redundan t base (dictionary). Applications are n umerous, ra nging fro m data-compression to high-resolution interpolatio n, and noise filtering [6 , 35]. There are man y formalizations of the problem, but one of the most co mmon is to cast it as min w k x 0 − Ψ w k 2 sub ject to k w k 0 ≤ m (32) where n < p ; x 0 ∈ R p ; Ψ ∈ R p × n is a matrix, who se co lumns r epresent a redundant base employ ed to approximate x 0 , and the “zero-nor m” (it is not a c tua lly a norm) k w k 0 := |{ i ∈ N | w i 6 = 0 }| (33) is defined by the num ber of non-zero en tries of a v ector w . It can b e sa id that w is a “simple” way to expres s x 0 as a linear combination of the columns o f Ψ , where the co ncept of “simplicity” is given by a constraint on the num ber o f non-zero entries of w . F o r each j = 1 , ..., n , define the following sets: W ( j ) = { W ( z ) ∈ 0 F 1 × n | W j ( z ) = 0 } , (34) where W j ( z ) deno tes the j -th co mpo nent of W ( z ) . F or any W ∈ W ( j ) , define the “zero- no rm” a s k W k 0 = { # of entries such that ∃ z ∈ C , W i ( z ) 6 = 0 } and define the random vector x = ( x 1 , ..., x n ) T . (35) The pro blem (31) can b e for mally cast as min W ∈W j k x j − W x k 2 sub ject to k W k 0 ≤ m , (36) which is , from a for mal p oint of view, equiv alent to the standar d l 0 problem as defined in (32). 5.4 Solution via sub optimal algorithms The pro blem o f netw o rk identification/complexit y reduction we hav e formulated in this pap er is equiv alent to the pro blem of determining a sparse Wiener filter, as explained in the pre v ious section, once a notio n of orthogo na lit y is int ro duced. This formal equiv a lence shows how the problem of identifying a top olog y can immediately inherit a set of pr actical to o ls already developed in the ar e a of compress ing sensing. Here we present, as illustrative examples, mo difications o f algorithms a nd strategies, well-kno wn in the Signa l Pro cessing communit y , which can be adopted to obtain sub optimal so lutions to the problem of iden tifying a net work. While formally identical to (32), the problem of iden tifying a top ology as cast in (36) still has its own characteristics. Since the “pr o jection” pro cedure in (36) is given by the estimation of a Wiener filter, it is computationally mor e expensive than the standa rd pro jection in the space of vectors of rea l num ber s. F o r this reaso n g reedy algorithms offer a go o d appro ach to ta c kle the problem since speed b ecomes a funda- men tal factor. Moreover, since the complexity of the netw o rk mo del is the final goal, greedy a lgorithms a re a natural solution allowing one to specify explicitly the connection degree m j of every no de x j . This feature is in general not provided b y other alg orithms. As a n alter na tiv e approach to gr eedy algo rithms w e a lso descr ibe a strategy based on iterated r eweigh ted optimizations as des crib ed in [4]. 16 SLE BUD KFT PEP KO PG MO SO AEP CPB HNZ MMM VZ T BMY MDT S PFE ATI COP IP WY AA DD SLB BHI HAL XOM EP WB HON CAT GE TYC COV BA RTN GD NSC BNI UPS DOW ROK CCU DIS BAX TWX CMCSA FDX HD CBS CVS TGT MCD INTC GM ALL WMT TXN MSFT CSCO F HIG AIG USB CI RF WFC BK BAC NYX AXP COF JPM C MS HPQ ORCL XRX AAPL DELL IBM EMC ABT JNJ MRK AES ETR EXC AMGN GOOG AVP CVX WMB UNH CL GS LEH MER UTX Figure 1: The tree str ucture obtained using a coher ence-based a nalysis. Every no de repr esents a sto ck and the color repres e n ts the business sector it b elongs to. The considered sectors are Basic Mater ial (y ellow), C onglom erates (white), H ealthc are (pink), Tran sporta tions (dar k blue), Te chnolo gy (red), Capit al Goods (or ange), Util ities (brown tin ts), Cons umer (vio let tints), Finan cial (g reen tints), Energy (grey tints) Servic es (light blue tints). Using the industry clas s ification given by Go o gle, the Finan cial sector has also been differentiated among Insurance Co mpanies (light green), Banks (a verage green) and In vestmen t Co mpanies (dark gr een); Servi ces have b een divided in Information T echnology (cyan) and Retail (aquama rina), Con sumer in F oo d (plum) a nd Personal-care (purple); Energy in Oil & Gas (dar k grey) and W ell Equipmen t (ligh t grey); Utilities in Electrica l (dar k brown) and Natural Gas (lig ht br own) . The complete sto ck list is rep or ted in T able 2. 6 Application to clustering and graphical m o del i ng 6.1 High-frequency sto c k mark et analysis A collection of 100 sto cks of the New Y ork Sto ck Exc hange has b een observed for four weeks (nineteen market days), in the lapse 03/03 /2008 - 03 /28/2 008 s ampling the prices every 2 min utes. The sto cks hav e b een chosen as the first 10 0 sto cks with highest trading v olume acco rding to the Standard & Poor Index at the first da y of observ ation. An a priori organization o f the mar ket ha s b een ass umed in accordance with the sector and industry gr oup classification pr ovided by Go o gle Finance R  , tha t is also the so urce of our data. W e under line that in this pap er we ar e mainly concerned with sectors, but in some cases we hav e also taken into a c c ount the industry group classification to refine the results. The whole observ ation horizon spans the month o f March. W e hav e not considere d reaso nable the h ypo thesis that the mutual influences among the c o mpanies are stationary during this time p erio d. Ho w ever, it is w orth observing that the market session data a re na tur a lly divided int o subp erio ds, namely w eeks and days and discontin uities a re present. Indeed, due to the pre- and po st-market sessio ns , there is a discontin uit y b etw een the end v alue o f a day and the op ening price of the next one. There fo re, in our analysis, we hav e follow e d the na tural approa ch of dividing the historical series in to nineteen subp erio ds corres p onding to sing le days. Mo reov er, we ha ve assumed the working h yp othesis that the relations betw een the companies can be a pproximated with stationary ones during a single da y . Hence, we have p erfor med the multiv ariate analys is for each subp erio d computing the coherence-based dista nce (11) among the stocks. Finally , we ha ve av er aged such dista nce s o ver the who le obser v ation horizon and the related results have b een exploited to ex tr act the MST, providing the corresp onding market s tructure. W e find useful to remar k that the computation of the dista nces for small data sets is als o computationally b etter per forming. The cor resp onding r esults a r e illustrated in Figure 1 . Every no de r e pr esents a sto ck a nd the color re presents the business sector o r industry it b elongs to. W e note that the sto cks are q uite sa tisfactorily group ed according to their bus ines s sectors. W e stress that the a prio r i classification in sector s is not a hard fact b y itself and we are not tr y ing to ma tc h it exactly . A company co uld well b e ca tegorized in a 17 sector becaus e of its business , but, at the same time, could show a b ehaviour similar to and explainable through the dynamics of other se c tors. Actually , we would b e very int erested into finding results o f this k ind. Indeed, in those very cases, our quantitativ e analysis would provide the gre a test contributions detecting in an ob jective w ay something which is “counter-in tuitiv e”. Thus, we just use such a priori classificatio n as a to ol to check if the final topo logy makes se ns e a nd if, at a ge ne ral level, the cohere nce appro ach perfor ms better. Despite this disclaim, it is worth noting that the Fina ncial (green tints), Cons umer (violet tint s), Basic Mate rials (yellow) , Energy (grey tints) and Tr anspor tation (dark blue) sectors are a ll p e r fectly group ed, with no exc eptions. W e note a sub c lusterization o f the F inanci al sector , as well. The C onsume r sector s hows another pr ominent s ub clusteriza tion in the Fo od (plum) and Perso nal/He althcare (purple) industries, while the Ener gy sector presents an evident sub clusterization in to the Oil & Gas (dar k gr ey) and Oil Well Equip ment (light grey). The close presence o f co mpanies of a different sector and industry , Utili ties/N atural Ga s (light brown ) is observed as well. The other Utilit ies/El ectricity co mpanies (dark brown) are, int erestingly , a different gro up. W e also o bserve a big cluster of companies classified as Servi ces (light blue tints). W e hav e differentiated them in the tw o industries Ret ail and Inf ormati on Techn ology us ing tw o slightly different colo rs, r espe c tively aqua marine and cyan. W e also note the presence of three Ser vices companies which are isolated from the other ones: V [V er izon], T [A T&T], and S [Sprint ]. All of them a re telephone companies. This might suggest that this industry sho uld show at least a slightly different dyna mics from the other ser v ice companies. Note also how the Techn ology s ector (red) is almost per fectly g roup ed a nd how IBM , a n IT compa ny , even though cla ssified as a Se rvices company , is lo cated in it. Finally , the t w o only automobile companies GM and F [F or d] happen to be linked together. The a nalysis of this four w eeks of the month of March cleanly shows a taxonomic arr angement of the sto cks ev en though the choice of a tr ee structure might have se emed quite reductive at first thought. 6.1.1 Benefits o f the dynamical mo deli ng In this section w e highlight the improv emen ts of the dynamical appro ach, i.e. the co herence-based distance (11), against the static pro cedure, i.e. the traditional correlation- based distance used, for example, in [18]. How ever, the application of the cor relation-base d a nalysis to the entire da ta set would be meaningless, since the sto cks a r e sampled at high frequency ov er a long observ a tio n horizo n. Hence, to av o id the no n-stationary problem, we p erform the pr o c e dure in the daily subp erio ds, av eraging the corresp onding distances. Figur e 2 illustrates a compar is on b etw ee n the “hea t map” cor resp onding to the co r relation-bas e d dis ta nce matrix and the coherence- based o ne. The energ y of the erro r a s so ciated to the dyna mical mo deling approa ch is exp ected to b e low er than the s ta tic one. Conversely , the cor relation-base d distance is exp e c ted to b e hig her than the cohere nce-based metric. Such a pr o per t y is highlighted in Figur e 2 by the darker lo o king of the map (b) against map (a). Indeed, it shows that the dynamical mo deling approach is able to describ e a lar ger v ariety of dependencies , es pecia lly rela tions inv o lving the presence o f delays. The Spear ma n index has b een ev aluated from these s et of data obtaining σ sp ≃ 0 . 158 . The MST derived from those corr elation-based distance s is illustrated in Fig ure 3. At a v ery first loo k , a lesser capa bilit y of grouping the stocks according to their sectors already app ears. Mo reov er, a further analysis r eveals, for instance, that the sub clusterizatio n of the Financ ial s e ctor highlighted in the dynamica l approach is no more present in the correla tion-based MST. T h e same happ ens fo r the sub clusterizations of the Con sumer a nd Ene rgy s ectors. The higher consistency of the coherence-based metric ca n b e also underlined considering the time evolution of the single distances. In Figure 4 we hav e co nsidered a test node b elonging to the Te chnolo gy sec tor a nd we hav e rep or ted alo ng the time hor izo n b oth its distances with resp ect to other Tec hnolog y sto cks and the distances r elated to the F inanci al o nes. As exp ected, the dynamical modeling appro a ch is more consistent and ro bust in detecting similar ities among sto ck time series b elonging to the same secto r s. 18 10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 (a) (b) Figure 2: A g raphical r epresentation of the averaged co rrelation (a) and coherence (b) dista nce matric e s . An y distance matrix has 100 × 100 en tries representing the distances o f a n y p ossible couple of the 100 sto c ks in our analy sis. These v alues have b een repres e nted using a grey scale: the lighter the spo t, the la rge the distance b etw een the tw o sto cks. The interpretation of b oth the distances in terms of a mo deling erro r shows, as exp ected, a b etter modeling capabilit y for the dynamic a pproach. This can be considered an evidence for the presence of propagative phenomena in the netw o rk, at the considered time scale o f 2 minutes. ALL HD HIG COF CI JPM WFC BAC RF AIG AXP MS BK F GM TWX MER NYX GE EMC MMM RTN S HON MDT HPQ ROK DOW TYC UTX BA CAT SLE KFT PEP BUD MO CBS CCU CPB KO HNZ CL BMY UNH ATI IP AA DD WY BNI CVX COP XOM SLB HAL EP SO VZ T BAX ABT MRK WB XRX GD NSC COV DIS CMCSA AMGN FDX CVS TGT INTC MSFT CSCO ORCL WMT MCD USB C GS TXN AAPL IBM PFE JNJ AES AEP GOOG AVP BHI UPS DELL PG ETR EXC LEH WMB Figure 3 : MST ass o ciated to the correla tion-based metr ic . The distances have b een computed in the daily subper io ds and then av eraged as in Figur e 1. An increased num b er of sto cks group ed with o thers b elonging to a different s ector is exp e r ienced. F or instance, HP Q [Hewlett-Pac k ard] is no mor e g roup ed with the other Techn ology sto cks. One might also observe that GM and F a re not directly linked or that the Cons umer sto cks (cyan) are almost spread over the whole gr aph. F ur thermore, the s ubclusteriza tion of the Financial se c to r, highlighted in the Figure 1, is not present. 19 0 2 4 6 8 10 12 14 16 18 20 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Day Distance 0 2 4 6 8 10 12 14 16 18 20 0.12 0.14 0.16 0.18 0.2 0.22 0.24 0.26 0.28 0.3 0.32 Day Distance (a) (b) Figure 4: The distances of 5 financia l s tocks (blue) and 5 technological sto cks (red) from a sample tec hnological sto ck hav e been computed on a daily base. The figur e repr esents bo th the corre lation distances (a) and the coherence ones (b). As it can b e noted, the corre la tion distance is mor e consistent a nd robust in detecting similarities among sto cks b elonging to the same sectors . 6.2 Graphical mo dels of Europ ean climatic regions In this section we provide an illustrative exa mple to s how the effectiveness of the pr o po sed g raphical mo deling approach (see [1 2] for a further reading on this c ase). T o this aim, we re p or t the results obtained by ana lyzing a set of real data. A cr ucial assumption lies in the adoptio n of linea r mo dels for the mo deling pro cess. Indeed, the application of our techni que to real data is exp ected to be a ble to provide information ab out the linear comp onent of the dynamics, while the nonlinear part will necessarily be embedded as no ise. In other words, the corr esp onding re s ults are exp ected to be as w orse as the nonlinear dynamics prev ail in the real net w ork. How ever, though a linear analysis may b e not sufficient to fully unders tand the pro cess dynamics, w e hig hlight that it provides a firs t step to der ive information ab o ut the system, especia lly ab out its internal connections. Indeed, we remind that no a prior i knowledge is assumed, th us ours is a fully blind appro ach. In general, it is p ossible to incor por ate a priori knowledge int o the mo deling proce dure if, for ex ample, the dynamics is known to hav e a cer tain structure. How ever, the development of similar pro cedur es go es b eyond the sco pe of this paper and it will b e matter of future research. In the following we address a collec tio n of meteor ological time series sampled in 50 sites lo ca ted in as muc h Europ ean towns, listed in T able 3. Our main goal is obtaining infor mation ab out the similarities among the temper a tures in thos e towns, or , rather, deriving a suitable top ologic a l mo del. The temper atures in the ab ov e loca tions are observed ov er a o ne year time hor izon cov er ing the whole 19 8 8. Suc h time s eries are hourly sampled, providing a high frequency da ta set 1 . Notice that consequently b oth daily and s easonal trends ma y be p ointed out from the observ ations. T o this re gard, any deterministic tr e nd has to b e r emov e d from the original time series. Hence, let { Y k ( n ) } n =1 ,..., 8760 , k = 1 , . . . , 50 , b e the hourly sampled temp eratures in the considered towns. Any trend different from the daily o ne c a n b e singled out just by observing the mean v alue ov er a 24 ho urs span, i.e. S k ( n ) = 1 24 11 X i = − 12 Y k ( n + i ) , n = 1 , . . . , 8 760 , k = 1 , . . . , 5 0 , assuming Y k ( n ) = 0 if n < 1 or n > 8760 as a working hypothesis. Ther efore, we provide the rejection o f such seas onal trends just b y defining the new time series X k ( n ) := Y k ( n ) − S k ( n ) , k = 1 , . . . , 50 . 1 In ternat ional Surface W eather Observ ations 1982-1997, V olume 3 (Europe) , 1998 Asheville, N.C 20 T a ble 2: Lis t of the c ompanies considere d in the analysis. Name Code Sector 3M Company MMM Conglomerates Abbott Laboratories ABT Healthcare Aes Corporation AES Utilities Alcoa Inc. AA Basic Materials Allegheny T echn ologies Inc. A TI Basic Materials Allstate Corporation ALL Financial Altria Group MO Consumer/Non-Cyclical American Electric Pow er AEP Utilities American Express AXP Financial American Interna tional Group AIG Financial Amgen Inc. AMGN Healthcare Anheuser Busch BUD Consumer/Non-Cyclical Apple Inc. AAPL T echnology A T& T T Services A von Products A VP Consumer/Non-Cyclical Baker H ughes Inc. BHI Energy Bank of America BAC Financial Bank of New Y ork Mellon BK Financial Baxter Inter national BA X Healthcar e Boeing BA Capital Goo ds Bristol Myer s Squibb BMY Healthcare Burlington Norther n Santa F e BNI T ransportation Campbell Soup CPB Consumer/Non-Cyclical Capital One Financial COF Financial Caterpillar Inc. CA T Capital Goods CBS CBS Services Chevron CVX Energy CIGNA CI Financial Cisco Systems CSCO T echno logy Citigroup Inc C Financ ial Clear Channel Communications CCU Services Coca-Cola KO Consum er/Non-Cyclical Colgate Palmolive CL Consumer/Non-Cyclical Comcast CMCSA Services Conoco Phill i ps COP Energy Covidien COV Hea lthcare CVS Caremark CVS Services Dell Inc DELL T echnology Dow Chemical Company DOW Basic Materials E.I. du Pont de Nemours DD Basic Materials El Paso EP Utili ties EMC EMC T echn ology Ente rgy ETR Utilities Exelon EXC Utilities Exxon Mobil XOM Energy F edEx FDX T ransportation F ord Motor F Consumer Cyclical General Dynamics G D Capital Goods General Electric GE Conglomerates General Motors GM Consumer Cyclical Name Code Sector Goldman Sachs Group GS Financial Google I nc. GOOG T echnology Halliburton HAL Energy Hartford Financial Services H IG Financ ial H. J. Heinz H N Z Consumer/Non-Cyclical Hewlett-Pac k ard H PQ T echnology Home Depot HD Services Honeywell I nte rnational HON Capital Goods Intel INTC T ec hnology Intern ational Business Machines IBM Services Intern ational Paper IP B asic Materials Johnson & Johnson JNJ H ealthcar e JPMorgan Chase JPM Financial Kraft F o ods KFT Consumer/Non-Cyclical Lehman Brothers Holding LEH Financial McDonald’s MCD Services Medtron ic MDT Healthcare Merc k MRK Healthcar e Merril Lynch MER Financial Microsoft MSFT T echnology Morgan Stanley MS Financial Norfolk Souther Group NSC T ransportation NYSE Euronext NYX Financial Oracle ORCL T ec hnology Pespi PEP Consumer/Non-Cyclical Pfizer Inc. PFE Healthcare Procter & G amble PG Consumer/Non-Cyclical Raythe on R TN C onglomerate s Regions Financial RF Financial Rockwell Automation ROK T echnology Sara Lee SLE Consumer/Non-Cyclical Schlumberg er Limited SLB Energy Southern SO Utilities Sprint Nextel S Services T arget TGT Services T exas Instruments Inc. TXN T echnology Time W arner TWX Services Tyco International TYC Conglomerates U. S. Bancorp USB Financial United Parc el Service UPS Tr ansportation United T echnologies UT X Cong lomerates UnitedHealth Group Inc. UNH Financial V erizon Communications VZ Services W achovia WB Financial W al-Mart Stores WMT Ser vices W alt Di sney DIS Services W ells F argo WFC Financial W eyerhaeuser Company WY Basic Materials Williams Companies WMB Utilities Xerox XRX T echnology This tec hnique to detrend data is a na logous to the pro cedure rep o rted in [28]. It is w orth noticing that the new time s eries mea ns are approximately eq ua l to zero, since the tempera tur e v ar iation over a 24 hours time span is fair ly distributed a round the mean v a lue, as exp ected by such time ser ies. Observe that also this daily trend s ho uld be ass umed as a deterministic p erio dic comp onent , but, due to the nature of the data, its range is expected to v a r y along the year. Theor etically , we could apply a “window ed” appr o ach to p oint out how the daily b ehavior c hanges, but the choice o f the related time length should b e derived by the seasona l trends. Therefore, for the sake of the simplicity , in this example we just neglect to remov e the daily trend, considering it a s part of the sto chastic v a r iation of the signal. In the following analysis we consider a given data set and we do not requir e any on-line computation. Moreov er, o nly the non-causal appro ach will be presented, s ince the causal one just differs for the numerical to ols used to compute the distance matrix. The coherence functions of the detrended time ser ies { X k } k =1 ,..., 50 are numerically ev aluated by employing the W elc h algo rithm [3 4]. Then, the corr esp o nding distances are computed a ccording to (11), pr oviding a quantit ative des cription of similarities and connections among the tempera tures in the consider e d to wns during the year 1988. Such infor mation can b e exploited to c ho ose the best sour c e of each time series according to the minim um mo deling error approach and to derive a connected tr e e top ology acco rding to the pro cedure summariz e d in T able 1. The corresp onding graphica l mo del is illustrated in Figure 5. Ev en though the appro ach is blind a nd it do es no t exploit a n y a prior i knowledge, the mo deled top ology absolutely turns out r easonable. T owns b elonging to the same clima tic r egion are co r rectly linked to gether, 21 Ab erdeen ABR Bordeaux BRX Berli n BRL Barce l ona BRC Belfast BLF Brest BRE Bremen BRE Granada GRN Birmingham BR M Chartres CHR Dre sd e n DRS Madri d MDR Bournemouth B MT Dijon DJN Dusseldorf DSS Malaga MLG Bristol BRS Ly on LYN F rankfurt Am Main FRN Salamanca SLM Cardiff CRD M arseille MAR Hamburg HMB Santiago De Comp ostela SDC Carlisle CRL Montpe l lier M NT Hannov er HAN Sevilla SVL Edinburgh EDN Nanc y NCY Hei delb erg HDL V alen c ia VAL Exeter EXE Nan tes N AN Kiel KIE Zaragoza ZAG Glasgo w GSW Nic e NIC Leipzig LPZ London LND Orl eans ORL M an nheim MNH Plymouth PLY Paris P RS Mun i ch MUN Manchester MAN Strasbourg STR Stu ttgart S TG Nottingham NOT T oul ouse TOU T a ble 3: The 50 Europ ean towns considered in Section 6.2 a nd their lab els. ABR BLF GSW EDN CRL MAN NOT BRM LND BMT CRD BRS BRE NAN ORL CHR PRS NCY DJN L YN BRX TOU MNT MAR NIC ZAG SLM SDC MDR SVL GRN VAL MLG BRC STR STG MNH FRN HDL MUN DRS BRL LPZ HAN DSS BRE HMB KIE PL Y EXE Figure 5: Solid plus dashed lines: the tree top o logy obtained by the application of the undir e cted pro cedure of T able 1 for the towns rep orted in T able 3. Solid lines: connections cor resp onding to the minima of ea ch row of the undirected dista nce matrix. also forming consistent clusters in term of the links asso ci a ted only to the minima o f the r ows of the distance matrix. Moreover, s ome connectio ns representing the pr opagative effects o f the weather are provided, as well. Hence, the co rresp onding top ology app ears consis tent with the problem of deriving a reduced optimal mo del of the whole system. F or instance , obser ve that, according to o ur a pproach, in order to derive a suitable reduce linear mo del to des crib e the temp eratures in Edinburgh, we should consider a MISO linear filter with the temp eratures o f Ab erdeen, Carlisle and Glasgow as inputs. References [1] Georg e Box, Gwilym M. Jenkins, and Greg ory Reinsel. Time Series Analysis: F or e c asting & Contr ol (3r d Edition) . Pren tice Hall, F ebruary 1 9 94. [2] A. Brovelli, M. Ding, A. Ledb erg, Y. Chen, R. Nak amura, a nd S. L . Bressler. Beta oscilla tions in a large-s c a le sensorimo tor co rtical netw ork: directional influences revealed by gr anger causality . Pr o c Natl 22 A c ad Sci USA , 1 01(26):98 49–98 54, June 2 004. [3] E. J. Candès a nd T erence T ao. Deco ding by linear progra mming. IEEE T r ansactions on Information The ory , 5 1(12):420 3–421 5, 200 5. [4] Emmanuel Candès, Michael W a kin, and Stephen Boyd. Enhancing sparsity by reweigh ted l1 minimiza- tion. Journal of F ourier Analysis and Applic ations , 14:87 7–905 , 20 08. [5] Marinaz z o D., P ellicoro M., a nd Strama glia S. Kernel metho d for nonlinea r gra nger causality . Phys. R ev. L ett. , 100:14 4103, 200 8. [6] I. Daubechies, R. De V ore, M. F or nasier, C. S&ido t, et al. Iter atively r eweigh ted least squares minimiza- tion for spar se recov ery . Communic ations on Pure and Applie d Mathematics , 63(1):1–38, 2009 . [7] R. Diestel. Gr aph The ory . Spring er-V erlag, Ber lin, Germany , 200 6. [8] Andre F ujita, Luciana Gomes, Joa o Sato, Rui Y amaguchi, Car los Thomaz, Mari Sogayar, and Satoru Miyano. Mu ltiv ariate gene expression analysis reveals functional connectivity changes betw een no r- mal/tumoral prosta tes. BMC Systems Biolo gy , 2 (1), 2008 . [9] Andre F ujita, J. .R. Sato, H. M. Garay-Malpartida, R. Y amaguchi, S. Miyano, M. C. Sogay ar, and C. E. F er reira. Mo deling g ene expression reg ulatory netw orks with the sparse vector autor egressive mo del. BMC Systems Biolo gy , 1 , 2007 . [10] G. Inno centi and D. Materassi. A mo deling appro ach to mul tiv ariate analy s is and clusterization theory . Journal of Physics A: Mathematic al and The or etic al , 41 :20510 1, 2008 . [11] G. Innocenti a nd D. Ma terassi. T op o logical prop erties in identification and mo deling techniques. Pr o c. of the 17th IF AC W orld Congr ess, Se oul, S out h Kor e a , 2 008. [12] G. Innocenti and D. Materassi. Modeling the top olog y of a dynamical netw ork via wiener filtering approach. Submitte d to Automatic a , 200 8 (first submission), 2011 (under review). [13] G. Inno centi and D. Materassi. T op olo gic al Pr op erties in Ident ific ation and Mo deling T e chniques . con- tribution to the bo ok Mo delling, Estimation and Co n trol o f Net w orked Complex Systems edited by A. Chiuso, L. F ortuna, M. F rasca, A. Rizz o , L. Schenato and S. Z ampieri for Spring er ed., 200 9. [14] T. K ailath, A.H. Say ed, and B. Hassibi. Line ar Est imation . Prentice-Hall, Inc., Upp er Saddle River, NJ, USA, 2 000. [15] D. Koller and N. F riedman. Pr ob abili stic Gr aphic al Mo dels: Principl es and T e chniques . MIT Press, 20 09. [16] D Lib erzon and D Nesic. A unified framework for des ign a nd analysis of netw orked and quantized control systems. IEEE T r ansactions on Automatic Contr ol , 54:732– 747, 2009 . [17] R. Madan, N.B. Mehta, A.F. Molisch, and Jin Zhang. Energy-efficient decentralized co o per ative r outing in wireless netw orks. IEEE T r ansactions on Automatic Contr ol , 54 :5 12–52 7, 2 009. [18] R.N. Mantegna a nd H.E. Stanley . An Intr o duction t o Ec onophysics: Corr elations and Complexity in Financ e . Cambridge Universit y Press, Cambridge UK, 200 0. [19] D. Mater assi and G. Inno centi. T op olog ical identification in netw orks of dy na mical systems. Pr o c. of the 47th IEEE Confer enc e on De cision and Contr ol, Cancun, Mexic o , 2008. [20] D. Matera ssi and G. Inno cent i. Un veiling the co nnectivit y structure of financial netw orks via high- frequency analys is . Physic a A , 388:3866 –3878 , 2009 . 23 [21] D. Matera ssi and G. Inno cent i. T o po logical iden tification in netw or ks o f dynamical systems. IEEE T r ansactions on Automatic Contr ol , 55(8):1860 –1871 , 20 1 0. [22] D. Matera s si, G. Inno centi, and L. Giarré . Reduced complexity mo dels in the ident ification o f dynamical net works: Links with spar sification problems. Pr o c. of the 48th IEEE Confer enc e on De cision and Contr ol, Shanghai, P.R. China , 2 009. [23] D. Nap oletani and T. Sauer . Reconstructing the top ology of spar sely connected dynamical net works. Phys. R ev. E , 7 7:0261 03, 200 8. [24] M.J. Naylora, L.C. Roseb, and B.J . Mo yle. T o po logy of fore ig n exchange markets using hierarchical structure metho ds. Physic a A , 382:199– 208, 2 007. [25] R. Olfati-Sab er. Distributed k alman filtering for sensor netw orks. In Pr o c. of IEEE CDC , pages 5492– 5498, New Or leans, 2007. [26] R. Olfati-Sab er and R. M. Mur ray . Consensus problems in netw orks o f agents with switching top olog y and time-delays. Automatic Contr ol, IEEE T r ansactions on , 49(9):152 0–153 3, 20 04. [27] Reza Olfati-Sab er. Flo cking for m ulti-agent dynamic systems: Algorithms and theor y . IEEE T r ansactions on Automatic Contr ol , 51:401– 420, 2006 . [28] B. Podobnik and H.E. Stanley . Detrended cr o ss-corr elation analysis: A new metho d for analyzing t w o non-stationar y time ser ies. Phys. R ev. L ett. , 100 :08410 2, 2 008. [29] A.H. Sayed and T. Kailath. A sur vey of sp ectral factoriza tion metho ds. Numeric al Line ar Algebr a with Applic ations , 8 :467–4 69, 200 1. [30] A.N. Shir yaev. Pr ob ability . Springer-V erlag, New Y o r k, 1995 . [31] M. Timme. Revealing netw o r k connectivity from re s po nse dynamics. Phys. Rev . L ett. , 9 8(22):224 101, 2007. [32] M. T umminello, C. Coronnello , F. Lillo, S. Miccic hé, and R.N. Mantegna. Spanning trees and bo o tstrap reliability estimation in corr elation-based net works. Int. J. of Bifur c ations & Chaos , 17:23 19–23 29, 2007. [33] D. Urba n a nd T. Keitt. La ndscap e c o nnectivit y: A graph- theoretic p ersp ective. Ec olo gy , 82(5):120 5– 1218, 20 01. [34] William T. V etterling and Brian P . Flannery . Numeric al R e cip es in C++: The Art of Scientific Com- puting (3r d e dition) . Cambridge University Press, F ebruar y 2007. [35] D. P . Wipf and S. Nag ara jan. Iterative reweigh ted l1 and l2 metho ds fo r finding sparse solutions. In Signal Pr o c essing with A dap tive Sp arse S tructur e d R epr esentations , Rennes, F rance, April 2009. 24