Influence Function and Robust Variant of Kernel Canonical Correlation Analysis

TL;DR

This paper develops a robust kernel canonical correlation analysis (CCA) method using influence functions and generalized loss functions, improving robustness to contaminated data and outliers in unsupervised kernel learning.

Contribution

It introduces a robust kernel covariance and cross-covariance operator, derives their influence functions, and proposes a robust kernel CCA method that outperforms standard kernel CCA on contaminated data.

Findings

Influence function of standard kernel CCA can identify outliers.

Robust kernel CCA performs better than standard kernel CCA on noisy data.

Proposed methods are applicable to other kernel-based unsupervised learning techniques.

Abstract

Many unsupervised kernel methods rely on the estimation of the kernel covariance operator (kernel CO) or kernel cross-covariance operator (kernel CCO). Both kernel CO and kernel CCO are sensitive to contaminated data, even when bounded positive definite kernels are used. To the best of our knowledge, there are few well-founded robust kernel methods for statistical unsupervised learning. In addition, while the influence function (IF) of an estimator can characterize its robustness, asymptotic properties and standard error, the IF of a standard kernel canonical correlation analysis (standard kernel CCA) has not been derived yet. To fill this gap, we first propose a robust kernel covariance operator (robust kernel CO) and a robust kernel cross-covariance operator (robust kernel CCO) based on a generalized loss function instead of the quadratic loss function. Second, we derive the IF for robust kernel CCO and standard kernel CCA. Using the IF of the standard kernel CCA, we can detect influential observations from two sets of data. Finally, we propose a method based on the robust kernel CO and the robust kernel CCO, called {\bf robust kernel CCA}, which is less sensitive to noise than the standard kernel CCA. The introduced principles can also be applied to many other kernel methods involving kernel CO or kernel CCO. Our experiments on synthesized data and imaging genetics analysis demonstrate that the proposed IF of standard kernel CCA can identify outliers. It is also seen that the proposed robust kernel CCA method performs better for ideal and contaminated data than the standard kernel CCA.

Tables

Table 1. Table 1: Mean and standard deviation of the measure of kernel covariance operator (Standard) and robust kernel covariance operator (Robust).

	Data	TCSD		SFSD
Measure	Kernel	Standard	Robust	Standard	Robust
	Poly-1	$0.9874\pm 0.0017$	$0.8963\pm 0.0069$	$0.3067\pm 0.1026$	$0.1669\pm 0.0626$
	Poly-2	$1.0000\pm 0.0000$	$0.9863\pm 0.0020$	$0.9559\pm 0.0622$	$0.5917\pm 0.1598$
${\Vert {\hat{\Sigma }}_{XX}\Vert }_F$	Poly-3	$1.0000\pm 0.0000$	$0.9996\pm 0.0001$	$0.9973\pm 0.0094$	$0.8793\pm 0.1067$
	Gaussian	$0.1153\pm 0.0034$	$0.1181\pm 0.0039$	$0.1174\pm 0.0266$	$0.1059\pm 0.0258$
	Laplacian	$0.1420\pm 0.0032$	$0.1392\pm 0.0035$	$0.1351\pm 0.0459$	$0.1280\pm 0.0366$
	Poly-1	$0.9993\pm 0.0001$	$0.9940\pm 0.0005$	$0.8074\pm 0.0838$	$0.6944\pm 0.1118$
	Poly-2	$1.0000\pm 0.0000$	$0.9996\pm 0.0001$	$0.9921\pm 0.0122$	$0.9070\pm 0.0703$
${\Vert {\hat{\Sigma }}_{XX}\Vert }_M$	Poly-3	$1.0000\pm 0.0000$	$1.0000\pm 0.0000$	$0.9994\pm 0.0020$	$0.9709\pm 0.0344$
	Gaussian	$0.1300\pm 0.0133$	$0.1028\pm 0.0038$	$0.1065\pm 0.0583$	$0.0735\pm 0.0370$
	Laplacian	$0.1877\pm 0.0053$	$0.1474\pm 0.0042$	$0.1065\pm 0.0583$	$0.0735\pm 0.0370$

Table 2. Table 2: Mean and standard deviation of the measures, η ρ subscript 𝜂 𝜌 \eta_{\rho} and η f subscript 𝜂 𝑓 \eta_{f} of the standard kernel CCA (Standard) and robust kernel CCA (Robust).

	Measure	${\eta }_{\rho }$		${\eta }_f$
Data	n	$Standard$	$Robust$	$Standard$	$Robust$
	$100$	$1.9114\pm 3.5945$	$1.2445\pm 3.1262$	$1.3379\pm 3.5092$	$1.3043\pm 2.1842$
MGSD	$500$	$1.1365\pm 1.9545$	$1.0864\pm 1.5963$	$0.8631\pm 1.3324$	$0.7096\pm 0.7463$
	$1000$	$1.1695\pm 1.6264$	$1.0831\pm 1.8842$	$0.6193\pm 0.7838$	$0.5886\pm 0.6212$
	$100$	$0.4945\pm 0.5750$	$0.3963\pm 0.4642$	$1.6855\pm 2.1862$	$0.9953\pm 1.3497$
SCFSD	$500$	$0.2581\pm 0.2101$	$0.2786\pm 0.4315$	$1.3933\pm 1.9546$	$1.1606\pm 1.3400$
	$1000$	$0.1537\pm 0.1272$	$0.1501\pm 0.1252$	$1.6822\pm 2.2284$	$1.2715\pm 1.7100$
	$100$	$0.6455\pm 0.0532$	$0.1485\pm 0.1020$	$0.6507\pm 0.2589$	$2.6174\pm 3.3295$
SMSD	$500$	$0.6449\pm 0.0223$	$0.0551\pm 0.0463$	$3.7345\pm 2.2394$	$1.3733\pm 1.3765$
	$1000$	$0.6425\pm 0.0134$	$0.0350\pm 0.0312$	$7.7497\pm 1.2857$	$0.3811\pm 0.3846$

Table 3. Table 3: Mean and standard deviation of the differences between the training and test correlation in 10 10 10 -fold cross-validation using standard kernel CCA (Standard) and robust kernel CCA (Robust), respectively.

Data		$\mathrm{Standard}$	$\mathrm{Robust}$
MGSD	ID	$0.4151\pm 0.210$	$0.3119\pm 0.09140$
	CD	$0.3673\pm 0.1196$	$0.2609\pm 0.09660$
SCFSD	ID	$0.0002\pm \mathrm{0.0.0001}$	$0.0002\pm 0.0001$
	CD	$0.0003\pm 0.0002$	$0.0002\pm 0.0001$

Table 4. Table 4: Mean and standard deviation of the differences between taring and test correlation of 10 fold cross-validation of MICI and TCGA data sets using standard kernel CCA (Standard) and robust kernel CCA (Robust).

Data			$\mathrm{Standard}$	$\mathrm{Robust}$
	SNP & fMRI	All	$0.8107\pm 0.1782$	$0.7867\pm 0.13012$
		Without outliers	$0.7361\pm 0.1494$	$0.7348\pm 0.1299$
MCIC	SNP & Methylation	All	$0.7337\pm 0.2000$	$0.7639\pm 0.1433$
		Without outliers	$0.6606\pm 0.1772$	$0.7852\pm 0.1776$
	fMRI &Methylation	All	$0.8424\pm 0.1803$	$0.7842\pm 0.1219$
		Without outliers	$0.7479\pm 0.1671$	$0.7617\pm 0.1759$
	OVSC	All	$0.1679\pm 0.0611$	$0.1792\pm 0.0631$
TCGA		Without outliers	$0.2335\pm 0.0482$	$0.1976\pm 0.0525$
	LISC	All	$0.0779\pm 0.0411$	$0.0713\pm 0.0369$
		Without outliers	$0.1349\pm 0.0547$	$0.0987\pm 0.0597$

Table 5. Table 5 : List of abbreviations

Abbreviation	Elaboration	Abbreviation	Elaboration
CC	Canonical correlation	CCA	Canonical correlation Analysis
CO	Kernel covariance operator	CCO	Kernel cross-covariance operator
CV	Canonical Variates	DE	Density estimation
DNA	Deoxyribonucleic acid	EIF	Empirical influence function
IF	Influence function	fMRI	Functional magnetic resonance imaging
IRWLS	Iteratively re-weighted least squares	KIRWLS	Kerneled iteratively re-weighted least squares
LUSC	lung squamous cell carcinoma	MCIC	The mind clinical imaging consortium
ME	Mean element	MGSD	Multivariate Gaussian structural data
NIH	National institutes of health	NSF	National Science Foundation
OVSC	Ovarian serous cystadenocarcinoma	PCA	Principal component analysis
PDK	Positive definite kernel	RKHS	Reproducing kernel Hilbert Space
SCFSD	Sine cosine function structural data	SFSD	Sine function of structural data
SNP	Single-nucleotide polymorphism	SZ	Schizophrenia
TCSD	Three circles structural data	TCGA	The cancer genome atlas

Table 6. Table 6 : List of symbols.

Symbol	Explanation	Symbol	Explanation
$ℝ$	The set of real numbers	$\Omega$	A sample space
$𝒜$	A set of events	$P$	A function from events to probabilities
$𝚽(\cdot )$	Feature map	$\widetilde{𝚽}(\cdot )$	Centered feature map
${\mathscr{H}}_X$	RKHS of X data	${\mathscr{H}}_X$	RKHS of Y data
$\mathscr{H}$	Hilbert space	${ℬ}_{\mathscr{H}}$	$\sigma$-filed of Borel sets
$k_X(X_i,X_j)$	PDK of X data	$k_Y(Y_i,Y_j)$	PDK of Y data
${\tilde{k}}_X(X_i,X_j)$	Centered PDK of X data	${\tilde{k}}_Y(Y_i,Y_j)$	Centered PDK of Y data
$k_X(\cdot ,X)$	${\mathscr{H}}_X$-valued random variable	$k_Y(\cdot ,Y)$	${\mathscr{H}}_Y$-valued random variable
$F_X$	Probability distribution of X	$F_Y$	Probability distribution of Y
$F_{XY}$	Joint probability distribution of (X,Y)	$F_{nXY}$	Empirical joint probability distribution.
$E_X$	Expectation of X	$E_Y$	Expectation of Y
${ℳ}_X$	Kernel mean element of X	${ℳ}_Y$	Kernel mean element of Y
${\widehat{ℳ}}_X$	Estimated kernel mean element of X	${\widehat{ℳ}}_R$	Robust kernel mean element
$\zeta (\cdot )$	Robust loss function	$\phi \left(t\right)$	$\frac{{\zeta }^{\prime }\left(t\right)}{t}$, Weight function
${𝐊}_X$	Gram matrix of X data	${𝐊}_Y$	Gram matrix of Y data
${𝐆}_X$	Centered Gram matrix of X data	${𝐆}_Y$	Centered Gram matrix of Y data
$\rho$	Robust kernel canonical correlation	$\hat{\rho }$	Estimate robust kernel canonical correlation
${𝐚}_X$	Canonical direction of X data	${𝐚}_Y$	Canonical direction of Y data
${𝐆}_{RX}$	Robust centered Gram matrix of X data	${𝐆}_{RY}$	Robust centered Gram matrix of Y data
${𝐛}_X$	Robust canonical direction of X data	${𝐛}_Y$	Robust canonical direction of Y data
${\rho }_{rkcc}$	Kernel canonical correlation	${\hat{\rho }}_{rkcc}$	Estimate kernel canonical correlation
${\Sigma }_{XX}$	Kernel CO	${\Sigma }_{XY}$	Kernel CCO
${\Sigma }_{RXX}$	Robust kernel CO	${\Sigma }_{RXY}$	Robust kernel CCO
${\hat{\Sigma }}_{XX}$	Estimate of the kernel CO	${\hat{\Sigma }}_{XY}$	Estimate of the kernel CCO
${\hat{\Sigma }}_{RXX}$	Estimate of the robust kernel CO	${\hat{\Sigma }}_{RXY}$	Estimate of the robust kernel CCO
$f_X$	Canonical projection along eigen-function of X data	$f_Y$	Canonical projection along eigen-function of Y data
$g_X$	Robust canonical projection along eigen-function of X data	$g_Y$	Robust canonical projection along eigen-function of Y data
${\eta }_{KCO}$	Performance measure on the kernel CO	${\eta }_{RKCO}$	Robust performance measure of the kernel CO
${\eta }_{\rho }$	Performance measure on the kernel CC	${\eta }_f$	Performance measure of the kernel CV

Equations

$${\mathcal{M}}_X(\cdot )={\mathbb{E}}_X[k_X(\cdot ,X)].$$

$$⟨{\mathcal{M}}_X,f{⟩}_{{\mathcal{H}}_X}=⟨{\mathbb{E}}_X[k_X(\cdot ,X)],f{⟩}_{{\mathcal{H}}_X}={\mathbb{E}}_X⟨k_X(\cdot ,X),f{⟩}_{{\mathcal{H}}_X}={\mathbb{E}}_X[f(X)],$$

$$⟨{\hat{\mathcal{M}}}_X,f{⟩}_{{\mathcal{H}}_X}=⟨\frac{1}{n}\sum _{i=1}^n{k}_X(\cdot ,X_i),f⟩=\frac{1}{n}\sum _{i=1}^{n}f(X_i).$$

$${\hat{\mathcal{M}}}_X=\underset{f\in {\mathcal{H}}_X}{\mathrm{argmin}}\sum _{i=1}^n\Vert \mathbf{\Phi }(X_i)-f{\Vert }_{{\mathcal{H}}_X}^2.$$

$${\hat{\mathcal{M}}}_R=\underset{f\in {\mathcal{H}}_X}{\mathrm{argmin}}\sum _{i=1}^n\zeta (\Vert \mathbf{\Phi }(X_i)-f{\Vert }_{{\mathcal{H}}_X}).$$

$$\zeta (t)=\{\begin{array}{l}{t}^2/2,0\le t\le c\\ ct-c^2/2,c\le t,\end{array}$$

$$\zeta (t)=\{\begin{array}{l}{t}^2/2,0\le t\le c_1\\ c_{1}t-c_1^2/2,c_1\le t<c_2\\ -\frac{c_1}{2(c_3-c_2)}(t-c_3{)}^2+\frac{c_1(c_2+c_3-c_1)}{2},c_2\le t<c_3\\ \frac{c_1(c_2+c_3-c_1)}{2},c_3\le t,\end{array}$$

$$\zeta (t)=\{\begin{array}{l}1-(1-(t/c{)}^2{)}^3,0\le t\le c\\ 1,c\le t,\end{array}$$

$${\hat{\mathcal{M}}}_R^{(h)}=\sum _{i=1}^n{w}_i^{(h-1)}{k}_X(\cdot ,X_i),$$

$${\tilde{K}}_{ij}^{test}=({\mathbf{K}}^{test}-1_t{\mathbf{w}}^T\mathbf{K}-{\mathbf{K}}^{test}\mathbf{w}{1}_n^T+1_t{\mathbf{w}}^T\mathbf{Kw}{1}_n^T{)}_{ij}$$

$${\Sigma }_{XY}=E_{XY}[\tilde{\mathbf{\Phi }}(X)\otimes \tilde{\mathbf{\Phi }}(Y)],$$

$$⟨f_X,{\Sigma }_{XY}{f}_Y{⟩}_{{\mathcal{H}}_X}$$

$${\hat{\Sigma }}_{XY}=\underset{\Sigma \in {\mathcal{H}}_X\otimes {\mathcal{H}}_Y}{\mathrm{argmin}}\sum _{i=1}^n\Vert \tilde{\mathbf{\Phi }}(X_i)\otimes \tilde{\mathbf{\Phi }}(Y_i)-\Sigma {\Vert }_{{\mathcal{H}}_X\otimes {\mathcal{H}}_Y}^2.$$

$${\hat{\Sigma }}_{XY}$$

$${\hat{\Sigma }}_{RXY}=\underset{\Sigma \in {\mathcal{H}}_X\otimes {\mathcal{H}}_Y}{\mathrm{argmin}}\sum _{i=1}^n\zeta (\Vert \tilde{\mathbf{\Phi }}(X_i)\otimes \tilde{\mathbf{\Phi }}(Y_i)-\Sigma {\Vert }_{{\mathcal{H}}_X\otimes {\mathcal{H}}_Y}).$$

$$J(\Sigma )=\sum _{i=1}^n\zeta (\Vert \tilde{\mathbf{\Phi }}(X_i)\otimes \tilde{\mathbf{\Phi }}(Y_i)-\Sigma {\Vert }_{{\mathcal{H}}_X\otimes {\mathcal{H}}_Y}).$$

$$\partial A(Z;\mathcal{T})=\underset{ϵ\to 0}{\lim }\frac{A(Z+ϵ\mathcal{T})+A(Z)}{ϵ}.$$

$$\delta J(\Sigma ,\mathcal{T})=-⟨S(\Sigma ),\mathcal{T}{⟩}_{{\mathcal{H}}_X\otimes {\mathcal{H}}_Y},$$

$$S(\Sigma )=\sum _{i=1}^n\phi (\Vert \tilde{\mathbf{\Phi }}(X_i)\otimes \tilde{\mathbf{\Phi }}(Y_i)-\Sigma {\Vert }_{{\mathcal{H}}_X\otimes {\mathcal{H}}_Y})\cdot (\tilde{\mathbf{\Phi }}(X_i)\otimes \tilde{\mathbf{\Phi }}(Y_i)-\Sigma ).$$

$${\hat{\Sigma }}_{RXY}(X,Y)=\sum _{i=1}^n{w}_i\tilde{k}(X,X_i)\tilde{k}(Y,Y_i)$$

$$w_i\propto \phi (\Vert \tilde{\mathbf{\Phi }}(X_i)\otimes \tilde{\mathbf{\Phi }}(Y_i)-{\hat{\Sigma }}_{RXY}{\Vert }_{{\mathcal{H}}_X\otimes {\mathcal{H}}_Y}).$$

$${\Sigma }^{(h)}=\sum _{i=1}^n{w}_i^{(h-1)}\tilde{\mathbf{\Phi }}(X_i)\otimes \tilde{\mathbf{\Phi }}(Y_i),$$

$$U=\{\Sigma \in {\mathcal{H}}_X\otimes {\mathcal{H}}_Y|S(\Sigma )=0\}$$

$$\Vert {\Sigma }^{(h)}-U{\Vert }_{{\mathcal{H}}_X\otimes {\mathcal{H}}_Y}\doteq \underset{\Sigma {\mathcal{H}}_X\otimes {\mathcal{H}}_Y}{\inf }\Vert {\Sigma }^{(h)}-\Sigma {\Vert }_{{\mathcal{H}}_X\otimes {\mathcal{H}}_Y}\to 0$$

$$T_n=T(F_n),$$

$$T_F^{\prime }(G-F)=\underset{ϵ\to 0}{\lim }\frac{T((1-ϵ)F+ϵG)-T(F)}{ϵ}.$$

$$IF(X^{\prime },T,F)=\underset{ϵ\to 0}{\lim }\frac{T(F^{ϵ})-T(F)}{ϵ}$$

$${\gamma }^{\ast }=su{p}_{X\in X}|IF(X,F,R)|.$$

$${\lambda }^{\ast }=su{p}_{X_1\ne X_2}\frac{|IF(X_1,F,T)-IF(X_2;F,T)|}{|X_1-X_2|}.$$

$${\rho }_{\mathcal{F}}^{\ast }=inf\{m>0;IF(X,F,T)=0,\mathrm{when}|X|>\mathrm{m}\}.$$

Full text

††journal: Journal of LaTeX Templates

1 Introduction

To accelerate the analysis of complex data, kernel based methods (i.e., the support vector machine, kernel ridge regression, multiple kernel learning, kernel dimension reduction in regression, and so on) have proved to be powerful techniques and have been actively studied over the last two decades due to their many flexibilities [1, 2, 3, 4, 5, 6]. Examples of unsupervised kernel methods include kernel principal component analysis (kernel PCA), kernel canonical correlation analysis (standard kernel CCA), and weighted multiple kernel CCA [7, 8, 9, 10, 11]. These methods have been extensively studied for decades in the use of unsupervised kernel methods. However, all of these approaches are not robust and are sensitive to the contaminated model. This paper introduces the robust kernel covariance operator (kernel CO) and kernel cross-covariance operator (kernel CCO) for unsupervised kernel methods such as kernel CCA.

Although many researchers have been studying the robustness issue in a supervised learning setting (e.g., the support vector machine for classification and regression [12, 13, 14]) there are generally few well-founded robust methods for kernel unsupervised learning. The robustness is an important and challenging issue in using statistical machine learning for multiple source data analysis. This is because outliers often occur in real data, which can wreak havoc when used in statistical machine learning methods. Since 1960s, many robust methods, which are less sensitive to outliers, have been developed to overcome this problem. The objective of robust statistics is to use the methods from the bulk of the data and detect the deviations from the original patterns [15, 16].

Recently, in the field of kernel methods, a robust kernel density estimator (robust kernel DE) based on robust kernel mean elements (robust kernel ME) has been proposed by [17], which is less sensitive to outliers than the kernel density estimator. Robust kernel DE is computed using a kernelized iteratively re-weighted least squares (KIRWLS) algorithm in a reproducing kernel Hilbert space (RKHS). In addition, two spatial robust kernel PCA methods have been proposed based on the weighted eigenvalue decomposition [18] and spherical kernel PCA [19], showing that the influence function (IF) of kernel PCA, a well-known measure of robustness, can be arbitrarily large for unbounded kernels.

The kernel methods explicitly or implicitly depend on the kernel CO or the kernel CCO. These operators are among the most useful tools in unsupervised kernel methods but have not yet been robustified. This paper shows that they can be formulated as an empirical optimization problem to achieve robustness by combining empirical optimization problems with the idea of Huber or Hampel on the M-estimation model [15, 16]. The robust kernel CO and robust kernel CCO can be computed efficiently via a KIRWLS algorithm.

In the past decade, CCA with a positive definite kernel has been proposed and is called standard kernel CCA. Several of its variants have also been proposed [20, 21, 22, 23]. Due to the use of simple eigen decomposition, they are still a well-used method for multiple source data analysis. An empirical comparison and sensitivity analysis for robust linear CCA and standard kernel CCA have also been discussed, which give a similar interpretation as kernel PCA but without any robustness measure (e.g., IF of standard kernel CCA) [24]. In addition, the author in [25] has proposed the IF of canonical correlation and canonical vectors of linear CCA. While the IF of an estimator can characterize its robustness, asymptotic properties and standard error, the IF of standard kernel CCA has not yet been proposed. In addition, a robust kernel CCA has not yet been studied. All of these considerations provide motivation to study the IF of kernel CCA and the robust kernel CCA in unsupervised learning.

The contribution of this paper is fourfold. First, we propose a robust kernel CO and robust kernel CCO based on a generalized loss function instead of the quadratic loss function. Second, we propose the IF of standard kernel CCA: kernel canonical correlation (kernel CC) and kernel canonical variates (kernel CV). Third, we propose a method for detecting the influential observations from multiple sets of data, by proposing a visualization method using the IF of kernel CCA. Finally, we propose a method based on robust kernel CO and robust kernel CCO, called robust kernel CCA, which is less sensitive than standard kernel CCA. Experiments on both synthesized data and imaging genetics analysis demonstrate that the proposed visualization and robust kernel CCA can be applied effectively to ideal and contaminated data.

The remainder of this paper is organized as follows. In the following section, we provide a brief review of positive definite kernel, kernel ME, robust kernel ME and kernel CCO. In Section 3 we present the definition, representer theorem, KIRWLS convergence, and a algorithm of robust kernel CCO. In Section 4, we discuss the basic notion of the IF, the IF of kernel ME, kernel CO, kernel CCO and robust kernel CCO. After a brief review of standard kernel CCA in Section 5.1, we propose the IF of standard kernel CCA (kernel CC and kernel CV) and the robust kernel CCA in Section 5.2 and in Section 5.3, respectively. In Section 6, we describe experiments conducted on both synthesized data and real imaging genetics analysis. In Section 7, concluding remarks and future research directions are presented. In the appendix, we discuss the detailed results.

2 Standard and robust kernel (cross-) covariance operator

The kernel ME, kernel CO, and kernel CCO with positive definite kernel have been extensively applied to nonparametric statistical inference through representing distributions in the form of means and covariance in the RKHS [26, 27, 28, 17, 29]. To define the kernel ME, robust kernel ME, kernel CO and kernel CCO, we need the basic notions of positive definite kernels and Reproducing kernel Hilbert space (RKHS), which are briefly addressed in the following [30, 31, 32].

2.1 Basic notion of kernel methods

Let $F_{X}$, $F_{Y}$ and $F_{XY}$ be probability measures on the given nonempty sets $\mathcal{X}$, $\mathcal{Y}$ and $\mathcal{X}\times\mathcal{Y}$, respectively, such that $F_{X}$ and $F_{Y}$ are the marginals of $F_{XY}$. Also let $X_{1},X_{2},\ldots,X_{n}$; $Y_{1},Y_{2},\ldots,Y_{n}$ and $(X_{1},Y_{1}),(X_{2},Y_{2}),\ldots,(X_{n},Y_{n})$ be the independent and identically distributed (IID) samples from the distribution $F_{X}$, $F_{Y}$ and $F_{XY}$, respectively. A symmetric kernel, $k(\cdot,\cdot)\colon\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$, defined on a space is called a positive definite kernel if the Gram matrix $(k(X_{i},X_{j}))_{ij}$ is positive semi-definite for all $i,j\in\{1,2,\cdots,n\}$. A RKHS is a Hilbert space with a reproducing kernel whose span is dense in the Hilbert space. We can equivalently define an RKHS as a Hilbert space of functions with all evaluation functionals bounded and linear. The Moore-Aronszajn theorem states that every symmetric, positive definite kernel defines a unique reproducing kernel Hilbert space [30]. The feature map is a mapping $\mathbf{\Phi}:x\to\mathcal{H}_{X}$ and defined as $\mathbf{\Phi}(\cdot)=k(\cdot,x),\forall\,x\in\mathcal{X}$). The vector $\mathbf{\Phi}(x)\in\mathcal{H}_{X}$ is called a feature vector. The inner product of two feature vectors can be defined as $\langle\mathbf{\Phi}(x),\mathbf{\Phi}(x^{\prime})\rangle_{\mathcal{H}_{X}}=k(x,x^{\prime})$ for all $x,x^{\prime}\in\mathcal{X}$. This is called the kernel trick. By the reproducing property, $f(x)=\langle f(\cdot),k(\cdot,x)\rangle_{\mathcal{H}_{X}}$, with $f\in\mathcal{H}_{X}$ and the kernel trick, the kernel can evaluate the inner product of any two feature vectors efficiently, without knowing an explicit form of either the feature map or the feature vector. Another great advantage is that the computational cost does not depend on the dimension of the original space after computing the Gram matrices [33, 10].

2.2 Standard kernel mean element

Let $k_{X}$ be a measurable positive definite kernel on $\mathcal{X}$ with $\mathbb{E}_{X}[\sqrt{k(X,X)}]<\infty$. The kernel mean, $\mathcal{M}_{X}$, of $X$ on $\mathcal{H}_{X}$ is an element of $\mathcal{H}_{X}$ and is defined by the mean of the $\mathcal{H}_{X}$-valued random variable $k_{X}(\cdot,X)$,

[TABLE]

The kernel mean always exists with arbitrary probability under the assumption that positive definite kernels are bounded and measurable. By the reproducing property, the kernel ME satisfies the following equality

[TABLE]

for all $f\in\mathcal{H}_{X}$.

The empirical kernel ME, $\widehat{\mathcal{M}}_{X}=\frac{1}{n}\sum_{i=1}^{n}\mathbf{\Phi}(X_{i})=\frac{1}{n}\sum_{i=1}^{n}k_{X}(\cdot,X_{i})$ is an element of the RKHS,

[TABLE]

The empirical kernel ME of the feature vectors $\mathbf{\Phi}(X_{i})$ can be regarded as a solution to the empirical risk optimization problem [17]

[TABLE]

2.3 Robust kernel mean element

As explained in Section 2.2, the kernel ME is the solution to the empirical risk optimization problem, which is a least square type of estimator. This type of estimator is sensitive to the presence of outliers in the feature, $\mathbf{\Phi}(X)$. To reduce the effect of outliers, we can use $M$-estimation. In recent years, the robust kernel ME has been proposed for density estimation [17]. The robust kernel ME, based on a robust loss function $\zeta(t)$ on $t\geq 0$, is defined as

[TABLE]

Examples of robust loss functions include Huber’s loss function, Hampel’s loss function, or Tukey’s biweight loss function. Unlike the quadratic loss function, the derivative of these loss functions are bounded [15, 34, 35]. The Huber’s function, a hybrid approach between squared and absolute error losses, is defined as:

[TABLE]

where c ($c>0$) is a tuning parameter. The Hampel’s loss function is defined as:

[TABLE]

where the non-negative free parameters $c_{1}<c_{2}<c_{3}$ allow us to control the degree of suppression of large errors. The Tukey’s biweight loss functions is defined as:

[TABLE]

where $c>0$.

The basic assumptions of the loss functions are; (i) $\zeta$ is non-decreasing, $\zeta(0)=0$ and $\zeta(t)/t\to 0$ as $t\to 0$, (ii) $\varphi(t)=\frac{\zeta^{\prime}(t)}{t}$ exists and is finite, where $\zeta^{\prime}(t)$ is the derivative of $\zeta(t)$, (iii) $\zeta^{\prime}(t)$ and $\varphi(t)$ are continuous, and bounded, and (iv) $\varphi(t)$ is Lipschitz continuous. All of these assumptions hold for Huber’s loss function as well as others [17]. Figure 1 presents the family of loss functions, $\zeta(t)$, $\zeta^{\prime}(t)$, $\varphi(t)$, and $\zeta^{\prime\prime}(t)$ (second derivative of $\zeta(t)$).

Essentially Eq. (2) does not have a closed form solution, but using KIRWLS, the solution of robust kernel mean is,

[TABLE]

where $w_{i}^{(h)}=\frac{\varphi(\|\tilde{\mathbf{\Phi}}(X_{i})-f^{(h)}\|_{\mathcal{H}_{X}})}{\sum_{b=1}^{n}\varphi(\|\mathbf{\Phi}(X_{b})-f^{(h)}\|_{\mathcal{H}_{X}})}\,,\rm{and}\,\varphi(x)=\frac{\zeta^{\prime}(x)}{x}.$

Given the weights of the robust kernel ME, $\mathbf{w}=[w_{1},w_{2},\cdots,w_{n}]^{T}$, of a set of observations $X_{i},\cdots,X_{n}$, the points $\tilde{\mathbf{\Phi}}(X_{i}):=\mathbf{\Phi}(X_{i})-\sum_{a=1}^{n}w_{a}\mathbf{\Phi}(X_{a})$ are centered and the centered robust Gram matrix is $\tilde{K}_{ij}=\langle\tilde{\mathbf{\Phi}}(X_{i}),\tilde{\mathbf{\Phi}}(X_{j})\rangle=(\mathbf{C}\mathbf{K}_{X}\mathbf{C}^{T})_{ij}$, where $\mathbf{K}_{X}=(k_{X}(X_{i},X_{j}))_{i=1}^{n}$ is a Gram matrix, $\mathbf{1}_{n}=[1_{1},1_{2},\cdots,1_{n}]^{T}$ and $\mathbf{C}=\mathbf{I}-\mathbf{1}_{n}\mathbf{w}^{T}$. For a set of test points $X^{t}_{1},X_{2}^{t},\cdots,X_{T}^{t}$, we define two matrices of order $T\times n$ as $K_{ij}^{test}=\langle\mathbf{\Phi}(X^{t}_{i}),\mathbf{\Phi}(X_{j})\rangle$ and $\tilde{K}_{ij}^{test}=\langle\mathbf{\Phi}(X^{t}_{i})-\sum_{b=1}^{n}w_{b}\mathbf{\Phi}(X_{b}),\mathbf{\Phi}(X_{j})-\sum_{d=1}^{n}w_{d}\mathbf{\Phi}(X_{b})\rangle$. Like the centered Gram matrix, the centered robust Gram matrix of test points, $K_{ij}^{test}$, in terms of the robust Gram matrix and $\mathbf{1}_{t}=[1_{1},1_{2},\cdots,1_{t}]^{T}$ is defined as,

[TABLE]

2.4 Standard kernel (cross-) covariance operator

In this section we study the covariance of two random feature vectors $k_{X}(\cdot,X)$ and $k_{Y}(\cdot,Y)$. As for the standard random vectors, the notion of kernel covariance is useful as the basis in describing the statistical dependence among two or more variables.

Let $(\mathcal{X},\mathcal{B}_{X})$ and $(\mathcal{Y},\mathcal{B}_{Y})$ be two measurable spaces and $(X,Y)$ be a random variable on $\mathcal{X}\times\mathcal{Y}$ with distribution $F_{XY}$. The kernel CCO (centered) is a linear operator $\Sigma_{XY}:=\mathcal{H}_{Y}\to\mathcal{H}_{X}$ defined as

[TABLE]

where $\tilde{\mathbf{\Phi}}(\cdot)=\mathbf{\Phi}(\cdot)-E[\mathbf{\Phi}(\cdot)]$ and $\otimes$ is a tensor product operator $((a_{1}\otimes b_{1})x=\langle x,b_{1}\rangle a_{1}\,\rm{and}\,\langle a_{1}\otimes a_{2},b_{1}\otimes b_{2}\rangle_{\mathcal{H}_{12}}=\langle a_{1},b_{1}\rangle_{\mathcal{H}_{1}}\langle a_{2},b_{2}\rangle_{\mathcal{H}_{2}},\,\forall\,a_{1},b_{1}\in\mathcal{H}_{1},\,\rm{and}\,a_{2},b_{2}\in\mathcal{H}_{2}$, where $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$ are Hilbert spaces) [36].

Given two $k_{X}$ and $k_{Y}$ measurable positive definite kernels with respective RKHS $\mathcal{H}_{X}$ and $\mathcal{H}_{Y}$. By the reproducing property, the kernel CCO, with $\mathbb{E}_{X}[k_{X}(X,X)]<\infty$, and $\mathbb{E}_{Y}[k_{Y}(Y,Y)]<\infty$ is satisfied

[TABLE]

for all $f_{X}\in\mathcal{H}_{X}$ and $f_{Y}\in\mathcal{H}_{Y}$. This is a bounded operator. As shown in Eq. (1), we can define kernel CCO as an empirical risk optimization problem as follows,

[TABLE]

The empirical kernel CCO is then

[TABLE]

where $\tilde{k}_{X}$ and $\tilde{k}_{Y}$ are centered kernels. For the special case, when $Y$ is equal to $X$, it gives a kernel CO.

3 Robust kernel (cross-) covariance operator

Because a robust kernel ME (see Section 2.3) is used, to reduce the effect of outliers, we propose to use $M$-estimation to find a robust sample covariance of $\mathbf{\Phi}(X)$ and $\mathbf{\Phi}(Y)$. To do this, we estimate kernel CO and kernel CCO based on robust loss functions, namely, robust kernel CO and robust kernel CCO, respectively. Eq. (3) can be written as

[TABLE]

3.1 Representation of robust kernel (cross-) covariance operator

In this section, we represent $\widehat{\Sigma}_{RXY}$ as a weighted combination of the product of two kernels $k_{X}(\cdot,X_{i})k_{Y}(\cdot,Y_{i})$. We will also address necessary and sufficient conditions for the robust kernel CCO. Eq (5) can be reformulated as $\hat{\Sigma}_{RXY}=\operatornamewithlimits{argmin}_{\Sigma\in{\mathcal{H}_{X}\otimes\mathcal{H}_{Y}}}J(\Sigma),$ where

[TABLE]

In order to optimize $J$ in a product RKHS, the necessary conditions are characterized through the Gâteaux differentials of $J$. Given a product vector space $\mathcal{X}\times\mathcal{Y}$ and a function $A:\mathcal{X}\times\mathcal{Y}\to[-\infty,\infty]$, the Gâteaux differential of $A$ at $Z=(X,Y)\in\mathcal{X}\times\mathcal{Y}$ with incremental $\mathcal{T}\in\mathcal{X}\times\mathcal{Y}$ is defined as

[TABLE]

The Gâteaux differential on a probability distribution is also defined in Section 4.

Based on the optimality principle [37], the Gâteaux differential is well defined for all $\mathcal{T}$ and a necessary condition for $A$ to have a minimum at $Z_{0}=(X_{0},Y_{0})\in\mathcal{X}\times\mathcal{Y}$ is that $\partial A(Z_{0};\mathcal{T})=0$. We can state the following lemma.

Lemma 3.1

Under the assumptions (i) and (ii) the Gâteaux differential of the objective function $J$ at $\Sigma\in\mathcal{H}_{X}\otimes\mathcal{H}_{Y}$ and incremental $\mathcal{T}\in\mathcal{H}_{X}\otimes\mathcal{H}_{Y}$ is

[TABLE]

where $S:\mathcal{H}_{X}\otimes\mathcal{H}_{Y}\to\mathcal{H}_{X}\otimes\mathcal{H}_{Y}$ is defined as

[TABLE]

A necessary condition for $\Sigma=\widehat{\Sigma_{RXY}}$, robust kernel CCO is $S(\Sigma)=\mathbf{0}.$

The key difference of Lemma 3.1 and Lemma $1$ of [17] is the RKHS. The latter lemma is based on a single RKHS $\mathcal{H}_{X}$ but the former one is on a product RKHS $\mathcal{H}_{X}\times\mathcal{H}_{Y}$. This is a generalization result.

Theorem 3.1

Under the same assumption of Lemma 3.1, the robust kernel CCO (centered) for any $(X,Y)\in\mathcal{X}\times\mathcal{Y}$ is then

[TABLE]

where $w_{i}\geq 0$, and $\sum_{i=1}^{n}w_{i}=1$. Furthermore,

[TABLE]

Representer Theorem 3.1 tells us that in the robust loss function, when $\varphi$ is decreasing the large value of $\|\tilde{\mathbf{\Phi}}(X_{i})\otimes\tilde{\mathbf{\Phi}}(Y_{i})-\widehat{\Sigma}_{RXY}\|_{\mathcal{H}_{X}\otimes\mathcal{H}_{Y}}$, $w_{i}$ will be small. Therefore, the robust kernel CCO is robust in the sense that it down-weights outlying points.

In order to state the sufficient condition for $\widehat{\Sigma}_{RXY}$ to be the minimizer of Eq. (5), we need an additional assumption on $J$.

Theorem 3.2

Under the assumptions (i), (ii), and $J$ is strictly convex, Eq. (7), Eq. (8) and $\sum_{i=1}^{n}w_{i}=1$ are sufficient conditions for the robust kernel CCO to be the minimizer of Eq. (5).

For a positive definite kernel, $J$ becomes strictly convex for the Huber loss function.

3.2 Algorithm for robust kernel (cross-) covariance operator

As explained in [17], Eq. (5) does not have a closed form solution, but using the kernel trick the standard IRWLS can be extended to a RKHS. The solution at hth iteration is then,

[TABLE]

where $w_{i}^{(h)}=\frac{\varphi(\|\tilde{\mathbf{\Phi}}(X_{i})\otimes\tilde{\mathbf{\Phi}}(Y_{i})-\Sigma^{(h)}\|_{\mathcal{H}_{X}\otimes\mathcal{H}_{Y}})}{\sum_{b=1}^{n}\varphi(\|\tilde{\mathbf{\Phi}}(X_{b})\otimes\tilde{\mathbf{\Phi}}(Y_{b})-\Sigma^{(h)}\|_{\mathcal{H}_{X}\otimes\mathcal{H}_{Y}})}\,,\rm{and}\,\varphi(x)=\frac{\zeta^{\prime}(x)}{x}.$

Theorem 3.3

Under the assumptions (i) - (iii) and $\varphi(t)$ is non-increasing. Let

[TABLE]

and $\{\Sigma^{(h)}\}_{h=1}^{\infty}$ be the sequence produced by the KIRWLS algorithm. Then $J(\Sigma^{(h)})$ decreases monotonically at every iteration and converges.

[TABLE]

as $h\to\infty$.

Theorem 3.3 sates that $\Sigma^{(h)}$ becomes close to the set of stationary points of $J$ by increasing the number of iterations. Under the assumptions of Theorem 3.3 and for a strictly convex set $J$, it is also granted that the $\{\Sigma^{(h)}\}_{h=1}^{\infty}$ converges to $\Sigma_{RXY}$ in the Hilbert-Schmidt norm and supremum norm.

The algorithm for estimating robust kernel CCO is given in Figure 2. The input of this algorithm is a robust kernel ME. The computational complexity of a robust kernel ME is $\mathcal{O}(n^{2})$ in each iteration, where $n$ is the number of data points. The algorithm that we have presented involves finding the robust kernel CCO with the dimension $n\times n$. A naive implementation of the algorithm in Figure 2 would show that both time and memory complexity are similar to $\mathcal{O}(n^{3})$ in each iteration. In practice, the required number of iterations is around $50$. A computational complexity with cubic growth in the number of data points would be a serious liability in application to large dataset. We are able to reduce the time complexity using the low-rank approximation of the Gram matrix [38]. We can also use the random features approach. Random Features provide a finite-dimensional alternative to the kernel trick by instead mapping the data to an equivalent randomized feature space [39].

4 Influence function of robust kernel and kernel (cross-) covariance operator

To define the robustness in statistics, different approaches have been proposed, for example, the minimax approach [40], the sensitivity curve [35], the IF [41, 34] and the finite sample breakdown point [42]. Due to its simplicity, the IF is the most useful approach in statistical supervised learning [13, 12]. In this section, we briefly introduce the definition of IF and the IF of kernel ME, kernel CO, and kernel CCO. We then propose the IF of robust kernel CO and the robust kernel CCO.

Let $X_{1},X_{2},\cdots,X_{n}\in\mathcal{X}$ is a IID sample from a population with distribution function $F$, its empirical distribution function is $F_{n}$, and $T_{n}=T_{n}(X_{1},X_{2},\cdots,X_{n})$ is a statistic. Also let $\mathcal{A}(\mathcal{X})$ be a class of all possible distributions containing $F_{n}$ for all $n\geq 1$ and $F$. We assume that there exists a functional $T:\mathcal{D}\to\mathbb{R}$, where $\mathcal{D}$ is the set of all probability distributions in $\mathcal{A}(\mathcal{X})$ for which $T$ is defined, such that

[TABLE]

where $T$ does not depend on $n$. $T$ is then called a statistical functional. If the domain of $T$ is a convex set containing all distributions, $\mathcal{D}$ and the data do not follow the model $F$ in $\mathcal{D}$ exactly but slightly going toward a distribution $G$. The Gâteaux derivative, $T_{F}^{\prime}$ of $T$ at $F$ is defined as

[TABLE]

The Gâteaux differentiability at $F$ ensures the directional derivative of $T$ exists in all directions that stay in $\mathcal{D}$.

Suppose $X^{\prime}\in\mathcal{X}$ and $G=\Delta_{X^{\prime}}\in\mathcal{D}$ is the probability measure which gives mass $1$ at the point $\{X^{\prime}\}$. Then, $F^{\epsilon}=(1-\epsilon)F+\epsilon\Delta_{X^{\prime}}$ is a $\epsilon-$ contaminated distribution. The influence function (special case of Gâteaux Derivative) of $T$ at $F$ is defined by

[TABLE]

provided that the limit exists. It can be intuitively interpreted as a suitably normalized asymptotic influence of outliers on the value of an estimate or test statistic. The IF exists with an even weaker condition than Gâteaux differentiability. The IF reflects the bias caused by adding a few outliers at the point $X^{\prime}$, standardized by the amount of contamination. Therefore a bounded IF accelerates the robustness of an estimator [34].

4.1 Influence function based robustness measures

The three metrics of the IF function that can be used for robustness measures of the functional $T$ are the gross error sensitivity, local shift sensitivity and rejection point. The gross error sensitivity of $T$ at $F$ is defined as

[TABLE]

The gross error sensitivity measures the worst effect that a small amount of contamination of fixed size can have on the estimator. The local shift sensitivity of $T$ at $F$ for all $X_{1},X_{2}\in\mathcal{X}$ is defined by

[TABLE]

$\lambda^{*}$ measures the worst effect of rounding error (small function in the observation). The rejection point of $T$ at $F$ is defined by

[TABLE]

The $\rho_{\mathcal{F}}^{*}$ is infinite if there exits no such $m$. We can reject those observations, which are farther away than $\rho_{\mathcal{F}}^{*}$. For a robust estimator, $\rho_{\mathcal{F}}^{*}$ will be finite.

4.2 Influence function of kernel (cross-) covariance operator

In kernel methods, every estimate is a function. For a scalar-valued estimate, we define the IF at a fixed point. But if the estimate is a function, we are able to express the change of the function value at every point. Suppose $T(\cdot,F)$ and $T(\cdot,F^{\epsilon})$ are two function estimates on the distribution $F_{X}$ and the contaminated distribution $F^{\epsilon}$ at $X^{\prime}$, respectively. The influence function for $T(\cdot,F)$ is defined by

[TABLE]

We can estimate the IF using the empirical distribution which is called empirical IF (EIF). Suppose a sample of size $n$ is drawn from the empirical distributions $F_{n}$. Also let $F^{\epsilon}_{n}$ be a contamination model with the empirical data. The empirical IF for $T(\cdot,F_{n})$ is defined as

[TABLE]

As a first example, let the kernel ME, $T(\cdot,F_{X})=\int k_{X}(\cdot,X)dF_{X}=E_{X}[k_{X}(\cdot,X)]$, where $X\sim F_{X}.$ The value of the parameter at the contamination model, $F^{\epsilon}=(1-\epsilon)F_{X}+\epsilon\Delta_{X^{\prime}}$ is

[TABLE]

Thus the IF of kernel ME at point $X^{\prime}$ is given by

[TABLE]

We can estimate the IF of the kernel ME with the empirical distribution, $F_{n}$, at the data points $X_{1},X_{2},\cdots,X_{n}\sim F_{n}$, at $X^{\prime}$ for every point $X$ as

[TABLE]

which is called the EIF of kernel ME.

As a second example, let the mean of the product of two random variables, $f(X)$ and $f(Y)$ with $(X,Y)\in\mathcal{X}\otimes\mathcal{Y}$, $T(X,Y,F_{XY})=\mathbb{E}_{XY}[f_{X}(X)f_{Y}(Y)]$, for all $f_{X}\in\mathcal{H}_{X}$, and $f_{Y}\in\mathcal{H}_{Y}$. The value of parameter at the contamination model at $(X^{\prime},Y^{\prime})\in\mathcal{X}\otimes\mathcal{Y}$ , $F_{XY}^{\epsilon}=(1-\epsilon)F_{XY}+\epsilon\Delta_{X^{\prime}Y^{\prime}}$ is given by

[TABLE]

Thus the IF of $T(X,Y,F_{XY})$ is given by

[TABLE]

We can find the IF for a combined statistic given the IF for the statistic itself. The IF of complicated statistics can be calculated with the chain rule, say $T(F)=\ell(T_{1}(F),\cdots,T_{s}(F))$, that is,

[TABLE]

For example, the IF of covariance of two random variables, $f_{X}(X)$ and $f_{Y}(Y)$ can be calculated using the above chain rule as

[TABLE]

for $f_{X}\in\mathcal{H}_{X}$, $f_{Y}\in\mathcal{H}_{Y}$, and $Z=(X,Y)\in\mathcal{X}\times\mathcal{Y}$.

Using Eq. (11) and the reproducing property, the IF of $T(z,F_{XY})$ with distribution, $F_{XY}$ at $Z^{\prime}=(X^{\prime},Y^{\prime})$ is given by

[TABLE]

Letting $f_{X}=k_{X}(\cdot,X)$ and $f_{Y}=k_{Y}(\cdot,Y),\,\forall\,X\in\mathcal{X}\,\rm{and}\,Y\in\mathcal{Y}$ be two random variables taking values in $\mathcal{H}_{X}$ and $\mathcal{H}_{Y}$, the IF of kernel CCO at $Z=(X^{\prime},Y^{\prime})$ is formulated as

[TABLE]

where $k_{X}(\cdot,X)$, $k_{X}(\cdot,X^{\prime})$, $k_{Y}(\cdot,Y)$, and $k_{Y}(\cdot,Y^{\prime})$ are random vectors in $\mathcal{H}_{X}$ and $\mathcal{H}_{Y}$, respectively.

Given data points $(X_{1}Y_{1}),(X_{2},Y_{2}),\cdots,(X_{n},Y_{n})\in\mathcal{X}\times\mathcal{Y}$ from the joint empirical distribution, $F_{nXY}$, for every point $(X_{i},Y_{i})$, we can estimate the IF of the kernel CCO, called EIF of kernel CCO as follows,

[TABLE]

In case of the outliers, the bounded kernels take the values in a range. Thus, the above IFs have the three properties: gross error sensitivity, local shift sensitivity and rejection point only for the bounded kernels. These properties are not true for the unbounded kernels, for example, linear and polynomial kernels. The unbounded kernels take the arbitrary values and IFs reflects the bias. We can make a similar conclusion for the kernel CO.

4.3 Influence function of robust kernel (cross-) covariance operator

To derive the IF of the robust kernel CCO like the robust kernel DE as shown in [17], we generalize the definition of robust kernel CCO to a joint general distribution $\mu_{XY}$,

[TABLE]

Let $\widehat{\Sigma}_{RXY}(X,Y;F_{XY})=\Sigma_{F_{XY}}(X,Y)$, the IF for the robust kernel CCO at $(X^{\prime},Y^{\prime})$ is

[TABLE]

Similarly for the definition of robust kernel CCO, we generalize the necessary condition $S(\widehat{\Sigma}_{RXY})$ to $S_{F{XY}}(\widehat{\Sigma}_{RXY})$. Besides the assumptions $(i)-(iv)$ in Section 2.3, assume that $\Sigma_{F_{XY}^{\epsilon}}\to\Sigma_{F_{XY}}$ as $\epsilon\to 0$. We need to find the Gâteaux differentiability of $S_{F{XY}}$ as in proof of Lemma 3.1 (in the appendix). If $\dot{\Sigma}_{F_{XY}}\doteq\lim_{\epsilon\to 0}\frac{\Sigma_{F_{XY}^{\epsilon}}-\Sigma_{F_{XY}}}{\epsilon}$ exists, the IF of robust kernel CCO is defined as

[TABLE]

where $\dot{\Sigma}_{F_{XY}}\in\mathcal{H}_{X}\otimes\mathcal{H}_{Y}$ satisfies

[TABLE]

where $q(t)=t\psi^{\prime}(t)-\psi(t)$. Unfortunately, Eq. (16) has no closed form solution. By considering the empirical joint distribution, $F_{n}=F_{nXY}$ instead of the joint distribution, $F_{XY}$, we can find $\dot{\Sigma}_{F_{n}}$ explicitly. To do this, besides the assumptions $(i)-(iv)$ we assume that $\Sigma_{F_{n}^{\epsilon}}\to\Sigma_{F_{n}}$ as $\epsilon\to 0$ (satisfied when $J$ is strictly convex) and the extended kernel matrices $\mathbf{K}_{X}^{\prime}$ and $\mathbf{K}_{Y}^{\prime}$ with $(X_{i},Y_{i})_{i=1}^{n}\cup(X^{\prime},Y^{\prime})$ are positive definite. Then, the IF of robust kernel CCO with $(X,Y)$ at $(X^{\prime},Y^{\prime})$ is defined as

[TABLE]

where $\alpha^{\prime}=n\frac{\varphi(\|\tilde{\mathbf{\Phi}}(X)\otimes\tilde{\mathbf{\Phi}}(Y)-\Sigma_{F_{n}}\|_{\mathcal{H}_{X}\otimes\mathcal{H}_{Y}})}{\gamma}$, $\gamma=\sum_{i=1}^{n}\varphi(\|\tilde{\mathbf{\Phi}}(X_{i})\otimes\tilde{\mathbf{\Phi}}(Y_{i})-\Sigma_{F_{n}}\|_{\mathcal{H}_{X}\otimes\mathcal{H}_{Y}})$ and $\alpha=[\alpha_{1},\alpha_{2},\cdots,\alpha_{n}]^{T}$ are the solution of the following system of linear equations:

[TABLE]

where $\mathbf{1}=[1,\cdots,1]^{T}$, $\mathbf{k}_{XY}=[k_{X}(X^{\prime},X_{1})k_{Y}(Y^{\prime},Y_{1}),\cdots,k_{X}(X^{\prime},X_{n})k_{Y}(Y^{\prime},Y_{n})]$, $\mathbf{I}_{n}$ is a $n$ ordered identity matrix, $\mathbf{Q}$ is a diagonal matrix with $Q_{ii}=\frac{q(\|\tilde{\mathbf{\Phi}}(X_{i})\otimes\mathbf{\Phi}(Y_{i})-\Sigma_{F_{n}}\|_{\mathcal{H}_{X}\otimes\mathcal{H}_{Y}})}{\|\tilde{\mathbf{\Phi}}(X_{i})\otimes\mathbf{\Phi}(Y_{i})-\Sigma_{F_{n}}\|^{3}_{\mathcal{H}_{X}\otimes\mathcal{H}_{Y}}}$, and $\mathbf{w}=[w_{1},\cdots,w_{n}]^{T}$ gives the weights as in robust kernel CCO. $\alpha^{\prime}$ captures the amount of contaminated data in the robust kernel CCO, which is given as

[TABLE]

For a standard kernel CCO, we have $\varphi\equiv 1$ and $\alpha^{\prime}=1$, which is in agreement with the IF of standard kernel CCO. The robust loss function $\zeta$, $\varphi(\|\tilde{\mathbf{\Phi}}(X^{\prime})\otimes\tilde{\mathbf{\Phi}}(Y^{\prime})-\Sigma_{F_{n}}\|_{\mathcal{H}_{X}\otimes\mathcal{H}_{Y}})$ can be regarded as a measure of “inlyingness”, with more inlying points having larger values than $\alpha^{\prime}<1$. Thus, the robust kernel CCO is less sensitive to outlying points than the standard kernel CCO.

5 Standard and robust kernel canonical correlation analysis

In this section, we review standard kernel CCA and propose the IF and empirical IF (EIF) of kernel CCA. After that we propose a robust kernel CCA method based on robust kernel CO and robust kernel CCO.

5.1 Standard kernel canonical correlation analysis

Standard kernel CCA has been proposed as a nonlinear extension of linear CCA [8, 43]. Researchers have extended the standard kernel CCA with an efficient computational algorithm, i.e., incomplete Cholesky factorization [9]. Over the last decade, standard kernel CCA has been used for various tasks [44, 45, 46, 23]. Theoretical results on the convergence of kernel CCA have also been obtained [20, 21].

The aim of the standard kernel CCA is to seek the sets of functions in the RKHS for which the correlation (Corr) of random variables is maximized. For the simplest case, given two sets of random variables $X$ and $Y$ with two functions in the RKHS, $f_{X}(\cdot)\in\mathcal{H}_{X}$ and $f_{Y}(\cdot)\in\mathcal{H}_{Y}$, the optimization problem of the random variables $f_{X}(X)$ and $f_{Y}(Y)$ is

[TABLE]

The optimizing functions $f_{X}(\cdot)$ and $f_{Y}(\cdotp)$ are determined up to scale.

Using a finite sample, we are able to estimate the desired functions. Given an i.i.d sample, $(X_{i},Y_{i})_{i=1}^{n}$ from a joint distribution $F_{XY}$, by taking the inner products with elements or “parameters” in the RKHS, we have features $f_{X}(\cdot)=\langle f_{X},\mathbf{\Phi}_{X}(X)\rangle_{\mathcal{H}_{X}}=\sum_{i=1}^{n}a_{X}^{i}k_{X}(\cdot,X_{i})$ and $f_{Y}(\cdot)=\langle f_{Y},\mathbf{\Phi}_{Y}(Y)\rangle_{\mathcal{H}_{Y}}=\sum_{i=1}^{n}a_{Y}^{i}k_{Y}(\cdot,Y_{i})$, where $k_{X}(\cdot,X)$ and $k_{Y}(\cdot,Y)$ are the associated kernel functions for $\mathcal{H}_{X}$ and $\mathcal{H}_{Y}$, respectively. The kernel Gram matrices are defined as $\mathbf{K}_{X}:=(k_{X}(X_{i},X_{j}))_{i,j=1}^{n}$ and $\mathbf{K}_{Y}:=(k_{Y}(Y_{i},Y_{j}))_{i,j=1}^{n}$. We need the centered kernel Gram matrices $\mathbf{G}_{X}=\mathbf{C}\mathbf{K}_{X}\mathbf{C}$ and $\mathbf{G}_{Y}=\mathbf{C}\mathbf{K}_{Y}\mathbf{C}$, where $\mathbf{C}=\mathbf{I}_{n}-\frac{1}{n}\mathbf{D}_{n}$ with $\mathbf{D}_{n}=\mathbf{1}_{n}\mathbf{1}^{T}_{n}$ and $\mathbf{1}_{n}$ is the vector with $n$ ones. The empirical estimate of Eq. (18) is then given by

[TABLE]

where

[TABLE]

and $\mathbf{W}$ is a diagonal matrix with elements $\frac{1}{n}$, and $\mathbf{a}_{X}$ and $\mathbf{a}_{Y}$ are the eigen-direction of $X$ and $Y$, respectively. The regularized coefficient $\kappa>0$. Solving the maximization problem in Eq. (19) is analogous to solving the following generalized eigenvalue problem:

[TABLE]

[TABLE]

It is easy to show that the eigenvalues of Eq. (20) and Eq. (21) are equal and that the eigenvectors for any equation can be obtained from the other. The square roots of the eigenvalues of Eq.(20) or Eq. (21) are the estimated kernel CC, $\hat{\rho}$. The $\hat{\rho}_{j}$ is the jth largest kernel CC and the jth kernel CVs are $\mathbf{a}_{X}^{T}\mathbf{G}_{X}$, and $\mathbf{a}_{Y}^{T}\mathbf{G}_{Y}$.

Standard kernel CCA can be formulated using kernel CCO, which makes the robustness analysis easier. As in [20], using the cross-covariance operator of (X,Y), $\Sigma_{XY}:\mathcal{H}_{Y}\to\mathcal{H}_{X}$ we can reformulate the optimization problem in Eq. (18) as follows:

[TABLE]

As with linear CCA [47], we can derive the solution of Eq. (22) using the following generalized eigenvalue problem.

[TABLE]

The eigenfunctions of Eq. (5.1) correspond to the largest eigenvalue, which is the solution to the kernel CCA problem. After some simple calculations, we reset the solution as

[TABLE]

It is known that the inverse of an operator may not exist. Even if it exists, it may not be continuous in general [20]. We can derive kernel CC using the correlation operator $\Sigma_{YY}^{-\frac{1}{2}}\Sigma_{YX}\Sigma_{XX}^{-\frac{1}{2}}$, even when $\Sigma_{XX}^{-\frac{1}{2}}$ and $\Sigma_{YY}^{-\frac{1}{2}}$ are not proper operators. The potential danger is that it might overfit, which is why introducing $\kappa$ as a regularization coefficient would be helpful. Using the regularized coefficient $\kappa>0$, the empirical estimators of Eq. (22) and Eq. (23) are

[TABLE]

and

[TABLE]

respectively.

Now we calculate a finite rank operator $\mathbb{B}_{YX}=(\hat{\Sigma}_{{YY}}+\kappa\mathbf{I})^{-\frac{1}{2}}\hat{\Sigma}_{{YX}}(\hat{\Sigma}_{{XX}}+\kappa\mathbf{I})^{-\frac{1}{2}}$. For $\kappa>0$, the square roots of the j-th eigenvalue of $\mathbb{B}_{YX}$ are the j-th kernel CC, $\rho_{j}$. The unit eigenfunctions of $\mathbb{B}_{YX}$ corresponding to the jth eigenvalues are $\hat{\nu}_{jX}\in\mathcal{H}_{X}$ and $\hat{\nu}_{j_{Y}}\in\mathcal{H}_{Y}$. The jth ($j=1,2,\cdots,n$) kernel CVs are

[TABLE]

where $\hat{f}_{jX}=(\hat{\Sigma}_{{XX}}+\kappa\mathbf{I})^{-\frac{1}{2}}\hat{\nu}_{jX}$ and $\hat{f}_{jY}=(\hat{\Sigma}_{{YY}}+\kappa\mathbf{I})^{-\frac{1}{2}}\hat{\nu}_{jY}.$

The generalized eigenvalue problem in Eq. (23) can be formulated as a simple eigenvalue problem. Using the j-th eigenfunction in the first equation of Eq. (23) we have

[TABLE]

where $e_{jX}=\Sigma_{XX}^{\frac{1}{2}}f_{jX}$.

5.2 Influence function of the standard kernel canonical correlation analysis

By using the IF of kernel PCA, linear PCA and linear CCA, we can derive the IF of kernel CCA (kernel CC and kernel CVs). For simplicity, let us define $\tilde{f}_{X}(X)=\langle f_{X},\tilde{k}_{X}(\cdot,X)$, $\mathbb{L}_{jX}=\Sigma_{XX}^{-\frac{1}{2}}(\Sigma_{XX}^{-\frac{1}{2}}\Sigma_{XY}\Sigma_{YY}^{-1}\Sigma_{YX}\Sigma_{XX}^{-\frac{1}{2}}-\rho^{2}_{j}\mathbf{I})^{-1}\Sigma_{XX}^{-\frac{1}{2}}$, and $\mathbb{L}_{jY}=\Sigma_{YY}^{-\frac{1}{2}}(\Sigma_{YY}^{-\frac{1}{2}}\Sigma_{YX}\Sigma_{XX}^{-1}\Sigma_{XY}\Sigma_{YY}^{-\frac{1}{2}}-\rho^{2}_{j}\mathbf{I})^{-1}\Sigma_{YY}^{-\frac{1}{2}}$.

Theorem 5.4

Given two sets of random variables $(X,Y)$ having the distribution $F_{XY}$ and the j-th kernel CC ( $\rho_{j}$) and kernel CVs ($f_{jX}(X)$ and $f_{jX}(Y)$), the influence functions of kernel CC and kernel CVs at $Z^{\prime}=(X^{\prime},Y^{\prime})$ are

[TABLE]

The above theorem has been proved on the basis of previously established ones, such as the IF of linear PCA [48, 49], the IF of linear CCA [25], and the IF of kernel PCA, respectively. To do this, we convert the generalized eigenvalue problem of kernel CCA into a simple eigenvalue problem. First, we need to find the IF of $\Sigma_{XX}^{-\frac{1}{2}}\Sigma_{XY}\Sigma_{YY}^{-1}\Sigma_{YX}\Sigma_{XX}^{-\frac{1}{2}}$, henceforth the IF of $\Sigma_{YY}^{-1},\Sigma_{XX}^{\frac{1}{2}}$ and $\Sigma_{XY}$.

Using the above result, we can establish some properties of kernel CCA: robustness, asymptotic consistency and its standard error. In addition, we are able to identify the outliers based on the influence of the data. All notations and proof are explained in the appendix.

The IF of inverse covariance operator exists only for the finite dimensional RKHS. For infinite dimensional RKHS, we can find the IF of $\Sigma_{XX}^{-\frac{1}{2}}$ by introducing a regularization term as follows

[TABLE]

which gives the empirical estimator.

Let $(X_{i},Y_{i})_{i=1}^{n}$ be a sample from the empirical joint distribution $F_{nXY}$. The EIF (IF based on empirical distribution) of kernel CC and kernel CVs at $(X^{\prime},Y^{\prime})$ for all points $(X_{i},Y_{i})$ are $\rm{EIF}(X_{i},Y_{i},X^{\prime},Y^{\prime},\rho_{j}^{2})=\widehat{\rm{IF}}(X^{\prime},Y^{\prime},\hat{\rho}_{j}^{2})$, $\rm{EIF}(X_{i},Y_{i},X^{\prime},Y^{\prime},f_{jX})=\widehat{\rm{IF}}(\cdot,X^{\prime},Y^{\prime},f_{jX})$, and $\rm{EIF}(X_{i},Y_{i},X^{\prime},Y^{\prime},f_{jY})=\widehat{\rm{IF}}(\cdot,X^{\prime},Y^{\prime},\widehat{f}_{jY})$, respectively.

For the bounded kernels, the IFs defined in Theorem 5.4 have three properties: gross error sensitivity, local shift sensitivity, and rejection point. But for unbounded kernels, say a linear, polynomial and so on, the IFs are not bounded. Consequently, the results of standard kernel CCA using the bounded kernels are less sensitive than the standard kernel CCA using the unbounded kernels. In practice, standard kernel CCA is sensitive to the contaminated data even with the bounded kernels [24].

5.3 Robust kernel canonical correlation analysis

In this section, we propose a robust kernel CCA method based on the robust kernel CO and the robust kernel CCO. While many robust linear CCA methods have been proposed to show that linear CCA methods cannot fit the bulk of the data and have points deviating from the original pattern for further investment [50, 24], there are no well-founded robust methods of kernel CCA. The standard kernel CCA considers the same weights for each data point, $\frac{1}{n}$, to estimate kernel CO and kernel CCO, which is the solution of an empirical risk optimization problem when using the quadratic loss function. It is known that the least square loss function is not a robust loss function. Instead, we can solve an empirical risk optimization problem using the robust least square loss function where the weights are determined based on KIRWLS. We need robust centered kernel Gram matrices of $X$ and $Y$ data. The centered robust kernel Gram matrix of $X$ is, $\mathbf{G}_{RX}=\mathbf{C_{R}}\mathbf{K}_{X}\mathbf{C_{R}}$ where $\mathbf{C}_{R}=\mathbf{I}_{n}-\mathbf{1}_{n}\mathbf{w}^{T}_{n}$, $\mathbf{1}_{n}$ is the vector with $n$ ones and $\mathbf{w}$ is a weight vector of robust kernel ME, $\mathcal{M}_{X}$. Similarly, we can calculate $\mathbf{G}_{RY}$ for $Y$. After getting robust kernel CO and kernel CCO, they are used in standard kernel CCA, which we call robust kernel CCA. The empirical estimate of Eq. (18) is then given by

[TABLE]

with for all $g_{X}\in{\mathcal{H}}_{X}$, $g_{Y}\in{\mathcal{H}}_{Y}$ and

[TABLE]

where $\mathbf{W}_{XY}$, $\mathbf{W}_{XX}$, and $\mathbf{W}_{YY}$ are diagonal matrices with elements corresponding to the weights of robust kernel CCO, and kernel COs, respectively. Also $\mathbf{b}_{X}$ and $\mathbf{b}_{Y}$ are the eigen-direction of $X$ and $Y$, respectively. As in Eq. (20), we can solve the maximization problem of Eq. (5.3) as an eigenvalue problem. Let $\Sigma_{RXY}$, $\Sigma_{RXX}$, and $\Sigma_{RYY}$ be the robust kernel CCO, robust kernel CO of $X$, and robust kernel CO of $Y$, respectively. Like standard kernel CCA, the robust empirical estimators of Eq. (22) and Eq. (23) are

[TABLE]

and

[TABLE]

respectively. Figure 3 presents a detailed algorithm of the proposed methods (all steps are similar to standard kernel CCA except the first one). This method is designed for contaminated data, and the principles we describe also apply to the kernel methods, which must deal with the issue of kernel CO and kernel CCO.

It is well-known that robust methods have higher time complexity than the standard methods. At each update of the robust kernel CO or robust kernel CCO, we need to store the $n\times n$ matrix. The memory complexity of robust kernel CCA is then $\mathcal{O}(n^{3})$. A naive implementation of the algorithm in Figure 3 would therefore require $\mathcal{O}(n^{3}h))$ operations (the time complexity), where $h$ is the number of iterations. The spectrum of Gram matrices tends to show rapid decay, and low-rank approximations of Gram matrices can often provide sufficient fidelity for the needs of kernel-based algorithms [38, 51, 9]. By assuming that the outliers have a similar effect on marginal distribution and the joint distribution, we can also reduce the memory complexity and time complexity. Under this assumption, we estimate the weight of kernel CCO and consider this weight for kernel CO of $X$ and $Y$ data.

6 Experiments

We conducted experiments on both the synthetic and real data sets. We generated two types of simulated data: ideal data and those with $5\%$ of contamination. The description of real data sets are in Sections 6.3. The five synthetic data sets are as follows:

Three circles of structural data (TCSD): Data are generated along three circles of different radii with small noise:

[TABLE]

where $r_{i}=1$, $0.5$ and $0.25$, for $i=1,\ldots,n_{1}$, $i=n_{1}+1,\ldots,n_{2}$, and $i=n_{2}+1,\ldots,n_{3}$, respectively, $\mathbf{Z}_{i}\sim U[-\pi,\pi]$ and $\epsilon_{i}\sim\mathcal{N}(0,0.01\,I_{2})$ independently for the ideal data and $\mathbf{Z}_{i}\sim U[-10,10]$ for the contaminated data.

Sine function of structural data (SFSD): 1500 data points are generated along the sine function with small noise:

[TABLE]

where $\mathbf{Z}_{i}\sim U[-2\pi,0]$ and $\epsilon_{i}\sim\mathcal{N}(0,0.01\,I_{10})$ independently for the ideal data and $\epsilon_{i}\sim\mathcal{N}(0,10\,I_{10})$ for the contaminated data.

**Multivariate Gaussian structural data (MGSD): ** Given multivariate normal data, $\mathbf{Z}_{i}\in\mathbb{R}^{12}\sim\mathbf{N}(\mathbf{0},\Sigma)$ ($i=1,2,\ldots,n$), where $\Sigma$ is the same as in [23]. We divide $\mathbf{Z}_{i}$ into two sets of variables ($\mathbf{Z}_{i1}$,$\mathbf{Z}_{i2}$), and use the first six variables of $\mathbf{Z}_{i}$ as $X$ and perform the $\log$ transformation of the absolute value of the remaining variables ($\log_{e}(|\mathbf{Z}_{i2}|))$) as $Y$. For the contaminated data $\mathbf{Z}_{i}\in\mathbb{R}^{12}\sim\mathbf{N}(\mathbf{1},\Sigma)$ ($i=1,2,\ldots,n$).

Sine and cosine function structural data (SCFSD): We use uniform marginal distribution, and transform the data by two periodic $\sin$ and $\cos$ functions to make two sets $X$ and $Y$, respectively, with additive Gaussian noise: $Z_{i}\sim U[-\pi,\pi],\,\eta_{i}\sim N(0,10^{-2}),~{}\,i=1,2,\ldots,n,X_{ij}=\sin(jZ_{i})+\eta_{i},\,Y_{ij}=\cos(jZ_{i})+\eta_{i},j=1,2,\ldots,100.$ For the contaminated model $\eta_{i}\sim N(1,10^{-2})$.

SNP and fMRI structural data (SMSD): Two data sets of SNP data X with $1000$ SNPs and fMRI data Y with 1000 voxels were simulated. To correlate the SNPs with the voxels, a latent model is used as in [52]). For simulation of contamination, we consider the signal level, $0.5$ and noise level, $1$ to $10$ and $20$, respectively.

In our experiments, first, we compare standard and robust kernel covariance operators. After that, the robust kernel CCA is compared with the standard kernel CCA. For the Gaussian kernel we use the median of the pairwise distance as a bandwidth and for the Laplacian kernel we set the bandwidth equal to 1. As shown in [45], we can optimize the regularization parameter but our goal is the robustness issue of different methods. Thus the regularization parameter in standard kernel CCA and robust kernel CCA is fixed as $\kappa=10^{-5}$. In robust methods, we consider Huber’s loss function with the constant, $c$, equal to the median of error.

6.1 Results of kernel CCO and robust CCO

We evaluate the performance of kernel CO and robust kernel CO in two different settings. First, we check the accuracy of both operators by considering the kernel CO (KCO) with large data (say a population kernel CO of size $N$) and kernel CO with small data (say a sample kernel CO of size $n$). Now, we can estimate the distance between sample kernel CO and population kernel CO as $\|\widehat{\Sigma}_{XX}-\Sigma_{XX}\|^{2}_{\mathcal{H}_{X}\otimes\mathcal{H}_{X}}=\|\widehat{\Sigma}_{XX}\|^{2}_{\mathcal{H}_{X}\otimes\mathcal{H}_{X}}-2\langle\widehat{\Sigma}_{XX},\Sigma_{XX}\rangle_{\mathcal{H}_{X}\otimes\mathcal{H}_{X}}+\|\Sigma_{XX}\|^{2}_{\mathcal{H}_{X}\otimes\mathcal{H}_{X}}$ and similarly for the robust kernel CO (RKCO). Thus, the performance measures of the kernel CO and robust kernel CO estimators are defined as

[TABLE]

and

[TABLE]

respectively.

In theory, the above two equations become zero for large population size, $N$, with the sample size, $n\to N$. To do this, we consider the synthetic data, TCSD with $N\in\{1500,3000,6000,9000\}$ and $n\in\{15,30,45,60,90,120,150,180,210,240,270,300\}$ ($n=n_{1}+n_{2}+n_{3}$). For each $n$, we take $5\%$ CD. We repeated the process for $100$ samples to confirm our findings. The results (mean with standard error) were plotted in Figure 4. These figures show that both estimators give similar performances for small sample sizes, but for large sample sizes the robust estimator (i.e., robust kernel CO) shows much better results than the kernel CO estimate at all population sizes.

In addition, we compare kernel CO and robust kernel CO estimators using five kernels: linear (Poly-1), a polynomial with degree $2$ (Poly-2) and polynomial with degree $3$ (Poly-3), Gaussian and Laplacian on two synthetic data sets: TCSD and SFSD. To measure the performance, we use two matrix norms: Frobenius norm (F) and maximum modulus of all the elements (M) [53]. We calculate the ratio between ideal model and contaminated model for the kernel CO. The ratio becomes zero if the estimator is not sensitive to contaminated data. For both estimators, kernel CO and robust kernel CO, we use the following performance measures,

[TABLE]

and

[TABLE]

respectively. We repeated the experiment for $100$ samples with sample size, $n=1500$. The results (mean $\pm$ standard deviation) for kernel CO (standard) and robust kernel CO (Robust) are tabulated in Table 1. From this table, it is clear that the robust estimator performs better than the standard estimator in all cases. Moreover, both estimators using Gaussian and Laplacian kernels are less sensitive than all polynomial kernels.

6.2 Visualizing influential subject using standard kernel CCA and robust kernel CCA

We evaluated the performance of the proposed method for three different settings. First, we compared robust kernel CCA with the standard kernel CCA using Gaussian kernel (same bandwidth and regularization). To measure the influence, we calculated the ratio of IF for kernel CC between ideal data and contaminated data. We also calculated a similar measure for the kernel CV. Based on these ratios, we defined two performance measures on kernel CC and kernel CVs

[TABLE]

respectively. For any method, that does not depend on the contaminated data, the above measures, $\eta_{\rho}$ and $\eta_{f}$, should be approximately zero. In other words, the best methods should give small values. To compare, we considered three simulated data sets: MGSD, SCFSD, SMSD with three sample sizes, $n\in\{100,500,1000\}$. For each sample size, we repeated the experiment for $100$ samples. Table 2 presents the results (mean $\pm$ standard deviation) of the standard kernel CCA and robust kernel CCA. From this table, we observed that the robust kernel CCA outperforms the standard kernel CCA in all cases.

Second, we considered a simple graphical display based on the EIF of kernel CCA, the index plots (the subject on the $x$-axis and the influence, $\eta_{\rho}$, on the $y$ axis), to assess the related influence in data fusion regarding EIF based on kernel CCA, $\eta_{\rho}$. To do this, we considered a simulated data set, SMSD. The index plots of the standard kernel CCA and robust kernel CCA using the SMSD are presented in Figure 5. The $1$st and $2$nd rows are for the ideal and contaminated, and $1$st and $2$nd columns are for the standard kernel CCA (Standard kernel CCA) and robust kernel CCA (Robust kernel CCA), respectively. These plots show that both methods have almost similar results for the ideal data. But for contaminated data, the standard kernel CCA is affected by the contaminated data significantly. We can easily identify the influence of observation using this visualization. On the other hand, the robust kernel CCA has almost similar results for the ideal and contaminated data.

Table 3 presents the mean and standard deviation for the difference between the training and test of correlation for $10$-fold cross-validation of MGSD and SCFSD simulated data, using standard kernel CCA and robust kernel CCA. From the table, we can conclude that standard kernel CCA is sensitive to the contamination for both data sets. On the other hand, the robust kernel CCA is not only less sensitive to the contaminated data, but also performs better than the standard kernel CCA.

6.3 Application to imaging genetics data from MCIC and TCGA

To demonstrate the application of the proposed methods, we used three data sets: the Mind Clinical Imaging Consortium (MCIC) and two data sets from the Cancer Genome Atlas (TCGA) project. The MCIC has collected three types of data: SNPs (723,404 loci), fMRI (51,056 voxels) and DNA methylation (9273 methylation profiles) from 208 subjects including $92$ schizophrenic patients (age: $34\pm 11$, $22$ females) and $116$ (age: $32\pm 11$, $44$ females) healthy controls. Without missing information, the number of subjects is reduced to $183$ ($79$ schizophrenia (SZ) patients and $104$ healthy controls). The detailed information of the MCIC data set is given in [54]. In addition, we consider ovarian serous cystadenocarcinoma (OVSC) and lung squamous cell carcinoma (LUSC) data sets from TCGA data portal. The RNA-Seq gene expression data and methylation profiles are selected from the OVSC and the LUSC patients. After merging the RNA-Seq and methylation data, the number of OVSC patients and LUSC patients are $294$ and $130$, respectively (https://tcga-data.nci.nih.gov/tcga/).

To detect influential subjects, we use the EIF of the kernel CC for the standard and robust kernel CCA methods. For robust kernel CCA, we consider robust kernel CC and kernel CVs as in Theorem 5.4. However, both standard and robust kernel CCA have identified a similar subject, but robust kernel CCA is less sensitive than standard kernel CCA. After getting the influence of the subject, we extracted the outlier subjects of each data set based on the ‘getOutliers’ function of “extremevalues” R packages. The outlier subjects of SNP and fMRI; SNP and Methylation; and fMRI and Methylation are

[TABLE]

[TABLE]

and

[TABLE]

respectively. We observed that the SZ patient number $58$ was common in all cases. In the clinical assessment, this patient has high current psychosis disorder diagnosis rate ($259.9$). For TCGA data, outlier patients of OVSC and LUSC are

[TABLE]

and

[TABLE]

respectively.

Finally, we investigated the difference between training and testing for correlations using $10$ fold cross-validation. Table 4 shows the results of all subjects and all but outlier subjects of MCIC and TCGA data sets using standard kernel CCA and robust kernel CCA. When comparing the subjects with the outliers to the subjects without the outliers the standard kernel CCA method produces drastically different results. Whereas when the robust kernel CCA method is used to compare the two, the results are similar.

7 Concluding remarks and future research

The robust estimator employs a robust loss function instead of a quadratic loss function for the analysis of contaminated data. The robust estimators are weighted estimators where smaller weights are given more outlying data points. The weights can be estimated efficiently using a KIRLS approach. In terms of accuracy and sensitivity, it is clear that the robust estimators (e.g., robust kernel CO and robust kernel CCO) perform better than standard estimators (e.g., kernel CO and kernel CCO). We propose the IF of kernel CCA (kernel CC and kernel CVs) and robust kernel CCA based on the robust kernel CO and robust kernel CCO. The proposed IF measures the sensitivity of kernel CCA, which shows that the standard kernel CCA is sensitive to the contamination, but the proposed robust kernel CCA is not. The visualization method can identify influential (outlier) data in both synthesized and real imaging genetics data.

While M- estimator based methods are robust with a high breakdown point, finding the theoretical IF of robust kernel CCA is a future research direction. Although the focus of this paper is on kernel CCA, we are able to robustify other kernel methods, which must deal with the issue of kernel CO and kernel CCO. In future work, it would be interesting to develop robust multiple kernel PCA and robust multiple weighted kernel CCA.

Appendix A: Proofs

We recall some definitions on Hilbert spaces. Hilbert-Schmidt operator and Hilbert-Schmidt norm will be used in the proofs of Lemma 3.1, Theorem 3.1 and 3.2.

Let $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$ be separable Hilbert spaces. A linear operator $\mathcal{T}:\mathcal{H}_{1}\to\mathcal{H}_{2}$ is called the Hilbert-Schmidt operator if $\sum_{i}\|\mathcal{T}\phi_{i}\|^{2}<\infty$ for an orthonormal basis $\{\phi_{i}\}_{i\in I}$ of $\mathcal{H}_{1}$ with index set $I$. The sum $\sum_{i\in I}\|\mathcal{T}\phi_{i}\|^{2}$ does not depend on the orthonormal basis $\{\phi_{i}\}_{i}$. The square root of this sum is the called Hilbert-Schmidt norm, $\|\mathcal{T}\|_{\rm{HS}}=\sqrt{\sum_{i}\|\mathcal{T}\phi_{i}\|^{2}}$.

7.1 Proof of Lemma 3.1

To prove, we need to calculate the Gâteaux differential of $J$. Let $\Sigma$ and $\mathcal{T}$ be two Hilbert-Schmidt operators. We consider the two cases: $\tilde{\mathbf{\Phi}}(X_{i})\otimes\tilde{\mathbf{\Phi}}(Y_{i})-(\Sigma+\epsilon\mathcal{T})=0$ and $\tilde{\mathbf{\Phi}}(X_{i})\otimes\tilde{\mathbf{\Phi}}(Y_{i})-(\Sigma+\epsilon\mathcal{T})\neq 0$. As in Section 2.3, $\tilde{\mathbf{\Phi}}(\cdot)$’s are centered feature maps, $\zeta(\cdot)$ is a robust loss function, $\zeta^{\prime}(\cdot)$ is the derivative of $\zeta(\cdot)$, and $\varphi(t)=\frac{\zeta^{\prime}(t)}{t}$.

Case 1: $\tilde{\mathbf{\Phi}}(X_{i})\otimes\tilde{\mathbf{\Phi}}(Y_{i})-(\Sigma+\epsilon\mathcal{T})\neq\mathbf{0}$

[TABLE]

[TABLE]

Case 2: $\tilde{\mathbf{\Phi}}(X_{i})\otimes\tilde{\mathbf{\Phi}}(Y_{i})-(\Sigma+\epsilon L)=\mathbf{0}$

[TABLE]

For $\mathcal{T}=\mathbf{0}$ and $\mathcal{T}\neq\mathbf{0}$ the above equation is equal to $\lim_{\delta\to 0}\frac{\zeta(0)}{\delta}$ and $\frac{\zeta(\|\delta\mathcal{T}\|_{\rm{HS}})}{\delta}\cdot\|\mathcal{T}\|_{\rm{HS}}$, receptively. Using the assumption (i) we have $\frac{\partial}{\partial\epsilon}\zeta(\|\tilde{\mathbf{\Phi}}(X_{i})\otimes\tilde{\mathbf{\Phi}}(Y_{i})-(\Sigma+\epsilon\mathcal{T})\|_{\rm{HS}})=0$. Since $\tilde{\mathbf{\Phi}}(X_{i})\otimes\tilde{\mathbf{\Phi}}(Y_{i})-(\Sigma+\epsilon\mathcal{T})=\mathbf{0}$ and $\varphi(0)$ is well-defined by the assumption (ii). Combining Eq. (7.1) and Eq. (7.1) we get

[TABLE]

Now it is clear that for any $\Sigma,\mathcal{T}\in\rm{HS}$,

[TABLE]

Therefore,

[TABLE]

The necessary condition for $\Sigma$ to be a minimizer of $J$, i.e., $\Sigma=\widehat{\Sigma}_{RXY}$, is that $\delta J(\Sigma;\mathcal{T})=0\,\forall\,\mathcal{T}\in\rm{HS}$, which leads to $S(\Sigma)=\mathbf{0}.$

7.2 Proof of Theorem 3.2

Using Lemma 3.1 we have $\varphi(\|\tilde{\mathbf{\Phi}}(X_{i})\otimes\tilde{\mathbf{\Phi}}(Y_{i})-\widehat{\Sigma}_{RXY})\|_{\rm{HS}})\cdot(\tilde{\mathbf{\Phi}}(X_{i})\otimes\tilde{\mathbf{\Phi}}(Y_{i})-\widehat{\Sigma}_{RXY}))=\mathbf{0}$ By solving $\widehat{\Sigma}_{RXY}$, we get $\widehat{\Sigma}_{RXY}=\sum_{i=1}^{n}w_{i}\tilde{\mathbf{\Phi}}(X_{i})\otimes\tilde{\mathbf{\Phi}}(Y_{i})$, where

[TABLE]

Since the function $\zeta$ is non-deceasing, $w_{i}\geq 0$ and $\sum_{i=1}^{n}w_{i}=1$.

7.3 Proof of Theorem 3.3

We will prove this theorem in three steps; (i) the monotone decreasing property of $J(\Sigma^{(h)})$, (ii) every limit point $\Sigma^{*}$ of $\{\Sigma^{(h)}\}_{h=1}^{\infty}\in U$, and (iii) by contradiction . We define a function

[TABLE]

where $c$ is a real number. As shown in [15], for non increasing $\varphi$, the function $p$ is a surrogate function of $\zeta$ with the following two properties

[TABLE]

and

[TABLE]

Now we define a bivariate function

[TABLE]

which is a continuous function in both arguments because both $\zeta^{\prime}$ and $\varphi$ are continuous functions.

Step (i): using Eq. (39) and Eq. (40), we have

[TABLE]

and

[TABLE]

Now,

[TABLE]

From Eq. (7.3), Eq. (7.3), and Eq. (7.3), we have

[TABLE]

Thus, it is proved that $J(\Sigma^{(h)})$ monotonically decreases at every iteration. In addition, since $J(\Sigma^{(h)})\geq 0$, for any $h\geq 1$ (the sequence is bounded below at 0), it converses.

Step (ii): as in [17], it is clear that $\{\Sigma^{(h)}\}_{h=1}^{\infty}$ has a convergent subsequence $\{\Sigma^{(h_{\ell})}\}_{\ell=1}^{\infty}$. Let $\Sigma^{*}$ be the limit of $\{\Sigma^{(h_{\ell})}\}_{\ell=1}^{\infty}$. Using Eq. (7.3), Eq. (7.3), Eq. (7.3), and the monotone decreasing property of $J(\Sigma^{(h)})$, we have also get

[TABLE]

Taking the limit on both sides of the above inequality, we have

[TABLE]

Therefore,

[TABLE]

and thus

[TABLE]

This implies $\Sigma^{*}\in U.$

Step (iii): suppose $\inf_{\Sigma\in U}\|\Sigma^{(h)}-\Sigma\|_{\rm{HS}}$ does not tend to [math]. Then there exists $\epsilon>0$ such that $\forall I\in\mathbb{N},\exists h>I$ with $\inf_{\Sigma\in U}\|\Sigma^{(h)}-\Sigma\|_{\rm{HS}}\geq\epsilon$. Thus we are able to regard a increasing sequence of indices such that the $\inf_{\Sigma\in U}\|\Sigma^{(h_{\ell})}-\Sigma\|_{\rm{HS}}\geq\epsilon$ for all $\ell=1,2,\cdots\cdots$. Since $\Sigma^{(h_{\ell})}$ lies in the compact subset of $\rm{HS}$, it has a subsequence converging to some $\Sigma^{\dagger}$, and we can choose $j$ such that $\|\Sigma^{(h_{j})}-\Sigma^{\dagger}\|_{\rm{HS}}<\epsilon/2$. Since $\Sigma^{\dagger}$ is also a limit point of $\{\Sigma^{(h)}\}_{h=1}^{\infty}$, $\Sigma^{\dagger}\in U$. This is a contradiction because

[TABLE]

7.4 Proof of Theorem 5.4

We present the derivation of the IF of standard kernel CCA in detail. Recall the generalized eigenvalue problem in Eq. (23). We can formulate this problem as a simple eigenvalue problem. Using the j-th eigenfunction of the first equation of Eq. (23) we have

[TABLE]

where $e_{jX}=\Sigma_{XX}^{\frac{1}{2}}f_{jX}$.

To establish the IF of kernel CCA, we convert the generalized eigenvalue problem of kernel CCA into a simple eigenvalue problem. Henceforth, we can use the results such as the IF of linear PCA analysis [48, 49], the IF of linear CCA [25] and the IF of kernel PCA (finite dimension and infinite dimension) [18], respectively. To do this, first we calculate the IF of $\Sigma_{XX}^{-\frac{1}{2}}\Sigma_{XY}\Sigma_{YY}^{-1}\Sigma_{YX}\Sigma_{XX}^{-\frac{1}{2}}$, henceforth, the IF of $\Sigma_{YY}^{-1},\Sigma_{XX}^{\frac{1}{2}}$ and $\Sigma_{XY}$.

Using simple algebra (as in Section4.2), we have the IF for the following operators at point $Z^{\prime}=(X^{\prime},Y^{\prime})$:

[TABLE]

For simplicity, let us define $\tilde{k}_{X}(\cdot,X^{\prime}):=k_{X}(\cdot,X^{\prime})-\mathcal{M}_{X},\,\tilde{k}_{Y}(\cdot,\mathbf{y}^{\prime}):=k_{Y}(\cdot,Y^{\prime})-\mathcal{M}_{Y}$. Also define $\mathbb{A}:=\Sigma_{XY}\Sigma_{YY}^{-1}\Sigma_{YX},\,\mathbb{B}:=\Sigma_{XX}^{-\frac{1}{2}}\mathbb{A}\Sigma_{XX}^{-\frac{1}{2}}$, and $\mathbb{L}=\Sigma_{XX}^{-\frac{1}{2}}(\Sigma_{XX}^{-\frac{1}{2}}\Sigma_{XY}\Sigma_{YY}^{-1}\Sigma_{YX}\Sigma_{XX}^{-\frac{1}{2}}-\rho^{2}\mathbf{I})^{-1}\Sigma_{XX}^{-\frac{1}{2}}$. Now

[TABLE]

Then,

[TABLE]

and

[TABLE]

The influence of $\mathbb{B}$ is then given by

[TABLE]

To define the IF of kernel CC ($\rho_{j}^{2}$) and kernel CVs ($f_{X}(X)$ and $f_{Y}(Y)$), we convert a generalized eigenvalue problem and use the Lemma 1 of [18] and Lemma 2 of [49]. Then the IF of kernel $\rho_{j}^{2}$ is defined as follows

[TABLE]

For simplicity, Eq. (49) can calculate in parts. The first part is derived as

[TABLE]

In the last equality, we used Eq. (7.4). The 2nd part of Eq. (49) is derived as

[TABLE]

In the last second equality, we used Eq. (23). Similarly we can write the 3rd term as

[TABLE]

where $\tilde{f}_{jX}(X^{\prime})=\langle f_{jX},\tilde{k}_{X}(\cdot,X^{\prime})\rangle$ and similar for $\tilde{f}_{jY}$. Therefore, substituting Eq. (50), (51) and (52) into Eq. (49), the IF of kernel CC is given by

[TABLE]

Now we derived the IF of kernel CVs. To this end, first we need to derive

[TABLE]

By the first term of Eq. (54) we have

[TABLE]

We derive each term of Eq. (55), respectively. The first term of Eq. (55) is given by

[TABLE]

The 2nd term of Eq. (55) is

[TABLE]

and the 3rd term of Eq. (55) is

[TABLE]

By substituting the above three equations into Eq. (55), we have

[TABLE]

The 2nd term of Eq. (54) is given

[TABLE]

Therefore, substituting Eq. (57) and Eq. (58) into Eq. (54) we get the j-th IF of kernel CV for the $X$ data:

[TABLE]

Similarly, for the $Y$ data we have,

[TABLE]

Appendix B: Abbreviations and symbols

Table 5 and 6 preset the list of abbreviations, and symbols, respectively.