Hyperharmonic analysis for the study of high-order information-theoretic signals

A popular way to analyze the interactions within $\text{bi}{X}^N$ is to represent them as networks, where each variable $X_0,\dots ,X_N$ is represented as a node, and edges between nodes represent the strength of their interaction. A simple way to build such a network is to calculate the correlation matrix $R_{\text{bi}{X}^N}:=\left[r_{i,j}\right]$ with components given by

and use it as an adjacency matrix — either binarising its components via a threshold, or considering weighted edges. This construction, however, only captures linear relationships between the variables. More encompassing analyses often consider non-linear measures of dependency such as Shannon’s mutual information and focus on the matrix $I_{\text{bi}{X}^N}:=[I(X_i;X_j)]$ of mutual information terms, which are computed as

It is important to realise that the joint probability distribution $p\left(\text{bi}{X}^N\right)$ may contain substantial information that is not assessed by any of the pairwise marginals $p(X_i,X_j)$. An elegant way of gaining insights into this is provided by information geometry, as presented in Amari (2001) (for related discussions, see also Rosas et al. (2016)). Consider the $k$-marginals that are obtained by marginalising $p\left(\text{bi}{X}^N\right)$ over $N-k$ variables. One can see that the $k$-marginals provide a more detailed description of the system than the $(k-1)$-marginals by noting that the latter can be directly computed from the former by marginalising the corresponding variables. In contrast, the process of marginalising involves irreversible information loss, as there are many $k-1$ marginals that are consistent with a given a set of $k$-marginals.

Perhaps the simplest example of high-order statistical dependency is given by the exclusive-or (xor) logic gate. To see this, consider two independent fair coins $X_0$ and $X_1$, and let

A quick calculation shows that the mutual information matrix $I_{\text{bi}{X}^2}$ of $(X_0,X_1,X_2)$ has all its off-diagonal elements equal to zero, making it indistinguishable from an alternative situation where $X_2$ is just another independent fair coin. Therefore, put simply, any “xor-like” effects is completely ignored by constructions based only on pairwise statistics — either networks or hypergraphs.

While commonly neglected, high-order statistics has been recently proven to be instrumental in a number of systems. For example, synergies have been shown to play a key role in distributed interdependent systems, including the most complex types of elementary cellular automata (Rosas et al., 2018). It has also been suggested that high-order interactions could drive thermalisation processes within closed systems (Lindgren and Olbrich, 2017). Additionally, high-order interactions have been argued to play a key role in neural information processing (Wibral et al., 2017) and high-order brain functions (Tononi et al., 1994), being at the core of popular metrics employed in consciousness science (Mediano et al., 2019). Furthermore, it has been recently shown that high-order statistics also play a crucial role in enabling emergent phenomena (Rosas et al., 2020a).

Shannon’s mutual information is limited to capture the dependencies of two groups of variables, but cannot directly assess triple or higher interactions. The most popular non-negative multivariate extensions of the mutual information are the Total Correlation (TC) (Watanabe, 1960) and the Dual Total Correlation (DTC) (Han, 1978), which are defined as

Above, $H\left(X_i\right)=-{\sum }_{x_i}p\left(x_i\right)\log p\left(x_i\right)$ corresponds to Shannon’s entropy, $H\left(X_i|X_j\right)=H(X_i,X_j)-H\left(X_j\right)$ is the conditional Shannon entropy, and $\text{bi}{X}_{-i}^N$ is the vector of all variables except $X_i$ (i.e., $\text{bi}{X}_{-i}^N=(X_0,\dots ,X_{i-1},X_{i+1},\dots ,X_N)$). Importantly, both TC and DTC are zero if and only if all variables $X_0,\dots ,X_N$ are jointly statistically independent — i.e. if $p\left(\text{bi}{X}^N\right)={\prod }_{i=0}^{N}p\left(X_i\right)$.

Unfortunately, both TC and DTC provide metrics for high-order interdependency which are difficult to analyse together. A recent approach, introduced in Rosas et al. (2019), proposes to employ a linear transform over these two metrics to obtain the O-information and the S-information, which have more intuitive interpretations.

The O-information can be seen as a revision of the measure of neural complexity proposed by Tononi, Sporns and Edelman in Tononi et al. (1994), which provides a mathematical construction that is closer to their original desiderata (Rosas et al., 2019). In effect, the O-information is a signed metric that satisfies the following key properties:

• It is zero for systems with only pairwise interdependencies.

• It is additive over non-interactive subsystems.

• It is maximised by redundant distributions, and minimised by synergistic (“xor-like”) distributions.

Hence, $\Omega \left(\text{bi}{X}^N\right)<0$ implies a predominance of statistical synergy within the system $\text{bi}{X}^N$. Conversely, $\Omega \left(\text{bi}{X}^N\right)>0$ implies that the system $\text{bi}{X}^N$ is redundancy-dominated.

On the other hand, the S-information is an over-encompassing account of interdependencies taking place at all orders, being sometimes described as a “very mutual information” (James et al., 2011). In fact, a quick calculation shows that the S-information can be decomposed as $\Sigma \left(\text{bi}{X}^N\right)={\sum }_{i=0}^{N-1}I(X_i;X_{-i}^N)$, which can be seen as a chain rule where the interdependencies involving each variable are sequentially addressed.

In summary, while the TC and DTC provide alternative representations to the same construct, $\Omega$ and $\Sigma$ provide a complementary account of the system: the latter addressing the overall strength of interdependencies, and the former qualitatively characterising their dominant nature.

Another popular metric of high-order interdependencies in the Interaction Information, first introduced in (McGill, 1954) for systems with three variables.¹¹1The Interaction Information is closely related to the I-measures (Yeung, 1991), the co-information (Bell, 2003), and the multi-scale complexity (Bar-Yam, 2004)¹ Building on an application of the inclusion-exclusion principle to entropies, the Interaction Information of $\text{bi}{X}^n$ is a signed metric given by

where the sum is performed over all subsets $\text{bi}\gamma \subseteq \{0,\dots ,N\}$, with $|\text{bi}\gamma |$ the cardinality of $\text{bi}\gamma$, and $\text{bi}{X}^{\text{bi}\gamma }$ the vector of all variables with indices in $\text{bi}\gamma$. While this measure has a direct interpretation as redundancy minus synergy for $N=2$, it no longer reflects this balance for larger system sizes (Williams and Beer, 2010, Section V). However, the Interaction Information can still be interpreted for arbitrary $N$ under a topological formulation of information, as described in Baudot and Bennequin (2015), Baudot et al. (2019).

Other well-known metrics of high-order effects include the Redundancy-Synergy Index (Chechik et al., 2002, Timme et al., 2014), the Connected Information (Amari, 2001, Schneidman et al., 2003b), and the Partial Entropy Decomposition (Ince, 2017). Furthermore, a detailed exploration of multivariate decompositions can be found in the Partial Information Decomposition framework (Williams and Beer, 2010) and its constantly growing associated literature (Faes et al., 2017, Finn and Lizier, 2018, Ay et al., 2019, Rosas et al., 2020b, Makkeh et al., 2020).

A hypergraph $S=(V,E)$ is determined by a set of vertices $V=\{0,\dots ,N\}$ with $N\in N$ and a set of hyper-edges $E\subseteq P\left(V\right)$, where $P\left(V\right)$ is the power set of $V$. Furthermore, a hypergraph $S$ is said to be a simplicial complex if it satisfies two conditions:

• all singletons $\left\{k\right\}$ with $k\in V$ are included in $E$, and

• if $\sigma \in E$ and $\rho$ is a subset of $\sigma$, then $\rho \in E$.

Please note that the passage to simplicial complexes does not restrict the theory of hypergraphs significantly, since to every hypergraph one can assign a canonical simplicial complex by downward closure. Explicitly, this is the smallest simplicial complex that contains the hypergraph. For an introduction to graphs and hypergraphs, we refer to Berge and Minieka (1973); for a comprehensive introduction to hypergraphs in the context of complex system analysis see Johnson (2013).

If $S=(V,E)$ is a simplicial complex, the elements of $E$ are referred to as simplices; and their dimension is defined as one less than their cardinality (i.e. the number of elements they connect). This shift can be understood noting that the natural dimension of a point, corresponding to a singleton, is $0$. The subset of $E$ consisting of simplices of dimension $n$ is denoted by $S_n$. The elements of $S_n$ are typically called “$n$-simplices”, with the $0$-simplices and $1$-simplices being informally called vertices and edges, respectively. If $v_0<\cdots <v_n$ are the vertices of a simplex, we denote this simplex by $[v_0,\dots ,v_n]$.

An $n$-dimensional signal on a simplicial complex $S$ is a function $\alpha :S_n\to R$, that is to say, an assignment of a real number to each $n$-simplex. We refer to $n$-signals with $n\ge 2$ as high-order signals.

For a simplicial complex $S$ and $n\in N$, we denote by $C_n\left(S\right)$ the vector space generated by all the $n$-simplices of $S$, i.e.,

Please note that there is a natural bijection between the set of $n$-signals on $S$ and $C_n\left(S\right)$ established by

This bijection is used implicitly when referring to the elements of $C_n\left(S\right)$ as $n$-signals.²²2In the mathematics literature, the elements of $C_n\left(S\right)$ are called “real-valued $n$-chains”, but we do not use this terminology.²

Let us now consider the linear map ${\partial }_n:C_n\left(S\right)\to C_{n-1}\left(S\right)$ defined on basis elements by

with ${\hat{v}}_i$ denoting the absence of $v_i$ from the simplex. One can visualise ${\partial }_n$ geometrically in terms of the boundary of a basis element. Up to a sign, the basis elements appearing as summands on ${\partial }_n([v_0,\dots ,v_n])$ are all simplices of dimension $n-1$ contained in $[v_0,\dots ,v_n]$. Furthermore, the sign is determined in terms of the orientations of the simplices, as illustrated in Figure 1. The maps ${\partial }_n$ play a central role in algebraic topology holding a significant amount of topological information. In particular, the difference between the dimension of the kernel of ${\partial }_n$ and the dimension of the image of ${\partial }_{n+1}$ is known as the $n$-Betti number of the simplicial complex, a powerful topological invariant generalising the Euler characteristic. For a systematic treatment of these ideas, please consult Hatcher (2002).

In this work we are not only interested in the topology of $S$, but also on the “geometric” structure encoded by weights assigned to its simplices. A weighted simplicial complex is a pair $(S,w)$ where $S$ is a simplicial complex and $w=\{w_n:S_n\to R_{>0}\}$ is a set of non-negative weight functions. In the following, a weighted standard simplex (see Example 1) is referred to as a structural simplex.

Given a weighted simplicial complex $(S,w)$, one can encode the structural information provided by $w_n$ (where the subscript denotes the dimension) as an inner product on the vector space of $n$-signals on $S$. Explicitly, for any pair of $n$-simplices we have

Importantly, this inner product allow us to introduce the adjoint operator of ${\partial }_n$, which is denoted by ${\delta }_n$. That is to say, the operator ${\delta }_n:C_n\left(S\right)\to C_{n+1}\left(S\right)$ is defined by the identity

which holds for any $\alpha \in C_{n+1}\left(S\right)$ and $\beta \in C_n\left(S\right)$. Note that ${\partial }_n$ does not depend on the weights $w$, but ${\delta }_n$ does.

The classical sine and cosine functions form the basis of the spectral representation of time-domain signals. In effect, classical Fourier analysis guarantees that a large class of functions are expressible in terms of linear combination of these functions. The coefficients of this linear combination are the Fourier coefficients, which correspond to a geometric projection of the original function over this new basis. Interestingly, numerous signals of practical relevance are more compactly represented in the spectral domain, and hence Fourier analysis is often employed as a principled way to perform dimensionality-reduction. For a more extensive exposition of these ideas, we refer the reader to Bracewell (1986).

Importantly, sine and cosine functions are also eigenfunctions of the Laplace operator on the circle — a one-dimensional manifold. Put differently, a spectral representation is equivalent to a diagonalisation of the Laplace operator. Mathematicians have generalized the Laplace operator to higher-dimensional manifolds via the Laplace-de Rham operators. In this more general context, functions on a smooth manifold lie at the bottom of a sequence of higher-dimensional objects on which the corresponding Laplace-de Rham operator acts. The eigenvectors of these operators play a central role in modern geometry — most notably through Hodge theory (Hodge and Atiyah, 1989). For an exposition of these ideas we refer the reader to Morita et al. (2001).

A discrete analogue of the higher-dimensional objects on which the Laplace-de Rham operators act is provided by high-order signals defined over a simplicial complex. In this correspondence, the Laplace-de Rham operators are represented by the discrete Laplace operators first introduced in Eckmann (1944). See also Horak and Jost (2013) where a weighted version of these is presented, and Parzanchevski and Rosenthal (2017) where the connection to random walks is explored. We now introduce the particular version of the discrete Laplace-de Rham operators that we use.

Notice that we are abusing notation by omitting $w$ from the notation referencing the Laplace operators since, although ${\partial }_n$ does not depend on $w_n$, ${\delta }_n$ and therefore ${\Delta }_n$ do. For the interested reader, we remark that this Definition 2 has a strong resemblance to the form the Laplace-de Rham operator adopts in Rimmanian geometry. In effect, by denoting by $d$ the exterior derivative of differential forms and $d^{\ast }$ its adjoint, the Laplace-de Rham operator in this context can be expressed as $d{d}^{\ast }+d^{\ast }d$. Furthermore, the geometry defined by the Riemannian metric is reflected in the Laplace-de Rham operator through the operator $d^{\ast }$ only. Consult Morita et al. (2001) for further details.

The well-known graph Laplacian, a central concept in spectral graph theory (Chung et al., 1997), is equivalent to the $0$-Laplace operator defined above — when a weighted graph is regarded as a weighted simplicial complex. Importantly, the $n$-Laplace operator is self-adjoint for any weighted simplicial complex $(S,w)$, i.e.

for all $\alpha ,\beta \in C_n\left(S\right)$. This implies that ${\Delta }_n\left(\alpha \right)$ is diagonalisable; that is to say, there exists a basis of $C_n\left(S\right)$ consisting of eigenvectors of ${\Delta }_n$.³³3This result is a more general form of the well-known fact that symmetric (i.e. self-adjoint) matrices are diagonalisable.³ An $n$-Fourier basis for $(S,w)$ is a maximal set of linearly independent orthonormal (with respect to $⟨\cdot ,\cdot {⟩}_w$) eigenvectors of ${\Delta }_n$. Please note that we speak of the Fourier basis, with the understanding that there is a sign choice for each of its elements.

Finally, the hyperharmonic representation of a high-order signal defined over a weighted simplicial complex is given by its change of bases from the canonical to the Fourier one. In the case of graphs, the use of this transformation as a dimensionality-reduction method was pioneered by Belkin and Niyogi (2003), and has served as an early application of the field of graph signal processing — see for example Shuman et al. (2013), Sandryhaila and Moura (2014) or Ortega et al. (2018) for an overview of this field. In recent years, some key ideas from graph signal processing have been adapted to hypergraphs, see for example Barbarossa and Sardellitti (2020) and Schaub and Segarra (2018); however, the applicability of Fourier analysis as a compression tool of high-order signals is, to a large extend, still unexplored territory.

The proposed pipeline consists in six steps.

• From data to a distribution: The starting point of the workflow is multivariate data, which typically takes the shape of a sequence of vectors of dimension $N+1$ (e.g. successive samples of time series of the form $\text{bi}s\left(t\right)=\left(s_0\right(t),\dots ,s_N(t\left)\right)$ at different values of $t\in R$). Using these data, the first step is to construct a joint probability distribution for a $(N+1)$-dimensional vector, denoted as $p(X_0,\dots ,X_N)$.

• The structural simplex: To build the structural simplex, the first step is to use $p$ to calculate the mutual information $I(X_i;X_j)$ for each pair of variables of the system. These values are then used to define the weights on a pairwise network, where each edge correspond to one of the variables $X_0,\dots ,X_N$. Subsequently, a weighted $n$-simplex is build for each $n=2,\dots ,N$ by taking the average value of all the edges that connect two elements within the simplex.

• High-order signals: For each subset of cardinality $n+1$ of $\{X_0,\dots ,X_N\}$, one computes the O-information and S-information and arrange them as high-order signals following Example 2. These signals are stored as $(\frac{N+1}{n+1})$-dimensional vectors $\text{bi}{o}^n$ and $\text{bi}{s}^n$, representing them in the canonical basis of simplices. Note that there is a standard order on the elements of this basis, which is given by the lexicographic principle (i.e. $[v_0,\dots ,v_n]<[v_0^{\prime },\dots ,v_n^{\prime }]$ if $v_j<v_j^{\prime }$ and $v_i=v_i^{\prime }$ for all $i<j$).

• The Laplace operator: Using the weights of the structural simplex constructed in Step (2), one then builds the corresponding $n$-Laplace operator as in Definition 2. Concretely, we use the canonical basis to represent the $d=$ $(\frac{N+1}{n+1})$ weights of order $n+1$ within a $d\times d$ diagonal matrix $W_n$, and represents the linear map ${\partial }_n$ as a matrix $B_n$ of the same dimensions (for an algorithmic description of the construction of $B_n$, see B). Then, the discrete $n$-Laplace operator is represented in the canonical basis as the matrix

  $$L_n=L_n^{\text{up}}+L_n^{\text{down}},$$

  (10)

  with $L_n^{\text{up}}$ and $L_n^{\text{down}}$ given by

  	$$L_n^{\text{up}}$$	$$=W_n^{-1}{B}_n^{⊺}{W}_{n+1}{B}_n,$$
  	$$L_n^{\text{down}}$$	$$=B_{n-1}{W}_{n-1}^{-1}{B}_{n-1}^{⊺}{W}_n,$$

  where $B_n^{⊺}$ is the transpose of $B_n$. To recapitulate, Eq. (10) is the matrix representation, in the canonical basis, of Eq. (9) with respect to the inner product given by Eq. (8).

• The Fourier basis: The $n$-Fourier basis of the structural simplex is constructed by choosing a maximal linearly independent set of eigenvectors of the $n$-Laplace operator that are orthonormal with respect to the inner product $⟨\cdot ,\cdot {⟩}_w$ defined by the weights of the structural simplex (see Eq. (8)). Concretely, one needs to find a matrix $F_n$ such that

  $$F_n{L}_n{F}_n^{-1}=D_n$$

  (11)

  with $D_n$ a diagonal matrix, which also satisfies

  $$(F_n^{-1}{)}^{⊺}{W}_n{F}_n^{-1}=I_n,$$

  where $I_n$ is the identity matrix of dimension $(\frac{N+1}{n+1})$.

• Hyperharmonic representation: As a final step, one calculates the Fourier transform of the high-order signals calculated in Step (3), denoted by ${\widehat{\text{bi}o}}^n$ and ${\widehat{\text{bi}s}}^n$, as follows:

  	$$F_n\text{bi}{o}^n$$	$$={\widehat{\text{bi}o}}^n,$$
  	$$F_n\text{bi}{s}^n$$	$$={\widehat{\text{bi}s}}^n.$$

This workflow has been designed with modularity and flexibility in mind, in order to facilitate its adaptation to the needs of diverse applications. This section highlights possible variations over the workflow that can better accommodate the specific needs of different scenarios.

First, note that our choice of $\Omega$ and $\Sigma$ in Step (3) as high-order information signals is based on their interpretability, and on their promising value for the analysis of complex systems. Nevertheless, other high-order measures (e.g. the ones described in 2.4) can also be analysed following the same pipeline.

Additionally, Step (2) suggests to build the structural simplex by propagating the underlying mutual information from edges to higher-dimensional simplices via averages. However, one could use other constructions: e.g. use the maximum or minimum mutual information instead of the average. Additionally, one could replace the mutual information with other non-negative metrics of similarity, including the total variation distance, the Wasserstein distance, or the absolute value of the Pearson correlation. Furthermore, one could also use non-negative high-order metrics (such as the TC or DTC, see Section 2.3) for building the structural simplex directly.

Finally, it is important to note that Step (1) is included mainly for pedagogical purposes, but is often omitted in practice. In effect, most modern techniques to estimate information-theoretic quantities from data avoid to build an explicit joint distribution, as this introduces additional biases. For discrete data, we recommend considering Bayesian estimators such as the ones discussed in Archer et al. (2013), and the software package DIT (James et al., 2018). For continuous data, our preferred choice are Kraskof estimators (Kraskov et al., 2004) that can be implemented via JIDT (Lizier, 2014). However, these estimators require substantial amounts of data, and a more flexible alternative are provided by estimators based on Gaussian Copulas (Ince et al., 2017).

Our analysis is based on electronic scores that are publicly available at http://kern.ccarh.org, from where we extracted the scores of Symphonies No. 93--94 and No. 99--104.⁴⁴4The data for Symphonies 95–98 was not available.⁴ All symphonies use the same basic instrumentation: flutes, oboes, bassoons, horns, trumpets, timpani, violins, viola, violoncello, and double bass; only some of these symphonies employ clarinets, which were therefore not included in the analysis. To avoid duplicates, our analyses consider only one instrument of each kind, which left an arrangement of nine parts (violoncello and double bass were assumed to be equivalent).

The scores of the four movements of the selected symphonies were pre-processed in Python 3.8.5 using the Music21 package (http://web.mit.edu/music21). Each movement was transformed into nine coupled time series taking 13 possible values — one for each note plus one for the silence, using a small rhythmic duration as time unit. With these data, the joint distribution of the values for the nine-note chords was estimated using their empirical frequency. Note that regularisation methods (such as Laplace smoothing) were not employed, as many configurations (e.g. highly dissonant chords) are never explored in the Classic repertoire.

Our subsequent analyses were restricted to high-order signals depending on the joint probability of no more than $6$ instruments. Given the structure of the dataset (8 symphonies with 4 movements each), we computed the structural simplex using a distribution calculated using all the data. In contrast, individual high-order signals were calculated for each movement of each symphony by using only the corresponding data.

To measure compressibility, we used the following metric defined for any basis. Consider a given $n$-dimensional high-order signal, whose coefficients on the given basis are $\alpha =\left\{{\alpha }_i\right\}$. Without loss of generality we assume that ${\alpha }_i^2\ge {\alpha }_j^2$ if $i\le j$. Then, we define the functions

with ${\text{EV}}_{\alpha }\left(k\right)$ being the (normalised) explained variance by the $k$-th strongest component, and ${\text{CEV}}_{\alpha }\left(k\right)$ the cumulative explained variance by the $k$ strongest components. This definition is motivated by Parseval’s theorem, which guarantees that the sum of the square of the Fourier coefficients is equal to the variance of the signal; hence, the $\text{CEV}\left(k\right)$ of the Fourier transform of a signal corresponds to the percentage of the variance that is accounted by the $k$ strongest components.

The hyperharmonic representation of the considered high-order signals ($\Omega$ and $\Sigma$) was found to be substantially more concentrated than their corresponding canonical representations. Figure 3 illustrates the curves of cumulative explained variance for various dimensions, and Table 1 presents the number of components required to fulfil various reconstruction levels. In particular, it was found that a small number of components suffices to account for most of the variance observed in the hyperharmonic representations of $\Omega$ and $\Sigma$. Moreover, this good performance is found to be stable across dimensions.

To verify the unique properties of the Laplace operators, an additional control was run were the same signals were transformed according to randomly-generated bases. Interestingly, these representations do not exhibit the degree of dimensionality-reduction shown by the hyperharmonic representations (see A).

Finally, our results also show that the canonical representation of the S-information is much less concentrated than the canonical representation of the O-information. This dissimilarity suggests that the O-information is capturing a more specific signature of the different modalities of interdependency that exist across the orchestra. Moreover, the fact that the dissimilarity is not seen in their hyperharmonic representation further suggests that both signals may actually carry an equivalent amount of information, which happens to be differently localised over the canonical basis.