Reconstruction of Multidimensional Signals from Irregular Noisy Samples

1 Reconstruction of Multidimensional Signals from Irre gular N oisy Samples Alessandro Nordio Member IEEE , Carla-Fabiana Chiasserini Member IEEE , Emanuele V iterbo Member IEEE Abstract W e focus on a multidimen sional field with uncorrelated spectrum, and study the quality of the reconstruc ted signal when the field sample s a re ir regularly spac ed a nd af fected by indep endent and identically distrib u ted noise. More specifically , we ap ply linear reconstru ction techn iques a nd take the mean s quare error (MSE) of the field estimate as a m etric to ev alu ate the signal reco nstruction quality . W e fin d th at the MSE a nalysis co uld be car ried ou t by using the closed-f orm expression o f th e eigenvalue distribution of the ma trix repr esenting the sam pling system. Unfortu nately , such distribution is still unknown. Thus, we first derive a closed-f orm expr ession o f the distribution moments, and we find that the eigenv alue distribution tends to the Mar ˇ cenko-Pastur distribution a s th e field dimension goes to infinity . Finally , b y u sing ou r ap proach , we derive a tight approx imation to the MSE of the reco nstructed field. I . I N T R O D U C T I O N W e address the important issue of recons tructing a m ultidimensional signal from a collection of sa mples that are noisy and not uniformly spa ced. As a case study , we c onsider a wireless sens or network for en vironmen tal mo nitoring, where the nodes sensing the physica l phenome non (hereina fter a lso called field) are randomly deployed over the area und er o bservation. The sens ors samp le a d -dimensional spatially finite p hysical field, where d ma y take into ac count spatial dimensions as well as the temporal dimension. Examp les of su ch fields are pressure or temperature, on a 4-dimens ion domain, i.e., three spatial coordinates plus the time dimension. A spatially fi nite p hysical field is not bandlimited, howe ver it admits a n infinite Fourier series expansion. Here, we co nsider a finite a pproximation of the physical fie ld obtained by tr u ncating su ch series, ass uming that the contrib ution of the truncated terms is negligible. This work was supported by MIUR through the MEADOW project Nov ember 6, 2021 DRAFT 2 In our c ase s tudy , we assume tha t the meas ured sa mples are transferred from the sensors to a common data-collecting unit, the so-called sink nod e, which is in cha rge of reconstructing the field. Distributed systems, where in-network processing, is performed are out of the sc ope of this work. W e do not dea l with issues related to information transport and, thus, we assume that all samples are co rrectly receiv ed at the sink node. The field samp les, howev e r , a re corrupted by additiv e i.i .d. noise, due to quantization, round-off errors or quality o f the sens ing device. Furthermore, the samp ling points are known at the sink node, becaus e (i) either sensors a re located at pre -defined p ositions or their position ca n b e estimated through a localization tec hnique [1], a nd (ii) the sa mpling time is either p eriodic or inc luded in the information sent to the sink. Several efficient and fast algo rithms have bee n proposed to numerically recons truct or approximate a signal in such setting, which amount to the solution of a linear syste m (see [2], [3] a nd reference s therein). A widely used technique consists in process ing the se nsors’ me asures by means of a linear filter , which is a fun ction of the system parameters known at the sink . W e observe tha t the following two major factors affect the linear reconstruction: (i) the given machine precision, which may prevent the rec onstruction algorithm from performing correctly an d may lead t o a non-negligible probab ility of recon struction failure [3], (ii) the noise level affecting the sensors’ measureme nts. In the latter cas e, a measure o f reconstruction acc uracy is gi ven by the mean square error (MSE) of the fie ld estimate. In [4], [5], we h ave found that thes e issues could be studied by us ing the eigen value distrib u tion of the reconstruction matrix; howe ver , obtaining such a distrib ution is still an open problem. In this work, we first extend the system model and the problem formulation presented in [5] to the ca se of multidimensional fields (Section III). Then, we deri ve a close d-form expres sion of the moments of the eigenv alue distributi o n, through as ymptotic analysis (Section V). By using the moments expression s, we prove that the eigen value distrib ution of the matrix repres enting the sa mpling system tends to the Mar ˇ cenko-Pastur d istrib ution [6] as the fi eld dimension d → ∞ (Sec tion VI). W e apply our results to the study of the MSE of the field estimate, when the sens ors measurements are noisy a nd the re construction at the sink is performed throug h linear filtering. W e ge neralize the MSE expressions to the multidimensional ca se (with fin ite d ), a nd we show tha t, by using the Mar ˇ ce nko-Pastur distributi on ins tead of the a ctual eigen value d istrib u tion, we obtain an approximation to the MSE of the reco nstructed fi eld which is very tight for d ≥ 3 (Section VII). Before providing a de tailed desc ription of o ur ana lysis, in the next section we discuss some related studies an d h ighlight our main contributi o ns with res pect to p revious work. Nov ember 6, 2021 DRAFT 3 I I . R E L A T E D W O R K A N D M A I N C O N T R IB U T I O N S Relev an t to our work is the literature on spectral analysis, where, however , se veral studies deal with regularly sampled signals (e.g ., [7] and references therein). An excellent guide to irregular sampling is [8], which p resents a la r g e number of tec hniques , algorithms, and a pplications. Reconstruction tec hniques for irregularly or rando mly samp led signals can be foun d in [3], [9], [10], just to n ame few . In particular , Feichtinger a nd Gr ¨ ochenig in [10] provide a n error analysis of an iterati ve rec onstruction algorithm taking into accou nt rou nd-off errors, jitters, truncation errors and aliasing. From the theoretical point of view , irregular sampling has b een studied in [3], [9]–[14] and referenc es therein. In the con text of se nsor ne tworks, efficient techniques for spatial sampling are propos ed in [15], [16]. In p articular , in [16], an a daptiv e sa mpling is described, w hich a llo ws the central data-co llector to v ary the number of ac ti ve sens ors, i.e., s amples, acc ording to the desired reso lution le vel. Data a cquisition is also studied in [17], where the authors conside r a unidimension al field, un iformly sampled at the Nyq uist frequency by low precision s ensors. The a uthors sh ow that the numbe r of samples can be trade d-of f with the p recision of sensors. Th e problem of the recons truction of a band limited signa l from an irregular se t of samp les at unknown locations is ad dressed in [18]. T here, different solution me thods are proposed, and the conditions for which there exist multiple solutions or a uniqu e solution are discus sed. Diff erently from [18], we assu me that the sink can either acq uire or estimate the sensor loca tions and that sensors are randomly de ployed. The field recons truction at the sink no de with spatial and temporal correlation among sen sor mea sures is s tudied, for instance, in [19]–[23]. Other interesting studies c an be found in [24], [25], wh ich address the perturbations of regular sampling in shift-in variant spac es [24] and the rec onstruction of irregularly sampled images in presence of measu rement noise [25]. W e point out that our ma in c ontributi o n with resp ect to p revious work on signal sampling and reconstruction is the probabilistic approac h we ad opt to analyze the quality lev e l o f a signal reconstructed from a set of irregular , noisy samples. Our analysis, howev e r , ap plies to sampling sys tems where the field recons truction is pe rformed in a ce ntralized manner . Finally , we high light that our previous work [5] assume s that senso rs are uniformly distrib u ted over the s patial observation interval a nd may be d isplaced around a known average loca tion. The e f fe cts of noisy measures an d jittered positions are analyze d when linear reconstruction techniques are employed. Howev e r , o nly the unidimensional case is studied and semi-analytical deri vations o f the MSE of the reconstructed fie ld a re obtained. In [26], instead, sensors are a ssumed to be fixed , and the objectiv e is to evaluate the p erformance of a linear recon struction Nov ember 6, 2021 DRAFT 4 technique in the presence of qua si-equally spaced se nsor layouts. A. Main r esults The goal of this work is to provide an analytical s tudy on the reconstruction qu ality of a multidi- mensional ph ysical field, with uncorrelated spectrum. Th e field samples a re (i) irregularly spaced, s ince they are gathe red by a randomly deployed sen sor network and (ii) affected b y i.i.d. noise. The sink node receives the field samples and runs the reconstruction a lgorithm in a ce ntralized ma nner . Our major contributi ons with res pect to p revious work are a s follows. 1. Given a d -dimensional prob lem formulation, w e obtain an alytical expressions for the mome nts o f the eigen value d istrib u tion of the rec onstruction matrix. Using the exp ressions o f the moments, we show that the eigen value d istrib ution ten ds to the Ma r ˇ cenko-Pastur distributi on [6] as the field dimension d → ∞ . 2. W e a pply our resu lts to the study of the quality of a rec onstructed field an d d eri ve a tight approxi- mation to the MSE of the estimated field. I I I . P R E L I M I N A R I E S W e first presen t the multidimensiona l formulation of o ur rec onstruction problem. Then, we give some backgrou nd on linear reconstruction techniques and generalize to the multidimensiona l case some results previously obtained in the unidimens ional ca se [5]. Finally , we highlight the main steps followed in our study . Notation: L ower case bo ld letters denote column vec tors, while up per c ase bo ld letters denote matrices. W e denote the ( h, k ) -th entry of the matrix X by ( X ) h,k , the transpos e o f x by x T , and the conjugate transpose of x by x † . The identity matrix is deno ted by I . Finally , E [ x ] is the average of x an d subs cripts to the a verage operator s pecify the variable with respe ct to which the a verage is tak en. A. Irr e g ular sampling of mu ltidimensional signals Let us consider a d -dimensional, s patially-finite physical field ( d ≥ 1 ), where r s ensors are located in the hype rcube H = { x | x ∈ [0 , 1) d } and measure the value of the field. W e ass ume that the senso r sampling points are known. At first, we c onsider that they a re deterministic, the n we will assume that they are i.i.d. rand om variables uniformly distrib uted in the h ypercube H . Nov ember 6, 2021 DRAFT 5 When obs erved over a finite region, a d -dimens ional phys ical field s ( x ) with finite energy E s admits an infi nite d -dimens ional Fourier s eries expa nsion with coefficients ˜ a ℓ , s uch tha t E s = P + ∞ ℓ 1 ,...,ℓ d = −∞ | ˜ a ℓ | 2 where ℓ = [ ℓ 1 , . . . , ℓ d ] is a vector o f integers a nd ℓ m , m = 1 , . . . , d represen ts the ind ex of the expansion along the m -th d imension. W e truncate the expansion to 2 M + 1 terms per dimension w here M is such that P − M − 1 ℓ 1 ,...,ℓ d = −∞ | ˜ a ℓ | 2 + P + ∞ ℓ 1 ,...,ℓ d = M +1 | ˜ a ℓ | 2 ≪ E s Therefore, one can think of M as the approximate one-sided bandwidth (per dimension) of the field, which can be approximated over the finite region H as s ( x ) = (2 M + 1) − d/ 2 X ℓ a ν ( ℓ ) e j2 π x T ℓ (1) where the term (2 M + 1) − d/ 2 is a normalization factor and P ℓ represents a d -dimensional sum over the vector ℓ , with ℓ m = − M , . . . , M . Also, a ν ( ℓ ) = ˜ a ℓ and the func tion ν ( ℓ ) = d X m =1 (2 M + 1) m − 1 ℓ m , − (2 M +1) d − 1 2 ≤ ν ( ℓ ) ≤ + (2 M +1) d − 1 2 maps the vector ℓ on to a scalar index. Note tha t, wh ile ˜ a ℓ has a vectorial index, a ν ( ℓ ) has a scalar index and it has been introduced to simplif y the notation. As an example, for d = 2 and M = 1 , we have ν ( ℓ ) = 3 ℓ 1 + ℓ 2 and s ( x 1 , x 2 ) = 1 3 P 1 ℓ 1 = − 1 P 1 ℓ 2 = − 1 a 3 ℓ 2 + ℓ 1 e j2 π ( x 1 ℓ 1 + x 2 ℓ 2 ) . Let X = { x 1 , . . . , x r } , wit h x q = [ x q , 1 , . . . , x q ,d ] T ∈ H , q = 1 , . . . , r , be the s et of sampling points, and s = [ s 1 , . . . , s r ] T , s q = s ( x q ) , the values o f the correspond ing field sa mples. Follo wing [3], we write the vector of field v a lues s as a function of the sp ectrum: s = G † d a (2) where a is a vector of size (2 M + 1) d , w hose ν ( ℓ ) -th entry is given by a ν ( ℓ ) , an d G d is the (2 M + 1) d × r matrix: ( G d ) ν ( ℓ ) ,q = (2 M + 1) − d/ 2 e − j2 π x T q ℓ (3) In general, the entries of a can be correlated with cov a riance matrix E [ aa † ] = σ 2 a C a , and T r { C a } = (2 M + 1) d . In the following, we restrict our a nalysis to the class o f fields c haracterized by E [ aa † ] = σ 2 a I . If the sensor me asurements , p = [ p 1 , . . . , p r ] T , are no isy , then the relation between sensors’ samples a nd field spec trum can be written as: p = s + n = G † d a + n (4) where the noise is represented by t he r -size , zero-mean random vector n , with cov ariance matrix E [ nn † ] = σ 2 n I r . W e de fine the signal-to-noise ratio on the measu re as: SNR m △ = σ 2 a /σ 2 n △ = 1 /α . Nov ember 6, 2021 DRAFT 6 B. Sampling rate Follo wing [5], we introduce the p arameter β define d as : β = (2 M + 1) d r (5) This paramete r represents the ratio between the n umber of harmonics used for the field reconstruction and the number o f sensors sampling the field. In the following, we consider β ∈ [0 , 1) . Notice that for fixed β and M , the numbe r r of samples exponen tially increas es with d . C. P r eviou s results on r eco nstruction qua lity Gi ven an e stimate ˆ a of the field spec trum a , the reco nstructed signal is: ˆ s ( x ) = (2 M + 1) − d/ 2 X ℓ ˆ a ν ( ℓ ) e j2 π x T ℓ (6) As reco nstruction performance metric, we conside r the MSE of the field es timate, which, for any given set of sa mpling points X , is defin ed as: MSE X = E a , n Z H | ˆ s ( x ) − s ( x ) | 2 d x = E a , n k ˆ a − a k 2 (2 M + 1) d (7) where the av erage is taken with respe ct to the subsc ripted random vectors. Note that (7) still assumes that the sampling po ints are de terministic; this as sumption will be remov e d later in the pape r . For li n ear models s uch as (4), several estimation techniques based on linear filtering are available in the literature [27]. W e emp loy a filter B such that the estimate of the field spec trum is giv en by the linear ope ration ˆ a = B † p (8) where B is an r × (2 M + 1) d matrix. In particular , we c onsider the linear filter providing the best performance in terms of MSE, i.e., the linear minimum MSE (LMMSE) filter 1 [27]: B = G † d ( R d + α I ) − 1 (9) where R d = G d G † d . From now on, we carry out our analysis u nder the assump tion that the elements of the set X are independ ent rando m vectors, w ith i.i.d. entries, uniformly distrib uted in the hypercu be H . In [5], we have shown that a simple and ef fecti ve tool to ev a luate the p erformance of lar g e finite systems is asymptotic an alysis. W e co mputed the MSE by letting the fi eld numb er of harmonics and the 1 Notice that when the cov ari ance matrix of a is known, the LMMSE fil ter generalizes to ( G † d C a G d + α I ) − 1 G † d C a . Nov ember 6, 2021 DRAFT 7 number of samples grow to infinity , while their ratio β = (2 M + 1) d /r is kept constant. W e observed the validity o f a symptotic ana lysis res ults, even for sma ll values of M and r . Similarly , here we consider as performance metric the asymptotic ave rage MS E, normalized to σ 2 a : MSE ∞ = lim M ,r → + ∞ β 1 σ 2 a E X [MSE X ] (10) where β belo w the limit d enotes the ratio which is kept cons tant. In (10), the av e rage is over all pos sible realizations of the set X . Using (7) –(9) and the above definitions, in Appendix E we s how that, in the case of the LMMSE: MSE ∞ = E λ d,β  αβ λ d,β + αβ  (11) where λ d,β is a random variable with p robability de nsity function (pdf) f d,β ( x ) , distributed as the asymptotic eigen values of T d = β R d = β G d G † d . The subscripts d and β o f λ indicate that the distrib u tion of the asymptotic eige n values of T d depend s on b oth the fie ld dimension d and the pa rameter β . The matrix T d plays an important role in our an alysis; in the follo wing , we introduce some of its properties. In the unidimen sional c ase ( d = 1 ), T 1 is a (2 M + 1) × (2 M + 1) He rmitian T oeplitz matrix giv e n b y T 1 = T † 1 =         t 0 t 1 · · · t 2 M t − 1 t 0 · · · t 2 M − 1 . . . t − 2 M t − 2 M +1 · · · t 0         where ( T 1 ) ℓ,ℓ ′ = t ℓ − ℓ ′ = 1 r P r q =1 exp( − j2 π ( ℓ − ℓ ′ ) x q ) , ℓ, ℓ ′ = − M , . . . , M . For d ≥ 2 , T d can be defined recursively as a (2 M + 1) d × (2 M + 1) d Hermitian Block T oeplitz matrix with non Hermitian T oeplitz blocks: T d =         B 0 B 1 · · · B 2 M B − 1 B 0 · · · B 2 M − 1 . . . . . . . . . B − 2 M B − 2 M +1 · · · B 0         where ( T d ) ν ( ℓ ) ,ν ( ℓ ′ ) = 1 r r X q =1 e − j2 π ( ℓ − ℓ ′ ) x q (12) and ℓ , ℓ ′ ∈ [ − M , . . . , M ] d . T hat is, the matrix T d is compo sed of (2 M + 1) 2 blocks B i of s ize (2 M + 1) d − 1 × (2 M + 1) d − 1 , each i ncluding (2 M + 1) 2 blocks of s ize (2 M + 1) d − 2 × (2 M + 1) d − 2 , Nov ember 6, 2021 DRAFT 8 and s o on. The smallest block s have s ize ( 2 M + 1) × (2 M + 1) ; they hav e the s ame s tructure as matrix T 1 in the unidimen sional case, howe ver only thos e on the main diag onal have a He rmitian structure. Proof of this is giv en in [28] for d = 2 ; the exten sion to the d - dimensional ca se is straightforward. I V . E S T I M A T IO N E R RO R C A L C U L AT I O N M E T H O D The ana lysis detailed in the n ext s ections cons ists o f the f ollo wing main s teps. (i) As a practical cas e, we consider the asy mptotic expre ssion of the LMMSE in (11) and notice that an analytical ev a luation of the asymptotic L MMSE co uld be obtained by exploiting the closed- form express ion of the eigen value distrib ution, f d,β ( x ) , of the reconstruction matrix. However , suc h expression is s till u nknown. Hence, as a first s tep we d eri ve a clos ed form expression of the moments of λ d,β , for a ny d and β , an d provide an algo rithm to comp ute them. (ii) W e s how that the v alue of the moments of the eigen value distrib ution d ecrease s as the field dimension d inc reases. (iii) W e prove that, as d → ∞ , the expres sion of the eigen value distributi on tends to the Mar ˇ c enko-Pastur distrib u tion. (iv) By u sing the Mar ˇ cenko-Pastur distrib ution, we are able to obtain a tight approximation for the LMMSE of the reconstructed field, wh ich holds for a ny finite value of d . V . C L O S E D F O R M E X P R E S S I O N O F T H E M O M E N T S O F T H E A S Y M P T OT I C E I G E N V A L U E P D F Ideally , we would like to ob tain the a nalytical express ion of the distrib u tion f d,β ( x ) of the asymptotic eigen value of T d , for a gi ven β . Unfortunately , such a ca lculation see ms to b e prohibitive and is still an open problem. Therefore, as a fi rst step, we comp ute the closed form expres sion of the mome nts E [ λ p d,β ] of λ d,β , for a ny positi ve integer p . In the limit for M and r growing to infi nity with constant β , the expression of E [ λ p d,β ] can be easily obtained from the powers of T d as in [29 ], [30], E [ λ p d,β ] = lim M ,r → + ∞ β 1 (2 M + 1) d T r E X  T p d  (13) In Section V -A, we show tha t E [ λ p d,β ] is a polynomial in β , of degree p − 1 (s ee (24)); the remaining subsec tions de scribe how to compute this polyno mial. Nov ember 6, 2021 DRAFT 9 A. P artitions Using (12), the term T r E X [ T p d ] in (13) can be written as T r E X  T p d  = E X   X ℓ 1 ( T p d ) ν ( ℓ 1 ) ,ν ( ℓ 1 )   = E X " X L ∈L d ( T d ) ν ( ℓ 1 ) ,ν ( ℓ 2 ) · · · ( T d ) ν ( ℓ p ) ,ν ( ℓ 1 ) # = 1 r p X q ∈Q X L ∈L d E X h e − j2 π P p i =1 x T q i ( ℓ i − ℓ [ i +1] ) i (14) where Q = { q | q = [ q 1 , . . . , q p ] } , q i = 1 , . . . , r L d = { L | L = [ ℓ 1 , . . . , ℓ p ] } , ℓ i = [ ℓ i 1 , · · · , ℓ i d ] T , ℓ i m = − M , . . . M and [ i + 1] =    i + 1 1 ≤ i < p 1 i = p In (14), the a verage is performed over the random set of positions X = { x 1 , . . . , x r } , with independent and uniformly distrib u ted elements. T o o btain a closed-form expression of the distribution mome nts, we rewr ite (14) by using set partitioning. Let P = { 1 , . . . , p } be the set of integers from 1 to p . W e observe t hat any giv e n vector q ∈ Q partitions the set P into 1 ≤ k ( q ) ≤ p disjoint non-empty sub sets P 1 ( q ) , . . . , P k ( q ) , where P j , j = 1 , . . . , k ( q ) , is the se t of ind ices of the entries of q ta king the s ame value γ j . Tha t is, P j ( q ) = { i ∈ P | q i = γ j } (15) and k ( q ) is the numbe r of distinct values γ j taken by the entries of vector q . Subs ets P j have the follo wing properties k ( q ) ∪ j =1 P j ( q ) = P , P j ( q ) ∩ P j ′ ( q ) = ∅ for , j 6 = j ′ . Also, we p oint out that, since r is the numbe r of values tha t the e ntries q i can take, there exist r ! / ( r − k ( q ))! vectors q ∈ Q generating a given partition of P mad e of k ( q ) s ubsets. In order to clarify the abov e conc epts, we provide an examp le b elow . Nov ember 6, 2021 DRAFT 10 Example 1: Let p = 6 , then P = { 1 , 2 , 3 , 4 , 5 , 6 } . Also, let q = [4 , 9 , 5 , 5 , 4 , 3] . Since the distinct values in q are γ j = 3 , 4 , 5 , 9 , we have k ( q ) = 4 . It follo ws that P is p artitioned into the following s ubsets: P 1 ( q ) = { 1 , 5 } ( q 1 = q 5 = 4) P 2 ( q ) = { 2 } ( q 2 = 9) P 3 ( q ) = { 3 , 4 } ( q 3 = q 4 = 5) P 4 ( q ) = { 6 } ( q 6 = 3) Hence, the partit ion of P induced by q is {{ 1 , 5 } , { 2 } , { 3 , 4 } , { 6 }} . Next, we introduce an ef fectiv e me thod to rep resent partitions of a set P , by building a tree o f depth p , as in Figure 1. Such a repres entation will allo w to simplify the notation in the f ollo wing analysis . T o buil d the tree of d epth p , we proceed a s follows. Ea ch node of the tree is assigned with a label from the s et P = { 1 , . . . , p } , s tarting from the root which is labe led by 1. Each node generates m + 1 leaves, labeled in increasing order from 1 to m + 1 , where m is the largest lab el on the pa th fr o m the root to such node . Note that, a t level p , any value in { 1 , . . . , p } is used to lab el the leaves a t least on ce. Then, gi ven a tree of depth p , we d efine ω = [ ω 1 , . . . , ω p ] as a path of length p from the tree’ s root to a leaf. W e obs erve that a vector q can b e repres ented as a p ath ω in the tree of d epth p . This is don e by assigning a label ( j = 1 , . . . , p ) in increasing order to every distinct value of q ; the v ector c ollecting the labels is the pa th ω co rresponding to the giv e n q . W e h ave that the path ω , correspo nding to a giv e n q , defi nes in the tree of depth p the same p artition of the se t P as the on e induced by q . Indeed, giv e n a pa rtition of P , the subse t P j defined in (15) c an be rewritten as P j ( ω ) = { i ∈ P | ω i = j } , (16) i.e., as the s et of integers correspond ing to the depths of the j - th label in the p ath. As a las t remark, consider the number of d istinct values in q (i.e., the numbe r of distinct labels in ω ) to be equ al to k ( q ) , and recall that the number of all possible values taken by the p elements o f q is equ al t o r . It follows that r ! / ( r − k ( q ))! dif ferent q ’ s yield the same vector ω . This is in agreemen t with the fact that, giv e n P , there are r ! / ( r − k ( q ))! dif feren t q ’ s generating the sa me partition consisting of k ( q ) sub sets. Aga in, for the s ake of clarity , we gi ve an example. Nov ember 6, 2021 DRAFT 11 p = 4 p = 3 p = 2 p = 1 1 2 1 2 3 1 2 3 3 2 4 3 2 2 1 2 1 1 2 3 1 1 1 Fig. 1. Partitions tree of depth p = 4 . The path ω = [1 , 2 , 1 , 1] employ ed in Example 2 is highlighted by using dashed lines Example 2: Let us consider p = 4 and q = [4 , 9 , 4 , 4] . The vec tor q c an be represen ted in the tr ee of d epth 4 as the path ω = [1 , 2 , 1 , 1] (i.e., the path highlighted with dash ed lines in Figure 1 ). In q there are two distinct values (name ly , 4 and 9), or , equiv alently , in the pa th ω there a re two labels (namely , 1 and 2); then k ( q ) = k ( ω ) = 2 . La bel 1 ap pears in ω at depths 1,3, and 4 ( P 1 = { 1 , 3 , 4 } ), while label 2 appears at de pth 2 ( P 2 = { 2 } ) . The partition of P = { 1 , 2 , 3 , 4 } induced by q or , equiv alen tly , by ω is therefore: {{ 1 , 3 , 4 } , { 2 }} . From the d iscuss ion a bove, it s hould be clea r that con sidering a partition of P is e quiv alent to considering a path ω in a tree of depth p . Hen ce, in the following an alysis, we will refer to a partition through its c orresponding path ω . W e now exploit set partitioning to rewri te (14). Since the random vectors x q are indepe ndent, gi ven Nov ember 6, 2021 DRAFT 12 q , the a verage operator in (14 ) factorizes into k ( q ) terms, i.e ., E X h e − j2 π P p i =1 x T q i ( ℓ i − ℓ [ i +1] ) i = k ( q ) Y j =1 E x γ j h e − j2 π x T γ j P i ∈P j ( q ) ℓ i − ℓ [ i +1] i (17) Each term de pends on a single random vector x γ j . Moreover , since the entries of x γ j are indepe ndent random variables u niformly distributed in [0 , 1) , we have: E x γ j h e − j2 π x T γ j P i ∈P j ( q ) ℓ i − ℓ [ i +1] i = d Y m =1 E x γ j ,m h e − j2 π x γ j ,m c jm i = d Y m =1 δ ( c j m ) (18) where x γ j ,m and ℓ i,m are the m -th entries of x γ j and ℓ i , resp ectiv ely , where the function δ ( · ) is the Kronecker’ s delta, and wh ere c j m = P i ∈P j ( q ) ℓ i,m − ℓ [ i +1] ,m . By substituting (17) and (18) in (14) and by expanding the s ummation P L ∈L d , we ob tain T r E X [ T p d ] = 1 r p X q ∈Q X ℓ 1 ∈L 1 · · · X ℓ d ∈L 1 k ( q ) Y j =1 d Y m =1 δ ( c j m ) = 1 r p X q ∈Q d Y m =1   X ℓ m ∈L 1 k ( q ) Y j =1 δ ( c j m )   = 1 r p X q ∈Q   X ℓ ∈L 1 k ( q ) Y j =1 δ ( c j )   d (19) where c j = P i ∈P j ( q ) ℓ i − ℓ [ i +1] . For any giv en q ∈ Q , the expression ζ M ( q ) = X ℓ ∈L 1 k ( q ) Y j =1 δ ( c j ) (20) is a polyn omial in M , sinc e it represents the number of points with integer co ordinates con tained in the hypercub e [ − M , . . . , M ] p and satisfying the k ( q ) constraints: c j = X i ∈P j ( q ) ℓ i − ℓ [ i +1] = 0 (21) for j = 1 , . . . , k ( q ) . In Appen dix A, we sho w that one of these constraints is alw a ys redundant and that the number o f linearly independ ent constraints is exactly equa l to k ( q ) − 1 . As a c onsequ ence, the polyno mial ζ M ( q ) ha s degree p − k ( q ) + 1 , and, for large values of M , we hav e ζ M ( q ) = O (( 2 M + 1) p − k ( q )+1 ) . Nov ember 6, 2021 DRAFT 13 Now , using (19) an d (20), the limit in (13) is gi ven by E [ λ p d,β ] = lim M ,r → + ∞ β X q ∈Q ζ M ( q ) d r p (2 M + 1) d (22) Equation (22) c an be further simplified by conside ring that there exist r ! / ( r − k ( q ))! vectors q ∈ Q generating a given partition o f P ma de of k ( q ) subse ts, or , equiv alently , a path ω of length p with k ( q ) distinct labels. Let Ω p be the s et of vectors ω , eac h c orresponding to a distinct partition of P . Also, let us write k ( q ) and ζ M ( q ) as functions of the partition indu ced by q , i.e., as functions o f ω . T hen, we obtain: E [ λ p d,β ] = lim M ,r → + ∞ β X ω ∈ Ω p X q ⇒ ω ζ M ( ω ) d r p (2 M + 1) d ( a ) = lim M ,r → + ∞ β X ω ∈ Ω p ζ M ( ω ) d r ! ( r − k ( ω ))! r p (2 M + 1) d (23) where • the nota tion P q ⇒ ω represents the sum over all vectors q ge nerating a given p ath ω , • the eq uality ( a ) ho lds beca use the number of vectors q ge nerating a given ω is r ! / ( r − k ( ω ) )! . Note that, for lar ge r , r ! / ( r − k ( ω ))! = r k ( ω ) + O ( r k ( ω ) − 1 ) . Also , since ζ M ( ω ) is a polynomial in M of degree p − k ( ω ) + 1 , for large values of M we have: ζ M ( ω ) = v ( ω )(2 M + 1) p − k ( ω )+1 + O ((2 M + 1) p − k ( ω ) ) , wh ere v ( ω ) is the coefficient of d egree p − k ( ω ) + 1 o f ζ M ( ω ) . Therefore, taking the limit, we obtain: E [ λ p d,β ] = X ω ∈ Ω p v ( ω ) d β p − k ( ω ) = p X k =1 β p − k X ω ∈ Ω p,k v ( ω ) d (24) where Ω p,k ⊆ Ω p is the su bset of Ω p containing pa ths with k ( ω ) distinct labels, a nd v ( ω ) = lim M → + ∞ ζ M ( ω ) (2 M + 1) p − k ( ω )+1 (25) Note that the coefficient v ( ω ) represen ts the volume of the con vex po lytope desc ribed by the cons traints in (21), wh en the v ariables ℓ i are c onsidered as rea l a nd limit ed to a p -dimens ional hy percube of volume 1 . As a co nsequ ence, we h ave: 0 ≤ v ( ω ) ≤ 1 . Equation (24) provides a c losed-form expression of the moment E [ λ p d,β ] , a s a polynomial in β of degree p − 1 . Again, for the sake o f clarity , we giv e an e xample below . Nov ember 6, 2021 DRAFT 14 Example 3: Let p = 6 and q = [4 , 9 , 5 , 5 , 4 , 3] . W e hav e ω = [1 , 2 , 3 , 3 , 1 , 4] , and the partition of P is {{ 1 , 5 } , { 2 } , { 3 , 4 } , { 6 }} . Then, the set o f k ( ω ) = 4 cons traints (as in (21)) a re g i ven by: ℓ 1 + ℓ 5 = ℓ 2 + ℓ 6 , ℓ 2 = ℓ 3 , ℓ 3 + ℓ 4 = ℓ 4 + ℓ 5 , ℓ 6 = ℓ 1 The last equation is redundant since can be obtained from the first three c onstraints. Simplifying, we obtain ℓ 1 = ℓ 6 , and ℓ 2 = ℓ 3 = ℓ 5 . Since each variable ℓ i ranges from − M to M , the nu mber of integer so lutions satisfying the cons traints is exac tly ζ M ( ω ) = (2 M + 1) 3 , and then v ( ω ) = lim M → + ∞ ζ M ( ω ) (2 M +1) p − k ( ω )+1 = 1 . Next, in order to compute E [ λ p d,β ] , we need: • to enume rate the partitions, i.e ., the vec tors ω ∈ Ω p,k , for e ach k = 1 , . . . , p (see Se ction V -B) ; • to compu te the coefficients v ( ω ) , for any ω ∈ Ω p,k and k = 1 , . . . , p (see Section V -C). B. P artitions enumeration W e notice tha t Ω p represents the set of partitions of a p -element s et, thus it has ca rdinality | Ω p | = B ( p ) , where B ( p ) is the p -th Bell number or exponen tial number [31]. Furthermore, the subset Ω p,k ⊆ Ω p has cardinality S ( p, k ) , w hich is a Stirling number of the s econd kind [32] gi ven by: S ( p, k ) = 1 k ! k X i =0 ( − 1) i  k i  ( k − i ) p with B ( p ) = P p k =1 S ( p, k ) . C. Co mputation of the coefficients v ( ω ) The last step required for the computation of E [ λ p d,β ] is the ev aluation of the c oefficients v ( ω ) , for ev ery ω ∈ Ω p . W e have the follo wing Lemma : Lemma 5.1: For any ω ∈ Ω p (or , equiv a lently , any partit ion of P ) a nd any arbitrary integer n , with n = 1 , . . . , k ( ω ) , the coefficient v ( ω ) in (25) is given by: v ( ω ) = Z R k ( ω ) − 1 p Y i =1 sinc  y ω i − y ω [ i +1]  | y n =0 d y n (26) where Nov ember 6, 2021 DRAFT 15 • y n = [ y 1 , . . . , y n − 1 , y n +1 , . . . , y k ( ω ) ] T , • sinc ( y ) = sin( π y ) π y . Pr oof: The p roof can be found in Appe ndix B. D. A practical method for the mom ents compu tation Equation (26) in Section V -C sho ws that the computation of E [ λ p d,β ] requires the evaluation of B ( p ) integrals. Howe ver , B ( p ) is very large even for small p , e.g. , B (10 ) = 115975 and B (20) ≈ 5 · 10 13 . The co mputational co mplexity c an be reduced b y recursively ap plying the simplification rules define d in the foll owi ng Lemma: Lemma 5.2: Let • ω = [ ω 1 , . . . , ω p ] be the p ath in a tree of dep th p , c orresponding to the pa rtition of P into k subs ets; • P 1 , . . . , P k be the subsets of P = { 1 , . . . , p } de fined a s in (16); • i ∈ P j : • ω ′ be the path o btained from ω by removing ω i . W e have the following rules: 1) if P j has cardinality 1 (i.e., P j is a sing leton ) or 2) if P j contains adjace ncies (in the circula r sens e), i.e ., both i and [ i + 1] ∈ P j , then v ( ω ) = v ( ω ′ ) Pr oof: The p roof is a direct consequen ce of Le mma 5.1 an d c an be fou nd in App endix C. T able I shows two examples of how the rules desc ribed in Le mma 5.2 can b e app lied. Examp le 1 in the T able assu mes p = 6 an d ω = [1 , 2 , 3 , 2 , 2 , 1] . At step 1, we no te that the third elemen t ( i = 3 ) o f ω is a singleton, then, by applying rule 1, we can remove it from the path. At step 2, we find that in ω the re are so me adjace ncies, hence we apply twice rule 2 (steps 2 and 3). At step 4, the second eleme nt of ω is a singleton, and we remove it by ap plying rule 1. Eventua lly , at ste p 7, the path ω is empty (i.e., has size p = 0 ) and , thu s, the corresponding co efficient is v ( ω ) = 1 ( E [ λ 0 d,β ] = E [1] = 1) . Example 2 in the T able a ssumes p = 6 and ω = [1 , 2 , 3 , 1 , 2 , 1] . After removing a singleton (step 1) and an ad jacency (step 2 ), the remaining pa th cannot be further reduce d. Then, to compu te the coefficient v ( ω ) , we need to apply directly Lemma 5.1 o n the p ath ω = [1 , 2 , 1 , 2] . W e obtain: v ( ω ) = Z + ∞ −∞ sinc( y 1 − y 2 ) sinc( y 2 − y 1 ) · sinc( y 1 − y 2 ) sinc( y 2 − y 1 ) | y 2 =0 d y 1 = Z + ∞ −∞ sinc( y 1 ) 4 d y 1 = 2 3 Nov ember 6, 2021 DRAFT 16 T ABLE I E XA M P L E O F C O M P L E X I T Y R E D U C T I O N U S I N G T H E R U L E S D E S C R I B E D I N L E M M A 5 . 2 Example 1 Example 2 Step ω Rule i ω Rule i 1 [1,2,3,2,2,1] 1 3 [1,2,3,1,2,1] 1 3 2 [1,2,2,2,1] 2 2 [1,2,1,2,1] 2 5 3 [1,2,2,1] 2 2 [1,2,1,2] 4 [1,2,1] 1 2 5 [1,1] 2 2 6 [1] 1 1 7 [] v ( ω ) = 1 v ( ω ) = 2 / 3 In the foll owi ng examp le, L emmas 5.1 and 5.2 a re exploited to explicitly co mpute E [ λ 4 d,β ] . Example 4: Le t us c onsider p = 4 . The total nu mber of pa rtitions of P = { 1 , 2 , 3 , 4 } is equal to B (4) = 15 . Conside ring the tree of de pth p , we apply to eac h pa th the rules of Lemma 5.2, and we find that 1 4 pa ths (pa rtitions) out of 1 5 redu ce to the empty path, thus contributi n g with v ( ω ) = 1 . The only path that cannot b e further reduced is ω = [1 , 2 , 1 , 2] . Thus , applying Lemma 5 .1 with n = 2 , we ob tain v ( ω ) = 2 / 3 . From (24 ) a nd c onsidering all contributions, we obtain: E [ λ 4 d,β ] = β 3 +  6 + (2 / 3) d  β 2 + 6 β + 1 . V I . C O N V E R G E N C E T O T H E M A R ˇ C E N KO -P A S T U R D I S T R I B U T I O N In Section V, we have shown that the moments o f the a symptotic eigen values o f T d are polynomials in β , gi ven by (24). In pa rticular , the p -th mome nt E [ λ p d,β ] has degree p − 1 an d is g i ven by the s um of B ( p ) positi ve contrib u tions of the f orm v ( ω ) d β p − k ( ω ) . Since 0 < v ( ω ) ≤ 1 an d β > 0 , for any d , the follo wing inequality holds: E [ λ p d +1 ,β ] ≤ E [ λ p d,β ] i.e., for any gi ven p and β , the moments of the asymptotic eigen values decrease as the field dimension increases . The series E [ λ p d,β ] , as a function of d , is positiv e and monotonically decreasing, thus it con verges Nov ember 6, 2021 DRAFT 17 to: E [ λ p ∞ ,β ] = lim d → + ∞ E [ λ p d,β ] (27) Lemma 6.1: The mome nts E [ λ p ∞ ,β ] a re the Na rayana polynomials, given by E [ λ p ∞ ,β ] = p X k =1 T ( p, k ) β p − k (28) where T ( p, k ) = 1 k  p − 1 k − 1  p k − 1  are the Naraya na numbers [34], [35]. Mo reover , the random variable λ p ∞ ,β follo ws the Mar ˇ cenko-Pastur d istrib ution [6] with pdf (see Figure 2): f ∞ ,β ( x ) = p ( c 1 − x )( x − c 2 ) 2 π xβ (29) where c 1 , c 2 = (1 ± √ β ) 2 , 0 < β ≤ 1 , c 2 ≤ x ≤ c 1 . Pr oof: The p roof is gi ven in Appen dix D. 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 0.5 1 1.5 2 2.5 3 3.5 4 f ∞,β (x) x β =0.2 β =0.4 β =0.6 β =0.8 β =1.0 Fig. 2. Mar ˇ cenko -Pastur distribution In the following, we ap ply our finding s to the stud y of the L MMSE of a rec onstructed multidimensional field; in particular , we exploit the Mar ˇ ce nko-Pastur distribution to compu te the expectation in (11). Nov ember 6, 2021 DRAFT 18 1e-6 1e-5 1e-4 1e-3 1e-2 1e-1 1 50 40 30 20 10 0 -10 -20 MSE SNR m β =0.1, MSE ∞ β =0.1, MP β =0.2, MSE ∞ β =0.2, MP β =0.4, MSE ∞ β =0.4, MP β =0.8, MSE ∞ β =0.8, MP Fig. 3. MS E of the reconstructed fi eld for d = 1 and va rying values of β . Comparison between the MSE asymptotic value (11) and the fully analytical expression deriv ed using the Mar ˇ cenko-Pastur distribution (30) V I I . S T U DY O F T H E R E C O N S T RU C T I O N Q UA L IT Y T H R O U G H T H E M A R ˇ C E N KO -P A S T U R D IS T R I B U T I O N Recall tha t the MSE pro vided by the L MMSE filter is [5]: MSE LMMSE ∞ = E λ d,β  αβ λ d,β + αβ  where λ d,β is distributed as the a symptotic eigenv alues of T d , with p df f d,β ( x ) . By us ing the Mar ˇ cen ko-Pastur distributi o n f ∞ ,β instead of f d,β , we ha ve: MSE ∞ = E λ ∞ ,β  αβ λ ∞ ,β + αβ  = Z c 1 c 2 αβ f ∞ ,β ( x ) x + αβ d x = 2 β − θ + p θ 2 − 4 β 2 β (30) where θ = 1 + β (1 + α ) . Equation (30) provides a n ap proximation to the MSE ∞ , which, as sho w n in the following p lots, ca n be exploited to d eri ve the quality o f the re constructed field, g i ven a finite d . W e first conside r d = 1 and compare in Figure 3 the expression of the MSE ∞ as in (11) (solid lines) with the one obtained by using the Mar ˇ cen ko-Pastur distrib ution (da shed lines). The results are Nov ember 6, 2021 DRAFT 19 1e-5 1e-4 1e-3 1e-2 1e-1 1 50 40 30 20 10 0 -10 -20 MSE SNR m d=1 d=2 d=3 MP Fig. 4. MSE of the reconstructed field, for β = 0 . 4 and d = 1 , 2 , 3 . Comparison between the MSE asymptotic v alue (11) and the fully analytical expression deriv ed using the Mar ˇ cenko-P astur distribution (30) presented as functions of the SNR m and for diff erent values of β . W e compu ted (11) by averaging over the e igen values of 200 rea lizations of the matrix T 1 , with M = 150 . The plot shows that, for small values of β , the Mar ˇ c enko-Pastur distrib ution (30) yields an excellent approximation to the MS E ∞ , already for d = 1 . Instead, for values of β g reater than 0 .2, the expres sion in (30) fails to provide a vali d approximation. Howe ver , it is interesting to notice that, for d > 1 , it is possible to obtain an ac curate approximation o f the MSE ∞ using the Mar ˇ cenko-Pastur distribution, even for lar ge values of β . This is shown by Figures 4 and 5, which plot the results ob tained throu gh (11) and (30) for β e qual to 0.4 and 0.8, respectiv ely . The results are presented as the SNR m varies and for dif ferent v alue s of the field dimension d . Looking a t Figure 4, we note that our approximation is tight for d ≥ 2 , wh ile Figure 5 shows that, when d = 3 , we s till get a fairly go od app roximation for β as lar ge as 0 . 8 . Nov ember 6, 2021 DRAFT 20 1e-5 1e-4 1e-3 1e-2 1e-1 1 50 40 30 20 10 0 -10 -20 MSE SNR m d=1 d=2 d=3 MP Fig. 5. MSE of the reconstructed field, for β = 0 . 8 and d = 1 , 2 , 3 . Comparison between the MSE asymptotic v alue (11) and the fully analytical expression deriv ed using the Mar ˇ cenko-P astur distribution (30) V I I I . C O N C L U S I O N S W e cons idered a large-scale wireless senso r network sampling a multidimensional field, and we in vestigated the mean square e rror (M SE) of the signal reconstructed a t the sink node. W e noticed that an analytical study of the quality of the reconstructed field could be carried out by using the eigen value distrib u tion of the matrix representing the samp ling sy stem. Sinc e s uch a distrib ution is un known, we first deri ved a c losed-form expression of the distrib ution momen ts. By us ing this expression , we were able to s how that the eigen value d istrib ution of the reconstruction matrix tends to the Mar ˇ cenko-Pastur distrib u tion as the fie ld dimension ten ds to infinity . W e applied our results to the study of the MSE of the reconstructed field, wh en linea r fi ltering is us ed at the sink no de. W e found that, b y using the Mar ˇ ce nko- Pastur dis trib ution ins tead of the a ctual eigen value distributi on, we obtain a close a pproximation to the MSE of the reco nstructed signal, which holds for field dimen sions d ≥ 2 . W e believe tha t our work is the basis for an analytical study of various aspec ts c oncerning the reconstruction qua lity of multidimension al sen sor fields, and , more generally , of irregularly sampled Nov ember 6, 2021 DRAFT 21 signals. A P P E N D I X A T H E C O N S T R A I N T S Let us consider a vector of integers q of size p partiti o ning the set P = { 1 , . . . , p } in k subse ts P j , 1 ≤ j ≤ k , and the set of k c onstraints (21). W e first show tha t one of these cons traints is a lw a ys redundan t. A. Redundant constraint Choose an integer j , 1 ≤ j ≤ k . Summing up a ll constraints exce pt for the j - th, we get: k X h =1 h 6 = j c h = k X h =1 c h − c j = X i ∈P ℓ i − ℓ [ i +1] − c j = − c j which giv es the j -th constraint since P i ∈P ℓ i − ℓ [ i +1] = 0 . Thus, one o f the constraints in (21 ) is always redundan t. Next, we s how that the remaining k − 1 constraints are linearly indep enden t. B. Linear indepe ndence The k c onstraints in (21) can be arranged in the form: W ℓ = 0 with ℓ = [ ℓ 1 , . . . , ℓ p ] T and W being a k × p matri x defi ned as W = W ′ − W ′′ (31) where ( W ′ ) j,i =    +1 i ∈ P j 0 otherw ise and W ′′ is obtaine d from W ′ by c ircularly s hifting the rows to the right by one po sition. Since one of the co nstraints (21) is redundant, the rank of W is: ρ ( W ) ≤ k − 1 . Now we prove that the rank of W is equ al to k − 1 . Since the subsets P j have empty interse ction, the rows of W ′ are linearly independen t; hence, W ′ has rank k . Also, W ′′ is obtaine d from W ′ by circularly s hifting the rows by one position to the right, thus W ′′ can be written as W ′′ = W ′ S where S is the p × p right-shift matrix [36], i.e. , the entries of the i -th row of S a re ze ros except for an entry equal to 1 at pos ition [ i + 1] . As a conseque nce, W = W ′ − W ′ S = W ′ ( I p − S ) w here the ro ws of the matrix I p − S are obtained b y circularly shifting Nov ember 6, 2021 DRAFT 22 the vector [+1 , − 1 , 0 , . . . , 0] and thus has rank ρ ( I p − S ) = p − 1 . Hence, using the prop erties of the rank of ma trix produc ts reported in [36], we have ρ ( W ) = ρ ( W ′ ( I p − S )) ≥ ρ ( W ′ ) + ρ ( I p − S ) − p = k − 1 W e reca ll that the system of linear equa tions Wz ha s a finite numbe r of integer solutions bounded in [ − M , . . . , M ] d . The number o f solutions dec reases as ρ ( W ) increase s. A P P E N D I X B P R O O F O F L E M M A 5 . 1 Pr oof: Using (20 ) and (25), we obtain: v ( ω ) = lim M → + ∞ 1 (2 M + 1) p − k ( ω )+1 X ℓ ∈L 1 k ( ω ) Y j =1 δ ( c j ) W e first notice tha t Q k ( ω ) j =1 δ ( c j ) = δ ( W ℓ ) where the k ( ω ) × p matrix W is de fined in Appendix A and δ ( W ℓ ) is a multidimensional Kronecker de lta [33]. Sinc e the rank of W is ρ ( W ) = k ( ω ) − 1 (see Appendix A), the n δ ( W ℓ ) define s a subspa ce of Z p with p − k ( ω ) + 1 dimens ions. The refore, co nsidering that ℓ is a v ector of inte g ers with e ntries ranging in the interval [ − M , . . . , M ] a nd taking the limi t for M → ∞ , we obtain v ( ω ) = R [ − 1 / 2 , 1 / 2) p δ d ( Wz ) d z where z ∈ R p and the fun ction δ d represents the Dirac delta. W e hav e that δ d ( Wz ) can be f actorized as δ d ( Wz ) = k ( ω ) Y j =1 δ d ( r T j z ) (32) where r T j is the j -th row of W . As a lready sh own in Appe ndix A, o ne o f the constraints (21 ) is redun dant and, hence, one of the factors in the right hand side of (32), say the n -th, must not be included in the product. Now , moving to the Fourier transform domain, we c an write: δ d ( r T j z ) = R + ∞ −∞ exp(j2 π y r T j z ) d y . Therefore, v ( ω ) = Z [ − 1 / 2 , 1 / 2) p k ( ω ) Y j =1 ,j 6 = n δ d ( r T j z ) d z = Z [ − 1 / 2 , 1 / 2) p k ( ω ) Y j =1 ,j 6 = n Z + ∞ −∞ e j2 π y j r T j z d y j d z = Z [ − 1 / 2 , 1 / 2) p Z R k − 1 e j2 π P k ( ω ) j =1 ,j 6 = n y j r T j z d y n d z Nov ember 6, 2021 DRAFT 23 where y n = [ y 1 , . . . , y n − 1 , y n +1 , . . . , y k ( ω ) ] T . Integrating first with respect to z , we ge t v ( ω ) = Z R k − 1 Z [ − 1 / 2 , 1 / 2) p e j2 π P p i =1 z i y T n w i d z d y n = Z R k − 1 Z [ − 1 / 2 , 1 / 2) p p Y i =1 e j2 π z i y T n w i d z d y n = Z R k − 1 p Y i =1 Z 1 / 2 − 1 / 2 e j2 π z i y T n w i d z i d y n = Z R k − 1 p Y i =1 e j π y T n w i sinc ( y T n w i ) d y n = Z R k − 1 e j π y T n P p i =1 w i p Y i =1 sinc ( y T n w i ) d y n where w i is the i -th column of W , taken after removing its n -th row . By defin ition, the j -t h rows of W ′ and of W ′′ contain both p − |P j | “0 ” and |P j | “+ 1”. Since W = W ′ − W ′′ , we h ave P p i =1 w i = 0 and v ( ω ) = R R k − 1 Q p i =1 sinc ( y T n w i ) d y n . Notice that, by d efinition of W (see (31)), y T n w i = y j − y j ′ | y n =0 if i ∈ P j and [ i + 1] ∈ P j ′ . Moreover , by the definition in (16), we ha ve y j = y ω i when i ∈ P j . Thus , v ( ω ) = Z R k − 1 p Y i =1 sinc ( y ω i − y ω [ i +1] ) | y n =0 d y n A P P E N D I X C P R O O F O F L E M M A 5 . 2 ( S I M P L I F I C A T IO N R U L E S ) Let P j be a singleton with P j = i and ω i = j . W e first notice that, since P j is a singleton, [ i − 1] , [ i +1] / ∈ P j . By a pplying Lemma 5 .1 with an arbit rary n 6 = j , we have v ( ω ) = Z R k − 1 p Y h =1 sinc ( y ω h − y ω [ h +1] ) | y n =0 d y n = Z R k − 1 p Y h =1 h 6 = i h 6 =[ i − 1] sinc ( y ω h − y ω [ h +1] ) sinc ( y ω [ i − 1] − y ω i ) · sinc ( y ω i − y ω [ i +1] ) | y n =0 d y n Nov ember 6, 2021 DRAFT 24 W e now integrate with respect to y j , with j = ω i and we obtain v ( ω ) = Z R k − 2 p Y h =1 h 6 = i h 6 =[ i − 1] sinc ( y ω h − y ω [ h +1] ) · Z R sinc ( y ω [ i − 1] − y j ) sinc ( y j − y ω [ i +1] ) | y n =0 d y j d y ′ n = Z R k − 2 p Y h =1 h 6 = i h 6 =[ i − 1] sinc ( y ω h − y ω [ h +1] ) · sinc ( y ω [ i − 1] − y ω [ i +1] ) | y n =0 d y ′ n = Z R k − 2 p − 1 Y h =1 sinc ( y ω ′ h − y ω ′ [ h +1] ) | y n =0 d y ′ n = v ( ω ′ ) where y ′ n and ω ′ have been o btained from y n and ω by removing their j -th and i -th element, respec ti vely . Obviously y ′ n has size k − 1 and ω ′ has size p − 1 . Let P j be such that: P j = i, [ i + 1] , i.e., ω i = ω [ i +1] = j . Then, v ( ω ) = Z R k − 1 p Y h =1 sinc ( y ω h − y ω [ h +1] ) | y n =0 d y n = Z R k − 1 p Y h =1 h 6 = i sinc ( y ω h − y ω [ h +1] ) · sinc ( y ω i − y ω [ i +1] ) | y n =0 d y n = Z R k − 1 p Y h =1 h 6 = i sinc ( y ω h − y ω [ h +1] ) sinc ( y j − y j ) | y n =0 d y n = Z R k − 1 p Y h =1 h 6 = i sinc ( y ω h − y ω [ h +1] ) | y n =0 d y n = Z R k − 1 p − 1 Y h =1 sinc ( y ω ′ h − y ω ′ [ h +1] ) | y n =0 d y n = v ( ω ′ ) where ω ′ has bee n obtained from ω by removing its i -th element. A P P E N D I X D P R O O F O F L E M M A 6 . 1 In order to prove Lemma 6.1, we first note that Ω p,k may contain both c rossing and n on-crossing partitions [37]. Nov ember 6, 2021 DRAFT 25 a) Non -cr oss ing partitions: Every non-crossing partition contains at lea st a singleton or a subse t with a djacenc ies, and therefore can b e redu ced by using the rules in L emma 5.2. After reduc tion, the resulting partition is still no n-crossing, thus it can be further reduced till the empty se t is rea ched. It follo ws that the non-cross ing pa rtition ω ∈ Ω p,k contributes to the expression of E [ λ p p,k ] with a coe fficient v ( ω ) = 1 . b) Crossing p artitions: Reca ll that, in ge neral, the coefficient v ( ω ) defined in (26) ca n be obtained by c ounting the solutions of the sy stem of equations: W ℓ = 0 where the k ( ω ) × p matrix W contains the co efficients of the k ( ω ) c onstraints in (21 ). If ω ∈ Ω p,k is a crossing p artition, then • k ( ω ) ≥ 2 (by d efinition, a partit ion with k ( ω ) = 1 is alw ay s non-crossing) • it con tains at least tw o subsets P j and P j ′ , with j 6 = j ′ , which a re cros sing. Some c rossing partit ions can be reduc ed by applying the rules in Le mma 5.2 but, even after reduction, they remain crossing. Let us now focu s on the crossing subse t P j of a partition ω which has bee n reduced by applying the rules in Lemma 5.2. W ithout loss o f generality , we a ssume that |P j | = h , i.e. the partition P j contains h elements with h ≥ 2 sinc e P j is no t a singleton. T hen, b y defin ition of the matrix W (see Appendix A) its j -t h ro w , r T j , contains h entries with value 1 , h en tries with value − 1 a nd p − 2 h ze ros. W e then buil d the 2 × p matrix f W as f W = [ r j , − r j ] T . Notice that f W has rank 1 and the sys tem of equations f W ℓ = 0 contains the con straints induced by a partition e ω = [1 , 2 , . . . , 1 , 2] with 2 h en tries. Since the system of equations f W ℓ = 0 co ntains a r educed set of cons traints with resp ect to W ℓ = 0 and, thus, a lar g er numb er of solutions , it follows that v ( ω ) ≤ v ( e ω ) . It is straightforward to s how that for a p artition such as e ω , with k ( e ω ) = 2 , the coefficient v ( e ω ) is giv e n b y Lemma 5.1 as v ( e ω ) = R R sinc ( y ) 2 h d y . T his is a decreasing function of h and since h ≥ 2 we have: v ( e ω ) = Z R sinc ( y ) 2 h d y ≤ Z R sinc ( y ) 4 d y = 2 3 Therefore, we c onclude that v ( ω ) ≤ v ( e ω ) < 1 . Nov ember 6, 2021 DRAFT 26 c) Crossing and non-cr os sing partitions: Le t Ω c p,k , Ω n p,k ⊂ Ω p,k be, respectively , the set of crossing and non -crossing partitions of Ω p,k , with Ω c p,k ∩ Ω n p,k = ∅ and Ω c p,k ∪ Ω n p,k = Ω p,k . The n, E [ λ p ∞ ,β ] = lim d → + ∞ p X k =1 β p − k X ω ∈ Ω p,k v ( ω ) d = lim d → + ∞ p X k =1 β p − k   X ω ∈ Ω c p,k v ( ω ) d + X ω ∈ Ω n p,k v ( ω ) d   ( a ) = p X k =1 β p − k | Ω n p,k | where the equality ( a ) is due to the fact that for non-cros sing partitions v ( ω ) = 1 , while for crossing partitions v ( ω ) < 1 an d, h ence, lim d → + ∞ v ( ω ) d = 0 . In [38 ], it can b e found tha t the number of non-crossing partit ion s of size k in a p -eleme nt se t is gi ven by the Narayan a numbers T ( p, k ) = | Ω n p,k | and therefore E [ λ p ∞ ,β ] = P p k =1 T ( p, k ) β p − k are the Narayana polynomials. In [6], it is shown that the Narayana polyn omials a re the mo ments of the Mar ˇ cenko-Pastur distrib ution. A P P E N D I X E P R O O F O F ( 1 1 ) W e show tha t wh en the LMMSE filter is us ed, the expres sion of the asy mptotic MSE is gi ven by (11). Indeed, by using (10), (7), (8), and (9) we hav e : MSE ∞ = lim M ,r → + ∞ β 1 σ 2 a (2 M + 1) d E X E a , n  k A − 1 d G d p − a k 2  where A d = R d + α I a nd R d = G d G † d . Substituting (4) in the above expression and ass uming E [ aa † ] = σ 2 a I and E [ nn † ] = σ 2 n I , we get 1 σ 2 a E a , n  k A − 1 d G d p − a k 2  = T r n ( A − 1 d R d − I )( A − 1 R d − I ) † + α A − 1 d R d A − 1 d o = T r  α ( R d + α I ) − 1  = T r  αβ ( T d + αβ I ) − 1  where T d = β R d . L et us consider a n analytic function g ( · ) in R + . L et X = UΛ U † be a random p ositi ve definite Hermitian n × n matrix, whe re U is the eigen vectors matrix of X and Λ is a diagonal matri x containing the eigen values of X . By using the res ult for symmetric ma trices in [39, Ch. 6] c ombined with the result in [40, pag. 481], we ha ve: lim n →∞ 1 n T r E [ g ( X )] = E λ [ g ( λ )] where the random v ariable Nov ember 6, 2021 DRAFT 27 λ is distri buted a s the as ymptotic eigen values of X . It follo ws that lim M ,r → + ∞ β T r E X  αβ ( T d + αβ I ) − 1  (2 M + 1) d = E λ β ,d  αβ λ β ,d + αβ  R E F E R E N C E S [1] D. Moore, J. Leonard, D. Rus, and S. T eller , “Robust distributed network localization with noisy r ange measurements, ” ACM Confer ence on Embedded Networked Sensor Systems (SenSys 2004) , B altimore, MD, pp. 50-61, Nov . 2004. [2] H. G. Feichtinger and K. Grchenig, “Theory and practice of irregular sampling, ” in W avelets: Mathematics and Applications, J. J. Benedetto and M. W . F razier , Eds. Boca Raton, FL, CRC P ress, pp. 305–363 , 1994. [3] H. G. Feichtinger , K. G r ¨ ochenig, and T . Strohmer , “Efficient numerical methods in non-uniform sampling theory , ” Numerische Mathematik , V ol. 69, pp. 423–440, 1995. [4] C.-F . Chiasserini, A. Nordio, and E . V iterbo, “On data acquisition and field reconstruction in wireless sensor networks, ” in Signals and Communication T echn ology , Distributed Cooper ative Laboratories: Networking, Instrumentation, and Measur ements , Eds. F . Davoli, S . Palazzo and S. Zappatore, Eds. Springer Berlin Heidelberg, pp. 161–171 , 2006. [5] A. Nordio, C.-F . Chiasserini, and E. V iterbo “Quality of field reconstruction in sensor networks, ” IEEE INFOCOM Mini- Symposium , Anchorage, AK, May 2007. [6] V . A. Mar ˇ cenko and L . A. Pastur , “Distribution of eigen v alues for some sets of random matrices, ” USSR Sbornik , V ol. 1, pp. 457–483, 1967. [7] P . Stoica and R. Moses, Intr oduction to Spectra l Analysis, Upper Saddle Riv er, NJ, Prentice-Hall , 2000. [8] F . A. Marvasti, “Nonuniform Sampling: Theory and Practice, ”, New Y ork: Kluwer Academic: Plenum, 2001. [9] H. Rauh ut, “Random sampling of sparse trigonometric p olynomials, ” online: http://home page.uni vie.ac.at/holger .rauhut/RandomSampling.pdf . [10] H.G. Feichtinger and K. G r ¨ ochenig, “Error analysis i n regular and irregular sampling theory , ” A pplicable Analysis, V ol. 50, pp. 167–189, 1993. [11] H. Rauhut, “Stability results for random sampling of sparse trigonometric polynomials, ” online: http://homepage.uni vie.ac.at/holger .rauhut/StabilityRS.pdf . [12] E. Candes and T . T ao, “The Dantzig selector: St atistical estimation when p i s much larger t han n , ” online: http://www .acm.caltech.edu/ ∼ emmanuel/pape rs/DantzigSelector .pdf . [13] A. Aldroubi and K. Gr ¨ o chenig, “Nonuniform sampling and reconstruction in shift-in variant spaces, ” SIAM Rev iew , Dec. 2001, V ol. 43, pp. 585-6 20. [14] K. Gr ¨ ochenig, “Irregular sampling, T oeplitz matrices, and the approximation of entire functions of exponential type, ” Mathematics of Computation, V ol. 68, No. 226, pp. 749–765, Apr . 1999. [15] M. Perill o, Z. Ignjatov ic, and W . Heinzelman, “ An energy conserv ati on method for wireless sensor networks employing a blue noise spatial sampling technique, ” IPSN, Apr . 2004. [16] R. W i llett, A. Martin, and R. Now ak, “Backcasting: adaptiv e sampling for sensor networks, ” IPSN, Apr . 2004. [17] P . Ishwar , A. Kumar , and K. Ramchandran, “Distributed sampling for dense sensor networks: a bit-conserv ation principle, ” IPSN, Apr . 2003. [18] P . Marziliano and M. V etterli, “Reconstruction of irre gularly sampled discrete-time bandlimited signals with unkno wn sampling locations, ” IEEE T ransactions on Signal Pr ocessing , V ol. 48, No. 12, pp. 3462–3471 , Dec. 2000. Nov ember 6, 2021 DRAFT 28 [19] R. Cristescu and M. V etterli, “On the optimal density for real-time data gathering of spatio-temporal processes i n sensor networks, ” IPSN, Los Angeles, CA, Apr . 2005. [20] Y . Sung, L. T ong, and H. V . Poor, “Sensor activ ati on and scheduling for field detection in large sensor arrays, ” IPSN, Los Angeles, CA, Apr . 2005. [21] D. Ganesan, S. Ratnasamy , H. W ang, and D . Estrin, “Coping wit h irregular spatio-temporal sampling in sensor networks, ” ACM SIGCOMM, pp. 125–130, Jan. 2004. [22] M. C. V uran, ¨ O. B. Akan, and I. F . Akyildiz “Spatio-temporal correlation: theory and applications f or wi reless sensor networks, ” Computer Networks, V ol 45, No. 3, pp. 245–259, June 2004. [23] Y . Rachlin, R. Negi, and P . Khosla, “Sensing capacity for discrete sensor network applications, ” IPSN, Los Angeles, CA, Apr . 2005. [24] P . Zhao, C. Zhao, and P . G. Casazza, “Perturbation of regular sampling in shift-inv ariant spaces for frames, ” IEEE T ransactions on Information Theory , V ol. 52, No. 10, pp. 4643–4648, Oct. 2006. [25] D. S. E arly and D. G. L ong, “Image reconstruction and enhance d resolution imaging from i rregular samples, ” IEEE T ransactions on Geoscience and Remote Sensing, V ol. 39, No. 2, pp. 291–302, Feb. 2001. [26] A. Nordio, C.-F . Chiasserini, and E. V it erbo, “The impact of quasi-equally spaced sensor layouts on field reconstruction, ” IPSN, Cambridge, MA, Apr . 2007. [27] S. M. Kay , Fundamentals of statistical signal pro cessing: Estimation theory , Prentice Hall, 1993. [28] T . Strohmer , “Computationally attractiv e reconstruction of bandlimited images from irregular samples, ” IEEE T ransactions on Image Processing , V ol. 6, No. 4, pp. 540–5 45, Apr . 1997. [29] P . Bill ingsley , Pr obability and measur e ( 3r d edition) , John W iley and Sons I nc, Ne w Y ork, 1995. [30] L. Li , A. M. Tulino and S. V erd ` u, “ Asymptotic eigen value moments for linear multiuser detection, ” Communications in Information and Systems, V ol. 1, No. 3, pp. 273–304, S ept. 2001. [31] E. W . W eisstein, “Bell number , ” from MathW orld – A W olfram W eb Resour ce , http://mathworld.wo l fram.com/BellNumber .html. [32] E. W . W eisstein, “S tirling number of the second kind, ” from MathW orld – A W olfram W eb Resour ce , http://mathworld.wo l fram.com/StirlingNumberoftheSecondKind .html . [33] Kr oneck er delta, http://en.wikipedia.org/wiki/Kroneck er delta. [34] “The on-line ency clopedia of integer sequences, ” http://www .research.att.com/ ∼ njas/sequences/A0 01263 . [35] I. Dumitriu and E. Rassart, “Path counting and random matrix theory , ” The E lectr onic Journ al of Combinatorics, V ol. 10, No. 1, 2003. [36] M. Brook es, T he matrix r efere nce manual, http:// www . ee.ic.ac.uk/hp/staf f /dmb/matrix/intro.html, 2005. [37] N. Dersho witz and S. Zaks, “Ordered trees and non-crossing partitions, ” Disc. Math, V ol. 62, pp.215–218, 1986. [38] R. P . Stanley , Enumerative Combinatorics, V ol. 2, Cambridge Univ ersity Press, Cambridge, 1999. [39] R. Bellman, Intr oduction to matrix analysis, 2nd Edition, McGraw-Hill, 1970. [40] W . H. Pr ess, S. A. T eukolsky , W . T . V etterling, and B. P . Flannery , Numerical rec ipes, 2nd Editi on, C ambridge University Press, 1997. Nov ember 6, 2021 DRAFT