\documentclass[11pt]{article}

\usepackage{epsfig}
\usepackage{multirow}




\renewcommand{\baselinestretch}{1.1}


\setlength{\topmargin}{-0.5in}
\setlength{\oddsidemargin}{-0.00in}
\setlength{\textwidth}{6.5in}
\setlength{\textheight}{9in}



\graphicspath{{./}{img_eps/}}
\newcommand{\Local}{./img_eps}

\def\I{{\bf I}}
\def\J{{\bf J}}
\def\B{{\bf B}}
\def\S{{\rm SUM}}
\def\M{{\rm MAX}}
\def\E{{\rm E}}


\begin{document}




\begin{center}
{\Huge Learning Active Basis Models \\ by EM-Type Algorithms}
\end{center}

\begin{center}
Zhangzhang Si$^1$, Haifeng Gong$^{1,2}$, Song-Chun Zhu$^{1,2}$, and
Ying Nian Wu$^1$\\
\vspace{3mm} \noindent $^1$ Department of
Statistics, University of
California, Los Angeles\\
$^2$ Lotus Hill Research Institute, Ezhou, China\\
$\{$zzsi, hfgong, sczhu, ywu$\}$@stat.ucla.edu\\
\vspace{3mm} September 3, 2009
\end{center}




\begin{center}
{\bf Abstract}
\end{center}
EM algorithm is a convenient tool for maximum likelihood model
fitting when the data are incomplete or when there are latent
variables or hidden states. In this review article, we explain
that EM algorithm is a natural computational scheme for learning
image templates of object categories where the learning is not
fully supervised. We represent an image template by an active
basis model, which is a linear composition of a selected set of
localized, elongated and oriented wavelet elements that are
allowed to slightly perturb their locations and orientations to
account for the deformations of object shapes. The model can be
easily learned when the objects in the training images are of the
same pose, and appear at the same location and scale. This is
often called supervised learning. In the situation where the
objects may appear at different unknown locations, orientations,
and scales in the training images, we have to incorporate the
unknown locations, orientations, and scales as latent variables
into the image generation process, and learn the template by
EM-type algorithms. The E-step imputes the unknown locations,
orientations, and scales based on the currently learned template.
This step can be considered self-supervision, which involves using
the current template to recognize the objects in the training
images. The M-step then relearns the template based on the imputed
locations, orientations, and scales, and this is essentially the
same as supervised learning. So the EM learning process iterates
between recognition and supervised learning. We illustrate this
scheme by several experiments.

\vspace{3mm} \noindent{\bf Keywords}: Generative models, object
recognition, wavelet sparse coding.


\newpage

%%%%%%%%% BODY TEXT
\section{Introduction: EM learning scheme}

The EM algorithm \cite{DLR} and its variations \cite{Meng} have been
widely used for maximum likelihood estimation of statistical models
when the data are incompletely observed or when there are latent
variables or hidden states. This algorithm is an iterative
computational scheme, where the E-step imputes the missing data or
the latent variables given the currently estimated model, and the
M-step re-estimates the model given the imputed missing data or
latent variables. Besides its simplicity and stability, a key
feature that distinguishes the EM algorithm from other numerical
methods is its interpretability: both the E-step and the M-step
readily admit natural interpretations in a variety of contexts. This
makes the EM algorithm rich, meaningful and inspiring.

In this review article, we shall focus on one important context
where the EM algorithm is useful and meaningful, that is, learning
patterns from signals in the settings that are not fully
supervised. In this context, the E-step can be interpreted as
carrying out the recognition task using the currently learned
model of the pattern. The M-step can be interpreted as relearning
the pattern in the supervised setting, which can often be easily
accomplished.

This EM learning scheme has been used in both speech and vision.
In speech recognition, the training of the hidden Markov model
\cite{Rabiner} involves the imputation of the hidden states in the
E-step by the forward and backward algorithms. The M-step computes
the transition and emission frequencies. In computer vision, we
want to learn models for different categories of objects, such as
horses, birds, bikes, etc. The learning is often easy when the
objects in the training images are aligned, in the sense that the
objects appear at the same pose, same location and same scale in
the training images, which are defined on a common image lattice
that is the bounding box of the objects. This is often called
supervised learning. However, it is often the case that the
objects appear at different unknown locations, orientations and
scales in the training images. In such a situation, we have to
incorporate the unknown locations, orientations and scales as
latent variables in the image generation process, and use the EM
algorithm to learn the model for the objects. In the EM learning
process, the E-step imputes the unknown location, orientation, and
scale of the object in each training image, based on the currently
learned model. This step uses the current model to recognize the
object in each training image, that is, where it is, at what
orientation and scale. The imputation of the latent variables
enables us to align the training images, so that the objects
appear at the same location, orientation and scale. The M-step
then relearns the model from the aligned images by carrying out
supervised learning. So the EM learning process iterates between
recognition and supervised learning. Recognition is the goal of
learning the model, and it serves as the self-supervision step of
the learning process. The EM algorithm has been used by Fergus,
Perona and Zisserman \cite{Perona} in training the constellation
model for objects.

In this article we shall illustrate EM learning or EM-like
learning by training an active basis model \cite{Wu, Wu1} that we
have recently developed for deformable templates \cite{Yuille,
Amit} of object shapes. In this model, a template is represented
by a linear composition of a set of localized, elongated and
oriented wavelet elements at selected locations, scales and
orientations, and these wavelet elements are allowed to slightly
perturb their locations and orientations to account for the shape
deformations of the objects. In the supervised setting, the active
basis model can be learned by a shared sketch algorithm, which
selects the wavelet elements sequentially. When a wavelet element
is selected, it is shared by all the training images, in the sense
that a perturbed version of this element seeks to sketch a local
edge segment in each training image. In the situations where
learning is not fully supervised, the learning of the active basis
model can be accomplished by the EM-type algorithms. The E-step
recognizes the object in each training image by matching the image
with the currently learned active basis template. This enables us
to align the images. The M-step then relearns the template from
the aligned images by the shared sketch algorithm.

We would like to point out that the EM algorithm for learning the
active basis model is different than the traditional EM
algorithm, where the model structure is fixed and only the
parameters need to be estimated. In our implementation of the EM
algorithm, the M-step needs to select the wavelet elements in
addition to estimating the parameters associated with the selected
elements. Both the selection of the elements and the estimation of
the parameters are accomplished by maximizing or increasing the
complete-data log-likelihood. So the EM algorithm is used for
estimating both the model structure and the associated parameters.

Readers who wish to learn more about the active basis model are
referred to our recent paper \cite{Wu1}, which is written for the
computer vision community. Compared to that paper, this review
paper is written for the statistical community. In this paper we
introduce the active basis model from an algorithmic perspective,
starting from the familiar problem of variable selection in linear
regression. This paper also provides more details about the
EM-type algorithms than \cite{Wu1}. We wish to convey to the
statistical audience that the problem of vision in general and
object recognition in particular is essentially a statistical
problem. We even hope that this article may attract some
statisticians to work on this interesting but challenging problem.

Section 2 introduces the active basis model for representing
deformable templates, and describes the shared sketch algorithm
for supervised learning. Section 3 presents the EM algorithm for
learning the active basis model in the settings that are not fully
supervised. Section 4 concludes with a discussion.



%------------------------------------------------------------------------
\section{Active basis model: an algorithmic tour}

The active basis model is a natural generalization of the wavelet
regression model. In this section we first explain the background
and motivation for the active basis model. Then we work through a
series of variable selection algorithms for wavelet regression,
where the active basis model emerges naturally.

\subsection{From wavelet regression to active basis}

\subsubsection{$p>n$ regression and variable section}

Wavelets have proven to be immensely useful for signal analysis and
representation \cite{Donoho_compression}. Various dictionaries of
wavelets have been designed for different types of signals or
function spaces \cite{Candes_Donoho, Simoncelli_shiftable}. Two key
factors underlying the successes of wavelets are the sparsity of the
representation and the efficiency of the analysis. Specifically, a
signal can typically be represented by a linear superposition of a
small number of wavelet elements selected from an appropriate
dictionary. The selection can be accomplished by efficient
algorithms such as matching pursuit \cite{Mallat_Zhang} and basis
pursuit \cite{Chen_basis_pursuit}.

From a linear regression perspective, a signal can be considered a
response variable, and the wavelet elements in the dictionary can
be considered the predictor variables or regressors. The number of
elements in a dictionary can often be much greater than the
dimensionality of the signal, so this is the so-called ``$p > n$''
problem. The selection of the wavelet elements is the variable
selection problem in linear regression. The matching pursuit
algorithm \cite{Mallat_Zhang} is the forward selection method, and
the basis pursuit \cite{Chen_basis_pursuit} is the lasso method
\cite{lasso}.



\subsubsection{Gabor wavelets and simple V1 cells}

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=gf.eps,height=2in}
\end{tabular}
\end{center}
\caption{A collection of Gabor wavelets at different locations,
orientations and scales. Each Gabor wavelet element is a sine or
cosine wave multiplied by an elongated and oriented Gaussian
function. The wave propagates along the shorter axis of the
Gaussian function.}
 \label{fig:gf}
\end{figure}

Interestingly, wavelet sparse coding also appears to be employed
by the biological visual system for representing natural images.
By assuming the sparsity of the linear representation, Olshausen
and Field \cite{Olshausen_Field} were able to learn from natural
images a dictionary of localized, elongated, and oriented basis
functions that resemble the Gabor wavelets. Similar wavelets were
also obtained by independent component analysis of natural images
\cite{Bell_Sejnowski}. From a linear regression perspective,
Olshausen and Field essentially asked the following question:
Given a sample of response vectors (i.e., natural images), can we
find a dictionary of predictor vectors or regressors (i.e., basis
functions or basis elements), so that each response vector can be
represented as a linear combination of a small number of
regressors selected from the dictionary? Of course, for different
response vectors, different sets of regressors may be selected
from the dictionary.

Figure (\ref{fig:gf}) displays a collection of Gabor wavelet
elements at different locations, orientations and scales. These
are sine and cosine waves multiplied by elongated and oriented
Gaussian functions, where the waves propagate along the shorter
axes of the Gaussian functions. Such Gabor wavelets have been
proposed as mathematical models for the receptive fields of the
simple cells of the primary visual cortex or V1 \cite{Daugman}.

The dictionary of all the Gabor wavelet elements can be very
large, because at each pixel of the image domain, there can be
many Gabor wavelet elements tuned to different scales and
orientations. According to Olshausen and Field
\cite{Olshausen_Field}, the biological visual system represents a
natural image by a linear superposition of a small number of Gabor
wavelet elements selected from such a dictionary.

\subsubsection{From generic classes to specific categories}

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|}
\hline
& \psfig{figure=deer1.eps,height= 0.5in}&
\psfig{figure=deer2.eps,height=0.5in}&
\psfig{figure=deer3.eps,height=0.5in}&
\psfig{figure=deer4.eps,height=0.5in}&
\psfig{figure=deer5.eps,height=0.5in}&
\psfig{figure=deer6.eps,height=0.5in}&
\psfig{figure=deer7.eps,height=0.5in}&
\psfig{figure=deer8.eps,height=0.5in}\\

 \psfig{figure=ex1_5_tsk50.eps,height=0.5in}&
\psfig{figure=deer_deformed_sk1.eps,height=0.5in}&
\psfig{figure=deer_deformed_sk2.eps,height=0.5in}&
\psfig{figure=deer_deformed_sk3.eps,height=0.5in}&
\psfig{figure=deer_deformed_sk4.eps,height=0.5in}&
\psfig{figure=deer_deformed_sk5.eps,height=0.5in}&
\psfig{figure=deer_deformed_sk6.eps,height=0.5in}&
\psfig{figure=deer_deformed_sk7.eps,height=0.5in}&
\psfig{figure=deer_deformed_sk8.eps,height=0.5in}\\

\psfig{figure=ex1_5_tsk_large15.eps,height=0.5in}&
\psfig{figure=deer_deformed_sk_large1.eps,height= 0.5in}&
\psfig{figure=deer_deformed_sk_large2.eps,height=0.5in}&
\psfig{figure=deer_deformed_sk_large3.eps,height=0.5in}&
\psfig{figure=deer_deformed_sk_large4.eps,height=0.5in}&
\psfig{figure=deer_deformed_sk_large5.eps,height=0.5in}&
\psfig{figure=deer_deformed_sk_large6.eps,height=0.5in}&
\psfig{figure=deer_deformed_sk_large7.eps,height=0.5in}&
\psfig{figure=deer_deformed_sk_large8.eps,height=0.5in}\\
\psfig{figure=multiscale_template.eps,height=0.5in}&
\psfig{figure=synthesized_deer1.eps,height= 0.5in}&
\psfig{figure=synthesized_deer2.eps,height=0.5in}&
\psfig{figure=synthesized_deer3.eps,height=0.5in}&
\psfig{figure=synthesized_deer4.eps,height=0.5in}&
\psfig{figure=synthesized_deer5.eps,height=0.5in}&
\psfig{figure=synthesized_deer6.eps,height=0.5in}&
\psfig{figure=synthesized_deer7.eps,height=0.5in}&
\psfig{figure=synthesized_deer8.eps,height=0.5in}\\
\hline
\end{tabular}
\end{center}
\caption{Active basis templates. Each Gabor wavelet element is
illustrated by a bar of the same length and at the same location
and orientation as the corresponding element. The first row
displays the training images. The second row displays the
templates composed of 50 Gabor wavelet elements at a fixed scale,
where the first template is the common deformable template, and
the other templates are deformed templates for coding the
corresponding images. The third row displays the templates
composed of 15 Gabor wavelet elements at a scale about twice as
large as those in the second row. In the last row, the template is
composed of 50 wavelet elements at multiple scales, where larger
Gabor elements are illustrated by bars of lighter shades. The rest
of the images are reconstructed by linear superpositions of the
wavelet elements of the deformed templates. }
 \label{fig:ab}
\end{figure}

Wavelets are designed for generic function classes or learned from
generic ensembles such as natural images, under the generic
principle of sparsity. While such generality offers enormous scope
for the applicability of wavelets, sparsity alone is clearly
inadequate for modeling specific patterns. Recently, we have
developed an active basis model for images of various object
classes \cite{Wu, Wu1}. The model is a natural consequence of
seeking a common wavelet representation simultaneously for
multiple training images from the same object category.

Figure (\ref{fig:ab}) illustrates the basic idea. In the first row
there are 8 images of deer. The images are of the same size of 122
$\times$ 120 pixels. The deer appear at the same location, scale
and pose in these images. For these very similar images, we want
to seek a common wavelet representation, instead of coding each
image individually. Specifically, we want these images to be
represented by similar sets of wavelet elements, with similar
coefficients. We can achieve this by selecting a common set of
wavelet elements, while allowing these wavelet elements to locally
perturb their locations and orientations before they are linearly
combined to code each individual image. The perturbations are
introduced to account for shape deformations in the deer. The
linear basis formed by such perturbable wavelet elements is called
an active basis.

This is illustrated by the second and third rows of Figure
(\ref{fig:ab}). In each row the first plot displays the common set
of Gabor wavelet elements selected from a dictionary. The
dictionary  consists of Gabor wavelets at all the locations and
orientations, but at a fixed scale. Each Gabor wavelet element is
symbolically illustrated by a bar at the same location and
orientation and with the same length as the corresponding Gabor
wavelet. So the active basis formed by the selected Gabor wavelet
elements can be interpreted as a template, as if each element is a
stroke for sketching the template. The templates in the second and
third rows are learned using dictionaries of Gabor wavelets at two
different scales, with the scale of the third row about twice as
large as the scale of the second row. The number of Gabor wavelet
elements of the template in the second row is 50, while the number
of elements of the template in the third row is 15. Currently, we
treat this number as a tuning parameter, although they can be
determined in a more principled way.

Within each of the second and third rows, and for each training
image, we plot the Gabor wavelet elements that are actually used
to represent the corresponding image. These elements are perturbed
versions of the corresponding elements in the first column. So the
templates in the first column are deformable templates, and the
templates in the remaining columns are deformed templates. Thus,
the goal of seeking a common wavelet representation for images
from the same object category leads us to formulate the active
basis, which is a deformable template for the images from the
object category.


In the last row of Figure (\ref{fig:ab}), the common template is
learned by selecting from a dictionary that consists of Gabor
wavelet elements at multiple scales instead of a fixed scale. The
number of wavelet elements is 50. In addition to Gabor wavelet
elements, we also include the center-surround difference of
Gaussian wavelet elements in the dictionary. Such isotropic
wavelet elements are of large scales, and they mainly capture the
regional contrasts in the images. In the template in the last row,
larger Gabor wavelet elements are illustrated by bars of lighter
shades. The difference of Gaussian elements are illustrated by
circles. The remaining images are reconstructed by such
multi-scale wavelet representations, where each image is a linear
superposition of the Gabor and difference of Gaussian wavelet
elements of the corresponding deformed templates.


\begin{figure}[!htb]
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|}
\hline \psfig{figure=ex1_5_tsk5.eps,height= 0.5in} &
\psfig{figure=ex1_5_tsk10.eps,height=0.5in} &
\psfig{figure=ex1_5_tsk20.eps,height=0.5in} &
\psfig{figure=ex1_5_tsk30.eps,height=0.5in} &
\psfig{figure=ex1_5_tsk40.eps,height=0.5in} &
\psfig{figure=ex1_5_tsk50.eps,height=0.5in}\\
\hline $n=5$ & $n=10$ & $n=20$ & $n=30$ & $n=40$ &
$n=50$ \\
\hline
\psfig{figure=ex1_5_tsk_large1.eps,height=0.5in} &
\psfig{figure=ex1_5_tsk_large2.eps,height=0.5in} &
\psfig{figure=ex1_5_tsk_large3.eps,height=0.5in} &
\psfig{figure=ex1_5_tsk_large5.eps,height=0.5in} &
\psfig{figure=ex1_5_tsk_large10.eps,height=0.5in} &
\psfig{figure=ex1_5_tsk_large15.eps,height=0.5in}
\\
\hline $n=1$ & $n=2$ & $n=3$ & $n=5$ & $n=10$ &
$n=15$ \\
\hline
\end{tabular}
\end{center}
\caption{Shared sketch process for learning the active basis
templates at two different scales.} \label{fig:ss}
\end{figure}


The active basis can be learned by the shared sketch algorithm
that we have recently developed \cite{Wu, Wu1}. This algorithm can
be considered a paralleled version of the matching pursuit
algorithm \cite{Mallat_Zhang}. It can also be considered a
modification of the projection pursuit algorithm \cite{Friedman}.
The algorithm selects the wavelet elements sequentially from the
dictionary. Each time when an element is selected, it is shared by
all the training images in the sense that a perturbed version of
this element is included in the linear representation of each
image. Figure (\ref{fig:ss}) illustrates the shared sketch process
for obtaining the templates displayed in the second and third rows
of Figure (\ref{fig:ab}).

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{cc}
\psfig{figure=deerfind_original.eps,height=1in}&
\psfig{figure=deerfind_sketch.eps,height=1in}
\end{tabular}
\end{center}
\caption{Left: Testing image. Right: Object is detected and
sketched by the deformed template.}
 \label{fig:find}
\end{figure}
While selecting the wavelet elements of the active basis, we also
estimate the distributions of their coefficients from the training
images. This gives us a statistical model for the images. After
learning this model, we can then use it to recognize the same type
of objects in testing images. See Figure (\ref{fig:find}) for an
example. The image on the left is the observed testing image. We
scan the learned template of deer over this image, and at each
location, we match the template to the image by deforming the
learned template.  The template matching is scored by the
log-likelihood of the statistical model. We also scan the template
over multiple resolutions of the image to account for the unknown
scale of the object in the image. Then we choose the resolution
and location of the image with the maximum likelihood score, and
superpose on the image the deformed template matched to the image,
as shown by the image on the right in Figure (\ref{fig:find}).
This process can be accomplished by a cortex-like architecture of
sum maps and max maps, to be described in Subsection
(\ref{sect:sum-max}). In machine learning and computer vision
literature, detecting or classifying objects using the learned
model is often called inference. The inference algorithm is often
a part of the learning algorithm. For the active basis model, both
learning and inference can be formulated as maximum likelihood
estimation problems.

\subsubsection{Local maximum pooling and complex V1 cells}

Besides wavelet sparse coding for V1 simple cells, another
inspiration to the active basis model also comes from
neuroscience. Riesenhuber and Poggio \cite{Riesenhuber_Poggio}
observed that the complex cells of the primary visual cortex or V1
appear to perform local maximum pooling of the responses from
simple cells. From the perspective of the active basis model, this
corresponds to estimating the perturbations of the wavelet
elements of the active basis template, so that the template is
deformed to match the observed image. Therefore, if we are to
believe Olshausen and Field's theory on wavelet sparse coding
\cite{Olshausen_Field} and Riesenhuber and Poggio's theory on
local maximum pooling, then the active basis model seems to be a
very natural logical consequence.

In the following subsections we shall describe in detail wavelet
sparse coding, the active basis model, and the learning and
inference algorithms.

\subsection{An overcomplete dictionary of Gabor wavelets}

The Gabor wavelets are translated, rotated, and dilated versions
of the following function:
\begin{eqnarray*}
G(x_1, x_2) \propto
\exp\{ -[(x_1/\sigma_1)^2 + (x_2/\sigma_2)^2]/2\} e^{i {x_1}},
\end{eqnarray*}
which is sine-cosine wave multiplied by a Gaussian function. The
Gaussian function is elongated along the $x_2$-axis, with
$\sigma_2 > \sigma_1$, and the sine-cosine wave propagates along
the shorter $x_1$-axis. We truncate the function to make it
locally supported on a finite rectangular domain, so that it has a
well defined length and width.

We then translate, rotate and dilate $G(x_1, x_2)$ to obtain a
general form of the Gabor wavelets:
\begin{eqnarray*}
B_{x_1, x_2, s, \alpha}(x_1', x_2') = G({\tilde{x}_1}/{s},
{\tilde{x}_2}/{s})/s^2,
\end{eqnarray*}
where
\begin{eqnarray*}
\tilde{x_1} =(x_1'-x_1)\cos\alpha +
(x_2'-x_2)\sin\alpha, \\
\tilde{x_2} = -(x_1'-x_1)\sin\alpha + (x_2'-x_2)\cos\alpha.
\end{eqnarray*}
Writing $x = (x_1, x_2)$, each $B_{x, s, \alpha}$ is a localized
function, where $x = (x_1, x_2)$ is the central location, $s$ is
the scale parameter, and $\alpha$ is the orientation. The
frequency of the wave propagation in $B_{x, s, \alpha}$ is $\omega
= 1/s$. $B_{x, s, \alpha}$  = ($B_{x, s, \alpha, 0}$, $B_{x, s,
\alpha, 1}$), where $B_{x, s, \alpha, 0}$ is the even-symmetric
Gabor cosine component, and $B_{x, s, \alpha, 1}$ is the
odd-symmetric Gabor sine component. We always use Gabor wavelets
as pairs of cosine and sine components. We normalize both the
Gabor sine and cosine components to have zero mean and unit
$\ell_2$ norm. For each $B_{x, s, \alpha}$, the pair $B_{x, s,
\alpha, 0}$ and $B_{x, s, \alpha, 1}$ are orthogonal to each
other.

The dictionary of Gabor wavelets is
\[
\Omega = \{B_{x, s, \alpha}, \forall (x, s, \alpha)\}.
\]
We can discretize the orientation so that $\alpha \in \{o \pi/O, o
= 0, ..., O-1\}$, that is, $O$ equally spaced orientations (the
default value of $O$ is 15 in our experiments). In this article we
mostly learn the active basis template at a fixed scale $s$. The
dictionary $\Omega$ is called ``overcomplete'' because the number
of wavelet elements in $\Omega$ is larger than the number of
pixels in the image domain, since at each pixel, there can be many
wavelet elements tuned to different orientations and scales.

For an image $\I(x)$, with $x \in D$, where $D$ is a set of
pixels, such as a rectangular grid, we can project it onto a Gabor
wavelet $B_{x, s, \alpha, \eta}$, $\eta = 0, 1$. The projection of
$\I$ onto $B_{x, s, \alpha, \eta}$, or the Gabor filter response
at $(x, s, \alpha)$,  is
\begin{eqnarray*}
\langle \I, B_{x, s, \alpha, \eta} \rangle = \sum_{x'} \I(x') B_{x, s, \alpha, \eta}(x').
\end{eqnarray*}
The summation is over the finite support of $B_{x, s, \alpha,
\eta}$. We write $\langle \I, B_{x, s, \alpha} \rangle$ =
($\langle \I, B_{x, s, \alpha, 0} \rangle$, $\langle \I, B_{x, s,
\alpha, 1} \rangle$). The local energy is
\begin{eqnarray*}
|\langle
\I, B_{x, s, \alpha} \rangle|^2 = \langle \I, B_{x,  s,
\alpha, 0} \rangle^2+\langle \I, B_{x, s, \alpha, 1}
\rangle^2.
\end{eqnarray*}
$|\langle \I, B_{x, s, \alpha} \rangle|^2$ is the local spectrum
or the magnitude of the local wave in image $\I$ at $(x, s,
\alpha)$.

Let
\begin{eqnarray*}
   \sigma^2_s = \frac{1}{|D|O} \sum_{\alpha} \sum_{x \in D}
   |\langle \I, B_{x, s, \alpha} \rangle|^2, \label{eq:normalize}
\end{eqnarray*}
where $|D|$ is the number of pixels in $\I$, and $O$ is the total
number of orientations. For each image $\I$, we normalize it to
$\I \leftarrow \I/\sigma_s$, so that different images are
comparable.

\subsection{Matching pursuit algorithm}

For an image $\I(x)$ where $x \in D$, we seek to represent it by
\begin{eqnarray}
  \I = \sum_{i=1}^{n} c_i B_{x_i, s, \alpha_i} + U, \label{eq:sc}
\end{eqnarray}
where $(B_{x_i, s, \alpha_i}, i = 1, ..., n) \subset \Omega$ is a
set of Gabor wavelet elements selected from the dictionary
$\Omega$, $c_i$ is the coefficient, and $U$ is the unexplained
residual image. Recall that each $B_{x_i, s, \alpha_i}$ is a pair
of Gabor cosine and sine components. So $B_{x_i, s, \alpha_i} =
(B_{x_i, s, \alpha_i, 0}, B_{x_i, s, \alpha_i, 1})$, $c_{i} =
(c_{i, 0}, c_{i, 1})$, and $c_{i} B_{x_i, s, \alpha_i} = c_{i, 0}
B_{x_i, s, \alpha_i, 0} + c_{i, 1} B_{x_i, s, \alpha_i, 1}$. We
fix the scale parameter $s$.



In the representation (\ref{eq:sc}), $n$ is often assumed to be
small, for example, $n = 50$. So the representation (\ref{eq:sc})
is called sparse representation or sparse coding. This
representation translates a raw intensity image with a huge number
of pixels into a sketch with only a small number of strokes
represented by $\B = (B_{x_i, s, \alpha_i}, i = 1, ..., n)$. Because of the sparsity, $\B$ captures the most visually meaningful elements in the image. The set of wavelet elements $\B = (B_{x_i, s, \alpha_i}, i = 1, ..., n)$ can be selected from $\Omega$ by the matching pursuit algorithm \cite{Mallat_Zhang}, which seeks to minimize $\|\I - \sum_{i=1}^{n}  c_i B_{x_i, s, \alpha_i}\|^2$ by a greedy scheme.\\

{\bf Algorithm 0: Matching pursuit algorithm}
\begin{itemize}
  \item[0] Initialize $i \leftarrow 0$, $U \leftarrow \I$.
  \item[1] Let $i \leftarrow i + 1$. Let $(x_i, \alpha_i) = \arg\max_{x, \alpha}|\langle U, B_{x, s, \alpha}\rangle|^2$.
  \item[2] Let $c_i = \langle U, B_{x_i, s, \alpha_i}\rangle$. Update $U \leftarrow U - c_i B_{x_i, s, \alpha_i}$.
  \item[3] Stop if $i = n$, else go back to 1.
\end{itemize}

In the above algorithm, it is possible that a wavelet element is selected more than once, but this is extremely rare for real images. As to the choice of $n$ or the stopping criterion, we can stop the algorithm if $|c_i|$ is below a threshold.

Readers who are familiar with the so-called ``large $p$ and small $n$'' problem in linear regression may have recognized that wavelet sparse coding is a special case of this problem, where $\I$ is the response vector, and each $B_{x, s, \alpha} \in \Omega$ is a predictor vector. The matching pursuit algorithm is actually the forward selection procedure for variable selection.

The forward selection algorithm in general can be too greedy. But
for image representation, each Gabor wavelet element only explains
away a small part of the image data, and we usually pursue the
elements at a fixed scale, so such a forward selection procedure
is not very greedy in this context.

\subsection{Matching pursuit for multiple images}

Let $\{\I_m, m = 1, ..., M\}$ be a set of training images defined
on a common rectangle lattice $D$, and let us suppose that these
images come from the same object category, where the objects
appear at the same pose, location and scale in these images. We
can model these images by a common set of Gabor wavelet elements,
\begin{eqnarray}
  \I_m = \sum_{i=1}^{n} c_{m, i} B_{x_i, s, \alpha_i} + U_m, \;\; m = 1, ..., M.  \label{eq:sc1}
\end{eqnarray}
$\B = (B_{x_i, s, \alpha_i}, i = 1, ..., n)$ can be considered a common template for these training images. Model (\ref{eq:sc1}) is an extension of model (\ref{eq:sc}).

We can select these elements by applying the matching pursuit algorithm on these multiple images simultaneously. The algorithm  seeks to minimize $\sum_{m=1}^{M}\|\I_m - \sum_{i=1}^{n}  c_{m, i} B_{x_i, s, \alpha_i}\|^2$ by a greedy scheme.\\

{\bf Algorithm 1: Matching pursuit on multiple images}
\begin{itemize}
  \item[0] Initialize $i \leftarrow 0$. For $m = 1, ..., M$, initialize $U_m \leftarrow \I_m$.
  \item[1] $i \leftarrow i + 1$. Select
  \[
   (x_i, \alpha_i) = \arg\max_{x, \alpha} \sum_{m=1}^{M} |\langle U_m, B_{x, s, \alpha}\rangle|^2.
   \]
  \item[2] For $m = 1, ..., M$, let $c_{m, i} = \langle U_m, B_{x_i, s, \alpha_i}\rangle$, and update $U_m \leftarrow U_{m} - c_{m, i} B_{x_i, s, \alpha_i}$.
  \item[3] Stop if $i = n$, else go back to 1.
\end{itemize}

Algorithm 1 is similar to algorithm 0. The difference is that, in Step 1, $(x_i, \alpha_i)$ is selected by maximizing the sum of the squared responses.

\subsection{Active basis and local maximum pooling}

The objects in the training images share similar shapes, but there can still be considerable variations in their shapes. In order to account for the shape deformations, we introduce the perturbations to the common template, and the model becomes
\begin{eqnarray}
  \I_m = \sum_{i=1}^{n} c_{m, i} B_{x_i+\Delta x_{m, i}, s, \alpha_i + \Delta \alpha_{m, i}} + U_m, \;\; m = 1, ..., M.  \label{eq:sc2}
\end{eqnarray}
Again, $\B= (B_{x_i, s, \alpha_i}, i = 1, ..., n)$ can be considered a common template for the training images, but this time, this template is deformable. Specifically, for each image $\I_m$, the wavelet element $B_{x_i, s, \alpha_i}$ is perturbed to $B_{x_i+\Delta x_{m, i}, s, \alpha_i + \Delta \alpha_{m, i}}$, where $\Delta x_{m, i}$ is the perturbation in location, and $\Delta \alpha_{m, i}$ is the perturbation in orientation. $\B_{m} = (B_{x_i+\Delta x_{m, i}, s, \alpha_i + \Delta \alpha_{m, i}}, i = 1, ..., n)$ can be considered the deformed template for coding image $\I_m$. We call the basis formed by $\B = (B_{x_i, s, \alpha_i}, i = 1, ..., n)$ the active basis, and we call $(\Delta x_{m, i}, \Delta \alpha_{m, i}, i = 1, ..., n)$ the activities or perturbations of the basis elements for image $m$. Model (\ref{eq:sc2}) is an extension of model (\ref{eq:sc1}).

Figure (\ref{fig:ab}) illustrates three examples of active basis
templates. In the second and third rows the templates in the first
column are $\B = (B_{x_i, s, \alpha_i}, i = 1, ..., n)$. The scale
parameter $s$ in the second row is smaller than the $s$ in the
third row. For each row, the templates in the remain columns are
the deformed templates $\B_{m} = (B_{x_i+\Delta x_{m, i}, s,
\alpha_i + \Delta \alpha_{m, i}}, i = 1, ..., n)$, for $m = 1,
..., 8$. The template in the last row should be more precisely
represented by $\B = (B_{x_i, s_i, \alpha_i}, i = 1, ..., n)$,
where each element has its own $s_i$ automatically selected
together with $(x_i, \alpha_i)$. In this article we focus on the
situation where we fix $s$ (default length of the wavelet element
is 17 pixels).

For the activity or perturbation of a wavelet element $B_{x, s,
\alpha}$, we assume that $\Delta x = (d \cos\alpha, d\sin\alpha)$,
with $d \in [-b_1, b_1]$. That is, we allow $B_{x, s, \alpha}$ to
shift its location along its normal direction. We also assume
$\Delta \alpha \in [-b_2, b_2]$. $b_1$ and $b_2$ are the bounds
for the allowed displacements in location and orientation (default
values: $b_1 = 6$ pixels, and $b_2 = \pi/15$). We define
\[
A(\alpha) = \left\{ (\Delta x = (d\cos\alpha, d\sin\alpha), \Delta
\alpha): d \in [-b_1, b_1], \Delta \alpha \in [-b_2, b_2]\right\}
\]
the set of all possible activities for a basis element tuned to
orientation $\alpha$.

We can continue to apply the matching pursuit algorithm to the multiple training images, the only difference is that we add a local maximum pooling operation in Steps 1 and 2. The following algorithm is a greedy procedure to minimize the least squares criterion
\begin{eqnarray}
    \sum_{m=1}^{M} \|\I_m -  \sum_{i=1}^{n} c_{m, i} B_{x_i+\Delta x_{m, i}, s, \alpha_i + \Delta \alpha_{m, i}}\|^2. \label{eq:LS}
\end{eqnarray}


{\bf Algorithm 2: Matching pursuit with local maximum pooling}

\begin{itemize}
  \item[0] Initialize $i \leftarrow 0$. For $m = 1, ..., M$, initialize $U_m \leftarrow \I_m$.
  \item[1] $i \leftarrow i + 1$. Select
  \[
   (x_i, \alpha_i) = \arg\max_{x, \alpha} \sum_{m=1}^{M} \max_{(\Delta x, \Delta \alpha)\in A(\alpha)} |\langle U_m, B_{x+\Delta x, s, \alpha+\Delta \alpha}\rangle|^2.
   \]
  \item[2] For $m = 1, ..., M$, retrieve
  \[
     (\Delta x_{m, i}, \Delta \alpha_{m, i}) = \arg\max_{(\Delta x, \Delta \alpha)\in A(\alpha_i)} |\langle U_m, B_{x_i+\Delta x, s, \alpha_i+\Delta \alpha}\rangle|^2.
  \]
  Let $c_{m, i} \leftarrow \langle U_m, B_{x_i+\Delta x_{m, i}, s, \alpha_i+\Delta \alpha_{m, i}}\rangle$, and update $U_m \leftarrow$ $U_m - c_{m, i} B_{x_i+\Delta x_{m, i}, s, \alpha_i+\Delta \alpha_{m, i}}$.
  \item[3] Stop if $i = n$, else go back to 1.
\end{itemize}

Algorithm 2 is similar to Algorithm 1. The difference is that we add an extra local maximization operation in Step 1: $\max_{(\Delta x, \Delta \alpha)\in A(\alpha)} |\langle U_m, B_{x+\Delta x, s, \alpha+\Delta \alpha}\rangle|^2$. With $(x_i, \alpha_i)$ selected in Step 1, Step 2 retrieves the corresponding maximal $(\Delta x, \Delta \alpha)$ for each image.

We can rewrite Algorithm 2 by defining $R_{m}(x, \alpha) = \langle U_m, B_{x, s, \alpha}\rangle$. Then instead of updating the residual image $U_m$ in Step 2, we can update the responses $R_{m}(x, \alpha)$. \\

{\bf Algorithm 2.1: Matching pursuit with local maximum pooling}
\begin{itemize}
  \item[0] Initialize $i \leftarrow 0$. For $m = 1, ..., M$, initalize $R_m(x, \alpha) \leftarrow \langle \I_m, B_{x, s, \alpha}\rangle$ for all $(x, \alpha)$.
  \item[1] $i \leftarrow i + 1$.
  Select
  \[
   (x_i, \alpha_i) = \arg\max_{x, \alpha} \sum_{m=1}^{M}
    \max_{(\Delta x, \Delta \alpha)\in A(\alpha)} |R_m(x+\Delta x, \alpha+\Delta \alpha)|^2.
   \]
  \item[2] For $m = 1, ..., M$, retrieve
  \[
     (\Delta x_{m, i}, \Delta \alpha_{m, i}) = \arg\max_{(\Delta x, \Delta \alpha)\in A(\alpha_i)} |R_m(x_i+\Delta x, \alpha_i+\Delta \alpha)|^2.
  \]
  Let $c_{m, i} \leftarrow R_m(x_i+\Delta x_{m, i}, \alpha_i+\Delta \alpha_{m, i})$, and update
  \[
  R_m(x, \alpha) \leftarrow R_m(x, \alpha)- c_{m, i} \langle B_{x, s, \alpha}, B_{x_i+\Delta x_{m, i}, s, \alpha_i+\Delta \alpha_{m, i}}\rangle.
  \]
  \item[3] Stop if $i = n$, else go back to 1.
\end{itemize}

\subsection{Shared sketch algorithm}

Finally, we come to the shared sketch algorithm that we actually used in the experiments in this paper. The algorithm involves two modifications to Algorithm 2.1. \\

{\bf Algorithm 3: Shared sketch algorithm}
\begin{itemize}
  \item[0] Initialize $i \leftarrow 0$. For $m = 1, ..., M$, initialize $R_m(x, \alpha) \leftarrow \langle \I_m, B_{x, s, \alpha}\rangle$ for all $(x, \alpha)$.
  \item[1] $i \leftarrow i + 1$.
  Select
  \[
   (x_i, \alpha_i) = \arg\max_{x, \alpha} \sum_{m=1}^{M}
    \max_{(\Delta x, \Delta \alpha)\in A(\alpha)} h(|R_m(x+\Delta x, \alpha+\Delta \alpha)|^2).
   \]
  \item[2] For $m = 1, ..., M$, retrieve
  \[
     (\Delta x_{m, i}, \Delta \alpha_{m, i}) = \arg\max_{(\Delta x, \Delta \alpha)\in A(\alpha_i)} |R_m(x_i+\Delta x, \alpha_i+\Delta \alpha)|^2.
  \]
  Let $c_{m, i} \leftarrow R_m(x_i+\Delta x_{m, i}, \alpha_i+\Delta \alpha_{m, i})$, and update $R_m(x, \alpha) \leftarrow$ 0 if
  \[
  {\rm corr}(B_{x, s, \alpha}, B_{x_i+\Delta x_{m, i}, s, \alpha_i+\Delta \alpha_{m, i}})>0.
  \]
  \item[3] Stop if $i = n$, else go back to 1.
\end{itemize}

The two modifications are as follows:

(1) In Step 1, we change $|R_m(x+\Delta x, \alpha+\Delta \alpha)|^2$ to $h(|R_m(x+\Delta x, \alpha+\Delta \alpha)|^2)$ where $h()$ is a sigmoid function, which increases from 0 to a saturation level $\xi$ (default: $\xi = 6$),
\begin{eqnarray}
h(r) = \xi \left[ \frac{2}{1 + e^{-2r/\xi}} - 1 \right].
\label{eq:sigmoid}
\end{eqnarray}
Intuitively, $\sum_{m=1}^{M} \max_{(\Delta x, \Delta \alpha)\in
A(\alpha)} h(|R_m(x+\Delta x, \alpha+\Delta \alpha)|^2)$  can be
considered the sum of the votes from all the images for the
location and orientation $(x, \alpha)$, where each image
contributes $\max_{(\Delta x, \Delta \alpha)\in A(\alpha)}
h(|R_m(x+\Delta x, \alpha+\Delta \alpha)|^2)$. The sigmoid
transformation prevents a small number of images from contributing
very large values. As a result, the selection of $(x, \alpha)$ is
a more ``democratic'' choice than in Algorithm 2, and the selected
element seeks to sketch as many edges in the training images as
possible. In the next section we shall formally justify the use of
sigmoid transformation by a statistical model.

(2) In Step 2, we update $R_m(x, \alpha) \leftarrow$ 0 if $B_{x,
s, \alpha}$ is not orthogonal to $B_{x_i+\Delta x_{m, i}, s,
\alpha_i+\Delta \alpha_{m, i}}$. That is, we enforce the
orthogonality of the basis $\B_{m} = (B_{x_i+\Delta x_{m, i}, s,
\alpha_i + \Delta \alpha_{m, i}}, i = 1, ..., n)$ for each
training image $m$. Our experience with matching pursuit is that
it usually selects elements that have little overlap with each
other. So for computational convenience, we simply enforce that
the selected elements are orthogonal to each other. For two Gabor
wavelets $B_1$ and $B_2$, we define their correlation as ${\rm
corr}(B_1, B_2) = \sum_{\eta_1 = 0}^{1}\sum_{\eta_2 = 0}^{1}
\langle B_{1, \eta_1}, B_{2, \eta_2}\rangle^2$, that is, the sum
of squared inner products between the sine and cosine components
of $B_1$ and $B_2$. In practical implementation, we allow small
correlations between selected elements, that is, we update $R_m(x,
\alpha) \leftarrow 0$ if ${\rm corr}(B_{x, s, \alpha},
B_{x_i+\Delta x_{m, i}, s, \alpha_i+\Delta \alpha_{m,
i}})>\epsilon$ (the default value of $\epsilon = .1$).

\subsection{Statistical modeling of images}

In this subsection we develop a statistical model for $\I_m$. A
statistical model is not only important for justifying Algorithm 3
for learning the active basis template, it also enables us to use
the learned template to recognize the objects in testing images,
because we can use the log-likelihood to score the matching
between the learned template and the image data.

The statistical model is based on the decomposition $\I_m =
\sum_{i=1}^{m} c_{m, i} B_{x_i+\Delta x_{m, i}, s, \alpha_i+\Delta
\alpha_{m, i}} + U_{m}$, where $\B_{m} = (B_{x_i+\Delta x_{m, i},
s, \alpha_i+\Delta \alpha_{m, i}}, i = 1, ..., n)$ is orthogonal,
and $c_{m, i} = \langle \I_m, B_{x_i+\Delta x_{m, i}, s,
\alpha_i+\Delta \alpha_{m, i}}\rangle$, so $U_m$ lives in the
subspace that is orthogonal to $\B_m$. In order to specify a
statistical model for $\I_m$ given $\B_m$, we only need to specify
the distribution of $(c_{m, i}, i = 1, ..., n)$ and the
conditional distribution of $U_m$ given $(c_{m, i}, i = 1, ...,
n)$.

The least squares criterion (\ref{eq:LS}) that drives Algorithm 2
implicitly assumes that $U_{m}$ is white noise, and $c_{m, i}$
follows a flat prior distribution. These assumptions are wrong.
There can be occasional strong edges in the background, but a
white noise $U_m$ cannot account for strong edges. The
distribution of $c_{m, i}$ should be estimated from the training
images, instead of being assumed to be a flat distribution.

In this work we choose to estimate the distribution of $c_{m, i}$
from the training images by fitting an exponential family model to
the sample $\{c_{m, i}, m = 1, ..., M\}$ obtained from the
training images, and we assume that the conditional distribution
of $U_m$ given $(c_{m, i}, i = 1, ..., n)$ is the same as the
corresponding conditional distribution in the natural images. Such
a conditional distribution can account for occasional strong edges
in the background, and it is the use of such a conditional
distribution of $U_m$ as well as the exponential family model for
$c_{m, i}$ that leads to the sigmoid transformation in Algorithm
3. Intuitively, a large response $|R_m(x+\Delta x, \alpha+\Delta
\alpha)|^2$ indicates that there can be an edge at $(x+\Delta x,
\alpha+\Delta \alpha)$. Because an edge can also be accounted for
by the distribution of $U_m$ in the natural images, a large
response should not be taken at its face value for selecting the
basis elements. Instead, it should be discounted by a
transformation such as $h()$ in Algorithm 3.



\subsection{Density substitution and projection pursuit}

More specifically, we adopt the density substitution scheme of
projection pursuit \cite{Friedman} to construct a statistical
model. We start from a reference distribution $q(\I)$. In this
article we assume that $q(\I)$ is the distribution of all the
natural images. We do not need to know $q(\I)$ explicitly beyond
the marginal distribution $q(c)$ of $c = \langle \I, B_{x, s,
\alpha}\rangle$ under $q(\I)$. Because $q(\I)$ is stationary and
isotropic, $q(c)$ is the same for different $(x, \alpha)$. $q(c)$
is a heavy tailed distribution because there are edges in the
natural images. $q(c)$ can be estimated from the natural images by
pooling a histogram of $\{\langle \I, B_{x, s, \alpha}\rangle,
\forall \I, \forall (x, \alpha)\}$ where $\{\I\}$ is a sample of
the natural images.

Given $\B_m = (B_{x_i+\Delta x_{m, i}, s, \alpha_i+\Delta
\alpha_{m, i}}, i = 1, ..., n)$, we modify the reference
distribution $q(\I_m)$ to a new distribution $p(\I_m)$ by changing
the distributions of $c_{m, i}$. Let $p_i(c)$ be the distribution
of $c_{m, i}$ pooled from $\{c_{m, i}, m = 1, ..., M\}$, which are
obtained from the training images $\{\I_m, m = 1, ..., M\}$. Then
we change the distribution of $c_{m, i}$ from $q(c)$ to $p_i(c)$,
for each $i = 1, ..., n$, while keeping the conditional
distribution of $U_{m}$ given $(c_{m, i}, i = 1, ..., n)$
unchanged. This leads us to
\begin{eqnarray}
    p(\I_m \mid \B_m = (B_{x_i+\Delta x_{m, i}, s, \alpha_i+\Delta \alpha_{m, i}}, i = 1, ..., n)) = q(\I_m) \prod_{i=1}^{n} \frac{p_i(c_{m, i})}{q(c_{m, i})}, \label{eq:MODEL}
\end{eqnarray}
where we assume that $(c_{m, i}, i = 1, ..., n)$ are independent
under both $q(\I_m)$ and $p(\I_m \mid \B_m)$, for orthogonal
$\B_m$. The conditional distributions of $U_m$ given $(c_{m, i}, i
= 1, ..., n)$  under $p(\I_m \mid \B_m)$ and $q(\I_m)$ are
canceled out in $p(\I_m \mid \B_m)/q(\I_m)$ because they are the
same. The Jacobians are also the same and are canceled out. So
$p(\I_m \mid \B_m)/q(\I_m) = \prod_{i=1}^{n} {p_i(c_{m,
i})}/{q(c_{m, i})}$.

The following are three perspectives to view  model (\ref{eq:MODEL}):

(1) Classification: we may consider $q(\I)$ as representing the negative examples, and $\{\I_m\}$ are positive examples. We want to find the basis elements $(B_{x_i, s, \alpha_i}, i = 1, ..., n)$ so that the projections $c_{m, i} = \langle \I_m, B_{x_i+\Delta x_{m, i}, s, \alpha_i+\Delta \alpha_{m, i}}\rangle$ for $i = 1, ..., n$ distinguish the positive examples from the negative examples.

(2) Hypothesis testing: we may consider $q(\I)$ as representing the null hypothesis, and the observed histograms of $c_{m, i}, i = 1, ..., n$ are the test statistics that are used to reject the null hypothesis.

(3) Coding: we choose to code $c_{m, i}$ by $p_i(c)$ instead of $q(c)$, while continuing to code $U_m$ by the conditional distribution of $U_m$ given $(c_{m, i}, i = 1, ..., n)$ under $q(\I)$.

For all the three perspectives, we need to choose $B_{x_i, s, \alpha_i}$ so that there is big contrast between $p_i(c)$ and $q(c)$. The shared sketch process can be considered as sequentially flipping dimensions of $q(\I_m)$ from $q(c)$ to $p_i(c)$ to fit the observed images. It is essentially a projection pursuit procedure, with an additional local maximization step for estimating the activities of the basis elements.


\subsection{Exponential tilting and saturation transformation}

While $p_i(c)$ can be estimated from $\{c_{m, i}, m = 1, ..., M\}$ by pooling a histogram, we choose to parametrize $p_i(c)$ with a single parameter so that it can be estimated from even a single image.

We assume $p_i(c)$ to be the following exponential family model:
\begin{eqnarray}
    p(c; \lambda) = \frac{1}{Z(\lambda)} \exp\{\lambda h(r)\} q(c),
    \label{eq:exponential}
\end{eqnarray}
where $\lambda>0$ is the parameter. For $c = (c_0, c_1)$, $r = |c|^2 = c_0^2 + c_1^2$.
\begin{eqnarray*}
Z(\lambda) = \int \exp\{\lambda h(r)\}q(c)dc = \E_{q}[\exp\{\lambda
h(r)\}]
\end{eqnarray*}
is the normalizing constant.  $h(r)$ is a monotone increasing
function. We assume $p_i(c) = p(c; \lambda_i)$, which accounts for
the fact that the squared responses $\{|c_{m, i}|^2 = |\langle
\I_m, B_{x_i+\Delta x_{m, i}, s, \alpha_i+\Delta \alpha_{m,
i}}\rangle|^2, m = 1, ..., M\}$ in the positive examples are in
general larger than those in the natural images, because
$B_{x_i+\Delta x_{m, i}, s, \alpha_i+\Delta \alpha_{m, i}}$ tends
to sketch a local edge segment in each $\I_m$. As mentioned
before, $q(c)$ is estimated by pooling a histogram from the
natural images.

We argue that $h(r)$ should be a saturation transformation in the
sense that as $r \rightarrow \infty$, $h(r)$ approaches a finite
number. The sigmoid transformation in (\ref{eq:sigmoid}) is such a
transformation. The reason for such a transformation is as
follows. Let $q(r)$ be the distribution of $r = |c|^2 =
|\langle\I, B\rangle|^2$ under $q(c)$ where $\I \sim q(\I)$. We
may implicitly model $q(r)$ as a mixture of $p_{\rm on}(r)$ and
$p_{\rm off}(r)$, where $p_{\rm on}$ is the distribution of $r$
when $B$ is on an edge in $\I$, and $p_{\rm off}$ is the
distribution of $r$ when $B$ is not on an edge in $\I$. $p_{\rm
on}(r)$ has a much heavier tail than $p_{\rm off}(r)$. Let $q(r) =
(1-\rho_0)p_{\rm off}(r) + \rho_0 p_{\rm on}(r)$, where $\rho_0$
is the proportion of edges in the natural images. Similarly, let
$p_i(r)$ be the distribution of $r = |c|^2$ under $p_i(c)$. We can
model $p_i(r) = (1-\rho_i)p_{\rm off}(r) + \rho_i p_{\rm on}(r)$,
where $\rho_i > \rho_0$, that is, the proportion of edges sketched
by the selected basis element is higher than the proportion of
edges in the natural images. Then, as $r \rightarrow \infty$,
$p_i(r)/q(r) \rightarrow \rho_i/\rho_0$, which is a constant.
Therefore, $h(r)$ should saturate as $r \rightarrow \infty$.


\subsection{Maximum likelihood learning and pursuit index}

Now we can justify the shared sketch algorithm as a greedy scheme
for maximizing the log-likelihood. With parametrization
(\ref{eq:exponential}) for the statistical model (\ref{eq:MODEL}),
the log-likelihood is
\begin{eqnarray}
\sum_{m=1}^{M} \sum_{i=1}^{n} \log \frac{p_i(c_{m, i})}{q(c_{m, i})} = \sum_{i=1}^{n} \left[\lambda_i \sum_{m=1}^{M} h(|\langle \I_m,  B_{x_i+\Delta x_{m, i}, s, \alpha_i+\Delta \alpha_{m, i}}\rangle|^2) - M \log Z(\lambda_i)\right]. \label{eq:loglik}
\end{eqnarray}
We want to estimate the locations and orientations of the elements
of the active basis, $(x_i, \alpha_i, i = 1, ..., n)$, the
activities of these elements, $(\Delta x_{m, i}, \Delta \alpha_{m,
i}, i = 1, ..., n)$, and the weights, $(\lambda_i, i = 1, ...,
n)$, by maximizing the log-likelihood (\ref{eq:loglik}), subject
to the constraints that $\B_m = (B_{x_i+\Delta x_{m, i}, s,
\alpha_i+\Delta \alpha_{m, i}}$, $i = 1, ..., n)$ is orthogonal
for each $m$.

First, we consider the problem of estimating the weight $\lambda_i$ given $\B_m$. To maximize the log-likelihood (\ref{eq:loglik}) over $\lambda_i$, we only need to maximize
\[
   l_i(\lambda_i) = \lambda_i\sum_{m=1}^{M} h(|\langle \I_m,  B_{x_i+\Delta x_{m, i}, s, \alpha_i+\Delta \alpha_{m, i}}\rangle|^2) - M \log Z(\lambda_i).
\]
By setting $l_i'(\lambda_i) = 0$, we get the well-known form of
the estimating equation for the exponential family model:
\begin{eqnarray}
   \mu(\lambda_i) = \frac{1}{M}\sum_{m=1}^{M} h(|\langle \I_m,  B_{x_i+\Delta x_{m, i}, s, \alpha_i+\Delta \alpha_{m, i}}\rangle|^2),
   \label{eq:MLe}
\end{eqnarray}
where the mean parameter $\mu(\lambda)$ of the exponential family model is
\begin{eqnarray}
   \mu(\lambda) = \E_{\lambda}[h(r)] = \frac{1}{Z(\lambda)} \int
h(r) \exp\{\lambda h(r)\} q(r)dr. \label{eq:ee}
\end{eqnarray}
The estimating equation (\ref{eq:MLe}) can be solved easily, because $\mu(\lambda)$ is a one-dimensional function. We can simply store this monotone function over a one-dimensional grid. Then we solve this equation by looking up the stored values, with the help of nearest neighbor linear interpolation for the values between the grid points. For each grid point of $\lambda$, $\mu(\lambda)$ can be computed by one-dimensional integration as in (\ref{eq:ee}). Thanks to the independence assumption, we only need to deal with such one-dimensional functions, which  relieves us from time consuming MCMC computations.

Next let us consider the problem of selecting $(x_i, \alpha_i)$, and estimating the activity $(\Delta x_{m, i}, \Delta \alpha_{m, i})$ for each image $\I_m$.  Let $\hat{\lambda}_i$ be the solution to the estimating equation (\ref{eq:MLe}). $l_i(\hat{\lambda}_i)$ is monotone in $\sum_{m=1}^{M} h(|\langle \I_m,  B_{x_i+\Delta x_{m, i}, s, \alpha_i+\Delta \alpha_{m, i}}\rangle|^2)$. Therefore, we need to find $(x_i, \alpha_i)$, and $(\Delta x_{m, i}, \Delta \alpha_{m, i})$, by maximizing $\sum_{m=1}^{M} h(|\langle \I_m,  B_{x_i+\Delta x_{m, i}, s, \alpha_i+\Delta \alpha_{m, i}}\rangle|^2)$. This justifies Step 1 of Algorithm 3, where $\sum_{m=1}^{M} h(|R_m(x+\Delta x, \alpha+\Delta \alpha)|^2)$ serves as the pursuit index.

\subsection{SUM-MAX maps for template matching  \label{sect:sum-max}}

After learning the active basis model, in particular, the basis elements $\B = (B_{x_i, s, \alpha_i}, i = 1, ..., n)$ and the weights $(\lambda_i, i = 1, ..., n)$, we can use the learned model to find the object in a testing image $\I$, as illustrated by Figure (\ref{fig:find}). The testing image may not be defined on the same lattice as the training images. For example, the testing image may be larger than the training images. We assume that there is one object in the testing image, but we do not know the location of the object in the testing image. In order to detect the object, we scan the template over the testing image, and at each location $x$, we can deform the template and match it to the image patch around $x$. This gives us a log-likelihood score at each location $x$. Then we can find the maximum likelihood location $\hat{x}$ that achieves the maximum of the log-likelihood score among all the $x$. After computing $\hat{x}$, we can then retrieve the activities of the elements of the active basis template centered at $\hat{x}$. \\

{\bf Algorithm 4: Object detection by template matching}
\begin{itemize}
\item[1] For every $x$, compute
\[
    l(x) = \sum_{i=1}^{n} \left[ \lambda_i \max_{(\Delta x, \Delta \alpha) \in A(\alpha_i)} h(|\langle \I, B_{x+x_i+\Delta x, s, \alpha_i + \Delta \alpha}\rangle|^2) - \log Z(\lambda_i) \right].
\]

\item[2] Select $\hat{x} = \arg\max_{x} l(x)$. For $i = 1, ..., n$, retrieve
    \[
       (\Delta x_i, \Delta \alpha_i) = \arg\max_{(\Delta x, \Delta \alpha) \in A(\alpha_i)} |\langle \I, B_{\hat{x}+x_i+\Delta x, s, \alpha_i + \Delta \alpha}\rangle|^2.
    \]

\item[3] Return the location $\hat{x}$, and the deformed template $(B_{\hat{x}+x_i+\Delta x_i, s, \alpha_i + \Delta \alpha_i}, i = 1, ..., n)$.
\end{itemize}

Figure (\ref{fig:find}) displays the deformed template $(B_{\hat{x}+x_i+\Delta x_i, s, \alpha_i + \Delta \alpha_i}, i = 1, ..., n)$, which is superposed on the image on the right.

Step 1 of the above algorithm can be realized by a computational architecture called sum-max maps. \\

{\bf Algorithm 4.1: sum-max maps}
\begin{itemize}
\item[1] For all $(x, \alpha)$, compute $\S1(x, \alpha) = h(|\langle \I, B_{x, s, \alpha}\rangle|^2)$.
\item[2] For all $(x, \alpha)$, compute $\M1(x, \alpha) = \max_{(\Delta x, \Delta \alpha) \in A(\alpha)}\S1(x+\Delta x, \alpha+\Delta \alpha)$.
\item[3] For all $x$, compute $\S2(x) = \sum_{i=1}^{n} [\lambda_i \M1(x+x_i, \alpha_i) - \log Z(\lambda_i)]$.
\end{itemize}
$\S2(x)$ is $l(x)$ in Algorithm 4.


The local maximization operation in Step 2 of Algorithm 4.1 has been hypothesized as the function of the complex cells of the primary visual cortex \cite{Riesenhuber_Poggio}. In the context of the active basis model, this operation can be justified as the maximum likelihood estimation of the activities. The shared sketch learning algorithm can also be written in terms of sum-max maps.

The activities $(\Delta x_{m, i}, \Delta \alpha_{m, i}, i = 1, ..., n)$ should be treated as latent variables in the active basis model. However, in both learning and inference algorithms, we treat them as unknown parameters, and we maximize over them instead of integrating them out. According to Little and Rubin \cite{Little}, maximizing the complete-data likelihood over the latent variables may not lead to valid inference in general. However, in natural images, there is little noise, and the uncertainty in the activities is often very small. So maximizing over the latent variables can be considered a good approximation to integrating out the latent variables.

\section{Learning active basis templates by EM-type algorithms}

The shared sketch algorithm in the previous section requires that
the objects in the training images $\{\I_m\}$ are of the same pose,
at the same location and scale, and the lattice of $\I_m$ is the
bounding box of the object in $\I_m$. It is often the case that the
objects may appear at different unknown locations, orientations and scales in $\{\I_m\}$. The unknown locations, orientations and scales can be incorporated into the image generation process as hidden variables. The template can still be learned by the maximum likelihood method.

\subsection{Learning with unknown orientations}

We start from a simple example of learning a horse template at the
side view, where the horses can face either to the left or to the
right. Figure (\ref{fig:horse}) displays the results of EM
learning. The three templates in the first row are the learned
templates in the first three iterations of the EM algorithm. The
rest of the figure displays the training images, and for each
training image, a deformed template is displayed to the right of
it. The EM algorithm correctly estimates the direction for each
horse, as can be seen by how the algorithm flips the template to
sketch each training image.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{*{3}{l}}
\psfig{figure=ex_flip_horse_template1.eps,height=.4in} &
\psfig{figure=ex_flip_horse_template2.eps,height=.4in} &
\psfig{figure=ex_flip_horse_template3.eps,height=.4in}
\\
\multicolumn{3}{c} { \psfig{figure=ex_flip_horse_canonical.eps,
width= 4in}
} \\
\multicolumn{3}{c} { \psfig{figure=ex_flip_horse_reflected.eps,
width= 4in} }
\end{tabular}
\end{center}
\caption{Template learned from images of horses facing two
different directions.  The first row displays the templates
learned in the first 3 iterations of the EM algorithm. For each
training image, the deformed template is plotted to the right of
it. The number of training images is $57$. The image size is
$120\times 100$ (width $\times$ height). The number of elements is
$40$. The number of EM iterations is $3$.  } \label{fig:horse}
\end{figure}


Let $\B = (B_{i} = B_{x_i, s, \alpha_i}, i = 1, ..., n)$ be the
deformable template of the horse, and $\B_m =$ $(B_{m, i} =
B_{x_i+\Delta x_{m, i}, s, \alpha_i + \Delta \alpha_{m, i}}$, $i =
1, ..., n)$ be the deformed template for $\I_m$. Then $\I_m$ can
either be generated by $\B_m$ or the mirror reflection of $\B_m$,
that is,  $(B_{R(x_i+\Delta x_{m, i}), s, -(\alpha_i + \Delta
\alpha_{m, i})}, i = 1, ..., n)$, where for $x = (x_1, x_2)$,
$R(x) = (-x_1, x_2)$ (we assume that the template is centered at
origin).  We can introduce a hidden variable $z_m$ to account for
this uncertainty, so that $z_m = 1$ if $\I_m$ is generated by
$\B_m$, and $z_m = 0$ if $\I_m$ is generated by the mirror
reflection of $\B_m$. More formally, we can define $\B_{m}(z_m)$,
so that $\B_m(1) = \B_m$, and $\B_m(0)$ is the mirror reflection
of $\B_m$. Then we can assume the following mixture model: $z_m
\sim {\rm Bernoulli}(\rho)$, where $\rho$ is the prior probability
that $z_m = 1$, and $[\I_m \mid z_m] \sim p(\I_m \mid \B_m(z_m),
\Lambda)$, where $\Lambda = (\lambda_i, i = 1, ..., n)$. We need
to learn $\B$, and estimate $\Lambda$ and $\rho$.

A simple observation is that $p(\I_m \mid \B_{m}(z_m)) =
p(\I_m(z_m) \mid \B_m)$, where $\I_m(1) = \I_m$ and $\I_m(0)$ is
the mirror reflection of $\I_m$. In other words, in the case of
$z_m = 1$, we do not need to make any change to $\I_m$ or $\B_m$.
In the case of $z_m = 0$, we can either flip the template or flip
the image, and these two alternatives will produce the same value
for the likelihood function.

In the EM algorithm, the E-step imputes $z_m$ for $m = 1, ..., M$
using the current template $\B$. This means recognizing the
orientation of the object in $\I_m$. Given $z_m$, we can change
$\I_m$ to $\I_m(z_m)$, so that $\{\I_m(z_m)\}$ become aligned with
each other, if $\{z_m\}$ are imputed correctly. Then in the
M-step, we can learn the template from the aligned images
$\{\I_m(z_m)\}$ by the shared sketch algorithm.

The complete data log-likelihood for the $m$th observation is
\begin{eqnarray*}
   \log p(\I_m, z_m \mid \B_m) &=&
   z_m \log p(\I_m \mid \B_m, \Lambda) + (1-z_m) \log p(\I_m (0) \mid
   \B_m, \Lambda)  \nonumber\\
   &+& z_{m} \log \rho + (1-z_m) \log (1-\rho),
\end{eqnarray*}
which is linear in $z_m$. So in the E-step, we only need to compute
the predictive expectation of $z_m$,
\begin{eqnarray*}
   \hat{z}_m &=& \Pr(z_m = 1 \mid \B_m, \Lambda, \rho) \nonumber \\
   &=& \frac{\rho p(\I_m \mid \B_m, \Lambda)}
   {\rho p(\I_m \mid \B_m, \Lambda) + (1-\rho) p(\I_m(0) \mid \B_m,
   \Lambda)}.
\end{eqnarray*}
Both $\log p(\I_m\mid \B_m, \Lambda)$ and $\log p(\I_m(0) \mid \B_m,
\Lambda)$ are readily available in the M-step.

The M-step seeks to maximize the expectation of the complete-data
log-likelihood,
\begin{eqnarray}
&& \sum_{i=1}^{n} \left[ {\lambda}_i \sum_{m=1}^{M}
\left(\hat{z}_{m} h(|\langle \I_m, B_{m,
i}\rangle|^2)+(1-\hat{z}_{m})h(|\langle \I_m(0), B_{m,
i}\rangle|^2)\right) - M \log Z({\lambda}_i)\right ]
\label{eq:L1}\\&+& \left[\log \rho \sum_{m=1}^{M} \hat{z}_m + \log
(1-\rho) (M - \sum_{m=1}^{M} \hat{z}_m)\right]. \label{eq:L2}
\end{eqnarray}
The maximization of (\ref{eq:L2}) leads to $\hat{\rho} =
\sum_{m=1}^{M} \hat{z}_m/M$. The maximization of (\ref{eq:L1}) can
be accomplished by the shared sketch algorithm, that is, Algorithm
3, with the following minor modifications:


(1) The training images become $\{\I_{m}, \I_{m}(0), m = 1, ..,
M\}$, that is, there are $2M$ training images instead of $M$ images.
Each $\I_m$ contributes two copies, the original copy $\I_m$ or
$\I_m(1)$, and the mirror reflection $\I_m(0)$. This reflects the
uncertainty in $z_m$. For each image $\I_m$, we attach a weight
$\hat{z}_m$ to $\I_m$, and a weight $1-\hat{z}_m$ to $\I_m(0)$.
Intuitively, a fraction of the horse in $\I_m$ is at the same orientation as the current template, and a fraction of it is at the opposite orientation --- a ``Schrodinger horse'' so to speak. We use $(\J_{k}, w_k, k = 1, ..., 2M)$ to represent these $2M$ images and their weights.

(2) In Step 1 of the shared sketch algorithm, we select $(x_i, \alpha_i)$ by
\begin{eqnarray*}
   (x_i, \alpha_i) = \arg\max_{x, \alpha} \sum_{k=1}^{2M} w_{k}
    \max_{(\Delta x, \Delta \alpha)\in A(\alpha)} h(|R_k(x+\Delta x, \alpha+\Delta \alpha)|^2).
\end{eqnarray*}

(3) The maximum likelihood estimating equation for
$\lambda_i$ is
\begin{eqnarray*}
   \mu(\lambda_i) = \frac{1}{M} \sum_{k=1}^{2M} w_k
    \max_{(\Delta x, \Delta \alpha)\in A(\alpha_i)} h(|R_k(x_i+\Delta x, \alpha_i+\Delta \alpha)|^2),
\end{eqnarray*}
where the right-hand side is the weighted average obtained from
the $2M$ training images.

(4) Along with the selection of $B_{x_i, s, \alpha_i}$ and the
estimation of $\lambda_i$, we should calculate the template
matching scores
\[
\log p(\J_k \mid \B_k, \Lambda) =
\sum_{i=1}^{n} \left[\hat{\lambda}_i \max_{(\Delta x, \Delta \alpha)\in A(\alpha_i)}h(|R_k(x_i+\Delta x, \alpha_i+\Delta \alpha)|^2) - \log Z(\hat{\lambda}_i)\right],
 \]
for $k = 1, ..., 2M$. This gives us $\log p(\I_m\mid \B_m,
\Lambda)$ and $\log p(\I_m(0) \mid \B_m, \Lambda)$, which can then
be used in the E-step.

We initialize the algorithm by randomly generating $\hat{z}_{m} \sim
{\rm Unif}[0, 1]$, and then iterate between the M-step and the
E-step. We stop the algorithm after a few iterations. Then we
estimate  $z_m = 1$ if $\hat{z}_{m}>1/2$ and $z_m = 0$ otherwise.

In Figure (\ref{fig:horse}) the results are obtained after 3
iterations of the EM algorithm. Initially, the learned template is
quite symmetric, reflecting the confusion of the algorithm
regarding the directions of the horses. Then the EM algorithm
begins a process of ``symmetry breaking'' or ``polarization.'' The
slight asymmetry in the initial template will push the algorithm
toward favoring for each image the direction that is consistent
with the majority direction. This process quickly leads to all the
images aligned to one common direction.

Figure (\ref{fig:pigeon}) shows another example where a template
of a pigeon is learned from examples with mixed directions.


\begin{figure}[!htb]
\begin{center}
\begin{tabular}{{l|lll}}
%\psfig{figure=ex_flip_p_template1.eps,height=.5in} &
%\psfig{figure=ex_flip_p_template2.eps,height=.5in} &
\psfig{figure=ex_flip_p_template3.eps,height=.4in}
&
\multicolumn{3}{l} { \psfig{figure=ex_flip_p_canonical.eps, height=
0.4in}
} \\ %\hline
&
\multicolumn{3}{l} { \psfig{figure=ex_flip_p_reflected.eps, height=
0.4in} }
\end{tabular}
\end{center}
\caption{Template learned from $11$ images of pigeons facing
different directions. The image size is $150\times 150$. The
number of elements is $50$. The number of iterations is $3$.  }
\label{fig:pigeon}
\end{figure}



We can also learn a common template when the objects are at more
than two different orientations in the training images. The
algorithm is essentially the same as described above. Figure
(\ref{fig:cap}) displays the learning of the template of a
baseball cap from examples where the caps turn to different
orientations. The E-step involves rotating the images by matching
to the current template, and the M-step learns the template from
the rotated images.


\begin{figure}[!htb]
\begin{center}
\begin{tabular}{l|l}
\psfig{figure=ex_trans_cap_sym.eps,height=1.2in} &
\psfig{figure=ex_trans_cap_results.eps,height=1.2in}
\end{tabular}
\end{center}
\caption{Template learned from 15 images of baseball caps facing
different orientations. The image size is $100\times 100$. The
number of elements is $40$. The number of iterations is $5$.  }
\label{fig:cap}
\end{figure}



\subsection{Learning from nonaligned images}

When the objects appear at different locations in the training
images, we need to infer the unknown locations while learning the
template. Figure (\ref{fig:bike}) displays the template of a bike
learned from the 7 training images where the objects appear at
different locations and are not aligned. It also displays the
deformed templates superposed on the objects in the training
images.

\begin{figure}[!htb]
\begin{tabular}{l}
\psfig{figure=ex5c_bike_sym0.eps,height=.6in}
\psfig{figure=ex5c_bike_sym1.eps,height=.6in}
\psfig{figure=ex5c_bike_sym2.eps,height=.6in}
\psfig{figure=ex5c_bike_sym3.eps,height=.6in} \\
\hline\\
\psfig{figure=ex5c_bike1.eps,height=.6in}
\psfig{figure=ex5c_bike2.eps,height=.6in}
\psfig{figure=ex5c_bike3.eps,height=.6in}
\psfig{figure=ex5c_bike4.eps,height=.6in}
\psfig{figure=ex5c_bike5.eps,height=.6in}
\psfig{figure=ex5c_bike6.eps,height=.6in}
\psfig{figure=ex5c_bike7.eps,height=.6in}
\end{tabular}
\caption{The first row shows the sequence of templates learned in
the first 3 iterations. The first one is the starting template,
which is learned from the first training image. The second row
shows the bikes detected by the learned template, where a deformed
template is superposed on each training image. The size of the
template is $225 \times 169$. The number of elements is $50$. The
number of iterations is $3$.} \label{fig:bike}
\end{figure}

In order to incorporate the unknown locations into the image
generation process, let us assume that both the learned template
$\B = (B_{x_i, s, \alpha_i}, i = 1, ..., n)$ and the training
images $\{\I_m\}$ are centered at origin. Then let us assume that
the location of the object in image $\I_m$ is $x^{(m)}$, which is
assumed to be uniformly distributed within the image lattice of
$\I_m$. Let us define $\B_m(x^{(m)}) = (B_{x^{(m)} + x_i+\Delta
x_{m, i}, s, \alpha_i + \Delta \alpha_{m, i}}, i = 1, ..., n)$ to
be the deformed template obtained by translating the template $\B$
from the origin to $x^{(m)}$ and then deforming it. Then the
generative model for $\I_m$ is $p(\I_m | \B_m(x^{(m)}), \Lambda)$.

Just like the example of learning the horse template, we can
transfer the transformation of the template to the transformation
of the image data, and the latter transformation leads to the
alignment of the images. Let us define $\I_m(x^{(m)})$ to be the
image obtained by translating the image $\I_m$ so that the center
of $\I_m(x^{(m)})$ is $-x^{(m)}$. Then $p(\I_m | \B_m(x^{(m)}),
\Lambda) = p(\I_m(x^{(m)})| \B_m, \Lambda)$. If we know $x^{(m)}$
for $m = 1, ..., M$, then the images $\{\I_m(x^{(m)})\}$ are all
aligned, so that we can learn a template from these aligned images
by the shared sketch algorithm. On the other hand, if we know the
template, we can use the template to recognize and locate the
object in each $\I_m$ by the inference algorithm, that is,
Algorithm 4, using the sum-max maps, and identify $x^{(m)}$. Such
considerations naturally lead to the iterative EM-type scheme.

The complete-data log-likelihood  is
\begin{eqnarray}
\sum_{i=1}^{n} \left[ {\lambda}_i \sum_{m=1}^{M} h(|\langle
\I_m(x^{(m)}), B_{x_i+\Delta x_{m, i}, s, \alpha_i + \Delta \alpha_{m, i}}\rangle|^2) - M \log Z({\lambda}_i)\right
]. \label{eq:La}
\end{eqnarray}

In the E-step, we perform the recognition task by calculating
\begin{eqnarray*}
  p_m(x) = \Pr(x^{(m)} = x \mid \B, \Lambda)  \propto p(\I_m(x) \mid \B_m, \Lambda), \; \forall x.
\end{eqnarray*}
That is, we scan the template over the whole image $\I_m$, and at
each location $x$, we evaluate the template matching between
the image $\I_m$ and the translated and deformed template $\B_m(x)$. $\log p(\I_m(x) \mid \B_m, \Lambda)$ is the $\S2(x)$ output by the sum-max maps in Algorithm 4.1. This gives us $p_m(x)$, which is the posterior or
predictive distribution of the unknown location $x^{(m)}$ within
the image lattice of $\I_m$. We can then use $p_m(x)$ to compute
the expectation of the complete-data log-likelihood (\ref{eq:La}) in
the E-step.

Our experience suggests that $p_m(x)$ is always highly peaked at
a particular position. So instead of computing the average of
(\ref{eq:La}), we simply impute $x^{(m)} = \arg\max_{x} p_m(x)$.

Then in the M-step, we maximize the complete data log-likelihood
(\ref{eq:La}) by the shared sketch algorithm, that is, we learn
the template $\B$ from $\{\I_m(x^{(m)})\}$. This step performs
supervised learning from the aligned images.

In our current experiment we initialize the algorithm by learning
$(\B, \Lambda)$ from the first image. In learning from this single
image, we set $b_1 = b_2 = 0$, that is, we do not allow the
elements $(B_{x_i, s, \alpha_i}, i = 1, ..., n)$ to perturb. After
that, we reset $b_1$ and $b_2$ to their default values, and
iterate the recognition step and the supervised learning step.

In addition to the unknown locations, we also allow the uncertainty
in scales. In the recognition step, for each $\I_m$, we search over a number of different resolutions of $\I_m$. We take $\I_m(x^{(m)})$ to be the optimal resolution that contains the maximum template matching score across all the resolutions.

In Figure (\ref{fig:bike}) the first row displays the templates
learned over the EM iterations. The first template is learned from
the first training image. Figures (\ref{fig:camel}) to
(\ref{fig:cat}) display more examples.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{l}
\psfig{figure=ex5c_camel_sym0.eps,height=.6in}
\psfig{figure=ex5c_camel_sym1.eps,height=.6in}
\psfig{figure=ex5c_camel_sym3.eps,height=.6in}
\psfig{figure=ex5c_camel_sym5.eps,height=.6in}\\
\hline\\
\psfig{figure=ex5c_camel1.eps,height=.6in}
\psfig{figure=ex5c_camel2.eps,height=.6in}
\psfig{figure=ex5c_camel6.eps,height=.6in}
\psfig{figure=ex5c_camel4.eps,height=.6in}
\psfig{figure=ex5c_camel9.eps,height=.6in}
\psfig{figure=ex5c_camel5.eps,height=.6in}\\
\psfig{figure=ex5c_camel3.eps,height=.6in}
\psfig{figure=ex5c_camel7.eps,height=.6in}
\psfig{figure=ex5c_camel8.eps,height=.6in}
\psfig{figure=ex5c_camel10.eps,height=.6in}
\psfig{figure=ex5c_camel11.eps,height=.6in}
\end{tabular}
\end{center}
\caption{The first row shows the sequence of templates learned in
iterations 0, 1, 3, 5. The second and third rows show the camel
images with superposed deformed templates. The size of the
template is $192\times 145$. The number of elements is $60$. The
number of iterations is $5$.} \label{fig:camel}
\end{figure}




\begin{figure}[!htb]
\begin{center}
\begin{tabular}{l}
\psfig{figure=ex5c_crane_sym0.eps,height=.6in}
\psfig{figure=ex5c_crane_sym1.eps,height=.6in}
\psfig{figure=ex5c_crane_sym3.eps,height=.6in}
\psfig{figure=ex5c_crane_sym5.eps,height=.6in}
\\
\hline
\\
\psfig{figure=ex5c_crane1.eps,height=.6in}
\psfig{figure=ex5c_crane2.eps,height=.6in}
\psfig{figure=ex5c_crane3.eps,height=.6in}
\psfig{figure=ex5c_crane4.eps,height=.6in}
\psfig{figure=ex5c_crane5.eps,height=.6in}\\
\psfig{figure=ex5c_crane6.eps,height=.6in}
\psfig{figure=ex5c_crane7.eps,height=.6in}
\psfig{figure=ex5c_crane8.eps,height=.6in}
\psfig{figure=ex5c_crane9.eps,height=.6in}
\psfig{figure=ex5c_crane10.eps,height=.6in}
\end{tabular}
\end{center}
\caption{The first row shows the sequence of templates learned in
iterations 0, 1, 3, 5. The other rows show the crane images with
superposed deformed templates. The size of the template is
$285\times 190$. The number of elements is $50$. The number of
iterations is $5$.} \label{fig:crane}
\end{figure}


\begin{figure}[!htb]
\begin{center}
\begin{tabular}{l}
\psfig{figure=ex5c_Whorse_sym0.eps,height=.5in}
\psfig{figure=ex5c_Whorse_sym2.eps,height=.5in}
\psfig{figure=ex5c_Whorse_sym4.eps,height=.5in}
\psfig{figure=ex5c_Whorse_sym6.eps,height=.5in}
\psfig{figure=ex5c_Whorse_sym8.eps,height=.5in}
\psfig{figure=ex5c_Whorse_sym10.eps,height=.5in}
\\
\hline\\
\psfig{figure=ex5c_Whorse1.eps,height=.5in}
\psfig{figure=ex5c_Whorse2.eps,height=.5in}
\psfig{figure=ex5c_Whorse3.eps,height=.5in}
\psfig{figure=ex5c_Whorse4.eps,height=.5in}
\psfig{figure=ex5c_Whorse5.eps,height=.5in}
\psfig{figure=ex5c_Whorse6.eps,height=.5in}
\psfig{figure=ex5c_Whorse7.eps,height=.5in}\\

\psfig{figure=ex5c_Whorse8.eps,height=.5in}
\psfig{figure=ex5c_Whorse9.eps,height=.5in}
\psfig{figure=ex5c_Whorse10.eps,height=.5in}
\psfig{figure=ex5c_Whorse11.eps,height=.5in}
\psfig{figure=ex5c_Whorse12.eps,height=.5in}
\psfig{figure=ex5c_Whorse13.eps,height=.5in}
\psfig{figure=ex5c_Whorse14.eps,height=.5in}\\
\psfig{figure=ex5c_Whorse15.eps,height=.5in}
\psfig{figure=ex5c_Whorse16.eps,height=.5in}
\psfig{figure=ex5c_Whorse17.eps,height=.5in}
\psfig{figure=ex5c_Whorse18.eps,height=.5in}
\psfig{figure=ex5c_Whorse19.eps,height=.5in}
\psfig{figure=ex5c_Whorse20.eps,height=.5in}
\end{tabular}
\end{center}
\caption{The first row shows the sequence of templates learned in
iterations 0, 2, 4, 6, 8, 10. The other rows show the horse images
with superposed deformed templates. We use the first 20 images of
the Weizmann horse data set \cite{Borenstein}, which are resized
to half the original sizes. The size of the template is 158
$\times$ 116. The number of elements is $60$. The number of
iterations is $10$. The detection results on the rest of the
images in this data set can be found in the reproducibility page.
} \label{fig:Whorse}
\end{figure}


\begin{figure}[!htb]
\begin{center}
\begin{tabular}{l}
\psfig{figure=ex5c_kitty_sym0.eps,height=.6in}
\psfig{figure=ex5c_kitty_sym1.eps,height=.6in}
\psfig{figure=ex5c_kitty_sym2.eps,height=.6in}
\psfig{figure=ex5c_kitty_sym3.eps,height=.6in}
\\
\hline\\
\psfig{figure=ex5c_kitty1.eps,height=.6in}
\psfig{figure=ex5c_kitty2.eps,height=.6in}
\psfig{figure=ex5c_kitty3.eps,height=.6in}
\psfig{figure=ex5c_kitty4.eps,height=.6in}
\psfig{figure=ex5c_kitty5.eps,height=.6in}\\
\psfig{figure=ex5c_kitty6.eps,height=.6in}
\psfig{figure=ex5c_kitty7.eps,height=.6in}
\psfig{figure=ex5c_kitty8.eps,height=.6in}
\psfig{figure=ex5c_kitty9.eps,height=.6in}
\end{tabular}
\end{center}
\caption{The size of the template is $240\times 180$. The number of elements is $60$. The number of iterations is $3$.} \label{fig:cat}
\end{figure}





\section{Discussion}

This paper experiments with EM-type algorithms for learning active
basis models from training images where the objects may appear at
unknown locations, orientations and scales. For more details on
implementing the shared sketch algorithm, the reader is referred
to \cite{Wu1} and the source code posted on the reproducibility
page.

We would like to emphasize two aspects of the algorithms that are
different from the usual EM algorithm. The first aspect is that
the M-step involves the selection of the basis elements, in
addition to the estimation of the associated parameters. The
second aspect is that the performance of the algorithms can rely
heavily on the initializations. In learning from nonaligned
images, the algorithm is initialized by training the active basis
model on a single image. Because of the simplicity of the model,
it is possible to learn the model from a single image. In
addition, the learning algorithm seems to converge within a few
iterations.

\subsection{Limitations}

The active basis model is a simple extension of the wavelet
representation. It is still very limited in the following aspects.
The model cannot account for large deformations, articulate
shapes, big changes in poses and view points, and occlusions. The
current form of the model does not describe textures and lighting
variations either. The current version of the learning algorithm
only deal with situations where there is one object in each image.
Also, we have tuned two parameters in our implementation. One is
the image resize factor that we apply to the training images
before the model is learned. Of course, for each experiment, a
single resize factor is applied to all the training images. The
other parameter is the number of elements in the active basis.

\subsection{Possible extensions}

It is possible to extend the active basis model to address some of
the above limitations. We shall discuss two directions of
extensions. One is to use active basis models as parts of the
objects. The other is to train active basis models by local
learning.

{\em Active basis models as part-templates:} The active basis
model is a composition of a number of Gabor wavelet elements. We
can further compose multiple active basis models to represent more
articulate shapes or to account for large deformations by allowing
these active basis models to change their overall locations,
scales and orientations within limited ranges. These active basis
models serve as part-templates of the whole composite template.
This is essentially a hierarchical recursive compositional
structure \cite{Geman_composition, Zhu_grammar}. The inference or
template matching can be based on a recursive structure of sum-max
maps. Learning such a structure should be possible by extending
the learning algorithms studied in this article. See \cite{Wu1}
for preliminary results. See also \cite{Ross, Sudderth} for recent
work on part-based models.

{\em Local learning of multiple prototype templates:} In each
experiment we assume that all the training images share a common
template. In reality, the training images may contain different
types of objects, or different poses of the same type of objects.
It is therefore necessary to learn multiple prototype templates.
It is possible to do so by modifying the current learning
algorithm. After initializing the algorithm by single image
training, in the M-step, we can relearn the template only from the
$K$ images with the highest template matching scores, that is, we
relearn the template from the $K$ nearest neighbors of the current
template. Such a scheme is consistent with the EM-clustering
algorithm for fitting mixture models. We can start the algorithm
from every training image, so that we learn a local prototype
template around each training image. Then we can trim and merge
these prototypes. See \cite{Wu1} for preliminary results.

\section*{Reproducibility}
All the experimental results reported in this paper can be
reproduced by the code that we have posted at
http://www.stat.ucla.edu/$\sim$ywu/ActiveBasis.html.


\section*{Acknowledgment}
We are grateful to the editor of the special issue and the two
reviewers for their valuable comments that have helped us improve
the presentation of the paper. The work reported in this paper has
been supported by NSF-DMS 0707055, NSF-IIS 0713652, ONR
N00014-05-01-0543, Air Force grant FA 9550-08-1-0489 and the Keck
foundation.


{\small
\begin{thebibliography}{9}

\bibitem{Amit}
Amit, Y. and Trouve, A. (2007) POP: patchwork of parts models for
object recognition. \emph{International Journal of Computer
Vision}, {\bf 75} 267-282.


\bibitem{Bell_Sejnowski}
Bell, A. and Sejnowski, T. J. (1997) The 'independent components' of
natural scenes are edge filters.  \emph{Vision Research}, {\bf 37}
3327-3338.

\bibitem{Borenstein}
Borenstein, E., Sharon, E. and Ullman, S. (2004) Combining
top-down and bottom-up segmentation. \emph{IEEE CVPR workshop on
Perceptual Organization in Computer Vision}, vol. 4, pp. 46--46.


\bibitem{Candes_Donoho}
Candes, E. J. and Donoho, D. L. (1999) Curvelets - a surprisingly
effective nonadaptive representation for objects with edges. {\em
Curves and Surfaces}, L. L. Schumakeretal. (eds).


\bibitem{Chen_basis_pursuit}
Chen, S., Donoho, D. and Saunders, M. A. (1999) Atomic decomposition by basis pursuit. {\em SIAM Journal on Scientific Computing}, {\bf 20} 33-61.

\bibitem{Daugman}
Daugman, J. (1985) Uncertainty relation for resolution in space,
spatial frequency, and orientation optimized by two-dimensional
visual cortical filters. {\em Journal of Optical Society of
America A: Optics, Image Science, and Vision}, {\bf 2} 1160-1169.

\bibitem{DLR}
Dempster, A. P., Laird, N. M. and Rubin, D. B. (1977) Maximum likelihood from incomplete data via the EM algorithm. {\em Journal of the Royal Statistical Society, B}, {\bf 39} 1-38.

\bibitem{Donoho_compression}
Donoho, D. L., Vetterli, M., DeVore, R. A. and Daubechie, I. (1998) Data compression and harmonic analysis. {\em IEEE Transactions on
Information Theory}, {\bf 6} 2435-2476.

\bibitem{Perona}
Fergus, R., Perona, P. and  Zisserman,  A. (2003)  Object class
recognition by unsupervised scale-invariant learning. {\em
Proceedings of Computer Vision and Pattern Recognition}, vol. 2,
pp. 264--271.


\bibitem{Friedman}
Friedman, J. H. (1987) Exploratory projection pursuit. {\em Journal
of the American Statistical Association}, {\bf 82} 249-266.

\bibitem{Geman_composition} Geman, S., Potter, D. F. and Chi, Z. (2002) Composition systems. \emph{Quartarly of Applied Mathematics}, {\bf 60} 707-736.


\bibitem{Little}
Little, R. J. A. and  Rubin, D. B. (1983) On jointly estimating parameters and missing data by maximizing the complete data likelihood. \emph{The American Statistician}, {\bf 37} 218-220.


\bibitem{Mallat_Zhang}
Mallat, S. and Zhang, Z. (1993)  Matching pursuit in a
time-frequency dictionary. {\em IEEE Transactions on Signal
Processing}, {\bf 41} 3397-415.


\bibitem{Meng}
Meng, X.-L. and van Dyk, D. (1997) The EM algorithm - an old folk-song sung to a fast new tune. {\em Journal of the Royal Statistical Society, B}, {\bf 59} 511-567.


\bibitem{Olshausen_Field}
Olshausen, B. A. and Field, D. J. (1996) Emergence of simple-cell
receptive field properties by learning a sparse code for natural
images. {\em Nature}, {\bf 381} 607-609.

\bibitem{Rabiner}
Rabiner, L. R. (1989) A tutorial on hidden Markov models and selected applications in speech recognition. {\em Proceedings of the IEEE},
{\bf 77} 257-286.


\bibitem{Riesenhuber_Poggio}
Riesenhuber, M. and Poggio, T. (1999) Hierarchical models of object
recognition in cortex. {\em Nature Neuroscience},  {\bf 2} 1019-1025.

\bibitem{Ross}
Ross, D. A. and Zemel, R. S. (2006) Learning parts-based representations of data. \emph{Journal of Machine Learning Research}, {\bf 7} 2369-2397.

\bibitem{Simoncelli_shiftable}
Simoncelli, E. P., Freeman, W. T., Adelson, E. H. and Heeger, D. J. (1992) Shiftable multiscale transforms. {\em IEEE Transactions on
Infomation Theory}, {\bf 38} 587-607.

\bibitem{Sudderth}
Sudderth, E. B., Torralba, A., Freeman, W. T. and Willsky, A. S. (2005) Learning hierarchical models of scenes, objects, and parts. {\em Proceedings of International Conference of Computer Vision}, {\bf 2} 1331-1338.

\bibitem{lasso}
Tibshirani, R. (1996) Regression shrinkage and selection via the lasso. {\em Journal of the Royal Statistical Society, B}, {\bf 58} 267-288.

\bibitem{Wu}
Wu, Y. N., Si, Z., Fleming, C. and Zhu, S. C. (2007) Deformable template as active basis. {\em Proceedings of International Conference of Computer Vision}.

\bibitem{Wu1}
Wu, Y. N., Si, Z., Gong, H. and Zhu, S. C. (2009) Learning active
basis model for object detection and recognition. {\em International
Journal of Computer Vision}, in press.


\bibitem{Yuille}
Yuille, A. L., Hallinan, P. W. and Cohen, D. S. (1992) Feature extraction from faces using deformable templates. {\em International Journal of Computer Vision}, {\bf 8} 99-111.

\bibitem{Zhu_grammar}
Zhu, S. C. and  Mumford, D. B. (2006) A stochastic grammar of images. {\em Foundations and Trends in Computer Graphics and Vision}, {\bf 2} 259-362.

\end{thebibliography}
}

\end{document}
