How do kernel-based sensor fusion algorithms behave under high dimensional noise?

In this subsection, we summarize some related results. Since the NCCA and AD matrices (3) are essentially products of transition matrices, we start from summarizing the results of the affinity and transition matrices when there is only one sensor. On one hand, in the noiseless setting, the spectral properties have been widely studied, for example, [2, 21, 22, 44, 18, 11], to name but a few. In summary, under the manifold model, researchers show that the Graph Laplacian(GL) converges to the Laplace–Beltrami operator in various settings with properly chosen bandwidth. On other hand, the spectral properties have been investigated in [4, 3, 8, 14, 13, 17, 28] under the null setup. These works essentially show that when $\mathcal{X}$ contains pure high-dimensional noise, the affinity and transition matrices are governed by a low-rank perturbed Gram matrix when the bandwidth $h_{1}=p_{1}$. Despite rich literature above about two extreme setups, limited results are available in the intermittent, or nonnull, setup [12, 6, 15]. For example, when the signal-to-noise ratio (SNR), which will be defined precisely later, is sufficiently large, the spectral properties of GL constructed from the noisy observation are close to that constructed from the clean signal. Moreover, the bandwidth plays an important role in the nonnull setup. For a more comprehensive review and sophisticated study on the spectral properties of the affinity and transition matrices for an individual point cloud, we refer the readers to [6, Sections 1.2 and 1.3].

For the NCCA and AD matrices, on one hand, in the noiseless setting, there have been several results under the common manifold model [32, 46]. On the other hand, under the null setup that both sensors only capture high dimensional white noise, its spectral property has been studied recently [7]. Specifically, except for a few larger outliers, when $h_{1}=p_{1}$ and $h_{2}=p_{2},$ the edge eigenvalues of $n^{2}\mathbf{A}$ or $n^{2}\mathbf{N}$ converge to some deterministic limit depending on the free convolution (c.f. Definition 2.3) of two Marchenko-Pastur (MP) laws [37]. However, in the nonnull setting when both sensors are contaminated by noise, to our knowledge, there does not exist any theoretical study, particularly under the high dimensional setup.

We now provide an overview of our results. The main contribution of this paper is a comprehensive study of NCCA and AD under the non-null case in the high dimensional setup. This result can be viewed as a continuation of the study under the null case [7]. We focus on the setup that the signal is modeled by a low dimensional manifold. It turns out that this problem can be recast as studying the algorithm under the commonly applied spiked model, which will be made clear later. In addition to providing a theoretical justification based on the kernel random matrix theory, we propose a method to choose the bandwidth adaptively. Moreover, peculiar and counterintuitive results will be presented when two sensors have different behavior, which emphasizes the importance of carefully applying these algorithms in practice. In Section 3, we investigate the eigenvalues of the NCCA and AD matrices when $h_{1}=p_{1}$ and $h_{2}=p_{2}$, which is a common choice in the literature. The behavior of the eigenvalues varies according to both SNRs of the point clouds. When both SNRs are small, the spectral behavior of $n^{2}\mathbf{N}$ and $n^{2}\mathbf{A}$ is like that in the null case, while both the number of outliers and the convergence rates rely on SNRs; see Theorem 3.1 for details. Furthermore, if one of the sensors has large SNR and the other one has small SNR, the eigenvalues of $n\mathbf{N}$ and $n\mathbf{A}$ provide limited information about the signal; see Theorem 3.2 for details. We emphasize that this result warns us that if we directly apply NCCA and AD without any sanity check, it may result in a misleading conclusion. When both SNRs are larger, the eigenvalues are close to the clean NCCA and AD matrices; see Theorem 3.3 for more details. It is clear that the classic bandwidth choices $h_{k}=p_{k}$ for $k=1,2$ are inappropriate when the SNR is large, since the bandwidth is too small compared with the signal strength. In this case $\mathbf{N}\approx\mathbf{I},\mathbf{A}\approx\mathbf{I}$, and we obtain limited information about the signal; see (42) for details. To handle this issue, in Section 4, we consider bandwidths that are adaptively chosen according to the dataset. With this choice, when the SNRs are large, NCCA and AD become non-trivial and informative; that is, NCCA and AD are robust against the high dimensional noise. See Theorem 4.1 for details.

The paper is organized as follows. In Section 2, we introduce the mathematical setup and some background in random matrix theory. In Section 3, we state our main results for the classic choice of bandwidth. In Section 4, we state the main results for the adaptively chosen bandwidth. In Section 5, we offer the technical proofs of the main results. In Appendix 5.1, we provide and prove some preliminary results which will be used in the technical proofs.

We focus on the following model for the datasets $\mathcal{X}$ and $\mathcal{Y}$. Assume that the first sensor i.i.d. sample $n$ clean signals from a sub-Gaussian vector $X:(\Omega,\mathcal{F},\mathbb{P})\to\mathbb{R}^{p_{1}}$, denoted as $\{\mathbf{u}_{ix}\}_{i=1}^{n}\in\mathbb{R}^{p_{1}}$, where $(\Omega,\mathcal{F},\mathbb{P})$ is a probability space. Similarly, assume that the second sensor also i.i.d. sample clean signals from a sub-Gaussian vector $Y$, denoted as $\{\mathbf{u}_{iy}\}_{i=1}^{n}\in\mathbb{R}^{p_{2}}$. Since we focus on the distance of pairwise samples, without loss of generality, we assume that

$$\mathbb{E}\mathbf{u}_{ix}=\mathbf{0}\quad\mbox{and}\quad\mathbb{E}\mathbf{u}_{iy}=\mathbf{0}.$$

(4)

Denote $S_{1}=\operatorname{Cov}(\mathbf{u}_{ix})$ and $S_{2}=\operatorname{Cov}(\mathbf{u}_{iy})$, and to simplify the discussion, we assume that $S_{1}$ and $S_{2}$ admit the following spectral decomposition

$$S_{k}=\operatorname{diag}\{\sigma_{k1}^{2},\cdots,\sigma_{kd}^{2},0,\cdots,0\}\in\mathbb{R}^{p_{k}\times p_{k}},\ k=1,2,$$

(5)

where $d_{1}$ and $d_{2}$ are fixed integers. We model the common information by assuming that there exists a bijection $\phi$ so that

$$X(\Omega)=\phi(Y(\Omega))\,;$$

(6)

that is, we have $\mathbf{u}_{ix}=\phi(\mathbf{u}_{iy})$ for any $i=1,\ldots,n$. In practice, the clean signals $\{\mathbf{u}_{ix}\}$ and $\{\mathbf{u}_{iy}\}$ are contaminated by two sequences of i.i.d. sub-Gaussian noise $\{\mathbf{z}_{i}\}\in\mathbb{R}^{p_{1}}$ and $\{\mathbf{w}_{i}\}\in\mathbb{R}^{p_{2}}$, respectively, so that the data generating process follows

$$\mathbf{x}_{i}=\mathbf{u}_{ix}+\mathbf{z}_{i}\ \ \mbox{and}\ \ \mathbf{y}_{i}=\mathbf{u}_{iy}+\mathbf{w}_{i}\,,$$

(7)

where

$$\mathbb{E}\mathbf{z}_{i}=\mathbf{0}_{p_{1}},\ \operatorname{Cov}(\mathbf{z}_{i})=\mathbf{I}_{p_{1}},\ \mathbb{E}\mathbf{w}_{i}=\mathbf{0}_{p_{2}},\ \operatorname{Cov}(\mathbf{w}_{i})=\mathbf{I}_{p_{2}}.$$

(8)

We further assume that $\mathbf{z}_{i}$ and $\mathbf{w}_{i}$ are independent with each other and also independent of $\{\mathbf{u}_{ix}\}$ and $\{\mathbf{u}_{iy}\}$. We are mainly interested in the high dimensional setting; that is, $p_{1}$ and $p_{2}$ are comparably as large as $n.$ More specifically, we assume that there exists some small constant $0<\gamma<1$ such that

$$\gamma\leq c_{1}:=\frac{n}{p_{1}}\leq\gamma^{-1}\ \ \mbox{and}\ \ \gamma\leq c_{2}:=\frac{n}{p_{2}}\leq\gamma^{-1}.\$$

(9)

The SNRs in our setting are defined as $\{\sigma_{1i}^{2}\}$ and $\{\sigma_{2i}^{2}\},$ respectively, so that for all $1\leq i\leq d_{1}$ and $1\leq j\leq d_{2}$,

$$\sigma_{1i}^{2}\asymp n^{\zeta_{1i}},\ \sigma_{2j}^{2}\asymp n^{\zeta_{2j}}\,,$$

(10)

for some constants $\zeta_{1i},\zeta_{2j}\geq 0$. To avoid repetitions, we summarize the assumptions as follows.

In view of (5), the model (7) for each sensor is related to the spiked covariance matrix models [27]. We comment that this seemingly simple model, particularly (5), includes the commonly considered nonlinear common manifold model. In the literature, the common manifold model means that two sensors sample simultaneously from one low dimensional manifold; that is, $\mathbf{u}_{ix}=\mathbf{u}_{iy}\in M$ and $\phi$ is an identity map, where $M$ is a low dimensional smooth and compact manifold embedded in the high dimensional Euclidean space. Since we are interested in the kernel matrices depending on pairwise distances, which is invariant to rotation, when combined with Nash’s embedding theory, the common manifold can be assumed to be supported in the first few axes of the high dimensional space, like that in (5). As a result, the common manifold model becomes a special case of the model (7). We refer readers to [6] for a detailed discussion of this relationship. A special example of the common manifold model is the widely considered linear subspace as the common component; that is, when $M=\mathbb{R}^{d}$ embedded in $\mathbb{R}^{p_{k}}$ for $k=1,2$. In this case, we could simply apply CCA to estimate the common component, and its behavior in the high dimensional setup has been studied in [1, 36].

We should emphasize that through the analysis of NCCA and AD under the common component model satisfying Assumption (2.1), we do not claim that we could understand the underlying manifold structure. The problem we are asking here is the nontrivial relationship between the noisy and clean affinity and transition matrices, while the problem of exploring the manifold structure from the clean datasets [18, 11] is a different one, which is usually understood as the manifold learning problem. To study the nontrivial relationship between the noisy and clean affinity and transition matrices, it is the spiked covariance structure that we focus on, but not the possibly non-trivial $\phi$. By establishing the nontrivial relationship between the noisy and clean affinity and transition matrices in this paper, when combined with the knowledge of manifold learning via the kernel-based manifold learning algorithm with clean datasets [31, 46, 43], we know how to explore the common manifold structure, which depends on $\phi$, from the noisy datasets.

In this subsection, we introduce some random matrix theory background and necessary notations. Let $\mathbf{Z}\in\mathbb{R}^{p_{1}\times n}$ be the data matrix associated with $\{\mathbf{z}_{i}\}$; that is, the $i$-th column $\mathbf{Z}$ is $\mathbf{z}_{i}$, and consider the scaled noise $s\mathbf{Z}$, where $s>0$ stands for the standard deviation of the scaled noise. Denote the empirical spectral distribution (ESD) of $\mathbf{Q}=\frac{s^{2}}{p_{1}}\mathbf{Z}^{\top}\mathbf{Z}$ as

$$\mu_{\mathbf{Q}}(x)=\frac{1}{n}\sum_{i=1}^{n}\mathbf{1}_{\{\lambda_{i}(\mathbf{Q})\leq x\}},\ x\in\mathbb{R}.$$

It is well-known that in the high dimensional regime (9), $\mu_{\mathbf{Q}}$ has the same asymptotic [29] as the so-called MP law [37], denoted as $\nu_{c_{1},s^{2}}$, satisfying

$$\nu_{c_{1},s^{2}}(I)=(1-c_{1})_{+}\chi_{I}(0)+\zeta_{c_{1},s^{2}}(I)\,,$$

(11)

where $I\subset\mathbb{R}$ is a measurable set, $\chi_{I}$ is the indicator function and $(a)_{+}:=0$ when $a\leq 0$ and $(a)_{+}:=a$ when $a>0$,

$$d\zeta_{c_{1},s^{2}}(x)=\frac{1}{2\pi s^{2}}\frac{\sqrt{(\lambda_{+,1}-x)(x-\lambda_{-,1})}}{c_{n}x}\mathrm{d}x\,,$$

(12)

$\lambda_{+,1}=(1+s^{2}\sqrt{c_{1}})^{2}$ and $\lambda_{-,1}=(1-s^{2}\sqrt{c_{1}})^{2}$. Denote

$$\tau_{1}\equiv\tau_{1}(\lambda_{1}):=2\left(\frac{\lambda_{1}}{p_{1}}+1\right),\ \tau_{2}\equiv\tau_{2}(\lambda_{2}):=2\left(\frac{\lambda_{2}}{p_{2}}+1\right),$$

(13)

and for $k=1,2$,

$$\varsigma_{k}\equiv\varsigma_{k}(\tau_{k}):=1-2\upsilon\exp(-\upsilon\tau_{k})-\exp(-\upsilon\tau_{k})\,.$$

(14)

For any constant $\mathsf{a}>0,$ denote $\mathrm{T}_{\mathsf{a}}$ be the shifting operator that shifts a probability measure $\nu$ defined on $\mathbb{R}$ by $\mathsf{a};$ that is

$$\mathrm{T}_{\mathsf{a}}\nu(I)=\nu(I-\mathsf{a}),$$

(15)

where $I-\mathsf{a}$ means the shifted set. Using the notation (11), for $k=1,2,$ denote

$$\nu_{k}:=\mathrm{T}_{\varsigma_{k}}\nu_{c_{k},\eta},\ \mbox{where }\ \eta=2\upsilon\exp(-2\upsilon).$$

(16)

Next, we introduce a $n$-dependent quantity of some probability measure. For a given probability measure $\mu$ and $n\in\mathbb{N},$ define $\gamma_{\mu}(j)$ as

$$\int_{\gamma_{\mu}(j)}^{\infty}\mu(\mathrm{d}x)=\frac{j}{n}.$$

(17)

Finally, we recall the following notion of stochastic domination [16, Chapter 6.3] that we will frequently use. Let $\mathsf{X}=\big{\{}\mathsf{X}^{(n)}(u):n\in\mathbb{N},\ u\in\mathsf{U}^{(n)}\big{\}}$ and $\mathsf{Y}=\big{\{}\mathsf{Y}^{(n)}(u):n\in\mathbb{N},\ u\in\mathsf{U}^{(n)}\big{\}}$ be two families of nonnegative random variables, where $\mathsf{U}^{(n)}$ is a possibly $n$-dependent parameter set. We say that $\mathsf{X}$ is stochastically dominated by $\mathsf{Y}$, uniformly in the parameter $u$, if for any small $\upsilon>0$ and large $D>0$, there exists $n_{0}(\upsilon,D)\in\mathbb{N}$ so that we have $\sup_{u\in\mathsf{U}^{(n)}}\mathbb{P}\Big{(}\mathsf{X}^{(n)}(u)>n^{\upsilon}\mathsf{Y}^{(n)}(u)\Big{)}\leq n^{-D}$, for a sufficiently large $n\geq n_{0}(\upsilon,D)$. We interchangeably use the notation $\mathsf{X}=\mathrm{O}_{\prec}(\mathsf{Y})$ or $\mathsf{X}\prec\mathsf{Y}$ if $\mathsf{X}$ is stochastically dominated by $\mathsf{Y}$, uniformly in $u$, when there is no danger of confusion. In addition, we say that an $n$-dependent event $\Omega\equiv\Omega(n)$ holds with high probability if for a $D>1$, there exists $n_{0}=n_{0}(D)>0$ so that $\mathbb{P}(\Omega)\geq 1-n^{-D}$, when $n\geq n_{0}.$

In this subsection, we summarize some preliminary results about free multiplication of random matrices from [5, 26]. Given some probability measure $\mu,$ its Stieltjes transform and $M$-transform are defined as

$$m_{\mu}(z)=\int\frac{1}{x-z}\mathrm{d}\mu(x)\ \ \mbox{and}\ \ M_{\mu(z)}=\frac{zm_{\mu}(z)}{1+zm_{\mu}(z)}\,,$$

where $z\in\mathbb{C}\backslash\mathbb{R}_{+}$, respectively. We next introduce the subordination functions utilizing the $M$-transform [26, 47]. For any two probability measures $\mu_{x}$ and $\mu_{y}$, there exist analytic functions $\Omega_{x}(z)$ and $\Omega_{y}(z)$ satisfying

$$zM_{\mu_{x}}(\Omega_{y}(z))=zM_{\mu_{y}}(\Omega_{x}(z))=\Omega_{x}(z)\Omega_{y}(z),\ \text{for}\ z\in\mathbb{C}\backslash\mathbb{R}_{+}.$$

(18)

Armed with the subordination functions, we now introduce the free multiplicative convolution of $\mu_{x}$ and $\mu_{y},$ denoted as $\mu_{x}\boxtimes\mu_{y}$, when $\mu_{x}$ and $\mu_{y}$ are compactly supported on $\mathbb{R}_{+}$ but not both delta measures supported on $\{0\}$; see Definition 2.7 of [5].

For $\nu_{1}$ and $\nu_{2}$ defined in (16), we have two sequences $a_{k}=\gamma_{\nu_{1}}(k-1)$ and $b_{k}=\gamma_{\nu_{2}}(k-1)$, where $1\leq k\leq n$. Note that we have

$$\int_{a_{k}}^{E_{+,1}}\mathrm{d}\nu_{1}(x)=\frac{k-1}{n},\ \int_{b_{k}}^{E_{+,2}}\mathrm{d}\nu_{2}(x)=\frac{k-1}{n},$$

(20)

where $E_{+,1},E_{+,2}$ are the right edges of $\nu_{1}$ and $\nu_{2},$ respectively. Denote two $n\times n$ positive definite matrices $\Sigma_{1}$ and $\Sigma_{2}$ as follows

$$\Sigma_{1}=\operatorname{diag}\{a_{1},a_{2},\cdots,a_{n}\},\ \Sigma_{2}=\operatorname{diag}\{b_{1},b_{2},\cdots,b_{n}\}.$$

(21)

Let $\mathbf{U}$ be an $n\times n$ Haar distributed random matrix in $O(n)$ and denote

$$\mathbf{H}=\Sigma_{2}\mathbf{U}\Sigma_{1}\mathbf{U}^{\top}.$$

The following lemma summarizes the rigidity of eigenvalues of $\mathbf{H}.$

In this subsection, we state the results when at least one sensor contains strong noise, or equivalently has a small SNR; that is, $0\leq\zeta_{2}<1$. In this case, the NCCA and AD will not be able to provide useful information or can only provide limited information for the underlying common manifold.

In this subsection, we state the results when both of the sensors have large SNR ($\zeta_{2}\geq 1$). Recall (29), (30) and denote $\widetilde{\mathbf{A}}_{2,s},\widetilde{\mathbf{A}}_{2,c}$ analogously for the point cloud $\mathcal{Y}.$ For some constant $C>0,$ denote

$$\mathsf{R}:=\begin{cases}C\log n,&\zeta_{2}=1\\ Cn^{\zeta_{2}-1},&1<\zeta_{2}<2.\end{cases}$$

(36)

Theorem 3.3 shows that when $h_{k}=p_{k}$, where $k=1,2,$ and both SNRs are large, the NCCA matrix from the noisy dataset could be well approximated by that from the clean dataset of the common manifold. The main reason has been elaborated in [6] when we have only one sensor. In the two sensors case, combining (37) and (38), we see that except the first $\mathsf{R}$ eigenvalues, the remaining eigenvalues are negligible and not informative. Moreover, (2) and (3) reveal important information about the bandwidth; that is, if the bandwidth choice is improper, like $h_{1}=p_{1}$ and $h_{2}=p_{2}$, the result could be misleading in general. For instance, when $\zeta_{1}$ and $\zeta_{2}$ are large, ideally we should have a “very clean” dataset and we shall expect to obtain useful information about the signal. However, this result says that we cannot obtain any useful information from NCCA; particularly, see (42). This however is intuitively true, since when the bandwidth is too small, the relationship of two distinct points cannot be captured by the kernel; that is, when $i\neq j$, $\exp(-\|\mathbf{x}_{i}-\mathbf{x}_{j}\|^{2}/p_{1})\approx 0$ with high probability (see proof below for a precise statement of this argument or [6]). This problem can be fixed if we choose a proper bandwidth. In Section 4, we will state the corresponding results when the bandwidths are selected properly, in which case this counterintuitive result is eliminated.