Stein Shrinkage for Cross-Covariance Operators\\ and Kernel Independence Testing
Aaditya Ramdas$^{*
Cross-covariance operators arise naturally in many applications using Reproducing Kernel Hilbert Spaces (RKHSs) and are typically estimated using an empirical plugin estimator, which we demonstrate are poor estimators of operator (eigen)spectra at low sample sizes. This paper studies the phenomenon of Stein shrinkage for infinite dimensional cross-covariance operators in RKHSs, as briefly initiated by Muandet et al (2014) who recently suggested two shrinkage estimators. We develop a third family of shrinkage estimators and undertake a study of how shrinkage improves estimation of operator spectra. We demonstrate an important and surprising application, that shrunk test statistics yield higher power for kernel independence tests and we provide insights into why they improve performance.
arxiv.tex
\documentclass{article}


\usepackage[left=1in,top=1in,right=1in,bottom=1in,head=1in]{geometry}
\usepackage{amsfonts,amsmath,amssymb,amsthm}
\usepackage{verbatim,float,url}
\usepackage{graphicx,subfigure,psfrag}
\usepackage{dsfont,bm}
\usepackage{natbib}

\newtheorem{algorithm}{Algorithm}
\newtheorem{theorem}{Theorem}
\newtheorem{lemma}{Lemma}
\newtheorem{corollary}{Corollary}

\theoremstyle{remark}
\newtheorem{remark}{Remark}
\theoremstyle{definition}
\newtheorem{definition}{Definition}

\newcommand{\argmin}{\mathop{\mathrm{argmin}}}
\newcommand{\argmax}{\mathop{\mathrm{argmax}}}
\newcommand{\st}{\mathop{\mathrm{subject\,\,to}}}

\newcommand{\Inocco}{\textrm{I}^{\textrm{NOCCO}}}
\newcommand{\fix}{\marginpar{FIX}}
\newcommand{\new}{\marginpar{NEW}}
\newcommand{\SigYX}{\Sigma_{XY}}
\newcommand{\SigXY}{\Sigma_{XY}}
\newcommand{\SXY}{S_{XY}}
\newcommand{\SXX}{S_{XX}}
\newcommand{\SigXX}{\Sigma_{XX}}
\newcommand{\SYY}{S_{YY}}
\newcommand{\SigYY}{\Sigma_{YY}}
\newcommand{\SYX}{S_{YX}}
\newcommand{\nep}{\epsilon}
\newcommand{\crossx}{\frac{1}{n}\sum_i \phi_{x_i}\otimes\varphi_{y_i}}
\newcommand{\crossy}{\frac{1}{n}\sum_i \varphi_{y_i}\otimes\phi_{x_i}}
\newcommand{\Hs}[1]{\mathcal{H}_{#1}}
\newcommand{\E}[1]{\mathbb{E}_{#1}}
\newcommand{\ex}[2]{\mathbb{E}_{#1}[#2]}
\newcommand{\exs}[2]{\mathbb{E}^2_{#1}[#2]}
\newcommand{\HSIC}{\textrm{HSIC}}
\newcommand{\tr}{\textrm{tr}}
\newcommand{\CX}{C_{XX}}
\newcommand{\tL}{\tilde L}
%\newcommand{\L}{\textrm{ L}}
\newcommand{\tK}{\tilde K}
\newcommand{\tM}{\tilde M}
%\newcommand{K}{\textrm{ K}}
\newcommand{\I}{\textrm{I}}
\newcommand{\KLW}{\tilde \K_{LW}}
\newcommand{\X}{\mathcal{X}}
\newcommand{\Y}{\mathcal{Y}}
%\newcommand{\ck}{\mathcal{k}}
%\newcommand{\l}{\mathcal{l}}
\newcommand{\tphi}{\tilde \phi}
\newcommand{\tpsi}{\tilde \psi}


\def\E{\mathbb{E}}
\def\P{\mathbb{P}}
\def\Cov{\mathrm{Cov}}
\def\Var{\mathrm{Var}}
\def\half{\frac{1}{2}}
\def\th{\mathrm{th}}
\def\tr{\mathrm{tr}}
\def\df{\mathrm{df}}
\def\dim{\mathrm{dim}}
\def\col{\mathrm{col}}
\def\row{\mathrm{row}}
\def\nul{\mathrm{null}}
\def\rank{\mathrm{rank}}
\def\nuli{\mathrm{nullity}}
\def\sign{\mathrm{sign}}
\def\supp{\mathrm{supp}}
\def\diag{\mathrm{diag}}
\def\aff{\mathrm{aff}}
\def\hy{\hat{y}}
\def\ty{\tilde{y}}
\def\hbeta{\hat{\beta}}
\def\tbeta{\tilde{\beta}}
\def\htheta{\hat{\theta}}
\def\halpha{\hat{\alpha}}
\def\hu{\hat{u}}
\def\hf{\hat{f}}
\def\lone{1}
\def\ltwo{2}
\def\linf{\infty}
\def\lzero{0}
\def\T{^\top}
\def\R{\mathbb{R}}
\def\cA{\mathcal{A}}
\def\cB{\mathcal{B}}
\def\cD{\mathcal{D}}
\def\cE{\mathcal{E}}
\def\F{\mathcal{F}}
\def\G{\mathcal{G}}
\def\H{\mathcal{H}}
\def\cM{\mathcal{M}}
\def\cN{\mathcal{N}}
\def\cS{\mathcal{S}}
\def\cT{\mathcal{T}}
\def\cW{\mathcal{W}}
\def\cZ{\mathcal{Z}}
\def\TV{\mathrm{TV}}

\title{Stein Shrinkage for Cross-Covariance Operators\\
and Kernel Independence Testing
}




\author{ 
Aaditya Ramdas$^{*}$ and Leila Wehbe\footnote{Equal contribution from both student authors.} \\ 
{Machine Learning Department,}\\
 {Carnegie Mellon University}\\
 {Pittsburgh, PA 15213, USA}\\
 \texttt{aramdas,lwehbe@cs.cmu.edu}
 }

\date{}

\begin{document}

\maketitle

\begin{abstract}
Cross-covariance operators arise naturally in many applications using Reproducing Kernel Hilbert Spaces (RKHSs) and are typically estimated using an empirical plugin estimator, which we demonstrate are poor estimators of operator (eigen)spectra at low sample sizes. This paper studies the phenomenon of Stein shrinkage for infinite dimensional  cross-covariance operators in RKHSs, as briefly initiated by  Muandet et al (2014) who recently suggested two shrinkage estimators. We develop a third family of shrinkage estimators and undertake a study of how shrinkage improves estimation of operator spectra. We demonstrate an important and \textit{surprising} application, that \textit{shrunk} test statistics yield higher power for kernel independence tests and we provide insights into why they improve performance. 
\end{abstract}
\section{Introduction}
\label{sec:intro}

Reproducing Kernel Hilbert Spaces (RKHSs) have become an indispensable part of machine learning. The long list of applications includes kernel PCA \cite{kernelpca}, dimension reduction \cite{kerneldimred}, structure discovery \cite{kernelstructdisc}, independence testing \cite{hsic},  Granger causality in time series \cite{kernelgranger}, learning graphical models \cite{graphicalmodels}, feature learning \cite{kernelfeature}, target or conditional shift \cite{kerneltargetshift} and kernel Bayes' rule \cite{kernelbayesrule}. \textit{All} the above methods rely on empirical covariance and cross-covariance operators. We derive, and analyze in detail, better estimators for these operators (than their corresponding plugin estimators) and demonstrate their usefulness in independence testing.
 

It is well known that the sample covariance \textit{matrix} is a poor estimator for estimating finite dimensional covariance \textit{matrices} from small samples, i.e. it results in highly inaccurate eigenvalues/vectors. Over the decades, there have been many estimators proposed for covariance matrices that outperform the sample covariance matrix in different situations - see \cite{ledoitwolf,covmat1,covmat2,covmat3}. The observation that led to the development of these estimators is a special case of a phenomenon known as the Stein effect \cite{stein} which has long been studied in the classical parametric statistics literature. The Stein estimator was originally used for the problem of estimating the mean of multivariate Gaussian distributions: it was shown that shrinking the sample mean towards any fixed point always had lower mean squared error than the sample mean. 

This phenomenon was only very recently proved to  occur in RKHSs by  \cite{kernelstein} who show that the sample kernel mean  estimator is inadmissible, i.e. there always exist an estimator that outperforms it. 
The cross-covariance operator is a specific case of a kernel mean estimator (see Appendix or \cite{kernelstein}) and hence its sample estimator is also inadmissible. This motivates us to find better estimators (possibly task-specific) that improve performance across the wide variety of aforementioned tasks. A first step was taken by \cite{kernelstein}, who extend their techniques for a more accurate kernel mean map to also mention two shrinkage estimators for cross-covariance operators. However, these shrunk operators were not given much special attention in their theoretical and experimental treatment of the topic. Here, we will introduce a third novel shrinkage estimator, and perform an in-depth analysis of all three, with special focus on their ability to improve power of independence testing.


\newpage
\textbf{A Novel Ledoit-Wolf Style (LW) Estimator}
The Ledoit-Wolf estimator for \textit{covariance matrices} \cite{ledoitwolf} is asymptotically optimal for estimation with respect to the quadratic risk (Frobenius norm, sum of squared singular values). We extend the techniques from \cite{ledoitwolf} into our more general infinite dimensional setting utilizing the fact that \textit{cross-covariance operators} are Hilbert-Schmidt (sum of squared singular values is finite). 
The shrinkage intensity is obtained in closed form and is quite interpretable - it reflects the relative improvement in average loss, it decreases with increasing sample size and if the variance of the sample covariance operator is small, then the shrinkage intensity is low. It also demonstrates the possibility of extending other shrinkage ideas from the rich covariance matrix literature.

\textbf{Kernel Independence Testing with HSIC} We introduce a potentially important application for covariance operator shrinkage, that of Kernel Independence Testing with HSIC. Kernel Independence Testing using the Hilbert-Schmidt Independence Criterion (HSIC) was recently  introduced \cite{hsic}. The HSIC is the Hilbert-Schmidt norm of the cross-covariance operator, and it is proven to be equal to 0 if the variables $x$ and $y$ are independent. We will see that the effect of shrinkage on HSIC value turns out to be similar to a soft-thresholding operator. 
It is not obvious, and possibly even counter-intuitive, that shrinking cross-covariance operators (and hence the test statistic HSIC) towards zero  could  improve the power of the test. 
However, if the shrinkage was monotone and solely a function of HSIC (as the soft-thresholding operator is), then the quantile of the test statistic during permutation testing would remain unchanged, leading to no improvement in power. Our expressions for shrunk HSIC reveal that this is in fact not the case. Experimentally, we demonstrate that the null distribution (HSIC for permuted samples) gets shrunk \textit{more} than the alternate distribution, and hence we often (but not always) see an increase in power but never a decrease. While experiments suggest that improvements with Laplace kernels are more than for Gaussian kernels, we leave a formal study of how improvement depends on the kernel spectral decay for future work. 


\textbf{Simple/Flexible COvariance Shrinkage Estimators (SCOSE, FCOSE)}
SCOSE and FCOSE were derived in \cite{kernelstein} from a related but slightly different viewpoint from \cite{ledoitwolf}, of minimizing a quadratic reconstruction loss with Tikhonov regularization on different objects. We give an alternate interpretation to the original LW estimator by noting that its expression is very similar to the SCOSE closed form expression with leave one out cross validation when the number of samples is large. FCOSE is more computationally intensive and quite different from LW/SCOSE as it is a non-linear shrinkage method. While SCOSE and LW have demonstrably lower quadratic risk, we show that FCOSE estimates singular functions better than LW, SCOSE, and hypothesize that this gives it an advantage in Kernel PCA and HSIC.


\paragraph{Notation}
Assume $\X$ and $\Y$ are measurable domains,  and $\F$ and  $\G$ are RKHSs defined on them, i.e. $\F$ is the set of functions $f: \X\rightarrow \R $ corresponding to a reproducing kernel $k : \X \times \X \rightarrow \R$ and feature map at each $x \in \X$ denoted by $\phi_x : \X \rightarrow \R$, where $\phi_x(y) = k(x,y)$. Similarly, $\G$ is the set  $g: \Y\rightarrow \R $ corresponding to kernel  $\ell : \Y \times \Y \rightarrow \R$ and and feature map   $\varphi_y : \Y \rightarrow \R$ for all $y \in \Y$. Given a joint  distribution $P_{XY}$ on $\X \times \Y$ with marginals $P_X,P_Y$, the cross-covariance operator   $\SigYX:\Y \rightarrow \X$ is defined  as:
 $$\SigXY := \E_{P_{XY}}[\phi_x \otimes \psi_y] - \mu_{P_X}\otimes \mu_{P_Y}$$
where $\otimes$ represents a tensor/outer product and $\mu_{P_X} : \X \rightarrow \R$  is the kernel mean map defined as $\mu_{P_X} := \E_P[\phi_x]$ whose empirical version is $\hat{\mu}_P := \frac1{n}\sum_{i=1}^n \phi_{x_i}$ (similarly $\mu_{P_Y} := \E_{P_Y}[\varphi_y]$,  $\hat{\mu}_Q := \frac1{n}\sum_{i=1}^n \varphi_{y_i}$).  The covariance operator $\SigXX:\X \rightarrow \X$ is defined as $\SigXX := \E_{P_{X}}[\phi_x \otimes \phi_x] - \mu_{P_X}\otimes \mu_{P_X}$. The most common estimator of $\SigYX $ is the  sample cross-covariance operator, a plugin estimator obtained by using the empirical joint and marginal distributions,
{\small $$\SXY := \frac1{n}\sum_{i=1}^n \phi_{x_i}\otimes \varphi_{y_i} - \hat{\mu}_P \otimes \hat{\mu}_Q.$$} \\

To illustrate that sample mean operators perform poorly at low sample sizes, we show two extremely varied examples where  the spectra of the sample operators are drastically far from the truth (Appendix Fig.\ref{fig:motivateApp}). Indeed, it shows that the problem (of plug-in operators being poor estimators) does indeed exist across a wide variety of kernels, types of operators, input distributions and underlying data dimensionality. This highlights the importance of the study of Stein shrinkage for kernel operators as taken up in this paper.

The paper outline is as follows. Section 2 introduces a new shrinkage estimator, and compares it to existing ones. Section 3 discusses the application of shrunk estimators to independence testing. Section 4 wraps up by showing promising future directions. 


\section{Ledoit-Wolf Style Shrinkage Estimation}

We now develop a new family of linear shrinkage estimators for the cross covariance operator $\SigXY$ (estimators for $\SigXX$ follows naturally) by extending a now classical technique described by \cite{ledoitwolf}. 
 

We derive the optimal (oracle) shrinkage towards the zero operator \textbf{0} (which maps each RKHS function to the zero function).
A common shrinkage target for covariance matrices is the identity matrix - it seems like a natural generalization would be to shrink the covariance operator towards the identity operator (which maps each function in the RKHS to itself). However, the identity operator is \textit{not} Hilbert-Schmidt (HS) - all its eigenvalues are unity, and hence the sum of squares is infinite - making any linear/convex combination with it also not Hilbert-Schmidt. We discuss another possible target in the future work section.

The oracle $\SXY^*$ is derived by finding the convex combination of $\SXY$ and $\mathbf{0}$ that minimizes quadratic risk:
\begin{eqnarray}
\SXY^* &:=& \arg\min_S \E\|S - \SigXY\|_{HS}^2  \label{eq:LW}\\
&& ~\mbox{ subject to }~ S = (1-\rho) \SXY + \rho \ \mathbf{0}.
\end{eqnarray}
First define the constants $\alpha^2 = \|\SigXY\|_{HS}^2$, $\beta^2 = \E\|\SXY - \SigXY\|_{HS}^2$ and $\delta^2 = \E\|\SXY\|_{HS}^2$. Since $\E[\SXY] = \SigXY$, 
 it is easy to verify that
$
\alpha^2  + \beta^2 = \delta^2.
$
Substituting and expanding the objective of the optimization problem, we get:
\begin{eqnarray}
\E\|S - \SigXY\|_{HS}^2 &=& \E\|-\rho \SXY + (\SXY - \SigXY)\|_{HS}^2 = \rho^2 \delta^2 + \beta^2  - 2\rho (\delta^2 - \alpha^2)\label{eq:bv}\\
&=& \rho^2 \alpha^2 + (1-\rho)^2 \beta^2 \nonumber.
\end{eqnarray}
Differentiating and equating to zero, the optimal value of $\rho$ is $\rho ~=~ \frac{\beta^2}{\delta^2}$, yielding\footnote{As often seen in Stein Shrinkage, the estimator $\left( 1 -\frac{\beta^2}{\delta^2} \right)_+$ is used to avoid negative multiples.}
$$
\SXY^* ~=~ \left( 1 -\frac{\beta^2}{\delta^2} \right) \SXY.
$$
The oracle operator $\SXY^*$ has the drawback that it depends on unknown constants $\beta^2,\delta^2$, and we now describe a bonafide estimator by suggesting simple sample estimators for these. $\delta^2$ has the simple plugin estimator
\begin{equation}\label{eq:d2}
d^2 ~=~ \|\SXY\|_{HS}^2 ~=~ \frac1{n^2}\sum_{i,j=1}^n \tK_{ij}\tL_{ij},
\end{equation}
where  $\tK = HKH$, $\tL = HLH$  are centered gram matrices.  $K,L$ are the $n\times n$ kernel matrices corresponding to $k,l$  i.e. $K_{ij}:=k(x_i,x_j), (y_i,y_j):=l(y_i,y_j)$ and $H = I - \frac1{n}\mathbf{11}^\top$, where $\mathbf{1}$ is the vector $[1, 1, ..., 1]^\top$. 

$\beta^2$ is the variance of the sample covariance operator, which can in turn be estimated by the sample variance of the sample covariance operator (note the extra factor of $1/n$ since the variance of the sample mean of $\tphi_{x_i} \otimes \tphi_{x_i}$ is $1/n$ times the variance of $\tphi_{x_i} \otimes \tphi_{x_i}$):
\begin{eqnarray}
b^2 ~=~ \frac{1}{n}\frac1{n} \sum_{k=1}^{n} || \tphi_{x_i} \otimes \tpsi_{x_i} - \SXY||^2 ~=~ \frac1{n} \left[ \frac{1}{n} \sum_{i=1}^{n}  \tK_{ii}\tL_{ii}  - \frac{1}{n^2} \sum_{i,j=1}^{n} \tK_{ij}\tL_{ij} \right].
\end{eqnarray}
Then we get a bona fide Ledoit-Wolf estimator that can be calculated from data:
 $$\SXY^{LW} := \left( 1 - \rho^{LW} \right ) \SXY ~\mbox{where}~ \rho^{LW} = \frac{b^2}{d^2}.$$
 We kept our notation consistent with Ledoit and Wolf \cite{ledoitwolf} to make our results visually similar to theirs. However, our constants $\alpha^2,\delta^2$ (and their corresponding estimators) are mathematically different from theirs due to a different shrinkage target and the use of kernel gram matrices. Note the use of centered kernel matrix $\tK$ and not $K$ in our formulae. The shrinkage intensity depends both on the kernel being used and on the underlying distribution of points, since the sample and true covariance operator themselves depend on these. However, it is presently unclear qualitatively how these quantities (eg: the decay of kernel eigenvalues) affect the shrinkage intensity and needs further investigation.
 
The shrinkage intensity $\beta^2/\delta^2$ is  precisely the relative improvement in average loss, i.e.
$$
\frac{\E\|\SXY - \SigXY\|_{HS}^2 - \E\|\SXY^* - \SigXY\|_{HS}^2}{\E\|\SXY - \SigXY\|_{HS}^2} = \frac{\beta^2}{\delta^2}.
$$
Intuitively, if the variance of the sample covariance operator is small, then we wouldn't need to shrink (since $\SXY$ would be both unbiased and have low variance). This is exactly the quantity $\beta^2$ in the numerator.

An important interpretation of the shrunk estimator is achieved via the bias-variance decomposition, since Eq.\eqref{eq:bv} can be alternately be expanded as
\begin{eqnarray}
\E\|S - \SigXY\|_{HS}^2 &=&  \|\E[S] - \SigXY\|_{HS}^2 + \E\|S-\E[S]\|_{HS}^2 = \mbox{bias}^2 + \mbox{variance}.
\end{eqnarray}
While $S=\SXY$ has zero bias and high variance, target $S=\mathbf{0}$ has high bias and zero variance, and $S=\SXY^*$ represents the optimal bias-variance tradeoff between the two extremes (a common idea in Stein shrinkage literature \cite{stein}). Below, we see an example of the singular values before and after shrinkage. Notice how, after shrinkage, the singular values are brought closer to the ``true" singular values (approximated with $n=5000$).

 

\begin{figure} [h]
\centering
\includegraphics[width=0.7\linewidth]{LW_example.pdf}
\caption{\small Leftmost: Plots of the dataset (hollow Gaussian) for different $N$. Middle: Unshrunk eigenvalues of $\SXY$. Right: Shrunk eigenvalues of $\SXY^{LW}$ with corresponding $\rho^{LW}$ given in legend. The purple line (5000) remains unchanged, while all the others are shrunk towards zero and are closer to the purple (``truth''), thought the light blue (n=20) is a bit overshrunk (which is expected at small $n$).}
\label{fig:LWeg}
\end{figure}
\subsection{Comparing LW to SCOSE and FCOSE}
We summarise some important details of shrunk estimators derived by \cite{kernelstein}, before comparing all three estimators qualitatively 
and quantitatively through detailed simulations.

\paragraph{SCOSE}
Since the sample covariance operator can be derived as  the solution to the optimization problem, $\min_{S \in \mathcal{HS}(\F,\G)} \frac1{2} \sum_{i=1}^n \left \|\tphi_{x_i}\otimes \tpsi_{y_i} - S \right\|_{HS}^2$, one can derive a regularized estimator by considering
$$
\min_{S \in \mathcal{HS}(\F,\G)} \frac1{2} \sum_{i=1}^n \left \|\tphi_{x_i}\otimes \tpsi_{y_i} - S \right\|_{HS}^2 + \lambda \|S\|_{HS}^2
$$
and the optimal solution  (called Simple COvariance Shrinkage Estimator) is
$$
\SXY^{S} := \left(1 - \frac{\lambda}{1+\lambda}\right) \SXY, %=: (1-\rho^S) \SXY
$$
where $\lambda$ (or the shrinkage intensity) is estimated by leave-one-out cross-validation as 

$$
\rho^S ~:=~ \left(\frac{\lambda^{CV}}{1+\lambda^{CV}}\right) ~=~ \frac{\left[ \frac{1}{n} \sum_{i=1}^{n}  \tK_{ii}\tL_{ii}  - \frac{1}{n^2} \sum_{i,j=1}^{n} \tK_{ij}\tL_{ij} \right]}{(n-2)\frac{1}{n^2} \sum_{i,j=1}^{n} \tK_{ij}\tL_{ij} +  \frac1{n} \frac{1}{n} \sum_{i=1}^{n}  \tK_{ii}\tL_{ii} }.
$$
SCOSE bears similarity to LW - by multiplying numerator and denominator by $1/n$, $\rho^S$ can be rewritten as
$$
\rho^S = \frac{b^2}{\frac{n-2}{n}d^2 + \frac{d^2}{n^2} + \frac{b^2}{n}},
$$ 
which looks very similar to $\rho^{LW} = b^2/d^2$ when $n$ is large (though both $\rho$s will be near zero then). Indeed, we shall see that they have very similar performance on most tasks. Hence SCOSE can be seen as an alternative way to interpret or derive the LW estimator, even for covariance matrices.

A subtle issue occurs when $b^2 > d^2$ and $(1 - \rho^{LW})_+ = 0$. This is often encountered at small $n$ when
$\frac1{n}\sum_{i=1}^n \tK_{ii}\tL_{ii} > (n+1) \frac1{n^2}\sum_{i,j=1}^n \tK_{ij}\tL_{ij}$. 
In this case, one can make the  observation that $\frac{n-2}{n}d^2 + \frac{d^2}{n^2} + \frac{b^2}{n} ~<~ \frac{n-2}{n}b^2 + \frac{b^2}{n} + \frac{b^2}{n} ~<~ b^2$ causing $\left(1-\rho^S\right) <0$. As with LW we should instead rightly use $(1-\rho^S)_+$, and we get 
$\SXY^{LW} = \SXY^S = 0.$ 
An alternate option
is to calculate $\lambda^{CV}$ to minimize LOOCV and if $\lambda^{CV} < 0$, to set $\lambda^{CV} =0$. This corresponds to setting  $\SXY^{S}$ to $\SXY$ and not to {\bf 0}, and is incorrect since  $\lambda^{CV} < 0$ when $\SXY$  is too unreliable. 


\paragraph{FCOSE}
By the representer theorem, the sample covariance operator can alternately be derived by solving $\min_{\beta \in \R^n} \frac1{2} \sum_{i=1}^n \left \|\tphi_{x_i}\otimes \tpsi_{y_i} - \sum_{i=1}^n \frac{\beta_i}{n} \tphi_{x_i}\otimes \tpsi_{y_i} \right\|_{HS}^2$, one can derive a regularized estimator as follows:
$$
\min_{\beta \in \R^n} \frac1{2} \sum_{i=1}^n \left \|\tphi_{x_i}\otimes \tpsi_{y_i} - \sum_{i=1}^n \frac{\beta_i}{n} \tphi_{x_i}\otimes \tpsi_{y_i} \right\|_{HS}^2 + \lambda \|\beta\|_2^2
$$
and the optimal solution (called Flexible COvariance Shrinkage Estimator) is given by 
\begin{eqnarray*}
\SXY^{F} &:=& \sum_{i=1}^n \frac{\beta_i}{n} \tphi_{x_i}\otimes \tpsi_{y_i} ~\mbox{where}~ \beta = (\tK \circ \tL + \lambda I)^{-1} \tK \circ \tL  \ \mathbf{1}.   \ \ \  ~\mbox{ ($\circ $ indicates  Hadamart product)}
\end{eqnarray*}
$\lambda$ is determined by leave one out cross-validation (LOOCV) -
a single eigenvalue decomposition of $\tK\circ\tL$ costing $O(n^3)$ can be done, following which evaluating LOOCV is only $O(n^2)$ per $\lambda$ (see \cite{kernelstein} for details). 
Even  though FCOSE is  more computationally intensive compared to SCOSE and LW due to cross-validation, since shrinkage is really useful at low sample sizes, it is probably the case that by the time $n$ gets large enough for eigenvalue decompositions to become prohibitively expensive, there is already little need for shrinkage.
Note FCOSE applies a nonlinear transformation to $\SXY$, while LW and SCOSE are linear shrinkage methods.

\paragraph{Quadratic Risk} Refer to Fig. 2 for a comparison of LW, SCOSE and FCOSE to unshrunk operators. Appendix Fig. \ref{fig:esterrorApp} compares two more kernels (linear, poly) where SCOSE, LW heavily outperform FCOSE.

\begin{figure} [h]
\centering
\includegraphics[width=0.45\linewidth]{extin_Laplace.pdf}
\includegraphics[width=0.45\linewidth]{extin_rbf.pdf}
\caption{\small Both panels show the quadratic risk $\E\|S-\SigXY\|_{HS}^2$ for $S\in\{\SXY,\SXY^{LW},\SXY^{S},\SXY^{F}\}$ for Laplace (left) and Gaussian (right) kernels. The same dataset as shown in Fig.\ref{fig:LWeg} was used. $\SigXY$ was approximated from 5000 points, and $\E$ was approximated by 1000 experiment repetitions. Takeaways: The improvement for Laplace is more than for Gaussian, the relative gain goes down with $n$ and LW/SCOSE outperform FCOSE.}
\label{fig:esterror}
\end{figure}


\section{Kernel Independence Testing}

One of the most important recent uses of the cross-covariance operator is in (kernel) independence testing, where we have samples $(x_i,y_i)$ from a joint distribution $P_{XY}$ and we want to test 
 $$ 
 H_0: P_{XY} = P_XP_Y ~\mbox{~ against ~}~ H_1 : P_{XY} \neq P_XP_Y
 $$ 
 \cite{hsic} introduced the test statistic Hilbert Schmidt Independence Criterion, $\HSIC=\|\SigXY\|_{HS}^2$, estimated as
$$
\HSIC_n = \|\SXY\|_{HS}^2 = \frac1{n^2}\tr(\tK\tL).
$$  
$H_0$ is rejected if the sample statistic is in the right tail of the statistic's null distribution determined by permutation testing, giving rise to a nonparametric test for independence. If a characteristic kernel (see \cite{hsic} for details) is used, then the test is consistent against any fixed alternative, i.e. $\HSIC = 0$ if and only if $X$ and $Y$ are independent. 
The experiments in \cite{hsic} show that HSIC has reasonable power against difficult alternatives. We will show that a better estimator of $\SigXY$ than $\SXY$ can improve the power of HSIC. This fact is not obvious - a monotone transformation of $\HSIC_n$ that is a function solely of $\HSIC_n$ will not change its quantile in permutation tests and hence cannot improve power. Indeed, we see below that the aforementioned shrinkage methods do lead to non-monotone transformations of HSIC.

\begin{figure} [h]
\centering
\includegraphics[width=0.45\linewidth]{scatter_hsic.pdf}
\hspace{0.4in}
\includegraphics[width=0.25\linewidth]{ratios.pdf}
\caption{\small  The scatterplots compare $\HSIC$ with $\HSIC^{LW}$ (left), $\HSIC^{S}$ (middle), $\HSIC^{F}$ (right). The top row corresponds to $n=20$, and the bottom row has $n=50$. Each cross mark  corresponds to the value of shrunk and unshrunk HSIC calculated during a single permutation of a permutation test. Visually, all three look like a soft-thresholding operator - if the unshrunk HSIC is small, it deems it as ``noise'' and shrinks it to zero. For larger $n$, it shrinks to zero in a smaller region, and the effect gets more monotone. Shrinking increases power because, under $H_1$, HSIC shrinks less than under $H_0$. This is demonstrated by the right panel, which shows the ratio between the mean real HSIC and the mean 95 percentile HSIC from the permutation test that is used as a threshold for rejecting $H_0$.
  In other words, for the shrunk estimators the unpermuted HSIC is differentiated better from (since it is shrunk less than)  the permuted HSICs.
   See Appendix Fig.\ref{fig:hsicscatterApp},\ref{fig:hsicpowerApp} for more examples. 
    }
\label{fig:hsicscatter}
\end{figure}


The familiar reader would have already noticed  that the $d^2$ term in $\rho^{LW}$ was nothing but HSIC. Indeed, even $b^2$ has an HSIC term lurking in it and the LW-shrunk test statistic $\|\SXY^{LW} \|_{HS}^2$ can be simplified as
\begin{eqnarray}
\HSIC^{LW}
&=& \left(1 - \frac{ \frac1{n}\sum_{i=1}^n\tK_{ii}\tL_{ii} - \HSIC_n}{n \HSIC_n} \right)^2 \HSIC_n.
\end{eqnarray}
Note that the larger the HSIC, the lesser the shrinkage. Also, while the shrinkage intensity depends on HSIC, it also depends on another term capturing the mean interaction between all elements of $K,L$ through $\tK,\tL$.
One can similarly derive expressions for HSIC when SCOSE and FCOSE are used for shrinkage:
\begin{eqnarray}
\HSIC^{S} &=& \left(1 - \frac{\frac1{n}\sum_{i=1}^n\tK_{ii}\tL_{ii} - \HSIC_n}{(n-2)\HSIC_n +  \tfrac{\frac1{n}\sum_{i=1}^n\tK_{ii}\tL_{ii}}{n} } \right)^2 \HSIC_n,\\
\HSIC^F &=& \frac1{n^2} \tr(\tK\circ\tL(\tK\circ\tL+\lambda I)^{-1}\tK\circ\tL(\tK\circ\tL+\lambda I)^{-1}\tK\circ\tL).
\end{eqnarray}
Refer to Fig.\ref{fig:hsicscatter} for visualizing the effect of shrinkage on HSIC in permutation testing (i.e. under $H_0$). Also note that power improves because the effect of shrinkage is different under $H_1$ (Fig.\ref{fig:hsicscatter} right panel, Fig.\ref{fig:hsicscatterApp}). 

\subsection{Simulations for Power}
We will draw  samples from a non-product distribution multiple (100) times, and for each test statistic we will calculate the power, which is the percentage of times it correctly rejects the null. Refer to the subtext of Fig.\ref{fig:hsicpower} for detailed analysis. We provide two more examples in the Appendix, Fig.\ref{fig:hsicpowerApp}. We do not plot SCOSE since we found that it overlaps with LW in most experiments. 
\begin{figure} [h!]
\centering
\includegraphics[width=0.85\linewidth]{HSIC2diff_extin.pdf}
\includegraphics[width=0.85\linewidth]{HSIC2diff_N_theta.pdf}
\caption{\small We present results from 100 repetitions of two independent examples (top, bottom). The leftmost panel shows the scatter plot of the distribution being used. In the first example, we fix the number of points and increase the radius of the empty circle in the middle. In the second, we increase the number of points, keeping the orientation angle fixed. The middle panels show the power of shrunk and unshrunk HSICs versus varying radius (top) or sample size (bottom). In each example, the top right panel shows how often shrunk HSIC gets it right when HSIC gets it right (correctly rejects the null) - higher is better and any drop from 1 shows a worsening of shrunk estimators. The bottom right panel shows how often shrunk HSIC gets it right when HSIC makes a mistake (incorrectly fails to reject the null) - higher is better and any rise from 0 shows an improvement for shrunk estimators. From the top right panels, one may conclude that the shrunk HSICs are almost never worse, and from the bottom right panels, we see that they can be significantly better. The middle panels show that you can get improvements in the low power as well as high power regime. Unlike for quadratic risk, FCOSE almost always outperforms LW, with LW  sometimes showing smaller improvement. }
\label{fig:hsicpower}
\end{figure}
\subsection{Understanding FCOSE}
\begin{figure} [h]
\centering
\includegraphics[width=0.75\linewidth]{eigenvectors.pdf}
\caption{\small The left panels shows 30, 1000 samples from the same distribution. In the middle panels, the ``true'' left and right singular functions $f_1^*,g_1^*$ of $\SigXY$ (Gaussian kernel) are estimated from 1000 points, and plotted in black. Then, from one example of 30 points, we plot the sample left and right singular functions, before ($f_1,g_1$, blue) and after ($f_1^F, g_1^F$, purple) shrinking. The right panels compare $f_1,f_1^F$ and $g_1,g_1^F$ for 100 different examples (100 different 30 point samples). 
The right panels show that over the 100 repeated experiments, the FCOSE error for the singular functions is on average smaller than the unshrunk error, demonstrating its better accuracy in estimating singular vectors. We verified this by calculating the MSE between the unshrunk/shrunk functions and the ``truth'', i.e. $\E\|f_1 - f^*_1\|_K^2$ and $\E\|f_1^F - f^*_1\|_K^2$. MSE of unshrunk estimators of $f$ is 0.5941 and of $g$ is 0.5704, which were larger than MSE of FCOSE for $f$ which was 0.5397, and for $g$ which was 0.5309. See Appendix for more examples. 
}
\label{fig:FCOSE}
\end{figure}
While LW is better than FCOSE for minimizing quadratic risk (sometimes far better, see Appendix Fig. \ref{fig:esterrorApp}), FCOSE seems to typically outperform SCOSE for Kernel PCA reconstruction error, as seen in the experiments of \cite{kernelstein} (LW performs similar to SCOSE). Our experimental results suggest that FCOSE is also better for HSIC power (it is never worse than $\HSIC$, and rarely worse than $\HSIC^{LW}$). 
Since HSIC is the sum of squared singular values of $\SXY$, it might seem that just shrinking the singular values ``optimally'' should yield the most improvement in power. However, we conjecture that estimating the singular functions accurately is also important. 
Unlike FCOSE, LW and SCOSE do \textit{not} change the eigen (singular) functions of the sample (cross) covariance operator, but only shrink the eigen (singular) values. Fig. \ref{fig:FCOSE} has preliminary experiments that suggests FCOSE could possibly result in more stable estimation of singular functions. 






\section{Conclusion}\label{sec:future}

In this paper we introduced and conducted a study of a novel third family of shrinkage estimators. We demonstrated a very important application of shrinkage, the ability of shrunk test statistics to yield higher power in independence testing. We clearly showed why this initially surprising effect occurs - the permuted test statistics shrunk more than the actual sample statistic. We also gave insights into why FCOSE seems to consistently outperform LW, SCOSE on Kernel PCA and HSIC, even though it has demonstrably worse quadratic risk. We end with some exciting future directions mentioned below.

\textbf{Non-linear shrinkage and kernel spectra.} It seems crucial to figure out the formal connection between the decay of eigenvalues of the kernel, and which kind (linear/nonlinear) of shrinkage works better. Even though shrinkage often helps and rarely hurts, there are indications that non-linear shrinkage methods work better. Alternately, one could find better shrinkage targets, like one with an exponential decay of singular values. 

\textbf{New Estimators}
It is of interest to investigate the properties of a Flexible Sparse COSE (FS-COSE) that uses $\|\beta\|_1$ instead of $\|\beta\|_2^2$.
Since the covariance operator is also trace class, another promising direction could be a Spectrally Sparse COSE (SS-COSE) that we define as:
$$
\min_{S \in \mathcal{HS}(\F,\F)} \frac1{2} \sum_{i=1}^n \left \|\tphi_{x_i}\otimes \tphi_{x_i} - S \right\|_{HS}^2 + \lambda \|S\|_*.
$$

\textbf{Shrinking gram matrices $K$.}
The LW technique can also be extended to shrinking kernel matrices to get more stable estimates. To see this, note that an estimator of the expected kernel matrix can be achieved by averaging the diagonals and off-diagonals separately. Then, using a similar LW strategy, one can shrink towards several targets, including the zero and the identity matrix, and the expected kernel matrix. 

\textbf{Large-scale Experiments.} COCO, NOCCO and HSNIC are  important alternatives to HSIC from \cite{coco}, \cite{nocco}. Our preliminary experiments suggest that shrinkage often improves and never worsens the power of  these related methods. While we see improvements to different degrees on various metrics, it is still unclear to us the extent to which shrinkage can improve performance of common kernelized ML tasks. A large scale empirical comparision with a variety of kernels, distributions and tasks is called for.


\subsection*{Acknowledgements}
We thank the authors of \cite{kernelstein} for sharing their code, Arthur Gretton for a great workshop at Columbia University, Jessica Chemali for tireless help and Rebecca Steorts for useful comments on a first draft.



\bibliographystyle{unsrt}
\bibliography{story} 

\newpage
\appendix


\section{Two Motivating Examples}

When $p$ is comparable or larger than $n$, the sample covariance matrix is quite a poor estimator of the true covariance matrix, even though it is the maximum likelihood estimator and is asymptotically (as $n \rightarrow \infty$ with $p$ fixed) optimal. Since covariance matrices are a special case of kernel covariance operators when a linear kernel is used, that will be our first example, shown in the left plot of Figure \ref{fig:motivateApp}.

We drew a varying number of points from a 50-dimensional standard Gaussian distribution ($n=20,50,100,5000$, $p=50$). For each sample size, we formed the sample covariance matrix, and we calculated the eigenvalues and plot them below (sorted)  as a curve for each $n$.  The main message is that at small sample sizes, the sample covariance matrices are poor estimators of the true covariance - the largest eigenvalues are biased upwards, and the smallest are biased downwards, and only at large values of $n$ (like $5000$) does the spectrum seem accurate.


For our second example, we drew a varying number of points from a joint distribution which is proportional to a 2-dimensional sinusoid of some arbitrarily chosen frequency ($n=20,50,100,5000$, $p=2$). For each dimension $x$ and $y$, we use the Gaussian kernel with the bandwidth parameter chosen by the well-known median heuristic (this detail is of little relevance to the main point of the example, and may be skipped by the unfamiliar reader). For each sample size, we formed the sample cross-covariance operator, calculated the singular values (this operator has at most $n$ non-zero singular values for each $n$), and plotted the logarithm of the top ten singular values as a curve for each $n$. In this example, all the leading singular values are heavily overestimated of factors of over 100, and the extremely small ones (beyond the top 10, not plotted) are estimated as zero.\\

\begin{figure} [h]
\centering
\includegraphics[width=0.6\linewidth]{motivation.pdf}
\caption{Left panel: $x$ axis is the eigenvalue count from 1 to 100 and $y$ axis is its magnitude. Right panel: $x$ axis is the singular value count from 1 to 10 and $y$ axis is logarithm of its magnitude. }
\label{fig:motivateApp}
\end{figure}

\section{The cross-covariance operator as a kernel mean}

In \cite{kernelstein}, the covariance operator is shown to be a kernel mean on $\F\otimes\G$. Since $\SYX := \frac1{n}\sum_{i=1}^n \tilde\phi_{x_i}\otimes \tilde\varphi_{y_i}$, where $\tilde\phi_{x_i}$ and $\tilde\varphi_{y_i}$ are centered versions of $\phi_{x_i}$ and $\varphi_{y_i}$, we can write $\SYX$ as the kernel mean of the kernel defined by the dot product:

$$
\langle \tilde\phi_{x_i} \otimes \tilde\varphi_{y_i} , \tilde\phi_{x_j}\otimes \tilde\varphi_{y_j} \rangle_{\F\otimes\G}  = \langle \tilde\phi_{x_i},  \tilde\phi_{x_j} \rangle_{\F} \langle  \tilde\varphi_{y_i} , \tilde\varphi_{y_j} \rangle_{\G}  = \tilde k(x_i,x_j)\tilde \ell(x_i,x_j).
$$

\newpage

\section{Quadratic Risk, more kernels}
\begin{figure} [h]
\centering
\includegraphics[width=0.48\linewidth]{extin_linear.pdf}
\includegraphics[width=0.48\linewidth]{extin_poly.pdf}
\caption{\small Both panels show the quadratic risk $\E\|S-\SigXY\|_{HS}^2$ for $S\in\{\SXY,\SXY^{LW},\SXY^{S},\SXY^{F}\}$ for linear (left) and polynomial (right) kernels. The same dataset as shown in Fig.\ref{fig:LWeg} was used. $\SigXY$ was approximated from 5000 points, and $\E$ was approximated by 1000 experiment repetitions. Takeaways: The improvement for the linear kernel is slightly more than for the polynomial, the relative gain goes down with $n$ and LW/SCOSE outperform FCOSE even more than in the Gaussian and Laplace case (Fig.\ref{fig:esterror}).}
\label{fig:esterrorApp}
\end{figure}


\section{Kernel Independence Testing, more examples}

\begin{figure} [h]
\centering
\includegraphics[width=0.28\linewidth]{ratios_extin_radius_6.pdf}
\hspace{0.2in}
\includegraphics[width=0.28\linewidth]{ratios_freq_2.pdf}
\hspace{0.2in}
\includegraphics[width=0.28\linewidth]{ratios_theta_angle_4.pdf}
\caption{\small  Illustration of how shrinkage improves the power of independence testing. Shrinking increases power because, under $H_1$, HSIC shrinks less than under $H_0$. This is demonstrated in three example experiments in the bar plots, which shows the ratio between the mean real HSIC and the mean 95 percentile HSIC from the permutation test that is used as a threshold for rejecting $H_0$. 
  For the shrunk estimators the unpermuted HSIC is differentiated better from the permuted HSICs than for the unshrunk one. [Left] This from the distribution in the experiment from Fig.4 (top), a hollow circle with radius 2.2. [Middle] This is from the distribution in Fig.9 (top), a 2D grid with frequency 1/3 and $n=15$. [Right] This is from the distribution in Fig.9 (bottom) with four gaussians tilted by $\theta=3/32\pi$. In Fig.3, the distribution shown was the one in Fig.4 (bottom) with four Gaussians tilted by $\theta = 1/16\pi$ and $n=30$. }
\label{fig:hsicscatterApp}
\end{figure}

\newpage



\section{HSIC Power, more examples}

\begin{figure} [h!]
\centering
\includegraphics[width=0.95\linewidth]{HSIC2diff_N_freq.pdf}
\includegraphics[width=0.95\linewidth]{HSIC2diff_theta.pdf}
\caption{\small We present results from 100 repetitions of two independent examples (top, bottom). The leftmost panel shows the scatter plot of the distribution being used. In the first example, we fix the spatial frequency (1/3) of a 2D grid and we increase the number of point. In the second, we keep the number of points fixed and tilt the four clusters with an increasing angle. The middle panels show the power of shrunk and unshrunk HSICs versus varying radius (top) or sample size (bottom). In each example, the top right panel shows how often shrunk HSIC gets it right when HSIC gets it right (correctly rejects the null) - higher is better and any drop from 1 shows a worsening of shrunk estimators. The bottom right panel shows how often shrunk HSIC gets it right when HSIC makes a mistake (incorrectly fails to reject the null) - higher is better and any rise from 0 shows an improvement for shrunk estimators. Gray rectangles indicates regions with too few samples to estimate the probability. From the top right panels, one may conclude that the shrunk HSICs are almost never worse, and from the bottom right panels, we see that they can be significantly better. The middle panels show that you can get improvements in the low power as well as high power regime. FCOSE and LW both outperform HSIC in the first case, while in the second only FCOSE with LW showing nearly no improvement. }
\label{fig:hsicpowerApp}
\end{figure}

\newpage

\section{Understanding FCOSE, more examples}

\begin{figure} [h]
\centering
\includegraphics[width=0.7\linewidth]{F1_G1_example2.pdf}
\caption{\small The left panels shows 30, 1000 samples from the same distribution. In the middle panels, the ``true'' left and right singular functions $f_1^*,g_1^*$ of $\SigXY$ (Gaussian kernel) are estimated from 1000 points, and plotted in black. Then, from one example of 30 points, we plot the sample left and right singular functions, before ($f_1,g_1$, blue) and after ($f_1^F, g_1^F$, purple) shrinking. The right panels compare $f_1,f_1^F$ and $g_1,g_1^F$ for 100 different examples (100 different 30 point samples). 
The right panels show that over the 100 repeated experiments, the FCOSE error for the singular functions is on average smaller than the unshrunk error, demonstrating its better accuracy in estimating singular vectors. We verified this by calculating the MSE between the unshrunk/shrunk functions and the ``truth'', i.e. $\E\|f_1 - f^*_1\|_K^2$ and $\E\|f_1^F - f^*_1\|_K^2$. MSE of unshrunk estimators of $f$ is 0.5488 and of $g$ is 0.7286, which were larger than MSE of FCOSE for $f$ which was 0.5207, and for $g$ which was 0.7336. 
}
\end{figure}
\begin{figure} [h!]
\centering
\includegraphics[width=0.7\linewidth]{F1_G1_example3.pdf}
\caption{\small The left panels shows 30, 1000 samples from the same distribution. In the middle panels, the ``true'' left and right singular functions $f_1^*,g_1^*$ of $\SigXY$ (Gaussian kernel) are estimated from 1000 points, and plotted in black. Then, from one example of 30 points, we plot the sample left and right singular functions, before ($f_1,g_1$, blue) and after ($f_1^F, g_1^F$, purple) shrinking. The right panels compare $f_1,f_1^F$ and $g_1,g_1^F$ for 100 different examples (100 different 30 point samples). 
The right panels show that over the 100 repeated experiments, the FCOSE error for the singular functions is on average smaller than the unshrunk error, demonstrating its better accuracy in estimating singular vectors. We verified this by calculating the MSE between the unshrunk/shrunk functions and the ``truth'', i.e. $\E\|f_1 - f^*_1\|_K^2$ and $\E\|f_1^F - f^*_1\|_K^2$. MSE of unshrunk estimators of $f$ is 0.1875 and of $g$ is 0.2429, which were larger than MSE of FCOSE for $f$ which was 0.0816, and for $g$ which was 0.1028. 
}
\end{figure}

\end{document}

arxiv.bbl
\begin{thebibliography}{10}

\bibitem{kernelpca}
Bernhard Scholkopf, Alexander Smola, and Klaus-Robert M{\"u}ller.
\newblock Kernel principal component analysis.
\newblock In {\em Advances in kernel methods-support vector learning}.
  Citeseer, 1999.

\bibitem{kerneldimred}
Kenji Fukumizu, Francis~R Bach, Michael~I Jordan, et~al.
\newblock Kernel dimension reduction in regression.
\newblock {\em The Annals of Statistics}, 37(4):1871--1905, 2009.

\bibitem{kernelstructdisc}
Matthew~B Blaschko and Arthur Gretton.
\newblock A hilbert-schmidt dependence maximization approach to unsupervised
  structure discovery.
\newblock In {\em International Workshop on Mining and Learning with Graphs},
  2008.

\bibitem{hsic}
Arthur Gretton, Olivier Bousquet, Alex Smola, and Bernhard Sch{\"o}lkopf.
\newblock Measuring statistical dependence with hilbert-schmidt norms.
\newblock In {\em Algorithmic learning theory}, pages 63--77. Springer, 2005.

\bibitem{kernelgranger}
Xiaohai Sun.
\newblock Assessing nonlinear granger causality from multivariate time series.
\newblock In {\em Machine Learning and Knowledge Discovery in Databases}, pages
  440--455. Springer, 2008.

\bibitem{graphicalmodels}
Francis~R Bach and Michael~I Jordan.
\newblock Learning graphical models with mercer kernels.
\newblock In {\em NIPS}, volume~15, pages 1009--1016, 2002.

\bibitem{kernelfeature}
Le~Song, Alex Smola, Arthur Gretton, Justin Bedo, and Karsten Borgwardt.
\newblock Feature selection via dependence maximization.
\newblock {\em The Journal of Machine Learning Research}, 98888(1):1393--1434,
  2012.

\bibitem{kerneltargetshift}
Kun Zhang, Krikamol Muandet, Zhikun Wang, et~al.
\newblock Domain adaptation under target and conditional shift.
\newblock In {\em Proceedings of the 30th International Conference on Machine
  Learning (ICML-13)}, pages 819--827, 2013.

\bibitem{kernelbayesrule}
Kenji Fukumizu, Le~Song, and Arthur Gretton.
\newblock Kernel bayes' rule.
\newblock In {\em NIPS}, pages 1737--1745, 2011.

\bibitem{ledoitwolf}
Olivier Ledoit and Michael Wolf.
\newblock A well-conditioned estimator for large-dimensional covariance
  matrices.
\newblock {\em Journal of multivariate analysis}, 88(2):365--411, 2004.

\bibitem{covmat1}
Michael~J Daniels and Robert~E Kass.
\newblock Shrinkage estimators for covariance matrices.
\newblock {\em Biometrics}, 57(4):1173--1184, 2001.

\bibitem{covmat2}
Juliane Sch{\"a}fer, Korbinian Strimmer, et~al.
\newblock A shrinkage approach to large-scale covariance matrix estimation and
  implications for functional genomics.
\newblock {\em Statistical applications in genetics and molecular biology},
  4(1):32, 2005.

\bibitem{covmat3}
Jianhua~Z Huang, Naiping Liu, Mohsen Pourahmadi, and Linxu Liu.
\newblock Covariance matrix selection and estimation via penalised normal
  likelihood.
\newblock {\em Biometrika}, 93(1):85--98, 2006.

\bibitem{stein}
Charles Stein et~al.
\newblock Inadmissibility of the usual estimator for the mean of a multivariate
  normal distribution.
\newblock In {\em Proceedings of the Third Berkeley symposium on mathematical
  statistics and probability}, volume~1, pages 197--206, 1956.

\bibitem{kernelstein}
Krikamol Muandet, Kenji Fukumizu, Bharath Sriperumbudur, Arthur Gretton, and
  Bernhard Schoelkopf.
\newblock Kernel mean estimation and stein effect.
\newblock In {\em Proceedings of The 31st International Conference on Machine
  Learning}, pages 10--18, 2014.

\bibitem{coco}
Arthur Gretton, Ralf Herbrich, Alexander Smola, Olivier Bousquet, and Bernhard
  Sch{\"o}lkopf.
\newblock Kernel methods for measuring independence.
\newblock {\em The Journal of Machine Learning Research}, 6:2075--2129, 2005.

\bibitem{nocco}
Kenji Fukumizu, Arthur Gretton, Xiaohai Sun, and Bernhard Sch{\"o}lkopf.
\newblock Kernel measures of conditional dependence.
\newblock In {\em NIPS}, volume~20, pages 489--496, 2007.

\end{thebibliography}
