extfieldsofstudy:
Computer Science
Mathematics
Physics

text:
Network depth: identifying median and contours in complex networks

Centrality descriptors are widely used to rank nodes according to specific concept(s) of importance. Despite the large number of centrality measures available nowadays, it is still poorly understood how to identify the node which can be considered as the `centre' of a complex network. In fact, this problem corresponds to finding the median of a complex network. The median is a non-parametric and robust estimator of the location parameter of a probability distribution. In this work, we present the most natural generalisation of the concept of median to the realm of complex networks, discussing its advantages for defining the centre of the system and percentiles around that centre. To this aim, we introduce a new statistical data depth and we apply it to networks embedded in a geometric space induced by different metrics. The application of our framework to empirical networks allows us to identify median nodes which are socially or biologically relevant.

Introduction
Nodes centrality is a very fundamental concept in network science, used to answer questions like "who are the most influential actors in a social network"? Or "which are the key nodes for the optimal functioning of a power grid"? Unfortunately centrality descriptors are defined upon a non-rigorous definition of importance. Consequently, there is a plethora of centrality measures and interpretations making difficult to understand which one could be used as a one-fits-all solution.
Similarly many statistical problems deal with the location, the centre of a population. Knowing very little about the underlying population, e.g. we may not know if the underlying distribution has finite first moment, we can still estimate the location through the median and set up distribution-free statistical tests for it. Nonparametric statistics largely rely on order statistic, which is fine in R, but gives some problems already in R 2 , where no natural order relation is defined. As we will see in Sec. 2, statistical data depths 2 of 26 G. BERTAGNOLLI, C. AGOSTINELLI AND M. DE DOMENICO enable the generalisation of order and order statistics, to multivariate spaces, functional spaces and, with this work, also to networks. Before moving to our new centrality, we have to summarise some key concepts about network centrality.
The need for a unified framework for centrality measures was already clear in the early days of network science, when Freeman proposed his three-fold categorisation of centralities [11]. More recently, Borgatti and Everett [4] approached centrality from a graph-theoretic perspective. They showed that all centralities evaluate the involvement of a node in the walkability of the network according to four features: Walk Type, Walk Property, Walk Position and Summary Type. Different Walk Types (e.g. shortest-paths or random walks) and restrictions on them (e.g. maximum length) might influence centrality as shown in Fig. 1. For instance the degree k i = ∑ j a i j of the i-th node can be seen as counting the paths of length 1 emanating from i. Allowing paths of length k leads to k-path centrality. Another example is given by the betweenness that, as the degree, belongs to the class of measures that count walks and have volume as the Walk Property, at variance with closeness that, instead, deals with paths length. The Walk Position can be radial, where walks emanate from a given node -as for degree and closeness -or medial, when the focus is on walks passing through that node (e.g. betweenness). Once the 3 of 26 information based on these features is gathered in a node-by-node matrix W ∈ R N×N , a summary statistic T : R N×N → R N is needed to obtain a centrality score vector C = T (W ). T represents the fourth feature: Summary Type. As observed in [4], the variation along this feature is small since usually the summary is a row sum T (W ) = W 1 or a (weighted) average. Table 1 places some popular centrality measures in this classification with respect to Walk Property and Position. We briefly recall the definitions of the five measures we use from now on, in bold font in Tab. 1, by using the same notation. In the following, let A = {a i j } N i, j=1 be the N × N adjacency matrix of graph G = (V, E). The degree of vertex i is C deg (i) = (W 1) i where W = A. Observing that in some occasions one's importance should depend on whether its neighbours are themselves important, many variations on the degree have been proposed. In 1972 Bonacich introduced the eigenvector centrality, which can be regarded as an elegant summary of Katzs, Hoedes and Hubbells measures [4]. Let x be the eigenvector corresponding to the largest eigenvalue of A, i.e. Ax = λ 1 x; then x is the eigenvector centrality vector and for each i ∈ 1, . . . , N ( 1.1) The entries of W are then w i j = a i j x j and C eig = T (W ) = λ −1 1 W 1. The closeness is obtained as the inverse of the marginals of W = D sp , the shortest-path distance matrix. The elements of its centrality score vector are computed as follows ∀i ∈ 1, · · · , N. (1.2) Observe that ∑ j =i d sp (i, j) is the total sp-distance from node i to all others. Sometimes the average distance [11] is considered to make the scores comparable among networks with different sizes. The betweenness centrality of node i, C btw (i), is the total number of shortest-paths between pairs of vertices that pass through node i.
Borgatti's Distance-weighted Fragmentation (DF) 1 [3] measures the extent to which the network cohesiveness depends on node i ∈ {1, . . . , N}. It depends both on the size of components of the network after the removal of i and on the relative cohesion of each component, in terms of total sp-distance between all pairs of nodes. Let d = max We have chosen this five measures, since they are representative of other centralities. The remainder of this paper is organised as follows. In Section 2 we present a new centrality built upon robust, nonparametric statistics, called statistical data depths. In particular, we start introducing the statistical background of data depths and the definition of our new multivariate depth function, the Projected Tukey depth (PTD). In Subsection 2.2 we make use of the graph-theoretical framework for centralities by Borgatti and Everett, to show that the PTD evaluated on a network embedded in space

Radial
Medial Volume Freeman degree, Bonacich eigenvector Betweenness Length closeness Borgatti DF Table 1: Table adapted from [4]. Walk Property and Walk Position dimensions and some popular centralities classified on the basis of this two features.
defines a new centrality in the sense of networks. We call this new centrality network depth. In addition, this framework allows differences and commonalities among the network depth and other centrality descriptors to be identified theoretically. In Subsection 2.3 we examine in detail the network embedding step in our procedure and study how it reflects on the network depth. In fact, networks can be embedded in different ways. We use Multidimensional scaling (MDS), a metric preserving embedding method that depends on two parameters: a distance (or dissimilarity) matrix and an embedding dimension. This step introduces a parameter space and, consequently, a family of network depths. We provide some applications to empirical networks in Section 3, together with the heuristics that should be followed to reproduce our results and to evaluate the network depth to real data. The conclusions in Section 4 summarise our contribution to the understanding of the concept of centre-outward ordering of nodes in a network.

Network Depth
A statistical data depth function, or simply a depth, can be seen as a measure of the centrality or outlyingness of a given sample with respect to its underlying probability distribution. In the analysis of multivariate data, it allows us to find an overall centre of the observations, which plays the role of the median. Further, depths induce a centre-outward and probability-based ordering of sample points, which results in the ranking of the observed units. Based on this ordering, quantitative and graphical univariate tools for data analysis have been generalised to the multivariate case [16]. The usefulness of depths is widely recognised, indeed they have been extended to functional data as well [19].
The idea underling data depths dates back to 1975 when John W. Tukey [24] proposed a new approach to the problem of ordering points in R p . Tukey's idea, which was further developed in [9], is based on the observation that, given a univariate random variable X with distribution function F X , a natural centre-outward ordering of points can be derived as follows: for each x ∈ R we can compute P(X x) = F X (x) and P(X x) = 1 − F X (x) and define D 1 (x, F X ) = min{P(X x), P(X x)}. The median of X -the value m such that P(X m), P(X m) 0.5 -is the most central (or the deepest) point, while points far away from it in both directions have smaller depth values, see Fig. 2. Then, the halfspace Tukey depth (HD) of a point x ∈ R p w.r.t. the distribution of a random vector X, is the least depth value of x in any one-dimensional projection of X.
Many depths and corresponding empirical estimators have been proposed since, but it was only in 2000 that statistical depth functions have been rigorously defined. In [27] Zuo and Serfling state four desirable properties a function should fulfil in order to be called a depth. These formalise the intuition that the concept of "centre" is relevant. Therefore, if a probability distribution is symmetric 2 around a point θ , that ought to be also the depth centre; further, the depth should take higher values near the FIG. 2: Centre-outward ordering of sample points. The multivariate median is the deepest point, while order statistics (quantiles) are generalised through depth-induced contours, which in the sample case can be estimated through convex hulls containing the most deep (central) (1 − α)100% sample points.
centre than towards the border. Properties (P2), (P3) and (P4) in the following definition are inspired precisely by this idea, while (P1) guarantees that common transformations (like standardisation) do not affect the depth of points. DEFINITION 2.1 (Depth function) Let x ∈ R p be a point and F X ∈ F , the probability distribution function of a random vector X on the Borel σ −algebra B(R p ). Furthermore, let the mapping D : R p × F → R be bounded, non-negative and satisfy P1) D(Ax + b, F AX+b ) = D(x, F X ) holds for any random vector X ∈ R p , any p × p non-singular matrix A and any vector b ∈ R p ; P2) D(θ θ θ , F X ) = sup x∈R p D(x, F X ) holds for any F X ∈ F having centre θ θ θ ; P3) for any F X ∈ F having deepest point θ θ θ , D(x, F X ) D(θ θ θ + α(x − θ θ θ ), F X ) holds for any 0 α 1; Then D(·, F X ) is called a statistical depth function.
A last key feature of depths, is the identification of a unique centre, regardless of the distribution being uni-or multi-modal. The reason why we emphasise it will be clear in Subsection 2.2.
Currently, there exist many multivariate depth functions, fulfilling all or just some properties of the definition 2.1, differing in robustness to outliers, ease of interpretation, and computational complexity among others. The halfspace Tukey depth is one of the most famous and well studied multivariate depth functions. Given that it is based on the "projection-pursuit" methodology, it behaves well overall compared with its competitors. Unfortunately, its exact computation is challenging (or prohibitive) in high-dimensional spaces. Indeed, given a point x ∈ R p and a p-variate random sample (X 1 , . . . , X n ) with empirical distributionF n , the (empirical) halfspace depth -HD n (x,F n ) = min n ∑ i=1 u t X i u t x : u = 1 -requires to consider projections along all possible u ∈ R p such that u = 1.
Our new depth is based on the original idea of Tukey and it addresses the main drawback of HD: the high computational cost.
2.1 Projected Tukey depth DEFINITION 2.2 (Projected Tukey depth) Let X be a random vector in R p with distribution F X , X 1 , X 2 being two independent copies and Σ the covariance matrix 3 of X (i.e. assume X ∈ L 2 (R p )). For a given x ∈ R p define the univariate random variable The Projected Tukey depth is a statistical depth function.
Proof. In the proof an equivalent 4 definition is also used: for each realisation Clearly PTD is bounded and non-negative. P1) Affine invariance: we have to prove that for any non singular matrix A ∈ R p×p and any vector b ∈ R p PTD(Ax + b, F AX+b ) = PTD(x, F X ).
= Y X ; the thesis follows since property P1 holds for HD. From now on we assume w.l.o.g. that Σ = Id.
P2) The proof of maximality at centre property involves probability distributions that are symmetric with around θ θ θ . There are several definitions of multivariate symmetry: spherical (most restrictive), elliptical, angular and halfspace (less restrictive) [26]; P2 can be proven in the latter, general case. Let the distribution of X be H-symmetric about a unique centre θ θ θ , i.e. in terms of the underlying probability measure P(H θ ) 1 2 for every closed halfspace H θ containing θ θ θ on its border. Assume θ θ θ = 0 thanks to P1. We ask if PTD(0, F X ) PTD(x, F X ) for all x ∈ R p . From (2.3) and [26,Th.2.3.2.] PTD(0, F X ) = inf HD 0, F Y 0 |X 2 1 2 3 Assume that Σ is invertible or use Moore-Penrose pseudo-inverse. 4 Equivalence can be easily proved writing the probability space (Ω , E , P) and X : Ω → R p with associated probability measure µ X (I) = P(X −1 (I)) for each Borel set I ⊂ R p . Then, ν = µ X 1 ⊗ µ X 2 is the product measure associated to (X 1 , X 2 ). In the following we write P(H) instead of µ X (H) for H ⊂ R p , with a negligible ambiguity, and F X (y) = P(X 1 y 1 , . . . , X p y p ) for y = (y 1 , . . . , y p ) ∈ R p . NETWORK DEPTH 7 of 26 Since the halfspace depth HD satisfies the maximality at centre property HD Π x,x 2 (0), F Y x |x 2 HD 0, F Y x |X 2 and HD 0, F Y x |X 2 1 2 for the uniqueness of the centre of symmetry assumption, it cannot happen that PTD(x, F X ) > 1 2 . P3) Monotonicity relative to deepest point: assuming w.l.o.g. θ θ θ = 0 we have to prove that PTD(x, F X ) PTD(αx, F X ) for α ∈ [0, 1]. The interesting case is when α ∈ (0, 1). The problem is that we have to compare two families of marginals Y (Ω )}. P3) has been proved only for spherical and elliptical symmetric distributions. In the first case we use a characterisation of spherical symmetry stating that rotations around the centre do not affect the distribution, formally X d = AX for orthogonal matrices A. Therefore, marginals on lines through the centre are all equal and the thesis follows immediately. Elliptic symmetry reduces to spherical when Σ = Id so by means of P1) we ho back to previous case.

P4) Vanishing at infinity -PTD
, 3π 2 then y x 0, otherwise y x 0. We want to show that for x sufficiently far away from the deepest point cos If X is bounded in R p , as soon as x ∈ R p is outside the ball of finite radius containing X(Ω ) the angle φ < π 2 and y x 0. If X is not bounded in R p , using Chebyshev's inequality on the univariate r.v. Y x|x 2 , i.e. ∀ε > 0 P Y x|x 2 ε 1 ε 2 we have that the probability of finding a point arbitrarily distant from E(X) is small. For x such that x → ∞ we can then conclude that P(φ < π 2 ) → 1 and the thesis follow. This closes the proof.
The sample version of the PTD is built around (2.3). It is similar to HD n with two important differences: first, only a subset of possible projecting directions are considered and exactly those indicated by that; secondly we do not take orthogonal but variance-scaled projections, akin in spirit to Mahalanobis distance. DEFINITION 2.4 (Sample PTD) Let (x 1 , . . . , x n ) be a random sample andF n its empirical distribution; ∀i, j = 1, · · · , n and x ∈ R p is the projection of the i-th sample point onto the line through x and x j . Collecting them for all j we get a univariate sample y · j = (y 1 j , · · · , y n j ) and we take the minimum probability mass carried by halflines with 0 on the boundary In the following we denote PT D n again by PT D.
In matrix form (2.5) becomesXΣ −1X t , whereX is obtained subtracting x from each row of X.
The set {PTD(x, F X ), ∀x ∈ X(Ω )} is called the depth space of X, or associated with distribution F X . In the empirical case the depth space is a vector whose entries are the depth value for each observation and this will be the centrality vector C PT D (G) of nodes in graph G. Another useful family of sets are the so-called depth trimmed (or central) regions R α , which generalise quantile-intervals, given α ∈ [0, 1] the depth region of order α is the set characterised by For instance, R α=0.99 is the depth central region containing the 1% deepest points. The set C α {x ∈ R p : PTD(x, P) = q α } is called depth contour of depth q α . The empirical version of the contour set is computed as the convex hull containing the most central fraction 1 − α of sample points and it is denoted by C α,n ; q α becomes q α,n and the depth region R α,n . Again, as we will always deal with empirical quantities, we drop the subscript n in the notation. The Projected Tukey depth can be computed exactly also in high-dimensional spaces. Given n points in R p , i.e. X ∈ R n×p , the algorithm of the PTD requires (i) the inversion of its sample covariance matrix Σ = cov(X), which costs O(p 3 ) in the worst case and (ii) the evaluation of the matrix productXΣ −1X t , which has complexity O(np 2 ) + O(pn 2 ) using the naive matrix multiplication algorithm. Usually p is not very high since we use dimensionality reduction techniques to embed our networks in space, so in the end the evaluation of the PTD of a point w.r.t. a cloud of n points has a complexity of O(pn 2 ) and of O(pn 3 ) for the whole data cloud.
In the next Subsection we show how this multivariate depth function can be used to generalise the definition of median to complex networks.

The Network Depth: a New Centrality
The key step linking the Projected Tukey depth to network centralities is network embedding. A network G = (V, E) embedded in some space R p is a cloud of points X(G) distributed accordingly to the connectivity of G. We map the problem of finding the centrality C(i, G) of the node i, to the problem of computing the depth PT D(x i ,F G ) of the point x i w.r.t. the (empirical) distributionF G of X(G). The idea is thatF G reflects the connectivity of graph G, adjacent vertices should lie nearer than nodes that do not share edges. Thus, the priority for our embedding is preserving distances.
On an undirected weighted connected network, we have a distance given by length of shortestpaths, D sp . If the graph is directed, D sp is no longer a metric and a non-metric embedding may be needed. Another distance can always be defined on (connected) networks, that is based on diffusive random walks, the diffusion distance D t [7]. Let p(t|i) be the vector with entries p k (t|i) representing the probability that a random walker starting at t = 0 in node i is found in k at time t. It is well known that this probability is linked to the network structure, in particular given the (random walk) normalised LaplacianL = Id − D −1 A, where D is the diagonal matrix of node degrees, the master equation of the "continuised" random walk isṗ (t) = −p(t)L.
The existence of a mapping from G to the geometric space induced by Markov dynamics, makes it possible to define the diffusion distance at time t, as shown in [5] and [7]  which expresses how nodes i and j are efficient in exchanging information. In other terms, D t (i, j; G) is small if there is a large transition probability from i to j along walks of length at most t. Observe that D t does not assume that the spreading takes place through shortest-paths and it is robust w.r.t. small fluctuations in the structure of the network. If the network is not connected, one has to regard each component separately or can consider a more general random walk dynamics with a teleportation term, as in diffusive processes used by PageRank [20] and Infomap [21]. Given a distance matrix describing structural and/or functional properties of our network we can recover (possibly with some approximations) the position of the N = |V | nodes in a space of dimension p ∈ N. p is always a free parameter of network embedding procedures. We call perfect embedding a configuration X ∈ R N×p of N points in R p for an appropriate p 1 such that x i − x j 2 = D(i, j) for all i, j = 1, . . . , N, where x i is the i-th row of matrix configuration X(G) and corresponds to the coordinates of node i in R p . For the moment, let us assume that, for a given distance matrix D, we have some p, such that the embedding of network G = (V, E) in R p is perfect (perfect embedding assumption). In statistical terms we do not lose information passing from G to X = X(G) and the centrality problem is equivalent to finding the depth space of X(G).
Let us denote by PTD(v; G, D sp , p) the PTD centrality of node v ∈ G with respect to geodesic embedding of G in R p and by PTD(G; D sp , p) the corresponding depth space, i.e. the network depth centrality vector with parameters D sp , p). We do not use the initial notation C PT D because the network depth, through the embedding step, depends on some free parameters. Similarly PTD(v; G, D t ,t, p) is the depth of node v in the diffusion embedding with parameters (t, p) and PTD(G; D t ,t, p) refers to the centrality vector ∀v ∈ G. If not ambiguous G is omitted. Using Borgatti's notation, the matrix W , containing the information on connectedness/social proximity, in our case is D or equivalently (under our perfect embedding assumption) X. Then the PTD is exactly the statistic summarising each node's share of total cohesion, as a matter of fact T (W ) = f X Σ −1X . Taking the Projected Tukey depth as a Summary Type makes the centrality, not only robust to outliers, but it also enables the extension to structured data of nonparametric statistical tools based on order statistics (e.g. location and scale estimators).
The other three features characterising centralities are Walk Type, Walk Property and Walk Position. Independently on the metric, PTD is a radial measure and, since it ignores multi-modality of distributions, it is not intended for the identification of local centres in a multi-community structure.
On the contrary Walk Type and Property depend on the distance. PTD(v; G, D sp , p), like Freeman's closeness and betweenness, takes into account only geodesics. In terms of walks, these measures assume the walker has a global knowledge of the network and it navigates it through the shortest possible paths 5 . Finally, PTD(v; G, D sp , p) is a length measure (Walk Property) since PTD(x, F X ) in R p can be interpreted as the "tailness" of point x with respect to the distribution F X . In other words, a point which lies deep inside the data cloud in the D sp embedding has small values of geodesic distance with all other points. On the other side, PTD(v; G, D t ,t, p) summarises structural information gathered through random walk exploration of the network. The diffusion time t sets a constraint on the length of otherwise unrestricted edges' sequences. As observed in [7], t plays the role of a scale parameter: for small values of t the micro-scale structure is revealed, see Fig. 4. In this case the network is embedded in a high dimensional space and the data cloud is very sparse with nodes with smaller D t lying deeper than all others, which are placed on a convex hull and have a depth of 1 N . When t increases diffusion distances display two different behaviours: D t among nodes exchanging efficiently information tends rapidly to zero, whereas D t (u, v) shrinks slowly if there is a small probability that two random walkers starting in u and v respectively meet somewhere in the network, following paths of length at most t. Since the Projected Tukey depth is not affected by the magnitude of the distances, but only on the relative position of points in space, hence after a certain amount of time, t > t * (G) the network depth spans a larger range and PTD(v; G, D t ,t, p) remains stable. Of course, in the limit for t → ∞, D t → 0 for all pairs of nodes, here small fluctuations due to non-perfect embedding might greatly affect the depth of points.
Finally, the Walk Property of the depth w.r.t. diffusion embedding is clearly volume, as a matter of fact i and j are near in the diffusion manifold if there is a high probability that two random walkers starting respectively in i and j meet somewhere in the network. In terms of connections, D t (i, j) is small if there are a lot of short paths connecting them.
The perfect embedding assumption is satisfied only if the original node distances D(i, j) = D E (i, j) are Euclidean distances in R p , then D E (i, j) = x i − x j and the configuration matrix X is reconstructed without errors. In general, x i − x j ≈ D(i, j) and matrix X is the result of a possibly non-metric Multidimensional Scaling (MDS). It is also not mandatory to have a distance matrix since MDS deals with general dissimilarity matrices ∆ = {δ i j } i, j . In the next subsection we will see the embedding step more in detail.

Network embedding
Multidimensional scaling is a collection of techniques aiming to find structures in multidimensional data, given a proximity matrix ∆ of observations and an embedding dimension p [6]. The MDS problem can be solved analytically or using iterative procedures, like the majorisation algorithm. The latter yields usually better results than the analytical solution, although its iterative nature may be time consuming for large N. When solving MDS by means of the majorization algorithm the loss function to be minimised FIG. 4: A Barabasi-Albert network with 100 nodes and its D t diffusion distance matrix for t = 3. The blocks in the distance matrix clearly reflect the groups in the network. Consider for instance v 1 , although it is linked to v 2 , which is a hub of the network, a random walker starting in v 1 will spend some time trapped in its branch of the network.
is Kruskal's stress-1: given a configuration X of N points in R p with 1 p N − 1 where δ i j are the original distances. If ∆ is not a proper distance matrix, δ i j are called dissimilarities and they may need to be transformed intod i j = f (δ i j ) before the scaling, for instance through a strictly increasing monotonic function f such that order is preserved δ i j < δ kl ⇒ f (δ i j ) < f (δ kl ). This case falls under the non-metric MDS methods. Furthermore, unfolding models allow us to deal with nonsymmetric matrices ∆ (which may be useful if in-and out-going geodesics have to be considered). The dimension p of the Euclidean space to which we map the network is then a free parameter of the network depth. If p is very small compared to N, the poor quality embedding will affect the depth too. Instead, if p is very large, the statistical information of the sample is so diluted in the high dimensional space, that all nodes happen to have the same depth of 1 N , i.e. they lie on a convex hull. Take for instance a star, a simple planar graph and the geodesic distance D sp . We may think that p = 2 should be the right embedding dimension but if we perform a metric MDS without any constraint we end up with a large stress value -since every permutation of nodes different from the star centre are equivalent -and a configuration X failing to match our expectations. A better approach is to use a constrained MDS on a 2D-sphere, maybe relaxing also the convergence tolerance.
In general, the most appropriate p is not known a-priori and we end up with an embedding for each p ∈ {1 p min , . . . , p max N −1} so that p is a free parameter for the depth. First of all, embeddings with "high" stress-1 values should be removed. What is "high" depends on the data. A typically approach is to look for the "legendary statistical elbow" in the stress curve (scree-plot) or to follow Kruskal's Rule of Thumb [14] stating that stress should be less than 0.2. Whether possible, it is warmly suggested -see [8] and [17] -to perform some additional goodness-of-fit analysis, like a permutation test to get a pvalue associated with the observed stress and diagnostic plots, e.g Stress-per-point plot to see if a subset of points accounts for large values of σ (X). As we already pointed out, for p → N − 1 the depth spans an increasingly small range since the sample of N points becomes more and more sparse in R p , hence the principle of parsimony plays a role in determining p too. Finally, as we will see in the next section there is an interval between p min and p max where the depth pattern (consequently the ranking) remains stable so that the PTD centrality score vector can be chosen among a {PTD(D, p)} in that interval. With the expression depth pattern we mean the ordering of nodes as a function of one free parameter. A depth pattern plot has the free parameter, e.g. t on the x-axis and the values of PTD(v, G; D,t) for each node v ∈ G on the y-axis. This plot enables us to graphically see how the depth range evolves, if the top ranking nodes remain approximately the same (stable pattern) or if the ranking changes considerably, indicating noise in the embedding.
The network depth PTD(D t , ·, ·) w.r.t. diffusion distance depends also on t and the scree-plot becomes a heatmap representing the stress surface, Fig.5, still, the same considerations apply for p, while t makes the PTD(D t , ·, ·) a multi-resolution centrality measure, as discussed in 2.2. To obtain a parameter-free centrality without having to fix the values for p and/or t, simply integrate over the parameter spacepossibly removing those parameters leading to the worst embeddings -and look at the aggregate network depth. In the case of diffusion embedding one may also consider the average diffusion distance matrix, as in [7], getting rid of one parameter.
We summarise the main steps the procedure as follows: 0) Choose a network metric, remembering that the shortest-path distance provides results which are similar to the closeness centrality and that if you choose the diffusion distance the diffusion time t plays the role of a scale parameter. Compute the distance matrix D (or set of matrices {D t } t 0 ). 1) Embed the network G = (V, E) in R p , using MDS (either solving the scaling problem analytically or through the majorisation algorithm); these statistical methods are preferred, since they give a feedback on the amount of information preserved/lost in the scaling process. Field-knowledge may suggest a dimension p, a-priori. As a rule of thumb, for graphs with more than 100 nodes, avoid p < 4, remember that p depends on the MDS algorithm one is using.
2) For each space configuration of points X(G) ∈ R N×p , compute Kruskal's stress function (2.9). It will be useful for selecting a range of reasonable dimensions p. It is not mandatory to evaluate the goodness-of-fit of the embedding, but it guarantess you that the subsequent findings are not the result of chance.
3) For each space configuration of points X(G), compute the depth space PTD(G; D, p), i.e. the network depth centrality vector, by means of equations (2.5)-(2.6). If the depth space reduces to just one value, 1 N , the dimension p is too large, so you can focus on smaller embedding spaces. To rule out values of p which are too small, plot the stress-1 values against p (scree plot), values larger than 0.2 rarely lead to good results. Additionally you can look at high order depth central regions R α , if there is a high variability in the nodes belonging to these regions, the depth pattern is surely not stable, so you may think of removing too small dimensions as well. 4) "Stable" depth patterns represent an indicator for good embedding dimensions p, where stable means that there are not nodes making big jumps across the depth ranking. In other terms, high order depth regions contain more or less the same nodes (possibly swapping their order inside the NETWORK DEPTH 13 of 26 region). You can either select a single p and have a single depth space or aggregate the depth score vectors, integrating over the (whole or part of) parameter space.

5)
Given the vector of depth values, quantiles q α are easily computed and so you can find depth central regions R α , while the network median is the node with maximum depth.

Applications to Empirical Networks
In this section we apply the network depth to some examples of social, biological and infrastructure networks and compare it with the measures in Tab. 1: degree, eigenvector, closeness, betweenness and distance-weighted fragmentation centrality 6 . As already mentioned we compute only these five centralities, since they are representative of other measures and the goal of our work is not to provide another comparative study on centralities applied to real networks, but to propose a statistically grounded centrality measure. The novelty of our work lies in the approach to the well-studied problem of node importance inside a network: the network depth is a generalisation of statistical depth function to structured data, through network embedding. On the other hand, depths induce multivariate order relations and enable the generalisation of order statistics (median and quantiles) to high dimensional spaces.

Social Networks
3.1.1 Zachary's Karate Club. The first example we propose is the famous Zachary's Karate Club network [25], it is a small (N = 34) social network where Actors are members of a university karate club and weighted edges represent interactions among them outside the club in a period of three years. It has a modular structure with two main communities, the factions with heads Mr. Hi, the president, and the coach John A, into which the group split after a conflict. The plot in Fig. 5 shows the stress-1, σ , for embeddings with parameters (D t ;t, p), which is in general low (stress-1 depends also on the sample size N, [17]), with higher values for small t, p. The pattern visible in the heatmap is found also for the other networks we analysed. Diffusion embeddings for short diffusion times tend to have a larger stress value since there are many nodes equally far away from each other according to D t and these are interchangeable in terms of configuration X(G) (similarly to a star graph). As p increases σ decreases, but remember that the sample becomes more and more diluted in high-dimensional spaces.
We then move to analyse the depth patterns plotting the depth value of each node as a function of t or p. These plots represent a well suited graphical tool to study the centrality range, scores and ranking, their stability and relation with the free parameters. Statistics is known to make a massive use of graphics, in particular our idea of depth pattern plots draws inspiration from forward plots in Forward Search analysis for outliers detection, see [2]. The key idea there is to monitor quantities of interest, e.g. test statistics, as the model is fitted to data subsets of increasing size. In our case, we study depth values PTD(i, G; . . . ) as space dimension p or diffusion time t increase.
We can find PTD(D t ;t, p = 3) as a function of t, in Fig. 6 For t = 1 the depth space has a small range 1 N , 5 N , while for t 7 it remains stable, with little variations on depth levels. Mr. Hi and John A. are the most central nodes until t = 4, after which Actors 20 and 33 gain the second position in the depth ranking. These two nodes have high centrality scores for betweenness and closeness (20) and degree (33). The last Fig. 7 shows the node of the Zachary Karate Club network that has the top PTD score Actors building the network into R p minimising the sum of squared differences (stress-1) between points in the Euclidean space and original diffusion distances D t on the graph, for varying p (x−axis) and t (y−axis). Main observations: (i) stress-1 is a decreasing function of p, implying that we cannot take it alone as a measure of embedding/centrality quality. (ii) 1D embeddings are poor quality in general, also for small networks.
(i.e. the generalised median), according to diffusion distance embedding (and its parameters t, p) and shortest-path distance embedding with parameter p.

2010 Network Scientists
Network. The second social network analysed in this work is the network of Network Scientists from 2010 [10]. It is a co-authorship network, where 522 scholars are linked by weighted edges representing the relation u ∼ v if and only if there is a paper with u, v among the authors. The weight of edge {u, v} depends on the number of papers co-authored by u, v.
Depth patterns in embeddings depending on D sp and D t (Fig. 8-9) are, as expected, very different. PTD(D sp ; p) summarises the information enclosed in the matrix D sp , as closeness centrality. As a matter of fact, if we look at the five nodes with top scores according to closeness -Ravasz, R. Jeong, H. Podani, J. Szathmary, E. and Barabasi A.L. -we already find a strong overlap. The correlation analysis in subsection 3.4 highlights this pattern.
On the other hand, the ranking based on PTD(D t ;t, p) may appear quite counterintuitive unless we have a clear picture of the information summarised by our network depth. Fig. 10, together with the four dimensions topology characterisation of our PTD(D t ;t, p), shades light on the results: the diffusion distance quantifies how difficult it is that two random walkers generated in two different nodes meet Lines are coloured if, for some t, the corresponding Actors belong to the depth trimmed region R α,n (2.7) of order α = 0.9, that is if the corresponding Actors have a degree of (depth) centrality of at least 90% with respect to the distribution induced by the network.
somewhere in the network. In terms of this co-authorship network, a node which is sufficiently near to all others corresponds to a scholar which has a multifaceted co-authoring history. Several blocks are clearly visible in the heatmap Fig. 10, but remember that depths ignore multi-modality, exactly as the usual median.

Drosophila's Connectome
The Drosophila connectome [22] is a directed weighted network with 49 nodes, 43 of which are local processing units (LPUs), while the remaining 6 are interconnecting units. Edge weights represent the strength of total connections in terms of neuronal fibres between neurons inside LPUs. LPUs have their own local interneurons populations (LNs) whose fibres are limited to that region and exchange information with other LPUs via bundled tracts. Interconnecting units lack an LNs population and appear to relay information to other LPUs unmodified by any local network interactions. In [22] the authors show that this has a modular structure with five functional clusters: four of them correspond to sensory modalities and the fifth (called "pre-motor centre") is the integration centre for visuo-locomotor behaviours, which is important for decision-making. The inter-  20, 9 and 33. Actor 20 in particular is also the third most central node according to betweenness and closeness, while Actor 33 scores three in the degree ranking. For small values of t, ignoring poor quality embeddings in R and R 2 , the most central vertex is Mr. Hi and the same is also true in the shortest-path embedding spaces. Missing tiles indicate that, in those embeddings, all points lie on a convex hull and it is hence meaningless to compute their depth. modular connectors are SPP, DMP, VLP-D, IDFB.
In Fig.11 we show the depth patterns in in-going and out-going geodesic embeddings with p 3. Since the directed shortest-path distance matrix is no more symmetric we get the space configuration of nodes through unfolding [6]. Lines carry a label if the corresponding node belongs to the depth region of order α = 0.90. We can compare the centrality patterns with the ranking based on node polarity, the difference between in-and out-strength; positive and negative values indicate "recievers" and "senders" respectively. Like degree, this is a very local measure, instead the PTD depends on the global distance induced space configuration. Observe that the aggregate PTD(D sp−out ) depth trimmed region of order α = 0.90 contains the following LPUs: PCB, AL, CMP, NOD and vlp-v; while for PTD(D sp−in ) the deepest LPUs are mb, og, NOD, AMMC, ammc. PCB, protocerebral bridge, and noduli(NOD), belong to the integration center for visuo-locomotor behaviours in all anthropods; antennal lobes (AL) and antennal mechano-sensory and motor center (AMMC) (with their respective counterparts in the opposite hemisphere) are first-layer sensory (resp. for olfaction and hearing) LPUs, while OG, Optic Glomerulus, is one of six interconnecting units in the Drosophila's brain, which appear to relay information to other LPUs unmodified by any local network interactions [22]. It is thus not surprising that OG/og is among the deepest node both w.r.t. in-coming and out-going connectivity. AMMC is a strong sender according to node polarity, but it also in the downstream of ammc and other senders (PAN/pan, AL/al, OG/og among others) and of LPUs with high degree/strength like DMP/dmp.
A different scenario appears if we look at depth patterns w.r.t. diffusion embeddings in Fig.12a: the top 10% of the depth is quite stable with dmp (Dorsomedial Protocerebrum) being the deepest node for most of the time. If we average over the whole parameter space, so that we do not even need to decide which pair (t, p) to use, and look at the aggregate PTD(D t ) the depth trimmed region of level α = 0.75 consists of: DMP, FB (Fanshaped Body), IDFP (Inferior Dorsofrontal Protocerebrum), SOG (Subesophageal Ganglion), SPP (Superpenduncular Protocerebrum), VLP-D (Ventrolateral Protocerebrum, Dorsal part), ccp (Caudalcentral Protocerebrum), dmp, idfp, sog, spp, vlp-d, vmp (Ventromedial Protocerebrum). Interestingly there is a large overlap with Shih's Global Centrality (GC), the result of a consensus among centralities: for each type of measure the top 12 out of 49 nodes (top quartile) earned 1 point and then by summing up one gets the node GC. One significant methodological difference between aggregate PTD(D t ) and GC is that we do not arbitrarily bestow points, we use sample depth quantiles, getting a continuous measure that can be interpreted in the same way as all other depths.

3.2.2
Protein-protein network Protein-protein interaction network of the yeast S. cervisiae [13]. This is a medium size network (1870 nodes, from which 1458 in the largest connected component), edges are undirected and unweighted, while nodes carry metadata on the protein lethality. According with lethality-centrality hypothesis, central nodes are important for cellular functioning. Nodes with highest degree (top 0.5%) are YDL100C (Non-Lethal), YLR423C (Non-Lethal), STE50 (Non-Lethal), SNP1 (Lethal), TEM1 (Unknown), LSM8 (Lethal), SEC17 (Lethal), Nup1 (Unknown). Fig. 13 shows the depth patterns for PTD(D t ,t, p = 5). Without a detailed field knowledge it is meaningless to discuss on the importance of these proteins, however from the data we can infer that the protein-protein interaction network is disassortative, which means that Lethal protein are not necessarily connected only to other Lethal proteins. This implies that the spatial distribution of points embedded according to diffusion (or shortest-path distance) does not reflect the Lethality-communities.

Infrastructure Networks
The Euro Roads infrastructure network [23] is the international E-road network, located mostly in Europe. It is an undirected network, with nodes representing cities and edges between nodes denoting E-road connections. We focus on its largest connected component, containing N = 1039 (over the total 1177) vertices.
This real network lives in R 3 , but as we try to embed it into R 3 , solving analytically the classical MDS problem with distance matrix D sp and p = 3, the result is poor with an explained variance of FIG. 10: Average diffusion distance matrix for the collaboration network of Network Scientists 2010. On the y-axis the 2.5% deepest nodes (depth region of order 97.5-percentile), on the x-axis the top 2.5% scholars w.r.t. betweenness centrality. We can see that there are several blue blocks across the diagonal representing groups of closely connected nodes. Deep nodes in the diffusion embedding, on the other hand, are not extremely close to anyone, but quite near to everyone. 74.9%. This agrees with our expectations, since the shortest-path distance is not equivalent to the Euclidean distance. Using the majorisation approach described in 2.3 and the non-metric MDS, a more appropriate embedding dimension according to stress evaluation is p = 5, which yields a stress-1 of 0.038. The cities in the depth region R α=0.99 (top 1%) based on PTD(D sp ; p = 5) are reported in Tab. 2, together with the top 1% of the common centrality measures used in this study.
We also embed the Euro-road network based on diffusion distance, with t ∈ {1, . . . , 20} and p{∈ 2, . . . , 30}. Fig. 14 shows the depth patterns for the network depth in R 5 and a subset of diffusion times. We can see a claer shift from "local" to "global" centrality of cities as t grows. For small t the nodes in the top 0.5% are Le Havre, Eupen, Arlon, Leeds, Tallinn and Paris (this remains true also for other values of p), with a high overlap with the eigenvector centrality. There is, however, a second stack of lines (nodes) with low depth values for small t, which gain depth as t grows. When t = 20 the depth central region of level α = 0.995 looks very different: the median nodes are Kingston upon Hull and Brussels, followed by Voss, Korczowa, Osnabrück and Olpe. This ranking is more similar to the rankings given by the betweenness, closeness and even PTD(D sp ; p = 5). The similarity between PTD(D t ;t, p) and PTD(D sp ; p) is a consequence of the sparsity of the network. The diffusion distance between two nodes is small if they are connected by many or by short   paths and since the network has a low edge density, the two metrics have similar patterns.

Correlations Analysis of Centrality Measures
In this subsection we perform a simple correlation analysis, by means of Spearman's rank correlation coefficients. For three of the previously studied networks and for each pair C 1 ,C 2 of centrality measures, we compute the Spearman's rank correlation coefficient r s ; additionally a two-sided statistical test of level α = 0.05 for H 1 : r s = 0, see [12,Sec. 8.5] for statistical details. Fig. 15 shows the correlation matrices for the three networks analysed in this paper, which share a common pattern. Indeed, blocks of strongly correlated measures emerge: for different values of parameters p and t the families {PTD(D sp ; p)} p and {PTD(D t ;t, p)} t,p have strong inter-family correlations and weak or even non-significant intra-families correlations. Further, closeness falls generally in the group {PTD(D sp ; p)} p while degree and eigenvector correlate more with {PTD(D t ;t, p)} t,p , as expected from the categorisation of centralities using the graph-theoretic framework of centrality measures [4]. Fig. 16 shows an extensive comparison in terms of mean Spearman's correlation between the network depth in different diffusion embeddings and common centralities, on more than 50 real networks (mostly from KONECT [15] repository). Correlations coefficients for which the p−value for H 0 : r s = 0 Nodes that lie in R α=0.995 for some value of t have a coloured line, while for each value of t the nodes in the central region for that t, R α=0.995 (t), are also labelled. We can see that median and highly deep nodes are almost always Non-Lethal.
is larger than 0.05 have been set to zero. The correlations do not show clear category-or size-based patterns. Correlations are generally moderate (below 0.7). One network stands out particularly: the interareal cortical network of the macaque [18].

Conclusions
We have shown that it is possible to extend the definition of median to structured data like networks, through a newly defined statistical data depth. This also provides a continuous centrality measure, the network depth, which depends on a network embedding step. This step introduce one parameter p -in case of an embedding based on geodesics -or two parameters p,t if we use the diffusion distance. The dimension of the embedding space p affects the quality of the network depth, but a reasonable p can be selected without much effort, as we have shown at the end of Sec. 2.3 and in the applications. The diffusion time t is more than a free-parameter, it is a resolution parameter, it allows us to move from the micro to the macro-scale, underlying locally central nodes, e.g. centres of communities, or community  connectors (see Drosophila's connectome in Sec. 3). The network depth arises from an unprecedented approach: the definition of the centrality does not follow the definition of the importance a node, it is the application of a pure statistical methodology, the generalisation of order relations through depth functions. For this reason and to prove that the network depth is a network centrality measure, we place it in the graph-theoretical taxonomy of centralities by Borgatti and Everett, theoretically inferring its properties. We have shown that the PTD is a radial measure and hence it looks for a global centre. Furthermore, which distance we choose to evaluate on the network has an important role in determining the distribution of points in space. Therefore, the depth-based node ranking has to be interpreted in light of the metric considered.
In this work we have analysed different social, biological and infrastructural systems to understand if the centrality score induced by the Projected Tukey depth reveals nodes with special or relevant meaning not detected by more traditional centrality measures, such as degree and eigenvector -radial, volume measures -, closeness -a radial, length centrality summarising information given by length of shortest paths -, or betweenness and Borgatti's fragmentation, two medial measures of centrality focusing on paths passing through each node. When diffusion geometry is used to embed the networks, we find our depth to unveil central nodes with a heterogeneous connectivity across the network. In other words, nodes that, instead of having a strong belonging to a community, are well-mixed inside the network structure.
Our results provide a grounded framework for the definition of median node and contours in complex networks.


Online Sensor Hallucination via Knowledge Distillation for Multimodal Image Classification

We deal with the problem of information fusion driven satellite image/scene classification and propose a generic hallucination architecture considering that all the available sensor information are present during training while some of the image modalities may be absent while testing. It is well-known that different sensors are capable of capturing complementary information for a given geographical area and a classification module incorporating information from all the sources are expected to produce an improved performance as compared to considering only a subset of the modalities. However, the classical classifier systems inherently require all the features used to train the module to be present for the test instances as well, which may not always be possible for typical remote sensing applications (say, disaster management). As a remedy, we provide a robust solution in terms of a hallucination module that can approximate the missing modalities from the available ones during the decision-making stage. In order to ensure better knowledge transfer during modality hallucination, we explicitly incorporate concepts of knowledge distillation for the purpose of exploring the privileged (side) information in our framework and subsequently introduce an intuitive modular training approach. The proposed network is evaluated extensively on a large-scale corpus of PAN-MS image pairs (scene recognition) as well as on a benchmark hyperspectral image dataset (image classification) where we follow different experimental scenarios and find that the proposed hallucination based module indeed is capable of capturing the multi-source information, albeit the explicit absence of some of the sensor information, and aid in improved scene characterization.

Introduction
The current generation has witnessed the abundance of satellite RS images acquired by sensors with varied resolution (both spatial and spectral). Additionally, it is nowadays possible to obtain images of a given area on the ground almost on a daily basis, thus assisting in (near) real-time monitoring of the Earth's surface. Notice that images acquired by sensors with different characteristics showcase complementary information regarding the underlying scene which, if fused properly, is expected to better analyze the ground objects. Historically, the notion of RS scene recognition is considered to be the prime ingredient of a number of RS applications which basically associates semantic meanings to the respective scenes. Precisely, such a semantic tagging invariably helps in applications related to military, urban planning, disaster mapping and many more.
The task of scene recognition ( [1]) is inherently different from the classical problem of image classification in RS which aims at predicting semantic land-cover classes at the pixel level ( [2], [3]). However, the focus on regions or sites is more important than analyzing the scenes at the pixel level in a number of applications and in particular for images with a very fine spatial resolution where very limited information is contained in the individual pixels. Nonetheless, both the tasks can be accomplished under the availability of different levels of supervision, namely, fully-supervised, semi-supervised, and unsupervised. The RS literature is traditionally rich in classification algorithms which, by definition, work on features extracted from the pixels/images. Mostly, such features or a combination of a number of base features are extracted explicitly from the image under consideration, thus giving rise to the notion of uni-modal feature extraction. Majority of the work in current literature follows this approach to address the classification problem ( [4], [5]).
On the other hand, with the availability of comparatively low-cost sensors and due to the successful launch of several international satellite missions (Sentinel 1 and 2, for example), multi-modal data are rather easier to obtain nowadays. Considering the complementary information captured by diverse sources (multi-spectral and PAN, for example), the performance of a given classifier system can substantially be boosted if this rich information can intuitively be combined since such a multi-modal fusion strategy is expected to enhance the discriminative capabilities of the input representations. As a naive example, while the panchromatic (PAN) images show better spatial resolution, it is understood that the corresponding multi-spectral (MS) images depicting the same ground area are rich in spectral properties. Hence, a weighted combination of both the modalities is always encouraged than considering each of them separately for a typical scene understanding task. Likewise, multi-spectral or hyperspectral imaging involves many measurements along the visible spectrum for each pixel and helps in coarse to fine level object analysis while synthetic aperture radar (SAR) observations are nearly insensitive to cloud cover and provide a completely different aspect regarding the back-scattering of the surface materials.
Nonetheless, it is always expected that the individual image features which are likely to be fused should possess enough generalization proficiency. While for typical RS image classification, the spectral bands are in general considered as the feature descriptors, the extraction of image level features for scene recognition is non-trivial. The feature extraction paradigms have evolved from the hand-crafted ad-hoc features to the recent data-driven deep feature descriptors. Amongst the ad-hoc feature descriptors, SIFT, HOG, based local features are encoded with techniques like bag-ofvisual-words, Fisher vector, VLAD etc. to obtain the image level representations ( [6], [7]). On the other hand, global scene level descriptors (GIST) can directly accumulate the responses for a number of filters into a unified feature representation. Besides, primitive spectral and texture statistics like Gray level co-occurrence matrix (GLCM) and Local binary patterns (LBP) are also popular for RS images ( [8], [9]). However, the advent of deep learning techniques has largely sidelined the ad-hoc descriptors for most of the visual inference tasks ( [10], [11], [12]).
In particular, the deep convolutional neural networks (CNN) based descriptors are capable of learning very generic to very specific abstract feature hierarchy from the images and thus possess improved discriminative capability. Deep CNN based techniques, which can deal with both the single as well as multi-modal data are prevalent in the literature ( [13], [14], [15], [16]). CNN based feature learning from a single modality is straight-forward as it involves the learning of a number of consecutive convolutional kernels in an end to end fashion. However, the design is non-trivial while dealing with multiple modalities. Apart from learning abstract features from multiple modalities in terms of fully-connected networks, attempts are made to augment multiple modalities prior to feeding the features to the classification networks. A number of recent endeavors have dealt with the problem of multi-modal feature fusion for RS images as mentioned in ( [17]).
However, the existing multi-modal fusion based CNN architectures require the presence of all the modalities during training and testing which may turn out to be a hard constraint for a number of applications particularly involving time-critical disaster management systems. Due to the variations in the temporal resolution of different sensors, it may not be possible to obtain all the different images immediately after a disaster takes place. Notwithstanding this limitation, the consideration of different images remains instrumental for efficient and quick response arrangement. One likely scenario in this regard can be formulated when all the modalities are present during training while one or more modalities are missing during testing or on field implementation.
Inspired by the aforementioned discussions, we propose a deep CNN based hallucination model for multi-modal image fusion where the model is trained with all the modalities where the image level abstract feature descriptors are learned in an end-to-end fashion. We basically work in the privileged information setting where we have additional data available during training procedure only but not while testing. However, we handle the absence of part of the input modalities during testing by separately modeling a hallucination network which is capable of approximating the missing modalities from the available ones on the fly. Ideally, the considered problem definition shares some ideological similarities with the classical missing value problem in RS. However, here we are interested in a more generic learning based scenario driven by the paradigms of modality hallucination and knowledge distillation jointly. To the best of our knowledge, there is an existing endeavor for RS modality hallucination ( [18]), but none on distilled hallucination. However as opposed to ( [18]), we judiciously incorporate the paradigm of knowledge distillation in our framework (see figure 4 for effect of distillation), apart from following a more intuitive training strategy. We summarize the noteworthy contributions of the paper in the following: • We propose a multi-modal fusion based scene recognition framework for RS data. Although we consider the presence of all the modalities during model training, some of the modalities are considered to be absent during testing. In order to approximate the missing modalities from the available ones, we design an intuitive deep hallucination network. • By incorporating the generalized distillation framework in the hallucination module, we demonstrate how an improved hallucination performance can be achieved. • We conduct extensive experiments on a large-scale dataset consisting of multi-spectral and PAN image pairs for the purpose of scene recognition and on a benchmark hyper-spectral dataset for pixel classification. We observe that the distilled hallucination network is capable of producing recognition accuracies comparable to a scenario when all the modalities are present during testing. In addition, we also show how the learning capacity of the hallucination network varies with changes in various model parameters.

Related work
Considering the focus of the paper, here we discuss briefly on i) multi-modal fusion based scene recognition in RS, and ii) knowledge distillation based approaches in machine learning, and iii) idea of hallucination networks.
The usefulness of employing data from multiple modalities was demonstrated by early works ( [19], [20]) well before deep neural networks became popular. This included thermal image fusion with LiDAR data for object recognition ( [21]) along with scene interpretation ( [22]) and multi-source image classification ( [23]) with optical and SAR data fusion to mention a few. In remote sensing, multi-modal fusion has been used for problems including urban mapping ( [24]), oil slick detection ( [25]), forest studies ( [26]) and disaster management ( [27]). There are also works on generating other modalities like Digital Surface Models (DSM) and Digital Elevation Models (DEM) ( [28]) to name a few.
The notion of two-stream networks is popular in fusing multi-modal information: RGB and optical flow image for video analysis ( [29]), Supervision transfer across RGB and depth ( [30]), spectralspatial data integration in RS ( [31]), RGB and depth images for scene understanding ( [32]), to name a few. By design, the two-stream network architecture is comprised of two separate neural networks and outputs from each of them are fused using an appropriate fusion strategy. As expected, each branch of the two-stream network, or a stream, typically deals with instances depicting only one modality. Recently, a number of endeavors can be observed in the area of hyperspectral image analysis for efficiently integrating both the spectral and spatial information ( [31]).
Besides, [33] propose a two-stream convolutional network for urban scene classification using PolSAR and hyperspectral data. [34] use two steam network architecture for pan-sharpening using panchromatic in one stream and multi-spectral data in the other and later fuse the outputs using a concatenation and fusion layer. There have been other fusion approaches such as max, sum and average but all these can be subsumed by a learning based fusion strategy. It is found that the learning based feature fusion technique offers improved performance than the fusion strategies following max or average pooling. Recently [35] propose a deep learning based multi-scale approach towards multi-modal semantic labeling with LIDAR and multi-spectral data.
There have been previous works to address the missing modality problem using neural networks. [36] propose a data/modality imputation approach and this offline imputed dataset is used for training a classifier. We on the other hand propose a to hallucinate missing modality on the fly in an online fashion. [37] use all the modalities for feature learning step but use one modality for supervised training step after which the testing with that modality missing. We in contrast use all available modalities for training and only one for testing. [38] propose a autoencoder based architectures for learning representations for data and in turn maximize the correlation of features across modalities. The method cannot be trained end to end and has a network for representation learning with SVM or classifier being modeled henceforth on the projections. Our method is end-to-end trainable and due to this nature can be quickly fine tuned to new data for improved hallucination performance as per need.
Before commencing the discussion on the existing knowledge distillation approaches, we briefly mention the notion of model training with side information for better inference. Precisely, training with privileged information consists of learning a student model in presence of additional information provided by a teacher model ( [39]). The idea is that the additional inputs provided by the teacher are expected to help the student to learn a better model. In the privileged information setting the additional inputs from the teacher are available to the student network only at the training time. At the test time the student network is expected to perform better inference pertaining to the initial teacher guided training process. [40] recently proposed a graph-based distillation method to learn in presence of side information for action detection and action recognition tasks in videos.
Knowledge distillation refers to the process of intelligently transferring the knowledge of an already trained, larger/complex learning module or an ensemble of several models to a rather smaller counterpart which is expected to be easier to comprehend and deploy. Knowledge distillation for machine learning models was first introduced by [41] while [42] presented a more general approach towards distillation within the scope of a feed-forward neural network and for the purpose of handwritten digit recognition. [41] considers the operations directly on the logits obtained from the softmax layer for learning the smaller model and minimize the squared difference between the logits from cumbersome and the smaller model. [42] scale the final logits by a temperature parameter to obtain soft thresholds and these are used to train the smaller model. They also show that matching the logits as done by [41], is a special case of the technique they propose.
However, from the perspective of the distillation framework mentioned in [42], the knowledge gathered by a deep network is not interpreted merely as the network parameters (weights and biases) obtained during training but also includes the final logits generated by the network corresponding to the inputs. Precisely, [42] propose a much general methodology for distillation of knowledge from a larger teacher network to a drastically smaller student network by using a parameter termed as the temperature. Essentially, the output logits of teacher network are normalized with respect to the temperature parameter and subsequently the softmax operation is carried out on them in order to obtain a softer threshold which leads to an improved and more general knowledge transfer to the student network. The outputs of neural networks are typically class probabilities which are obtained by using the softmax function which converts the output logits z i , computed for each class into a probability, by comparing it to other logits.
Consider a scenario where we have data from multiple modalities available during training, but while testing only one or a subset of these are available. Hallucination networks are trained to learn to mimic the modality that is absent by training with the hallucination loss. This paradigm can readily be applied to RS applications considering that the missing data problem is quite frequent in many scenarios. However, there exists only one previous work on hallucination networks by [18] which particularly hallucinates the missing Digital Surface Model (DSM) data in a remote sensing dataset and demonstrates improvement over an RGB based standard classification approach. However, the model is validated on a rather constrained scenario with limited samples and there has been no effort on knowledge distillation.
As opposed to the previous endeavor, ours is one of the initial approaches introducing the idea of knowledge distillation for missing data prediction in RS. Combining the concept of distillation and privileged information we demonstrate a two-stream network that can be trained with two modalities and can hallucinate the missing modality at test time. Note that the data from the missing modality during testing is what we consider as the so-called privileged information in the training stage. We transfer the knowledge learned by the stream that is trained on the missing modality at test time to another network which is referred to as hallucination. The approach presented in this paper is influenced by the notion of generalized distillation framework [32] which is considered for the purpose of human action recognition combining RGB and depth sensors. Apart from introducing the idea to the RS community, we specifically update their distillation framework with additional loss measure, which is showcased to boost the classification performance substantially.
3 Proposed Approach

Preliminaries and Problem definition
Let us consider a scenario where we have training data from multiple modalities: where M denotes the total number of sensory modalities. In the classical supervised learning setting, one would prefer to design a classifier system by training the model with a subset or all of these M features and subsequently evaluate the model on similar data items so as to preserve the consistency of the feature space. As opposed to this scenario, we aim to leverage data from all the M available modalities to train a robust classifier but subsequently ensure that consistent or even improved generalization performance is guaranteed while only a subset of the M modalities is available during testing. Precisely, let X i denote the available data samples from the i th modality with Y being the corresponding class label while X j represents the data from missing sensor modality during testing considering that we deal with two different data sources in a given scenario. In this circumstances, we are interested in learning three functional mappings: i) f comb : . . , C} scene themes where each image is captured by two separate sensors showcasing two completely different perspectives: X i and X j may denote the multi-spectral and PAN images for a given area or a pair of non-overlapping subsets of the original set of bands in hyperspectral image. For the next section, since we experiment with two data modalities, we assume X 1 and X 2 to represent panchromatic and multi-spectral data respectively and the network stream working with them are termed as PAN-Net and MS-Net for brevity. Note that we describe the proposed model using this setup, however, needless to mention, X i and X j may represent any different information sources in general. In the following, we describe the training process followed under generalized distillation for the proposed pipeline.

Training in detail
We follow a four-step modular procedure for an intuitive training of our two-stream network. In particular, in the first two stages of training, we learn the teacher network which is used in the following two steps in order to model the hallucination network. In the typical distillation setting, a teacher network is the particular network whose outputs essentially guide the training of the second network which is termed as the student network. Notice that with both X 1 and X 2 data streams being present, we initially learn separate network streams for each of the modalities followed by a joint learning combining both X 1 and X 2 in order to model f comb . Subsequently, the stream for data modality that is known to be missing during testing is considered to be the teacher network and is deployed to learn the hallucination stream, f hall efficiently in conjunction with the generalized distillation loss. The trained f hall , in turn, is used in the hallucinated two-stream network for classification which essentially has only the available modality as input and approximate the missing modality through F hall as if the testing is performed using both X 1 and X 2 . The individual loss measures for each of the component networks are described in the following and Figure 1 depicts a flow for the training stages.

Step-1: Learning the individual streams
As aforementioned, we train each of the PAN-Net and MS-Net streams separately at the beginning. Essentially, these are designed as typical classification networks which perform classification solely based on the individual information source X 1 and X 2 , respectively. Both of these networks are trained using cross entropy loss(L CE ) for classification as given below: where, Y k andŶ k are the original and predicted labels in the typical one-hot vector fashion, respectively. Step-3: Learning the hallucination network using the L GD and soft targets from frozen MS-Net, and finally (d) Step-4: Fine-tuning the hallucinated two-stream network with weights of the Hall-Net being frozen. The final step also represents how the network is deployed while testing. The striped networks represent that their weights are frozen in that step.

3.2.2
Step-2: Learning the two-stream network In the second phase of training, we build our two-stream network (Figure 2), f comb which is trained using data from both X 1 and X 2 . The first stream uses the architecture of the PAN-Net and the second stream uses the architecture of MS-Net, respectively. We fuse the two streams at the end and the new representations are classified following a learned late fusion strategy.
Specifically, the final class probabilities from both the streams are eventually concatenated and fed into a fusion layer. We employ a learned late fusion strategy, suppose p 1 , p 2 ∈ R N ×1 are the probability score from the two streams respectively, we learn a set of fusion weights w ∈ R N ×2N , such that w * p = w 1 * p 1 + w 2 * p 2 , where p ∈ R 2N ×1 is the concatenated probability vector. We want to learn a fusion layer such that, where Y are the final labels and γ is the regularization parameter. The symbol * represents the element-wise multiplication and . is the standard l 2 norm.
For efficiently train the two-stream network, we first initialize the weights of the individual stream with the respective weights learned in Step-1 and the whole network is subsequently trained using cross-entropy loss. This strategy aids in better convergence in comparison to the random initialization. This completes the learning of our teacher network for subsequent learning of hallucination network, which is the MS-Net for this description.

Step-3: Learning the hallucination network
In this step we train our hallucination stream, f hall , using the generalized distillation process. Assume for instance that we wish to hallucinate the X 2 data stream while testing given data from X 1 . In this  case, we use the MS-Net architecture with data from X 2 as input which will act as a teacher network for the hallucination net which we alternatively term as the Hall-net. The Hall-Net or f hall follows the previous PAN-Net architecture and uses data from X 1 to generate corresponding the MS-Net activations. As shown in Figure 1c, first we initialize the teacher network with the weights of MS-Net form Step-2 and freeze these weights. Next, we train the Hall-Net stream using the generalized distillation loss, to generate the multi-spectral counterparts for a given PAN instance.
In a standard knowledge distillation setting, we train the student model using a softer target distribution for the data samples via corresponding soft targets from the teacher network with the same high temperature in its softmax as the softmax of the student model. Let z i and s i be the logits and soft thresholds respectively, and T is the temperature parameter. The outputs of the final layer in a neural network, before taking a softmax, are called the logits. The soft thresholds are obtained from the logits as follows, In the usual classification setting, the temperature parameter is typically set to 1 to obtain the final probability distribution. But using a higher temperature value produce a softer probability distribution over classes which is employed during the knowledge distillation process.
Finally, this network is trained using the generalized distillation loss which comprises of the knowledge distillation loss term (L KD ) and hallucination loss term(L H ) as given below. The hallucination loss is fundamentally a mean square error loss which minimizes the l 2 distance between the two soft thresholds, s X2 and s Hall , of teacher and hallucination networks, respectively, which are obtained using a specific temperature.
The knowledge distillation loss, on the other hand, consists of two parts. Where first being the Kullback-Leibler divergence (D KL ) between the soft thresholds, s X2 and s Hall . This is employed since we want to bring the two soft target distribution closer, allowing the student model to achieve better generalization. And the second part is a cross entropy loss between the hard labels, z X2 and z Hall , from each of the two networks as we also want the hard target distribution to be as close as possible. The KL Divergence loss (L KL ) between the soft thresholds can be written as follows: The complete knowledge distillation loss (L KD ) can now be formulated using the hyper-parameter λ as follows: Combining the hallucination and knowledge distillation loss we write the final generalized distillation loss (L GD ) as below: where α is chosen to weigh in both these components suitably.
Note that the performance of the initial MS-Net stream is what we would ideally want the hallucination network to achieve eventually.

Step-4: Learning the hallucinated(pseudo) two-stream network
As the final phase of training and after having trained the f hall network to hallucinate the missing data stream X 2 , we use it in the two-stream network, f comb from Step-2. We replace the missing modality stream of f comb with f hall to learn the final f class network. In particular, we initialize the first stream with weights from PAN-Net in Step-2 and the second stream with weights from Hall-Net learned in the previous step. Notice that the weights of the hallucination stream are frozen and subsequently we only fine-tune this hallucinated two-stream network with cross-entropy loss. This network is considered to be ready to perform classification with missing X 2 during decision making. Again, note that the performance of the initial two-stream network form Step-2 will be what we would ideally want to achieve with this hallucinated two-stream model.

Testing
The testing is straightforward given the hallucinated two-stream model from Step-4. We input X 1 in the PAN-net and f hall (X 1 ) in the HALL-net, respectively. The final classification scores are obtained by following the aforementioned late-fusion strategy mentioned in our original two-stream model (Step-2).

Dataset description
We use two datasets for analyzing the efficacy of the proposed hallucination model: the first one containing panchromatic images of aerial scenes along with the corresponding multi-spectral data. This dataset consists of a total of 80000 PAN-MS pairs from eight different land-cover categories [43]. The images are acquired by the GF-1 panchromatic and GF-1 multi-spectral sensors, respectively. The panchromatic images are of size 128 × 128, with a spatial resolution of 2m and single spectral band. The respective multi-spectral data has a resolution of 64 × 64, with a resolution of 8m, and consists of a total of four spectral bands. The classes considered here are aquafarm, cloud, forest, high building, low building, farmland, river, and water.
The second is the Pavia University hyperspectral datacube acquired by the ROSIS sensor 1 . It is taken over the Pavia region in northern Italy and consists of 103 spectral reflectance bands with 610 × 340 pixels each post discarding bands with unimportant information contents. The datacube bands have a geometric resolution of 1.3 meters and the area contains pixels from nine land-cover classes namely, Asphalt, Meadows, Gravel, Trees, Painted metal sheets, Bare Soil, Bitumen, Self-Blocking Bricks and Shadows

Model architectures
Here we summarize the typical experimental setup considered for both the datasets. For experiments on the PAN-MS dataset, we alternatively consider the PAN and MS modality to be missing during testing. For the hyper-spectral data, we consider a subset of bands to be missing and aim at hallucinating them given the remaining bands.
For a given set of experiments conducted on the PAN-MS data for scene classification, we train one of the streams of our initial two-stream network with the panchromatic band and the second with multi-spectral data and provide only the panchromatic data for testing. Needless to mention, we focus to hallucinate the multi-spectral data for each test sample in this case. Notice that both of the network streams considered here are based on convolutional neural networks as our inputs are image objects. By design principle, the PAN-Net consists of four convolutional layers followed by four fully connected layers. Note that the structure is fixed heuristically. On the other hand, the MS-net contains three convolutional layers and three subsequent fully-connected layers. In a separate set of experiments, we hallucinate the PAN images given the MS data-cubes.
Our network is inspired from the alexnet architecture [11] and since our dataset is relatively small we use a rather minimalistic version of the same. We tested the method with 3,4, and 5 convolution layers and find that the performances of the models do not drastically vary. Hence, we report the performance of the model with 3 hidden convolution layers. The late fusion scheme is standard in multi-stream deep learning models. One of our goals here was to keep the model simple and yet delivering good performance.
For the experiments concerning the hyper-spectral data, we use first 52 bands (according to indices) of Pavia University datacube as our first sensor modality (X 1 ) and the goal is to hallucinate the subsequent 53 bands during testing. Likewise in another set of experiments, we reverse the interpretation of both the streams and hallucinate the first 52 bands while testing with only the last 53.
As far as the network design is concerned, we use a 3 layer network consisting of fully-connected layers for modeling the two network streams. As opposed to the PAN-MS case where the input data dimensionality are different, both the streams, in this case, follows identical architecture.

Experiments protocols
For the case of PAN-MS data set, in Step-1, the PAN and MS network streams are individually trained with a learning rate or 3 × 10 −4 . The same parameter values are also used to train the hallucination stream as well. For Step-2 where the initial two-stream network is trained and Step-4 where the hallucinated two-stream network is trained, the learning is performed with a learning rate of 1.5 × 10 −4 . Each of the networks is trained for 100 iterations while the initial two-stream network is trained for 200 iterations. For experiments on HS data, each of the two streams and hallucination stream is learned with a learning rate of 6 × 10 −4 and the initial and hallucinated two-stream networks are learned using a learning rate of 3 × 10 −4 .
The performance is measured using the overall and class-wise accuracies of classification. We would like to highlight that we would ideally like the performance of Hallucination network to be as close to the performance of networks stream for the missing modality. Similarly, we would like the performance of the hallucinated two-stream network to be higher than that of the available modality stream but of course as close as possible to the original two-stream network trained with both the data modalities.
We also extensively study the effects of various parameters in modeling our network, namely, α, λ and T along with the network depth.

Results on PAN-MS Data
We perform our experiments on the panchromatic and multi-spectral pair dataset with the convolutional network architecture as shown in figure 2. We use different training ratios for our experiments and study the performance variations with respect to them. The results are shown in Table 1 and Table  2 for hallucinating multi-spectral and panchromatic streams, while considering 50% of the entire data for training, from which 5% is used for validation. In these tables, the PAN-Net and MS-Net refer to the individual network streams for each of the modalities as in Step-1. The Two-stream Net is the network from Step-2, where the two individual streams are fused. The MS Hall-Net and PAN Hall-Net refer to the Hallucination networks from Step-3 for hallucinating Panchromatic and multi-spectral data streams respectively. The proposed Hallucinated two-stream Net is the two-stream network from Step-4 comprising of the hallucination stream in place of the missing modality. The table 2 shows overall and class-wise accuracies of each of the streams with hallucination, from where it is evident that the multi-spectral hallucination network has indeed learned to approximate the missing modality and achieves a comparable accuracy as its teacher network which is the original multi-spectral stream. The results presented are mean and standard deviation of the of the values.
And also the hallucinated two-stream network constructed with this hallucination stream gives us a higher overall accuracy and class-wise accuracies in comparison to the panchromatic stream. Similar observations can be made from table 2 about the panchromatic hallucination network which gives higher overall accuracy than the multi-spectral stream alone. Note that the performance of the two-stream network trained in Step-2 acts as the upper bound to what the hallucinated two-stream network is expected to achieve.
We can observe from results in table 2 that for multi-spectral data hallucination, out of the eight classes, the Cloud class is best classified by the each of the streams with an accuracy going over 99% in each case. The Water class, on the other hand, is best hallucinated by the Hallucination Net and the hallucinated two-stream network achieves a performance close to the original two-stream accuracy for it. Also, we can observe that the class farmland is the most difficult to classify and is extremely difficult to hallucinate as well with its classification accuracy being the least among the available classes but still has a performance jump of around +1% in the hallucinated two-stream network in comparison to the PAN-Net. For panchromatic data hallucination from results in table 2, the Cloud class is still best classified by all the networks, while the class River remains most difficult to hallucinate having a lower classification accuracy than even the Farm land class. The classification of the Farmland class suffers in this case too as the performance of the initial PAN and MS network streams is far from being good for it in particular. For this case, the jump in classification accuracy of Farmland class is around 4% from MS-Net to hallucinated two-stream model.
We further present the error and accuracy curves of the four networks while training namely, PAN-net, MS-net, Two-Stream Net and the hallucinated two-stream net, for 100, 100, 200 and 100 training epochs, respectively. We consider the best performing network snapshots over validation set during each of these training processes for deployment and use in subsequent stages. It can in Figure 3 be observed that the final network converges quickly since it merely requires slight fine-tuning which is performed at a lower learning rate in comparison to the other networks.

Effect of varying network parameters
We subsequently conduct experiments to study the effect of changing various hyper-parameters of our model and also consider different training-test splits in order to assess the learning capability of our network under different levels of supervision. We defined α to be the weighing factor of knowledge distillation loss in the overall loss term. Figure 4 shows how by increasing the parameter α, the accuracy of the hallucination network increases first and then decreases for higher values of α. From this, we can easily infer the advantage of the inclusion of the knowledge distillation loss term in our objective function. In particular, α = 0 implies only the hallucination loss term is used and α = 1 is for when only the knowledge distillation loss term is prioritized. We can observe from figure 4 that    the knowledge distillation loss results in an improved hallucination network. As desired, a weighted combination of these two loss measures provides an improved performance in terms of classification accuracy in comparison to just using one of them. Essentially, we observe a substantial improvement in the learning capability of the hallucination network by using the proposed generalized distillation technique instead of resorting to only the hallucination loss.
In figure 5 we can observe how the increase in the value of λ parameter changes the learning of our hallucination stream. As a part of the knowledge distillation loss, the parameter λ weighs the constituent losses of the L KD . If λ = 0, only soft threshold are used for knowledge distillation and in the case of λ = 1, only the hard targets are used for the distillation process. We can observe that soft targets give us a better hallucination performance in comparison to only the hard targets. In this case, too, we obtain a much better hallucination network by using a weighted combination of these losses instead of just resorting to using hard targets for distillation loss.

Effect of varying temperature
Although the temperature parameter is difficult to tune we found a trend in the same on averaging the results across multiple runs. We repeated the proposed MS hallucination process for multiple temperature values ranging from 1 to 100 while α and λ were each set to 0.5. The results of this are presented in Fig. 6. We can observe from the figure that the soft threshold obtained at T = 15 give best hallucination performance for this particular case. Another observation is that the performance initially increases with increase in temperature but quickly starts dropping with further increase.

Effect of varying stream depth
The next experiment involves varying the capacity of each our streams. Specifically, we reduce the number of layers in each network stream and observe how this affects the hallucination and hallucinated two-stream network performance. The corresponding performance measures are tab-    3. We observe that reducing the network depth affects the learning capacity of the hallucination network along with all of the other networks. The results are shown for the case of hallucinating multi-spectral data stream. We observe that the performance of initial MS net drop by 1.5% in comparison to the results from table 1 and that of the hallucination net drops by 3% due to the reduction in network parameters.

Effect of varying training data size
For this study, we use three different training-test splits to carry out the training process and report the observed changes in final classification performance. In particular, we use training sets comprising of 20%, 50% and 80% samples from the entire dataset. Table 4 shows the accuracies for varying training ratios for each of our networks, where we hallucinate the multi-spectral data stream. As expected, it can be observed that the accuracies of each of the networks increase with an increase in training set size. The hallucination network also improves with an increase in data and as a result, the hallucinated two-stream network always showcase a higher accuracy measure than considering solely the PAN-stream.

Results on Hyperspectral Data
We divide the Pavia University hyperspectral datacube into two parts and used them to train each of the streams in order to hallucinate the other half of the bands in the datacube. The same train-test From the tables, it can be noted that these two sub-cubes offer fairly complementary contributions to the final classification accuracy. Similar to the case of PAN-MS dataset, we can see that the hallucination network has successfully learned to mimic the missing bands of the datacube to a considerable extent. Besides. we can also observe that the hallucinated two-stream network offers an improved performance in comparison to the single stream networks trained on either half the datacube. This indicates that we can obtain an improved performance with a lesser number of bands using our hallucination architecture than a unimodal classifier.

Conclusions
We propose a hallucination network model in this paper which can deal with the problem of missing data modality during testing in a typical remote sensing classification task. It is imperative that information fusion from multiple sources provide complementary insights to a geographical area which helps in better characterization. However, in many situations, some of the modalities may remain absent. Our approach is capable of approximating such missing modalities during the inference stage. Ideally, we train our model with all the available modalities, thus giving rise to the notion of learning with privileged information. We validate the proposed framework on a large-scale PAN-MS dataset for scene recognition and a hyperspectral image for pixel classification where the effectiveness of model hallucination is established through extensive experimentation. We plan to extend this model to handle multi-temporal data as a future endeavor.