\documentstyle[11pt]{article}

%\documentstyle[times,art10,twocolumn,latex8]{article}

\pagestyle{empty}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{corollary}[theorem]{Corollary}
\def\Section{\section}
\def\SubSection{\subsection}
\def\var{\hbox{Var}}
\def\stdnor{e^{-|x|^2}}
\def\gennor{e^{-xAx^T}}
\def\unicond{{{l^2 \sigma}\over{2\sqrt{n\pi}}}}
\def\stdcond{{1\over{20 e \sqrt{\pi} n^{3/2}}}}
\def\gencond{{{\lambda_n}\over{\lambda_1 32\sqrt{2}(1+\sqrt{3}) e
n\ln n}}}
\def\orthant{\{(x_1,x_2,...,x_n)\in R^n:\, x_i\geq 0\}}
\def\vol{{\rm vol}}
\def\dist{{\rm dist}}
\def\mindim{10}
\def\mixone{(20 e \sqrt{\pi} n^{3/2})^2}
\def\mixtwo{({{\lambda_1}\over{\bar{\rho_n}\lambda_n}}32\sqrt{2}(1+\sqrt{3}) e
n\ln n )^2}
\def\pr{\noindent{\bf Proof.\hbox{ }}}
\def\prob{\hbox{Prob}}

\begin{document}  

\title{Sampling according to the multivariate normal density}

\author{Ravi Kannan\\School of Computer Science\\Carnegie-Mellon
University\\Pittsburgh, PA 15213\thanks{Supported in 
part by NSF grant CCR 9528973.}\\{\bf ravi.kannan@cs.cmu.edu}
\and 
Guangxing Li\\Department of Mathematics\\Carnegie-Mellon
University\\Pittsburgh, PA 15213 }

\maketitle
\thispagestyle{empty}

\begin{abstract}
This paper \footnote{To appear in the Proceedings of the IEEE Symposium
on the Foundations of Computer Science, 1996}
deals with the normal density of $n$ dependent random
variables. This is a function of the form : $$c e^{-x^TAx}$$
where $A$ is an $n\times n$ positive definite matrix, $x$ is the $n-$vector of the random
variables and $c$ is a suitable constant. 

The first problem we consider is
the (approximate) evaluation of the integral of this function over the
positive orthant 
$$\int_{x_1=0}^\infty \int_{x_2=0}^\infty \ldots \int_{x_n=0}^\infty  
                c e^{-x^TAx}.$$
This problem has a long history and a substantial literature 
. Related to it is the problem of drawing a sample from
the positive orthant with probability density (approximately) equal to
$c e^{-x^TAx}$. We solve both these problems here in polynomial time
using rapidly mixing Markov Chains. For proving rapid convergence of the
chains to their stationary distribution, we use a geometric property 
called the Isoperimetric Inequality. Such an inequality has been 
the subject of recent papers for general log-concave functions. We use these techniques,
but the main thrust of the paper is to exploit the special property 
of the normal density to prove a stronger inequality than for general
log-concave functions. 

We actually consider first the problem of drawing a sample according to the
normal density with $A$ equal to the identity matrix from a convex set
$K$ in ${\bf R}^n$ which contains the unit ball. This problem is motivated
by the problem of computing the volume of a convex set in a way we explain
later. Also, the methods used in the solution of this and the orthant problem 
are similar. 

\end{abstract}

\Section{Introduction}\label{sec:intro}
The problem of integrating the multivariate 
normal density over certain convex sets, especially orthants is a
classical
problem. It is related to other problems in probability and statistics. 
The approximation of orthant probabilities of the normal distribution has
an extensive
literature, see for example the results and bibliographies in 
\cite{mar}, \cite{satish}, \cite{moran}
and \cite{evans}. The problem of integration is well-known to be related to the
problem of sampling and many methods are known for estimating the integral based
on a number of samples. (For example, see \cite{mar} , \cite{evans}). However,
the difficulty with the known methods is that one needs a number of samples
that grows exponentially with the number of dimensions $n$. For this reason,
it has been noted (\cite{mar},\cite{evans}) that as a rule, the estimation
of the integral is only feasible for very modest $n$.

There is an even more fundamental problem
: that of sampling according to 
$e^{-x^TAx}$ in, say, the positive orthant quickly. (I.e., in time bounded by a polynomial
in $n$.)

[We use the phrase ``sampling according to a (nonnegative real valued function)
$F(x)$ in a region $R$'' to mean : pick $x\in R$ at random from the probability
distribution over $R$ with density $F(x)/\int_R F$.]

It is easy to
sample in the whole space according to $e^{-x^TAx}$ : one simply draws a sample
according to the standard normal density $ e^{-|x|^2}$ by any one of a number
of schemes and applies a (suitable) linear transformation to it. (We omit the
details.) But, then if one wishes to sample in the positive orthant according to this
function, by say Rejection Sampling, we expect to draw 
$${\int_{ {\bf R}^n} e^{-x^TAx}\over  \int_ { {\bf R}^n_+ } e^{-x^TAx} }$$
samples before one lands in the positive orthant ${\bf R}^n_+$. This ratio in general
will grow exponentially with $n$. 

We use instead a Markov Chain with the desired steady state density - 
$e^{-x^TAx} /\int_{ {\bf R}^n_+} e^{-x^TAx}$ in ${\bf R}_+^n$. 
The chain itself is a natural
``Metropolis'' algorithm. The Metropolis algorithm can be used in a wide range
of contexts - see \cite{ds} for a recent survey. In our cases, the chain is easily
described : Suppose we wish to sample from a set $R$ according to function
$F$. We first choose a ``step size'' $\delta  $ suitably. (To be described later.)
Suppose the Chain is at a point $x\in R$. Then we choose a random point $y$ from 
a ball of radius $\delta $ around $x$ uniformly. If $y$ is not in $R$, we stay at $x$.
Otherwise, if $F(y)\geq F(x)$, we move to $y$, else, we move to $y$ with probability
$F(y)/F(x)$. Formally, the probability density $P_{xy}$ of moving from
$x$ to $y$ is given by :
$$P_{xy}= \left\{
\begin{array}{ll}
{1\over \vol {\delta  B} }\min \{ 1, {F(y)\over F(x)} \} & \hbox{   if } x\not= y\in R;\; 
 |x-y|\leq \delta  ;\\ 
0 & \hbox{\rm otherwise }.\\
\end{array}
\right. %}
$$
$$P_{xx}=1-\int_{y \not= x} P_{xy} \ dy.$$

We call this the {\bf Metropolis Chain according to $F$}. It is easily see that
its steady state distribution is $F(x)/\int_R F$. 

Our proof that this chain
``mixes rapidly'' (i.e., converges to the steady state fast) uses techniques
developed recently. These techniques were developed to solve the problem of
estimating the volume of a convex set (this problem itself motivates another 
sampling question in this paper) in \cite{dfk},\cite{ls},\cite{ak},\cite{df}
(see also \cite{k} for a survey.) Using these techniques, one relates the convergence
to a purely geometric property called the Isoperimetric Inequality. In this paper,
we exploit the special nature of the normal density (specifically the fact that its
logarithm declines quadratically in every direction) to prove a sharper Isoperimetric
Inequality than the one for general log-concave functions proved in \cite{ak},
\cite{polson}, \cite{df} and other papers. 
We remind the reader that a function is said to be log-concave if its logarithm is
concave. 

The Markov Chain yields samples (nearly) according to $F(x)$. We show how to use a polynomial
(in $n$ and an error parameter) number of samples to estimate the integral of $F$
over the positive orthant.
To avoid the problem of high variance, we use a trick often used in volume
estimation algorithms -  the integral is estimated as a product of polynomially
many quantities, each of low variance. 

Our sampling algorithm uses $O^*(\Lambda^2 n^3)$ steps (of the random walk)
for the first sample and
$O^*(\Lambda^2 n^2)$
steps for each subsequent nearly independent sample, where $\Lambda$ is the
condition number
of $A$, namely the ratio of the largest eigenvalue to the smallest. 
We also show that we can estimate the orthant integral with relative error $\varepsilon$
with $O^*(\Lambda^2 n^2/\varepsilon ^2 )$  samples. 

Before tackling the orthant problem, 
we first consider the problem of drawing a sample according to
the standard normal (we use this phrase here to mean the case of $A=I$)
from a convex set $K$ containing the unit ball. For this problem, our
algorithm finds
the first sample in $O^*(n^4)$ steps and each subsequent (nearly independent)
sample in 
$O^*(n^3)$ steps. [The $O^*$ notation hides constant factors, log factors and factors
depending 
upon the accuracy needed.]

A crucial parameter in the algorithm is the ``step size'' $\delta $. In the present 
paper, we use $\delta  = \Theta ({1\over n})$ for the latter problem
(of sampling from a convex set containing the unit ball according to the
standard normal density) to achieve the bounds on the number
of steps claimed. By using the techniques of \cite{kls}, we hope to show that
a step size 
$\delta  = \Theta ({1\over \sqrt n})$ may be used for the subsequent samples, 
which would  then improve the number of steps needed for each subsequent sample
in this problem 
to $O^*(n^2)$ steps. Now by direct integration, it is not difficult to show that 
we may discard the part of $K$ beyond a ball of radius $O(\sqrt n)$ from the origin
as this has negligible measure. Thus in this case we have that $B\subseteq K\subseteq
c\sqrt n B$. The number of steps needed in such a situation for the algorithm in
$\cite{kls}$ to draw each subsequent sample is $O^*(n^3)$ (see Theorem 2.2 of \cite{kls}).
So, if the improvement holds, then at least under the normal density, samples come
faster than the best known for the uniform density. This may ultimately 
help improve the complexity of the difficult problem of estimating the volume of
a convex set. 

\Section{Isoperimetric Inequalities}
\label{sec:iso}
Isoperimetric inequalities play an important role in analysis of the
mixing time
of Markov
chains. This is because they can be related to  the conductance of the chain, a
notion 
introduced by Jerrum and Sinclair \cite{js} which is very useful in the analysis of convergence. 

Roughly speaking an isoperimetric inequality in this context says that for any
partition of
a convex set $K$ into two parts, the measure (according to the desired steady
state function
$F$) of the dividing surface is not too small compared with the minimum of the
measures of the
two parts. The version we use here is slightly different - here, we consider a
division into three
parts $K_1,K_2,K_3$ where $K_1,K_2$ correspond to the two main parts described
above and 
$K_3$ may be thought of as the set of points within a certain distance $\sigma
/2$ of 
the dividing surface above.
We now state precisely our isoperimetric inequality;
this is the central new theorem we prove and use in the paper. One important
feature
of this is that there is no dependence on the diameter of the convex set
$K$. 

\begin{theorem}
\label{theorem:iso1}
Let $F(x)= \stdnor g(x)$, where $g(x)$ is an arbitrary log-concave function.
Suppose $K$ is
a convex set in ${\bf R}^n$. $K_1$ and $K_2$ are two measurable subsets of
$K$, and the
distance between $K_1$ and $K_2$ (i.e., $\inf_{x\in K_1,y\in K_2} |x-y|$)
is $\sigma$. $K_3 = K\backslash (K_1\cup K_2)$. 
The following inequality holds.

\begin{equation} 
\label{iso1}
\int_{K_3} F(x)\, dx \geq 
{{2\sigma e^{-\sigma^2}}\over{\sqrt{\pi}}} \min \{
\int_{K_1} F(x)\, dx,\, \int_{K_2} F(x)\, dx\}
\end{equation}

\end{theorem}

{\bf Idea of Proof} : 
We use a powerful technique due to Lov\'asz and Simonovits \cite{ls} to prove this
result. This technique uses their Localization Lemma to reduce the $n-$ dimensional
case to essentially a 1 dimensional problem. In this situation, this lemma implies that
the assertion in the Theorem holds provided for each line segment $[a\ b]$ and each
linear function $\ell (\cdot )$, we have the following :
$$
\begin{array}{l}
\int_{K_3\cap [a\ b]} F\ell^{n-1} \geq {{2\sigma e^{-\sigma^2}}\over{\sqrt{\pi}}} \\ 
\min \{
\int_{K_1\cap [a\ b]} F\ell^{n-1} dx, \int_{K_2\cap [a\ b]} F\ell^{n-1} \} \ \ 
(\ref{iso1}')
\end{array}
$$
where now all the integrals involved are just one dimensional integrals. (Intuitively, it says
that the inequality asserted in the Theorem holds provided it holds for every 
pointed cone with an infinitismal cross section.) As we move along the line $[a\ b]$,
$x$ is a function of one parameter $t$. We abuse notation and let all the functions
involved be functions of $t$. Also, we let $G(t) = g(t) \ell (t) ^{n-1}$ and note that
$G$ is log-concave. Then, we have that $F(t) = e^{-t^2 + \delta_1 t + \delta_2} G(t)$ for 
some constants $\delta_1 ,\delta_2$. 
In this sketch of the proof, we also deal only with the case
when $K_3\cap [a\ b]$ is an interval; the more general case can be dealt with the same
ideas and the full proof is deferred to the final paper. By the hypothesis of the
Theorem, we have in this case that $K_3\cap [a\ b]$ is an interval of length at least
$\sigma $. It will suffice to prove that for every real number
$\nu$ :
\begin{eqnarray}
\int_\nu^{\nu + \sigma} e^{-t^2 + \delta_1 t + \delta_2} G(t) dt
\geq {{2\sigma e^{-\sigma^2}}\over{\sqrt{\pi}}} \nonumber\\ 
\min \{
\int_{-\infty}^{\nu} e^{-t^2 + \delta_1 t + \delta_2} G(t) dt, 
\int_{\nu + \sigma}^{\infty} e^{-t^2 + \delta_1 t + \delta_2} G(t) dt\}.\label{iso3}
\end{eqnarray}


There are positive reals $c_1$ and $c_2$ such that $G(\nu) = c_1 e^{c_2 \nu}$ and 
$G(\nu +\sigma) = c_1 e^{c_2(\nu + \sigma)}$. By using the log-concavity of
$G(t)$, we can see that the following is true:

$$
G(t) 
\left\{
\begin{array}{ll}
\leq c_1 e^{c_2 t} & t\in (-\infty, \nu)\cup (\nu +\sigma, \infty)\\
\geq c_1 e^{c_2 t} & t\in (\nu, \nu+\sigma)
\end{array}
\right.
$$

So we can replace $G(t)$ in (\ref{iso3}) by $e^{c_2 t}$ and prove the resulting 
inequality :
$$
\begin{array}{l}
\int_\nu^{\nu + \sigma} e^{-t^2 + \delta_1' t + \delta_2'}dt
\geq {{2\sigma e^{-\sigma^2}}\over{\sqrt{\pi}}} \\ 
\min \{
\int_{-\infty}^{\nu} e^{-t^2 + \delta_1' t + \delta_2'} dt,\, 
\int_{\nu + \sigma} e^{-t^2 + \delta_1' t + \delta_2'} dt\}\ \  (\ref{iso3}')
\end{array}
$$
After cancelling common factors, it suffices to prove the following. 

\begin{lemma}
The following inequality holds:
\begin{equation}
\label{iso4}
\int_\nu^{\nu + \sigma} e^{-t^2}dt
\geq {{2\sigma e^{-\sigma^2}}\over{\sqrt{\pi}}}\min \{
\int_{-\infty}^{\nu} e^{-t^2} dt,\, 
\int_{\nu + \sigma}^\infty e^{-t^2} dt\}
\end{equation}
\end{lemma}

\pr Without loss of generality, we assume that $\nu > 0$. The other case can be
dealt with exactly symmetrically. Consider the following function of $x$: 
$$\phi (x) ={{\int_x^{x+\sigma} e^{-t^2} dt}\over{\int_{x+\sigma}^{\infty}
 e^{-t^2} dt}}$$
This function is increasing for $x>0$, as seen by differentiating it. 
This means that it suffices to consider the case when $\nu = 0$. In this case, we have
$$\int_0^{\sigma} e^{-t^2} dt \geq \sigma e^{-\sigma^2}.$$
Also, we have
$$\int_{\sigma}^{\infty} e^{-t^2} dt \leq {{\sqrt{\pi}}\over 2}$$
Therefore the lemma follows and the theorem is proved.

For the general normal density function, we have similar result and the proof of it 
is similar to that of the standard case.
\begin{theorem}
\label{theorem:iso2}
Let $F(x)= e^{-xAx^T} g(x)$, where $A$ is a positive definite matrix and 
$g(x)$ is an arbitrary log-concave function. $K$ is
a convex set in $R^n$. $K_1$ and $K_2$ are two subsets of $K$. The distance between
$K_1$ and $K_2$ is $\sigma$. $K_3 = K\backslash (K_1\cup K_2)$. 
The following inequality holds, where $\lambda_1$ is the smallest
eigenvalue of $A$ :  
\begin{equation}
\label{iso5}
\int_{K_3} F(x) dx \geq {{2\sqrt{\lambda_1}\sigma e^{-\sigma^2}}\over{\sqrt{\pi}}} \min \{
\int_{K_1} F(x) dx,\, \int_{K_2} F(x) dx\}
\end{equation}
\end{theorem}
\noindent{\bf Proof.} Again, with the Localization Lemma, we can reduce it to
the one dimensional case. Let $v$ be the direction of this $[a\ b]$, then
\begin{eqnarray*}
F(t) &=& e^{-(a+tv)A(a+tv)^T} g(a+tv)\\
&=& e^{-(t^2\sigma_1 + t\sigma_2 + \sigma_3)} g(a+tv)
\end{eqnarray*}
where $\sigma_1 = vAv^T\geq \lambda_1$, $\sigma_2 = aAv^T + vAa^T$ and $\sigma_3 = aAa^T$.
Note that $g$ is again log-concave. By following the same procedure
as in the proof of the last Theorem, we can see that it suffices to prove the following
$$
\int_\nu^{\nu + \sigma} e^{-\sigma_1 t^2}dt
\geq {{2\sqrt{\lambda_1}\sigma e^{-\sigma^2}}\over{\sqrt{\pi}}}\min \{
\int_{-\infty}^{\nu} e^{-\sigma_1 t^2} dt,\, 
\int_{\nu + \sigma} e^{-\sigma_1 t^2} dt\}
$$
With a simple variable substitution $u=\sqrt{\sigma_1}t$  and (\ref{iso4}), as well as the fact
that $\sigma_1 \geq \lambda_1$ we can derive (\ref{iso5}). This completes the proof.

\Section{Rapid Mixing}

\SubSection{Local conductance, step size and conductance}
\label{sec:lsc}

A central quantity of interest in rapid mixing proofs is the notion of
conductance
introduced by Jerrum and Sinclair. 
Suppose we have a Markov Chain in $K$ with steady state probability density
$\pi (\cdot )$. The conductance of the chain is defined 
as the infimum over all (measurable) partitions of $K$ into
two sets
$S_1,S_2$ of the quantity below : 
\begin{equation} \label{def:cond}
{\int_{x\in S_1}\int_{y\in S_2} \pi (x) P_{xy} \ dy\ dx
       \over
\min\{\int_{S_1} \pi(x)\, dx,\ \int_{S_2}\pi (y)\, dy\} }.
\end{equation}

If $\pi (S_1)\leq 1/2$, then quantity is the conditional 
probability of
going from $S_1$ to $S_2$ in one step conditioned on starting in the
steady state
in $S_1$. The central result of Jerrum and Sinclair is that under mild
conditions on the chain (which hold for us), high conductance implies rapid convergence. 
Below we will use a more directly relevant version of this from \cite{ls}. 

In the analysis of convergence rates of geometric random walks, a lower bound on
conductance has often been proved using an Isoperimetric Inequality. (See \cite{dfk},
\cite{ls},\cite{ak}.) The connection between conductance and Isoperimetry  
is intuitively not difficult to see. But to carry out the argument, we need another
quantity, namely the local conductance.

The local conductance at a point $x$ in space is defined as the
probability
of making a proper move from $x$ in one step. For example, if only an
exponentially small
volume of the $\delta $ ball around $x$ is contained in the convex set, then the
local
conductance of the Metropolis random walk is exponentially small. In this case,
we have
to make exponentially many steps before moving away from $x$, so we cannot
guarantee
rapid convergence. 
[In a sense, one may view local
conductance as
conductance specialized to single element sets.] 

We establish the following result about the conductance in terms of local
conductance and 
step size. 
In a manner similar to the proof of Theorem 3.2 in \cite{ls}, we 
relate conductance to Isoperimetric Inequality, but in 
order to achieve better coefficient,
we will carry out the estimation more carefully. 
Then we use our Isoperimetric Inequality (Theorem \ref{theorem:iso2} ). 
The full proof is deferred to the
final paper. 
\begin{theorem}
\label{theorem:lsc}
Suppose 
$$F(x)=e^{-x^TAx}g(x),$$
where $g$ is a log-concave function and $A$ is a positive definite
matrix with least eigenvalue $\lambda_1$. Suppose 
we are doing Metropolis
random walk
with respect
to $F(x)$. Also, suppose the step size is $\delta $ and  the local conductance
everywhere is at least a certain $l>0$. Then the conductance is at least
$$\Phi \geq {l^2\delta \sqrt\lambda_1\over 2\sqrt {n\pi }}.$$

\end{theorem}

\SubSection{M-distance}
Now we relate conductance to convergence. To do so, we use a quantity 
called the ``M-distance'' which measures how far apart two distributions are. 
The M-distance between two distributions $P$ and $Q$ is defined as (it is really
not a distance, but is however a convenient measure of the difference) :
$$M(P,Q)=\sup_{S: 0<Q(S)\leq 1/2}
{{|P(S)-Q(S)|}\over {\sqrt{Q(S)}}}.
$$
See \cite{kls} for a note about why the $M$ distance is convenient to use. 
The more common measure used is the total variation distance between $P$ and $Q$
which is denoted $|P-Q|_{\hbox{tv}}$ and defined as
$$\sup_S |P(S)-Q(S)|.$$
Clearly the M-distance is an upper bound on the total variation distance. 

In the context of this paper, we run a Markov chain starting from distribution
$Q_0$ and achieve distribution $Q_k$ at step $k$. Our steady state distribution 
is denoted $Q$. Therefore we want to know the bound for $M(Q_k,Q)$. 
Our main tool is the following theorem from \cite{ls}.
\begin{theorem}
\label{theorem:mixing}
Suppose $\Phi$ is the conductance of the Markov chain, we have
\begin{equation}
M(Q_k, Q) \leq M(Q_0, Q) \left(1-{{\Phi^2}\over 2}\right)^k
\label{equation:mixing}
\end{equation}
\end{theorem} 

In many contexts, it is not sufficient to get one sample. We like many
independent samples all from (nearly) the same distribution $Q$. To get
this, one may of course restart the random walk from a different (independently
chosen) point. We will see that this is not necessary. In fact, in general, for
geometric chains as in this paper, it turns out that it is easy to start from
distributions $Q_0$ which have $M(Q_0,Q)\leq e^{O^*(n)}$ (as we show later)
and so applying the above
Theorem, one sees that in $O^*(n \log (1/\varepsilon )/\Phi^2)$ steps, we can get a 
sample from a distribution with M-distance at most $\varepsilon$ from the desired
$Q$. Intuitively, now for the second sample, it is better to start from the sample
we already have, but this introduces dependence. However, we will argue below that
after $O(1/\Phi^2)$ steps, we will have a second sample which is ``nearly'' independent
of the first and is from a distribution with M-distance at most $\varepsilon$ again to
$Q$. This improves steps needed for subsequent samples by a factor of $O^*(n)$. 
Something akin to this happens for general Markov Chains, as first pointed out in
\cite{aldous} for discrete chains. Here, we will follow the simple proof in
\cite{kls}. First, we need a definition of a particular weak form of independence
that turns out to be sufficient for us. 

{\bf Definition } : Suppose $X,Y$ are two random variables (with values in two possibly
different $\sigma-$ algebras) and $\varepsilon $ any positive real. We say that 
$X,Y$ are $\varepsilon-$independent if
$$|\prob (X\in A\ , \ Y\in B) - \prob (X\in A)\prob (Y\in B) |\leq\varepsilon$$
for all measurable sets $A,B$. 

\begin{corollary}\label{subsequent}
Suppose we start from a point $w_0$ drawn according to distribution
$Q_0$ and run the Metropolis random walk with conductance $\Phi$ for
$t$ steps to get a point $w_t$. Then $w_o$ and $w_t$ are $\tau-$independent
where $\tau = [2M(Q_0,Q)+1] \left( 1-{1\over\Phi^2}\right)^t$.
\end{corollary}

\pr To see the $\tau$-independence, we argue as follows. Let $A,B$ be any
two measurable sets. Let $f(A,B)$ be defined as 
\begin{eqnarray}{l}\label{INDEP1}
\big|\prob (w_0\in A,w_t\in B)-\prob (w_0\in A)\prob (w_t\in
B)\big|\nonumber\\
= Q_0(A)\big|\prob (w_t\in B\ | \ w_0\in A) - Q_t (B)\big|
\end{eqnarray}
Let $Q_0'$ be the distribution of $w_0$ conditioned on it being in $A$,
i.e., $Q_0'(S) = Q_0(S\cap A)/Q_0(A)$ for any measurable $S$.
Then, $\prob (w_t\in B | w_0\in A)$ is the distribution of $w_t$ if we start
with $w_0$ drawn from $Q_0'$. So, applying the
theorem above, we get
\begin{eqnarray}\label{INDEP2}
|\prob (w_t\in B\ | \ w_0\in A) - Q (B)|\nonumber \\ 
\leq \sqrt { Q(B) } M(Q_0', Q)
\left( 1-{1\over\Phi^2}\right)^t.
\end{eqnarray}
For any measurable $S$,
$$
Q(S)=  {1\over Q_0(A)} Q(S\cap A) - \left({1\over Q_0(A)}-1\right)
Q(S\cap A) + Q(S\setminus A),
$$
and hence
\begin{eqnarray*}
|Q_0'(S)-Q(S)| &\le & {1\over Q_0(A)} | Q_0(S\cap A) -Q(S\cap A)|\\
&& + \left({1\over Q_0(A)}-1\right) Q(S\cap A) + Q(S\setminus A).
\end{eqnarray*}
So, we have $M(Q_0',Q)\leq {1\over Q_0(A)}[M(Q_0,Q) +1]$
and now using (\ref{INDEP1}) and (\ref{INDEP2}), and the 
theorem, we get the claimed $\tau$-independence.

In the rest of this write-up, we will make assertions of the form ``each subsequent
nearly independent sample takes ..... number of steps''. The specifications of how nearly
independent are postponed to the final paper. But such assertions will all be based
on the Corollary above and independence will mean the $\varepsilon$ independence
we have above. 

\Section{The standard normal density case}
\label{sec:std}
Now we consider the case: sampling according to the standard normal density
$\stdnor$ in a convex set $K$ with $B(0,1)\subseteq K$. This is split to
several
steps:
penalty function, step size and local conductance, conductance and sampling.

\SubSection{Penalty function}
\label{sec:std:penalty}
If we run the Metropolis random walk in $K$, according to just $F$, we run into the
following problem : $K$ may have very sharp corners and the local conductance may be
exponentially small in such a corner. 

One way to handle sharp boundary is introducing penalty function of the convex
set as used in \cite{ak},\cite{ls} etc. The idea
is: instead of restricting the random walk to $K$, we allow it to 
go outside of $K$, but ``dampen'' the probability outside more as we go away from $K$
. We dampen by using some appropriate penalty function to
penalize
being outside.
\begin{definition}
Given a convex set $K$ with the origin in its interior
, the gauge function of it is defined as follows:

$$\varphi_K(x) = \left\{
\begin{array}{ll}
1 & \hbox{\rm if } x\in K;\\
\min\{t: x\in tK\} & \hbox{\rm if } x\not\in K.
\end{array}
\right.
$$
\end{definition}

It is well-known that 
$g(x) = e^{-c(\varphi_K(x) - 1)} $ is a log-concave function for any positive
constant $c$. We use this as the penalty function. 

We will do a Metropolis random walk in all of ${\bf R}^n$ according to the function
$$F(x) = e^{-|x|^2 - c(\varphi_K(x) - 1)}$$
for a certain number of steps to be specified and reject if the resulting sample
is not in $K$. 
We will choose $c$ big enough to guarantee that the probability of being
outside of $K$ is small, so that we do not reject too many times before we accept a
sample. 

It can be shown that with $c={3\over 2} n$, the probability outside
is at most $3$ times of that of inside. [This is shown by showing that
it holds in every infinitismal pointed cone emanating from the origin.]
The proof is deferred to the final paper.

From now on, we sample in the whole of $R^n$ according to the following
function:
\begin{equation}
\label{equation:std}
F(x) = e^{-|x|^2 - {3\over 2}n(\varphi_K(x) - 1)}
\end{equation}
\begin{proposition}
\label{prop:std:inout}
The probability inside $K$ of distribution defined by (\ref{equation:std}) 
is at least $0.25$.
\end{proposition}

\SubSection{Step size and local conductance}
\label{sec:std:local}
Choosing appropriate step size is tricky. Usually, small step size
gives large
local conductance. But if it is too small, it lowers the lower bound on conductance 
that we use (namely from Theorem \ref{theorem:lsc}). So, there is 
a balance between small and large step size. Here it turns out
that we want
to choose step size with
order $1\over n$. In order to achieve the best bound for conductance, we use a 
parameter $1\leq c\leq 3$, and choose $1\over {cn}$ as the step size. $c$ will
be identified when
we discuss conductance. We prove in the final paper that : 
\begin{lemma}
\label{lemma:local}
With step size $\delta ={1\over {cn}}$, the local conductance at every point is at least
$${1\over 2}e^{-{1\over{c^2n^2}}-{3\over {2c}}}.$$
\end{lemma}

\SubSection{Conductance and convergence rate}
\label{sec:std:con}
We can plug in the above value for $l$ (and $\delta = 1/(cn)$) into Theorem \ref{theorem:lsc}
and choose the best value for $c$. It is convenient to ignore the $e^{-{1\over{c^2n^2}} }$ term
and then to differentiate with respect to a new variable $\lambda = 1/c$. This tells us to
choose $c=3$.

Now to describe the random walk, one only needs to specify the starting distribution, for which
we just use the uniform distribution over the unit ball. Then we get the following :

{\bf Algorithm } Starting from a point $w_0$ picked according to the uniform density from
the unit ball, do the Metropolis random walk in all of ${\bf R}^n$ according to the function
$F(x)$ defined in (\ref{equation:std}) using step size $\delta = 1/(3n)$. [We do the walk
for a number of steps chosen to ensure that the $M$ distance to the steady state is within
a desired bound (see below); if the resulting point is not in $K$, we reject and repeat.]

In the final paper, we show the following facts about this walk. 

\begin{theorem}
The conductance of the above Markov chain is at least 
\begin{equation}
\label{equation:std:cond}
\stdcond 
\end{equation}
\end{theorem}

\begin{proposition}
$$M(Q_0,Q)\leq 1 + e \pi^{n/2} n!.$$
\end{proposition}

From (\ref{equation:std:cond}) and (\ref{equation:mixing}), we get
that we can get
the first
sample in number of steps 
$$O^*(n^4).$$
Also, from Corollary \ref{subsequent}
the number of steps for each subsequent (nearly independent) sample is 
$$O^*(n^3).$$

\Section{Orthant case with general normal distribution}
Now we consider the general normal distribution $f(x) = e^{-x^TAx}$ 
in the orthant
$\Omega = \{ x: x_i\geq 0, i = 1,2,...,n\}$. $A$ is a positive definite 
matrix with eigenvalues 
$0< \lambda_1\leq\lambda_2\leq\cdots\leq\lambda_n$. Let $\Lambda = \lambda_n/\lambda_1$. 

\SubSection{Penalty function}
\label{sec:gen:penalty}
Again we introduce a penalty function to overcome the sharp boundary problem.
It is :
\begin{definition}
Function $\varphi(x)$ is defined as follows:
$$
\varphi(x) =
\left\{ 
\begin{array}{ll}
0 & \hbox{if $x_i\geq 0$ for each $i$};\\
\max \{ -x_i\sqrt{n}:\, x_i<0\} & \hbox{otherwise.}
\end{array}
\right.
$$
\end{definition}
In other words, for $x$ not in the orthant, if we draw a line starting from it 
and along direction $\tau=({1\over{\sqrt{n}}}, {1\over{\sqrt{n}}}, ...,
{1\over{\sqrt{n}}})$, then $\varphi(x)$ is the distance between $x$ and the 
intersecting point of this line with the boundary of the orthant. It is easy to
show that this function is convex.

The proof of the following proposition is technically complicated and long
and is deferred to the final paper. 

\begin{proposition}
Suppose that  $c= 4\alpha\lambda_n\sqrt{n}/\sqrt{\lambda_1}$ 
and $\alpha\geq 1$. Then with density function proportional to
the following $F$ in all of ${\bf R}^n$,
the probability inside the orthant is at least $0.25$.
\begin{equation}
\label{equation:gen:f} 
F(x) = e^{-xAx^T - c\varphi(x)}
\end{equation}
\end{proposition}

\SubSection{Step size, local conductance and conductance}
\label{sec:gen:local}
From now on, we sample in the whole $R^n$ according to the following function:
$$F(x) =
e^{-xAx^T-(4\alpha\lambda_n\sqrt{n}/\sqrt{\lambda_1})\varphi(x)}.$$The
step size used is 
$$\delta ={{\sqrt{\lambda_1}}\over{\gamma\lambda_n\sqrt{n}\ln n}}$$
where $\gamma\geq 2$ is a positive real.  It can be proved that
with this step size, the local conductance is at least
$$({1\over 2} - {1\over{4\gamma\sqrt{n} \ln n}})\omega$$
where $\omega$ is a positive real. The proof involves many details, and thus
will be presented in the final paper. 

We then optimize the parameters as in the other problem. Here is the algorithm with the
parameters chosen (nearly) optimally.

{\bf Algorithm } Start at $w_0$ chosen according to the uniform density on the unit
ball intersected with the positive orthant. Then do the Metropolis random walk in all of
${\bf R}^n$ according to the function 
$$F(x) = e^{-x^T Ax -5.5 {\sqrt n \lambda_n\over \sqrt\lambda_1}\varphi (x) }$$
and with step size 
$$\delta = {\sqrt\lambda_1\over 16 \lambda_n \sqrt n\log n}.$$
We do the random walk for a number of steps dictated by the desired $M$ distance
to steady state as determined by Proposition \ref{start2}, Theorem \ref{theorem:gencond}
, Theorem \ref{theorem:mixing} and Corollary \ref{subsequent}.

We prove the following about this walk. 

\begin{proposition}\label{start2}
$$M(Q_0,Q)\leq 1+ 2^n {e^{\lambda_n}\over\lambda_1} \pi^{n/2} n!.$$
\end{proposition}

\begin{theorem}
\label{theorem:gencond}
The conductance of this chain is at least
\begin{equation}
\label{equation:gen:cond}
{1\over 900 \Lambda n\log n}.
\end{equation}
\end{theorem}

Just like the standard normal case that we discussed in the last
section, we can give bounds for the first sample mixing time and
the subsequent sample mixing time. The first
sample takes $O^*(n^3\Lambda^3)$ steps. 
Each subsequent (nearly independent) sample
takes number of steps 
$$O^*(n^2\Lambda^2).$$

\Section{From Sampling to Integration}

In this section we show how the ability to sample from ${\bf R}^n_+$
according to functions of the form $e^{- x^TAx}$ helps us estimate the
integrals of such functions over the orthant. 

The connection between sampling and integration is in a sense a
natural and old one, but the point here is to describe a scheme
which given a polynomial time sampling algorithm, can estimate the
integral in polynomial time. 

For $p \in [0\ \ 1]$, define 
$$ f(p) = \int_{ {\bf R}_+^n } e^{-p x^TAx - (1-p) |x|^2 } \ dx
  = \int_{ {\bf R}^n_+}e^{-x^TMx} \ dx ,$$
where $M = p A + (1-p)I$ is a positive definite matrix. 
$f(0) $ is known and we are interested in $f(1)$.

Let $X$ be a ($n-$vector) random variable with values in ${\bf R}^n_+$
and with density 
$$e^{-p x^TAx - (1-p) |x|^2 } / f(p).$$

Let $\nu $  be a positive real number to be chosen later. 
Let $Y= e^{-\nu X^TAX +\nu |X|^2 }$. Then, clearly, we have
$$E(Y) = {f(p+\nu )\over f(p) }.$$

We will use this to estimate ${f(p+\nu )\over f(p) }$ and we will then 
find $f(1)/f(0)$ as the telescoping product of several such ratios. First,
let us bound the number of samples of $X$ needed.

We have 
$$\var (Y)\leq E(Y^2) = {f(p+2\nu )\over f(p) }.$$
So, we have 
$${\var (Y) \over (E(Y))^2} \leq {f(p)f(p+2\nu )\over f(p+\nu )f(p+\nu )}.$$
We will use the following Lemma to bound this ratio.

{\bf Lemma } Assuming $\lambda_1 = 1$ after rescaling if
necessary, we have 
$$\left| {d(\log f)\over dp}\right| \leq {n\over 2}
                \Lambda\qquad \forall p.$$

\pr We see by differentiating that 
$${df\over dp} = \int_{ {\bf R}^n_+ } [-x^TAx + |x|^2 ]e^{-x^TM x} \ dx,$$
Since 
$x^TAx\leq \lambda_n |x|^2$, we have 
$$\left| {df\over dp} \right|
  \leq \lambda_n
       \int_{ {\bf R}^n_+ } |x|^2 e^{-x^TM x} \ dx.$$
Define a function $g : {\bf R}_+\rightarrow {\bf R}$ by 
$$g(t) = \int_ { {\bf R}^n_+} e^{ - t x^T Mx}\ dx.$$
Then with the change of variables $y_i = \sqrt t x_i$, we see that
$g(t) = t^{-n/2} g(1)$. Differentiating we get $dg/dt = - n/ 2$ at $t=1$.
But differentiating under the integral we get (noting that the eigenvalues of
$M$ are all also at least 1), 
$$-{dg\over dt}  = \int_ { {\bf R}^n_+} (x^T M x) e^{ - t x^T Mx}\ dx
\geq  \int_ { {\bf R}^n_+} |x|^2 e^{ - t x^T Mx}\ dx.    $$
So we get $$|{df\over dp}|\leq \Lambda f(p),$$ from which the lemma
follows.

From the Lemma we get 
$$\left| \log \left( {f(p+2\nu )\over f(p+\nu )} \right)\right| 
       = \int_{q=p+\nu }^{p+2\nu } \left| {d (\log f)\over dp}\right| dq\leq 
           n\nu \Lambda / 2.$$
The same bound holds for $\left| \log \left( {f(p+\nu )\over f(p)} \right)
\right|$. So, choosing
$$\nu = {1\over n\Lambda},$$
we see that 
$${\var (Y) \over (E(Y))^2}\quad \leq 3 .$$
Thus then the ratio $f(p+\nu )/ f(p)$ can be estimated to 
relative error at most $\delta $ with $O(1/\delta ^2)$ samples.
Let $m=n\Lambda$. 

We use the notation introduced in this
section to describe the algorithm, stating the number of 
samples used etc.

{\bf Algorithm to estimate the integral}

\begin{itemize}

\item Given $A; 0\leq \varepsilon\leq 1/8$ and $\delta >0$, the algorithm 
will find a real number $T$ such that with
probability at least $1-\delta $, we have :
$$Te^{-\varepsilon }\leq \int_ { {\bf R}^n_+ }e^{-x^TAx } \leq Te^{\varepsilon }.$$

\item For $p=0,\nu ,2\nu ,\ldots $ (upto 1), we do the following :
\begin{itemize}
     \item Draw $s=20m/(\varepsilon ^2\delta )$ (pairwise nearly independent)
           samples $X_1,X_2,\ldots X_s$ (nearly) according to the density
               $$e^{-p x^TAx - (1-p) |x|^2 } / f(p).$$
      \item Let $$T_p =  {1\over s} \sum_{i=1}^s e^{-\nu |X_i^T AX_i| +\nu |X_i|^2 }.$$
\end{itemize}
\item Return the product of all the $T_p$ 's as $T$.
\end{itemize}

{\bf Remark } : The samples drawn for a particular value of $p$ will be totally independent
of samples drawn for other values of $p$.
           
To give the essential idea of the analysis of the algorithm, we will assume here that

(i) $E(T_p) = {f(p+\nu )\over f(p)}$ and

(ii) $\var (T_p) \leq \varepsilon^2\delta /(4m)$.

Both these are valid if the $s$ samples used for finding $T_p$ are pairwise completely
independent (i.e., 0-independent) 
and each exactly according to the desired density. This is not quite true,
but the errors cause only minor technical problems and the final paper will contain
the full analysis. Under the assumptions, we have first that
$$\prob (T_p\leq E(T_p)/\sqrt 2)\leq \delta/(16m).$$
So with probability at least $1-(\delta /16)$, we have that for every $p$ : $T_p\geq E(T_p)/\sqrt 2$. This technical fact is only used in the second case below. 


Now, consider the random variable 
$$t=\log {T\over \prod_p (E(T_p))} = \sum_{p} \log {T_p\over E(T_p)}.$$
First consider the case when $t\geq 0$. 
Since $\log x\leq x-1$ for all real $x$, we have
$$t = 
\sum_{p} \log {T_p\over E(T_p)}\leq \sum_p {T_p-E(T_p)\over E(T_p) }
.$$
The random variables ${T_p-E(T_p)\over E(T_p) }\; ;\; p = 0,\nu ,\ldots $ are independent,
so the variance of the sum above is just the sum of the variances. Thus we have by Markov
inequality that 
$$\prob (t\geq \varepsilon  ) \leq \delta /4.$$
The second case when $t<0$ is a little more complicated. In this case, we use the
inequality $\log x\geq (x-1)-(x-1)^2$ which is valid for all $x\geq 1/\sqrt 2$. We will defer the
full argument to the final paper. 

We saw that we can produce the first sample for each $p$ in $O^*(\Lambda^2 n^3)$ steps
and subsequent samples in $O^*(\Lambda^2 n^2)$ steps per sample.
So the total number of steps needed to estimate the integral to relative error
$\epsilon$ is at most 
$$O^*(\Lambda^4 n^4).$$

\Section{Open Problems}

Several interesting open problems remain. The first will be to improve the number of
steps needed. In this connection, it should be possible to use the results of \cite{kls},
especially the notion of proper moves. A proper move of the Metropolis random walk is what the
name indicates - it is a move where we move from $x$ to $y\not= x$. If the local conductance
is low at a point, we have to make many moves before we make a proper move. But in certain
situations, like sampling from the orthant, it may be easy to run a different chain whose
transition probability density function is : 
$$P_{xy} = {1\over \int_{ {\bf R}^n_+\cap (x+\delta B)} F} \min \left( 1, {F(y)\over F(x)}\right).$$
In other words, we may be
able to make proper moves with the correct probabilities without having to ``try'' out 
(possibly) improper moves. In \cite{kls}, better convergence bounds are proved in such cases and 
it will be interesting to exploit this.
There is a hope that for such an improvement for the case of the convex set containing the
unit ball (as in section \ref{sec:std}) may help improve the volume algorithm. 

Our arguments do not carry over for the case when the orthant is shifted - i.e., when the 
region is of the form : $\{ x : x_i\geq p_i\}$ where not all the $p_i$ are zero. The main
difficulty is the cumbersome penalty function and the argument that the positive orthant
has sizable (here at least 1/4) measure. Perhaps this can be overcome and a better penalty 
function or a better argument can be found. 

It is also possible that these methods here extend to other log-concave densities. 

{\bf Acknowledgements } : We thank Persi Diaconis for bringing the problem to
our attention. Thanks to Satish Iyengar for helpful discussions. 

\begin{thebibliography}{99}
\bibitem{aldous} 
D. Aldous,
``On the Markov chain simulation method for uniform combinatorial distributions
and simulated annealing,''
{\em Probability in the Engineering and Informational Sciences 1}, pp. 33-46,
1987.
\bibitem
D. Aldous, ``Random Walks on Finite Groups and Rapidly
Mixing Markov Chains'', {\it S\'{e}minaire de Probabilit\'{e}s XVII}, 
Springer Lecture Notes in Mathematics 986, pp. 243-297, 1981$/$82.
%\bibitem{ad} D. Aldous and P. Diaconis, ``Shuffling Cards and Stopping
Times'',
%{\it American Mathematical Monthly 93}, pp. 333-348, 1986.
\bibitem{ak} D. Applegate and R. Kannan, ``Sampling and integration of near
log-concave
functions'', {\it Proceedings of the 23rd ACM Symposium on Theory of
Computing}, IEEE, 1991.

\bibitem{bro} A. Z. Broder, ``How hard is it to marry at random? (On the
approximation
of the permanent)'', {\it Proceedings of the 18th ACM Symposium
on Theory of Computing}, pp. 50-88, 1986.

\bibitem{ds} P. Diaconis and Saloff-Coste, ``What do we know about Metropolis
Algorithm'', 
{\it Symposium on the Theory of Computing}, ACM (1995). 

\bibitem{df} M. Dyer and A. Frieze, ``Computing the volume of convex bodies: a
case
where randomness provably helps'', {\it Probabilistic Combinatorics and Its
Applications} {\bf 44}, 1991.

\bibitem{dfk} M. Dyer, A. Frieze and R. Kannan, ``A random polynomial time
algorithm for approximating the volume of convex bodies'', {\it Journal
of the ACM} {\bf 38}, 1991.

%\bibitem{mount} M. Dyer, R. Kannan and J. Mount, ``Sampling Contingency
Tables'',
%{\it unpublished manuscript}, 1995.

\bibitem{evans} M. Evans, and T. Swartz, ``Monte Carlo Computation of Some Multivariate 
Normal Probabilities'', 
{\it J. Statist. Comput. Simul.}, Vol. 30, pp. 117-128, 1988.

\bibitem{polson} A. Frieze, R. Kannan and N. Polson, ``Sampling from
Log-concave
Distributions'', {\it Annals of Applied Probability}, Vol. 4, No. 3, (1994) pp 812-837. 

\bibitem{satish} S. Iyengar, Evaluation of Normal Probabilities of Symmetric
 Regions, {\it SIAM J. Sci. Stat. Comput.}, Vol.9, No.3, May 1988.

%\bibitem{jerrum} M. Jerrum, ``A Very Simple Algorithm for Estimating
%the Number
%of $k-$colourings of a Low-degree Graph'' {\it unpublished manuscript}, 1995.

\bibitem{js} M. Jerrum and A. Sinclair, ``Approximating the permanent'', {\it 
SIAM J. Comp.} {\bf 18}, 1989.

\bibitem{k} R. Kannan, Markov Chains and Polynomial Time Algorithms, {\it
Proceedings of the 35th Annual Symposium on Foundations of Computer Science},
pp.656-671,
1994.

\bibitem{kls} R. Kannan, L. Lov\'{a}sz and M. Simonovits, {\it Random walks
and an $O^*(n^5)$ volume algorithm for convex sets} Preprint (1996).

\bibitem{ls} L. Lov\'{a}sz and M. Simonovits, ``Random walks in a convex
body and
an 
improved volume algorithm'', {\it Random Structures and Algorithm} {\bf
4}, 1993.

\bibitem{mar} G. Martynov, Evaluation of the Normal Distribution Function, 
{\it J. Soviet Math.}, 77 (1980), pp. 1857-1875.

\bibitem{moran} P. Moran, ``Calculation of Multivariate Normal 
Probabilities---Another Special Case'', {\it Austral. J. Statist.}, 
27 (1), pp. 60-67, 1985.
\end{thebibliography}
\end{document}


