\documentclass[10pt,twocolumn,letterpaper]{article}

\usepackage{cvpr}
\usepackage{times}
\usepackage{epsfig}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{cancel}
\usepackage{algorithm}
\usepackage{algorithmic}

% Include other packages here, before hyperref.

% If you comment hyperref and then uncomment it, you should delete
% egpaper.aux before re-running latex.  (Or just hit 'q' on the first latex
% run, let it finish, and you should be clear).
\usepackage[pagebackref=true,breaklinks=true,letterpaper=true,colorlinks,bookmarks=false]{hyperref}


% \cvprfinalcopy % *** Uncomment this line for the final submission

\def\cvprPaperID{****} % *** Enter the CVPR Paper ID here
\def\httilde{\mbox{\tt\raisebox{-.5ex}{\symbol{126}}}}

% Pages are numbered in submission mode, and unnumbered in camera-ready
\ifcvprfinal\pagestyle{empty}\fi
\begin{document}

%%%%%%%%% TITLE
\title{Particle Filter Estimation of Posterior Densities}

\author{First Author\\
Institution1\\
Institution1 address\\
{\tt\small firstauthor@i1.org}
% For a paper whose authors are all at the same institution,
% omit the following lines up until the closing ``}''.
% Additional authors and addresses can be added with ``\and'',
% just like the second author.
% To save space, use either the email address or home page, not both
\and
Second Author\\
Institution2\\
First line of institution2 address\\
{\small\url{http://www.author.org/~second}}
}

\maketitle
% \thispagestyle{empty} % *** Uncomment this line for the final submission

%%%%%%%%% ABSTRACT
\begin{abstract}
\end{abstract}

%%%%%%%%% BODY TEXT





\section{Background}\label{sec:background}

Our goal is to compute a posterior distribution $p(X_1, \dots, X_m \ | \ Z )$,
where $(X_1, \dots, X_m)$ is a vector of random variables (RVs)
and $Z=\{ z_1, \dots, z_m \}$ is a set of observations.
This will allow us to find value assignments $X_t=\widehat{x}_t$ for $t=1, \dots, m$
to RVs that maximize this posterior:
\begin{align}
& \widehat{x}_{1:m} = \underset{x_{1:m}}{\operatorname{argmax}} \ p(X_1=x_{1}, \dots, X_m=x_{m} \ | \ Z),
\label{eqG2}
\end{align}
where $x_{1:m}= (x_1, \dots, x_m) \in \mathcal{X}^m$ is a state space vector.
As it is commonly the case, we will abbreviate
$$p(X_1=x_{1}, \dots, X_m=x_{m} \ | \ Z)=p(x_{1:m} \ | \ Z).$$
We will achieve our goal by approximating the posterior distribution with a final number of samples
in the framework of Bayesian Importance Sampling (BIS).
Since it is usually difficult to draw samples from the pdf $p(x_{1:m}|Z)$,
we will draw samples $x_{1:m}^{(i)} \sim q(x_{1:m}|Z)$ for $i=1, \dots, N$
from a proposal pdf $q$, form which samples are easily generated.
Then approximation to the density $p$ is given by
\begin{align}\label{eq:PFn1}
p(x_{1:m}|Z)& \approx \sum_{i=1}^N w^{(i)} \delta(x_{1:m} - x_{1:m}^{(i)}),
\end{align}
where $\delta$ is the Dirac delta function and
\begin{align}\label{eq:PFn2}
w^{(i)} = \frac{p(x_{1:m}^{(i)}|Z)}{q(x_{1:m}^{(i)}|Z)}
\end{align}
are normalized weights (so that they sum to one).

Since it is still hard to draw samples from $q$ due to high dimensionality of $x_{1:m}$,
Sequential Importance Sampling (SIS) is usually utilized.
Following the order of dimensions in the vector of RVs $X=( X_1, \dots, X_t )$
samples are generated
\begin{align}\label{eq:sampling}
x_{t}^{(i)} \sim q(x_t|x_{1:t-1}, Z) = q(x_t|x_{1:t-1}, z_{1:t})
\end{align}
for $t = 2, \dots m$, and the particles are built sequentially $x_{1:t}^{(i)} =(x_{1:t-1}^{(i)}, x_{t}^{(i)})$ for $i=1, \dots, N$.
Since $q$ factorizes as
\begin{align}\label{eq:factor1}
& q(x_{1:m} | Z) = q(x_1 | Z) \prod_{t=2}^m q(x_t | x_{1:t-1}, Z),
\end{align}
we obtain that $x_{1:m}^{(i)} \sim q(x_{1:m} | Z)$.


The weights are recursively updated according to
\begin{align}\label{eq:PFn3}
w(x_{1:t}^{(i)}) = \frac{p(z_t|x_{1:t}^{(i)},z_{1:t-1}) p(x_{t}^{(i)}|x_{1:t-1}^{(i)})}
{q(x_{t}^{(i)}|x_{1:t-1}^{(i)},z_{1:t})} \ w(x_{1:t-1}^{(i)}) .
\end{align}
We show now that by recursively updating the weights according to \eqref{eq:PFn3} for $t=2, \dots, m$,
the weight $w(x_{1:m}^{(i)})$ of particle $(i)$ is equal to $w^{(i)}$.
\begin{align}\label{eq:PFn3.5}
&w(x_{1:t}^{(i)}) = \frac{p(x_{1:t}^{(i)}|Z)}{q(x_{1:t}^{(i)}|Z)} =
\frac{ p(x^{(i)}_{t}|x_{1:t-1}^{(i)},Z) \ p(x_{1:t-1}^{(i)}|Z) }
{ q(x^{(i)}_{t}|x_{1:t-1}^{(i)},Z) \ q(x_{1:t-1}^{(i)}|Z) } \nonumber \\
&=\frac{ p(x^{(i)}_{t}|x_{1:t-1}^{(i)},Z) }{ q(x^{(i)}_{t}|x_{1:t-1}^{(i)},Z) } \ w(x_{1:t-1}^{(i)}) \\
&= \frac{ p(Z|x_{1:t}^{(i)}) } { p(Z|x_{1:t-1}^{(i)})  }  \
\frac{ p(x^{(i)}_{t}|x_{1:t-1}^{(i)}) } {  q(x^{(i)}_{t}|x_{1:t-1}^{(i)},Z) } \ w(x_{1:t-1}^{(i)}) \label{eq:PFn3.6}
\end{align}
Eq.~\eqref{eq:PFn3.6} follows from \eqref{eq:PFn3.5} by
Bayes rule interchanging $x^{(i)}_{t}$ and $Z$ in $p(x^{(i)}_{t}|x_{1:t-1}^{(i)},Z)$.
It remains to show that the first fraction in \eqref{eq:PFn3.6} is equal to
$p(z_t|x_{1:t}^{(i)},z_{1:t-1})$.
This is true, since by keeping in mind that $Z=z_{1:t}$ and by factorizing we obtain
\begin{align}\label{eq:PFn3.8}
& p(Z | x_{1:t}^{(i)}) = p(z_t|x_{1:t}^{(i)},z_{1:t-1})  \ p(z_{1:t-1} |x_{1:t-1}^{(i)}) \nonumber \\
& = p(z_t|x_{1:t}^{(i)},z_{1:t-1})  \ p(Z |x_{1:t-1}^{(i)})
\end{align}


We have just derived the SIS theorem:

\noindent {\bf Theorem. }
By sampling particles according to \eqref{eq:sampling} and weighting them according \eqref{eq:PFn3},
we obtain a set of weighted samples form $p(x_{1:m}|Z)$ once $t=m$.
Consequently, we can approximate $p(x_{1:m}|Z)$ with any precision
if the number of particles $N$ is sufficiently large. Thus, we can write
\begin{align}\label{eq:PFSIS}
p(x_{1:m}|Z)& \approx \sum_{i=1}^N w(x_{1:t}^{(i)}) \delta(x_{1:m} - x_{1:m}^{(i)}),
\end{align}

A common assumption underlying the equations \eqref{eq:sampling} and \eqref{eq:PFn3} is that there exist two known functions
$f$ and $h$ such that
\begin{align}\label{eq:PFa1}
x_{t} &= f(x_{t-1}) + u_t\\
z_{t} &= h(x_{t}) + v_t \label{eq:PFa2}
\end{align}
where both $u_t$ and $v_t$ are mutually independent and identically distributed sequences with known probability density functions (pdfs). Often they are assumed to be Gaussians that model state prediction noise $u_t$
and observation noise $v_t$.
(The assumption in \eqref{eq:PFa1} can be replied by a weaker one: $x_{t} = f(x_{1:t-1}) + u_t$,
i.e., that we can determine $x_{t}$ if we know all previous states $x_{1:t-1}$.
For the simplicity of presentation we will use \eqref{eq:PFa1}.


The equation (\ref{eq:PFn3}) can be simplified by making a common assumption that
$q(x_{t}^{(i)}|x_{1:t-1}^{(i)},z_{1:t}) = p(x_{t}^{(i)}|x_{1:t-1}^{(i)})$, which yields
\begin{align}\label{eq:PFn4}
w(x_{1:t}^{(i)}) = w(x_{1:t-1}^{(i)}) p(z_t|x_{1:t}^{(i)},z_{1:t-1}),
\end{align}
and, consequently, the samples are generated from
\begin{align}\label{eq:sampling2}
x_{t}^{(i)} \sim p(x|x_{1:t-1}^{(i)}),
\end{align}
i.e., the current observations $z_{1:t-1}$ have no direct influence on the samples.
The influence of observations on the samples is introduced when particles are resampled according to their weights, which is a common step in PF algorithms,
since particle weights depend on observations.

As a summery of this section, we outline a {\bf standard PF algorithm}.
For every time step $t=1, \dots, m$ and for every particle $i =
1, \dots, N$ execute the following three steps:

\noindent 1) \textbf{Importance sampling / proposal:}
Sample followers of particle $(i)$ according to \eqref{eq:sampling} and set
$x_{1:t}^{(i)} = (x_{1:t-1}^{(i)}, x_{t}^{(i)})$.

\noindent 2) \textbf{Importance weighting / evaluation:}
An importance weight is assigned
to each particle $x_{1:t}^{(i)}$ according to \eqref{eq:PFn3}.

\noindent 3) \textbf{Resampling:}
Sample with replacement $N$ new particles form the current set of particles according to the weights.
We obtain a set of new particles $x_{1:t}^{(i)}$ for $i=1, \dots, N$ all with weight of $1/N$.
This procedure is called Sampling Importance Resampling (SIR) approach <cite>.


\section{New Approach}\label{sec:new}

We illustrate our key ideas on an example of multi robot localization.
A team of $m$ robots obtained $m$ observations $Z=\{ z_1, \dots, z_m \}$ by exploring their environment,
where observation $z_t$ comes from robot $t$.
For example, if the robots are equipped with laser range scanners,
then each $z_t$ could be a vector of laser range readings representing distances to the closest obstacles
from the robot $t$.
Each RV $X_t$ describes robot poses,
i.e., its values $x_t$ represent coordinates and the heading direction of the robot $t$.
Our goal is to determine the state vector of robot poses $(x_1, \dots, x_m )$ in a given top view map of the environment that maximizes the posterior distribution $p(X_1, \dots, X_m | Z)$.
In other words, $(x_1, \dots, x_m )$ is a vector of most likely robot poses given the measurements $Z$.

It is possible to apply the classical PF robot localization by utilizing the order of
the observations $Z=\{ z_1, \dots, z_m \}$, which follows the numbering of the robots.
Hence the follower for each particle $(i)$ is determined by importance sampling from the proposal distribution, i.e.,
sample $x_t^{(i)} \sim p(x|x_{1:t-1}^{(i)})$ and the particle weight is updated based on the evaluation of the observation $z_t$
according to the recursive formula in Eq.~\ref{eq:PFn4}.
However, by doing so, we would have selected an arbitrary order,
and in particular, if the robot localization task fails, it could be due to the selected order.
Would we have selected a different order, the localization task could have been successful.
Moreover, the observations $Z=\{ z_1, \dots, z_m \}$ are collected simultaneously at the same time.
i.e., after each robot completed its exploration.
Consequently, there is no reason to favor any particular order without utilizing further information.

In the proposed approach, the order of the observations is not predetermined,
in particular, we do not follow the order of indices of the observations in $Z$.
Our key idea is to utilize the PF framework to determine the most informative order of the observations.
This way we are able to simultaneously find the most informative order and
to utilize the observations in the order of their informativeness.
Intuitively, it make sense, for example,
if the first robot took its laser readings in the middle of a along a long corridor and
the second robot at the entrance to a room with many distinctive features,
then our approach will first process the laser readings obtained by the second robot,
since they are more informative.


We stress that the SIS in Eq.~\ref{eq:sampling}
and particle evaluation in Eq.~\ref{eq:PFn3} utilize the sequential order of the RVs
reflected in the order of dimensions in the state space $(x_1, \dots, x_m )$.
In many applications, this order is determined naturally by the time stamp of the observations,
e.g., a single robot is collecting laser measurements at consecutive time points,
in which case $x_t$ denotes the robot pose at time $t$.
The goal of our work is to extend SIS to applications in which there is no natural order of observations like the case of multi robot localization.

The key idea of the proposed approach is not to utilize the fix order of the dimensions,
but instead compute the best possible order of the dimensions $(x_{i_1}, \dots, x_{i_m})$ (or equivalently RVs)
so that the corresponding sequence of observations $Z=(z_{i_1}, \dots, z_{i_m})$ is most informative.
%This leads to faster convergence with a small number of particles $N$.
For this we extend the underlying assumption \eqref{eq:PFa1} to the existence of a sequence of functions
\begin{align}\label{eq:PFa3}
x_{s} &= f_{s|k} (x_{{k}}) + u_{s,k}  \ \ \mathrm{for all} \ \  s,k \in \{1, \dots, m \},
\end{align}
which means that all RVs depend on each other, and each function $f_{s|k}$
can determine the value of variable $s$ for any given value $x_k$ of variable $k$.
In our multi robot mapping example, $f_{s|k}$ allows to determine the pose of robot $s$
when we know the pose of robot $k$.
Consequently, \eqref{eq:PFa3} means that each robot knows the relative pose of the other robots.
In comparison, in the standard assumption \eqref{eq:PFa1}, we have $f=f_{t|t-1}$,
i.e., we only can determine the value of the next state.

We observe that now the sequence of states visited before time $t$ is not a sequence of consecutive numbers
$(1, \dots, t-1)$ but any subsequence $(i_1, \dots, i_{t-1})$ formed by $t-1$ different numbers in
$\{1, \dots, m \}$, and
the function $f_{s|k}$ allows us to
determine not only the value of the next state but values of all remaining states.
Due to noise factor $u_{s,k}$, each value $x_s$ is an estimate or prediction of a true unknown value.
Since a given state $x_s$ can be determined based on all the states in the current sequence,
we can combine all the predictions $f_{s|i_k}(x_{i_k}) + u_{s,i_k}$ for $ k =1, \dots, t-1$
and improve the accuracy of the estimation of the state $x_s$.




The proposed sampling is as follows.
At iteration $t$, the assumption \eqref{eq:PFa3} allows us to generate $m-t+1$ samples from
\begin{align}\label{eq:sampling2.5}
x_{i_s}^{(i)} \sim p_{i_s|i_{t-1}}(x|x_{i_{1:t-1}}^{(i)}) \ \ \mathrm{for} \ \  s \in \{t, \dots, m \}.
\end{align}
Hence at iteration $t$ particle $(i)$ has $m-t+1$ followers.
For example, if $m=5$ and $x_{i_{1:3}}^{(i)} = x_{1, 5, 3}^{(i)}$,
then particle $(i)$ will have two followers $x_{2}^{(i)}$ and $x_{4}^{(i)}$.
With reference to our multi robot example, when we determined the locations of robots $1, 3, 5$,
we consider the two remaining robots 2 and 4 as the next robot whose pose we want to determine by particle (i).
The poses of robots 2 and 4 are determined with reference to the pose of robot 3
as the last element of the sequence 1, 5, 3.
Of course, we repeat an analogous process for each particle (i) for $i=1, \dots, N$,
where $N$ is the number of particles.

In contrast, in the standard application of rule \eqref{eq:sampling2},
at each iteration $t$ particle $(i)$ has one follower.
Even when sometimes each particle $(i)$ has many followers, all followers are in the same dimension,
which means that we only determine possible locations of say robot 2
by the followers and do not consider locations of robot 4 for particle (i), since a strict order of the state dimensions is followed in the classical setting.





\section{Particle Filter with Static Observations}\label{sec:PF}
Our derivation is analog to the
PF derivation, but it differs fundamentally, since unlike the
standard PF framework, the observations $Z$ do not arrive
sequentially, but are available at once.
To simplify the notation we replace the double indexing of the state variables $x^{(i)}_{i_s}$
with a bijection (onto and one-to-one function)
$<\cdot>^{(i)}: \{1, \dots, m \} \to \{1, \dots, m \}.$
Although we may have a different bijection for each particle,
we will drop the index $(i)$ from $<1:t>^{(i)}$, since the state variables carry the particle index.
For example, we will denote $(x^{(i)}_{i_{1}},  x^{(i)}_{i_{2}}, x^{(i)}_{i_{3}}) = x^{(i)}_{i_{1:3}}$ as $x^{(i)}_{<1:3>}$.


We derive now the recursive weight update formula for the static observations $Z$.
For every $t$ from 2 to $m$, we have
\begin{align}\label{EQ:pf3}
&w(x_{<1:t>}^{(i)})
=\frac{ p(x_{<1:t>}^{(i)}  \ | \ Z) }{ q(x_{<1:t>}^{(i)}  \ | \ Z) } \nonumber \\
&=\frac{ p(x^{(i)}_{<t>}|x_{<1:t-1>}^{(i)},Z) \ p(x_{<1:t-1>}^{(i)}|Z) }
{ q(x^{(i)}_{<t>}|x_{<1:t-1>}^{(i)},Z) \ q(x_{<1:t-1>}^{(i)}|Z) } \nonumber \\
&=\frac{ p(x^{(i)}_{<t>}|x_{<1:t-1>}^{(i)},Z) }{ q(x^{(i)}_{<t>}|x_{<1:t-1>}^{(i)},Z) } \ w(x_{<1:t-1>}^{(i)}) \nonumber \\
&= \frac{ p(Z|x_{<1:t>}^{(i)}) p(x^{(i)}_{<t>}|x_{<1:t-1>}^{(i)}) }
{ p(Z|x_{<1:t-1>}^{(i)}) q(x^{(i)}_{<t>}|x_{<1:t-1>}^{(i)},Z) } \ w(x_{<1:t-1>}^{(i)})
\end{align}
To obtain the last equation, we apply Bayes rule to decompose
$p(x^{(i)}_{<t>}|x_{<1:t-1>}^{(i)},Z)$ that interchanges $x^{(i)}_{<t>}$ and $Z$.

As it is often the case in PF applications, we assume that
$q(x^{(i)}_{<t>}|x_{<1:t-1>}^{(i)},Z) = p(x^{(i)}_{<t>}|x_{<1:t-1>}^{(i)})$.
Using this simple exploration based proposal the weight recursion in
\eqref{EQ:pf3} becomes:
\begin{align}\label{eqn:recur-wgt}
&w(x_{<1:t>}^{(i)})=
&=w(x_{<1:t-1>}^{(i)}) \ \frac{p(Z|x_{<1:t>}^{(i)})}{p(Z|x_{<1:t-1>}^{(i)})}
\end{align}

By recursive substitution of weights in \eqref{eqn:recur-wgt},
i.e., by applying \eqref{eqn:recur-wgt} to $w(x_{<1:t-1>}^{(i)}), w(x_{<1:t-2>}^{(i)}), \dots, w(x_{<1:2>}^{(i)})$,
we obtain
\begin{align}\label{eqn:ind-wgt}
&w(x_{<1:t>}^{(i)})
=w(x_{<1:t-2>}^{(i)}) \ \frac{\cancel{p(Z|x_{<1:t-1>}^{(i)})}}{p(Z|x_{<1:t-2>}^{(i)})} \ \frac{p(Z|x_{<1:t>}^{(i)})}{\cancel{p(Z|x_{<1:t-1>}^{(i)})}} \nonumber \\
&= \dots =w(x_{<1>}^{(i)}) \ \frac{p(Z|x_{<1:t>}^{(i)})}{p(Z|x_{<1>}^{(i)})}
\end{align}
Under the assumption that all particles have the same initial weight $w(x_{<1>}^{(i)})$
and the same initial observation probability $p(Z|x_{<1>}^{(i)})$ for $i = 1, \dots, N$, we obtain
\begin{align}\label{eqn:partweight}
&w(x_{<1:t>}^{(i)})= p(Z|x_{<1:t>}^{(i)})
\end{align}


Since for $t=m$ we have $w(x_{<1:m>}^{(i)})=w(x_{1:m}^{(i)})$,
we obtain that the weights computed by the recursive formulas \eqref{eqn:recur-wgt} -
\eqref{eqn:partweight} are equal to the weights in Eq.~\eqref{eq:PFn2}.





For comparison, the corresponding weight update in the standard PF
framework (\cite{ThrunBook05}) is
\begin{align}\label{eq:PFstandard}
w(x_{1:t}^{(i)})&=
w(x_{1:t-1}^{(i)}) \ p(z_t | x_{1:t-1}^{(i)}, x_t),
\end{align}
where $z_t$ denotes the new observations obtained at time $t$.
Because our observations $Z$ do not have any natural order,
$Z$ cannot be expressed as a sequence of observations.
We do not make any Markov assumption in the proposed formula \eqref{eqn:partweight},
i.e., the new state $x^{(i)}_{<t>}$ is dependent on all previous states $x_{<1:t-1>}^{(i)}$ for each particle $(i)$.




We outline now the proposed PF algorithm with static observations (PFSO).
For every time step $t=1, \dots, m$ and for every particle $i =
1, \dots, N$ execute the following three steps:

\noindent 1) \textbf{Importance sampling / proposal:}
Sample followers of particle $(i)$ for $s \in \{1, \dots, m \} \setminus <1:t-1>$
\begin{align}\label{eq:sampling2.5.1}
x_{s}^{(i)} \sim p_{s|<t-1>}(x|x_{<1:t-1>}^{(i)}) %\ \ \mathrm{for} \ \
\end{align}
and set $x_{<1:t-1>,s}^{(i)} = (x_{<1:t-1>}^{(i)}, x_{s}^{(i)})$.
This step is possible by assumption \eqref{eq:PFa3}.
We sample one follower for each RV $X_s$ whose dimension index $s$ is not in $<1:t-1>$.
Consequently, we have $m-t+1$ followers at step $t$.
In the first iteration ($t=1$) we generate $m$ samples
\begin{align}\label{eq:sampling2.7}
x_{<1>}^{(i)}=x_{s}^{(i)} \sim p_s(x) \ \ \mathrm{for} \ \  s \in \{1, \dots, m \}.
\end{align}

\noindent 2) \textbf{Importance weighting/evaluation:}
An individual importance weight is assigned
to each follower of each particle according to Eq.~\ref{eqn:partweight}.

\noindent 3) \textbf{Resampling:}
At the sampling step we sample more followers than the number of particles.
Thus we have a larger
set of particles $x_{<1:t-1>,s}^{(i)}$ for $i=1, \dots, N$ and
$$s \in \{1, \dots, m \} \setminus <1:t-1>$$
from which we sub-sample $N$ particles and assign qual weights to all of them as in the standard
Sampling Importance Resampling (SIR) approach <cite>.
We obtain a set of new particles $x_{<1:t>}^{(i)}$ for $i=1, \dots, N$.
The resampling is not performed in the last step, i.e., when $t=m$. \\


\noindent \textbf{Theorem.} \
The PFSO algorithm computes an approximation of Eq.~\ref{eq:PFn1},
i.e., for $t=m$ the particles $x_{<1:t>}^{(i)}$ with $i=1, \dots, N$
provide an approximation to the posterior $p(x_{1:m}|Z)$ for sufficiently large $N$, i.e.,
\begin{align}\label{eq:PFn1.1}
p(x_{1:m}|Z)& \approx \sum_{i=1}^N w(x_{1:m}^{(i)}) \delta(x_{1:m} - x_{1:m}^{(i)}).
\end{align}
\noindent  \textbf{Proof.} \
Our proof is based on Eq.~\eqref{eq:PFn1.1} with weights defined in Eq.~\eqref{eq:PFn2}.
We observe that the weights computed by the recursive formulas
\eqref{eqn:recur-wgt} - \eqref{eqn:partweight} are equal to the weights in Eq.~\eqref{eq:PFn2},
since for $t=m$ we have $w(x_{<1:m>}^{(i)})=w(x_{1:m}^{(i)})$.

The Sampling Importance Resampling (SIR) replaces weighted particles
with $N$ particles with the weight equal to $1/N$, which provides an approximation to the same target pdf.
This proves the theorem.





The fact that we can consider more than one follower of each particle and
reduce the number of followers by resampling is known in the PF literature and
is referred to as prior boosting \cite{gordon93novel,carpenter99building}.
It is used to capture multi-modal likelihood regions.
We stress that the resampling plays in our framework an additional and a very crucial role.
It selects the the most informative random variables (i.e., state space dimensions) as followers of particles.
Since the weight of $x_{<1:t-1>,s}^{(i)}$ is determined by the observations $Z$,
and the resampling uses the weights to selects a follower $x_{<t>} =x_s$
from not yet considered dimensions
$$s \in \{ 1, \dots, m \} \setminus <1:t-1>,$$
the resampling determines the order of RVs, i.e., the bijection $<t>$ for $t=1, \dots m$.
Consequently, the order of RVs is heavily determined by $Z$,
and this order may be different for each particle $(i)$.
This is in strong contrast to the classical PF, where observations $Z$ have no influence on the order of RVs,
which is fixed.


The weight \eqref{eqn:partweight} of each particle is based on the evaluation
how the predicted observations $h(x_{<1:t-1>,s}^{(i)})$ differ from the current observations $z_{<1:t-1>,s}$.
Consequently, in the proposed approach, the weights of different particles
are evaluated with respect to different observations.
Thus, in the proposed approach the value of the weight depends on two factors,
\begin{itemize}
\item how descriptive a given vector of observations $z_{<1:t-1>,s}$ is and
\item how good the prediction $h(x_{<1:t-1>,s}^{(i)})$ of observation $z_{<1:t-1>,s}$ is.
\end{itemize}





%
%In order to be able to execute the proposed algorithm, we need to define the proposal distribution
%$q(x^{(i)}_{<t>}|x_{<1:t-1>}^{(i)},Z)=p(x^{(i)}_{<t>}|x_{<1:t-1>}^{(i)})$ and the posterior distribution
%$p(Z|x_{<1:t>}^{(i)})$. We describe their constructions below.
%Now we show that if we keep assumption \eqref{eq:PFa2} unchanged,
%then $p(Z|x_{<1:t>}^{(i)})$ can be factorized.
%We recall that $Z=\{ z_1, \dots, z_m \}$ is a set of observations, and that
%we know the observation $z_t$ corresponds to state $x_t$ by \eqref{eq:PFa2}.
%By recursive factorization we obtain
%\begin{align}\label{eq:f1}
%&p(Z|x_{<1:t>}^{(i)}) = p( z_{1:m} |x_{<1:t>}^{(i)}) \nonumber \\
%&=p(z_{1:m-<t>} | z_{<t>}, x_{<1:t>}^{(i)}) p(z_{<t>} | x_{<1:t>}^{(i)}) \nonumber \\
%&=p(z_{1:m-<t>} | x_{<1:t>}^{(i)}) p(z_{<t>} | x_{<1:t>}^{(i)}) \nonumber \\
%&=p(z_{1:m-<t,t-1>} | x_{<1:t>}^{(i)}) p(z_{<t-1>} | x_{<1:t>}^{(i)}) p(z_{<t>} | x_{<1:t>}^{(i)}) \nonumber \\
%&= \dots = p(z_{1:m-<1:t>} ) \prod_{k=1}^t p(z_{<k>} | x_{<1:t>}^{(i)}) \nonumber \\
%& \propto \prod_{k=1}^t p(z_{<k>} | x_{<1:t>}^{(i)}),
%\end{align}
%where the subscript $1:m - <t>$ denotes the sequence of indices $(1, \dots, m )$ without the index $<t> = i_t$
%and we assume that $t < m$.
%The last equation follows from the fact that at a given iteration $t$ of the particle filter
%$p(z_{1:m-<1:t>} )$ is a constant that does not depend on the state vector $x_{<1:t>}^{(i)}$.
%Hence we can drop $p(z_{1:m-<1:t>} )$.
%If $t = m$, then $p(z_{1:m-<1:t>} )$ disappears and we have an equality in \eqref{eq:f1}.
%
%Since the proposed PF framework performs sequential filtering, there are two important issues that need to be addressed:
%setting the initial particles $x^{(i)}_{<1>}$ for each $i=1, \dots, N$
%(Section \ref{sec:propdist})
%and the number of particles, which will be determined experimentally.
%We only mention here that the proposed approach is in some sense robust to the initial correspondences,
%since it does not matter with which correspondence we start as long as we start at some element of
%the optimal set of correspondences $\widehat{x}_{1:m} = (\widehat{x}_1, \dots, \widehat{x}_m)$.
%




\section{Rao-Blackwellized Particle Filter}

The following derivation is inspired by Rao-Blackwellized particle filter often used in SLAM \cite{Grisetti07}.
Our goal is to derive the posterior $ p(x_{1:m}, \mu |Z)$ not only over the set of stats $x_{1:m}$ but also over possible shape models $\mu$.
To make this estimation possible, we use the following factorization
\begin{align}\label{eq:RB1}
& p(x_{1:t}, \mu |Z) = p(\mu | x_{1:t}, Z)  \ p(x_{1:t} | Z)
\end{align}
for $t=1, \dots, m$.
This factorization allows us to first estimate only the set of states
and then compute the shape model.
The key assumption here is that the shape model is a deterministic function of the set of states.
Hence, the posterior over models $p(\mu | x_{1:t}, Z)$ can be computed analytically
since $x_{1:t}, Z$ are known.
This technique is often referred to as Rao-Blackwellization.

By substituting Eq.~\eqref{eq:RB1} to Eq.~\eqref{eq:PFn1.1}, we obtain
\begin{align}\label{eq:RB2}
& p(x_{1:m}, \mu |Z) \approx p(\mu | x_{1:m}, Z)  \sum_{i=1}^N w(x_{1:m}^{(i)}) \delta(x_{1:m} - x_{1:m}^{(i)}) \nonumber \\
&=\sum_{i=1}^N p(\mu | x_{1:m}, Z) w(x_{1:m}^{(i)}) \delta(x_{1:m} - x_{1:m}^{(i)})
\end{align}
By denoting $ u(x_{1:m}^{(i)}) = p(\mu | x_{1:m}^{(i)}, Z) w(x_{1:m}^{(i)})$, we obtain
\begin{align}\label{eq:RB2}
& p(x_{1:m}, \mu |Z) \approx \sum_{i=1}^N u(x_{1:m}^{(i)}) \delta(x_{1:m} - x_{1:m}^{(i)}),.
\end{align}
Consequently, if we modify the weight \eqref{eqn:partweight} of each particle to
\begin{align}\label{eqn:RBpartweight}
&u(x_{<1:t-1>}^{(i)}) = p(\mu | x_{<1:t-1>}^{(i)}, Z) p(Z|x_{<1:t-1>}^{(i)}),
\end{align}
we obtain a PF algorithm that estimates the posterior $p(x_{1:t}, \mu |Z)$.








\section{Improved Importance Sampling from MRFs}
In this section we will develop the key motivational story behind our work i.e. sampling from the posterior distribution of a Markov Random Field. As is well known in computer vision community many low level problems have been modeled using MRFs. A limited work is done in modeling mid-level and high-level problems using MRFs not because they cannot be modeled in such a way but the inference process becomes computationally hard. In fact for a general MRF, inference becomes NP-hard. This is because inference in MRFs is closely related to assignment problems. Except for simple costs of assignments the assignment problem quickly turns out to be intractable. Our contribution is to be able to infer to a reasonable degree on the problem instances of contour grouping and object recognition where we can exploit certain results from importance sampling theory to effectively navigate the exponential space of assignments. this has applications in clustering, object recognition, contour grouping and any general applications where inferences can be modeled using Markov Random Fields or Conditional Random Fields.
\section{Markov Random Fields}\label{sec:mrf}
Let $S=\{s_1,s_2,\ldots,s_m\}$ be a family of random variables (RVs), which take the values $e_i\in E=\{e_1,e_2,\ldots,e_n\}$. $E$ is called a set of labels. $s_i=e_i$ denotes the random event that the RV $s_i$ gets $e_i$ assigned. If $f$ is a short hand for $(s_1=e_1,s_2=e_2,\ldots,s_m=e_m)$ and $S$ forms a MRF as per definitions in \cite{Stan08}:
\begin{align}
p(f)=\frac{1}{Z}e^{-\frac{1}{T}U(f)}
\end{align}
where $Z$ is the normalizing constant and $T$ is the temperature parameter determining the sharpness of the distribution: high-temperature makes all configurations equally likely since the effect of $U(f)$ goes down. Simply put $e^{-\frac{1}{T}U(f)}\rightarrow 1$ as $T\rightarrow \infty$.

At given $T$ and particular $U$ we can sample ``patterns'' of assignments by sampling from $p$.

The state space to be explored to understand the patterns of the posterior $p$ is $n^m$. At lower $T$s it finding Maximum A Posteriori (MAP) estimates of $p$ has important applications in high-level vision problems. Thus one of the goals of sampling from MRF is:
\begin{align}
& \widehat{f} = \underset{f \in E^S }{\operatorname{argmax}} \ p(f)
\label{eqG1}
\end{align}
There have been many sampling algorithms like Gibbs sampler, Hot Coupling (\cite{hot_coupling_nips2005}), Tree sampling, Swendsen-Wang sampling etc. But most of them assume restrictive conditional independences. Recently Hamze et. al. proposed a very generic importance sampling method called Large Flip Importance Sampling (LFIS) to sample from the posterior \cite{hamze_freitas_uai2007}. The main motivation for their approach comes from N-Fold Way (NFW, \cite{bortz_etal:jcp1975}) and Tabu search (\cite{glover:jc1989}) where they use heuristics to improve the sampling of the exponential state space using memory and heuristics to design good moves in the state space. Since the moves are no-longer MCMC in the traditional sense they introduce importance weights to the distinct states visited by $N$ copies of the sampler. Independently there has been an application of similar strategy using particle filters with static observations in \cite{luetal_iccv2009}. In this paper we combine the strengths of both the approaches and present an improved sampler that employs better weigthing scheme and navigational strategy to explore state space so as to compute MAP in an efficient way.

We first present a brief overview of both the approaches and then combine both into an integrated approach.
\subsection{Large Flip Importance Sampling}
Sampling algorithms used in computer vision often tend to use different terminology which might makes things a bit difficult to understand convergence complexity issues. Hence, even though the Gibbs sampler is the simplest MCMC sampler for MRF we would like to present the algorithm below so as to make the connections between different approaches explicit.
\begin{algorithm}
\caption{Gibbs sampler}\label{alg:gibbs} %\algsetup{indent=2em}
\begin{algorithmic}[1]
\STATE Initialize a sample: $f^{(0)}=\left(s_1^{(0)}=e_1^{(0)},s_2^{(0)}=e_2^{(0)},\ldots,s_m^{(0)}=e_m^{(0)}\right)$.
\FOR{$i=1$ to $T$}
\STATE Draw $f_1^{(i)}\sim p(f_1|f_2^{(i-1)},f_3^{(i-1)},\ldots,f_m^{(i-1)})$.
\STATE Draw $f_2^{(i)}\sim p(f_2|f_1^{(i)},f_3^{(i-1)},\ldots,f_m^{(i-1)})$.
\vdots
\STATE Draw $f_m^{(i)}\sim p(f_m|f_1^{(i)},f_2^{(i)},\ldots,f_{m-1}^{(i)})$.
\ENDFOR
\STATE Return $\{f^{(1)},f^{(2)},\ldots,f^{(T)}\}$.
\end{algorithmic}
\end{algorithm}
where $f_i$ indicates the random event $s_i=e_i$.
The above algorithm shows that the invariance kernel at a particular iteration above, for MCMC type convergence is defined as:
\begin{equation}
 K(f,f')=\frac{1}{m}\sum_{i=1}^m p(f_i|f_{M\setminus i})
\end{equation}
The convergence condition says
\begin{equation}
 \|K^{(T)}(f^{(0)},f')-p(f')\|\rightarrow 0,~\textrm{as}~T\rightarrow \infty
\end{equation}
As can be seen from above there are two computational complexity issues in a sampling algorithm viz. (1) simulation complexity, (2) optimization complexity. Simulation complexity is because of the computations needed in drawing the samples while the optimization complexity involves $T$. The longer one runs the algorithm the closer one gets to the optimum. Of course in the most general case $T$ can be exponentially large. Most practical algorithm designs involve in reducing $T$ so as to visit non-trivial and important states of $E^S$ as soon as possible without wasting computing resources. Such design algorithms are called \emph{event driven} MCMC approaches. Motivated from NFW (algorithm below) LFIS was developed to address the ``cycling'' problem of NFW.
\begin{algorithm}
\caption{N-Fold Way}\label{alg:nfw} %\algsetup{indent=2em}
\begin{algorithmic}[1]
\STATE Initialize a sample: $\widehat{f}^{(0)}=\left(s_1^{(0)}=e_1^{(0)},s_2^{(0)}=e_2^{(0)},\ldots,s_m^{(0)}=e_m^{(0)}\right)$.
\FOR{$i=1$ to $T$}
\STATE Draw $\tau\sim \textrm{Geometric}(p_{\textrm{flip}}(\widehat{f}_1^{(i-1)},\widehat{f}_2^{(i-1)},\ldots,\widehat{f}_m^{(i-1)}))$.
\STATE Draw $\widehat{f}^{(i)}\sim \nu(j,e,\widehat{f}^{(i-1)})$.
\STATE Set $\Theta_i=\Theta_i+\tau$.
\ENDFOR
\STATE Return $\{\widehat{f}^{(1)},\widehat{f}^{(2)},\ldots,\widehat{f}^{(T)}\}$.
\end{algorithmic}
\end{algorithm}
where $\widehat{f}^{(i)}=f^{(\Theta_i\ldots\Theta_{i+1}-1)}$ i.e. the above algorithm avoids visiting duplicate states and the actual mixing time of the equivalent Gibbs sampler is $\Theta_T$. Hence NFW effectively reduced the mixing time by cleverly covering more states in less time. $\textrm{Geometric}(p_{\textrm{flip}}(\widehat{f}_1^{(i-1)},\widehat{f}_2^{(i-1)},\ldots,\widehat{f}_m^{(i-1)}))$ is the probability that the $\widehat{f}^{(i-1)}$ changes over the next iteration. $\nu(j,e,\widehat{f}^{(i-1)})$ is the discrete posterior probability of the joint event "variable $j$ assumes value $e$" given that a change in state occurred, i.e. $\widehat{f}^{(i)}\neq \widehat{f}^{(i-1)}$. Essentially the NFW simulates $T-1$ \emph{flips}.

One reason for longer mixing times for some posteriors is the conditionals involved in the kernel can make the state not change i.e. $f=f'$ for a long time. NFW exploits the clever perspective of ``simulating'' the long waiting time if the conditionals for change of state space visited can be computed. The main problem involved is in computing this ``change conditionals'' which if not computed using exhaustive flips can result in ``cycling'' of the states i.e. same states get visited over and over again instead of visiting new states. LFIS avoids this by having explicit memory and running $N$ copies of the sampler. After the run it computes importance weights to each of the unique states visited as follows:
\begin{equation}
 w^{(i)}=\frac{p(\widetilde{f}^{(i)})}{\sum_{i=1}^{S_L}p(\widetilde{f}^{(i)})}\label{eq:lfis_wgt}
\end{equation}
where $S_L\subset E^S$. These importance samples represent the posterior $p(f|S_L)$ instead of $p(f)$. Obviously if the algorithm is run long enough $S_L\rightarrow E^S$ and the samples represent the posterior $p(f)$.
\begin{algorithm}
\caption{Large Flip Importance Sampling}\label{alg:lfis} %\algsetup{indent=2em}
\begin{algorithmic}[1]
\STATE Initialize a sample: $\widehat{f}^{(0)}=\left(s_1^{(0)}=e_1^{(0)},s_2^{(0)}=e_2^{(0)},\ldots,s_m^{(0)}=e_m^{(0)}\right)$.
\STATE $k=1$.
\STATE $\Gamma_k\sim\mathcal{U}[\Gamma_{\textrm{min}},\Gamma_{\textrm{max}}]$
\STATE $n=0$.
\STATE $F_{0}^{k}=\phi$.
\FOR{$i=1$ to $T$}
\STATE Draw $\widehat{f}^{(i)}\sim \nu(j,e,\widehat{f}^{(i-1)})$.
\STATE $\mathcal{F}_n^{k}=\mathcal{F}_{n-1}^k\cup \widehat{f}^{(i)}$.
\STATE $n=n+1$.
\IF{$n=\Gamma_k$}
\STATE $k=k+1$.
\STATE $\Gamma_k\sim\mathcal{U}[\Gamma_{\textrm{min}},\Gamma_{\textrm{max}}]$
\STATE $n=0$.
\STATE $F_{0}^{k}=\phi$.
\ENDIF
\ENDFOR
\STATE Pick all the \emph{unique} states, $\{\widetilde{f}^{(i)}\}_{i=1}^T$, from $\{\mathcal{F}_{\Gamma_k}^k\}$ for all available $k$.
\STATE Assign importance weights to those unique states as per Eq. \eqref{eq:lfis_wgt}.
\STATE Return $\{(\widetilde{f}^{(1)},w^{(1)}),(\widetilde{f}^{(2)},w^{(2)}),\ldots,(\widetilde{f}^{(T)},w^{(T)})\}$.
\end{algorithmic}
\end{algorithm}
\subsection{Particle Filters with Static Prior}
Another recent work that uses importance sampling for sampling from MRF is in \cite{luetal_iccv2009}. Owing to the similarity to the notations used in robot mapping, they call it particle filters with static prior. There the idea is to sample states \emph{sequentially} and in such a way that samples build-up the assignments bottom-up using clever conditionals. Below we summarize that approach using the terminology used in this paper. $f_{<1:t>}$ denotes $t$ random variables are instantiated not necessarily $s_1$ through $s_t$. Also $p(f_{<t>}|f_{<1:t-1>})$ denotes the conditional probability distribution of instantiating \emph{a} random variable which has not been instantiated so far.

Assume there are $N$ particles. Let's take a closer look at the journey of one particle.
\begin{algorithm}
\caption{Particle filter with static prior}\label{alg:pf-static} %\algsetup{indent=2em}
\begin{algorithmic}[1]
\FOR{$i=1$ to $N$}
\STATE Draw $f_{<1>}^{(i)}\sim p(f_{<1>})$.
\STATE Set $w_{<1>}^{(i)}=p(f_{<1>}^{(i)})$.
\STATE Draw $f_{<2>}^{(i)}\sim p(f_{<2>}|f_{<1>}{(i)})$.
\STATE Set $f_{<1:2>}^{(i)}=\{f_{<2>}^{(i)},f_{<1>}^{(i)}\}$.
\STATE Set $w_{<1:2>}^{(i)}=w_{<1>}^{(i)}*p(f_{<2>}^{(i)}|f_{<1>}{(i)})$.
\vdots
\STATE Draw $f_{<m>}^{(i)}\sim p(f_{<m>}|f_{<1:m-1>}^{(i)})$.
\STATE Set $f_{<1:m>}^{(i)}=\{f_{<m>}^{(i)},f_{<1:m-1>}^{(i)}\}$.
\STATE Set $w_{<1:m>}^{(i)}=w_{<1:m-1>}^{(i)}*p(f_{<m>}^{(i)}|f_{<1:m-1>}^{(i)})$.
\ENDFOR
\STATE Return $\{(f_{<1:m>}^{(1)},w_{<1:m>}^{(1)}),(f_{<1:m>}^{(2)},w_{<1:m>}^{(2)}),\ldots,(f_{<1:m>}^{(N)},w_{<1:m>}^{(N)})\}$.
\end{algorithmic}
\end{algorithm}
By definition $f_{<1:m>}\equiv f$ since $f$ has $m$ random variables. Since there are only $m$ random variables the equivalent of number of iterations ($T$) in Gibbs sampler needed is only $m$. The importance samples $\{(f_{<1:m>}^{(i)},w_{<1:m>}^{(i)})\}$ represent the posterior $p(f_{<m>}|f_{<1:m-1>}^{(i)})$. But there are three key observations:
\begin{itemize}
 \item If $N$ is large enough and the conditionals are designed properly the samples would represent non-trivial patterns in the posterior and since in \cite{luetal_iccv2009} the goal was to identify good matching between model and image contours the samples were able to serve the purpose.
\item The samples are weighted \emph{incrementally} using recursive importance weighting unlike LFIS where the samples are weighted at the end of the sampling process. For derivations of the recursive importance weighting please refer to \cite{luetal_iccv2009}.
\item Although the samples are not generated in a standard MCMC fashion like Gibbs sampler the importance weighting and the conditional distributions result in useful samples. This is the key basis for LFIS also.
\end{itemize}
%One observation about the algorithms is that algorithms \ref{alg:gibbs}, \ref{alg:nfw} give out \emph{random} samples from the posterior
%(i.e. their weights are all equal) while algorithms \ref{alg:lfis}, \ref{alg:pf-static} give \emph{importance} samples.

\section*{Acknowledgments}
This work was supported in part by NSF IIS-0812118 and by DOE DE-FG52-06NA27508 grants.
N. Adluru is supported by CIBM and MIR at the UW-Madison.
{\small
%\footnotesize
\bibliographystyle{ieee}
%\bibliography{bib/shape,bib/EMrobot,bib/contours,bib/seg,bib/Symmetry,bib/symmskel,bib/shape08,bib/proposal,bib/shape08Jianbo,bib/shapeMRF,bib/ChengEn}
\bibliography{../../Bib/RobotMapping/EMrobot,../../Bib/RobotMapping/FFS,../../Bib/Shape/shape,../../Bib/Shape/contours,../../Bib/BMI/bmi,../../Bib/BMI/mrfs,../../embs2009,../../Bib/Shape/seg}
}

\end{document} 