Determining the Dimension and Structure of the Subspace Correlated Across Multiple Data Sets

TL;DR

This paper presents a method to fully identify the correlation structure of components across multiple data sets by analyzing eigenvalues and eigenvectors, advancing beyond just estimating the number of correlated components.

Contribution

It introduces a novel approach using eigenanalysis and bootstrap hypothesis testing to determine the complete correlation structure across data sets.

Findings

Eigenvalues and eigenvectors characterize correlation structures

Bootstrap hypothesis testing effectively identifies correlated components

Method improves understanding of relationships among data set components

Abstract

Detecting the components common or correlated across multiple data sets is challenging due to a large number of possible correlation structures among the components. Even more challenging is to determine the precise structure of these correlations. Traditional work has focused on determining only the model order, i.e., the dimension of the correlated subspace, a number that depends on how the model-order problem is defined. Moreover, identifying the model order is often not enough to understand the relationship among the components in different data sets. We aim at solving the complete modelselection problem, i.e., determining which components are correlated across which data sets. We prove that the eigenvalues and eigenvectors of the normalized covariance matrix of the composite data vector, under certain conditions, completely characterize the underlying correlation structure. We use these results to solve the model-selection problem by employing bootstrap-based hypothesis testing.

Tables

Table 1. TABLE I : Example with the correlation structure in Fig. 1 with three data sets each with five signal components. The entries are the correlation coefficients between signal components of different pairs of data sets.

	${\rho }_{12}^{\left(i\right)}$	${\rho }_{13}^{\left(i\right)}$	${\rho }_{23}^{\left(i\right)}$
$i=1$	0.5	0.6	0.6
$i=2$	0.7	0	0
$i=3$	0	0	0.8
$i=4$	0	0.4	0
$i=5$	0	0	0

Table 2. TABLE II : Example of 4 4 4 data sets each with 4 4 4 signal components.

	${\rho }_{12}^{\left(i\right)}$	${\rho }_{13}^{\left(i\right)}$	${\rho }_{14}^{\left(i\right)}$	${\rho }_{23}^{\left(i\right)}$	${\rho }_{24}^{\left(i\right)}$	${\rho }_{34}^{\left(i\right)}$
$i=1$	0	0	0	0.7	0.2	0.8
$i=2$	0	0.6	0.4	0	0	0.5
$i=3$	0.5	0	0	0	0	0
$i=4$	0.5	0	0	0	0	0

Table 3. TABLE III : Correlation structure of the three correlated components in four data sets used in the first simulation setup.

	${\rho }_{12}^{\left(i\right)}$	${\rho }_{13}^{\left(i\right)}$	${\rho }_{14}^{\left(i\right)}$	${\rho }_{23}^{\left(i\right)}$	${\rho }_{24}^{\left(i\right)}$	${\rho }_{34}^{\left(i\right)}$
$i=1$	0.63	0.78	0.69	0.81	0.64	0.91
$i=2$	0.62	0.67	0.74	0.71	0.82	0.91
$i=3$	0.84	0.81	0.72	0.57	0.71	0.62

Equations

$${\mathbf{x}}_p={\mathbf{A}}_p{\mathbf{s}}_p,p=1,2,\dots ,P,$$

$$E[{\mathbf{s}}_p]$$

$${\mathbf{R}}_{s_p{s}_p}$$

$${\mathbf{R}}_{s_p{s}_q}=\text{diag}({\rho }_{pq}^{(1)},{\rho }_{pq}^{(2)},\dots ,{\rho }_{pq}^{(n)}),$$

$$\mathbf{x}=[{\mathbf{x}}_1^T,\dots ,{\mathbf{x}}_P^T{]}^T,$$

$$\mathbf{C}={\mathbf{R}}_D^{-\frac{1}{2}}\mathbf{R}{\mathbf{R}}_D^{-\frac{1}{2}}.$$

$$\mathbf{C}=\left[\begin{array}{cccc}\mathbf{I}& {\mathbf{C}}_{12}& \cdots & {\mathbf{C}}_{1P}\\ {\mathbf{C}}_{21}& \mathbf{I}& \cdots & {\mathbf{C}}_{2P}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathbf{C}}_{P1}& {\mathbf{C}}_{P2}& \cdots & \mathbf{I}\end{array}\right],$$

$${\mathbf{C}}_{pq}$$

$$=({\mathbf{A}}_p{\mathbf{R}}_{s_p{s}_p}{\mathbf{A}}_p^T{)}^{-\frac{1}{2}}{\mathbf{A}}_p{\mathbf{R}}_{s_p{s}_q}{\mathbf{A}}_q^T({\mathbf{A}}_q{\mathbf{R}}_{s_q{s}_q}{\mathbf{A}}_q^T{)}^{-\frac{1}{2}},$$

$$=({\mathbf{A}}_p{\mathbf{A}}_p^T{)}^{-\frac{1}{2}}{\mathbf{A}}_p{\mathbf{R}}_{s_p{s}_q}{\mathbf{A}}_q^T({\mathbf{A}}_q{\mathbf{A}}_q^T{)}^{-\frac{1}{2}},$$

$${\mathbf{C}}_{pq}$$

$$\mathbf{C}=\mathbf{F}{\mathbf{R}}_{ss}{\mathbf{F}}^T,$$

$$\mathbf{C}$$

$$=\mathbf{F}{\mathbf{P}}^T{\tilde{\mathbf{R}}}_{ss}\mathbf{P}{\mathbf{F}}^T,$$

$${\tilde{\mathbf{R}}}_{ss}=\text{blkdiag}({\mathbf{R}}^{(1)},\dots ,{\mathbf{R}}^{(n)}),$$

$${\mathbf{R}}^{(i)}=\left[\begin{array}{cc}{\mathbf{B}}^{(i)}& 0\\ 0& \mathbf{I}\end{array}\right].$$

$${\mathbf{H}}^{(i)}=\left[\begin{array}{cc}0& {\rho }_{pq}^{(i)}\\ {\rho }_{pq}^{(i)}& 0\end{array}\right],$$

$${\mathbf{H}}^{(i)}=\left[\begin{array}{ccc}0& {\rho }_{pq}^{(i)}& {\rho }_{pr}^{(i)}\\ {\rho }_{pq}^{(i)}& 0& {\rho }_{qr}^{(i)}\\ {\rho }_{pr}^{(i)}& {\rho }_{qr}^{(i)}& 0\end{array}\right].$$

$${\mathbf{R}}^{(i)}=\left[\begin{array}{cc}{\mathbf{B}}^{(i)}& 0\\ 0& \mathbf{I}\end{array}\right],$$

$${\mathbf{R}}^{(i)}=\left[\begin{array}{cc}\mathbf{I}+{\mathbf{H}}^{(i)}& 0\\ 0& \mathbf{I}\end{array}\right],$$

$${\mathbf{H}}^{(i)}=\left[\begin{array}{cccc}0& {\rho }^{(i)}& \cdots & {\rho }^{(i)}\\ {\rho }^{(i)}& 0& & ⋮\\ ⋮& & \ddots & ⋮\\ {\rho }^{(i)}& \cdots & \cdots & 0\end{array}\right],$$

$${\mathbf{R}}^{(i)}=\left[\begin{array}{cc}{\mathbf{B}}^{(i)}& 0\\ 0& \mathbf{I}\end{array}\right],$$

$${\tilde{\mathbf{v}}}^{(i)}=\left[\begin{array}{c}{\mathbf{v}}^{(i)}\\ 0\end{array}\right],$$

$${\mathbf{u}}^{(i)}=\mathbf{F}{\mathbf{P}}^T{\tilde{\mathbf{u}}}^{(i)}.$$

$$H_0$$

$$H_1$$

$$H_0:{\lambda }^{(s+1)}={\lambda }^{(s+2)}=\dots ={\lambda }^{(s+P)}=1.$$

$$T(s)=\sum _{i=s+1}^{s+P}({\lambda }^{(i)}-1{)}^2,$$

$$H_0$$

$$H_1$$

$$T=\Vert {\mathbf{u}}_p^{(i)}{\Vert }^2,$$

$$\text{SNR (dB)}=10\log_{10}\bigg{(}\frac{\sigma_{s}^{2}}{\sigma_{n}^{2}}\bigg{)},$$

$$k\le R(\underset{ctimes}{\underbrace{G_{j+1},G_{j+1},\dots ,G_{j+1}}}).$$

Full text

I Introduction

Canonical correlation analysis (CCA), multiset CCA (mCCA), and their variants have been widely used for jointly analyzing the relationships among different sets of data [1, 2, 3, 4, 5]. These tools extract components from each data set, referred to as canonical variables, that are highly correlated across different data sets. The canonical variables can then be used in a plethora of applications including data fusion in biomedicine [6, 7, 8], feature extraction in machine learning [9], finding coupled patterns in oceanography [10], extracting gene clusters in genomics [11], etc. Despite the abundance of interest in extracting correlated components, one key question often remains overlooked. Is the estimated correlation of these components actually present in the underlying system, or is it an artifact due to noise or an insufficient number of samples? Sometimes the answer to this question is assumed to be known a priori from domain-specific knowledge. However, outside of that limited realm, most applications employ methods for thresholding the correlation coefficients as a way to determine the significant correlated components. These solutions are generally heuristic and often fail for one of two main reasons: 1) when the number of samples is limited, the correlation coefficient among the estimated components is overestimated, even to the point of identifying nonexistent, spurious correlations, 2) the number of possible correlation structures among the extracted components combinatorially increases as the number of data sets increases.

Some techniques in the past have aimed to solve this as a model-order selection problem. In signal processing, “model order” is the term used for the dimension of a parameter vector, i.e., the number of parameters of the data model [12]. Thus, estimating the number of correlated components can be posed as a model-order selection problem. For two data sets, the signal components are either correlated or uncorrelated across both sets, and the model-order selection problem is well-defined. In this case, counting the number of correlated components also answers the question of which components are correlated. Some of the techniques for estimating the model order for two data sets are [13, 14, 15]. However, for more than two data sets, the model-order selection problem is ambiguous. It is possible for the components to be correlated across no data sets, all data sets, or some subset of the collection. In this context, many generalizations of model-order selection are valid, and the problem must be more precisely defined. Fig. 1 illustrates an example of correlation structure between the latent signal components of three data sets, ${\bf s}_{1}$, ${\bf s}_{2}$, and ${\bf s}_{3}$. In most of cases, however, we observe linear mixtures of these latent components instead of observing them directly. Nonetheless, Fig. 1 can be used as a reference example to differentiate between the possible definitions of model-order selection for multiple data sets. Each column of the figure indicates the components of one data set and thus components of each column have the same subscript. Each row represents the individual signal components that can be correlated between different pairs of data sets. In this example each data set contains five signal components that are mutually uncorrelated within their set. The set $\mathbf{s}_{1}$ is repeated again in the last column to illustrate the correlation between the first and the third data set. Here, the first signal component of each data set is correlated with the first component of each other data set (the correlation is indicated with red arrows), and the next three components are correlated only between a pair of data sets (indicated with black arrows). The fifth component of each data set is uncorrelated with all other data sets.

One variation of the model-order selection problem for multiple data sets is to determine the number of signal components correlated across every pair of data sets. In Fig. 1, only the first signal components demonstrate this property, so the model order by this definition would be one. Model orders of this type can be identified using the methods described in [16, 17, 18, 19, 20]. The methods in [16] and [17] require the assumption that if signal components are correlated across any data sets then they must be correlated across all data sets, and these methods fail to estimate the model order when this assumption is not met. This assumption is relaxed in [18] and generalized to arbitrary correlation structures in [19] and [20].

Alternatively, signal components correlated across a subset of the collection of data sets might also be of interest. For instance, when tracking an object in videos recorded by spatially separated cameras, the object might not be visible in every frame of each camera [3]. Thus if multiset canonical correlation is used to measure the similarity of frames from the different views, it would be pertinent for the model order to represent the number of components correlated across all data sets or a subset of data sets. Similarly in brain imaging, estimating the number of signal components activated in the fMRI data of all the subjects is useful for multi-subject analysis [6, 21]. However, some brain regions may not appear active for some subjects, due to noise or other factors, even if biological intuition suggests that they should be. Knowing that correlations exist among the signal components corresponding to these regions is useful even if they are only present in a subset of subjects. These scenarios suggest that an appropriate definition of model order should count the signals that demonstrate correlation across all or a subset of data sets. By this definition the model order of the example in Fig. 1 is four. In [3], the authors propose a test statistic to estimate this model order and showed when this number can be correctly estimated for different settings of the signal-to-noise ratio (SNR) and the number of samples.

In the end, determining only the model order is insufficient to completely characterize the correlation structure in multiple data sets. This summary statistic only provides the knowledge that the components exhibit correlation. This knowledge, although sufficient for two data sets, is incomplete for multiple data sets as it is also required to determine which components are correlated across which data sets. For the example in Fig. 1, the complete solution is not just determining that the first four components are correlated but also that the first component in each data set is correlated, and the successive components are correlated between data sets $1$ and $2$, $2$ and $3$, and $1$ and $3$, respectively. Therefore, we formulate and solve a more general model-selection problem, which includes model-order selection as a subproblem. In this work, we propose a technique for solving this model-selection problem for multiple data sets using the eigenvalues and eigenvectors of the coherence matrix (normalized covariance matrix [22]) of the composite data set. The main contributions of this work are:

• •

We prove that, under fairly general conditions, the correlation structure in multiple data sets can be fully characterized from the eigenvector decomposition of the coherence matrix of the composite data.

• •

Using this theoretical result, we develop an algorithm that identifies the correlation structure effectively in practical scenarios.

A competing solution to the model-selection problem was recently presented by [23]. The technique in [23] applies a series of hypothesis tests to pairs of components extracted using CCA. It first determines the number of pairwise correlations, two data sets at a time, and amalgamates these detections using joint information from the complete collection of data sets. The technique proposed in this work has both practical and theoretical benefits over [23] in the right circumstances. First, this work provides conditions under which the correlation structure can be identified, and justifies these conditions theoretically. Second, the proposed method relies solely on joint information from all of the data sets, and as a result, it demonstrates superior detection accuracy for components correlated across more data sets. Finally, for $P$ number of data sets, the proposed solution requires significantly fewer hypothesis tests per component; running only $P$ tests compared to the $2^{P}$ needed for [23].

Our program for this paper is as follows. After defining the problem in Section II, we derive the main result in Section III. In Section IV, we present two hypothesis testing procedures to determine the underlying correlation structure based on the results derived in Section III. We use the bootstrap to estimate the distribution of the test statistic under the null hypothesis. Finally, in Section V simulations show that our technique reliably estimates the complete correlation structure of multiple data sets, and is competitive with the existing state-of-the-art approaches.

II Problem Formulation

We consider $P$ data sets consisting of zero-mean, real-valued random vectors, $\mathbf{x}_{1},\ldots,\mathbf{x}_{P}$, each of dimension $n$. The data sets are generated by an unknown linear mixing of underlying signal vectors, $\mathbf{s}_{1},\ldots,\mathbf{s}_{P}\in\mathbb{R}^{n}$. The generating data model is

[TABLE]

where $\mathbf{A}_{p}\in\mathbb{R}^{{n}\times n}$ is an unknown but fixed mixing matrix with full rank111All the results in this work can be easily extended to the general case where the data sets have different dimensions and the mixing matrices are non-square with full column rank. We omit the general case as the extension is trivial when the inverses are replaced by pseudo-inverses, because the more cumbersome bookkeeping distracts from the actual result. . We make the following additional assumptions:

• •

The signal vectors each contain $n$ signal components. The $i$th signal component of the $p$th data set is denoted by $s_{p}^{(i)}$. These components are assumed to be zero-mean, unit variance, and uncorrelated within each data set, i.e.,

[TABLE]

where $\mathbf{0}$ is a zero-vector, and $\mathbf{I}$ is an identity matrix.

• •

Only signal components with the same component number can have nonzero correlation. As such, the signal cross-covariance matrix between data sets $p$ and $q$ ($p\neq q$) is

[TABLE]

where diag(.) denotes a diagonal matrix with specified diagonal elements, and $\rho_{pq}^{(i)}$ represents the unknown (possibly zero) correlation coefficient between their $i$th components. When analyzing the correlations between two data sets only, ${\bf R}_{s_{p}s_{q}}$ can be assumed to be diagonal without loss of generality. For more than two sets, diagonal cross-covariance matrices are a restriction on the problem, as they do not represent all possible correlation structures. However, this assumption makes the multiset problem tractable and is commonplace in the literature [5, 24].

Our aim is to completely determine the underlying correlation structure among the signal components ${\bf s}_{p}$ based on $M$ independent and identically distributed (i.i.d.) joint samples of the data vectors ${\bf x}_{p}$, $p=1,\ldots,P$. We solve this problem by first estimating the model order $d$, which represents the total number of components that demonstrate nonzero correlation i.e., $d=|\{i|\exists p,q\text{ for which }\rho_{pq}^{(i)}\neq 0\}|$. Then we estimate the data sets across which these $d$ components are correlated, i.e., we identify all data sets $p$, $q$ for which $\rho_{pq}^{(i)}\neq 0$.

To compare to existing techniques, two scalar model orders are also defined. We denote the number of signal components correlated pairwise between data sets $p$ and $q$ as $d_{pq}$. The number of signals correlated across all data sets is denoted by $d_{\textrm{all}}$, i.e., $d_{\textrm{all}}=|\{i|\rho_{pq}^{(i)}\neq 0\text{ }\forall p,q|$. Prior works have either focused on determining $d_{pq}$ for two data sets [15] or on estimating only $d_{\text{all}}$ for more than two data sets [16, 17, 18, 19, 20]. However, as discussed earlier, even knowing both $d_{pq}$ and $d_{\text{all}}$ is not enough to completely determine the underlying correlation structure (except with very special types of correlation structures).

Examples: Let us revisit the example in Fig. 1. Table I provides one example of correlation coefficients that match the structure presented in Fig. 1. The entries are the correlation coefficients between signal components of different pairs of data sets. In this case $d_{pq}=2$ for all choices of $p$ and $q$, $d_{\text{all}}=1$, and $d=4$. Now consider an example of $4$ data sets each with $4$ signal components as shown in Table II. The first component of data sets $2,3,$ and $4$ and the second component of data sets $1,3,$ and $4$ are correlated. The third and fourth components of data sets $1$ and $2$ are correlated. Hence, $d_{\text{all}}=0$, $d=4$ and $d_{pq}$ is the number of nonzero entries in the corresponding column. In both these examples, the techniques in [15, 16, 17, 18, 19, 20] provide solutions for either $d_{pq}$ or $d_{\text{all}}$. Our technique, however, aims to identify which of the entries in Tables I and II are nonzero.

III Correlated Subspace in Multiple Data Sets

The problem as stated in Section II requires joint knowledge of the relationships between all data sets. In this section, our main results demonstrate that the pertinent information can be found in the eigenvector decomposition of the whitened covariance matrix of the concatenation of all data sets. Consider the composite data vector $\mathbf{x}$ obtained by vertically concatenating the individual data vectors,

[TABLE]

and the composite covariance matrix $\mathbf{R}=E[\mathbf{x}\mathbf{x}^{T}]$. Similarly, the composite signal vector is defined as ${\bf s}=[\mathbf{s}_{1}^{T},\ldots,\mathbf{s}_{P}^{T}]^{T}$, and the composite signal covariance matrix as ${\bf R}_{ss}=E[{\bf s}{\bf s}^{T}]$. The definition of the coherence matrix for two data sets [22] can be generalized in a natural way for this composite data as

[TABLE]

Here ${\bf R}_{D}=\text{blkdiag}(\mathbf{R}_{11},\ldots,\mathbf{R}_{PP})$ is a block-diagonal matrix with ${\bf R}_{pp}=E[{\bf x}_{p}{\bf x}_{p}^{T}]$, and exponent ${-\frac{1}{2}}$ on ${\bf R}_{D}$ denotes the inverse of the matrix square root. The composite coherence matrix, $\mathbf{C}$, can be written in a block structure given as

[TABLE]

where the diagonal blocks are identity matrices and the off-diagonal blocks are the coherence matrices of two data sets. For each pair of data sets ${\bf x}_{p}$ and ${\bf x}_{q}$, the coherence matrix can be decomposed as

[TABLE]

since ${\bf R}_{s_{p}s_{p}}={\bf R}_{s_{q}s_{q}}={\bf I}$. Let ${\bf F}_{p}=({{\bf A}}_{p}{{\bf A}}_{p}^{T})^{-1/2}{{\bf A}}_{p}$ and similarly ${{\bf F}}_{q}=({{\bf A}}_{q}{{\bf A}}_{q}^{T})^{-1/2}{{\bf A}}_{q}$ so that ${{\bf F}}_{p}{{\bf F}}_{p}^{T}={{\bf F}}_{q}{{\bf F}}_{q}^{T}={{\bf I}}$, and we have

[TABLE]

Using (7) and (9), the composite coherence matrix ${\bf C}$ can be written as

[TABLE]

where ${\bf F}=\text{blkdiag}({\bf F}_{1},\ldots,{\bf F}_{P})$. Based on the assumption that ${\bf R}_{s_{p}s_{q}}$ is diagonal, the elements of ${\bf R}_{ss}$ can be permuted to form a block-diagonal matrix whose $i$th block is the covariance matrix formed by the $i$th components of each data set. That is, there exists a permutation ${\bf P}$ where

[TABLE]

such that $\widetilde{{\bf R}}_{ss}={\bf P}{\bf R}_{ss}{\bf P}^{T}$ is a block-diagonal matrix defined as

[TABLE]

and ${\bf R}^{(i)}\in\mathbb{R}^{P\times P}$ is the covariance matrix of the $i$th components of each data set.

III-A Eigenvalues of ${\bf C}$

We now illustrate an explicit relationship between the dimension of the correlated subspace, $d$, and the eigenvalues of ${\bf C}$. It can be observed from (7) that ${\bf C}$ is an identity matrix when all signal components are uncorrelated, and thus all eigenvalues of ${\bf C}$ are one. However, when some components are correlated, ${\bf C}$ has eigenvalues that are different from one. More specifically, when data sets demonstrate correlations across $d$ components, ${\bf C}$ has at least $d$ eigenvalues greater than one. A key question then is: when is the dimension of the correlated subspace, $d$, exactly equal to the number of eigenvalues of ${\bf C}$ greater than one?

The answer is not as straightforward as one would hope. However, we can identify a set of sufficient conditions under which this property holds. The property also often holds when these conditions are not met and the algorithm proposed in this work is robust to the violation of these assumptions. The proof of the sufficiency of these conditions relies on three things:

1. i)

The composite coherence matrix ${\bf C}$ is similar (through a similarity transformation) to a block diagonal matrix where each block contains a non-identity matrix corresponding to a collection of data sets whose $i$th components are all correlated with each other and an identity matrix corresponding to data sets whose $i$th components are uncorrelated. 2. ii)

If the $i$th components of four or more data sets are correlated with each other, all nonzero correlations are greater than a prescribed threshold. 3. iii)

All correlations are transitive. That is, if a signal component is correlated between data sets $p$ and $q$, and between data sets $q$ and $r$, then it is also correlated between data sets $p$ and $r$.

For simplicity, we prove the result in Theorem 1 when there is only one block of the form described in i) for the $i$th components of all data sets, but of course the result holds for any number of such blocks. We will discuss the requirement of ii) and iii) during the proof of Theorem 1.

Theorem 1**.**

Let ${\bf C}$ be the composite coherence matrix of $P$ data sets constructed according to the linear mixing model in (1) with pairwise diagonal signal cross-covariance matrices. Let $k^{(i)}$ be the number of data sets whose $i$th components are correlated. Assume that correlations are transitive, and for $k^{(i)}\geq 4$, each correlation coefficient is either $\rho_{pq}^{(i)}=0$ or $\rho_{pq}^{(i)}>\epsilon^{(i)}=(\frac{k^{(i)}-1}{k^{(i)}})^{2}$ for all $p,q$. Let $\mathcal{I}=\{1,\ldots,n\}$ be an index set for the signals. ${\bf C}$ has exactly $d$ eigenvalues greater than one if and only if there exists a subset of signals $\mathcal{D}\subseteq\mathcal{I}$ with $|\mathcal{D}|=d$, and for each $i\in\mathcal{D}$ there exists a $p\neq q$ such that $s_{p}^{(i)}$ and $s_{q}^{(i)}$ are correlated.

Proof.

We begin by showing that if there exists a subset of signals $\mathcal{D}\subseteq\mathcal{I}$ with $|\mathcal{D}|=d,$ and for each $i\in\mathcal{D}$ there exists a $p,q$ with $s_{p}^{(i)}$ and $s_{q}^{(i)}$ correlated, ${\bf C}$ has $d$ eigenvalues greater than one. Let $\widetilde{{\bf R}}_{ss}$ be an $nP\times nP$ matrix with the structure defined in (12). The diagonal blocks of $\widetilde{{\bf R}}_{ss}$ can be indexed by the associated signal component. That is, ${\bf R}^{(i)}\in\mathbb{R}^{P\times P}$ is the covariance matrix of the $i$th components of each data set with the form

[TABLE]

${\bf B}^{(i)}\in\mathbb{R}^{k^{(i)}\times k^{(i)}}$ is a symmetric matrix with diagonal elements equal to one and off-diagonal elements equal to the correlation coefficients between the correlated $i$th components. The indices associated with the subset of data sets whose $i$th components are correlated is $\mathcal{K}^{(i)}\subseteq\mathcal{P}=\{1,\ldots P\}$, and the dimensions of ${\bf B}^{(i)}$ are the size of this subset, $|\mathcal{K}^{(i)}|=k^{(i)}$. As $\widetilde{{\bf R}}_{ss}$ is block-diagonal, its eigenvalues are equal to the eigenvalues of the blocks ${\bf R}^{(i)}$. Since ${\bf F}{\bf P}^{T}$ is an orthogonal matrix, the eigenvalues of ${\bf C}$ are equal to the eigenvalues of $\widetilde{{\bf R}}_{ss}$ and therefore, equal to the eigenvalues of ${\bf R}^{(i)}$.

Let ${\bf B}^{(i)}$ be decomposed as ${\bf B}^{(i)}={\bf I}+{\bf H}^{(i)}$ for each $i$. ${\bf H}^{(i)}$ is a hollow (with zeros on the diagonal) symmetric matrix whose off-diagonal elements are the nonzero correlation coefficients corresponding to the $i$th component. We show that the matrix ${\bf H}^{(i)}$ has exactly one positive eigenvalue for each $i$ as follows:

Case 1 ($k^{(i)}=2$): ${\bf H}^{(i)}$ has exactly one positive eigenvalue for all values of $\rho_{pq}^{(i)}>0$, because

[TABLE]

which has eigenvalues $\{\rho_{pq}^{(i)},-\rho_{pq}^{(i)}\}$.

Case 2 ($k^{(i)}=3$): When the $i$th components of exactly three data sets, $p,q$ and $r$, are correlated, ${\bf H}^{(i)}$ is

[TABLE]

The characteristic polynomial of ${\bf H}^{(i)}$ is $y=-\lambda^{3}+\lambda(\rho_{pq}^{(i)2}+\rho_{pr}^{(i)2}+\rho_{qr}^{(i)2})+2\rho_{pq}^{(i)}\rho_{pr}^{(i)}\rho_{qr}^{(i)}$. Using Descartes’ rule of sign change, $y$ has only one positive root for any $\rho_{pq}^{(i)},\rho_{pr}^{(i)},\rho_{qr}^{(i)}>0$ [25]. Therefore, ${\bf H}^{(i)}$ has only one positive eigenvalue.

Case 3 ($k^{(i)}\geq 4$): [26, Theorem 3.5] can be used to show that ${\bf H}^{(i)}$ has exactly one positive eigenvalue if all the off-diagonal elements of ${\bf H}^{(i)}$ are greater than $\epsilon^{(i)}=(\frac{k^{(i)}-1}{k^{(i)}})^{2}$. This result is demonstrated in Appendix VI. Without the assumption of transitive correlations, ${\bf H}^{(i)}$ cannot be guaranteed to have all the positive off-diagonal elements as required. Thus, ${\bf H}^{(i)}$ has exactly one positive eigenvalue in each case as desired.

Let $\mathcal{D}\subseteq\mathcal{I}$ be the $d$ values of $i$ for which correlation exists. For each $i\in\mathcal{D}$, ${\bf H}^{(i)}$ has one positive eigenvalue and $k^{(i)}-1$ non-positive eigenvalues. Let ${\bf H}^{(i)}={\bf U}\mathbf{\Lambda}^{(i)}{\bf U}^{T}$ be the eigenvalue decomposition (EVD) of ${\bf H}^{(i)}$ with ${\bf U}{\bf U}^{T}={\bf I}$. The EVD of ${\bf B}^{(i)}$ can be written as ${\bf B}^{(i)}={\bf U}({\bf I}+\mathbf{\Lambda}^{(i)}){\bf U}^{T}$ so that, for each $i\in\mathcal{D}$, ${\bf B}^{(i)}$ has one eigenvalue greater than one and the remaining $k^{(i)}-1$ eigenvalues less than or equal to one. Using (13), the maximum eigenvalue of ${\bf B}^{(i)}$ is also the maximum eigenvalue of ${\bf R}^{(i)}$, implying that ${\bf R}^{(i)}$ has exactly one eigenvalue greater than one. Hence ${\bf C}$ has exactly $d$ eigenvalues greater than one as desired.

We now show that the converse is also true. Let ${\bf C}$ has $d$ eigenvalues greater than one. There exists a permutation as described by (III) where

[TABLE]

and ${\bf B}^{(i)}$ is the sum of an identity matrix and a hollow symmetric matrix of correlation coefficients, i.e., ${\bf B}^{(i)}={\bf I}+{\bf H}^{(i)}$, as described above. When no data sets are correlated for a given $i$, ${\bf R}^{(i)}={\bf I}$. Let $\mathcal{D}^{\prime}\subseteq\mathcal{I}$ be the indices for which nonzero correlation exists (and by assumption is greater than $\epsilon^{(i)}$ for $k^{(i)}\geq 4$). Thus,

[TABLE]

for $i\in\mathcal{D}^{\prime}$ and ${\bf R}^{(i)}={\bf I}$ for $i\in\mathcal{I}\textbackslash\mathcal{D}^{\prime}$. Clearly ${\bf R}^{(i)}$ has no eigenvalues greater than one for $i\in\mathcal{I}\textbackslash\mathcal{D}^{\prime}$.

Suppose that $|\mathcal{D}^{\prime}|<d$. Then there exists an $i$ for which ${\bf H}^{(i)}$ has more than one eigenvalue greater than zero. However, this contradicts our proof that ${\bf H}^{(i)}$ has exactly one positive eigenvalue for each $i$. Similarly, suppose that $|\mathcal{D}^{\prime}|>d$. Then there exists an $i$ for which ${\bf H}^{(i)}$ has no eigenvalues greater than zero. This also contradicts our proof, and thus $|\mathcal{D}^{\prime}|=d$ and there are exactly $d$ values of $i$ for which $\rho_{pq}^{(i)}\neq 0$ for some $p\neq q$. ∎

Theorem 1 guarantees that, when $d$ eigenvalues of ${\bf C}$ are greater than one, $d$ signal components are correlated across at least a pair of data sets. Therefore, the number of correlated signal components can be determined by testing for the number of eigenvalues of ${\bf C}$ that are greater than one. Such a test is formulated in Section IV-A using bootstrap-based hypothesis testing.

One of the assumptions in Theorem 1 that needs further discussion is that if the $i$th components of four or more data sets are correlated (i.e., $k^{(i)}\geq 4$), then the correlation coefficient between any pair of $i$th signal components must be either zero or greater than $\epsilon^{(i)}=(\frac{k^{(i)}-1}{k^{(i)}})^{2}$. This assumption guarantees that only one eigenvalue of ${\bf C}$ corresponding to the $i$th component is greater than one. The threshold, $\epsilon^{(i)}$, is derived in Appendix VI and is a restrictive threshold since $\lim_{k^{(i)}\rightarrow\infty}\epsilon^{(i)}=1$. However, the proof does not claim to represent all matrices with the desired eigenvalue structure. That is, there is a nonempty set of real positive hollow symmetric matrices that have exactly one positive eigenvalue but do not meet this element-wise threshold. One example is the following: Suppose the nonzero correlation coefficients associated with the $i$th component are equal for all $k^{(i)}(\geq 4)$ data sets. That is, $\rho_{pq}^{(i)}=\rho^{(i)}$ $\forall p,q\in\mathcal{K}^{(i)}$, where $\mathcal{K}^{(i)}$ is the subset of indices associated with the data sets whose $i$th components are correlated. In this case,

[TABLE]

as defined in the proof of Theorem 1. This can be simplified as ${\bf H}^{(i)}=\rho^{(i)}\mathbf{1}\mathbf{1}^{T}-\rho^{(i)}{\bf I}$, where $\mathbf{1}\in\mathbb{R}^{k^{(i)}}$ is a vector with all elements equal to one. The maximum eigenvalue of ${\bf H}^{(i)}$ is $(k^{(i)}-1)\rho^{(i)}$ and the remaining $k^{(i)}-1$ eigenvalues are $-\rho^{(i)}$. Therefore, ${\bf H}^{(i)}$ has one positive eigenvalue for any $\rho^{(i)}>0$. In this example, the relationship between the eigenvalues of ${\bf C}$ and the number of signals with nonzero correlations described by Theorem 1 holds true for $0<\rho_{pq}^{(i)}\leq 1$.

In the general case, even though $\epsilon^{(i)}$ is restrictive, an element-wise threshold like this is perhaps the best that can be hoped for without imposing further constraints on the structure of the correlation among the components. As noted in [26], for any $k\in\mathbb{N}$ with $k\geq 3$, there exists a positive hollow symmetric ${\bf H}\in\mathbb{R}^{k\times k}$ such that ${\bf H}$ has only two nonpositive eigenvalues. That is to say, without the element-wise constraint there will always be feasible correlation structures for which ${\bf C}$ has more eigenvalues greater than one than signals with nonzero correlations. Moreover, as we will see in our numerical examples later, our hypothesis-test based techniques presented in Section IV may still perform satisfactorily even in cases where the assumptions of Theorem 1 are violated.

As an immediate consequence of Theorem 1, any eigenvalue of ${\bf C}$ that is equal to the maximum possible value, $P$, identifies a signal component where all $P$ data sets are perfectly correlated. This is shown in Corollary 1.

Corollary 1**.**

If any eigenvalue of ${\bf C}$, $\lambda^{(i)}$, is equal to $P$, there exists an $i\in\mathcal{D}$ such that the correlation between $s_{p}^{(i)}$ and $s_{q}^{(i)}$ is one for all $p,q=1,\ldots,P$.

Proof.

By Theorem 1, if $\lambda^{(i)}>1$ is an eigenvalue of ${\bf C}$, there exists an ${\bf R}^{(i)}\neq{\bf I}$ whose largest eigenvalue is equal to $\lambda^{(i)}$. The diagonal elements of $\mathbf{R}^{(i)}$ are equal to one by definition, so $\text{trace}({\bf R}^{(i)})=P=\sum_{j=1}^{P}\lambda_{j}^{(i)}$, where $\lambda_{j}^{(i)}$ is the $j$th largest eigenvalue of ${\bf R}^{(i)}$. Since the largest eigenvalue of ${\bf R}^{(i)}$, $\lambda_{1}^{(i)}=\lambda^{(i)}=P$, we have $\lambda^{(i)}_{2}=\ldots=\lambda^{(i)}_{P}=0,$ and thus, the rank of ${\bf R}^{(i)}$ is one.

Let ${\bf R}^{(i)}={\bf w}{\bf w}^{T}$ be a rank-one matrix, where ${\bf w}=[w_{1},\ldots,w_{P}]^{T}\in\mathbb{R}^{P}$. The diagonal elements of ${\bf R}^{(i)}$ are equal to one, implying that $w_{p}^{2}=1$ for all $p$. Since the off-diagonal elements are bounded by zero and one, $w_{p}$ is positive for $p=1,\ldots,P$. Thus, $w_{1}=w_{2}=\ldots=w_{P}=1$ is the only solution for ${\bf w},$ and ${\bf R}^{(i)}=\mathbf{1}\mathbf{1}^{T}$. Therefore, $\rho^{(i)}_{pq}=1\forall p,q$ and the $i$th component of each data set is perfectly correlated with the $i$th component of all other data sets. ∎

III-B Eigenvectors of ${\bf C}$

The eigenvalues of ${\bf C}$ provide information about the dimension of the correlated subspace, but identifying exactly which data sets demonstrate correlation in a particular component requires more information than this summary contains. The eigenvectors of ${\bf C}$, on the other hand, contain as their elements the coefficients for constructing the correlated signals from each data set that correspond to the associated eigenvalue (that is greater than one). Data sets connected to the nonzero elements of an eigenvector are then the ones whose components are correlated among the associated group of components.

Theorem 2**.**

Let ${\bf C}$ be the composite coherence matrix of $P$ data sets constructed according to the linear mixing model in (1) with pairwise diagonal signal cross-covariance matrices. Let $k^{(i)}$ be the number of data sets whose $i$th components are correlated. Assume that correlations are transitive, and for $k^{(i)}\geq 4$, each correlation coefficient is either $\rho_{pq}^{(i)}=0$ or $\rho_{pq}^{(i)}>\epsilon^{(i)}=(\frac{k^{(i)}-1}{k^{(i)}})^{2}$ for all $p,q$. Let ${\bf C}{\bf u}^{(i)}=\lambda^{(i)}{\bf u}^{(i)}$ such that $\lambda^{(i)}>1$ is an eigenvalue with algebraic multiplicity of one, and let the eigenvector ${\bf u}^{(i)}$ be partitioned into $P$ subvectors, ${\bf u}^{(i)}=[{\bf u}_{1}^{(i)T},{\bf u}_{2}^{(i)T},...,{\bf u}_{P}^{(i)T}]^{T},$ where ${\bf u}_{p}^{(i)}\in\mathbb{R}^{n}$ contains the elements of ${\bf u}^{(i)}$ associated with the $p$th data set. Then the $i$th signal component in the $p$th data set is among the group of correlated $i$th components if and only if ${\bf u}_{p}^{(i)}\neq\mathbf{0}$.

Proof.

Let the $i$th components of $k^{(i)}$ data sets be correlated, and let $\mathcal{K}^{(i)}\subseteq\mathcal{P}=\{1,\ldots,P\}$ be the subset of indices associated with these correlated data sets so that $|\mathcal{K}^{(i)}|=k^{(i)}$. Suppose the $i$th signal component in the $p$th data set is among the correlated components, i.e., $p$ is an element of $\mathcal{K}^{(i)}$. Then there exists a permutation as described by (III) where ${\bf R}^{(i)}\in\mathbb{R}^{P\times P}\neq{\bf I}$ such that

[TABLE]

where ${\bf B}^{(i)}\in\mathbb{R}^{k^{(i)}\times k^{(i)}}$ is a symmetric and element-wise positive matrix that contains the correlation coefficients between the correlated $i$th components. According to Theorem 1, the largest eigenvalue of ${\bf B}^{(i)}$ is greater than one and satisfies ${\bf B}^{(i)}{\bf v}^{(i)}=\lambda_{\text{max}}^{(i)}{\bf v}^{(i)},$ where ${\bf v}^{(i)}$ is the eigenvector of ${\bf B}^{(i)}$ associated with $\lambda_{\text{max}}^{(i)}$. According to the Perron-Frobenius theorem, all the entries of ${\bf v}^{(i)}$ are positive [27]. Since ${\bf R}^{(i)}$ is block-diagonal with blocks ${\bf B}^{(i)}$ and ${\bf I}$, $\lambda_{\text{max}}^{(i)}$ is also the largest eigenvalue of ${\bf R}^{(i)}$ with the eigenvector

[TABLE]

where $\mathbf{0}$ is a zero vector of dimensions $P-k^{(i)}$.

${\bf R}^{(i)}$ is a block of $\widetilde{{\bf R}}_{ss}$ as defined in (12), and thus $\lambda_{\max}^{(i)}$ is also an eigenvalue of $\widetilde{{\bf R}}_{ss}$. Each eigenvector, $\tilde{{\bf u}}^{(i)}$, of $\widetilde{{\bf R}}_{ss}$ can be partitioned into $n$ subvectors where the elements of the $i$th subvector correspond to the $i$th block of $\widetilde{{\bf R}}_{ss}$. Since $\lambda_{\text{max}}^{(i)}$ has an algebraic multiplicity of one, the eigenvector of $\widetilde{{\bf R}}_{ss}$ associated with $\lambda_{\text{max}}^{(i)}$ has $\widetilde{{\bf v}}^{(i)}$ at the $i$th position and zeros everywhere else and can be written as $\tilde{{\bf u}}^{(i)}=[\mathbf{0},\ldots,\tilde{{\bf v}}^{(i)T},\ldots,\mathbf{0}]^{T}$. Using (III), the eigenvector of ${\bf C}$ associated with $\lambda_{\text{max}}^{(i)}$ is related to $\tilde{{\bf u}}^{(i)}$ by

[TABLE]

For each $p\in\mathcal{K}^{(i)}$, the element of $\tilde{{\bf v}}^{(i)}$ corresponding to the $p$th data set is strictly positive. Therefore, from the definition of $\tilde{{\bf u}}^{(i)}$ and (21), ${\mathbf{u}}^{(i)}_{p}\neq\mathbf{0}$ as desired.

To show the other implication, suppose that the $i$th signal component in the $p$th data set is uncorrelated among the $i$th group of components. Thus, $p$ is not an element of $\mathcal{K}^{(i)}$. Therefore, the $p$th data set is not represented in ${\bf v}^{(i)}$ but rather corresponds to one of the elements in the zero vector of $\tilde{{\bf v}}^{(i)}$. Using the definition of $\tilde{{\bf u}}^{(i)}$ and (21), it is easy to see that the $p$th part of ${\bf u}^{(i)},{\bf u}_{p}^{(i)}=\mathbf{0}$. This implies that if the $p$th part of ${\bf u}^{(i)}$ associated with $\lambda^{(i)}>1$, ${\bf u}_{p}^{(i)}\neq\mathbf{0}$, then the $i$th signal component in the $p$th data set is among the group of correlated $i$th components as desired. ∎

Theorem 2 assumes that any eigenvalue of ${\bf C}$ that is greater than one has an algebraic multiplicity of one. This is not a necessary but a sufficient condition. Appendix -B discusses its sufficiency and also scenarios in which the correlation structure can be completely determined using Theorem 2 even when this assumption is not true.

As a result of Theorems 1 and 2, if the $i$th eigenvalue of ${\bf C}$ is greater than one and is unique as an eigenvalue, the existence of correlation associated with the $i$th component of the $p$th data set can be determined by testing the hypothesis $\mathbf{u}^{(i)}_{p}=\mathbf{0}$. A bootstrap-based hypothesis test for this purpose is proposed in Section IV-B.

IV Bootstrap-based Tests for Eigenvalues and Eigenvectors of ${\bf C}$

IV-A Test for detecting the eigenvalues of ${\bf C}$ representing the correlated subspace

Theorem 1 allows us to infer that the number of eigenvalues of ${\bf C}$ greater than one is equal to the dimension of the correlated subspace, $d$. In practice, however, the composite coherence matrix, ${\bf C}$, is unknown and has to be estimated from the samples. As a result, the number of eigenvalues of the sample composite coherence matrix that are greater than one will often not equal the dimension of the correlated subspace, $d$. This inconsistency is addressed in the related model-order selection literature by setting a threshold for the eigenvalues that is determined with an information-theoretic criterion (ITC) or hypothesis testing. The algorithm presented here fits in the category of hypothesis testing.

A popular approach to solve the problem of model-order selection using hypothesis testing is to perform a series of binary hypothesis tests and stop when a certain condition is met. In this context, this means starting with a counter $s=0$ and performing the following binary test of null hypothesis $H_{0}$ and alternative $H_{1}$

[TABLE]

If $H_{0}$ is rejected, $s$ is incremented and another test of $H_{0}$ vs. $H_{1}$ is run. This is repeated until $H_{0}$ is retained or $s$ reaches its maximum possible value.

The binary test in (IV-A) requires a statistic whose (asymptotic) distribution under $H_{0}$ is theoretically known or estimated from the samples. The distribution of the statistic can be derived assuming a particular correlation structure among the components [17], but for arbitrary correlation structures the distribution is unknown. In this work, we use the bootstrap technique to estimate this distribution.

Let the eigenvalues of ${\bf C}$ be arranged with ordering ${\lambda}^{(1)}\geq{\lambda}^{(2)}\ldots\geq{\lambda}^{(nP)}$. To estimate $d$, we propose a statistic that is based on testing whether the eigenvalues of ${\bf C}$ following its $d$ largest eigenvalues are equal to one. If a signal component in one data set is not correlated with the signal components of any other data set, at least one eigenvalue of C is one. In order to increase the power and stability of the test, we make the stronger assumption that each data set has at least one signal component that is completely uncorrelated in this manner, i.e., $d_{pq}<n$ for all $p$,$q$. Therefore, ${\bf C}$ has at least $P$ eigenvalues equal to one and the null hypothesis for each test is

[TABLE]

Note that we cannot include all $nP-d$ eigenvalues following the $d$ largest eigenvalues of ${\bf C}$ since some of them are less than one. The null hypothesis is rejected when the proposed test statistic,

[TABLE]

is sufficiently far from zero. The distribution of $T(s)$ under $H_{0}$ is estimated using the bootstrap [28].

In Algorithm 1, we present an algorithm for estimating $d$ using $T{(s)}$ given $M$ joint samples of data sets ${\bf x}_{1},\ldots,{\bf x}_{P}$. The $M$ samples of each data set form the columns of the data matrices $\mathbf{X}_{1},\ldots,\mathbf{X}_{P}$. The $P$ data matrices along with the number of bootstrap resamples, $B$, and the probability of false alarm, $P_{\text{fa}}$, are the inputs of Algorithm 1.

IV-B Test for eigenvectors of ${\bf C}$ corresponding to correlated components

In addition to identifying the dimension of the correlated subspace, $d$, our stated goal is to estimate the structure of the correlations between the collection of data sets. As a result of Theorem 2, we need to identify the values of $i$ and $p$ for which the subvector ${\bf u}_{p}^{(i)}=\mathbf{0}$ in order to determine which data sets have an uncorrelated $i$th signal component. However, we still do not have direct access to the composite coherence matrix. When ${\bf C}$ is estimated from samples, these subvectors will not be exactly zero. Thus we propose a novel method for identifying multiset correlation structure that uses a bootstrap-based test to detect zero subvectors in the eigenvectors of ${\bf C}$.

Assuming ${d}$ correlated components (which can be estimated via Algorithm 1), the technique tests the ${d}$ eigenvectors associated with the $d$ largest eigenvalues of the sample composite coherence matrix. For $i=1\ldots{d}$ and $p=1\ldots P$ we test the hypotheses

[TABLE]

The null hypothesis is rejected when the squared Euclidean norm of $\mathbf{u}_{p}^{(i)}$,

[TABLE]

is sufficiently far from zero. The distribution of the statistic $T$ under the null hypothesis is estimated using a bootstrap akin to the one in Algorithm 1. The full procedure is given in Algorithm 2.

Combining Algorithms 1 and 2 leads to an effective method for determining the correlation structure among multiple data sets. Algorithm 1 determines how many signal components, ${d}$, have nonzero correlations, and Algorithm 2 reveals the data sets across which these ${d}$ components are correlated. The final output is a binary matrix, $\mathbf{Z}$, that is similar to Tables I and II except that nonzero correlation coefficients are represented by ones. We refer to this as a correlation map, an example of which can be seen in Fig. 6.

V Numerical Results

In this section, we use Monte Carlo simulations to demonstrate the performance of the proposed technique that combines Algorithms 1 and 2. Initially, we compare our technique with those reported by [16, 17, 19, 20] which aim to estimate the number of components correlated across all data sets, $d_{\text{all}}$. To estimate $d_{\text{all}}$ for our technique, we run Algorithms 1 and 2 and then count the number of components that are correlated across all the data sets. Next, we investigate the behavior of the proposed technique when the pairwise correlation coefficients are not above the threshold necessary for the proof of Theorem 1. We show that the method is robust to the violation of this assumption and that the accuracy remains high for many correlation structures. Finally, we compare the performance of our technique with an approach in [23], which also determines the complete correlation structure in multiple data sets. This comparison highlights the quantitative and qualitative differences between the two techniques.

We present results with four different correlation structures. Some simulation settings are common to all four. The signal components in each data set are Gaussian distributed with variance of 1. The mixing matrices are randomly generated orthogonal matrices. Each data set is corrupted by additive white Gaussian noise (AWGN)222The Gaussian distribution of signal components and noise, and the orthogonality of mixing matrices is not a requirement for the technique proposed in this work. MATLAB code for the proposed method and competing approaches can be found at

https://github.com/SSTGroup/Correlation-Analysis-in-High-Dimensional-Data and can be tested with different settings from the ones mentioned here.. The variance of noise components is chosen according to the signal-to-noise-ratio (SNR) which is defined per component as

[TABLE]

where $\sigma_{s}^{2}$ is the variance of the signal and $\sigma_{n}^{2}$ is the variance of the noise. The SNR is the same for all data sets. The number of bootstrap resamples is $B=1000$ and the probability of false alarm is $P_{\text{fa}}=0.05$. The performance plots are shown as a function of SNR, which is varied from $-10$ dB to $15$ dB. The results are averaged over $500$ independent trials. The performance of each method for determining model order is measured by the mean accuracy (number of correct estimates divided by number of trials). The performance in estimating the complete correlation structure is measured using precision, i.e., the percentage of correctly detected correlations among all the detected correlations, and recall, i.e., the percentage of correctly detected correlations among all actual correlations. The four different scenarios are the following:

1. i)

Evaluation of model-order selection with special correlation structure, for $P=4$ data sets with $d=d_{\text{all}}=3$: Each data set is of dimension $n=7$ and the number of samples is $M=350$. The components are either correlated across all data sets or are uncorrelated. Thus, the number of components that are correlated between at least a pair of data sets, $d$, is equal to the number of components correlated across all data sets, $d_{\text{all}}$. This type of correlation structure satisfies the assumption in [16, 17]. The pairwise correlation coefficients for the three correlated components are shown in Table III, all of which exceed the $\epsilon^{(i)}=(3/4)^{2}$ threshold as assumed by Theorem 1.

Fig. 2 shows the mean accuracy of $\hat{d}_{\text{all}}$ as a function of the SNR for the proposed and the competing techniques. All the techniques correctly estimate $d_{\text{all}}$ when the SNR is high. When the SNR is low, the proposed approach outperforms all competing techniques.

2. ii)

Evaluation of model-order selection with arbitrary correlation structure, for $P=4$ data sets with $d=3,d_{\text{all}}=1$: In this setting all the parameters are the same as in the previous scenario except that the first component of each data set is correlated with all other data sets but the other two components are only correlated across a subset of data sets. The second component is correlated across all except the first data set and the third component is correlated between data sets two and four. The methods in [16, 17] are not evaluated as they are inapplicable in this setting.

Fig. 3 shows that the proposed technique works better than the techniques in [19, 20, 23] in estimating the model order $d_{\text{all}}$ for low values of SNR. It is also worth noting that the techniques in [19, 20] estimate only $d_{\text{all}}$ while the proposed method and [23] also detect the components correlated across subsets of the data sets along with their correlation structure.

3. iii)

Performance of the proposed method when the element-wise threshold is not met, for $P=5$ data sets with $d=d_{\text{all}}=2$: We also investigate the performance of the proposed technique for determining the number of correlated components when some of the pairwise correlation coefficients do not meet the threshold required for Theorem 1. For this, we assume that the first two components in each data set are correlated across all data sets. Thus, the threshold for the pairwise correlation coefficients given by Theorem 1 is $\epsilon^{(1)}=\epsilon^{(2)}=\epsilon=(\frac{4}{5})^{2}=0.64$. We keep some of the pairwise correlation coefficients above the threshold $\epsilon$ and vary the remaining ones. More precisely, we set $\rho^{(i)}_{pq}=0.7$ for $p>q$, $p=2,3,4,5$ and $q=3,4,5$, which exceeds the threshold $\epsilon$, and jointly vary the remaining correlation coefficients $\rho_{12}^{(i)}=\rho_{13}^{(i)}=\rho_{14}^{(i)}=\rho_{15}^{(i)}=\rho$, for $i=1,2$. To demonstrate the relative robustness of the method against violating the assumption in Theorem 1, we show the accuracy of $\hat{d}$ as a function of $\rho$ in Fig. 4 for different values of SNR.

For $\rho<\epsilon$, the performance depends on the SNR: For low SNR, it becomes increasingly difficult to correctly estimate $d$ for only weakly correlated components. On the other hand, as long as the SNR is high enough, violating the threshold in Theorem 1 does not present a problem, and $d$ can still be correctly determined.

4. iv)

Evaluation of complete correlation structure, for $P=5$ data sets with $d=3$, $d_{\text{all}}=1$: Finally, we compare the performance of the proposed method and the technique in [23] for determining the complete correlation structure. Each data set has $n=4$ components and the number of samples is $M=250$. Of the $d=3$ correlated components, $d_{\text{all}}=1$ and the first component of each data set is correlated with the first component of each other data set. The second component of each data set is correlated across all the data sets except the fourth data set. Finally, the third components of data sets one, four and five are correlated. Each pairwise correlation coefficient is $0.7$, thus exceeding the threshold as required in Theorem 1.

Fig. 5(a) shows the mean accuracy of estimating ${d}$ for both techniques as a function of SNR. Both methods perform well for medium and high SNR values when estimating $d$, with the proposed method slightly outperforming the method in [23] for low values of SNR. Fig. 5(b) shows the mean precision and recall values for determining the complete correlation structure. The precision of the method in [23] is better at low SNR while the recall for the proposed method is better at both low and medium SNRs.

The main reason for this difference is tied to the information that each method uses for hypothesis testing. The method in [23] performs mCCA to extract the sets of signal components that are highly correlated via a deflationary approach. The components are extracted jointly from all data sets. Then, hypothesis tests are conducted on the extracted signal components from pairs of data sets to detect the underlying correlation structure. At low SNR values some of the correlations are missed while testing an individual pair of data sets, so its average recall is small in this regime.

On the other hand, the proposed method detects the components and their correlation structure by applying hypothesis tests to the eigenvalues and eigenvectors of the composite coherence matrix directly. The eigenvalues and the eigenvectors are functions of all the pairwise correlation coefficients associated with the component and thus, the proposed method tests on this joint information.

This is further illustrated in Fig. 7, which shows heat maps of the average accuracy for the two methods in estimating the complete correlation structure for this scenario. The true correlation structure is visualized in Fig. 6 using a binary map. This map mirrors the structure of Tables I, II and III but represents the nonzero correlation coefficients with white blocks and the zero correlation coefficients with black blocks. This binary map can be compared to heat maps of the simulation results to assess qualitative differences between the proposed method and the technique in [23]. Ideally, these heat maps should look exactly like the binary map in Fig. 6.

In the low SNR regime (shown in Fig. 7a and b at SNR = $-7$dB), both techniques detect very few correlations as illustrated by the dark color of the boxes. Some of the boxes in the second row corresponding to $i=2$ in Fig. 7a that should be black are dark brown, indicating that the proposed method detects a few false positives, and therefore has low precision value compared to the method in [23].

The main differences between performance of the two methods start to appear from SNR of $-4$dB. Fig. 7c,d and 7e,f show the heat maps for SNR$=-4$dB and SNR$=-1$dB, respectively. In Fig. 7d and 7f, boxes corresponding to the nonzero correlations of same component number have different colors. This is expected because [23] conducts tests on pairs of data sets so it is possible to detect the correlation between one pair and miss it between another. On the other hand, the proposed technique tests the eigenvalues and eigenvectors of the composite coherence matrix, which provides a summary of all the pairwise correlations. Therefore, the boxes of the nonzero correlations corresponding to components with same component number are more uniform in color in Fig. 7c and 7e.

At low SNR values, the number of data sets correlated across a particular component affects the accuracy of the proposed method. The more data sets that are correlated along a given component, the better the proposed method performs. This can be observed in Fig. 7c. The boxes of the first row ($i=1$) corresponding to nonzero correlations are significantly brighter than those of second and third rows ($i=2,3$), indicating a higher accuracy when detecting the first component. Similarly in Fig. 7e, there is a contrast between the high accuracy when detecting the correlations of first and second components and the relatively lower accuracy of detecting the third component. This is because the eigenvalue associated with the component correlated across more data sets is significantly greater than one and thus makes its detection possible even when the noise power is high. This stands in contrast with the method from [23] where no advantage is gained by the number of data sets across which a component is correlated. In Fig. 7d and 7f we can see that boxes corresponding to the nonzero correlations of $i=1,2,3$ have less variation in color, indicating similar accuracy across the board.

VI Conclusion

We have presented a technique that solves a model-selection problem to determine the complete correlation structure in multiple sets of data, i.e., identifying which components are correlated across which data sets. We proved that the eigenvalues and the eigenvectors of the composite coherence matrix completely characterize this information. Using these results, we developed a bootstrap-based hypothesis testing technique that estimates the complete correlation structure, which has shown competitive performance and broad applicability in various simulation scenarios. The proposed technique utilizes the joint information from all the data sets provided by the eigenvalues and the eigenvectors and thus performs well in detecting the components that are correlated across a large number of data sets even at low SNR. This is promising for various fields, e.g., biomedicine, genomics and remote sensing, where a common phenomenon is observed using different modalities.

-A Number of positive eigenvalues of a hollow symmetric matrix

Crucial to the goal of identifying correlated signals via the eigenvalues of ${\bf C}$ is a requirement that the correlated subspace of signal components associated with the same component number correspond to a single eigenvalue. As we see in Theorem 1, this requirement necessitates identifying classes of matrices that have exactly one positive eigenvalue. The most general result in this domain comes from [26] and characterizes the eigenvalues of hollow (zero-diagonal) symmetric nonnegative matrices.

To leverage this result, we must first define a generalized Ramsey number. Let $\{G_{a_{1}},G_{a_{2}},\ldots,G_{a_{c}}\}$ be a collection of simple graphs where $a_{i}$ is the number of vertices of the $i$th graph. Suppose we wish to color the edges of a complete graph $G$ with $c$ colors. The generalized Ramsey number, $R(G_{a_{1}},G_{a_{2}},\ldots,G_{a_{c}})$, is the minimum number of vertices of the complete graph $G$ such that for any $c$-coloring of $G$, there exists an $i\in\{1,2,\ldots,c\}$ such that $G_{a_{i}}$ is an induced subgraph of $G$ with all edges of color $i$. By Ramsey’s Theorem [29] such a number always exists.

Theorem 3**.**

(Charles, Farber, Johnsons, and Kennedy-Schaffer [[26], Theorem 3.5]).* Let $k$ and $2\leq j\leq k-1$ be two positive integers, let $G_{j+1}$ be a complete graph on $j+1$ vertices, and let $c>1$ be the smallest integer for which*

[TABLE]

Let $\epsilon=(\frac{j}{j+1})^{c}.$ Then all hollow symmetric nonnegative matrices of order at least $k$ and with off-diagonal entries from $(\epsilon,1]$ have at least $j$ nonpositive eigenvalues.

As an immediate consequence of this theorem we have the following relevant result.

Corollary 2**.**

Let $\epsilon=(\frac{k-1}{k})^{2}$. If ${\bf H}\in\mathbb{R}^{k\times k}$ is a hollow symmetric matrix with off-diagonal elements from $(\epsilon,1],$ then the largest eigenvalue of ${\bf H}$ is positive and the other $k-1$ eigenvalues are nonpositive.

Proof.

In the special case of $j=k-1$, we are trying find the smallest value of $c>1$ for which the Ramsey number of $c$ copies of $G_{k}$ is greater than $k$, that is, the smallest value of $c$ for which

[TABLE]

For $c=2$, the Ramsey number of $\{G_{k},G_{k}\}$ is the minimum number of vertices needed for a complete graph such that any coloring results in an isomophic copy of $G_{k}$ whose edges are all monochromatic (one color). When $k=2$, any 2-coloring of $G_{2}$ contains a monochromatic copy of $G_{2}$, therefore $k=R(G_{2},G_{2})$. For $k>2$, there exists a 2-coloring of $G_{k}$ that does not contain a monochromatic isomorphic copy of $G_{k}$. For example, any 2-coloring that is not monochromatic will not contain a copy of $G_{k}$. Thus the number of vertices needed must be greater than $k$ and $k\leq R(G_{k},G_{k})\forall k$. This implies that ${\bf H}$ will have at least $k-1$ nonpositive eigenvalues when the off-diagonal elements are chosen from the interval $(\epsilon,1]$ with $\epsilon=(\frac{k}{k+1})^{2}$. Since the trace of a symmetric matrix is the sum of its eigenvalues, the eigenvalues of ${\bf H}$ must sum to zero. This means that $k-1$ eigenvalues are nonpositive, but not all identically zero, implying that the largest eigenvalue must be positive. Thus ${\bf H}$ has exactly one positive eigenvalue as desired. ∎

-B Identifiability of the underlying correlation structure

Theorem 1 and 2 state the conditions that allow us to determine the correlated components along with their correlation structure using the eigenvalue decomposition of ${\bf C}$. One additional assumption in Theorem 2 is that the eigenvalues of ${\bf C}$ greater than one are distinct, i.e., have algebraic multiplicity of one. In this section, we will briefly discuss why this assumption is needed. We will also mention the scenarios in which the correlation structure can still be completely determined using Theorem 2 even if the assumption is not true.

Let $\lambda^{(i)}$ and $\lambda^{(j)}$ be the two eigenvalues of ${\bf C}$ with $\lambda^{(i)}>1$ and $\lambda^{(j)}>1$. Let ${\bf u}^{(i)}$ and ${\bf u}^{(j)}$ be the eigenvectors associated with $\lambda^{(i)}$ and $\lambda^{(j)}$, respectively. Let ${\bf u}=a{\bf u}^{(i)}+b{\bf u}^{(j)}$ be a vector formed by a linear combination of ${\bf u}^{(i)}$ and ${\bf u}^{(j)}$, and $a$ and $b$ are scalars. If $\lambda^{(i)}=\lambda^{(j)}$, any linear combination of ${\bf u}^{(i)}$ and ${\bf u}^{(j)}$ is an eigenvector of $\lambda^{(i)}$ or $\lambda^{(j)}$. In this case, if the $i$th and $j$th group of components are correlated across different data sets, their correlation structure, i.e., across which data sets the components are correlated, cannot always be determined using Theorem 2. For instance, if the $i$th components are correlated across all data sets except the $p$th data set, then according to Theorem 2, the $p$th part of ${\bf u}^{(i)}$, ${\bf u}_{p}^{(i)}=\mathbf{0}$. Similarly, ${\bf u}_{q}^{(j)}=\mathbf{0}$ if the $j$th components are correlated across all data sets except the $q$th data set. When $\lambda^{(i)}=\lambda^{(j)}$, then $a{\bf u}^{(i)}+b{\bf u}^{(j)}$ can also be an eigenvector of $\lambda^{(i)}$ or $\lambda^{(j)}$ for any $a,b$. Therefore, ${\bf u}_{p}^{(i)}$ or ${\bf u}_{q}^{(j)}$ are not necessarily equal to zero.

However, if the $i$th and $j$th components are correlated across the same subset of data sets, even when $\lambda^{(i)}=\lambda^{(j)}$, their correlation structure can be determined using Theorem 2. This is due to the fact that the zeros in ${\bf u}=a{\bf u}^{(i)}+b{\bf u}^{(j)}$ are at the same positions as those of ${\bf u}^{(i)}$ and ${\bf u}^{(j)}$ for any $a,b$.

To conclude, Theorem 2 can completely identify the correlation structure of the components when the eigenvalues associated with the components that are correlated across different subset of data sets are distinct.