\documentclass[twoside,11pt]{article}

% Any additional packages needed should be included after jmlr2e.
% Note that jmlr2e.sty includes epsfig, amssymb, natbib and graphicx,
% and defines many common macros, such as 'proof' and 'example'.
%
% It also sets the bibliographystyle to plainnat; for more information on
% natbib citation styles, see the natbib documentation, a copy of which
% is archived at http://www.jmlr.org/format/natbib.pdf

\usepackage{jmlr2e}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{latexsym}
\usepackage{epsfig}
\usepackage{color}
\usepackage{graphicx}
\usepackage{float}
\usepackage{natbib}
%\usepackage{amsthm}
%\usepackage{chicago}

% ----------------------------------------------------------------

% THEOREMS -------------------------------------------------------
\newtheorem{thm}{Theorem}
\newtheorem{cor}[thm]{Corollary}
\newtheorem{lma}[thm]{Lemma}
\newtheorem{prop}[thm]{Proposition}

% MATH -----------------------------------------------------------
\newcommand{\norm}[1]{\left\Vert#1\right\Vert}
\newcommand{\abs}[1]{\left\vert#1\right\vert}
\newcommand{\set}[1]{\left\{#1\right\}}
\newcommand{\Real}{\mathbb R}
\newcommand{\eps}{\varepsilon}
\newcommand{\To}{\longrightarrow}
\newcommand{\BX}{\mathbf{B}(X)}
\newcommand{\A}{\mathcal{A}}
\newcommand{\celsius}{${}^{\circ}${\kern-0.5mm}C}
\newcommand{\farenheit}{${}^{\circ}${\kern-0.5mm}F}
% ----------------------------------------------------------------
\floatstyle{ruled}
\newfloat{algorithm}{hhh}{loa}
\floatname{algorithm}{Algorithm}

\DeclareMathAlphabet{\mathitbf}{OT1}{cmr}{bx}{it}
\newcommand{\vx}{\mathitbf{x}}
\newcommand{\vz}{\mathitbf{z}}
\newcommand{\vv}{\mathitbf{v}}
\newcommand{\va}{\mathitbf{a}}
\newcommand{\vb}{\mathitbf{b}}
\newcommand{\vw}{\mathitbf{w}}
\newcommand{\vu}{\mathitbf{u}}
\newcommand{\C}{{\mathbb{C}}}
\newcommand{\R}{{\mathbb{R}}}
\newcommand{\N}{{\mathbb{N}}}
\newcommand{\cF}{{\cal F}}
\newcommand{\cG}{{\cal G}}
\newcommand{\cH}{{\cal H}}
\newcommand{\indep}{{\bot\negthickspace\negthickspace\bot}}
\newcommand{\hK}{\hat{K}}
\newcommand{\hS}{\hat{\Sigma}}
\newcommand{\la}{\langle}
\newcommand{\ra}{\rangle}

\newcommand{\vbeta}{\boldsymbol{\beta}}
\newcommand{\eq}[1]{Equation~(\ref{#1})}
\newcommand{\mysec}[1]{Section~\ref{sec:#1}}

\jmlrheading{5}{2004}{73-99}{5/03; Revised 10/03}{1/04}{Kenji
Fukumizu, Francis R. Bach, and Michael I. Jordan}
\ShortHeadings{Kernel Dimensionality Reduction}{Fukumizu, Bach and Jordan}
\firstpageno{73}

\begin{document}

\title{Dimensionality Reduction for Supervised Learning
with Reproducing Kernel Hilbert Spaces}

\author{\name Kenji Fukumizu
\email fukumizu@ism.ac.jp \\
\addr Institute of Statistical Mathematics \\
4-6-7 Minami-Azabu, Minato-ku, Tokyo 106-8569, Japan \AND \name
Francis R. Bach
\email fbach@cs.berkeley.edu \\
\addr
Computer Science Division\\
University of California\\
Berkeley, CA 94720, USA
\AND
\name Michael I. Jordan
\email jordan@cs.berkeley.edu \\
\addr
Computer Science Division and Department of Statistics\\
University of California\\
Berkeley, CA 94720, USA
}

\editor{Chris Williams}


\maketitle



\begin{abstract}
We propose a novel method of dimensionality reduction for
supervised learning problems.  Given a regression or
classification problem in which we wish to predict a response
variable $Y$ from an explanatory variable $X$, we treat the
problem of dimensionality reduction as that of finding a
low-dimensional ``effective subspace'' for $X$ which retains the
statistical relationship between $X$ and $Y$. We show that this
problem can be formulated in terms of conditional independence. To
turn this formulation into an optimization problem we establish a
general nonparametric characterization of conditional independence
using covariance operators on reproducing kernel Hilbert spaces.
This characterization allows us to derive a contrast function for
estimation of the effective subspace. Unlike many conventional
methods for dimensionality reduction in supervised learning, the
proposed method requires neither assumptions on the marginal
distribution of $X$, nor a parametric model of the conditional
distribution of $Y$.  We present experiments that compare the
performance of the method with conventional methods.

\end{abstract}


\begin{keywords}
regression, dimensionality reduction, variable selection, feature
selection, kernel methods, conditional independence.
\end{keywords}


\section{Introduction}

Many statistical learning problems involve some form of dimensionality
reduction, either explicitly or implicitly.  The goal may be one of
\emph{feature selection}, in which we aim to find linear or nonlinear
combinations of the original set of variables, or one of \emph{variable
selection}, in which we wish to select a subset of variables from the
original set.  The setting may be unsupervised learning, in which a
set of observations of a random vector $X$ are available, or supervised
learning, in which desired responses or labels $Y$ are also available.
Developing methods for dimensionality reduction requires being clear on
the goal and the setting, as methods developed for one combination of
goal and setting are not generally appropriate for another.  There are
additional motivations for dimensionality reduction that it is also
helpful to specify, including: providing a simplified explanation of a
phenomenon for a human (possibly as part of a visualization algorithm),
suppressing noise so as to make a better prediction or decision, or
reducing the computational burden.  These various motivations are
often complementary.

In this paper we study dimensionality reduction in the setting of
supervised learning.  Thus, we consider problems in which our data
consist of observations of $(X,Y)$ pairs, where $X$ is an
$m$-dimensional explanatory variable and where $Y$ is an
$\ell$-dimensional response. The variable $Y$ may be either
continuous or discrete.  We refer to these problems generically as
``regression'' problems, which indicates our focus on the
conditional probability density function $p_{Y|X}(y \mid x)$.  In
particular, our framework includes discriminative approaches to
classification problems, where $Y$ is a discrete label.

We wish to solve a problem of feature selection in which the
features are linear combinations of the components of $X$.  In
particular, our modeling framework posits that there is an
$r$-dimensional subspace $S \subset \R^m$ such that
\begin{equation*}
 p_{Y|X}(y \mid x) = p_{Y|\Pi_S X}(y \mid \Pi_S x),
\end{equation*}
for all $x$ and $y$, where $\Pi_S$ is the orthogonal projection of
$\R^m$ onto $S$.  The subspace $S$ is called the \emph{effective
subspace for regression}. Based on a set of observations of
$(X,Y)$ pairs, we wish to recover a matrix whose columns span the
effective subspace.  The effective subspace can help to provide
explanation of the statistical relation between $X$ and $Y$, by
isolating the feature vectors that capture that relation.  Also, 
finding such a space can suppress noise, to the extent that the 
orthogonal direction to the effective subspace is noisy vis-a-vis
the prediction of $Y$.

We approach the problem as a \emph{semiparametric} statistical problem;
in particular, we make no assumptions regarding the conditional distribution
$p_{Y|\Pi_S X}(y \mid \Pi_S x)$, nor do we make any assumptions regarding
the marginal distribution $p_X(x)$.  That is, we wish to estimate a
finite-dimensional parameter (a matrix whose columns span the effective
subspace), while treating the distributions $p_{Y \mid \Pi_S X}(y \mid \Pi_S x)$
and $p_X(x)$ nonparametrically.

Having found an effective subspace, we may then proceed to build a
parametric or nonparametric regression model on that subspace.  Thus
our approach is an explicit dimensionality reduction method for supervised
learning that does not require any particular form of regression model,
and can be used as a preprocessor for any supervised learner.  This can
be compared to the use of methods such as principal components analysis
(PCA) in regression, which also make no assumption regarding the subsequent
regression model, but fail to make use of the response variable $Y$.

There are a variety of related approaches in the literature, but most
of them involve making specific assumptions regarding the conditional distribution
$p_{Y|\Pi_S X}(y \mid \Pi_S x)$, the marginal distribution $p_X(x)$, or both.
For example, classical two-layer neural networks involve a linear transformation
in the first ``layer,'' followed by a specific nonlinear function and
a second layer~\citep{Bishop}.  Thus, neural networks can be seen as
attempting to estimate an effective subspace based on specific assumptions
about the regressor $p_{Y|\Pi_S X}(y \mid \Pi_S x)$.  Similar comments
apply to projection pursuit regression~\citep{Friedman81PP},
ACE~\citep{Breiman85ACE} and additive models~\citep{Hastie86},
all of which provide a methodology for dimensionality reduction in which
an additive model $E[Y \mid X]=g_1(\beta_1^TX)+ \cdots + g_K(\beta_K^TX)$
is assumed for the regressor.

Bayesian approaches to variable selection include the Automatic
Relevance Determination (ARD) method proposed by~\cite{NealBayesian96}
for neural networks.  \cite{Vivarelli99} have adapted this approach
to feature selection in the setting of Gaussian process regression.
These methods also depend on a specific parametric model, on which
a Bayesian prior distribution is assumed.

Canonical correlation analysis (CCA) and partial least
squares~\citep[PLS,][]{Holkuldsson88,Helland88} are classical
multivariate statistical methods that can be used for
dimensionality reduction in
regression~\citep{Fung_etal02,NguyenRocke02}.  These methods are
based on a linearity assumption for the regressor, however, and
thus are quite strongly parametric.

The line of research that is closest to our work has its origin
in a technique known as sliced inverse regression
\citep[SIR,][]{KerChauLi91}.  SIR is a semiparametric method for
finding effective subspaces in regression. The basic idea is that
the range of the response variable $Y$ is partitioned into a set
of ``slices,'' and the sample means of the observations $X$ are
computed within each slice.  This can be viewed as a rough
approximation to the inverse regression of $X$ on $Y$.  For
univariate $Y$ the method is particularly easy to implement.
Noting that the inverse regression must lie in the effective
subspace if the forward regression lies in such a subspace,
principal component analysis is then used on the sample means to
find the effective subspace.  \cite{KerChauLi91} has shown that
this approach can find effective subspaces, but only under strong
assumptions on the marginal distribution $p_X(x)$---in particular,
the marginal distribution must be elliptically symmetric.

Further developments in the wake of SIR include principal Hessian
directions~\citep[pHd,][]{KerChauLi92}, and sliced average variance
estimation~\citep[SAVE,][]{CookWeisberg91, CookYin01}.
These are all semiparametric methods in that they make no assumptions about the
regressor~\citep[see also][]{Cook}. However, they again place strong restrictions on
the probability distribution of the explanatory variables. If these assumptions
do not hold, there is no guarantee of finding the effective subspace.

There are also related nonparametric approaches that estimate the derivative of
the regressor to achieve dimensionality reduction, based on the fact that the
derivative of the conditional expectation $E[y \mid B^T x]$ with respect to $x$
belongs to the effective subspace~\citep{Samarov93,Hristache01}.
However, nonparametric estimation of derivatives is quite challenging
in high-dimensional spaces.

There are also dimensionality reduction methods with a semiparametric flavor
in the area of classification, notably the work of \cite{torkkola}, who has
proposed using nonparametric estimation of the mutual information between
$X$ and $Y$, and subsequent maximization of this estimate of mutual information
with respect to a matrix representing the effective subspace.

In this paper we present a novel semiparametric method for dimensionality
reduction that we refer to as \emph{Kernel Dimensionality Reduction (KDR)}.
KDR is based on the estimation and optimization of a particular class of
operators on reproducing kernel Hilbert spaces~\citep{Aronszajn50}.
Although our use of reproducing kernel Hilbert spaces is related to their
role in algorithms such as the support vector machine and kernel
PCA~\citep{SVM92,SVregression,KernelPCA}, where the kernel function
allows linear operations in function spaces to be performed in a
computationally-efficient manner, our work differs in that it cannot
be viewed as a ``kernelization'' of an underlying linear algorithm.
Rather, we use reproducing kernel Hilbert spaces to provide characterizations
of general notions of independence, and we use these characterizations to
design objective functions to be optimized.  We build on earlier work by
\cite{KernelICA}, who showed how to use reproducing kernel Hilbert spaces to
characterize \emph{marginal independence} between pairs of variables,
and thereby design an objective function for independent component analysis.
In the current paper, we extend this line of work, showing how to characterize
\emph{conditional independence} using reproducing kernel Hilbert spaces.
We achieve this by expressing conditional independence in terms of covariance
operators on reproducing kernel Hilbert spaces.

How does conditional independence relate to our dimensionality
reduction problem?  Recall that our problem is to find a projection
$\Pi_S$ of $X$ onto a subspace $S$ such that the conditional
probability of $Y$ given $X$ is equal to the conditional probability
of $Y$ given $\Pi_S X$. This is equivalent to finding a projection
$\Pi_S$ which makes $Y$ and $(I-\Pi_S)X$ conditionally independent
given $\Pi_S X$.  Thus we can turn the dimensionality reduction
problem into an optimization problem by expressing it in terms of
covariance operators.

In a presence of a finite sample, we need to estimate the covariance
operator so as to obtain a sample-based objective function that we
can optimize.  We derive a natural plug-in estimate of the covariance
operator, and find that the resulting estimate is identical to the
\emph{kernel generalized variance} that has been described earlier by
\cite{KernelICA} in the setting of independent component analysis.
In that setting, the goal is to measure departures from independence,
and the minimization of the kernel generalized variance can be viewed as
a surrogate for minimizing a certain mutual information.  In the
dimensionality reduction setting, on the other hand, the goal is to
measure \emph{conditional} independence, and minimizing the kernel
generalized variance can be viewed as a surrogate for \emph{maximizing} a
certain mutual information.  Not surprisingly, the derivation that leads
to the kernel generalized variance that we present here is quite different
from the one presented in the earlier work on kernel ICA.  Moreover, the
argument that we present here can be viewed as providing a rigorous foundation
for other, more heuristic, ways in which the kernel generalized variance has been
used, including the model selection algorithms for graphical models
presented by \cite{BachJordanNIPS02}.

The paper is organized as follows.  In \mysec{setup}, we introduce
the problem of dimensionality reduction for supervised learning,
and describe its relation with conditional independence and mutual
information. \mysec{kdr} derives the objective function for
estimation of the effective subspace for regression, and describes
the KDR method. All of the mathematical details needed for the
results in \mysec{kdr} are presented in the Appendix, which also
provides a general introduction to covariance operators in
reproducing kernel Hilbert spaces. In \mysec{experiments}, we
present a series of experiments that test the effectiveness of our
method, comparing it with several conventional methods.
\mysec{varsel} describes an extension of KDR to the problem of
variable selection.  \mysec{conclusions} presents our conclusions.


\section{Dimensionality Reduction for Regression}
\label{sec:setup}

We consider a regression problem, in which $Y$ is an
$\ell$-dimensional random vector, and $X$ is an $m$-dimensional
explanatory variable. (Note again that we use ``regression'' in a
generic sense that includes both continuous and discrete $Y$).
The probability density function of $Y$ given $X$ is denoted
by $p_{Y|X}(y \mid x)$. Assume that there is an $r$-dimensional
subspace $S \subset \R^m$ such that
\begin{equation}
\label{eq:cond_def}
 p_{Y|X}(y \mid x) = p_{Y|\Pi_S X}(y \mid \Pi_S x),
\end{equation}
for all $x$ and $y$, where $\Pi_S$ is the orthogonal projection of
$\R^m$ onto $S$.  The subspace $S$ is called the {\it effective
subspace for regression}.

The problem that we treat here is that of finding the subspace $S$
given an \emph{i.i.d.}~sample $\{ (X_1,Y_1), \ldots, (X_n,Y_n) \}$
from $p_X$ and $p_{Y|X}$.  The crux of the problem is that we assume
no \emph{a priori} knowledge of the regressor, and place no assumptions
on the conditional probability $p_{Y|X}$.

As in the simpler setting of principal component analysis, we make
the (generally unrealistic) assumption that the dimensionality $r$
is known and fixed.  While choosing the dimensionality is a very
important problem, we separate it from the task of finding the
best subspace of a fixed dimensionality, and mainly focus on the
latter problem.  In \mysec{conclusions}, we briefly discuss
various approaches to the estimation of the dimensionality.

The notion of effective subspace can be formulated in terms of
conditional independence.  Let $(B,C)$ be the $m$-dimensional
orthogonal matrix such that the column vectors of $B$ span the
subspace $S$, and define $U=B^TX$ and $V=C^T X$.  Because $(B, C)$
is an orthogonal matrix, we have
\begin{equation}
p_X(x) = p_{U,V}(u,v), \qquad p_{X,Y}(x,y) = p_{U,V,Y}(u,v,y),
\label{eq:var_trans}
\end{equation}
for the probability density functions.  From
\eq{eq:var_trans}, \eq{eq:cond_def} is equivalent to
\begin{equation*}
p_{Y|U,V}(y \mid u,v) = p_{Y|U}(y \mid u).
\end{equation*}
This shows that the effective subspace $S$ is the one which makes
$Y$ and $V$ conditionally independent given $U$ (see Figure~\ref{fig:cond_indep}).

Mutual information provides another point of view on the equivalence
between conditional independence and the existence of the effective
subspace. From \eq{eq:var_trans}, it is straightforward to
see that
\begin{equation}
\label{eq:mutual}
I(Y,X) = I(Y,U) + E_{U}\bigl[ I(Y|U,  V|U ) \bigr],
\end{equation}
where $I(Z,W)$ denotes the mutual information defined by
\begin{equation*}
I(Z,W) := \int\int p_{Z,W}(z,w)\log \frac{p_{Z,W}(z,w)}{p_Z(z)
p_W(w) }dz\;dw.
\end{equation*}
Because \eq{eq:cond_def} means $I(Y,X) = I(Y,U)$, the
effective subspace $S$ is characterized as the subspace which retains the
mutual information of $X$ and $Y$ by the projection onto that subspace,
or equivalently, which gives $I(Y|U,V|U) = 0$. This is again the
conditional independence of $Y$ and $V$ given $U$.

The expression in \eq{eq:mutual} can be understood in terms of the
decomposition of the mutual information according to a tree-structured
graphical model---a quantity that has been termed the \emph{T-mutual information}
by \cite{TCAjmlr}.  Considering the tree $Y-U-V$ in Figure~\ref{fig:cond_indep}(b),
we have that the T-mutual information $I^T$ is given by
\begin{equation*}
I^T = I(Y,U,V) - I(Y,U) - I(U,V).
\end{equation*}
This is equal to the KL-divergence between a probability distribution
on $(Y,U,V)$ and its projection onto the family of distributions that
factor according to the tree; that is, the set of distributions that
verify $Y \indep V \mid U$.  Using \eq{eq:var_trans}, we can easily see
that $I(Y,U,V)=I(Y,X)+I(U,V)$, and thus we obtain
\begin{equation}
\label{eq:TMI2}
I^T = I(Y,X) - I(Y,U) = E_U[ I(Y|U, V|U)].
\end{equation}
Then, dimensionality reduction for regression can be viewed as the problem
of minimizing the T-mutual information for the fixed tree structure
in Figure~\ref{fig:cond_indep}(b).


\begin{figure}
  \begin{center}
  \includegraphics[height=4cm]{ps/cond_indep.eps}\\
  (a) \hspace{4.8cm} (b) \hspace{1cm}
  \caption{Graphical representation of dimensionality reduction for regression.
  The variables $Y$ and $V$ are conditionally independent given $U$, where $X=(U,V)$.}
  \label{fig:cond_indep}
  \end{center}
\end{figure}

From \eq{eq:mutual} and \eq{eq:TMI2}, we see that there are two
approaches to solve the problem of dimensionality reduction for
regression; one is to maximize the mutual information $I(Y,U)$,
and the other is to make $Y$ and $V$ conditionally independent
given $U$. We choose the latter in deriving our method.  Direct
evaluation of mutual information is not a straightforward task in
general, because it requires an explicit form for the probability
density functions. Although assuming a specific probability model
leads readily to a solution, it implies a restriction of the
range of problems to which the method can be applied, and does
not satisfy our goal of developing a general semiparametric method.
Nonparametric estimation of the probabilities also provides an
approach to evaluation of mutual information. However, nonparametric
estimation and numerical integration do not give accurate estimates
for high-dimensional variables, and it is a challenge to make a
nonparametric method viable~\citep{torkkola}.  An alternative
semiparametric approach was presented by \cite{KernelICA}, who
showed that the kernel generalized variance could serve as a
surrogate for the mutual information in the setting of independent
component analysis.  Our approach to dimensionality reduction also
makes use of the kernel generalized variance.  However, we do not
take the mutual information as our point of departure, because we
expect to be far from the setting of mutual independence in which
\cite{KernelICA} showed that the mutual information and kernel
generalized variance are closely related.  Instead, we take an entirely
different path, presenting a rigorous characterization of conditional
independence using reproducing kernels, and showing that this
characterization leads once again to the kernel generalized variance.

\section{Kernel Method for Dimensionality Reduction in Regression}
\label{sec:kdr}

In this section we present our kernel-based method for dimensionality
reduction.  We discuss the basic definition and properties of
cross-covariance operators on reproducing kernel Hilbert spaces,
derive an objective function for characterizing conditional independence
using cross-covariance operators, and finally present a sample-based
objective function based on this characterization.

\subsection{Cross-Covariance Operators on Reproducing Kernel Hilbert
Spaces} \label{sec:cov_op}

We use cross-covariance operators on reproducing kernel Hilbert
spaces to derive an objective function for dimensionality
reduction. While cross-covariance operators are generally defined
for random variables in Banach spaces~\citep{ProbBanach,Baker73},
the theory is much simpler for reproducing kernel Hilbert spaces.
We summarize only basic mathematical facts in this subsection, and
defer the details to the Appendix.  Let $(\cH,k)$ be a reproducing
kernel Hilbert space of functions on a set $\Omega$ with a
positive definite kernel $k: \Omega\times\Omega \to \R$.  The
inner product of $\cH$ is denoted by $\la \cdot, \cdot \ra_\cH$.
We consider only real Hilbert spaces for simplicity.  The most
important aspect of reproducing kernel Hilbert spaces is the
reproducing property:
\begin{equation*}
\la f, k(\cdot,x) \ra_\cH = f(x) \quad\text{ for all } x\in\Omega
\text{ and } f \in \cH.
\end{equation*}
Throughout this paper we use the Gaussian kernel
\begin{equation*}
k(x_1,x_2) = \exp\bigl(-\|x_1 - x_2\|^2 /\sigma^2\bigr),
\end{equation*}
which corresponds to a Hilbert space of smooth functions.

Let $(\cH_1,k_1)$ and $(\cH_2, k_2)$ be reproducing kernel Hilbert
spaces over measurable spaces $(\Omega_1, {\cal B}_1)$ and
$(\Omega_2, {\cal B}_2)$, respectively, with $k_1$ and $k_2$
measurable. For a random vector $(X, Y)$ on $\Omega_1\times
\Omega_2$, the cross-covariance operator from $\cH_1$ to $\cH_2$
is defined by the relation
\begin{equation}
\la g, \Sigma_{YX}f \ra_{\cH_2} = E_{XY}[ f(X)g(Y) ] -
E_X[f(X)]E_Y[g(Y)]
 \label{eq:cov_op}
 \end{equation}
for all $f\in \cH_1$ and $g\in \cH_2$.  \eq{eq:cov_op} implies
that the covariance of $f(X)$ and $g(Y)$ is given by the action of
the linear operator $\Sigma_{YX}$ and the inner product. (See the
Appendix for a basic exposition of cross-covariance operators.)

Covariance operators provide a useful framework for discussing
conditional probability and conditional independence.  As we show
in Corollary \ref{cor:cond} of the Appendix, the following
relation holds between the conditional expectation and the
cross-covariance operator, given that $\Sigma_{XX}$ is
invertible:\footnote{Even if $\Sigma_{XX}$ is not invertible, a
similar fact holds.  See Corollary \ref{cor:cond}.}
\begin{equation}
E_{Y|X}[ g(Y) \mid X ] = \Sigma_{XX}^{-1} \Sigma_{XY} g \quad
\text{ for all } g \in \cH_2,
\label{eq:cond_exp}
\end{equation}
 \eq{eq:cond_exp} can be
understood by analogy to the conditional expectation of
Gaussian random variables. If $X$ and $Y$ are Gaussian random
variables, it is well known that the conditional expectation is
given by
\begin{equation}
\label{eq:Gauss_cond}
    E_{Y|X}[ a^T Y\mid X=x] = x^T \Sigma_{XX}^{-1}\Sigma_{XY}a,
\end{equation}
for an arbitrary vector $a$, where $\Sigma_{XX}$ and $\Sigma_{XY}$
are the variance-covariance matrices in the ordinary sense.


\subsection{Conditional Covariance Operators and Conditional Independence}
\label{sec:cond_cov}

We derive an objective function for characterizing conditional independence
using cross-covariance operators.  Suppose we have random variables $X$
and $Y$ on $\R^m$ and $\R^\ell$, respectively. The variable $X$ is
decomposed into $U\in\R^r$ and $V\in\R^{m-r}$ so that $X=(U,V)$.
For the function spaces corresponding to $Y$, $U$ and $V$, we consider the
reproducing kernel Hilbert spaces $(\cH_1, k_1)$, $(\cH_2, k_2)$,
and $(\cH_3,k_3)$ on $\R^\ell$, $\R^r$, and $\R^{m-r}$,
respectively, each endowed with Gaussian kernels.  We define the
conditional covariance operator $\Sigma_{YY|U}$ on $\cH_1$ by
\begin{equation}
\label{eq:def_cond_var}
\Sigma_{YY|U} := \Sigma_{YY}-
\Sigma_{YU}\Sigma_{UU}^{-1}\Sigma_{UY},
\end{equation}
where $\Sigma_{YY}$, $\Sigma_{UU}$, $\Sigma_{YU}$ are the
corresponding covariance operators.  As shown by Proposition
\ref{prop:cond_indep_pre} in the Appendix, the operator
$\Sigma_{YY|U}$ captures the conditional variance of a random
variable in the following way:
\begin{equation}
\label{eq:cond_var} \la g , \Sigma_{YY|U}g\ra_{\cH_1} = E_{U}\bigl[
\text{Var}_{Y|U}[ g(Y) \mid U] \bigr],
\end{equation}
where $g$ is an arbitrary function in $\cH_1$.  As in the case of
\eq{eq:Gauss_cond}, we can make an analogy to Gaussian variables.
In particular, Equations~(\ref{eq:def_cond_var}) and (\ref{eq:cond_var})
can be viewed as the analogs of the following well-known equality
for the conditional variance of Gaussian variables:
\begin{equation*}
\text{Var}[ a^T Y \mid U ] = a^T (\Sigma_{YY} -
\Sigma_{YU}\Sigma_{UU}^{-1} \Sigma_{UY} )a.
\end{equation*}

It is natural to use minimization of $\Sigma_{YY|U}$ as the basis of
a method for finding the most informative direction $U$.  This
intuition is justified theoretically by Theorem
\ref{thm:cond_indep1} in the Appendix.  That theorem shows
that
\begin{equation}
\label{eq:cond_cov_comp}
\Sigma_{YY|U} \geq \Sigma_{YY|X} \qquad\text{for any }\;\; U,
\end{equation}
and
\begin{equation*}
\Sigma_{YY|U}-\Sigma_{YY|X}=0 \qquad \Longleftrightarrow \qquad
Y\indep V \,|\,U,
\end{equation*}
where, in \eq{eq:cond_cov_comp}, the inequality should be
understood as the partial order of self-adjoint operators. From
these relations, the effective subspace $S$ can be characterized
in terms of the solution to the following minimization problem:
\begin{equation}
\label{eq:crit_RKHS} \min_{S} \Sigma_{YY|U},\quad \text{subject to } U=\Pi_S X.
\end{equation}
In the following section we show how to turn this population-based
criterion into a sample-based criterion that can be optimized in
the presence of a finite sample.

\subsection{Kernel Generalized Variance for Dimensionality Reduction}
\label{sec:method}

To derive a sample-based objective function from
\eq{eq:crit_RKHS}, we have to estimate the conditional covariance
operator from data, and choose a specific way to evaluate the size
of self-adjoint operators.  While there are many possibilities for
approaching these two problems, we make a specific choice, adopting
an approach which has been used successfully for independent component
analysis~\citep{KernelICA}.

We describe a regularized empirical estimate of the cross-covariance
operator.  Define $\tilde{k}_{1}^{(i)}\in \cH_1$ by
$\tilde{k}_1^{(i)}=k_1(\cdot,Y_i) - \frac{1}{n}\sum_{j=1}^n
k_1(\cdot,Y_j)$, and $\tilde{k}_{2}^{(i)}\in\cH_2$ similarly using
$U_i$. By replacing the expectation with the empirical average,
the covariance is estimated by
\[
\bigl\la f, \Sigma_{YU} g \bigr\ra_{\cH_1} \approx
\frac{1}{n}\sum_{i=1}^n \bigl\la \tilde{k}_{1}^{(i)}, f \bigr\ra
\bigl\la \tilde{k}_{2}^{(i)}, g \bigr\ra.
\]
Let $\hK_Y$ be the centralized Gram
matrix~\citep{KernelICA,KernelPCA}, defined by
\begin{equation}
\label{eq:hatK}
 \hK_Y = \bigl( I_n - \tfrac{1}{n}{\mathbf
1}_n{\mathbf 1}_n^T \bigr) G_Y \bigl( I_n - \tfrac{1}{n}{\mathbf
1}_n{\mathbf 1}_n^T \bigr),
\end{equation}
where $(G_Y)_{ij}=k_1(Y_i, Y_j)$ is the Gram matrix and ${\mathbf
1}_n = (1,\ldots,1)^T$ is the vector with all elements equal to
$1$.  The matrix $\hK_U$ is defined similarly using
$\{U_i\}_{i=1}^n$.  Then, it is easy to see that $\bigl\la
\tilde{k}_1^{(i)}, \tilde{k}_1^{(j)}\bigr\ra =
\bigl(\hK_Y\bigr)_{ij}$ and $\la \tilde{k}_2^{(i)},
\tilde{k}_2^{(j)}\bigr\ra = \bigl(\hK_U\bigr)_{ij}$. If we
decompose $f$ and $g$ as $f=\sum_{\ell=1}^n a_i
\tilde{k}_1^{(\ell)} + f^\perp$ and $g=\sum_{m=1}^n b_m
\tilde{k}_2^{(m)} + g^\perp$, where $f^\perp$ and $g^\perp$ are
orthogonal to the linear hull of $\{ \tilde{k}_1^{(i)} \}_{i=1}^n$
and $\{ \tilde{k}_2^{(i)} \}_{i=1}^n$, respectively, we see that
the covariance is approximated by
\[
    \bigl\la f, \Sigma_{YU} g\bigr\ra_{\cH_1} \approx \sum_{\ell=1}^n \sum_{m=1}^n
    a_\ell b_m \bigl(\hat{K}_Y \hat{K}_U)_{lm}.
\]
Thus, by restricting the covariance operator $\Sigma_{YU}$ to the
$n$-dimensional subspaces spanned by $\{ k_1^{(i)} \}_{i=1}^n$ and
$\{ k_2^{(i)} \}_{i=1}^n$, we can estimate the operator by
\[
\hS_{YU} = \hat{K}_Y \hat{K}_U.
\]
In estimating $\Sigma_{YY}$ and $\Sigma_{UU}$, we use the same
regularization technique as \cite{KernelICA}.  The empirical
conditional covariance matrix $\hS_{YY|U}$ is then defined by
\begin{equation}
\label{eq:empcondcov}
 \hS_{YY|U} :=
\hS_{YY}-\hS_{YU}\hS_{UU}^{-1}\hS_{UY} = (\hK_Y + \varepsilon
I_n)^2 - \hK_Y\hK_U (\hK_U + \varepsilon I_n)^{-2}\hK_U \hK_Y,
\end{equation}
where $\varepsilon > 0$ is a regularization constant.

The size of $\hS_{YY|U}$ in the ordered set of positive definite
matrices can be evaluated by its determinant. Although there are other
choices for measuring the size of $\hS_{YY|U}$, such as the
trace and the largest eigenvalue, we focus on the determinant
in this paper. Using the Schur decomposition
$\text{det} (A-B C^{-1} B^T) = \text{det} \Bigl(
\begin{smallmatrix} A & B \\ B^T & C \end{smallmatrix} \Bigr)/ \text{det} C$,
the determinant of $\hS_{YY|U}$ can be written as follows:
\begin{equation*}
\det \hat{\Sigma}_{YY|U}  = \frac{ \det \hS_{[YU][YU]}}{\det
\hS_{UU} },
\end{equation*}
where $\hS_{[YU][YU]}$ is defined by
\begin{equation*}
\hS_{[YU][YU]} = \begin{pmatrix} \hS_{YY} & \hS_{YU} \\ \hS_{UY} &
\hS_{UU} \end{pmatrix} = \begin{pmatrix} (\hK_Y + \varepsilon
I_n)^2 & \hK_Y\hK_U \\ \hK_U \hK_Y & (\hK_U + \varepsilon I_n)^2
\end{pmatrix}.
\end{equation*}
We symmetrize the objective function by dividing by the constant
$\det\hS_{YY}$, which yields the objective function
\begin{equation}
\label{eq:KMI} \frac{ \det \hS_{[YU][YU]} }{ \det \hS_{YY} \det
\hS_{UU} }.
\end{equation}
We refer to the problem of minimizing this function with respect to the
choice of subspace $S$ as \emph{Kernel Dimensionality Reduction (KDR)}.

\eq{eq:KMI} has been termed the ``kernel generalized variance''
by~\cite{KernelICA}, who used it as a contrast function for
independent component analysis.  In that setting, the goal is to
\emph{minimize} a mutual information (among a set of recovered
``source'' variables), in the attempt to obtain independent
components. \cite{KernelICA} showed that the kernel generalized
variance is in fact an approximation of the mutual information of
the recovered sources, when this mutual information is expanded
around the manifold of factorized distributions.  In the current
setting, on the other hand, our goal is to \emph{maximize} the
mutual information $I(Y,U)$, and we certainly do not expect to be
near a manifold in which $Y$ and $U$ are independent.  Thus the
argument for the kernel generalized variance as an objective
function in the ICA setting does not apply here.  What we have
provided in the previous section is an entirely distinct argument
that shows that the kernel generalized variance is in fact an
appropriate objective function for the dimensionality reduction
problem, and that minimizing the kernel generalized variance in
\eq{eq:KMI} can be viewed as a surrogate for maximizing the mutual
information $I(Y,U)$,
while the value of $I(Y,U)$ may not be explicitly related to the
value of the kernel generalized variance.

For the optimization of \eq{eq:KMI}, we use gradient descent with
line search in our experiments.  In a straightforward
implementation, the matrix $B$ is updated iteratively according to
\begin{align}
B(t+1) & = B(t) - \eta \frac{\partial \log \det \hS_{YY|U}}
{\partial B} \nonumber \\
& = B(t) - \eta \text{Tr}\Bigl[ \hS_{YY|U}^{-1} \frac{\partial
\hS_{YY|U}}{\partial B} \Bigr], \label{eq:update}
\end{align}
where $\eta$ is optimized through golden section
search.  The trace term in \eq{eq:update} is rewritten by
\[
    \text{Tr}\Bigl[\hS_{YY|U}^{-1} \frac{\partial
\hS_{YY|U}}{\partial B} \Bigr] = 2 \varepsilon \text{Tr} \Bigl[
\hS_{YY|U}^{-1} \hat{K}_Y \bigl( \hat{K}_U+\varepsilon
I_n\bigr)^{-1} \frac{\partial \hat{K}_U}{\partial B}\bigl(
\hat{K}_U+\varepsilon I_n\bigr)^{-2} \hat{K}_U \hat{K}_Y \Bigr].
\]
All of these matrices are calculated directly from the definitions
in Equations~(\ref{eq:hatK}) and (\ref{eq:empcondcov}).

Given that the numerical task that must be solved in KDR is the
same as the numerical task that must be solved in kernel ICA, we
can import all of the computational techniques developed
by~\cite{KernelICA} for minimizing kernel generalized variance in
the KDR setting.  In particular, we can exploit incomplete
Cholesky decomposition, which approximates an $n\times n$ positive
semidefinite matrix $K$ by $K\approx A A^T$, where $A$ is an
$n\times \ell$ matrix for $\ell< n$.   Application of this
decomposition reduces the $n\times n$ matrices $\hat{K}_Y$ and
$\hat{K}_U$ required in \eq{eq:update} to low-rank approximations.
This diminishes the computational cost associated with multiplying
and inverting large matrices, especially for a large $n$.  For an
exposition of incomplete Cholesky decomposition, see \cite{KernelICA}.
Another computational issue, which may arise in minimizing \eq{eq:KMI},
is the problem of local minima, because the objective function is not
a convex function. To cope with this problem, we make use of an annealing
technique, in which the scale parameter $\sigma$ for the Gaussian
kernel is decreased gradually during the iterations of
optimization.  For a larger $\sigma$, the contrast function is
smoother with fewer local optima, which makes optimization easier.
The search becomes more accurate as $\sigma$ is decreased.


\section{Experimental Results}
\label{sec:experiments}

We study the effectiveness of the new method through experiments,
comparing it with several conventional methods: SIR, pHd, CCA, and
PLS.  For the experiments with SIR and pHd, we use an implementation
in R due to \cite{Weisberg02}.

In all of our experiments, we use a fixed value $0.1$ for the
regularization coefficient $\varepsilon$; empirically, the
performance of the algorithm is robust to small variations in 
this coefficient. This coefficient could also be chosen using
cross-validation.


\subsection{Synthetic Data}
The data sets Data I and Data II comprise one-dimensional $Y$ and
two-dimensional $X=(X_1,X_2)$.  One hundred \emph{i.i.d.}~data
points are generated by
\begin{align*}
\text{I}:  \quad Y & \sim 1/( 1+\exp(-X_1) ) + Z, \\
\text{II}: \quad Y & \sim 2 \exp( - X_1^2 ) + Z,
\end{align*}
where $Z\sim N(0,0.1^2)$, and $X=(X_1,X_2)$ follows a normal
distribution and a normal mixture with two components for Data I 
and Data II, respectively.  The effective subspace is spanned by
$B_0=(1,0)^T$ in both cases.  The data sets are depicted in
Figure~\ref{fig:hut}.


\begin{figure}
\begin{center}
\includegraphics[height=3.5cm]{ps/sigm1_X1Y.eps}\hspace{1cm}
\includegraphics[height=3.5cm]{ps/sigm_XY.eps}\\
\text{I}: $(X_1,Y)$ \hspace{2.7cm} I: $(X_1, X_2)$\\
\vspace{0.5cm}
  \includegraphics[height=3.5cm]{ps/plot_hut_x1y_b.eps}\hspace{1cm}
  \includegraphics[height=3.5cm]{ps/plot_hut_xy_b.eps}\\
  \text{II}: $(X_1,Y)$ \hspace{2.7cm} II: $(X_1, X_2)$
  \caption{Data I and II. One-dimensional $Y$ depends only on $X_1$ in $X=(X_1,X_2)$.}
    \label{fig:hut}
  \end{center}
\end{figure}


Table \ref{tbl:hut} shows the angles between $B_0$ and the
estimated direction.   For Data I, all the methods except PLS
yield a good estimate of $B_0$.  Data II is surprisingly difficult
for the conventional methods, presumably because the distribution
of $X$ is not spherical and the regressor has a strong
nonlinearity.  This indicates the weakness of model-based methods;
if data do not fit the model, the obtained result may not be
meaningful.  The KDR method succeeds in finding the correct
direction for both data sets.


Data III has 300 samples of 17 dimensional $X$ and one dimensional
$Y$, which are generated by
\begin{equation*}
\text{III}: \quad Y \sim 0.9 X_1 + 0.2 \frac{1}{1 + X_{17}} + Z,
\end{equation*}
where $Z \sim N(0,0.01^2)$ and $X$ follows a uniform
distribution on $[0,1]^{17}$.  The effective subspace is given by
$\vb_1=(1,0,\ldots,0)$ and $\vb_2=(0,\ldots,0,1)$.  We compare the
KDR method with SIR and pHd only---CCA and PLS cannot find a
two-dimensional subspace, because $Y$ is one-dimensional.  To
evaluate the accuracy of the results, we use the multiple
correlation coefficient
\begin{equation*}
    R(\vb)= \max_{\vbeta\in B} \frac{\vbeta^T \Sigma_{XX}\vb }
    {\sqrt{\vbeta^T \Sigma_{XX}\vbeta \cdot \vb^T \Sigma_{XX}\vb} }, \qquad
    (\vb \in B_0),
\end{equation*}
which is used in \cite{KerChauLi91}.  As shown in Table \ref{tbl:MG},
the KDR method outperforms the others in finding the weak contribution
of the second direction.


\begin{table}
\begin{center}
\begin{tabular}{c||c|c|c|c|c}
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
  \hline
  & SIR & pHd & CCA & PLS & Kernel \\
  \hline
  I: angle (rad.) & 0.0087 & -0.1971 & 0.0099 & 0.2736 & -0.0014 \\
  \hline
  II: angle (rad.) & -1.5101 & -0.9951 & -0.1818 & 0.4554 & 0.0052 \\
  \hline
\end{tabular}
\caption{Angles between the true and the estimated subspaces for
Data I and II.}\label{tbl:hut}
\end{center}
\end{table}


\begin{table}
\begin{center}
\begin{tabular}{c||c|c|c|c|c|c}
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
  \hline
  & SIR(10) & SIR(15) & SIR(20) & SIR(25) & pHd & Kernel  \\
  \hline
  $R(\vb_1)$ & 0.987 & 0.993 & 0.988 & 0.990 & 0.110 & 0.999\\
  $R(\vb_2)$ & 0.421 & 0.705 & 0.480 & 0.526 & 0.859 & 0.984 \\
  \hline
\end{tabular}
\caption{Correlation coefficients for Data III. SIR($m$) indicates
SIR with $m$ slices.}\label{tbl:MG}
\end{center}
\end{table}




\subsection{Real Data: Classification}
In this section we apply the KDR method to classification problems.
Many conventional methods of dimensionality reduction for regression
are not suitable for classification.  In particular, in the case
of SIR, the dimensionality of the effective subspace must be less
than the number of classes, because SIR uses the average of $X$ in
slices along the variable $Y$.  Thus, in binary classification,
only a one-dimensional subspace can be found, because at most
two slices are available.  The methods CCA and PLS have a similar
limitation on the dimensionality of the effective subspace;
they cannot find a subspace of larger dimensionality than that
of $Y$.  Thus our focus is the comparison between KDR and pHd,
which is applicable to general binary classification problems.
Note that \cite{CookLee99} discuss dimensionality reduction methods for
binary classification, and propose the \emph{difference of covariance
(DOC)} method.  They compare pHd and DOC theoretically, and show that
these methods are the same in binary classification if the population
ratio of the classes is 1/2, which is almost the case in our
experiments.

In the first experiment, we demonstrate the visualization capability of
the dimensionality reduction methods.  We use the {\it Wine} data set
in the UCI machine learning repository~\citep{UCI} to see how the
projection onto a low-dimensional space realizes an effective
description of data.  The wine data consist of 178 samples with
13 variables and a label of three classes.  We apply the KDR
method, CCA, PLS, SIR, and pHd to these data.   Figure~\ref{fig:wine} shows
the projection onto the two-dimensional subspace estimated by each
method.  The KDR method separates the data into three classes most
completely, while CCA also shows perfect separation. We can see
that the data are nonlinearly separable in the two-dimensional
space. The other methods do not separate the classes completely.

Next we investigate how much information on $Y$ is preserved in
the estimated subspace.  After reducing the dimensionality, we use
the support vector machine (SVM) method to build a classifier in
the reduced space, and compare its accuracy with an SVM trained
using the full dimensional vector $X$.\footnote{In our experiments
with the SVM, we used the Matlab Support Vector Toolbox by
S.~Gunn; see http://www.isis.ecs.soton.ac.uk/resources/svminfo.}
Although we should theoretically use the (unknown) optimum
classifier to evaluate the extent of the information preserved in the
subspace, we use an SVM classifier whose parameters were chosen by
exhaustive grid search. We use the {\it Heart-disease} data
set,\footnote{We use the {\it Cleveland} data set, created by Dr.
Robert Detrano of V.A. Medical Center, Long Beach and Cleveland
Clinic Foundation. Although the original data set has five
classes, we use only ``no presence'' (0) and ``presence'' (1-4)
for the binary class labels. Samples with missing values are
removed in our experiments.} {\it Ionosphere}, and {\it
Wisconsin-breast-cancer} from the UCI repository. A description of
these data is presented in Table \ref{tbl:binary_class}.


\begin{table}
  \centering
  \begin{tabular}{c|c|c|c}
    % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
    \hline
    Data set & dim. of $X$  & training sample & test sample \\
    \hline
    Heart-disease & 13 & 149 & 148  \\
    Ionosphere & 34 & 151 & 200 \\
    Breast-cancer-Wisconsin & 30 & 200 & 369 \\ \hline
  \end{tabular}
  \caption{Data description for the binary classification problem.}\label{tbl:binary_class}
\end{table}


Figure~\ref{fig:svm} shows the classification rates for the test
set in subspaces of various dimensionality.   We can see that KDR
yields good separation even in low-dimensional subspaces, while
pHd is much worse in low dimensions.  It is noteworthy that in the
Ionosphere data set the classifier in dimensions 5, 10, and 20
outperforms the classifier in the full dimensional space.  This is
presumably due to the suppression of noise irrelevant to the
prediction of $Y$. 
These results show that
the kernel method successfully finds a subspace which preserves
the class information even when the dimensionality is reduced
significantly.





\begin{figure}
\begin{center}
  \includegraphics[height=5cm]{ps/kernel_b.eps}\hspace{0.5cm}
  \includegraphics[height=5cm]{ps/cca_b.eps}\\
  (a) KDR \hspace{4cm} (b) CCA\\
  \includegraphics[height=5cm]{ps/pls_b.eps}\hspace{0.5cm}
  \includegraphics[height=5cm]{ps/sir_b.eps}\\
  (c) PLS \hspace{4cm} (d) SIR\\
  \includegraphics[height=5cm]{ps/phd_b.eps}\hspace{5.5cm}
  \hbox{}\\
  \hspace{1cm} (e) pHd \hspace*{6.5cm} \\
  \caption{Wine data.  Projections onto the estimated two-dimensional space.
  The symbols `+', `$\Box$', and gray `$\bigcirc$' represent the three
  classes.}
  \label{fig:wine}
  \end{center}
\end{figure}


\begin{figure}
\begin{center}
    (a) Heart-disease \\
    \includegraphics[height=5.5cm]{ps/heart_res_b.eps}\\
  \vspace{0.2cm}
  (b) Ionosphere \\
  \includegraphics[height=5.5cm]{ps/result_ion_gr_b.eps}\\
  \vspace{0.2cm}
  (c) Wisconsin Breast Cancer\\
  \includegraphics[height=5.5cm]{ps/breast_svm_b.eps}\\
  \caption{Classification accuracy of the SVM for test data after dimensionality reduction.}
  \label{fig:svm}
  \end{center}
\end{figure}



\section{Extension to Variable Selection}
\label{sec:varsel}

In this section, we describe an extension of the KDR method to
the problem of variable selection. Variable selection is different
from dimensionality reduction; the former involves selecting a subset of
the explanatory variables $\{ X_1, \ldots, X_m\}$ in order to obtain
a simplified prediction of $Y$ from $X$, while the latter involves
finding linear combinations of the variables.  However, the objective
function that we have presented for dimensionality reduction can be
extended straightforwardly to variable selection.  In particular, given
a fixed number of variables to be selected, we can compare the KGV for
subspaces spanned by combinations of this number of selected variables.
This gives a reasonable way to select variables, because for a
subset $ W=\{X_{j_1},\ldots,X_{j_r}\}\subset \{X_1,\ldots,X_m\}$,
the variables $Y$ and $W^C$ are conditionally independent given
$W$ if and only if $Y$ and $\Pi_{W^c}X$ are conditionally
independent given $\Pi_W X$, where $\Pi_W$ and $\Pi_{W^C}$ are the
orthogonal projections onto the subspaces spanned by $W$ and
$W^C$, respectively.  If we try to select $r$ variables from
among $m$ explanatory variables, the total number of evaluations
is ${m \choose r}$.

When ${m \choose r}$ is large, we must address the computational cost
that arises in comparing large numbers of subsets.  As in most other
approaches to variable selection~\citep[see, e.g.,][]{Guyon03}, we
propose the use of a greedy algorithm and random search for this
combinatorial aspect of the problem.  (In the experiments presented
in the current paper, however, we confine ourselves to small problems
in which all combinations are tractably evaluated).

We apply this kernel-based method of variable selection to the
\emph{Boston Housing} data \citep{Harrison78} and the \emph{Ozone}
data \citep{Breiman85ACE}, which have been often used as testbed
examples for variable selection.  Tables \ref{tbl:Boston} and
\ref{tbl:Ozone} give the detailed description of the data sets.
There are 506 samples in the Boston Housing data, for which the
variable MV, the median value of house prices in a tract, is
estimated by using the 13 other variables. We use the corrected
version of the data set given by \cite{GilleyPace_Boston96}.  In
the Ozone data in which there are 330 samples, the variable UPO3
(the ozone concentration) is to be predicted by 9 other variables.

Table \ref{tbl:varsel_BosOz} shows the best three sets of four
variables that attain the smallest values of the kernel
generalized variance. For the Boston Housing data, RM and LSTAT
are included in all the three of the result sets in Table
\ref{tbl:varsel_BosOz}, and PTRATIO and TAX are included in two of
them. This observation agrees well with the analysis using
alternating conditional expectation (ACE) by \cite{Breiman85ACE},
which gives RM, LSTAT, PTRATIO, and TAX as the four major
contributors. The original motivation in the study was to
investigate the influence of nitrogen oxide concentration (NOX) on
the house price \citep{Harrison78}.  In accordance with the
previous studies, our analysis shows a relatively small
contribution of NOX. For the Ozone data, all three of the result
sets in the variable selection method include HMDT, SBTP, and
IBHT.  The variables IBTP, DGPG, and VDHT are chosen in one of the
sets.  This shows a fair accordance with earlier results by
\cite{Breiman85ACE} and \cite{LiLueChen00}; the former concludes
by ACE that SBTP, IBHT, DGPG, and VSTY are the most influential,
and the latter selects HMDT, IBHT, and DGPG using a pHd-based
method.



\begin{table}
\centering
  \begin{tabular}{c c l}
  \hline
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
  Variable & & Description \\
  \hline
  MV & \rule[0.5ex]{3ex}{0.5pt} & median value of owner-occupied home \\
  CRIM & \rule[0.5ex]{3ex}{0.5pt} & crime rate by town \\
  ZN  & \rule[0.5ex]{3ex}{0.5pt}& proportion of town's residential land zoned for lots \\
      &  & \hspace{5ex}greater than 25,000 square feet \\
  INDUS  & \rule[0.5ex]{3ex}{0.5pt} &  proportion of nonretail business acres per town \\
  CHAS  & \rule[0.5ex]{3ex}{0.5pt} & Charles River dummy \\
  & & \hspace{5ex}($=1$ if tract bounds the Charles River, 0 otherwise) \\
  NOX  & \rule[0.5ex]{3ex}{0.5pt} & nitrogen oxide concentration in pphm \\
  RM  & \rule[0.5ex]{3ex}{0.5pt} & average number of rooms in owner units \\
  AGE  & \rule[0.5ex]{3ex}{0.5pt} & proportion of owner units build prior to 1940 \\
  DIS  & \rule[0.5ex]{3ex}{0.5pt} & weighted distances to five employment centers \\
  &  & \hspace{5ex} in the Boston region \\
  RAD  & \rule[0.5ex]{3ex}{0.5pt} & index of accessibility to radial highways \\
  TAX  & \rule[0.5ex]{3ex}{0.5pt} & full property tax rate (\$/\$10,000) \\
  PTRATIO  & \rule[0.5ex]{3ex}{0.5pt} & pupil-teacher ratio by town school district \\
  B  & \rule[0.5ex]{3ex}{0.5pt} & black proportion of population \\
  LSTAT & \rule[0.5ex]{3ex}{0.5pt} & proportion of population that is lower status \\
  \hline
\end{tabular}
  \caption{Boston Housing Data}\label{tbl:Boston}
\end{table}


\begin{table}
  \centering
\begin{tabular}{c c l}
  \hline
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
  Variable & & Description \\
  \hline
  UPO3 & \rule[0.5ex]{3ex}{0.5pt} & upland ozone concentration (ppm) \\
  VDHT & \rule[0.5ex]{3ex}{0.5pt} & Vandenburg 500 millibar height (m)\\
  HMDT & \rule[0.5ex]{3ex}{0.5pt} & himidity (percent)\\
  IBHT & \rule[0.5ex]{3ex}{0.5pt} & inversion base height (ft.) \\
  DGPG & \rule[0.5ex]{3ex}{0.5pt} & Daggett pressure gradient (mmhg)  \\
  IBTP & \rule[0.5ex]{3ex}{0.5pt} & inversion base temperature (\farenheit) \\
  SBTP & \rule[0.5ex]{3ex}{0.5pt} & Sandburg Air Force Base temperature (\celsius) \\
  VSTY & \rule[0.5ex]{3ex}{0.5pt} & visibility (miles) \\
  WDSP &  \rule[0.5ex]{3ex}{0.5pt} & wind speed (mph) \\
  DAY  &  \rule[0.5ex]{3ex}{0.5pt} & day of the year\\
  \hline
\end{tabular}
  \caption{Ozone data}\label{tbl:Ozone}
\end{table}




\begin{table}
  \centering
  \begin{tabular}{c|c|c|c}
    \hline
    % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
    Boston & 1st & 2nd & 3rd \\
    \hline
   CRIM &   &   X &  \\
    ZN  &  & & \\
  INDUS &  &  & \\
  CHAS  &  &  & \\
  NOX   &  &  & \\
  RM    & X & X & X \\
  AGE   &  &  & \\
  DIS  &  &   &X \\
  RAD  &  &  & \\
  TAX  & X &   & X \\
  PTRATIO  & X & X & \\
  B     &  & &  \\
  LSTAT &  X & X & X \\  \hline
   KGV &  .1768  & .1770 & .1815 \\
    \hline
  \end{tabular}
  \hspace*{0.5cm}
\begin{tabular}{c|c|c|c}
\multicolumn{4}{c}{ $ $} \\
\multicolumn{4}{c}{$ $} \\
\multicolumn{4}{c}{ $ $} \\
\multicolumn{4}{c}{ $ $} \\
\hline
    Ozone & 1st & 2nd & 3rd \\
\hline
  VDHT &   &   & X  \\
  HMDT & X & X & X  \\
  IBHT & X & X & X  \\
  DGPG &   & X &    \\
  IBTP & X &   &    \\
  SBTP & X & X & X  \\
  VSTY &   &   &    \\
  WDSP &   &   &    \\
  DAY  &   &   &    \\
  \hline
  KGV  & .2727 & .2736 & .2758 \\
  \hline
  \end{tabular}
  \caption{Variable selection using the proposed kernel method.}
  \label{tbl:varsel_BosOz}
\end{table}



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

We have presented KDR, a new kernel-based approach to dimensionality
reduction for regression and classification.   One of the most notable
aspects of this method is its generality---we do not impose any strong
assumptions on either the conditional or the marginal distribution.
This allows the method to be applicable to a wide range of problems,
and gives it a significant practical advantage over existing methods
such as CCA, PPR, SIR, and pHd.  These methods all impose significant
restrictions on the conditional probability, the marginal distribution,
or the dimensionality of the effective subspaces.

Our experiments have shown that the KDR method can provide many of
the desired effects of dimensionality reduction: it provides data
visualization capabilities, it can successfully select important
explanatory variables in regression, and it can yield classification
performance that is better than the performance achieved with the
full-dimensional covariate space.  We have also discussed the extension
of the KDR method to variable selection.  Experiments with classical
data sets has shown an accordance with the previous results on these
data sets and suggest that further study of this application of KDR
is warranted.

The theoretical basis of KDR lies in the nonparametric characterization
of conditional independence that we have presented in this paper.
Extending earlier work on the kernel-based characterization of
independence in ICA~\citep{KernelICA}, we have shown that conditional
independence can be characterized in terms of covariance operators on
a reproducing kernel Hilbert space.  While our focus has been on the
problem of dimensionality reduction, it is also worth noting that there
are many other possible applications of this characterization.  In
particular, conditional independence plays an important role in the
structural definition of probabilistic graphical models, and our
results may have applications to model selection and inference in
graphical models.

There are several statistical problems which need to be addressed in
further research on KDR.  First, a basic analysis of the statistical
consistency of the KDR-based estimator---the convergence of the estimator
to the true subspace when such a space really exists---is needed.
We expect that, to prove consistency, we will need to impose a
condition on the rate of decrease of the regularization coefficient 
$\varepsilon$ as the sample size $n$ goes to infinity.
Second, and most significantly, we need rigorous methods for choosing
the dimensionality of the effective subspace.  If the goal is that of
achieving high predictive performance after dimensionality reduction,
we can use one of many existing methods (e.g., cross-validation,
penalty-based methods) to assess the expected generalization as
a function of dimensionality.  Note in particular that by using KDR as
a method to select an estimator given a fixed dimensionality, we have
substantially reduced the number of hypotheses being considered,
and expect to find ourselves in a regime in which methods such as
cross-validation are likely to be effective.  It is also worth
noting, however, that the goals of dimensionality reduction are
not always simply that of prediction; in particular, the search
for small sets of explanatory variables will need to be guided
by other principles.  Finally, asymptotic analysis may provide
useful guidance for selecting the dimensionality; an example of
such an analysis that we believe can be adopted for KDR has been
presented by \citet{KerChauLi91} for the SIR method.


\acks{
This work was done while the first author was visiting the
University of California, Berkeley.  We thank the action editor and
reviewers for their valuable comments.  In particular, the 
proof of Theorem \ref{thm:Gauss_pd} was suggested by one of the
reviewers.  We also thank Dr.\ Noboru Murata of Waseda University
and Dr.\ Motoaki Kawanabe of Fraunhofer, FIRST for their helpful
comments on the early version of this work. We wish to acknowledge
support from JSPS KAKENHI 15700241, ONR MURI N00014-00-1-0637, NSF
grant IIS-9988642, and a grant from Intel Corporation.}


\appendix

\section{Cross-Covariance Operators on Reproducing Kernel Hilbert Spaces
and Independence of Random Variables}
In this appendix, we present additional background and detailed proofs
for results relating cross-covariance operators to marginal and
conditional independence between random variables.
\subsection{Cross-Covariance Operators}
While cross-covariance operators are generally defined for random
variables on Banach spaces~\cite{ProbBanach,Baker73}, they are
more easily defined on reproducing kernel Hilbert spaces (RKHS).
In this subsection, we summarize some of the basic mathematical
facts used in Sections \ref{sec:cov_op} and \ref{sec:method}.
While we discuss only real Hilbert spaces, extension to the
complex case is straightforward.
\begin{thm}
\label{thm:cov_op} Let $(\Omega_1, {\cal B}_1)$ and $(\Omega_2,
{\cal B}_2)$ be measurable spaces, and let $(\cH_1,k_1)$ and
$(\cH_2, k_2)$ be reproducing kernel Hilbert spaces on $\Omega_1$
and $\Omega_2$, respectively, with $k_1$ and $k_2$ measurable.
Suppose we have a random vector $(X, Y)$ on $\Omega_1 \times
\Omega_2$ such that $E_X[ k_1(X,X) ]$ and $E_Y[ k_2(Y,Y)]$ are
finite. Then, there exists a unique operator $\Sigma_{YX}$ from
$\cH_1$ to $\cH_2$ such that
\begin{equation}
\label{eq:cov_def} \la g, \Sigma_{YX}f \ra_{\cH_2} = E_{XY}[
f(X)g(Y) ] - E_X[f(X)]E_Y[g(Y)]
\end{equation}
holds for all $f\in \cH_1$ and $g\in \cH_2$.  This is called the
\emph{cross-covariance operator}.
\end{thm}
\begin{proof}
Obviously, the operator is unique, if it exists.   From Riesz's
representation theorem~\citep[see][Theorem II.4, for
example]{ReedSimon}, the existence of $\Sigma_{YX}f \in \cH_2$ for
a fixed $f$ can be proved by showing that the right hand side of
\eq{eq:cov_def} is a bounded linear functional on $\cH_2$. The
linearity is obvious, and the boundedness is shown by
\begin{align}
& \bigl|E_{XY}[ f(X)g(Y) ] - E_X[f(X)]E_Y[g(Y)]\bigr| \nonumber \\
& \leq E_{XY}\bigl|\la k_1(\cdot, X), f\ra_{\cH_1} \la
k_2(\cdot,Y),g \ra_{\cH_2} \bigr| + E_X \bigl| \la k_1(\cdot, X),
f \ra_{\cH_1}\bigr|  \cdot E_Y \bigl| \la k_2(\cdot,Y),g\ra_{\cH_2} \bigr|  \nonumber \\
& \leq E_{XY}\bigl[ \|k_1(\cdot, X)\|_{\cH_1} \|f\|_{\cH_1}
\|k_2(\cdot,Y)\|_{\cH_2} \|g\|_{\cH_2}  \bigr] \nonumber \\
& \qquad + E_X\bigl[ \|k_1(\cdot, X)\|_{\cH_1} \|f\|_{\cH_1}\bigr]
E_Y \bigl[ \|k_2(\cdot,Y)\|_{\cH_2} \|g\|_{\cH_2} \bigr]   \nonumber \\
& \leq \bigl\{ E_X[ k_1(X,X)]^{1/2} E_Y[k_2(Y,Y)]^{1/2} + E_X[
k_1(X,X)^{1/2}] E_Y[k_2(Y,Y)^{1/2}]\bigr\}
\|f\|_{\cH_1}\|g\|_{\cH_2}  \nonumber \\
& \leq 2 E_X\bigl[ k_1(X,X) \bigr]^{1/2} E_Y\bigl[
k_2(Y,Y)\bigr]^{1/2} \|f\|_{\cH_1}\|g\|_{\cH_2}  .
\label{eq:cov_bound}
\end{align}
For the second last inequality, $\| k(\cdot,x) \|_\cH^2 = k(x,x)$
is used. The linearity of the map $\Sigma_{YX}$ is given by the
uniqueness part of Riesz's representation theorem.
\end{proof}

\noindent From \eq{eq:cov_bound}, $\Sigma_{YX}$ is bounded,
and by definition, we see $\Sigma_{YX}^* = \Sigma_{XY}$, where
$A^*$ denotes the adjoint of $A$.  If the two RKHS are the same,
the operator $\Sigma_{XX}$ is called the \emph{covariance
operator}. A covariance operator $\Sigma_{XX}$ is bounded,
self-adjoint, and trace-class.


\medskip
In an RKHS, conditional expectations can be expressed by
cross-covariance operators, in a manner analogous to
finite-dimensional Gaussian random variables.

\begin{thm}
\label{thm:cond} Let $(\cH_1, k_1)$ and $(\cH_2, k_2)$ be RKHS on
measurable spaces $\Omega_1$ and $\Omega_2$, respectively, with
$k_1$ and $k_2$ measurable, and $(X,Y)$ be a random vector on
$\Omega_1\times \Omega_2$. Assume that $E_X[k_1(X,X)]$ and
$E_Y[k_2(Y,Y)]$ are finite, and for all $g\in \cH_2$ the
conditional expectation $E_{Y|X}[ g(Y) \mid X=\cdot]$ is an
element of $\cH_1$. Then, we have for all $g\in \cH_2$
\begin{equation*}
\Sigma_{XX} E_{Y|X}[ g(Y) \mid X =\cdot] = \Sigma_{XY} g,
\end{equation*}
where $\Sigma_{XX}$ and $\Sigma_{XY}$ are the covariance and
cross-covariance operator.
\end{thm}
\begin{proof}
For any $f\in \cH_1$, we have
\begin{align}
& \la f, \Sigma_{XX} E_{Y|X}[g(Y) \mid X=\cdot] \ra_{\cH_1}  \nonumber \\
& \qquad = E_{X}\bigl[ f(X) E_{Y|X}[ g(Y) \mid X]\bigr] -
E_{X}[f(X)] E_X\bigl[
E_{Y|X}[g(Y) \mid X]\bigr]     \nonumber \\
& \qquad = E_{XY}[ f(X)g(Y)] - E_{X}[f(X)] E_{Y}[g(Y)] \quad =
\quad \la f, \Sigma_{XY}g \ra_{\cH_1}.  \nonumber
\end{align}
This completes the proof.
\end{proof}
\begin{cor}
\label{cor:cond} Let $\tilde{\Sigma}_{XX}^{-1}$ be the right
inverse of $\Sigma_{XX}$ on $(Ker \Sigma_{XX})^{\perp}$. Under the
same assumptions as Theorem \ref{thm:cond}, we have
\begin{equation*}
\la f, \tilde{\Sigma}_{XX}^{-1} \Sigma_{XY}g \ra = \la f,
E_{Y|X}[g(Y) \mid X=\cdot] \ra
\end{equation*}
for all $f\in(\mathrm{Ker} \Sigma_{XX})^{\perp}$ and $g\in\cH_2$.
In particular, if $\mathrm{Ker} \Sigma_{XX} = 0$, we have
\begin{equation*}
\Sigma_{XX}^{-1} \Sigma_{XY}g  =  E_{Y|X}[g(Y) \mid X=\cdot].
\end{equation*}
\end{cor}
\begin{proof}
Note that the product $ \tilde{\Sigma}_{XX}^{-1} \Sigma_{XY}$ is
well-defined, because $\overline{\mathrm{Range} \Sigma_{XY}}
\subset \overline{\mathrm{Range} \Sigma_{XX}} =
(\mathrm{Ker}\Sigma_{XX})^\perp $. The first inclusion is shown
from the expression
$\Sigma_{XY}=\Sigma_{XX}^{1/2}V\Sigma_{YY}^{1/2}$ with a bounded
operator $V$~\citep[][Theorem 1]{Baker73}, and the second equation
holds for any self-adjoint operator.  Take $f = \Sigma_{XX}h \in
\text{Range}\Sigma_{XX}$. Then, Theorem \ref{thm:cond} yields
\begin{align*}
\la f, \tilde{\Sigma}_{XX}^{-1} \Sigma_{XY} g \ra & = \la h,
\Sigma_{XX}\tilde{\Sigma}_{XX}^{-1}\Sigma_{XX}E_{Y|X}[g(Y) \mid
X=\cdot]
\ra \\
& = \la h, \Sigma_{XX}E_{Y|X}[g(Y) \mid X=\cdot] \ra  = \la
f,E_{Y|X}[g(Y) \mid X=\cdot] \ra.
\end{align*}
This completes the proof.
\end{proof}


The assumption $E_{Y|X}[ g(Y) \mid X=\cdot] \in \cH_1$ in Theorem
\ref{thm:cond} can be simplified so that it can be checked without
reference to a specific $g$.
\begin{prop}
Under the condition of Theorem \ref{thm:cond}, if there exists
$C>0$ such that
\begin{equation*}
E_{Y|X}[ k_2(y_1, Y) \mid X=x_1]E_{Y|X}[ k_2(y_2, Y) \mid X=x_2]
\leq C k_1(x_1,x_2)k_2(y_1,y_2)
\end{equation*}
for all $x_1, x_2 \in \Omega_1$ and $y_1,y_2\in \Omega_2$, then
for all $g\in\cH_2$ the conditional expectation $E_{Y|X}[ g(Y)
\mid X=\cdot]$ is an element of $\cH_1$.
\end{prop}
\begin{proof}
See Theorem 2.3.13 in \cite{Alpay01}.
\end{proof}
From this proposition, it is obvious that $E_{Y|X}[g(Y) \mid
X=\cdot] \in \cH_1$ holds, if the range of $X$ and $Y$ are
bounded.


For a function $f$ in an RKHS, the expectation of $f(X)$ can be
formulated as the inner product of $f$ and a fixed element. Let
$(\Omega, \cal{B})$ be a measurable space, and $(\cH, k)$ be an
RKHS on $\Omega$ with $k$ measurable. Note that for a random
variable $X$ on $\Omega$, the linear functional $f \mapsto
E_X[f(X)]$ is bounded if $E_X[k(X,X)]$ exists.  By Riesz's
theorem, there is $u\in \cH$ such that $\la u,f \ra_\cH
=E_X[f(X)]$ for all $f\in \cH$. If we define $E_X[k(\cdot,X)] \in
\cH$ by this element $u$, we formally obtain the equality
\begin{equation*}
\la E_X[k(\cdot, X)], f\ra_\cH = E_X[ \la k(\cdot, X), f \ra_\cH
],
\end{equation*}
which looks like the interchangeability of the expectation by $X$
and the inner product.  While the expectation $E_X[k(\cdot,X)]$
can be defined, in general, as an integral with respect to the
distribution on $\cH$ induced by $k(\cdot,X)$, the element
$E_X[k(\cdot,X)]$ is formally obtained as above in a reproducing
kernel Hilbert space.


\subsection{Conditional Covariance Operator and Conditional Independence}
We define the conditional (cross-)covariance operator, and derive
its relation with the conditional covariance of random variables.
Let $(\cH_1,k_1)$, $(\cH_2,k_2)$, let $(\cH_3,k_3)$ be RKHS on
measurable spaces $\Omega_1$, $\Omega_2$, and $\Omega_3$,
respectively, and let $(X,Y,Z)$ be a random vector on
$\Omega_1\times\Omega_2\times\Omega_3$. The \emph{conditional
cross-covariance operator of $(X,Y)$ given $Z$} is defined by
\begin{equation*}
\Sigma_{YX|Z} := \Sigma_{YX} -
\Sigma_{YZ}\tilde{\Sigma}_{ZZ}^{-1}\Sigma_{ZX}.
\end{equation*}

Because $\mathrm{Ker}\Sigma_{ZZ} \subset \mathrm{Ker}\Sigma_{YZ}$
from the fact $\Sigma_{YZ} = \Sigma_{YY}^{1/2} V
\Sigma_{ZZ}^{1/2}$ for some bounded operator $V$ \citep[][Theorem
1]{Baker73}, the operator $\Sigma_{YZ}\Sigma_{ZZ}^{-1}\Sigma_{YX}$
can be uniquely defined, even if $\Sigma_{ZZ}^{-1}$ is not unique.
By abuse of notation, we write
$\Sigma_{YZ}\Sigma_{ZZ}^{-1}\Sigma_{ZX}$, when cross-covariance
operators are discussed.


The conditional cross-covariance operator is related to the
conditional covariance of the random variables.
\begin{prop}
\label{prop:cond_indep_pre} Let $(\cH_1, k_1)$, $(\cH_2, k_2)$,
and $(\cH_3,k_3)$ be reproducing kernel Hilbert spaces on
measurable spaces $\Omega_1$, $\Omega_2$, and $\Omega_3$,
respectively, with $k_i$ measurable, and let $(X,Y,Z)$ be a
measurable random vector on $\Omega_1\times \Omega_2\times
\Omega_3$ such that $E_X[k_1(X,X)]$, $E_Y[k_2(Y,Y)]$, and
$E_Z[k_3(Z,Z)]$ are finite.  It is assumed that $E_{X|Z}[f(X) \mid
Z=\cdot]$ and $E_{Y|Z}[g(Y) \mid Z=\cdot]$ are elements of $\cH_3$
for all $f\in\cH_1$ and $g\in\cH_2$.  Then, for all $f\in\cH_1$
and $g\in\cH_2$, we have
\begin{align}
\la g, \Sigma_{YX|Z} f \ra_{\cH_2} & = E_{XY}[ f(X)g(Y)]
-E_{Z}\bigl[ E_{X|Z}[ f(X) \mid Z] E_{Y|Z}[ g(Y) \mid Z] \bigr]
 \nonumber \\
& = E_Z\bigl[ \mathrm{Cov}_{XY|Z}\bigl( f(X),g(Y) \mid Z \bigr)
\bigr].
\end{align}
\end{prop}
\begin{proof}
From the decomposition $\Sigma_{YZ} = \Sigma_{YY}^{1/2} V
\Sigma_{ZZ}^{1/2}$, we have $\Sigma_{ZY}g \in
(\mathrm{Ker}\Sigma_{ZZ})^\perp$.  Then, by Corollary
\ref{cor:cond}, we obtain
\begin{multline*}
\la g, \Sigma_{YZ}\tilde{\Sigma}_{ZZ}^{-1} \Sigma_{ZX} f \ra  =
\la \Sigma_{ZY}g, \tilde{\Sigma}_{ZZ}^{-1} \Sigma_{ZX} f\ra
= \la \Sigma_{ZY}g, E_{X|Z}[f(X) \mid Z] \ra  \\
= E_{YZ}\bigl[ g(Y)E_{X|Z}[ f(X) \mid Z ] \bigr] -
E_X[f(X)]E_Y[g(Y)].
\end{multline*}
From this equation, the theorem is proved by
\begin{align}
\la g, \Sigma_{YX|Z} f \ra & = E_{XY}[
f(X)g(Y)]-E_{X}[f(X)]E_Y[g(Y)] \nonumber \\
& \qquad - E_{YZ}\bigl[ g(Y)E_{X|Z}[ f(X) \mid
Z ] \bigr] + E_X[f(X)]E_Y[g(Y)] \nonumber \\
& = E_{XY}[ f(X)g(Y)] -E_{Z}\bigl[ E_{X|Z}[ f(X) \mid Z] E_{Y|Z}[
g(Y) \mid Z] \bigr].
\end{align}
\end{proof}


The following definition is needed to state our main
theorem.  Let $(\Omega,\mathcal{B})$ be a measurable space, let
$(\cH,k)$ be a RKHS over $\Omega$ with $k$ measurable and bounded,
and let $\mathcal{S}$ be the set of all the probability measures
on $(\Omega,\mathcal{B})$.  The RKHS $\cH$ is called {\it
probability-determining}, if the map
\begin{equation*}
\mathcal{S} \ni P \quad\mapsto\quad (f\mapsto E_{X\sim P}[f(X)])
\in \cH^*
\end{equation*}
is one-to-one, where $\cH^*$ is the dual space of $\cH$. From
Riesz's theorem, $\cH$ is probability-determining if and only if
the map
\[
\mathcal{S} \ni P \quad\mapsto\quad E_{X\sim P}[ k(\cdot,X) ] \in
\cH
\]
is one-to-one.  For Gaussian kernels, the following theorem can be
proved by an argument similar to that used in the proof of Theorem 2 
in \cite{KernelICA} and the uniqueness of the characteristic
function.  For completeness, we present another simple proof here.
\begin{thm}
\label{thm:Gauss_pd}
For an arbitrary $\sigma >0$, the reproducing
kernel Hilbert space $\cH$ with Gaussian kernel $k(x,y)=\exp(-\| x
- y\|^2/\sigma^2 )$ on $\R^m$ is probability-determining.
\end{thm}
\begin{proof}
Suppose $P$ and $Q$ are different probabilities on $\R^m$ such
that $E_{Z\sim P}[f(Z)] = E_{Z\sim Q}[f(Z)]$ for all $f\in \cH$.
Let $y_1$ and $y_2$ be two different vectors in $\R^m$, and $Y$ be
a random variable with probability $1/2$ for each of $Y=y_1$ and
$Y=y_2$. Define a random variable $X$ so that the probability of
$X$ given $Y=y_1$ and $Y=y_2$ are $P$ and $Q$, respectively.
Noting that the marginal distribution of $X$ is $(P+Q)/2$,  we
have for all $f,g\in\cH$,
\begin{align*}
    & E_{X,Y}[f(X)g(Y)] - E_X[f(X)]E_Y[g(Y)] = E_Y\bigl[E_{X|Y}[f(X)
    |Y ] g(Y)\bigr] - E_X[f(X)]E_Y[g(Y)]   \\
    & = \frac{1}{2}g(y_1)E_{X|Y}[f(X)|Y=y_1]
    + \frac{1}{2}g(y_2)E_{X|Y}[f(X)|Y=y_2] \\
    & \hspace{2cm} - \frac{1}{2}\bigl( E_{Z\sim
    P}[f(Z)] + E_{Z\sim Q}[f(Z)] \bigr) \frac{g(y_1)+ g(y_2)}{2}
    \\
    & = 0.
\end{align*}
From Theorem 2 in \cite{KernelICA}, $X$ and $Y$ must be
independent, which contradicts the construction of $X$ and $Y$.
\end{proof}


Recall that for two RKHS $\cH_1$ and $\cH_2$ on $\Omega_1$ and
$\Omega_2$, respectively, the direct product $\cH_1\otimes\cH_2$
is the RKHS on $\Omega_1\times\Omega_2$ with the positive definite
kernel $k_1 k_2$~\citep[see][]{Aronszajn50}.  Note that if the two
RKHS have Gaussian kernels, their direct product is also a RKHS
with Gaussian kernel, and thus probability-determining.  The
relation between conditional independence and the conditional
covariance operator is given by the following theorem:

\begin{thm}
\label{thm:cond_indep1} Let $(\cH_{11},k_{11})$, $(\cH_{12},
k_{12})$, and $(\cH_2, k_2)$ be reproducing kernel Hilbert spaces
on measurable spaces $\Omega_{11}$, $\Omega_{12}$, and $\Omega_2$,
respectively, with continuous and bounded kernels.  Let
$(X,Y)=(U,V,Y)$ be a random vector on
$\Omega_{11}\times\Omega_{12}\times\Omega_2$, where $X=(U,V)$, and
let $\cH_1 = \cH_{11}\otimes\cH_{12}$ be the direct product.  It
is assumed that $E_{Y|U}[g(Y) \mid U=\cdot]\in\cH_{11}$ and
$E_{Y|X}[g(Y) \mid X=\cdot]\in\cH_1$ for all $g\in\cH_2$. Then, we
have
\begin{equation}
\label{eq:var_ineq} \Sigma_{YY|U} \geq \Sigma_{YY|X},
\end{equation}
where the inequality refers to the order of self-adjoint
operators, and if further $\cH_2$ is probability-determining, the
following equivalence holds
\begin{equation}
\label{eq:cond_thm} \Sigma_{YY|X} = \Sigma_{YY|U} \qquad
\Longleftrightarrow \qquad Y \indep V \mid U   .
\end{equation}
\end{thm}
\begin{proof}
The right hand side of \eq{eq:cond_thm} is equivalent to $P_{Y|X}
= P_{Y|U}$, where $P_{Y|X}$ and $P_{Y|U}$ are the conditional
probability of $Y$ given $X$ and given $U$, respectively.  Taking
the expectation of the well-known equality
\begin{equation*}
\text{Var}_{Y|U}[g(Y) \mid U] = E_{V|U}\bigl[ \text{Var}_{Y|U,V}[
g(Y) \mid U,V] \bigr] + \text{Var}_{V|U} \bigl[ E_{Y|U,V}[g(Y)
\mid U,V] \bigr]
\end{equation*}
with respect to $U$, we derive
\begin{equation}
\label{eq:yyz_yyx}
 E_U\bigl[ \text{Var}_{Y|U}[ g(Y) \mid U] \bigr] = E_X\bigl[
\text{Var}_{Y|X}[ g(Y) \mid X] \bigr] + E_U\bigl[
\text{Var}_{V|U}[ E_{Y|X}[g(Y) \mid X]]\bigr].
\end{equation}
Since the last term of \eq{eq:yyz_yyx} is nonnegative, we obtain
\eq{eq:var_ineq} from Proposition \ref{prop:cond_indep_pre}.

Equality holds if and only if $\text{Var}_{V|U}[ E_{Y|X}[g(Y) \mid
X]]=0$ for almost every $U$, which means $E_{Y|X}[g(Y) \mid X]$
does not depend on $V$ almost surely.  This is equivalent to
\begin{equation*}
 E_{Y|X}[ g(Y) \mid X ] = E_{Y|U}[ g(Y) \mid U ]
\end{equation*}
for almost every $V$ and $U$.  Because $\cH_2$ is
probability-determining, this means $P_{Y|X} = P_{Y|U}$.
\end{proof}




\subsection{Conditional Cross-Covariance Operator and Conditional Independence}
Theorem \ref{thm:cond_indep1} characterizes conditional
independence using the conditional covariance operator.  Another
formulation is possible with a conditional cross-covariance
operator.

Let $(\Omega_1,\mathcal{B}_1)$, $(\Omega_2,\mathcal{B}_2)$, and
$(\Omega_3,\mathcal{B}_3)$ be measurable spaces, and let $(X,Y,Z)$
be a random vector on $\Omega_1\times\Omega_2\times\Omega_3$ with
law $P_{XYZ}$.  We define a probability measure $E_{Z}[
P_{X|Z}\otimes P_{Y|Z} ]$ on $\Omega_1\times\Omega_2$ by
\begin{equation*}
E_{Z}[ P_{X|Z}\otimes P_{Y|Z} ] (A\times B) = E_Z \bigl[ E_{X|Z}[
\chi_A|Z] \, E_{Y|Z}[ \chi_B \mid Z] \bigr],
\end{equation*}
where $\chi_A$ is the characteristic function of a measurable set
$A$. It is canonically extended to any product-measurable sets in
$\Omega_1\times\Omega_2$.
\begin{thm}
\label{thm:cond_indep2} Let $(\Omega_i,\mathcal{B}_i)$ ($i=1,2,3$)
be a measurable space, let $(\cH_i, k_i)$ be a RKHS on $\Omega_i$
with kernel measurable and bounded, and let $(X,Y,Z)$ be a random
vector on $\Omega_1\times \Omega_2\times \Omega_3$.  It is assumed
that $E_{X|Z}[f(X) \mid Z=\cdot]$ and $E_{Y|Z}[g(Y) \mid Z=\cdot]$
belong to $\cH_3$ for all $f\in\cH_1$ and $g\in\cH_2$, and that
$\cH_1\otimes\cH_2$ is probability-determining. Then, we have
\begin{equation}
\label{eq:cond_indepw} \Sigma_{YX|Z} = O \qquad
\Longleftrightarrow \qquad P_{XY} = E_{Z}[ P_{X|Z}\otimes P_{Y|Z}
].
\end{equation}
\end{thm}
\begin{proof}
The right-to-left direction is trivial from Theorem
\ref{prop:cond_indep_pre} and the definition of $E_{Z}[
P_{X|Z}\otimes P_{Y|Z} ]$.  The left-hand side yields $E_Z[
E_{X|Z}[f(X) \mid Z]E_{Y|Z}[g(Y) \mid Z]] = E_{XY}[f(X)g(Y)]$ for
all $f\in\cH_1$ and $g\in\cH_2$.   Because $\cH_1\otimes \cH_2$ is
defined as the completion of all the linear combinations of
$f_i(x)g_i(y)$ for $f_i\in\cH_1$ and $g_i\in\cH_2$, we have
$E_{(X',Y')\sim Q}[h(X',Y')] = E_{XY}[h(X,Y)]$ for every such
linear combination, and thus every $h\in\cH_1\otimes\cH_2$ as a
limit, where $Q= E_Z[ P_{X|Z}\otimes P_{Y|Z}]$.  This implies the
right-hand side, because $\cH_1\otimes\cH_2$ is
probability-determining.
\end{proof}

The right-hand side of \eq{eq:cond_indepw} is weaker than the
conditional independence of $X$ and $Y$ given $Z$.  However, if
$Z$ is a part of $X$, we obtain conditional independence.

\begin{cor}
Let $(\cH_{11},k_{11})$, $(\cH_{12}, k_{12})$, and $(\cH_2, k_2)$
be reproducing kernel Hilbert spaces on measurable spaces
$\Omega_{11}$, $\Omega_{12}$, and $\Omega_2$, respectively, with
kernels measurable and bounded.  Let $(X,Y)=(U,V,Y)$ be a random
vector on $\Omega_{11}\times\Omega_{12}\times\Omega_2$, where
$X=(U,V)$, and let $\cH_1 = \cH_{11}\otimes\cH_{12}$ be the direct
product.  It is assumed that $E_{X|U}[f(X) \mid U=\cdot]$ and
$E_{Y|U}[g(Y) \mid U=\cdot]$ belong to $\cH_{11}$ for all
$f\in\cH_1$ and $g\in\cH_2$, and $\cH_1\otimes\cH_2$ is
probability-determining. Then, we have
\begin{equation}
\label{eq:cond_indep3} \Sigma_{YX|U}=O \qquad \Longleftrightarrow
\qquad Y\indep V \,|\,U.
\end{equation}
\end{cor}
\begin{proof}
For any measurable sets $A\subset \Omega_{11}$, $B\subset
\Omega_{12}$, and $C\subset\Omega_2$, we have, in general,
\begin{align}
& E_{U}\bigl[E_{X|U}[\chi_{A\times B}(U,V) \mid U]
E_{Y|U}[\chi_{C}(Y) \mid U] \bigr] - E_{XY}[\chi_{A\times
B}(U,V)\chi_{C}(Y)]
\nonumber \\
& = E_U\bigl[ E_{V|U}[\chi_{B}(V) \mid U] \chi_{A}(U)
E_{Y|U}[\chi_{C}(Y) \mid U]\bigr]  - E_U\bigl[
E_{VY|U}[\chi_{B}(V)\chi_C(Y) \mid U] \chi_{A}(U) \bigr] \nonumber \\
&= \int_A \bigl\{ P_{V|U}(B \mid u)P_{Y|U}(C \mid u)  - P_{VY \mid
U}(B\times C
 \mid  u)\bigr\} dP_U(u). \label{eq:intA}
\end{align}
From Theorem \ref{thm:cond_indep2}, the left-hand side of
\eq{eq:cond_indep3} is equivalent to $E_U[P_{X|U}\otimes
P_{Y|U}]=P_{XY}$, which implies that the last integral in
\eq{eq:intA} is zero for all $A$. This means $P_{V|U}(B \mid
u)P_{Y|U}(C \mid u)  - P_{VY \mid U}(B\times C \mid u)=0$ for
almost every $u$-$P_U$.  Thus, $Y$ and $V$ are conditional
independent given $U$. The converse is trivial.
\end{proof}
Note that the left-hand side of \eq{eq:cond_indep3} is not
$\Sigma_{YV|U}$ but $\Sigma_{YX|U}$, which is defined on the
direct product $\cH_{11}\otimes\cH_{12}$.

\bibliography{fukubib}





\end{document}
