\documentclass{article} % For LaTeX2e
\usepackage{nips11submit_e,times}
%\documentstyle[nips10submit_09,times,art10]{article} % For LaTeX 2.09

\usepackage{graphicx} % more modern
\usepackage{epstopdf}

% For algorithms
\usepackage{algorithm}
\usepackage{algorithmic}

%my settings
\usepackage{amsmath}
\usepackage{color}
\usepackage{natbib}
\usepackage{multirow}
\newcommand{\dist}{\mathrm{dist}}
\newcommand{\sump}{\mathrm{sp}}
\newcommand{\ssp}{\mathrm{ssp}}
\newcommand{\vecx}{{\mathbf{x}}}
\newcommand{\vecy}{{\mathbf{y}}}
\newcommand{\veca}{{\mathbf{a}}}
\newcommand{\vecb}{{\mathbf{b}}}
\newcommand{\vecq}{{\mathbf{q}}}
\newcommand{\nvecx}{{\hat{\vecx}}}
\newcommand{\nvecy}{{\hat{\vecy}}}
\newcommand{\nvecq}{{\hat{\vecq}}}
\newcommand{\na}{{\hat{a}}}
\newcommand{\nb}{{\hat{b}}}
\newcommand{\vareps}{{\varepsilon}}
\newcommand{\msl}{{\text{\sc msl}}}
\newcommand{\msu}{{\text{\sc msu}}}
% For algorithms
\usepackage{algorithm}
\usepackage{algorithmic}

\newtheorem{thm}{Theorem}
\newtheorem{lem}{Lemma}
\newtheorem{cor}{Corollary}

\newcommand{\complain}[1]{\textcolor{red}{#1}}

\title{Convergent % Incremental %Upper and Lower 
Bounds on the Euclidean Distance}

\author{Yoonho Hwang \hspace{1cm} Hee-Kap Ahn\\
Department of Computer Science and Engineering\\
Pohang University of Science and Technology\\
POSTECH, Pohang, Gyungbuk, Korea(ROK)\\
\texttt{\{cypher,heekap\}@postech.ac.kr}\\
}

% The \author macro works with any number of authors. There are two commands
% used to separate the names and addresses of multiple authors: \And and \AND.
%
% Using \And between authors leaves it to \LaTeX{} to determine where to break
% the lines. Using \AND forces a linebreak at that point. So, if \LaTeX{}
% puts 3 of 4 authors names on the first line, and the last on the second
% line, try using \AND instead of \And before the third author name.

\newcommand{\fix}{\marginpar{FIX}}
\newcommand{\new}{\marginpar{NEW}}


\nipsfinalcopy % Uncomment for camera-ready version

\begin{document}

\maketitle

\begin{abstract}
Given a set $V$ of $n$ vectors in $d$-dimensional space, we provide
an efficient method for computing quality upper and lower bounds of
the Euclidean distances between a pair of vectors in $V$.
For this purpose, we define a distance measure, called
the {\em MS-distance}, by using the mean and the standard
deviation values of vectors in $V$.
% We call it MS(Mean-Standard deviation) distance.
Once we compute the mean and the standard
deviation values of vectors in $V$ in $O(dn)$ time,
the MS-distance provides upper and lower bounds
of Euclidean distance between any pair of vectors in $V$
in constant time.
Furthermore, these bounds can be refined further in such a way 
to converge monotonically to the exact Euclidean distance within 
$d$ {\em refinement steps}.
An analysis on a random sequence of refinement steps shows that
the MS-distance provides very tight bounds in only a few refinement
steps.
% of a typical sequence.
The MS-distance can be used to various applications where the Euclidean
distance is used to measure the proximity or similarity between objects.
%various proximity searches. 
We provide experimental results on
the nearest and the farthest neighbor searches.
\end{abstract}

\section{Introduction}
\label{sec:intro}
The Euclidean distance between two vectors $\vecx$ and $\vecy$ in 
$d$-dimensional
space is a typical distance measure that reflects their proximity in the 
space. Measuring the Euclidean distance is a fundamental operation in 
computer science, including the areas 
of database, computational geometry, computer vision and computer 
graphics. In machine learning, 
the Euclidean distance, denoted by $\dist(\vecx,\vecy)$, or 
it's variations(for example, $e^{||\vecx-\vecy||}$) are widely 
used to measure data similarity for clustering~\cite{kmeans}, 
classification~\cite{svm} and so on.

A typical problem is as follows. Given two sets $X$ and $Y$ of 
vectors in $d$-dimensional space, our goal is to find a pair $(\vecx, \vecy)$,
for $\vecx\in X$ and $\vecy\in Y$,
such that $\dist(\vecx,\vecy)$ is the optimum (minimum or maximum) over all 
such pairs. For the nearest or farthest neighbor searches, $X$ is the set
consisting of a single query point while $Y$ consists of all candidate 
data points. % that are candidates of the nearest (or the farthest) point 
% from the query.
If the dimension is low, a brute-force computation would be fast enough.
However, data sets in areas such as optimization, computer vision, machine
learning or statistics often live in spaces of dimensionality in the
order of thousands or millions.
In $d$ dimensional space, a single distance computation already takes $O(d)$ 
time, thus the cost for finding the nearest or farthest neighbor becomes
$O(dnm)$ time, where $n$ and $m$ are the cardinalities of $X$ and $Y$, 
respectively.

Several techniques have been proposed to reduce computation cost 
for computing distance. Probably PCA (principal component analysis) 
is the most frequently used technique
for this purpose~\cite{pca}, in which we use an orthogonal 
transformation based on PCA to convert a set of given data so that 
the dimensionality of the transformed data is reduced. Then 
it computes distances between pairs of transformed data efficiently.
% to use a dimensionality 
%reduction technique so that 
However, this transformation does not preserve the pairwise distances 
of data in general, therefore there is no guarantee on the computation 
results.

If we restrict ourselves to the nearest neighbor search,
some methods using space partitioning trees such as KD-tree~\cite{kdtree}, 
R-tree~\cite{rtree}, or their variations  have been widely used.
However, they become impractical for high dimensions because of their 
poor performance in constructing data structures for queries.
Recently, cover tree~\cite{covertree} has been used for high dimensional
nearest neighbor search, but its construction time increases drastically
as the dimension increases~\cite{intrinsic}.
% and they do not seem 
% useful for related application problems other than the nearest neighbor 
% search.

Another approach that has attracted some attention is to compute a good 
bound of the exact Euclidean distance efficiently such that it
can be used to filter off some unnecessary computation, for example,
the distance computation between two vectors that are far apart from each 
other in nearest neighbor search. 
One of such methods is to compute a distance bound using the inner 
product approximation~\cite{inner}. This method, however, requires 
the distribution of the input data to be known in advance, and works 
only on data in some predetermined distribution. Another method is
to compute a distance bound using bitwise operations~\cite{vafile}. 
% also present a 
But this method works well only on uniformly distributed vectors, 
and requires $O(2^d)$ bitwise operations in $d$ dimension. A method
using an index structure~\cite{idist} provides an effective filtering method
based on the triangle inequality. But this works well only when data 
are well clustered.

In this paper, we
define a distance measure, called
the {\em MS-distance}, by using the mean and the standard
deviation values of vectors in $V$.
% We call it MS(Mean-Standard deviation) distance.
Once we compute the mean and the standard
deviation values of vectors in $V$ in $O(dn)$ time,
the MS-distance provides tight upper and lower bounds
of Euclidean distance between any pair of vectors in $V$
in constant time. 
Furthermore, these bounds can be refined further in such a way to
converge monotonically to the exact Euclidean distance within 
$d$ refinement steps. Each refinement step takes constant time.

We provide an analysis on a random sequence
of $k$ refinement steps for $0\leq k\leq d$, which shows a
good expectation on the lower and
upper bounds. This can justify
that the MS-distance provides very tight bounds in 
a few refinement steps of a typical sequence.
We also show that the MS-distance can be used in 
fast filtering.
Note that we do not use any assumption on data distribution.

The MS-distance can be used to various applications where the Euclidean
distance is a measure for proximity or similarity between objects.
Among them, we provide experimental results on
the nearest and the farthest neighbor searches.

\section{An Upper and A Lower Bounds of the Euclidean Distance}
For a $d$-dimensional vector $\vecx=[x_1,x_2,\ldots,x_d]$, we denote
its mean by $\mu_\vecx = \frac{1}{d}\sum_{i=1}^d x_i$ and its variance by
$\sigma_\vecx^2 = \frac{1}{d}\sum_{i=1}^d (x_i-\mu_\vecx)^2$.  For a pair
of vectors $\vecx$ and $\vecy$, we can reformulate the squared
Euclidean distance between $\vecx$ and $\vecy$ as follows. 
Let $\veca=[a_1,a_2,\ldots,a_d]$ and $\vecb=[b_1,b_2,\ldots,b_d]$ such that $a_i = x_i-\mu_\vecx$ and $b_i = y_i-\mu_\vecy$.
% For a query vector $\vecy\notin V$, we similar notation to $\vecx$,
% and $\mu_\vecy$,$\sigma_\vecy$ are computed once in runtime.  For a
% little trick on the formulation, we define $a_i = x_i-\mu_\vecx$,
% $b_i = y_i-\mu_\vecy$.  Then the squared Eculidean distance between
% $\vecx$ and $\vecy$ formulated by setting $a_i = x_i-\mu_\vecx$ and
% $b_i = y_i-\mu_\vecy$ as follows.
\begin{align}
  \dist(\vecx,\vecy)^2 &= \sum_{i=1}^d(x_i-y_i)^2 \nonumber\\
  &= \sum_{i=1}^d((\mu_\vecx+a_i)-(\mu_\vecy+b_i))^2\nonumber\\
  &= \sum_{i=1}^d(\mu_\vecx^2+2a_i\mu_\vecx+a_i^2
  +\mu_\vecy^2+2b_i\mu_\vecy+b_i^2
  -2(\mu_\vecx\mu_\vecy+a_i\mu_\vecy+b_i\mu_\vecx+a_ib_i))
  \label{eqn:dist.expand}\\
  &= \sum_{i=1}^d(\mu_\vecx^2 -2\mu_\vecx\mu_\vecy +\mu_\vecy^2 +a_i^2
  +b_i^2 -2a_ib_i)
  \label{eqn:dist.middle}\\
  &= d\left((\mu_\vecx-\mu_\vecy)^2+(\sigma_\vecx+\sigma_\vecy)^2\right) 
    - 2d\sigma_\vecx\sigma_\vecy - 2\sum_{i=1}^d a_ib_i
  \label{eqn:dist.upper.pre}\\
  &= d\left((\mu_\vecx-\mu_\vecy)^2+(\sigma_\vecx-\sigma_\vecy)^2\right) 
    + 2d\sigma_\vecx\sigma_\vecy - 2\sum_{i=1}^d a_ib_i.\label{eqn:dist.lower.pre}
\end{align}
By the definitions of $a_i$ and $b_i$, we have $\sum_{i=1}^d a_i =\sum_{i=1}^d b_i=0$, and $\frac{1}{d}\sum_{i=1}^d a_i^2 = \sigma_\vecx^2$. 
By the first properties, equation (\ref{eqn:dist.expand}) is 
simplified to (\ref{eqn:dist.middle}), 
and by the second property, equations (\ref{eqn:dist.middle}) becomes 
(\ref{eqn:dist.upper.pre}) and (\ref{eqn:dist.lower.pre}). 

Note that equations (\ref{eqn:dist.upper.pre}) and 
(\ref{eqn:dist.lower.pre}) % of the squared Euclidean distance 
are composed of the mean and variance values 
(their products and squared values, multiplied by $d$) 
of $\vecx$ and $\vecy$, except the last summations. 
Thus, once we preprocess $V$ of $n$ vectors 
such that both $\mu_\vecx$ and $\sigma_\vecx$ for all $\vecx\in V$ 
are computed in $O(dn)$ time and stored in a table of size $O(n)$,
this sum can be computed in constant time for any pair of vectors, regardless
of the dimension.

% This is great advantage, because we can compute much non-negative term $d((\mu_\vecx-\mu_\vecy)^2+\sigma_\vecx^2+\sigma_\vecy^2)$ in constant time, regardless of it's dimension.
The last summation, $\sum_i^d a_ib_i$, is the inner product $\langle\veca,\vecb\rangle$, 
and therefore by applying the Cauchy-Schwarz inequality we get
\begin{align}
|\langle\veca,\vecb\rangle|=|\sum_{i=1}^d a_ib_i| \le \sqrt{(\sum_{i=1}^d a_i^2)(\sum_{i=1}^d b_i^2)} = d \sigma_\vecx \sigma_\vecy. \label{eqn:cauchy}
\end{align}
This gives us the following upper and lower bounds of the squared Euclidean 
distance from equations (\ref{eqn:dist.upper.pre}) and 
(\ref{eqn:dist.lower.pre}).
\begin{lem}
For two $d$-dimensional vectors $\vecx,\vecy$, the followings hold.
\begin{align}
\dist(\vecx,\vecy)^2\ge d\left((\mu_\vecx-\mu_\vecy)^2+(\sigma_\vecx-\sigma_\vecy)^2\right)\label{eqn:dist.lower}\\
\dist(\vecx,\vecy)^2\le d\left((\mu_\vecx-\mu_\vecy)^2+(\sigma_\vecx+\sigma_\vecy)^2\right)\label{eqn:dist.upper}
\end{align}
\end{lem}
% The equality in the Cauchy-Schwarz inequality holds if and only if two 
% vectors are linearly dependent. In our case, equalities in (~\ref{eqn:dist.bound}) hold if and only if $a$=$b$ or $a$=$-b$,
% and these are simply realizable.
% Therefore, for any arbitrarily small $\delta>0$, if we add $\delta$ ($-\delta$ for upper bound) or multiply $1+\delta$($1-\delta$) on the lower bound, we loose correctness because of $a$=$b$ ($a$=$-b$) case.
% In this sence, the bounds in (~\ref{eqn:dist.bound}) are \emph{tight}.

% Note that these bounds (~\ref{eqn:dist.bound}) are consist of constant number of terms.
% Therefore, we can compute the upper and lower bound of distance in constant time, regardless of vector's dimension.

\section{The MS-distance}
\label{sec:ms-distance}
The lower and upper bounds in inequalities (\ref{eqn:dist.lower}) and (\ref{eqn:dist.upper}) can be computed in constant time once we compute
the mean and standard variance values of each vector in $V$ in the
preprocessing.
However, in some applications these bounds may not be tight enough.
In this section, we introduce the {\em MS-distance} which not only 
provides lower and upper bounds of the Euclidean distance 
in constant time, but also could be refined further in such a way
to converge to the exact Euclidean distance within $d$ steps.
% Furthermore, if we refine this bound with $O(d)$ more computation,
% it becomes identical to the Euclidean distance.\\

% Then we have 
% \[
% \dist^2(\vecx,\vecy)=d((\mu_\vecx-\mu_\vecy)^2+\sigma_\vecx^2+\sigma_\vecy^2)-2\langle \veca,\vecb\rangle.
% \]
To do this, we reformulate the last term of equations (\ref{eqn:dist.upper.pre}) 
and (\ref{eqn:dist.lower.pre}),
that is, the inner product $\langle\veca,\vecb\rangle$.
If the norms $||\veca||=\sqrt{\sum_{i=1}^da_i^2}$ or $||\vecb||=\sqrt{\sum_{i=1}^db_i^2}$ are zero, then $\sum_{i=1}^d a_ib_i=0$, thus the upper and lower 
bounds become the same. This implies that we can compute the exact Euclidean
distance in constant time. So from now on, we assume that both $||\veca||$ and 
$||\vecb||$ are non-zero. We reformulate the inner product $\langle\veca,\vecb\rangle$. 
% Let $\nveca=[\na_1,\ldots,\na_d]$ and $\nvecb=[\nb_1,\ldots,\nb_d]$ be the normalized vectors of $\veca$ and $\vecb$, that is, 
% $\na_i=a_i/||\veca||$ and $\nb_i=b_i/||\vecb||$ for all $1\leq i\leq d$.
% Because $||\veca||=\sqrt{\sum_{i=1}^d a_i^2}= \sigma_\vecx\sqrt{d}$,
% % by definition of $a_i$, and $||a||$,$\sigma_\vecx$,$d$ are all positive, 
% we can rewrite $\na_i=a_i/(\sigma_\vecx\sqrt{d})$ and 
% $\nb_i=b_i/(\sigma_\vecy\sqrt{d})$.
% % If we apply the Cauchy-Schwarz inequality~\ref{eqn:cauchy} again on 
% % $\langle\nveca,\nvecb\rangle$,
% % equality holds if and only if they are linearly dependent. 
% % Because of 
% % normalization, $\langle\nveca,\nvecb\rangle$ is maximized when 
% %$\nveca=\nvecb$, and becomes strictly smaller when $\nveca\neq\nvecb$.
% % To see how much smaller $\langle\nveca,\nvecb\rangle$ is 
% % when $\nveca\neq\nvecb$, we define the difference $\vareps_i=\nb_i-\na_i$.
% Since $\sum_{i=1}^d \na_i^2 = \sum_{i=1}^d \nb_i^2 = 1$, the following holds.
% \begin{align}
% \langle\nveca,\nvecb\rangle=\sum_{i=1}^d \na_i\nb_i = -\frac{1}{2}(\sum_{i=1}^d (\nb_i-\na_i)^2-\sum_{i=1}^d\nb_i^2-\sum_{i=1}^d \na_i^2) =1-\frac{1}{2}\sum_{i=1}^d (\nb_i-\na_i)^2.
% \label{eqn:normal.ab}
% \end{align}
%\sum_i^d \nb_i^2 &= \sum_i^d (\na_i+\vareps_i)^2 = \sum_i^d (\na_i^2 + 2\na_i\vareps_i +\vareps_i^2) = 1 \label{eqn:normal.b}\\
%\end{align}
% By using the above property, the inner product $\langle\nveca,\nvecb\rangle$
% can be formulated as follows.

% \begin{align}
% \sum_{i=1}^d \na_i\nb_i &= \sum_{i=1}^d \na_i(\na_i+\vareps_i)\nonumber\\
% &=1+\sum_{i=1}^d \na_i\vareps_i\nonumber\\
% &= 1-\frac{1}{2} \sum_{i=1}^d \vareps_i^2\nonumber\\
% &= 1-\frac{1}{2} \sum_{i=1}^d (\nb_i - \na_i)^2 \label{eqn:nzero}
% \end{align}

\begin{align}
\sum_{i=1}^d a_ib_i
&= d\sigma_\vecx\sigma_\vecy - d\sigma_\vecx\sigma_\vecy +\sum_{i=1}^d a_ib_i\nonumber\\
&= d\sigma_\vecx\sigma_\vecy - \frac{\sigma_\vecx\sigma_\vecy}{2}\left(2d-\sum_{i=1}^d\frac{2a_ib_i}{\sigma_\vecx\sigma_\vecy}\right)\nonumber\\
&= d\sigma_\vecx\sigma_\vecy - \frac{\sigma_\vecx\sigma_\vecy}{2}\left(\sum_{i=1}^d\left(\frac{a_i}{\sigma_\vecx}\right)^2+\sum_{i=1}^d\left(\frac{b_i}{\sigma_\vecy}\right)^2-\sum_{i=1}^d\frac{2a_ib_i}{\sigma_\vecx\sigma_\vecy}\right)\label{eqn:inner.squared}\\
% &=d\sigma_\vecx\sigma_\vecy (\sum_{i=1}^d \na_i\nb_i)\nonumber\\
% &=d\sigma_\vecx\sigma_\vecy \left(1-\frac{1}{2}\sum_{i=1}^d (\nb_i - \na_i)^2\right)\nonumber\\
% &=d\sigma_\vecx\sigma_\vecy \left(1-\frac{1}{2}\sum_{i=1}^d \frac{1}{d}(\frac{b_i}{\sigma_\vecy}-\frac{a_i}{\sigma_\vecx})^2\right)\nonumber\\
&=d\sigma_\vecx\sigma_\vecy - \frac{\sigma_\vecx\sigma_\vecy}{2}\sum_{i=1}^d {(\frac{b_i}{\sigma_\vecy}-\frac{a_i}{\sigma_\vecx})^2}\label{eqn:inner.upper}\\
&=-d\sigma_\vecx\sigma_\vecy + \frac{\sigma_\vecx\sigma_\vecy}{2}\sum_{i=1}^d {(\frac{b_i}{\sigma_\vecy}+\frac{a_i}{\sigma_\vecx})^2}\label{eqn:inner.lower}
\end{align}
Equation (\ref{eqn:inner.squared}) is because of 
$\sum_{i=1}^da_i^2=d\sigma_\vecx^2$ and $\sum_{i=1}^db_i^2=d\sigma_\vecy^2$.
We can also get equation (\ref{eqn:inner.lower}) by switching the roles of
the term $-d\sigma_\vecx\sigma_\vecy$ and the term $d\sigma_\vecx\sigma_\vecy$ in the above 
equations. 

\paragraph{Definition.} Now we define the MS-distance between $\vecx$ and $\vecy$ 
in its lower bound form, denoted by $\msl(\vecx,\vecy,k)$,
by replacing the last term of equation (\ref{eqn:dist.upper.pre}) with equation 
(\ref{eqn:inner.upper}), and in its upper bound form, denoted by
$\msu(\vecx,\vecy,k)$
by replacing the last term of equation (\ref{eqn:dist.lower.pre}) with equation 
(\ref{eqn:inner.lower}).
The MS-distance makes use of the nonincreasing intermediate values for
its upper bound and the nondecreasing intermediate values for its lower
bound. We let $a_0=b_0=0$.
\begin{align}
\msl(\vecx,\vecy,k)=d\left((\mu_\vecx-\mu_\vecy)^2+(\sigma_\vecx-\sigma_\vecy)^2\right) + \sigma_\vecx\sigma_\vecy\sum_{i=0}^k {\left(\frac{b_i}{\sigma_\vecy}-\frac{a_i}{\sigma_\vecx}\right)^2} \label{eqn:ms.lower}\\
\msu(\vecx,\vecy,k)=d\left((\mu_\vecx-\mu_\vecy)^2+(\sigma_\vecx+\sigma_\vecy)^2\right) - \sigma_\vecx\sigma_\vecy\sum_{i=0}^k {\left(\frac{b_i}{\sigma_\vecy}+\frac{a_i}{\sigma_\vecx}\right)^2} \label{eqn:ms.upper}
\end{align}

\paragraph{Properties.}
Note that equation (\ref{eqn:ms.lower}) is nondecreasing
and equation (\ref{eqn:ms.upper}) is nonincreasing while $i$ 
increases from $0$ to $d$, because $d$, $\sigma_\vecx$, and $\sigma_\vecy$ 
are all nonnegative, and $(\frac{b_i}{\sigma_\vecy}-\frac{a_i}{\sigma_\vecx})^2$ and $(\frac{b_i}{\sigma_\vecy}+\frac{a_i}{\sigma_\vecx})^2$ are also nonnegative 
for all $i$.
This is very useful because, in equation (\ref{eqn:ms.lower}), 
the first term, $\msl(\vecx,\vecy,0)$, 
is already a lower bound of $\dist(\vecx,\vecy)^2$
by inequality (\ref{eqn:dist.lower}) , and
%equal to upper bound from the Cauchy-Schwarz inequality~\ref{eqn:cauchy}.
the lower bound can be refined further
nondecreasingly over the summation in the second term.
% computing $i$th summation.
If we stop the summation at $i=k$, for $k<d$, the intermediate result
is also a refined lower bounds of $\dist(\vecx,\vecy)^2$.
% because every remaining terms ($i>k$) are positive.
Similarly, in equation (\ref{eqn:ms.upper}), the first term, 
$\msu(\vecx,\vecy,0)$, 
is already an upper bound of $\dist(\vecx,\vecy)^2$ by inequality 
(\ref{eqn:dist.upper}) , and
%equal to upper bound from the Cauchy-Schwarz inequality~\ref{eqn:cauchy}.
the upper bound can be refined further
nonincreasingly over the summation in the second term.
% Note that MS-distance include bounds from inequalities (\ref{eqn:dist.lower}) 
% and (\ref{eqn:dist.upper}) in the first term.
This means we can stop the summation as soon as we find a bound good 
enough for the application under consideration. If we need the exact Euclidean distance, we can get it by continuing to the full summation. We summarize the above
properties in the following.
\begin{lem} (Monotone Convergence)
Let $\msl(\vecx,\vecy,k)$ and $\msu(\vecx,\vecy,k)$ be the lower and upper
bounds of MS-distance as defined above, respectively. Then the following
properties hold.
\begin{itemize}
\item $\msl(\vecx,\vecy,0)\leq\msl(\vecx,\vecy,1)\leq\cdots\leq\msl(\vecx,\vecy,d-1)\leq\msl(\vecx,\vecy,d)=\dist(\vecx,\vecy)^2$.
\item $\msu(\vecx,\vecy,0)\geq\msu(\vecx,\vecy,1)\geq\cdots\geq\msu(\vecx,\vecy,d-1)\geq\msu(\vecx,\vecy,d)=\dist(\vecx,\vecy)^2$.
\item $\msl(\vecx,\vecy,k)=\msl(\vecx,\vecy,k+1)$  if and only if $b_{k+1}/\sigma_\vecy=a_{k+1}/\sigma_\vecx$.
\item $\msu(\vecx,\vecy,k)=\msu(\vecx,\vecy,k+1)$  if and only if $b_{k+1}/\sigma_\vecy=-a_{k+1}/\sigma_\vecx$.
\end{itemize}
\end{lem}
\begin{lem}
For $0\leq k<d$, we can update $\msl(\vecx,\vecy,k)$ to $\msl(\vecx,\vecy,k+1)$,
and $\msu(\vecx,\vecy,k)$ to $\msu(\vecx,\vecy,k+1)$ in constant time.
\end{lem}
% \section{Geometric Property}
\paragraph{Fast Filtering.}
We must emphasize that $\msl(\vecx,\vecy,0)$ and $\msu(\vecx,\vecy,0)$ can
be used for fast filtering. Let $\phi$ denote a threshold for filtering
defined in some proximity search problem under consideration. If 
$\phi<\msl(\vecx,\vecy,0)$ in case of nearest search or 
$\phi>\msl(\vecx,\vecy,0)$ in case of farthest search, 
we do not need to consider this pair $(\vecx,\vecy)$ as a candidate, 
thus we can save time from computing their exact Euclidean 
distance.
% Before There is a property of MS-distance which can be used for fast 
% filtering.
% in previous sections are about to reducing computation time related to $d$.
% To check whether distance from a base vector(or query vector) to a target vector is over threshold or not,
% MS distance requires at least $O(1)$ computation.
% Therefore, if we have $n$ vectors to check, at least $O(n)$ computation 
% is required.
%
% In this section, we present a geometric property of MS-distance to avoid $O(n)$ bound computation.
% Note that dimension $d$ is positive and same for all $n$ vectors,
% Therefore, we can also use first bound of MS-distance divided by $d$ as following:

Precisely speaking, we map each $d$-dimensional vector 
$\vecx=[x_1,x_2,\ldots,x_d]$ 
into a pair of points, $\nvecx^+$ and $\nvecx^-$, in the 2-dimensional plane
such that $\nvecx^+=[\mu_\vecx,\sigma_\vecx]$ and 
$\nvecx^-=[\mu_\vecx,-\sigma_\vecx]$. Then

\begin{align}
\dist(\nvecx^+,\nvecy^+)^2&=\msl(\vecx,\vecy,0)/d\\
\dist(\nvecx^+,\nvecy^-)^2&=\msu(\vecx,\vecy,0)/d.
\end{align}

To see why it is useful in fast filtering, consider the case of finding
the nearest vector. For $d$-dimensional vectors in $V$ of size $n$,
we have $n$ pairs of points in the plane as in Figure~\ref{fig:geo}.
Since $\sigma_\vecx$ is nonnegative, exactly $n$ points lie on
or below $\mu$-axis. Let $\vecq$ be a query vector, and let $\hat{\vecq}^+$
denote the point mapped in the plane as defined above. Among these $n$ 
points lying on or below $\mu$-axis, 
let $\hat{\vecx}_i^-$ be the point that is nearest to
$\hat{\vecq}^+$. Note that the closest point from the query can be 
computed efficiently in $2$-dimensional space, for example, after 
constructing some space partitioning structures such as kd-trees
or R-trees, each query can be answered in poly-logarithmic search
time.

Then we can ignore all $d$-dimensional vectors $\vecx$
whose mapped point $\hat{\vecx}^+$ lies outside the circle centered 
at $\hat{\vecq}^+$ and of radius $\dist(\hat{\vecq}^+, \hat{\vecx}_i^-)$ 
in the plane, because they are strictly farther than $\vecx_i$ from $\vecq$.

\begin{figure}[ht]
\centering
\includegraphics[width=3in]{msdist.eps}
\caption{Fast filtering using $\msl(\vecx,\vecy,0)$ and $\msu(\vecx,\vecy,0)$. 
All $d$-dimensional 
vectors $\vecx$ whose mapped point $\hat{\vecx}^+$ lies outside the circle 
are strictly farther than $\vecx_i$ from $\vecq$.}
\label{fig:geo}
\end{figure}

\section{Estimating the Expected Difference Between Two Bounds}
We now turn to estimating the expected difference between 
$\msl(\vecx,\vecy,k)$ and $\msu(\vecx,\vecy,k)$.
Observe that $\msl(\vecx,\vecy,k)$ is almost the same as 
$\msl(\vecx,\vecy,k-1)$ if $b_{k}/\sigma_\vecy\approx a_{k}/\sigma_\vecx$.
Hence, in the worst case, $\msl(\vecx,\vecy,0)=\msl(\vecx,\vecy,d-1)<\msl(\vecx,\vecy,d)=\dist(\vecx,\vecy)^2$ when 
$b_{k}/\sigma_\vecy=a_{k}/\sigma_\vecx$ for
all $k=0,1,\ldots,d-1$, except $k=d$. Therefore, if we need a lower bound strictly 
better than $\msl(\vecx,\vecy,0)$, then we need to go through all $d$ 
refinement steps, which takes $O(d)$ time.
It is not difficult to see that this also applies to the case of 
$\msu(\vecx,\vecy,k)$.

However, this is unlikely to happen. Consider a random order for the last term 
in equation $\msl(\vecx,\vecy,k)$ and for the last term in 
equation $\msu(\vecx,\vecy,k)$. We show below that their expected
values increase and decrease linearly, respectively, as $k$
increases from $0$ to $d$.
Formally, let $(a_{\gamma(i)}, b_{\gamma(i)})$ denote the $i$th pair
in the random order. We measure the expected quality of the bounds
by the difference between the bounds, that is, 
$\msu(\vecx,\vecy,k)-\msl(\vecx,\vecy,k)$ as follows.

\begin{align}
\msu(\vecx,\vecy,k)-\msl(\vecx,\vecy,k)&=4d\sigma_\vecx\sigma_\vecy-2\sigma_\vecx\sigma_\vecy
\sum_{i=0}^k\left(\left(\frac{a_{\gamma(i)}}{\sigma_\vecx}\right)^2+\left(\frac{b_{\gamma(i)}}{\sigma_\vecy}\right)^2\right)\label{eqn:rgap.init}\\
%&=4d\sigma_\vecx\sigma_\vecy-2\sigma_\vecx\sigma_\vecy \frac{1}{d!} \sum_j^{d!}
%\sum_i^k((\frac{a_{\gamma_j(i)}}{\sigma_\vecx})^2+(\frac{b_{\gamma_j(i)}}{\sigma_\vecy})^2)\label{eqn:rgap.simp}\\
%&=4d\sigma_\vecx\sigma_\vecy-2\sigma_\vecx\sigma_\vecy \frac{1}{d!} \sum_j^{(d-1)!k}
%\sum_i^d((\frac{a_i}{\sigma_\vecx})^2+(\frac{b_i}{\sigma_\vecy})^2)\label{eqn:rgap.nrand}\\
&=4d\sigma_\vecx\sigma_\vecy-2\sigma_\vecx\sigma_\vecy \frac{k}{d} 
\sum_{i=0}^d\left(\left(\frac{a_i}{\sigma_\vecx}\right)^2+\left(\frac{b_i}{\sigma_\vecy}\right)^2\right)\label{eqn:rgap.res}\\
&=4d\sigma_\vecx\sigma_\vecy-4k\sigma_\vecx\sigma_\vecy\label{eqn:rgap.simplified}\\
&=4\sigma_\vecx\sigma_\vecy(d-k)\label{eqn:rgap.final}
\end{align}
Let us explain how we get Equation (\ref{eqn:rgap.res}) from 
(\ref{eqn:rgap.init}). Let $N$ denote the set of all pairs, and let $N^k$ 
denote the set of first $k$ pairs in the random order. 
Since each pair in $N$ is treated equally, $N^k$ is a random subset 
of $N$ of size $k$. Therefore, $\sum_{i=1}^k(a_{\gamma(i)}/\sigma_\vecx)^2$ 
is equivalent to take the total sum of $(a_i/\sigma_\vecx)^2$ with $i$ 
from $1$ to $d$ and divide it by $d/k$. We can also show this for 
$\sum_{i=1}^k(b_{\gamma(i)}/\sigma_\vecy)^2$ by a similar augment.

Equations (\ref{eqn:rgap.simplified}) and (\ref{eqn:rgap.final}) are because 
$\sum_{i=1}^d a_i^2=d\sigma_\vecx^2$ and $\sum_{i=1}^d b_i^2=d\sigma_\vecy^2$ 
by definitions of $a_i$ and $b_i$. By replacing each squared sum with $d$, 
that is , by applying $\sum_{i=1}^d (a_i/\sigma_\vecx)^2=\sum_{i=1}^d (b_i/\sigma_\vecy)^2=d$, we have Equation (\ref{eqn:rgap.final}). 
\begin{lem}
The expected value of $\msu(\vecx,\vecy,k)-\msl(\vecx,\vecy,k)$ is 
$4\sigma_\vecx\sigma_\vecy(d-k)$.\label{thm:gap}
\end{lem}
Because $\dist(\vecx,\vecy)^2$ always lies in between the two bounds,
the following also holds.
\begin{cor}
Both expected values of $\msu(\vecx,\vecy,k)-\dist(\vecx,\vecy)^2$ and
$\dist(\vecx,\vecy)^2-\msl(\vecx,\vecy,k)$ are 
at most $4\sigma_\vecx\sigma_\vecy(d-k)$.
\end{cor}
This shows a good theoretical expectation on the lower and upper bounds. 
This can justify that the MS-distance provides very tight bounds in a few 
refinement steps of a typical sequence.

\section{Applications : Proximity Searches}
The MS-distance can be used to application problems where the Euclidean 
distance is a measure for proximity or similarity of objects.
As a case study, we implemented the nearest neighbor search (NNS) and 
the farthest neighbor search (FNS) using the MS-distance. 
% in proximity searches: .experimental results

Given a set $X$ of $d$-dimensional vectors $\vecx_i$, for $i=1,\ldots,n$,
and a $d$-dimensional query vector $\vecq$, we use the following 
simple randomized algorithm for NNS. Initially, we set $\phi$ to the threshold 
given from the application under consideration or computed from 
the fast filtering in 2-dimension in Section~\ref{sec:ms-distance}.
\begin{enumerate}
% \item Construct a random sequence $\vecx_1,\vecx_2,\ldots,\vecx_n$ of the 
% vectors in $X$.
\item Consider the vectors in $X$ one at a time according to this sequence. 
  At the $i$th stage, we do the followings.\\
\begin{quote}
\begin{tabular}{ l }
% Construct a random sequence $\vecx_1,\vecx_2,\ldots,\vecx_d$ of the 
% vectors in $X$.\\
% \textbf{for} ~$i=1,2,...,n:$\\
%\hspace*{0.5cm}$\text{lowerbound}=\msl(\vecq,\vecx_i,0)$\\
\hspace*{0.5cm}\textbf{if} ~$\msl(\vecq,\vecx_i,0) < \phi:$\\
\hspace*{1.0cm}\textbf{for} ~$j=1,2,...,d:$\\
%\hspace*{1.5cm}$\text{lowerbound} = \msl(\vecq,\vecx_i,j)$\\
\hspace*{1.5cm}\textbf{if} ~$\msl(\vecq,\vecx_i,j) > \phi:$\\
\hspace*{2.0cm}\textbf{break;}\\
\hspace*{1.0cm}\textbf{if} ~$j=d$:\\
\hspace*{1.5cm}{$\phi=\msl(\vecq,\vecx_i,d)$;}\\
\hspace*{1.5cm}{NN = $i$;}
\end{tabular}
\end{quote}
\item return NN as the nearest neighbor of $\vecq$ with the squared Euclidean
distance $\phi$.
\end{enumerate}
Note that the first line of the pseudocodes filters out the vectors whose
distance to $\vecq$ is larger than $\phi$ as in the fast filtering 
in Section~\ref{sec:ms-distance}. In the \textbf{for} loop, we compute
$\msl(\vecq,\vecx_i,j)$ from $\msl(\vecq,\vecx_i,j-1)$ in
constant time.
From the last two lines of the pseudocodes, 
we update $\phi$ to the exact Euclidean 
distance between $\vecq$ and $\vecx_i$ and store the index as the \emph{current}
nearest neighbor (NN).
The algorithm for the farthest neighbor search is similar to this one, 
except that it uses $\msu(\vecx_i,\vecy,j)$ and maintains the maximum distance.

For empirical comparison, we implemented a linear search algorithm that simply 
computes distances from $\vecq$ to every $\vecx_i$ and chooses the one 
with the minimum distance. We also used the implementation of the cover tree~\cite{covertree}.
A cover tree is a data structure that 
supports fast nearest neighbor queries given a fixed 
intrinsic dimensionality~\cite{intrinsic}.

We tested these implementations on data sets from UCI machine learning 
archive~\cite{uci}.
We selected data sets $D$ from various dimensions (from $10$ to $100,000$),
and randomly selected $30$ queries points $Q\subset D$, and queried them 
on $D\setminus Q$.
We labelled the data set on $d$-dimension as ``D$d$''.
The data sets D500, D5000, D10000, D20000, D100000 were used in 
\emph{NIPS 2003 challenge on feature selection}~\cite{feature}.
% We include them to test performance of MS-distance based proximity 
% search algorithm on very high dimensions.
The test machine has one CPU, Intel Q6600 with 2.4GHz, 3GB memory,
and 32bit Ubuntu 10 operating system running on the machine.



\begin{figure}[ht]
\centering
\parbox[t]{2.6in}{
\includegraphics[width=2.70in]{filter.pdf}
\caption{Data filtered off in percentage.}
\label{fig:filter}
} \hspace{0.2in}
\parbox[t]{2.6in}{
\includegraphics[width=2.70in]{nnprep.pdf}
\caption{Preprocessing time for nearest neighbor search in log-scaled second.}
\label{fig:nnprep}
}
\end{figure}

Figure~\ref{fig:filter} shows the percentage of data filtered off.
For the data sets on relaxed dimensions,
the MS-distance filtered off over $95\%$ of data without lose of 
accuracy. For high dimensional data, MS-distance failed to filter off 
many data. Probably this is because the distances from queries to 
their nearest vectors tend to converge to the distances to their 
farthest vectors as described in~\cite{meaningful}.
This makes it hard to decrease (or increase in FNS) the 
threshold $\phi$ for the MS-distance enough to filter off many data.
However, on such high dimensions, both the linear search and the 
cover tree algorithm also show poor performance.

Figure~\ref{fig:nnprep} shows the preprocessing time of the 
MS-distance and the cover tree for NNS.
The time axis is log-scaled second.
This shows that the preprocessing time of the MS-distance is up to 
1000 times faster than the one in the cover tree.
This is because for the MS-distance it requires only $O(dn)$ time 
to compute the mean and the standard deviation values.

\begin{figure}[ht]
\centering
\parbox[t]{2.6in}{
\includegraphics[width=2.70in]{nnsearch.pdf}
\caption{Relative running time for the nearest neighbor search queries, normalized by linear search time.}
\label{fig:nnsearch}
} \hspace{0.2in}
\parbox[t]{2.6in}{
\includegraphics[width=2.70in]{fnsearch.pdf}
\caption{Relative running time for the farthest neighbor search queries, normalized by linear search time.}
\label{fig:fnsearch}
}
\end{figure}


Figure~\ref{fig:nnsearch} shows the time spent for NNS queries.
The graph shows the query time that is normalized by the 
linear search time. It is clear that the filtering algorithm based
on the MS-distance beats the linear search algorithm, even on high 
dimensional data in the results. The cover tree, which is 
designed exclusively for NNS, shows slightly better query performance 
than ours. However, the MS-distance is more general and flexible: it supports
addition of a new vector to the data set (our data structure) in $O(d)$ time
for computing the mean and the standard deviation values of the vector.
Deletion of a vector from the data set can be done in constant time.
Furthermore, the data structure for NNS can also be used for FNS.

Figure~\ref{fig:fnsearch} shows the time spent for FNS queries.
This is outstanding compared to the linear search algorithm. We hardly
know any other previous work achieving better performance than this.
% aspect of the MS-distance, because most tree based 
% algorithms hardly achieve such a result.
% On FNS, we have an example that our algorithm is slightly slower than 
% linear search.
% However, we already know that MS-distance becomes exact Euclidean distance after $d$ refinement steps.
% Therefore, we can say that even in worst case scenario that every filtering failed,
% the theoritical running time is still bounded as same as linear distance computation.

\begin{table}[ht]
\caption{Data sets}
\label{dtable}
\begin{center}
\begin{tabular}{l l r | l l r }
%\multicolumn{1}{c}{\bf }  &\multicolumn{1}{c}{\bf DESCRIPTION}
\\ \hline \\
Data Label & Name & $\#$ of vectors & Data Label & Name & $\#$ of vectors \\
\\ \hline \\
%D1000 & 1000 & 10000 & D1000 & 1000 & 10000 \\
D10 & Page Blocks & 5473 & D64 & Optical Recognition & 5620 \\
D11 & Wine Quality & 6497  & D86 & Insurance Company & 5822 \\
D16 & Letter Recognition & 20000 & D90 & YearPredictionMSD & 515345 \\
D19 & Image Segmentation & 2310 & D167 & Musk2  & 6597 \\
D22 & Parkinsons Tel & 5875 & D255 & Semeion & 1593 \\
D27 & Steel Plates Faults & 1941 & D500 & Madelon & 4400 \\
D37 & Statlog Satellite& 6435 & D617 & ISOLET & 7795 \\
D50 & MiniBooNE & 130064 & D5000 & Gisette & 13500 \\
D55 & Covertype & 581012 & D10000 & Arcene & 900 \\
D57 & Spambase & 4601 & D20000 & Dexter & 2600 \\
D61 & IPUMS Census & 233584 & D100000 & Dorothea & 1950 \\
\\ \hline \\
\end{tabular}
\end{center}
\end{table}

\section{Conclusion}
We introduce a fast distance bounding technique, called the MS-distance, 
by using the mean and the standard deviation values. 
The MS-distance between two vectors 
provides upper and lower bounds of Euclidean distance between them
in constant time, and these bounds converge monotonically to 
the exact Euclidean distance over iteration.
The MS-distance can be used to application problems where the Euclidean 
distance is a measure for proximity or similarity of objects. The 
experimental results show that our method is efficient enough even to 
replace the best known algorithms for proximity searches.

\subsection*{Acknowledgments}
This work was supported by the National Research Foundation of Korea
Grant funded by the Korean Government (MEST) (NRF-2010-0009857).

\small{
\bibliographystyle{unsrt}
\bibliography{ref}
}
\end{document}
