\documentclass[twoside,11pt]{article}
\usepackage{jmlr2e}
\usepackage{amsmath}            % higher mathematics
\usepackage{amsfonts}           % more math fonts
\usepackage{array}              % extensions to tabular and array 
\usepackage{tabularx}           % support correct-width tables
\usepackage{times}              % typeset in Times Roman
\usepackage{color}      % color support

\newif\ifpdf
\ifx\pdfoutput\undefined
  \pdffalse
  \newcommand{\mypdfoption}{{}}
\else
  \pdfoutput=1
  \pdftrue
  \newcommand{\mypdfoption}{pdftex}
  \pdfcompresslevel=9
\fi

\ifpdf
  \DeclareGraphicsExtensions{.jpg,.png,.pdf,.mps} % .png before .pdf ==> okay!
\else
  \DeclareGraphicsExtensions{.jpg,.eps,.png,.mps} % .eps before .png ==> okay!
\fi

\newcommand{\rset}{\ensuremath{\mathbb{R}}}

%%% The set of non-negative reals

\newcommand{\trainset}{{\mathcal{D}}}
\newcommand{\testset}{{\mathcal{T}}}

\newcommand{\defequal}{\stackrel{\Delta}{=}}

\newcommand{\Var}{\mbox{$\mathrm{Var}$}}
\newcommand{\Cov}{\mbox{$\mathrm{Cov}$}}
\newcommand{\Bias}{\mbox{$\mathrm{Bias}$}}


\newcommand{\PE}{\mbox{$\mathrm{PE}$}}
\newcommand{\EPE}{\mbox{$\mathrm{EPE}$}}
\newcommand{\hatPE}{\mbox{$\widehat{\PE}$}}
\newcommand{\hatEPE}{\mbox{$\widehat{\EPE}$}}
\newcommand{\CV}{\mbox{$\mathrm{CV}$}}
\newcommand{\JK}{\mbox{$\mathrm{JK}$}}

% \newtheorem{corollary}{Corollary}
% \newtheorem{lemma}{Lemma}
% \newtheorem{theorem}{Theorem}
\newtheorem{experiment}{Experiment}

% Heading arguments are {volume}{year}{pages}{submitted}{published}{author-full-names}

\jmlrheading{5}{2004}{1089--1105}{05/03; Revised 9/03}{9/04}{Yoshua Bengio and Yves Grandvalet}

% Short headings should be running head and authors last names

\ShortHeadings{Variance of K-Fold Cross-Validation}{Bengio and Grandvalet}

\firstpageno{1089}

\begin{document}
\sloppy
\title{No Unbiased Estimator of the Variance of K-Fold Cross-Validation}

\author{\name Yoshua Bengio \email bengioy@iro.umontreal.ca \\
        \addr Dept. IRO, Universit{\'e} de Montr{\'e}al \\
              C.P. 6128, Montreal, Qc, H3C 3J7, Canada
        \AND      
        \name Yves Grandvalet \email Yves.Grandvalet@utc.fr \\
        \addr Heudiasyc, UMR CNRS 6599 \\
              Universit\'e de Technologie de Compi\`egne, France}

\editor{Dana Ron}

\maketitle

%\tableofcontents

\begin{abstract}
  
Most machine learning researchers perform quantitative experiments to estimate generalization error and compare the performance of different algorithms (in particular, their proposed algorithm). In order to be able to draw statistically convincing conclusions, it is important  to estimate the uncertainty of such estimates.
This paper studies the very commonly used K-fold cross-validation estimator of generalization performance. The main theorem shows that there exists no universal (valid under all distributions) unbiased estimator of the variance of K-fold cross-validation. The analysis that accompanies this result is based on the eigen-decomposition of the covariance matrix of errors, which has only three different eigenvalues corresponding to three degrees of freedom of the matrix and three components of the total variance. This analysis helps to better understand the nature of the problem and how it can make naive estimators (that don't take into account the error correlations due to the overlap between training and test sets) grossly underestimate variance. This is confirmed by numerical experiments in which the three components of the variance are compared when the difficulty of the learning problem and the number of folds are varied.

\end{abstract}

\begin{keywords}
cross-validation, variance estimators, k-fold cross-validation,
statistical comparisons of algorithms
\end{keywords}


%------------------------------------------------------------------------------%
\section{Introduction}
%------------------------------------------------------------------------------%


In machine learning, the standard measure of accuracy for trained
models is the prediction error (\PE), i.e. the expected loss on future
examples.  Learning algorithms themselves are often compared according to their average
performance, which is formally defined as the expected value of prediction error
(\EPE) over training sets.

When the data distribution is unknown, \PE\ and \EPE\ cannot be
computed.  If the amount of data is large enough, \PE\ can be
estimated by the mean error over a hold-out test set.  The usual
variance estimates for means of independent samples can then be
computed to derive error bars on the estimated prediction error, and
to assess the statistical significance of differences between models.

The hold-out technique does not account for the variance with respect
to the training set, and may thus be considered inappropriate for the
purpose of algorithm comparison~\citep{Dietterich98}.  Moreover, it
makes an inefficient use of data which forbids its application to
small sample sizes.  In this situation, one rather uses computer
intensive resampling methods such as cross-validation or bootstrap to
estimate \PE\ or \EPE.

We focus here on K-fold cross-validation.
While it is known that cross-validation provides an unbiased estimate of \EPE, it is also known that its variance may be very large~\citep{Breiman96}. 
This variance should be estimated to provide faithful confidence intervals on \PE\ or \EPE, and to test the significance of observed differences between algorithms.
This paper provides theoretical arguments showing the difficulty of this estimation.

The difficulties of the variance estimation have already been addressed \citep{Dietterich98,Kohavi95,Nadeau03}.
Some distribution-free bounds on the deviations of cross-validation are available, but they are specific to some locally defined decision rules, such as nearest neighbors \citep{Devroye96}.
This paper builds upon the work of \citet{Nadeau03}, 
which investigated in detail the theoretical and practical merits of several estimators of the variance of cross-validation.
Our analysis departs from this work in the sampling procedure defining the cross-validation estimate.
While \citet{Nadeau03} consider K independent training and test splits, 
we focus on the standard K-fold cross-validation procedure, 
where there is no overlap between test sets: 
each example of the original data set is used once and only once as a test example. 

This paper is organized as follows.  Section \ref{sec:framework}
defines the measures of performance for algorithms, their estimation
by K-fold cross-validation and similar procedures such as delete-$m$
jackknife.  Our theoretical findings are summarized in Sections
\ref{sec:decomposition}--\ref{sec:possiblevalues}.  They are followed
in Section~\ref{sec:experiments} by experiments illustrating the
effect of experimental conditions on the total variance and its
decomposition in three components, and confirming the underestimation
of variance obtained by the naive estimator commonly used by
researchers.

%------------------------------------------------------------------------------%
\section{General Framework}\label{sec:framework}
%------------------------------------------------------------------------------%

In machine learning, the performance measure differs according to the
experimenter's viewpoint.  In applications, we are interested in
finding the best algorithm for solving the particular task at hand,
specified by one particular training set and some information about
the data generating process.  In algorithm evaluation, we want to
compare several learning algorithms for different learning tasks, and
we care about the sensitivity of the learning algorithm to the
choice of training examples.

%------------------------------------------------------------------------------%
\subsection{Measures of Performance}
%------------------------------------------------------------------------------%

Formally, we have a training set
$D=\{\mathbf{z}_1,\ldots,\mathbf{z}_n\}$, with $\mathbf{z}_i \in \cal
Z$, independently sampled from an unknown distribution $P$.  We also
have a learning algorithm $A$, which maps a data set of (almost)
arbitrary size to a function $A: {\cal Z}^* \rightarrow {\cal F}$.
Throughout this paper, we consider symmetric algorithms, i.e. $A$ is
insensitive to the ordering of examples in the training set $D$.  The
discrepancy between the prediction and the observation $\mathbf{z}$ is
measured by a loss functional $L: {\cal F} \times {\cal Z}\rightarrow
\rset$.  Typically, $L$ is the quadratic loss in regression 
($L(f,(x,y)) = (f(x)-y)^2$) and the
misclassification $\{0,1\}$-loss in classification ($(L(f,(x,y)) = 1_{f(x)\neq y}$).

Let $f=A(D)$ be the function returned by algorithm $A$ on the training set $D$. 
In application based evaluation, the goal of learning is usually stated as the minimization of the prediction error, 
i.e. the expected loss on future test examples
%
\begin{equation}
\PE(D) = E[L(f,\mathbf{z})], 
\end{equation}
%
where the expectation is taken with respect to $\mathbf{z}$ sampled from $P$.\footnote{
Note that we are using the same notation for random variables and their realization.
The intended meaning will be specified when not clear from the context.
}

In algorithm based evaluation, we are not really interested in performances on a specific training set;
we would like comparisons on a more general basis.
In this context, the lowest level of generality can be stated as ``training sets of size $n$ sampled from $P$'', 
and the performance of learning algorithm $A$ can be measured by the expected performance of the functions returned in this situation
%
\begin{equation}
\EPE(n) = E[L(A(D),\mathbf{z})] , 
\end{equation}
%
where the expectation is taken with respect to $D$ sampled from $P^n$ and $\mathbf{z}$ independently sampled from $P$.

Note that other types of performances measure can be proposed, 
based for example on parameters, 
or defined by the predictability in other frameworks,
such as the prequential analysis \citep{Dawid97}.

When the data distribution is unknown, \PE\ and \EPE\ cannot be computed. 
They have to be estimated, and it is often crucial to assess the uncertainty attached to this estimation:
\begin{itemize}
  \item
  in application-oriented experiments, to give a confidence interval on \PE;
  \item
  in algorithm-oriented experiments, to take into account the stability of a given algorithm.  
  For comparisons between algorithms, it is essential to assess the statistical significance of observed differences in the estimate \hatEPE.
\end{itemize}
Although this point is often overlooked, estimating the variance of the estimates \hatPE\ and \hatEPE\ requires caution. 

%------------------------------------------------------------------------------%
\subsection{Hold-Out Estimates of Performance}\label{sec:framework:holdout}
%------------------------------------------------------------------------------%

If the amount of data is large enough, \PE\ can be estimated by the mean error over a hold-out test set, 
and the usual variance estimate for means of independent variables can then be computed. 
However, even in the ideal situation where several independent training and test sets would be available, 
this estimate should not be applied to compute the variance of \hatEPE:
even though training and test examples are independent, the test errors are correlated,
since many test errors are computed for each training set, now considered as a random variable.

Figure \ref{fig:illustration2} illustrates how crucial it is to take these correlations into account.
The mean of two variance estimators is reported {\em vs.} the empirical variance of the hold-out estimate,
in an ideal situation where 10 independent training and test sets are available.
The variance of $\hatEPE(n)$ (estimated on 100,000 independent experiments) is displayed for reference by the dotted line.
The average of $\widehat{\theta_1}$, the variance estimator ignoring correlations, shows that this estimate is highly biased,
even for large sample sizes,
whereas the variance estimator $\widehat{\theta_2}$, taking into account correlations, is unbiased.
The details of this experiment are given below.


\begin{figure}[htb]
 \begin{center} 
   \includegraphics[scale = 0.6]{illustration0.ps} 
  \caption{Estimates of the variance of $\hatEPE(n)$ vs. empirical variance of $\hatEPE(n)$ (shown by bold curve) on 100,000 experiments.
           The average of the variance estimators $\widehat{\theta_1}$ (ignoring correlations, dashed curve) and $\widehat{\theta_2}$ 
           (taking into account correlations, dotted curve) are displayed for different training sample size $n$.}\label{fig:illustration2} 
 \end{center} 
\end{figure}

\begin{experiment}\label{exp:hold-out}
Ideal hold-out estimate of \EPE.

We have $K=10$ independent training sets  $D_1,\ldots,D_K$ of $n$ independent examples $\mathbf{z}_i = (\mathbf{x}_i,y_i)$, 
where $\mathbf{x_i} = (x_{i1},\ldots,x_{id})'$ is a $d$-dimensional centered, unit covariance Gaussian variable ($d=30$), 
$y_i = \sqrt{3/d}\sum_{k=1}^d x_{ik} + \varepsilon_i$ with $\varepsilon_i$ being independent, centered, 
unit variance Gaussian variables.\footnote{
The $\sqrt{3/d}$ factor provides an $R^2$ of approximately $3/4$.
}
We also have $K$ independent test sets $T_1,\ldots,T_K$ of size $n$ sampled from the same distribution. 

The learning algorithm consists in fitting a line by ordinary least squares, 
and the estimate of \EPE\ is the average quadratic loss on test examples
$\hatEPE\ = \bar{L} = \frac{1}{K} \sum_{k=1}^K  \frac{1}{n} \sum_{\mathbf{z}_i\in T_k} L_{ki}$, 
where $L_{ki}=L(A(D_k),\mathbf{z}_i)$.

The first estimate of variance of \hatEPE\ is
$\widehat{\theta_{1}} = \frac{1}{K n (K n-1)} \sum_{k=1}^K \sum_{i} (L_{ki}-\bar{L})^2$,
which is unbiased provided there is no correlation between test errors.
The second estimate is
$\widehat{\theta_{2}} = \frac{1}{K (K-1) n^2} \sum_{k=1}^K \sum_{i,j} (L_{ki}-\bar{L})(L_{kj}-\bar{L})$, 
which takes into account correlations between test errors. 
\end{experiment}

Looking at Figure~\ref{fig:illustration2} suggests that asymptotically
the naive estimator of variance converges to the true variance. This
can be shown formally by taking advantage of the results in this paper,
as long as the learning algorithm converges as the amount of training
data goes to infinity (i.e. as $n \rightarrow \infty$ the function
$A(D)$ obtained does not depend on the particular training set $D$).
In that limit, the correlations between test errors converge to 0. The
rate of convergence will depend on the stability of the learning
algorithm as well as on the nature of the data distribution (e.g., the
presence of thick tails and outliers will slow down convergence).

The hold-out technique makes an inefficient use of data which forbids its application in most 
real-life applications with small samples.
Then, K-fold cross-validation can provide estimates of \PE\ or \EPE. 

%------------------------------------------------------------------------------%
\subsection{K-Fold Cross-Validation Estimates of Performance}\label{sec:framework:kfoldxval}
%------------------------------------------------------------------------------%

Cross-validation is a computer intensive technique, using all available examples as training and test examples.
It mimics the use of training and test sets by repeatedly training the algorithm $K$ times with a fraction $1/K$ of training examples left out for testing purposes. 
This kind of hold-out estimate of performance lacks computational efficiency due to the repeated training,
but the latter are meant to lower the variance of the estimate \citep{Stone74}.


In practice, the data set $D$ is first chunked into $K$ disjoint subsets (or {\em blocks}) of the same 
size\footnote{
To simplify the analysis below we assume that $n$ is a multiple of $K$.
}
$m \defequal n/K$.
Let us write $T_k$ for the $k$-th such block, and $D_k$ the training set obtained by removing the elements in $T_k$ from $D$. 
The cross-validation estimator is defined as the average of the errors on test block $T_k$ obtained when the training set is deprived from $T_k$:
\begin{equation}
  \CV(D) = \frac{1}{K} \sum_{k=1}^K \frac{1}{m} \sum_{\mathbf{z}_i \in T_k} L(A(D_k),\mathbf{z}_i) .
\label{eq:xv1}
\end{equation}

% So what ? 

Does \CV\ estimate \PE\ or \EPE? 
Such a question may seem pointless considering that $\PE(D)$ is an estimate of $\EPE(n)$, 
but it becomes relevant when considering the variance of \CV:
does it inform us of the uncertainty about \PE\ or \EPE?

On the one hand, only one training set, $D$, enters the definition of \CV, 
which can be, up to an approximation, an unbiased estimate of $\PE(D)$ \citep{Hastie90}.\footnote{
More precisely, following~\citet{Hastie90}, 
when $L$ is the quadratic loss, and writing $f=A(D)$, $f^{-k}=A(D_k)$,
assuming that for $(\mathbf{x}_i,y_i) = \mathbf{z}_i \in T_k$, $\frac{1}{K}\sum_{k=1}^K f^{-k}(\mathbf{x}_i) \approx f(\mathbf{x}_i)$ (which is weaker than  $f^{-k} \approx f$) yields $E[\CV] \approx E[\frac{1}{n}\sum_{i=1}^n (f(\mathbf{x}_i)-y_i)^2]$, where the expectation is taken with respect to $y_1,\ldots,y_n$.
} 
Some distribution-free bounds on the expected deviations of $|\CV(D)-PE(D)|$ are available for leave-one-out cross-validation applied to specific algorithms $A$ such as nearest neighbors \citep{Devroye96}.
In a more general context, it has also been proved that, under suitable stability assumptions on the algorithm $A$, 
$\CV(D)$ estimates $\PE(D)$ at least as accurately as the training error \citep{Kearns96,Anthony98}. 
A more appealing result states that \CV\ is a more accurate estimate of \PE\ than hold-out testing \citep{Blum99}.
However, this statement does not apply to $\PE(D)$, but to the prediction error of a randomized algorithm picking solutions uniformly within $\{A(D_k)\}_{k=1}^K$.

On the other hand, \CV\ is explicitly defined from the learning algorithm $A$,
and not from the function $f=A(D)$.
The inner average in the definition of \CV\ (\ref{eq:xv1}) is an average test loss for $A(D_k)$ which thus estimates unbiasedly $\PE(D_k)$.
The training sets $D_1,\ldots,D_K$ are clearly not independent, but they are sampled from $P^{n-m}$. 
Hence, the outer average of (\ref{eq:xv1}) estimates unbiasedly $\EPE(n-m)$.\footnote{
Note that leave-one-out cross-validation is known to fail to estimate \EPE\ for unsmooth statistics \cite[e.g.][]{Breiman96,Efron93}.
This failure is due to the similarity of the training sets $D_1,\ldots,D_K$ which are far from being representative samples drawn from $P^{n-m}$.
}
Here, following \citet{Dietterich98} and \citet{Nadeau03}, we will adopt this latter point of view.

The variance estimate of \hatEPE\ provided by the hold-out estimate has to account for test error dependencies
due to the choice of training set, which  cannot be estimated using a single training/test
experiment.
Here, the situation is more complex, since there are additional dependencies due to the overlapping training sets $D_1,\ldots,D_K$.
Before describing this situation in detail and summarizing the results of our theoretical analysis in Sections \ref{sec:decomposition}--\ref{sec:possiblevalues},
we detail some procedures similar to K-fold cross-validation, for which the forthcoming analysis will also hold. 

%------------------------------------------------------------------------------%
\subsection{Other Estimates of the K-Fold Cross-Validation Type}
%------------------------------------------------------------------------------%

% Extensions

One of the main use of variance estimates of \hatEPE\ is to compare learning algorithms.
The analysis presented in this paper also applies to the version of cross-validation dedicated to this purpose:
if we want to compare the performances of algorithms $A_1$ and $A_2$, cross-validation with matched pairs should be the method of choice
%
\begin{equation}\label{eq:xv2}
  \Delta\CV(D) = \frac{1}{K} \sum_{k=1}^K \frac{1}{m} \sum_{\mathbf{z}_i \in T_k} L(A_1(D_k),\mathbf{z}_i) - L(A_2(D_k),\mathbf{z}_i) .
\end{equation}
%
Compared to the difference of two independent cross-validation estimates, $\Delta\CV$ avoids the additional variability due to train/test splits.

In application oriented experiments, we would like to estimate $\PE(D)$, the expected
error when training with the given $D$.
We have seen in Section \ref{sec:framework:kfoldxval} that under stability assumptions,
\CV\ can be used to estimate \PE.
Alternatively, 
we may resort to the jackknife or the delete-$m$ jackknife (see e.g. \citet{Efron93}) to estimate 
% the bias of the training error (i.e the sample mean error of $A(D)$ on training examples) used as an estimate of $\PE(D)$.
the optimism (i.e. the bias of the mean error on training examples, when the latter is used to estimate $\PE(D)$).
Ideally, the estimate of optimism should be an average over all subsets of size $n-m$, 
but a less computationally intensive alternative is
%
\begin{equation}
%  \JK = \frac{1}{n} \sum_{i=1}^nL(A(D),\mathbf{z}_i) - 
  (K-1) \left( \frac{1}{K(n-m)} \sum_{k=1}^K 
  \sum_{\mathbf{z}_i\in D_k} L(A(D_k),\mathbf{z}_i) -
  \frac{1}{n} \sum_{i=1}^n L(A(D),\mathbf{z}_i) \right)  .
\label{eq:opt}
\end{equation}

The link with cross-validation is exhibited more clearly by the following expression of the (debiased) jackknife estimate of \PE
%
\begin{equation}
  \JK = \CV + \frac{1}{n} \sum_{k=1}^K \sum_{i=1}^n \left( L(A(D),\mathbf{z}_i) - L(A(D_k),\mathbf{z}_i) \right) .
\label{eq:jk2}
\end{equation}
%
For additional information about jackknife estimates and clues on the derivation of (\ref{eq:opt}) and (\ref{eq:jk2}),
the reader is referred to \citet{Efron93}.

%------------------------------------------------------------------------------%
\subsection{Generic Notations}\label{sec:framework:notations}
%------------------------------------------------------------------------------%

This paper studies the variance of statistics such as \CV, $\Delta\CV$ or \JK. 
In what follows, these statistics will be denoted by $\hat{\mu}$, a generic notation for means of observations $e_i$ split in $K$ groups.

%
\begin{eqnarray*}
  \hat{\mu} & = & \frac{1}{n} \sum_{i=1}^n e_i  \\
            & = & \frac{1}{K} \sum_{k=1}^K \frac{1}{m} \sum_{i \in T_k} e_i ,
\end{eqnarray*}
%
where, slightly abusing notation, $i \in T_k$ means $\mathbf{z}_i\in T_k$ and
%
\begin{equation*}
  \forall i \in T_k,\  e_i = \left\{\begin{array}{ll}
  L(A(D_k),\mathbf{z}_i) & \mbox{for} \ \hat{\mu}  = \CV , \\
  L(A_1(D_k),\mathbf{z}_i) - L(A_2(D_k),\mathbf{z}_i)& \mbox{for} \ \hat{\mu}  = \Delta\CV , \\
  K L(A(D),\mathbf{z}_i) - \sum_{\ell\neq k}L(A(D_\ell),\mathbf{z}_i) & \mbox{for} \ \hat{\mu}  = \JK .
  \end{array}\right.
\end{equation*}
%

% Estimating $\theta \defequal \Var[\hat{\mu}]$ provides a useful means
% to build a confidence interval or to test a hypothesis on $\mu$. 
% Furthermore, as 
Note that $\hat{\mu}$ is the average of identically distributed
(dependent) variables. 
Thus, it asymptotically converges to a normally distributed variable,
which is completely characterized by its expectation $E[\hat{\mu}]$
% (i.e. \EPE\ for \CV\ and an approximation of \PE\ for \JK) 
and its variance $\Var[\hat{\mu}] = E[\hat{\mu}^2]-E[\hat{\mu}]^2$. 
% The following Section describes the actual structure of the
% covariance of $e_{i}$, which should be taken into account in the
% estimation process.

%------------------------------------------------------------------------------%
\section{Structure of the Covariance Matrix}\label{sec:decomposition}
%------------------------------------------------------------------------------%

The variance of $\hat{\mu}$ is defined as follows
\[
 \theta % \Var[\hat{\mu}] 
 = \frac{1}{n^2}\sum_{i,j} \Cov(e_i,e_j) ,
\]
where $\Cov(e_i,e_j)=E[e_i e_j]-E[e_i]E[e_j]$ is the covariance between variables $e_i$ and $e_j$.

By using symmetry arguments over permutations of the examples in $D$, 
we show that many distributions on $e_i$ and pairwise joint distributions on $(e_i,e_j)$ are identical. 
As a result, the covariance matrix $\boldsymbol{\Sigma}$ has a very particular block structure, 
with only three possible values for $\Sigma_{ij} = \Cov(e_i,e_j)$, and the expression of $\theta$ % $\Var[\hat{\mu}]$ 
is thus a linear combination of these three values. 

\begin{lemma}
\label{lemma:distribution}

Using the notation introduced in section~\ref{sec:framework:notations}, 
\begin{enumerate}
\item all $e_i$ are identically distributed:  

there exists $f$ such that, $\forall i, \ P(e_i=u) = f(u)$.
\item all pairs $(e_i,e_j)$ belonging to the same test block are jointly identically distributed:

there exists $g$ such that, $\forall (i,j) \in T_k^2:j\neq i, \ P(e_i=u,e_j=v) = g(u,v)$.

\item all pairs $(e_i,e_j)$ belonging to different test blocks are jointly identically distributed:

there exists $h$ such that, $\forall i \in T_k,\ \forall j\in T_\ell:\ell \neq k, \ P(e_i=u,e_j=v) = h(u,v)$.

\end{enumerate}

\noindent
{\bf Proof}

\noindent
These results are derived immediately from the permutation-invariance of $P(D)$ and the symmetry of $A$.
\begin{itemize}
\item invariance with respect to permutations within test blocks:
   \begin{enumerate}
   \item  $\forall (i,i')\in T_k^2, \ P(e_i=u) = P(e_{i'}=u) = f_k(u)$; 
   
   \item []
   $\forall (i,i')\in T_k^2, \ \forall j\in T_\ell$:
   
   $P(e_i=u,e_j=v) = P(e_{i'}=u,e_j=v)$
   
   hence:
      
   \item 

   $\forall (i,j) \in T_k^2:j\neq i, \ P(e_i=u,e_j=v) = g_k(u,v)$.
   
   \item 
   
   $\forall i \in T_k,\ \forall j\in T_\ell:\ell \neq k, \ P(e_i=u,e_j=v) = h_{k\ell}(u,v)$.
   \end{enumerate}
\item invariance with respect to permutations between test blocks.
   \begin{enumerate}
   \item $\forall (k,k'), \ f_k(u)= f_{k'}(u)=f(u)$;
   \item $\forall (k,k'), \ g_k(u,v)=g_{k'}(u,v)=g(u,v)$;
  \item $\forall (k,k'), \ \forall (\ell,\ell'): \ell\neq k, \ell\neq k', \ell'\neq k, \ell'\neq k', \ h_{k\ell}(u,v)=h_{k\ell'}(u,v)=h_{k'\ell'}(u,v)=h_{k'\ell}(u,v)=h(u,v)$.
   \end{enumerate}
\end{itemize}

\flushright
Q.E.D.
\end{lemma}


\begin{corollary}
\label{corollary:covariances}
The covariance matrix $\boldsymbol{\Sigma}$ of cross-validation errors $\mathbf{e}=(e_1,\ldots,e_n)'$ has the simple block structure depicted in Figure~\ref{fig:Covariance}:
\begin{enumerate}
\item all diagonal elements are identical

$\forall i, \ \Cov(e_i,e_i) = \Var[e_i] = \sigma^2$;
\item all the off-diagonal entries of the $K$ $m \times m$ diagonal blocks are identical

$\forall (i,j) \in T_k^2:j\neq i,\ \Cov(e_i,e_j) = \omega$;
\item all the remaining entries are identical

$\forall i \in T_k,\ \forall j\in T_\ell:\ell \neq k, \ \Cov(e_i,e_j) = \gamma$.
\end{enumerate}

\end{corollary}


\begin{figure}[htb]
 \begin{center} 
  \setlength{\unitlength}{0.8cm}%
    \begin{picture}(5,5)(0,0)
       {\color[rgb]{0,0,0.75}\put(0,0){$\rule{5\unitlength}{5\unitlength}$}
                            \put(0,5){$\overbrace{\rule{5\unitlength}{0\unitlength}}^{n}$}}%
       {\color[rgb]{0.6,1,0.6}\put(0,4){$\rule{1\unitlength}{1\unitlength}$}                  
                           \put(0,4){$\underbrace{\rule{1\unitlength}{0\unitlength}}_{m}$}
                           \multiput(1,3)(1,-1){4}{\rule{1\unitlength}{1\unitlength}}}%
       {\color[rgb]{1,0.25,0.25}\multiput(0,4.9)(0.1,-0.1){50}{\rule{0.1\unitlength}{0.1\unitlength}}}%
   \end{picture}
   \caption{Structure of the covariance matrix.}\label{fig:Covariance}
 \end{center} 
\end{figure}


\begin{corollary} 
\label{corollary:variance}
The variance of the cross-validation estimator is a linear combination of three moments: 
% on the cross-validation errors:
\begin{eqnarray}
 %\Var[\hat{\mu}]
 \theta         &=& \frac{1}{n^2} \sum_{i,j} \Cov(e_i,e_j) \nonumber \\
%                &=& \frac{1}{n^2} (n \sigma^2  + K (m^2 - m) \omega) + (n^2 - K m^2) \gamma \nonumber\\
                &=& \frac{1}{n} \sigma^2 + \frac{m-1}{n} \omega + \frac{n - m}{n} \gamma 
\label{eq:v}
\end{eqnarray}
\end{corollary}

Hence, the problem of estimating $\theta$ % $\Var[\hat{\mu}]$ 
does not involve estimating $n(n+1)/2$ covariances, but it cannot be reduced to that of estimating a single variance parameter. 
Three components intervene, which may be interpreted as follows when $\hat{\mu}$ is the K-fold cross-validation estimate of \EPE:
\begin{enumerate} 
\item
The variance $\sigma^2$ is the average (taken over training sets) variance of errors for ``true'' test examples when algorithm $A$ is fed with training sets of size $m(K-1)$.
\item
The within-block covariance $\omega$ would also apply to ``true'' test examples; 
it arises from the dependence of test errors stemming from the common training set. 
\item
The between-blocks covariance $\gamma$ is due to the dependence of training sets (which share $n(K-2)/K$ examples) and the fact that test block $T_k$ appears in all the training sets $D_\ell$ for $\ell\neq k$. 
% This component is specific to K-fold cross-validation, and has no counterpart in hold-out validation.
\end{enumerate} 
The forthcoming section makes use of this structure to show that there is no universal unbiased estimator of $\theta$. 
% $\Var[\hat{\mu}]$.
 
%------------------------------------------------------------------------------%
\section{No Unbiased Estimator of $\Var[\hat{\mu}]$ Exists}\label{sec:unbiasedness}
%------------------------------------------------------------------------------%

Consider a generic estimator $\hat{\theta}$ that depends on the sequence of cross-validation errors 
$\mathbf{e}=(e_1,e_2,\ldots,e_n)'$.
Let us assume that $\hat{\theta}$ is an analytic function of the errors, so that
we can write its Taylor expansion
\begin{equation}
 \hat{\theta} = \alpha_0 + \sum_i \alpha_{1}(i) e_i + \sum_{i,j} \alpha_{2}(i,j) e_i e_j + \sum_{i,j,k} \alpha_{3}(i,j,k) e_i e_j e_k + \ldots .
\label{eq:theta}
\end{equation}
We first show that for unbiased variance estimates (i.e. $E[\hat{\theta}]=\Var[\hat{\mu}]$),
all the $\alpha_i$ coefficients must vanish except for the second order coefficients $\alpha_{2,i,j}$.

\begin{lemma}
\label{lemma:order2}
There is no universal unbiased estimator of $\Var[\hat{\mu}]$ that involves the $e_i$ in a non-quadratic way.

% All unbiased estimators of $\Var[\hat{\mu}]$ only have second order coefficients in their Taylor expansion.

\noindent
{\bf Proof}

\noindent
Take the expected value of $\hat{\theta}$ expressed as in~(\ref{eq:theta}), and equate it with $\Var[\hat{\mu}]$ (\ref{eq:v}):
\begin{equation*}
\left\{\begin{array}{r@{\,}l}
E[\hat{\theta}] = & \displaystyle{\alpha_0 + \sum_i \alpha_{1}(i) E[e_i] + \sum_{i,j} \alpha_{2}(i,j) E[e_i e_j] + \sum_{i,j,k} \alpha_{3}(i,j,k) E[e_i e_j e_k] + \ldots} \\
\theta          = & \frac{1}{n}\sigma^2 + \frac{m-1}{n} \omega + \frac{n - m}{n} \gamma  .
\end{array}\right.
\end{equation*}
                
For having $E[\hat{\theta}]=\theta$ for all possible values of the moments of $\mathbf{e}$,
one must have $\alpha_0=0$ because $\theta$ has no such constant term,
not depending on any of the moments of $\mathbf{e}$. 
Similarly, $\alpha_{1}(\cdot)$ must be zero because $\theta$ has no term in $E[e_i]=\mu$. 
Finally, the third and higher order coefficients $\alpha_{\ell}(\ldots)$, $\ell>2$ must also be zero because 
$\theta$ has only quantities depending on the second order moments $\sigma^2$, $\omega$ and $\gamma$.

\flushright
Q.E.D.
\end{lemma}

Since estimators that include moments other than the second moments in their expectation
are biased, we now focus on the class of estimators which are quadratic forms of the errors, i.e.
\begin{equation}
 \hat{\theta} = \mathbf{e}' \mathbf{W} \mathbf{e} = \sum_{i,j} W_{ij} e_i e_j .
\label{eq:quadratic}
\end{equation}

\begin{lemma}
\label{lemma:estimators-classes}
The expectation of quadratic estimators $\hat{\theta}$ defined as in (\ref{eq:quadratic}) is a linear combination
of only three terms
%
\begin{equation}\label{E(theta-hat)}
 E[\hat{\theta}] = a (\sigma^2 + \mu^2) + b (\omega + \mu^2) + c (\gamma + \mu^2) ,
\end{equation}
%
where $(a,b,c)$ are defined as follows:
%
\begin{equation*}
\left\{\begin{array}{lcl}
 a &\defequal& \sum_{i=1}^n  W_{ii} , \\
 b &\defequal& \sum_{k=1}^K \sum_{i \in T_k} \sum_{j \in T_k:j\neq i} W_{ij} , \\
 c &\defequal& \sum_{k=1}^K \sum_{\ell\neq k} \sum_{i \in T_k} \sum_{j \in T_\ell} W_{ij} .
\end{array}\right.
\end{equation*}
%

A ``trivial'' representer of estimators with this expected value is
\begin{equation}
 \hat{\theta} = a s_1 + b s_2 + c s_3,
\label{eq:theta-s}
\end{equation}
where $(s_1,s_2,s_3)$ are the only quadratic statistics of $\mathbf{e}$ that are invariants to the within blocks and between blocks permutations described in Lemma~\ref{lemma:distribution}:
%
\begin{equation}
\label{eq:s}
\left\{\begin{array}{lcl}
  s_1 &\defequal& \displaystyle{ \frac{1}{n} \sum_{i=1}^n e_i^2 ,}  \\
  s_2 &\defequal& \displaystyle{ \frac{1}{n(m-1)} \sum_{k=1}^K \sum_{i \in T_k} \sum_{j \in T_k:j\neq i} e_i e_j,}  \\
  s_3 &\defequal& \displaystyle{ \frac{1}{n(n-m)} \sum_{k=1}^K \sum_{\ell\neq k} \sum_{i \in T_k} \sum_{j \in T_\ell} e_i e_j .}
\end{array}\right.
\end{equation}
%
\noindent
{\bf Proof}

\noindent
This result is obtained exploiting Corollary~\ref{corollary:covariances} and grouping
the terms of $\hat{\theta}$ in Equation~(\ref{eq:quadratic}) that have the
same expected values. 
%
\begin{eqnarray}
 E[\hat{\theta}] &=& \sum_{k=1}^K \sum_{i \in T_k} \left( W_{ii} E[e_i^2] + 
 \sum_{j \in T_k:j\neq i} W_{ij} E[e_i e_j]  + 
 \sum_{\ell\neq k} \sum_{j \in T_\ell}  W_{ij} E[e_i e_j] \right) \nonumber \\
 &=& (\sigma^2+\mu^2) \sum_{i=1}^n W_{ii} + 
 (\omega + \mu^2) \sum_{k=1}^K \sum_{i \in T_k} \sum_{j \in T_k:j\neq i}  W_{ij} + \nonumber \\
 & & (\gamma + \mu^2) \sum_{k=1}^K \sum_{\ell\neq k} \sum_{i \in T_k} \sum_{j \in T_\ell} W_{ij} \nonumber \\
 &=& a (\sigma^2+\mu^2) + b (\omega + \mu^2) + c (\gamma + \mu^2) \nonumber \\
 &=& a E[s_1] + b E[s_2] + c E[s_3]  , \nonumber 
\end{eqnarray}
%
which is recognized as the expectation of the estimator defined in Equation~(\ref{eq:theta-s}).

\flushright
Q.E.D.
\end{lemma}

We now use Lemma~\ref{lemma:estimators-classes} to prove that there is no {\em universally} unbiased estimator
of $\Var[\hat{\mu}]$, i.e. there is no estimator $\hat{\theta}$ such that $E[\hat{\theta}] = \Var[\hat{\mu}]$ for all possible distributions of $\mathbf{e}$.

\begin{theorem}
\label{theorem:no-unbiased}

There exists no universally unbiased estimator of $\Var[\hat{\mu}]$.

\noindent
{\bf Proof}

\noindent
Because of Lemma~\ref{lemma:order2} and~\ref{lemma:estimators-classes},
it is enough to prove the result for estimators that are quadratic forms
expressed as in Equation~(\ref{eq:theta-s}). To obtain unbiasedness, the
expected value of that estimator must be equated with $\Var[\hat{\mu}]$
(\ref{eq:v}):
\begin{equation}\label{eq:no-unbiased:proof0}
 a (\sigma^2+\mu^2) + b (\omega + \mu^2) + c (\gamma + \mu^2) = \frac{1}{n} \sigma^2 + \frac{m-1}{n} \omega + \frac{n - m}{n} \gamma .
\end{equation}
For this equality to be satisfied for all distributions of cross-validation errors, it must be satisfied
for all admissible values of $\mu$, $\sigma^2$, $\omega$, and $\gamma$.
This imposes the following unsatisfiable constraints on $(a,b,c)$:
%
\begin{equation}\label{eq:no-unbiased:proof}
 \left\{\begin{array}{lcl}
   a &=& \frac{1}{n} ,   \\
   b &=& \frac{m-1}{n} , \\
   c &=& \frac{n-m}{n} , \\
   a + b + c &=& 0 . 
 \end{array}\right.
\end{equation}
%
\flushright
Q.E.D.
\end{theorem}


%------------------------------------------------------------------------------%
\section{Eigenanalysis of the Covariance Matrix}\label{sec:eigenanalysis}
%------------------------------------------------------------------------------%

One way to gain insight on the origin of the negative statement of Theorem \ref{theorem:no-unbiased} is via the eigenanalysis of $\boldsymbol{\Sigma}$, the covariance of $\mathbf{e}$.
This decomposition can be performed analytically thanks to the very particular block structure displayed in Figure \ref{fig:Covariance}.
% It provides another viewpoint on what can and cannot be estimated unbiasedly in $\boldsymbol{\Sigma}$. 


\begin{lemma}\label{lemma:eigendecomposition}
Let $\mathbf{v}_k$ be the binary vector indicating the membership of each example to test block $k$.
The eigensystem of $\boldsymbol{\Sigma}$ is as follows:
\begin{itemize}
\item 
$\lambda_1 = \sigma^2 - \omega$ with multiplicity $n-K$ and eigenspace defined by the orthogonal of basis $\{\mathbf{v}_k\}_{k=1}^K$;
\item 
$\lambda_2 = \sigma^2 + (m-1) \omega - m \gamma$ with multiplicity $K-1$ and eigenspace defined in the orthogonal of $\mathbf{1}$ by the basis $\{\mathbf{v}_k\}_{k=1}^K$;
\item 
$\lambda_3 =  \sigma^2 + (m-1) \omega + (n-m) \gamma$ with eigenvector $\mathbf{1}$.
\end{itemize}

\noindent
{\bf Proof}

\noindent
From Corollary \ref{corollary:covariances}, the covariance matrix 
$\boldsymbol{\Sigma} = E[\mathbf{e}\mathbf{e}']-E[\mathbf{e}]E[\mathbf{e}]'$ 
can be decomposed as
\[
\boldsymbol{\Sigma} = (\sigma^2-\omega) \boldsymbol{\Sigma}_1 + m (\omega-\gamma) \boldsymbol{\Sigma}_2 + n \gamma \boldsymbol{\Sigma}_3 , 
\]
where $\boldsymbol{\Sigma}_1 = \mathbf{I}$, 
$\boldsymbol{\Sigma}_2 = \frac{1}{m} \left( \mathbf{v}_{1} \ldots \mathbf{v}_{K}\right) 
\left( \mathbf{v}_{1} \ldots \mathbf{v}_{K} \right)'$
and $\boldsymbol{\Sigma}_3 = \frac{1}{n} \mathbf{1} \mathbf{1}'$.

$\boldsymbol{\Sigma}_1$, $\boldsymbol{\Sigma}_2$ and $\boldsymbol{\Sigma}_3$ share the same eigenvectors, with eigenvalues being equal either to zero or one:
\begin{itemize}
\item
the eigenvector $\mathbf{1}$ has eigenvalue $1$ for $\boldsymbol{\Sigma}_1$, $\boldsymbol{\Sigma}_2$ and $\boldsymbol{\Sigma}_3$; 
\item
the eigenspace defined in the orthogonal of $\mathbf{1}$ by the basis $\{\mathbf{v}_k\}_{k=1}^K$ defines $K-1$ eigenvectors with eigenvalues $1$ for $\boldsymbol{\Sigma}_1$ and $\boldsymbol{\Sigma}_2$ and $0$ for $\boldsymbol{\Sigma}_3$; 
\item
all remaining eigenvectors have eigenvalues $1$ for $\boldsymbol{\Sigma}_1$ and $0$ for $\boldsymbol{\Sigma}_2$ and $\boldsymbol{\Sigma}_3$.
\end{itemize}

\flushright
Q.E.D.
\end{lemma}

Lemma \ref{lemma:eigendecomposition} states that the vector
$\mathbf{e}$ can be decomposed into three uncorrelated parts: $n-K$
projections to the subspace orthogonal to $\{\mathbf{v}_k\}_{k=1}^K$,
$K-1$ projections to the subspace spanned by
$\{\mathbf{v}_k\}_{k=1}^K$ in the orthogonal of $\mathbf{1}$, and one
projection on $\mathbf{1}$. A single vector example with $n$
independent elements can be seen as $n$ independent examples.
Similarly, these projections of $\mathbf{e}$ can be equivalently
represented by respectively $n-K$, $K-1$ and one uncorrelated
one-dimensional examples, corresponding to the coordinates of
$\mathbf{e}$ in these subspaces.

In particular, for the projection on $\mathbf{1}$, with only a single one-dimensional point, 
the sample variance is null, resulting in the absence of an unbiased
variance estimator of $\lambda_3$.  The projection of $\mathbf{e}$ on
the eigenvector $\frac{1}{n}\mathbf{1}$ is precisely $\hat{\mu}$.
Hence there is no unbiased estimate of
$\Var[\hat{\mu}]=\frac{\lambda_3}{n}$ when we have only one realization
of the vector $\mathbf{e}$.  For the same reason, even with simple
parametric assumptions on $\mathbf{e}$ (such as $\mathbf{e}$
Gaussian), the maximum likelihood estimate of $\theta$ is not defined.
Only $\lambda_1$ and $\lambda_2$ can be estimated unbiasedly.  Note
that this problem cannot be addressed by performing multiple K-fold
splits of the data set.  Such a procedure would not provide
independent realizations of $\mathbf{e}$.

%------------------------------------------------------------------------------%
\section{Possible Values for $\omega$ and $\gamma$}\label{sec:possiblevalues}
%------------------------------------------------------------------------------%

Theorem \ref{theorem:no-unbiased} states that no estimator is unbiased, 
and in its demonstration, it is shown that the bias of any quadratic estimator is a linear combination of  
$\mu^2$, $\sigma^2$, $\omega$ and $\gamma$.
Regarding estimation, it is thus interesting to see what constraints restrict the possible range of these quantities. 

\begin{lemma}\label{lemma:sog}
For $\hat{\mu}=\CV$ and $\hat{\mu}=\Delta\CV$, the following inequalities hold:
%
\begin{equation*}%\label{eq:sog:proof}
 \begin{array}{ll}
 & \left\{\begin{array}{ccccc}
   0 &\leq& \omega &\leq& \sigma^2 \\
   -\frac{1}{n-m}(\sigma^2+(m-1)\omega)  &\leq& \gamma &\leq& \frac{1}{m} (\sigma^2+(m-1)\omega)  
 \end{array}\right. \\
 \Rightarrow & \left\{\begin{array}{ccccc}
   0 &\leq& \omega &\leq& \sigma^2  \\
   -\frac{m}{n-m}\sigma^2 &\leq& \gamma &\leq& \sigma^2 . 
 \end{array}\right.
 \end{array}
 \end{equation*}
%
The shape of the admissible $(\omega,\gamma)$ region corresponding to the first set of (tighter) inequalities is displayed in Figure \ref{fig:admissible}.


\begin{figure}[htb]
 \begin{center} 
  \begin{tabular}{cc}
    \setlength{\unitlength}{1cm}%
    \begin{picture}(4.2,4.2)(0,0)
      \put(0.75,0.5){\includegraphics[scale = 0.6]{admissibility2.ps}}
     
      \put(2.3,3.9){\makebox(0,0.5)[b]{$K=2$}}
      
      \put(0.9,0){\footnotesize{\makebox(0.3,0.5)[t]{$-\sigma^2$}}}
      \put(2.3,0){\footnotesize{\makebox(0.3,0.5)[t]{$0$}}}
      \put(3.6,0){\footnotesize{\makebox(0.3,0.5)[t]{$\sigma^2$}}}

      \put(2.3,0){$\omega$}

      \put(0.25,0.9){\footnotesize{\makebox(0.5,0.5)[rb]{$-\sigma^2$}}}
      \put(0.25,2.1){\footnotesize{\makebox(0.5,0.5)[rb]{$0$}}}
      \put(0.25,3.4){\footnotesize{\makebox(0.5,0.5)[rb]{$\sigma^2$}}}
      \put(0,2.1){\rotatebox{90}{$\gamma$}}
    \end{picture} &  
    \setlength{\unitlength}{1cm}%
    \begin{picture}(4.2,4.2)(0,0)
      \put(0.75,0.5){\includegraphics[scale = 0.6]{admissibility5.ps}}
     
      \put(2.3,3.9){\makebox(0,0.5)[b]{$K=5$}}
      
      \put(0.9,0){\footnotesize{\makebox(0.3,0.5)[t]{$-\sigma^2$}}}
      \put(2.3,0){\footnotesize{\makebox(0.3,0.5)[t]{$0$}}}
      \put(3.6,0){\footnotesize{\makebox(0.3,0.5)[t]{$\sigma^2$}}}

      \put(2.3,0){$\omega$}

      \put(0.25,0.9){\footnotesize{\makebox(0.5,0.5)[rb]{$-\sigma^2$}}}
      \put(0.25,2.1){\footnotesize{\makebox(0.5,0.5)[rb]{$0$}}}
      \put(0.25,3.4){\footnotesize{\makebox(0.5,0.5)[rb]{$\sigma^2$}}}
      \put(0,2.1){\rotatebox{90}{$\gamma$}}
    \end{picture} \\ 
    \setlength{\unitlength}{1cm}%
    \begin{picture}(4.2,4.2)(0,0)
      \put(0.75,0.5){\includegraphics[scale = 0.6]{admissibility10.ps}}
     
      \put(2.3,3.9){\makebox(0,0.5)[b]{$K=10$}}
      
      \put(0.9,0){\footnotesize{\makebox(0.3,0.5)[t]{$-\sigma^2$}}}
      \put(2.3,0){\footnotesize{\makebox(0.3,0.5)[t]{$0$}}}
      \put(3.6,0){\footnotesize{\makebox(0.3,0.5)[t]{$\sigma^2$}}}

      \put(2.3,0){$\omega$}
      \put(0.25,0.9){\footnotesize{\makebox(0.5,0.5)[rb]{$-\sigma^2$}}}
      \put(0.25,2.1){\footnotesize{\makebox(0.5,0.5)[rb]{$0$}}}
      \put(0.25,3.4){\footnotesize{\makebox(0.5,0.5)[rb]{$\sigma^2$}}}
      \put(0,2.1){\rotatebox{90}{$\gamma$}}
    \end{picture} &  
    \setlength{\unitlength}{1cm}%
    \begin{picture}(4.2,4.2)(0,0)
      \put(0.75,0.5){\includegraphics[scale = 0.6]{admissibility100.ps}}
     
      \put(2.3,3.9){\makebox(0,0.5)[b]{$K=100$}}
      
      \put(0.9,0){\footnotesize{\makebox(0.3,0.5)[t]{$-\sigma^2$}}}
      \put(2.3,0){\footnotesize{\makebox(0.3,0.5)[t]{$0$}}}
      \put(3.6,0){\footnotesize{\makebox(0.3,0.5)[t]{$\sigma^2$}}}

      \put(2.3,0){$\omega$}

      \put(0.25,0.9){\footnotesize{\makebox(0.5,0.5)[rb]{$-\sigma^2$}}}
      \put(0.25,2.1){\footnotesize{\makebox(0.5,0.5)[rb]{$0$}}}
      \put(0.25,3.4){\footnotesize{\makebox(0.5,0.5)[rb]{$\sigma^2$}}}
      \put(0,2.1){\rotatebox{90}{$\gamma$}}
    \end{picture} 
   \end{tabular}
   \caption{Possible values of $(\omega,\gamma)$ according to $\sigma^2$ for  $n=200$ and $K=\{2,5,10,100\}$.}\label{fig:admissible}
 \end{center} 
\end{figure}
\noindent
{\bf Proof}

\noindent
The constraints on $\omega$ result from the Cauchy-Schwartz inequality which provides $\Cov(u,v)^2 \leq \Var[u]\Var[v]$, hence
%
\begin{equation*}
   -\sigma^2 \leq \omega \leq \sigma^2 .
\end{equation*}
%
Moreover, the following reasoning shows that, 
for $\hat{\mu}=\CV$ and $\hat{\mu}=\Delta\CV$, 
$\omega$ is non-negative:
$\omega$ is the covariance of (differences in) test errors for training sets of size $n-m$ and test sets of size $\ell=m$.
The variance of the average test error is given by the mean of covariances $\frac{1}{\ell}(\sigma^2+(\ell-1)\omega)$.
The variance $\sigma^2$ and covariance $\omega$ of test errors are not affected by $\ell$, 
and the variance of the average test error should be non-negative for any test set size $\ell$. 
Hence $\omega$ is bound to be non-negative.
When this type of reasoning cannot be used, as for $\hat{\mu}=\JK$, 
$\omega$ can only be proved to be greater than $-\sigma^2/(m-1)$.

The constraints on $\gamma$ simply rephrase that the eigenvalues $\lambda_2$ and $\lambda_3$ of the covariance matrix $\boldsymbol{\Sigma}$ should be non-negative.
The simpler (and looser) form is obtained by using $\omega \leq \sigma^2$.
\flushright
Q.E.D.
\end{lemma}

The admissible $(\omega,\gamma)$ region obtained in Lemma \ref{lemma:sog} is very large.
Furthermore, there is no constraint linking $\mu$ and $\sigma^2$, the mean and variance of $e_i$. 
Hence we cannot propose a variance estimate with universally small bias.

%------------------------------------------------------------------------------%
\section{Experiments}
\label{sec:experiments}
%------------------------------------------------------------------------------%

We already mentioned that the bias of any quadratic estimator is a linear combination of  
$\mu^2$, $\sigma^2$, $\omega$ and $\gamma$.
The admissible values provided in the preceding section suggest that $\omega$ and $\gamma$ cannot be proved to be negligible compared to $\sigma^2$.
This section illustrates that in practice, 
the contribution to the variance of $\hat{\mu}$ due to $\omega$ and $\gamma$ (see Equation (\ref{eq:v}))
can be of same order than the one due $\sigma^2$.
It therefore suggests that the estimators of $\theta$ should indeed take into account the correlations of $e_i$.

\begin{experiment}\label{exp:truevar}
True variance of K-fold cross-validation.

We repeat the experimental setup of Experiment \ref{exp:hold-out}, except that now, 
we are in the more realistic situation where only one sample of size $n$ is available.
Since cross-validation is known to be sensitive to the instability of algorithms, 
in addition to this standard setup, we also consider another one with outliers:

The input $\mathbf{x_i} = (x_{i1},\ldots,x_{id})'$ is still $30$-dimensional, 
but it is now a mixture of two centered Gaussian variables: 
let $t_i$ be a binary variable, with $P(t_i=1)=p=0.95$;
when $t_i=1$, $x_i\sim{\cal N}(0,\mathbf{I})$;
when $t_i=0$, $x_i\sim{\cal N}(0,100\mathbf{I})$;
$y_i = \sqrt{3/(d(p+100(1-p)))}\sum_{k=1}^d x_{ik} + \varepsilon_i$ 
with $\varepsilon_i\sim{\cal N}(0,1/(p+100(1-p)))$   when $t_i=1$
and  $\varepsilon_i\sim{\cal N}(0,100/(p+100(1-p)))$ when $t_i=0$.
\end{experiment}

We now look at the variance of K-fold cross-validation ($K=10$), 
and decompose in the three orthogonal components $\sigma^2$, $\omega$ and $\gamma$.
The results are shown in Figure \ref{fig:variance}.
% Note that in these experiments, $\gamma$ was always estimated positive.


\begin{figure}[htb]
 \begin{center} 
  \begin{tabular}{ccc}
   \includegraphics[scale = 0.6]{LinTrueRepK10noOut.ps} &  
   \includegraphics[scale = 0.6]{LinTrueRepK10badOut.ps} \\
   no outliers &  outliers
   \end{tabular}
   \caption{Bar plots of the contributions to total variance $\Var[\CV]$
   due to $\sigma^2$, $\omega$ and $\gamma$
   {\em vs.} the number of training examples $n-m$ for Experiment \ref{exp:truevar}.}\label{fig:variance}
 \end{center} 
\end{figure}

When there are no outliers, the contribution of $\gamma$ is very important for small sample sizes. 
For large sample sizes, the overall variance is considerably reduced and is mainly caused by $\sigma^2$.
In these situations, the learning algorithm returns very similar answers for all training sets.
When there are outliers, $\omega$ has little effect, but the contribution of $\gamma$ is of same order as the one of $\sigma^2$, even when the ratio of examples to free parameters is large (here up to 20).
Thus, in difficult situations, where $A(D)$ varies according to the realization of $D$,
neglecting the effect of $\omega$ and $\gamma$ can be expected to introduce a bias of the order of the true variance.
  
It is also interesting to see how these quantities are affected by the number of folds $K$.
The decomposition of $\theta$ in $\sigma^2$, $\omega$ and $\gamma$ (\ref{eq:v}) does not imply that $K$ should be set either to $n$ or to $2$ (according to the sign of $\omega-\gamma$) in order to minimize the variance of $\hat{\mu}$.
Modifying $K$ affects $\sigma^2$, $\omega$ and $\gamma$ through the size and overlaps of the training sets $D_1,\ldots,D_K$, as illustrated in Figure \ref{fig:variance2}.
For a fixed sample size, the variance of $\hat{\mu}$ and the contribution of $\sigma^2$, $\omega$ and $\gamma$ effects varies smoothly with $K$.\footnote{
Of course, the mean of $\hat{\mu}$ is also affected in the process.
}
The experiments with and without outliers illustrate that there is no general trend either in variance or decomposition of the variance in its $\sigma^2$, $\omega$ and $\gamma$ components. 
The minimum variance can be reached for $K=n$ or for an intermediate value of $K$.

\begin{figure}[htb]
 \begin{center} 
  \begin{tabular}{ccc}
   \includegraphics[scale = 0.6]{LinTrueRepN120noOut.ps} &  
   \includegraphics[scale = 0.6]{LinTrueRepN120badOut.ps} \\
   no outliers &  outliers
   \end{tabular}
   \caption{Bar plots of contributions of $\sigma^2$, $\omega$ and $\gamma$ to $\theta$
   {\em vs.} $K$ for $n=120$ for Experiment \ref{exp:truevar}.}\label{fig:variance2}
 \end{center} 
\end{figure}

We also report an experiment illustrating that the previous observations also apply to classification on real data.
The variance of K-fold cross-validation ($K=10$), 
decomposed in the three orthogonal components $\sigma^2$, $\omega$ and $\gamma$ is displayed in Figure \ref{fig:letter}.

\begin{experiment}\label{exp:letter}
Classification with trees on the Letter data set.

The Letter data set comprises 20,000 examples described by  16 numeric features. 
The original setup considers 26 categories reprensenting the letters of the roman alphabet. 
Here, we used a simplified setup with 2 classes (A to M) {\em vs.} (N to Z) in order to obtain sensible results for small sample sizes.

Accurate estimates of $\sigma^2$, $\omega$ and $\gamma$ require many independent training samples.
This was achieved by considering the set of 20,000 examples to be the population, from which independent training samples were drawn by uniform sampling with replacement.
\end{experiment}

% Note that in these experiments, $\gamma$ was always estimated positive.


\begin{figure}[htb]
 \begin{center} 
   \includegraphics[scale = 0.6]{LetterK10.ps}  
   \caption{Bar plots of the contributions to total variance $\Var[\CV]$
   due to $\sigma^2$, $\omega$ and $\gamma$
   {\em vs.} the number of training examples $n$  for Experiment \ref{exp:letter}.}\label{fig:letter}
 \end{center} 
\end{figure}

Here again, the variance of \CV\ is mainly due to $\sigma^2$ and $\gamma$. 
According to the number of training examples, $\sigma^2$ is only responsible for 50 to 70\,\% of the total variance,
so that a variance estimate based solely on $\sigma^2$ has a negative bias of the order of magnitude of the variance itself.
 
%------------------------------------------------------------------------------%
\section{Special Cases}
%------------------------------------------------------------------------------%
This section addresses how our main result can be transposed to
hold-out estimates of generalization error.  We also detail how it
applies to two specific instances of the general K-fold
cross-validation scheme: two-fold and leave-one-out cross validation.
\subsection{Hold-Out Estimate of \EPE}

When having $K$ {\em independent training and test sets}, the
structure of hold-out errors resemble the one of cross-validation
errors, except that we know (from the independence of training and
test sets) that $\gamma =0$.  This knowledge allows to build the
unbiased variance estimate $\widehat{\theta_2}$ described in
\ref{sec:framework:holdout}.  This can be seen directly in the proof
of Theorem \ref{theorem:no-unbiased}: knowing that $\gamma =0$ removes
the third equation in the linear system (\ref{eq:no-unbiased:proof}).
In practice, one is often restricted to $K=1$ (ordinary hold-out
test), which allows to estimate the variance due to the finite test
set but not due to the particular choice of training set.

\subsection{Two-Fold Cross-Validation}

Two-fold cross-validation has been advocated to perform hypothesis testing \citep{Dietterich98,Alpaydin99}.
It is a special case of K-fold cross-validation since the training blocks are mutually independent since they do not overlap.
However, this independence does not modify the structure of $\mathbf{e}$ in the sense that $\gamma$ is not null.
The between-block correlation stems from the fact that the training block $D_{1}$ is the test block $T_{2}$ and vice-versa.

\subsection{Leave-One-Out Cross-Validation}

Leave-one-out cross-validation is a particular case of K-fold cross-validation, where $K=n$.
The structure of the covariance matrix is simplified, without diagonal blocks: 
$\boldsymbol{\Sigma} = (\sigma^2-\gamma) \boldsymbol{\Sigma}_1 + n \gamma \boldsymbol{\Sigma}_3$.
The estimation difficulties however remain: even in this particular case, there is no unbiased estimate of variance.
From the definition of $b$ (Lemma~\ref{lemma:estimators-classes}), we have $b=0$, and with $m=1$ the linear system~(\ref{eq:no-unbiased:proof}) reads
\[
 \left\{\begin{array}{lcl}
   a &=& \frac{1}{n} ,   \\
   c &=& \frac{n-1}{n} , \\
   a +  c &=& 0 , 
 \end{array}\right.
\]
which still admits no solution.

%------------------------------------------------------------------------------%
\section{Conclusions}\label{sec:conclusion}
%------------------------------------------------------------------------------%

It is known that K-fold cross-validation may suffer from high variability, 
which can be responsible for bad choices in model selection and erratic behavior in the estimated expected prediction error.

In this paper, 
we show that estimating the variance of K-fold cross-validation is difficult.
Estimating a variance can be done from independent realizations or from dependent realizations whose correlation is known.
K-fold cross-validation produces dependent test errors. 
Our analysis shows that although the correlations are structured in a very simple manner,
their values cannot be estimated unbiasedly.
Consequently, there is no unbiased estimator of the variance of K-fold cross-validation.

Our experimental section shows that in very simple cases, 
the bias incurred by ignoring the dependencies between test errors will be of the order of the variance itself.
These experiments illustrate thus that the assessment of the significance of observed differences in cross-validation scores should be treated with much caution.
The problem being unveiled, the next step of this study consists in building and comparing variance estimators dedicated to the very specific structure of the test error dependencies.
  
%------------------------------------------------------------------------------%
\bibliography{var-kfoldxv}
%------------------------------------------------------------------------------%

% \appendix


\end{document}
