
%%%%%%%%%%%%%%%%%%%%%%% file template.tex %%%%%%%%%%%%%%%%%%%%%%%%%
%
% This is a general template file for the LaTeX package SVJour3
% for Springer journals.          Springer Heidelberg 2010/09/16
%
% Copy it to a new file with a new name and use it as the basis
% for your article. Delete % signs as needed.
%
% This template includes a few options for different layouts and
% content for various journals. Please consult a previous issue of
% your journal as needed.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% First comes an example EPS file -- just ignore it and
% proceed on the \documentclass line
% your LaTeX will extract the file if required
\begin{filecontents*}{example.eps}
%!PS-Adobe-3.0 EPSF-3.0
%%BoundingBox: 19 19 221 221
%%CreationDate: Mon Sep 29 1997
%%Creator: programmed by hand (JK)
%%EndComments
gsave
newpath
  20 20 moveto
  20 220 lineto
  220 220 lineto
  220 20 lineto
closepath
2 setlinewidth
gsave
  .4 setgray fill
grestore
stroke
grestore
\end{filecontents*}
%
\RequirePackage{fix-cm}
%
%\documentclass{svjour3}                     % onecolumn (standard format)
%\documentclass[smallcondensed]{svjour3}     % onecolumn (ditto)
%\documentclass[smallextended]{svjour3}       % onecolumn (second format)
\documentclass[twocolumn]{svjour3}          % twocolumn
%
\smartqed  % flush right qed marks, e.g. at end of proof
%
\usepackage{graphicx}
%
% \usepackage{mathptmx}      % use Times fonts if available on your TeX system
%
% insert here the call for the packages your document requires
%\usepackage{latexsym}
% etc.
%
% please place your own definitions here and don't use \def but
% \newcommand{}{}
%
% Insert the name of "your journal" with
% \journalname{myjournal}
%

\usepackage{times}
\usepackage{url}
\usepackage{epsfig}
\usepackage{graphicx}
\usepackage{subfig}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{mathrsfs}
\usepackage{algorithm}
\usepackage{algorithmicx}
\usepackage{array}
\usepackage{bm}
\usepackage{algpseudocode}
\usepackage[pagebackref=true,breaklinks=true,letterpaper=true,colorlinks,bookmarks=false]{hyperref}
\usepackage[misc]{ifsym}
\renewcommand{\algorithmicrequire}{\textbf{Input:}}
\renewcommand{\algorithmicensure}{\textbf{Output:}}

\def\E{{\rm E}}
\def\Var{{\rm Var}}
\def\Cov{{\rm Cov}}
\def\I{{\bf I}}
\def\J{{\bf J}}
\def\tI{\tilde{\bf I}}
\def\B{{\bf B}}
\def\bB{\bar{\bf B}}
\def\bb{\bar{B}}
\def\R{{\bf R}}


\begin{document}

\title{Learning Sparse FRAME Models for Natural Image Patterns}
%\subtitle{Do you have a subtitle?\\ If so, write it here}

%\titlerunning{Short form of title}        % if too long for running head

\author{Jianwen Xie        \and
        Wenze Hu \and
        Song-Chun Zhu \and
        Ying Nian Wu
}

\authorrunning{J. Xie et al.}



\institute{J. Xie $\cdot$ W. Hu $\cdot$ S.-C. Zhu $\cdot$ Y. N. Wu (\Letter) \\
Department of Statistics, UCLA\\
Los Angeles, CA, USA \\
             % Tel.: +123-45-678910\\
             % Fax: +123-45-678910\\
              \email{ywu@stat.ucla.edu}           %  \\
%             \emph{Present address:} of F. Author  %  if needed
         %  \and
         %  S. Author \at
         %     second address
}
%

\date{Received: 1 February 2014 / Accepted: 13 August 2014}
% The correct dates will be entered by the editor

\journalname{International Journal of Computer Vision (IJCV)}

\maketitle

\begin{abstract}
It is well known that natural images admit sparse representations by redundant dictionaries of basis functions such as Gabor-like wavelets. However, it is still an open question as to what the next layer of representational units above the layer of wavelets should be. We address this fundamental question by proposing a sparse FRAME (Filters, Random field, And Maximum Entropy) model for representing natural image patterns. Our sparse FRAME model is an inhomogeneous generalization of the original FRAME model. It is a non-stationary Markov random field model that reproduces the observed statistical properties of filter responses at a subset of selected locations, scales and orientations.  Each sparse FRAME model is intended to represent an object pattern and can be considered a deformable template. The sparse FRAME model can be written as a shared sparse coding model, which motivates us to propose a two-stage algorithm for learning the model. The first stage selects the subset of wavelets from the dictionary by a shared matching pursuit algorithm. The second stage then estimates the parameters of the model given the selected wavelets. Our experiments show that the sparse FRAME models are capable of representing a wide variety of object patterns in natural images and that the learned models are useful for object classification. 

\keywords{Generative models \and Markov random fields \and Shared sparse coding}
% \PACS{PACS code1 \and PACS code2 \and more}
% \subclass{MSC code1 \and MSC code2 \and more}
\end{abstract}

\section{Introduction}

\subsection{Background and motivation}

Sparsity underlies various  types of data arising from different scientific disciplines. For natural images, the seminal work of Olshausen and Field (1996) \cite{Olshausen_Field} showed that natural image patches admit sparse linear representations by an over-complete or redundant dictionary of basis functions that resemble Gabor wavelets. The dictionary of basis functions can be learned from training image patches by minimizing the reconstruction error with a sparsity-inducing penalty as in the original work of Olshausen and Field. It can also be learned by pursuit-based algorithms such as K-SVD \cite{KSVD}. The learned dictionaries prove to be useful for tasks such as image recovery \cite{sparsity} \cite{Denoising} and image classification \cite{4} \cite{Fergus}. 

During the past decade, a rich literature has been developed on learning dictionaries of wavelets for sparse coding. See, for example, the recent book by Elad \cite{Elad} and the references therein. However, it remains an open question as to what the next layer of representational units above the layer of wavelets should be. The goal of this article is to address this fundamental question by proposing a class of statistical models for representing natural image patterns based on wavelets sparse coding. In this class of models, each model is composed of a subset of wavelets selected from the dictionary of wavelets. The model assumes that the image intensities are generated by a linear superposition of the selected wavelets, and the model implies a probability distribution on the coefficients of the selected wavelets. 

\subsection{Model and algorithm}

One technical difficulty with modeling the coefficients of the selected wavelets explicitly is that it is difficult to specify the multi-dimensional joint distribution of the coefficients, and yet it is unrealistic to assume that the coefficients are statistically independent. To get around the difficulty, we choose instead to model  the responses of the image to the selected wavelets by adopting the mathematical form of the FRAME (Filters, Random field, And Maximum Entropy) model of Zhu, Wu, and Mumford (1997) \cite{Minimax} . 

The FRAME model, which was originally proposed for stochastic texture patterns, is a spatially stationary or homogeneous Markov random field model. Furthermore, it  is the maximum entropy distribution that reproduces the observed marginal histograms of responses from a bank of filters, where for each filter tuned to a specific scale and orientation, the marginal histogram is spatially pooled over all the pixels in the image domain. 

By modifying the original FRAME model, we propose a sparse FRAME model for representing natural image patterns. It is a generative model with a well-defined probability distribution on the image intensities. Unlike the original FRAME model for texture patterns, each sparse FRAME model is intended to model an object pattern, and can be considered a deformable template for this pattern. It is spatially non-stationary or inhomogeneous, and it is the maximum entropy distribution that reproduces statistical properties of filter responses at a  subset of selected locations, scales and orientations. 

The sparse FRAME model can be written as a shared sparse coding model, where the observed images are represented by a commonly shared subset of wavelets at selected locations, scales and orientations, subject to local perturbations to account for shape deformations. The sparse FRAME model implicitly assumes a joint probability distribution on the coefficients of the selected wavelets. 

We then propose a two-stage algorithm to learn the sparse FRAME model from roughly aligned image patches. The first stage selects a subset of wavelets to simultaneously reconstruct all the observed images, while allowing the selected wavelets to perturb their locations and orientations to represent each individual image. The second stage of the algorithm then estimates the parameters of the model given the selected wavelets. This stage  implicitly estimates the probability distribution of the coefficients of the  subset of selected wavelets. The computation of the second stage can be accomplished by stochastic gradient ascent \cite{Younes}, which seeks to reproduce the observed statistical properties of the responses of the selected wavelets.  Our experiments show that the learned model or template can synthesize realistic images and can be used to detect similar object patterns in the testing images. The learning algorithm can also be used for clustering images of different patterns. 

The above two-stage algorithm can learn a single sparse FRAME model from a training set of aligned images. It can also be employed to  learn a codebook of sparse FRAME models or templates from non-aligned training images, so that each image can be represented by a small number of spatially translated, rotated and scaled versions of the templates selected from the learned codebook. We use an unsupervised learning algorithm that iterates the following two steps: (1) Image encoding: Given the current codebook, select templates to encode the training images using a template matching pursuit algorithm. (2) Codebook re-learning: Given the encoding of the training images, re-learn the templates from the training image patches by the two-stage learning algorithm. Our experiments show that it is possible to learn codebooks of sparse FRAME models and the learned models are useful for image classification. 

In this paper, we assume that the bank of filters is given. For example, the filter bank contains Gabor filters and Difference of Gaussian (DoG) filters as in the original FRAME model. In other words, we assume that there already exists a dictionary of wavelets that gives sparse representations of the observed images. Presumably this dictionary can itself be learned if there are enough training data. We shall not pursue this issue in this paper, and shall focus on learning the generative models based on the given dictionary of wavelets. 

\subsection{Related work}

The two stages of the learning algorithm for training the sparse FRAME model naturally connect two major frameworks in image representation and modeling, namely the sparse coding framework with its root in harmonic analysis and the Markov random field framework with its root in statistical physics. There have been vast literatures on both themes of research. In the following, we shall review and compare with some of the papers that are most relevant to our work. 

{\em Markov random field models.} Models in this class are also called energy-based models \cite{Energy} \cite{BM}, exponential family models, and Gibbs distributions depending on the context. Examples include the FRAME model \cite{Minimax} as well as its inhomogeneous extension for face shape data \cite{LiuCe}, field of experts \cite{FoE}, product of experts \cite{PoE}, product of t model \cite{PoT}, restricted Boltzmann machine \cite{RBM} \cite{Hinton} and its many recent generalizations such as those found in \cite{Ranzato} and the references therein. A Markov random field model  is defined by an energy function and may involve latent variables or hidden units. If the latent variables are conditionally independent given the observed data or visible units,  the latent variables can be integrated out in closed form, resulting in a marginal energy-based model. These models usually assume fixed energy functions and do not involve explicit feature selection. In addition, these models seek to approximate the probability distributions of the training images but do not attempt to reconstruct individual training images. Compared to these models, the sparse FRAME model performs feature selection via a linear additive model and can reconstruct the training images. 

{\em Sparse coding models.} The sparse FRAME model selects the wavelets via a shared sparse coding scheme. Such a scheme has been studied in harmonic analysis and signal processing \cite{Chen_Huo} \cite{Tropp} under the name of simultaneous sparse coding, and in statistics and machine learning \cite{Supportunion} \cite{Multitask} under the names of multi-task learning and support union recovery. These methods seek to reconstruct the observed signals but do not attempt to approximate the probability distributions of these signals. In contrast, the sparse FRAME model defines an explicit probability distribution on image intensities and can synthesize new images by sampling from this distribution. 

There has also been work on learning dictionaries that are sparse combinations of wavelets from a base dictionary, such as \cite{doubleSparsity}, where the coefficients of sparse linear combinations are fixed. In our model, the coefficients are allowed to vary according to a certain probability distribution. In addition, the wavelets are also allowed to perturb their locations and orientations. Our work is also related to subspace clustering \cite{subspace}, where each cluster is spanned by a subset of wavelets or basis functions. In subspace clustering, the distributions of the coefficients of the basis functions that span the subspaces are not modeled. In our work, we seek to model the coefficients or responses of the selected basis functions.

The proposed model (and the original FRAME model) is related to ``analysis priors'' developed in the last few years by the sparse modeling community \cite{AnalysisModel} \cite{AnalysisSyn}. Our learning algorithm can be viewed as a principled statistical method to learn the parameters or weights of the analysis prior models via maximum likelihood. Here our focus is on statistical modeling and random sampling, whereas the focus of the analysis priors is on optimization and reconstruction tasks. 

{\em Deep learning}. Our work is not exactly within the domain of deep learning, but is closely related to it. In particular, we try to understand the layer of representational units or nodes above the layer of sparse coding wavelets. Our proposal is that each node  is a sparse FRAME model, which is selectively and sparsely connected to a subset of wavelets selected for this model. The connection weights are the parameters of the model. The sparse connections are obtained by seeking the shared sparse coding of a collection of image patches of similar patterns. The advantage of seeking the shared sparse coding is that the time-consuming explaining-away computation in sparse coding can be memorized in the learning stage by sparse connections, so that explaining-away sparsification does not need to be recalculated on-line in the inference stage. In other words, we believe that ``sparse connectivities = shared sparse activities'' or ``sparse wiring = shared sparse coding''.  In contrast, the current methods of deep learning are mostly based on stacking restricted Boltzmann machines (RBM) \cite{Hinton}  \cite{LeeNg}  or auto-encoders \cite{Bengio}.  They do not pursue explicit sparse representations and sparse connections.  Our work is also related to the deconvolution network of \cite{Fergus}. It appears that in the deconvolution model, given the values of the top layer units, the values of the units at lower layers are fixed by the learned weights, and are not allowed to vary. In our model, we allow the coefficients to vary according to the learned probability distribution. 

Our method can be extended to learning hierarchical models. After learning one layer of sparse FRAME models, we can treat these models as re-usable parts, and continue to compose them into higher layers of sparse FRAME models. 

{\em Compositional models}. The sparse FRAME model is a special case of compositional models advocated by S. Geman et al. for vision \cite{Geman_composition}. In particular, it follows the And-Or grammar studied by Zhu and Mumford  \cite{Zhu_grammar}, where the composition of the wavelets forms an And-node, and the perturbation of each selected wavelet and the variation of its coefficient form an Or-node. Our model is also related to \cite{LeoZhu} \cite{Fidler} , which are about compositions of edgelets but which are not based on explicit generative models as in our work.  

Compared to our own previous work,  this paper can be considered a fusion of the original FRAME model  \cite{Minimax}  and the active basis model \cite{Wu_active_basis_IJCV}  \cite{ABC}. While the active basis model focuses on the ``sketching'' aspect, this paper adds the ``painting'' aspect. In order to avoid MCMC computation in learning, the active basis model makes the simplifying assumptions that the selected wavelets are orthogonal and their coefficients are statistically independent. In this paper, we do not make such simplifying assumptions, and thus our model is more rigorously defined and is capable of synthesizing realistic image patterns. This paper is an expanded version of our conference paper \cite{iFRAME}. 

\subsection{Contributions}

The following are the main contributions of this paper. (1) We propose an inhomogeneous dense FRAME model for object patterns, and we show that it can model a wide variety of objects in natural scenes. (2) We propose a sparse FRAME model and connect it to the shared sparse coding model. We then propose a two-stage algorithm for learning the sparse FRAME model. (3) We show that it is possible to learn codebooks of sparse FRAME models from non-aligned and unannotated images. 


\section{Inhomogeneous FRAME model} \label{sect:inhomogeneous}

\begin{figure}[h]
\begin{center}
\includegraphics[height=0.15\linewidth]{hummingbird/hummingbird/gray/h1.jpg}
\includegraphics[height=0.15\linewidth]{hummingbird/hummingbird/gray/h6.jpg}
\includegraphics[height=0.15\linewidth]{hummingbird/hummingbird/gray/h8.jpg}
\includegraphics[height=0.15\linewidth]{hummingbird/hummingbird/gray/h9.jpg}
\includegraphics[height=0.15\linewidth]{hummingbird/hummingbird/gray/h10.jpg}\\ \vspace{.5mm}
\includegraphics[height=0.15\linewidth]{cat/cat/gray/1.jpg}
\includegraphics[height=0.15\linewidth]{cat/cat/gray/2.jpg}
\includegraphics[height=0.15\linewidth]{cat/cat/gray/3.jpg}
\includegraphics[height=0.15\linewidth]{cat/cat/gray/4.jpg}
\includegraphics[height=0.15\linewidth]{cat/cat/gray/5.jpg}
\includegraphics[height=0.15\linewidth]{cat/cat/gray/8.jpg}
\end{center}
\caption{The inhomogeneous FRAME  is a generative model that seeks to represent and generate object patterns as shown above. }
\label{fig:patterns}
\end{figure}


This section presents a dense version of the inhomogeneous FRAME model to lay the foundation for the next section which will focus on the sparsified version. 

\subsection{Model and learning algorithm}

{\em Notation}. We start by modeling roughly aligned images of object patterns from the same category, such as the images in Figure 1. Let $\{\I_m, m = 1, ..., M\}$ be a set of training images defined on an image domain ${\cal D}$. We use the notation $B_{x, s, \alpha}$ to denote a basis function such as a Gabor wavelet centered at pixel $x$ (which is a two-dimensional vector), and tuned to scale $s$ and orientation $\alpha$. $B_{x, s, \alpha}$ is also an image on ${\cal D}$,  although it is non-zero only within a local range. We assume that $B_{x, s, \alpha}$ are translated, dilated and rotated versions of each other. We assume that $s$ and $\alpha$ take values within a finite and properly discretized range. The inner product $\langle \I, B_{x, s, \alpha}\rangle$ can be considered the filter response of $\I$ to a filter of scale $s$ and orientation $\alpha$ at pixel $x$. We assume that $B_{x, s, \alpha}$ are all normalized to have unit $\ell_2$ norm. 

{\em Model}. The inhomogeneous FRAME model is a probability distribution defined on $\I$, 
\begin{eqnarray} 
   p(\I; \lambda) = \frac{1}{Z(\lambda)} \exp \left( \sum_{x, s, \alpha} \lambda_{x, s, \alpha}(\langle \I, B_{x, s, \alpha}\rangle) \right) q(\I), \label{eq:denseFRAME}
\end{eqnarray} 
where $q(\I)$ is a known reference distribution or a null model, $\lambda_{x, s, \alpha}()$ are one-dimensional functions that depend on $(x, s, \alpha)$, $\lambda = \{\lambda_{x, s, \alpha}, \forall x, s, \alpha\}$, and 
\begin{eqnarray}
Z(\lambda) &=& \int  \exp \left( \sum_{x, s, \alpha} \lambda_{x, s, \alpha}(\langle \I, B_{x, s, \alpha}\rangle)\right) q(\I) d\I \\
&=& \E_{q}\left[\exp\left( \sum_{x, s, \alpha} \lambda_{x, s, \alpha}(\langle \I, B_{x, s, \alpha}\rangle)\right) \right]
\end{eqnarray}
 is the normalizing constant, where the notation $\E_{q}$ means the expectation with respect to the probability distribution $q$. 

In the original FRAME model for stochastic textures \cite{Minimax}, $\lambda_{x, s, \alpha}()$ is assumed to be independent of $x$ (but dependent of $s$ and $\alpha$, which index the scale and orientation of the filter). So the model is spatially stationary. For modeling object patterns that are not spatially stationary, $\lambda_{x, s, \alpha}()$ must depend on $x$, in addition to $s$ and $\alpha$. 

In the original homogeneous FRAME, the potential functions $\lambda_{s, \alpha}()$ (we drop the subscript $x$ due to stationarity) are estimated non-parametrically as step functions. In the inhomogeneous FRAME, we have to estimate $\lambda_{x, s, \alpha}()$ for each individual $x$. With small data sets, we may not afford estimating $\lambda_{x, s, \alpha}()$ non-parametrically. We therefore decide to parametrize 
\begin{eqnarray}
\lambda_{x, s, \alpha}(r) = \lambda_{x, s, \alpha} |r|, \label{eq:parametrization}
\end{eqnarray}
where $r = \langle \I, B_{x, s, \alpha}\rangle$, and with slight abuse of notation, $\lambda_{x, s, \alpha}$ on the right hand side of the above equation becomes a constant (instead of a function as on the left hand side). The parametrization (\ref{eq:parametrization}) is inspired by the Laplacian distribution that can account for heavy tails in the distributions of filter responses. It is possible to replace the function $|r|$ by other classes of parametrized functions to encourage heavy tails of the responses, and we shall investigate this issue in future work. 

In many Markov random field models including the original FRAME model, the reference measure $q(\I)$ is simply the uniform measure. In our work,  we assume $q(\I)$ to be the Gaussian white noise model, under which the image intensities follow independent ${\rm N}(0, \sigma^2)$ distributions. So 
\begin{eqnarray}
q(\I) = \frac{1}{(2\pi\sigma^2)^{|{\cal D}|/2}} \exp\left( - \frac{1}{2\sigma^2} \sum_{x} \I(x)^2\right),
\end{eqnarray}
where $|{\cal D}|$ is the number of pixels in the image domain. $q(\I)$ itself is a maximum entropy model relative to a uniform measure, and it reproduces the marginal mean and variance of the image intensities. In our work, we normalize the observed images to have marginal mean 0 and a fixed variance, and we set $\sigma^2 = 1$. This $q(\I)$ can be considered an initial model or a model of the background residual image with the foreground object removed. As a result, $p(\I; \lambda)$ in equation (\ref{eq:denseFRAME}) can be written as an exponential family model relative to the uniform measure. 

{\em Maximum likelihood learning}. The inhomogeneous version of the  FRAME model is a special case of the exponential family model, and the parameter $\lambda = (\lambda_{x, s, \alpha}, \forall x, s, \alpha)$ can be estimated from the training images $\{\I_m, m = 1, ..., M\}$ by maximum likelihood. The log-likelihood function is
\begin{eqnarray}
L(\lambda)&=&\frac{1}{M} \sum_{m=1}^{M}\log p(\I_m;\lambda)\\
&=&\frac{1}{M} \sum_{m=1}^{M}\sum_{x,s,\alpha} \lambda_{x,s,\alpha}|\langle \I_m, B_{x, s, \alpha}\rangle|  -\log Z(\lambda) \\ 
&&+ \frac{1}{M}\sum_{m=1}^{M} \log q(\I_m). \nonumber 
\end{eqnarray} 
The maximization of $L(\lambda)$ can be accomplished by gradient ascent. The gradient is 
\begin{eqnarray}
\frac{\partial L(\lambda)}{\partial \lambda_{x,s,\alpha}}=\frac{1}{M}\sum_{m=1}^{M} |\langle \I_m, B_{x, s, \alpha}\rangle|&& - \E_{p(\I;\lambda)}\left[|\langle \I, B_{x, s, \alpha}\rangle|\right], \nonumber \\&&\forall x,s,\alpha,
\end{eqnarray} 
where $\E_{p(\I;\lambda)}[|\langle \I, B_{x, s, \alpha}\rangle|]$ is the expectation of $|\langle \I, B_{x, s, \alpha}\rangle|$ with $\I$ following the distribution $p(\I;\lambda)$. $\E_{p(\I;\lambda)}[|\langle \I, B_{x, s, \alpha}\rangle|]$ is the derivative of $\log Z(\lambda)$. 

The gradient ascent algorithm then becomes  
\begin{eqnarray}
\lambda^{(t+1)}_{x, s, \alpha} = \lambda^{(t)}_{x, s, \alpha} + 
\gamma_t \Bigg(\frac{1}{M} \sum_{m=1}^{M} |\langle \I_m, B_{x, s, \alpha}\rangle| - \nonumber \\ 
 \E_{p(\I;\lambda^{(t)})}\left[|\langle \I, B_{x, s, \alpha}\rangle|\right]\Bigg), \label{eq:learning0}
\end{eqnarray}
where $\gamma_t$ is the step size. The analytic form of the expectation under the current model at step $t$, $\E_{p(\I;\lambda^{(t)})}[|\langle \I, B_{x, s, \alpha}\rangle|]$, is not available, so we approximate it from a sample set of synthesized images $\{\tI_{m}, m = 1, ..., \tilde{M}\}$ generated from $p(\I; \lambda^{(t)})$:
\begin{eqnarray}
\E_{p(\I;\lambda)}[|\langle \I, B_{x, s, \alpha}\rangle|] \approx  \frac{1}{\tilde{M}} \sum_{m=1}^{\tilde{M}} |\langle \tI_m, B_{x, s, \alpha}\rangle|. \label{eq:w1}
\end{eqnarray}
The synthesized images $\{\tI_m\}$ can be sampled from $p(\I; \lambda^{(t)})$ by Hamiltonian Monte Carlo (HMC) \cite{HMC}. Unlike the Gibbs sampler \cite{Gibbs}, HMC makes use of the gradient of the energy function, and it is particularly natural for our model. The computation of HMC involves a bottom-up convolution step followed by a top-down convolution step. Both steps can be efficiently implemented in Matlab by GPU. With HMC and warm start, $\{\tI_m\}$ are produced by $\tilde{M}$ parallel chains. All the synthesized images presented in the figures of this paper are generated by the HMC algorithm along the learning process. More details about simulation by the HMC algorithm are presented in Section \ref{sect:HMC} in the appendix. 

With $\E_{p(\I;\lambda)}[|\langle \I, B_{x, s, \alpha}\rangle|] $ approximated according to (\ref{eq:w1}), we arrive at the stochastic gradient algorithm analyzed by Younes (1999)
\cite{Younes} : 
\begin{eqnarray}
\lambda^{(t+1)}_{x, s, \alpha} = \lambda^{(t)}_{x, s, \alpha} + 
\gamma_t \Bigg(\frac{1}{M} \sum_{m=1}^{M} |\langle \I_m, B_{x, s, \alpha}\rangle| - \nonumber \\ 
\frac{1}{\tilde{M}} \sum_{m=1}^{\tilde{M}} |\langle \tI_m, B_{x, s, \alpha}\rangle|\Bigg).  \label{eq:w2}
\end{eqnarray}
This is the algorithm we use for maximum likelihood estimation of $\lambda$. 

{\em Computing normalizing constants.} Thanks to HMC, we can simulate from $p(\I;\lambda)$ without knowing its normalizing constant, thus estimating $\lambda$ by MLE. Nevertheless, computing the normalizing constant $Z(\lambda)$ is still required in situations such as fitting a mixture model or learning a codebook of models. The ratio of the normalizing constants at two consecutive steps is 
\begin{eqnarray}
\frac{Z(\lambda^{(t+1)})}{Z(\lambda^{(t)})} = \E_{p(\I;\lambda^{(t)})}\Bigg[\exp \Big(\sum_{x, s, \alpha} (\lambda^{(t+1)}_{x, s, \alpha} - \lambda^{(t)}_{x, s, \alpha}) \nonumber\\
\times |\langle \I, B_{x, s, \alpha}\rangle|\Big) \Bigg] \label{eq:Z}
\end{eqnarray} 
which can be approximated by averaging over the sampled images $\{\tI_{m}\}$ as an application of importance sampling \cite{bridgeSampling}:
\begin{eqnarray}
\frac{Z(\lambda^{(t+1)})}{Z(\lambda^{(t)})} \approx \frac{1}{\tilde{M}}\sum_{m=1}^{\tilde{M}}\Bigg[\exp \Big(\sum_{x, s, \alpha} (\lambda^{(t+1)}_{x, s, \alpha} - \lambda^{(t)}_{x, s, \alpha}) \nonumber\\
\times |\langle \tI_m, B_{x, s, \alpha}\rangle|\Big) \Bigg]. \label{eq:Z2}
\end{eqnarray} 
Starting from $\lambda^{(0)} = 0$ and $\log Z(\lambda^{(0)}) = 0$, we can compute $\log Z(\lambda^{(t)})$ along the learning process by iteratively updating its value as follows:
\begin{eqnarray}
\log Z(\lambda^{(t+1)}) = \log Z(\lambda^{(t)}) + \log \frac{Z(\lambda^{(t+1)})}{Z(\lambda^{(t)})}.
\end{eqnarray}

The calculation of $Z$ is based on running parallel Markov chains for a sequence of distributions $p(\I; \lambda^{(t)})$. The setting is similar to annealed importance sampling \cite{ImportanceSampling} and bridge sampling \cite{bridgeSampling}. We shall explore these methods in future work. 

\subsection{Summary of the learning algorithm} \label{sect:iFRAMEalgorithm}

Pseudocode of the algorithm for learning the inhomogeneous FRAME model is shown in Algorithm \ref{code:FRAME}. The algorithm stops when the gradient of the log-likelihood is close to 0, i.e., when the statistics of the synthesized images closely match those of the observed images. Figure \ref{fig:synthesized} displays the synthesized images $\{\tI_{m}\}$ generated by the models learned from training images shown in Figure \ref{fig:patterns}  (a separate model is learned from each training set). Figure \ref{fig:learningprocess} illustrates the learning process by showing the synthesized images with $\lambda$ being updated by the algorithm. The synthesized image starts from Gaussian white noises sampled from $q(\I)$,  then gradually gets similar to the observed images in the overall shape and appearance. 

The computational complexity of  Algorithm \ref{code:FRAME} is of the order $O(U \times \tilde{M} \times L \times K \times H_B \times W_B)$ with $U$ being the number of updating steps for $\lambda$, $\tilde{M}$ the number of synthesized images, $L$ the number of leapfrog steps in HMC, $K$ the number of filters, and $H_B$ and $W_B$ the average window sizes (height and width) of the filters. As to the actual running time, for the cat example, each iteration of a single chain takes about 2 seconds on a current PC, with $L=30, K=240100, H_B=12$, and $W_B=12$. 

\begin{figure}[h]
\begin{center}
\includegraphics[height=0.15\linewidth]{denseFRAME/h1.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/h2.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/h3.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/h4.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/h5.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/h6.png}\\ \vspace{.5mm}
\includegraphics[height=0.15\linewidth]{denseFRAME/c1.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/c2.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/c3.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/c4.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/c5.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/cat.png}
\end{center}
\caption{Synthesized images generated by the inhomogeneous FRAME models learned separately from a training set of 5 hummingbird images and another training set of 12 cat images. Some of the images are displayed in Fig. \ref{fig:patterns}. The sizes of the images are $70 \times 70$.  }
\label{fig:synthesized}
\end{figure}

\begin{figure*}
\begin{center}
\includegraphics[width=1\linewidth]{learningprocess/cat_bird.jpg}
\end{center}
\caption{Learning sequence by inhomogeneous FRAME. The sizes of the images are $70 \times 70$. A separate model is learned from each training set shown in Fig 1 (Top: Hummingbird. The number of training images is 5 as shown in  Fig. \ref{fig:patterns}. Bottom: Cat. The number of training images is 12, with 6 of them shown in Fig. \ref{fig:patterns}.) Synthesized images are generated in iterations $t = 1, 4, 7, 10, 13, 20, 50, 100, 200, 300, 400$, and $500$.}
\label{fig:learningprocess}
\end{figure*}

\begin{algorithm}
\caption{Learning algorithm for dense FRAME}
\label{code:FRAME}
\begin{algorithmic}[1]
\Require 
\Statex training images $\{\I_m, m=1,...,M\}$ 
\Ensure
\Statex $\lambda=\{\lambda_{x,s,\alpha}, \forall x,s,\alpha\}$ and $\log Z(\lambda)$ 
\Statex
\State Create a filter bank $\{B_{x,s,\alpha}, \forall x,s,\alpha \}$
\State Initialize $\lambda^{(0)}_{x,s,\alpha} \leftarrow 0, \forall x,s,\alpha$.
\State Calculate observed statistics:

$H^{obs}_{x,s,\alpha} \leftarrow \frac{1}{M}\sum_{m=1}^{M} |\langle \I_m, B_{x, s, \alpha}\rangle|, \forall x,s,\alpha.$ 
\State Initialize synthesized images $\tI_m$ as Gaussian white noise images
\State Initialize $\log Z(\lambda^{(0)})\leftarrow 0$ 
\State Let $t\leftarrow 0$
\Repeat 
\State Generate $\{\tI_m,m=1,...,\tilde{M}\}$ from $p(\I;\lambda^{(t)})$ by HMC  

\State Calculate synthesized statistics:

$   H^{syn}_{x,s,\alpha} \leftarrow \frac{1}{\tilde{M}}\sum_{m=1}^{\tilde{M}} |\langle \tI_m, B_{x, s, \alpha}\rangle|, \forall x,s,\alpha.$
\State Update $\lambda^{(t+1)}_{x,s,\alpha} \leftarrow \lambda^{(t)}_{x,s,\alpha}+ \gamma_t ( H^{obs}_{x,s,\alpha} - H^{syn}_{x,s,\alpha}) $, $\forall x,s,\alpha$.

\State Compute $Z$ ratio $\frac{Z(\lambda^{(t+1)})}{Z(\lambda^{(t)})}$ by Eq. (\ref{eq:Z2})
\State Update $\log Z(\lambda^{(t+1)}) \leftarrow \log Z(\lambda^{(t)}) + \log \frac{Z(\lambda^{(t+1)})}{Z(\lambda^{(t)})}$ 
\State Let $t \leftarrow t+1$
\Until $\sum_{x,s,\alpha} |H^{obs}_{x,s,\alpha} - H^{syn}_{x,s,\alpha}| \leq \epsilon$
\end{algorithmic}
\end{algorithm}

\begin{figure}[h]
\begin{center}
\includegraphics[height=0.15\linewidth]{denseFRAME/cat.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/lion.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/lionese.png}
%\includegraphics[height=0.15\linewidth]{denseFRAME/cougar.png}
%\includegraphics[height=0.15\linewidth]{denseFRAME/snowleopard.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/tiger.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/wolf.png}
%\includegraphics[height=0.15\linewidth]{denseFRAME/bear.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/panda.png}\\ \vspace{.5mm}
\includegraphics[height=0.15\linewidth]{denseFRAME/horsehead.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/deerface.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/deer.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/cow.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/pony.png}\\ \vspace{.5mm}
\includegraphics[height=0.14\linewidth]{denseFRAME/dog.png}
\includegraphics[height=0.14\linewidth]{denseFRAME/dog1.png}
%\includegraphics[height=0.15\linewidth]{denseFRAME/camel.png}\\ \vspace{.5mm}
\includegraphics[height=0.15\linewidth]{denseFRAME/elephant.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/zebra.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/leopard.png}\\ \vspace{.5mm}
\includegraphics[height=0.15\linewidth]{denseFRAME/pigeon.png} 
\includegraphics[height=0.15\linewidth]{denseFRAME/seagull.png}
\includegraphics[height=0.12\linewidth]{denseFRAME/swan.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/goose.png}
\includegraphics[height=0.18\linewidth]{denseFRAME/egret.png}
%\includegraphics[height=0.18\linewidth]{denseFRAME/crane.png}
\includegraphics[height=0.18\linewidth]{denseFRAME/flamingo.png}\vspace{.5mm}
\includegraphics[height=0.24\linewidth]{denseFRAME/giraffe.png}
\includegraphics[height=0.21\linewidth]{denseFRAME/woman.png} 
\includegraphics[height=0.21\linewidth]{denseFRAME/women3.png}
\includegraphics[height=0.10\linewidth]{denseFRAME/butterfly.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/sunflower.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/susan.png}
\\ \vspace{.5mm}
\includegraphics[height=0.15\linewidth]{denseFRAME/star.png}
%\includegraphics[height=0.15\linewidth]{denseFRAME/five.png}
%\includegraphics[height=0.15\linewidth]{denseFRAME/fivepetal.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/forgetmenot.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/springbeauty.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/fivepetal.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/sixpetal.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/eight.png}
\\ \vspace{.5mm}
\includegraphics[height=0.15\linewidth]{denseFRAME/peach.png} 
\includegraphics[height=0.15\linewidth]{denseFRAME/tomato.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/teapot.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/cup.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/watch.png}
\includegraphics[height=0.15\linewidth]{denseFRAME/bag.png}  \\ \vspace{.5mm}
\includegraphics[height=0.145\linewidth]{denseFRAME/motorcycle.png}
\includegraphics[height=0.145\linewidth]{denseFRAME/bike.png}
\includegraphics[height=0.145\linewidth]{denseFRAME/car.png}
\includegraphics[height=0.145\linewidth]{denseFRAME/carfront.png}
\end{center} 
   \caption{Synthesis by dense FRAME. Images generated by the dense FRAME models learned from different categories of objects. The training images are collected from the internet and are cropped so that the training images for each category are roughly aligned. Typical number of training images for each category is around 10. }
\label{fig:denseFRAME}
\end{figure}

{\em Experiment 1: Learning dense FRAME}. Figure \ref{fig:denseFRAME} displays some images generated by the dense models learned from roughly aligned training images. We run a single chain in the learning process, i.e., $\tilde{M} = 1$ in this experiment. The learned models can generate a wide variety of natural image patterns. Typical sizes of the images are $70 \times 70$. 


In the appendix, Section \ref{sect:app1} gives a justification of the inhomogeneous FRAME model by the maximum entropy  principle. 

\section{Sparse FRAME model}  \label{sect:sparse}

{\em Sparsification}. In model (\ref{eq:denseFRAME}), the $(x, s, \alpha)$ in $\sum_{x, s, \alpha}$ is over all the pixels $x$ and all the scales $s$ and orientations $\alpha$. We call such a model the dense FRAME. It is possible to sparsify the model by selecting only a small set of $(x, s, \alpha)$ so that $\sum_{x, s, \alpha}$ is restricted to this selected subset. More explicitly, we can write the sparsified model as 
\begin{eqnarray} 
   p(\I; \B, \lambda) = \frac{1}{Z(\lambda)} \exp\left( \sum_{i=1}^{n} \lambda_i|\langle \I, B_{x_i, s_i, \alpha_i}\rangle|\right) q(\I), \label{eq:sparseFRAME}
\end{eqnarray} 
where $\B = (B_{x_i, s_i, \alpha_i}, i = 1, ..., n)$ are the selected basis functions, and  $\lambda = (\lambda_i, i = 1, ..., n)$ collects the parameters. Given the selected basis functions $\B$, the model can still be trained by maximum likelihood as in the previous section, and properties such as maximum entropy still hold. The following are the reasons why a sparsified model is desirable. (1) It makes the computation faster. (2) It leads to more reliable parameter estimates because it involves a much smaller number of parameters. Estimation efficiency or  accuracy is an important aspect of statistical modeling. (3) The MCMC sampling may converge faster if the selected basis functions are not heavily correlated. (4) It is connected to the linear additive sparse coding model for image reconstruction. (5) It allows the selected basis functions to perturb their locations and orientations to account for shape deformations. 

{\em Deformation}. To be more specific about the above point (5), we may treat $p(\I; \B, \lambda)$ as a deformable template, so that when it is fitted to each training image $\I_m$,  we may allow the basis functions in $\B = (B_{x_i, s_i, \alpha_i}, i = 1, ..., n)$ to perturb their locations and orientations so that $\B$ is deformed to $\B_m = (B_{x_i +  \Delta x_{m, i}, s_i, \alpha_i + \Delta  \alpha_{m, i}}, i = 1, ..., n)$, where $(\Delta x_{m, i}, \Delta \alpha_{m, i})$ are the perturbations of the location and orientation of the $i$-th basis function $B_{x_i, s_i, \alpha_i}$. Both $\Delta x_{m, i}$ and  $\Delta \alpha_{m, i}$ are assumed to vary within limited ranges (default setting:  $\Delta x_{m, i} \in [-3, 3]$ pixels along the normal direction of the Gabor wavelet, and $\Delta \alpha_{m, i} \in \{-1, 0, 1\} \times \pi/16$). When we fit the model $p(\I; \B, \lambda)$ to $\I_m$, we model $\I_m$ by $p(\I_m; \B_m, \lambda)$, in which $B_{x_i, s_i, \alpha_i}$ in (\ref{eq:sparseFRAME}) is changed to $B_{x_i +  \Delta x_{m, i}, s_i, \alpha_i + \Delta  \alpha_{m, i}}$. 


\subsection{Shared sparse coding}

We can select $B_{x_i, s_i, \alpha_i}$ or $(x_i, s_i, \alpha_i)$ in model (\ref{eq:sparseFRAME}) sequentially using a procedure like projection pursuit \cite{Friedman_projection_pursuit} or filter pursuit \cite{Minimax}, but the computational speed of such a sequential procedure can be slow. In this article, we choose to employ a different strategy by exploring the connection between sparse FRAME and shared sparse coding. 

For simplicity, let us temporarily ignore the issue of deformation. For the sparse model in equation (\ref{eq:sparseFRAME}), the number of selected basis functions $n$ is always much smaller than the number of pixels $|{\cal D}|$. We can then project each $\I_m \sim p(\I; \lambda)$ onto the subspace spanned by the selected basis functions $\B = (B_{x_i, s_i, \alpha_i}, i = 1, ..., n)$, so that 
\begin{eqnarray}
   \I _m= \sum_{i=1}^{n} c_{m, i} B_{x_i, s_i, \alpha_i} + \epsilon_m,   \label{eq:nodeform}
\end{eqnarray}
where $c_{m, i}$ are the least squares reconstruction coefficients of the linear projection, and $\epsilon_m$ is the resulting residual image that resides in the $|{\cal D}| - n$ dimensional residual subspace that is orthogonal to the subspace spanned by $\B = (B_{x_i, s_i, \alpha_i}, i = 1, ..., n)$. Equation (\ref{eq:nodeform}) is a shared sparse coding model, where the small set of basis functions $\B = (B_{x_i, s_i, \alpha_i}, i = 1, ..., n)$ is shared by the training images $\{\I_m, m = 1, ..., M\}$. 

With the Gaussian white noise background model $q(\I)$ where each $\I(x) \sim {\rm N}(0, \sigma^2)$ independently, the sparse model (\ref{eq:sparseFRAME}) implies that $C_m = (c_{m,i}, i = 1, ..., n)$ follows a certain distribution $p_C(C; \lambda)$,  $\epsilon_m$ is the projection of the Gaussian white noise image onto the $|{\cal D}|-n$ residual subspace, and $C_m$ and $\epsilon_m$ are independent of each other. The log-likelihood of $\I_m$ can be decomposed into the log-likelihood of $C_m$ and the log-likelihood of the projected Gaussian white noise $\epsilon_m$. While the former depends on $\lambda$, the latter only depends on the squared norm of the residual image  $\|\epsilon_m\|^2 = \|\I_m - \sum_{i=1}^{n} c_{m, i} B_{x_i, s_i, \alpha_i}\|^2$. 

The above consideration suggests a two-stage learning algorithm for fitting the sparse FRAME model. In the first stage, we selected $\B = (B_{x_i, s_i, \alpha_i}, i = 1, ..., n)$ by minimizing the overall least squares reconstruction error $\sum_{m=1}^{M} \|\epsilon_m\|^2$. In the second stage, we then estimate $\lambda$ given the selected $\B$.

In the appendix, Section \ref{sect:app2} gives a more detailed explanation of the connections between the sparse FRAME model and the shared sparse coding model. In particular, it shows that the sparse FRAME model is equivalent to the shared sparse coding model with an implied joint distribution on the coefficients of the selected basis functions. 

Now let us consider the issue of deformation. Since model (\ref{eq:sparseFRAME}) is deformable, we can also make the sparse coding model  (\ref{eq:nodeform}) deformable by allowing the shared basis functions to perturb their locations and orientations to account for the shape deformation in each image. This leads to the deformable shared sparse coding first proposed in our previous work on active basis  \cite{Wu_active_basis_IJCV}
\begin{eqnarray}
  \I_m = \sum_{i=1}^{n} c_{m, i} B_{x_i +  \Delta x_{m, i}, s_i, \alpha_i + \Delta  \alpha_{m, i}} + \epsilon_m, \label{eq:deform}
\end{eqnarray} 
where $(\Delta x_{m, i}, \Delta \alpha_{m, i})$ are the perturbations of the location and orientation of the $i$-th basis function. 
  
\subsection{The two-stage learning algorithm} \label{sect:learning}

\begin{algorithm}[h]
\caption{Stage 1:  Deformable shared matching pursuit} \label{alg:MP}
\begin{algorithmic}[1]
\Require 
\Statex training images $\{\I_m, m=1,...,M\}$ 
\Ensure
\Statex selected basis functions $\B=\{B_{x_i,s_i,\alpha_i}, i=1,...,n\}$
\Statex
\State
Initialize $i \leftarrow 0$. For $m = 1, ..., M$, initialize the residual image $\epsilon_m \leftarrow \I_m$.

\State Let $i \leftarrow i + 1$. Then we select 
\[
(x_i, s_i, \alpha_i) = \arg\max_{x, s, \alpha} \sum_{m=1}^{M} \max_{\Delta x, \Delta \alpha}
|\langle \epsilon_m, B_{x+\Delta x, s, \alpha+\Delta \alpha}\rangle|^2,
\]
where $\max_{\Delta x, \Delta \alpha}$ is  local max pooling 
within the small ranges of $\Delta x_{m, i}$ and  $\Delta \alpha_{m,
i}$.

\State For each $m$, given $(x_i, s_i, \alpha_i)$, infer the
perturbations in locations and orientations by retrieving the arg-max
in the local max pooling of step 2:
\[
(\Delta x_{m, i}, \Delta \alpha_{m, i}) = \arg\max_{\Delta x, \Delta \alpha} |\langle \epsilon_m,
B_{x_i+\Delta x, s_i, \alpha_i+\Delta \alpha}\rangle|^2. \label{eq:argMax}
\]
Let the coefficient
\[
c_{m, i} \leftarrow \langle \epsilon_m, B_{x_i+\Delta x_{m, i}, s_i, \alpha_i+\Delta \alpha_{m, i}}\rangle,
\]
and update the residual image by explaining away:
\[
\epsilon_m \leftarrow \epsilon_m - c_{m, i} B_{x_i+\Delta x_{m, i}, s_i, \alpha_i+\Delta \alpha_{m, i}}.
\label{eq:subtraction}
\]

\State Stop if $i = n$, else go back to step 2.
\end{algorithmic}
\end{algorithm}



\begin{figure}[h]
\begin{center}
\includegraphics[height=0.17\linewidth]{hummingbird/output/LearnedTemplate1.png}
\includegraphics[height=0.17\linewidth]{hummingbird/output/LearnedTemplate2.png}
\includegraphics[height=0.17\linewidth]{hummingbird/output/LearnedTemplate3.png}
\includegraphics[height=0.17\linewidth]{hummingbird/output/LearnedTemplate4.png}
\includegraphics[height=0.17\linewidth]{hummingbird/output/LearnedTemplate6.png}\\ \vspace{.5mm}
\includegraphics[height=0.18\linewidth]{hummingbird/output/ReconstructedImage1.png}\\ \vspace{.5mm}
\includegraphics[height=0.18\linewidth]{hummingbird/output/ReconstructedImage2.png}\\ \vspace{.5mm}
\includegraphics[height=0.18\linewidth]{hummingbird/output/ReconstructedImage3.png}\\ \vspace{.5mm}
%\includegraphics[height=0.18\linewidth]{hummingbird/output/ReconstructedImage4.png}\\ \vspace{.5mm}
\includegraphics[height=0.18\linewidth]{hummingbird/output/ReconstructedImage5.png}\\ \vspace{.5mm}
\includegraphics[height=0.22\linewidth]{hummingbird/EachChainReconSyn/reconstruct_chain-0001.png}
\includegraphics[height=0.22\linewidth]{hummingbird/EachChainReconSyn/reconstruct_chain-0004.png}
\includegraphics[height=0.22\linewidth]{hummingbird/EachChainReconSyn/reconstruct_chain-0031.png}
\includegraphics[height=0.22\linewidth]{hummingbird/EachChainReconSyn/reconstruct_chain-0033.png}
\end{center}
\caption{Reconstruction and synthesis by sparse FRAME model (hummingbirds).  The number of selected wavelets is 300. The first row contains symbolic sketches of selected Gabor wavelets at different scales, where each Gabor wavelet is illustrated by a bar. The first 4 sketches correspond to 4 different scales. The last one is the superposition of the 4  scales, where smaller scales appear darker.  The next 4 rows display examples of the training images,  the deformed sketches, the reconstructed images, and the residual images. The last row displays examples of synthesized images generated by the learned model. The number of training images is 5 as shown in  Fig. \ref{fig:patterns}. The sizes of images are scaled to 100 $\times$ 100.}
\label{fig:sparseFRAME}
\end{figure}



\begin{figure}[h]
\begin{center}
\includegraphics[height=0.17\linewidth]{cat/output/LearnedTemplate1.png}
\includegraphics[height=0.17\linewidth]{cat/output/LearnedTemplate2.png}
\includegraphics[height=0.17\linewidth]{cat/output/LearnedTemplate3.png}
\includegraphics[height=0.17\linewidth]{cat/output/LearnedTemplate4.png}
\includegraphics[height=0.17\linewidth]{cat/output/LearnedTemplate6.png}\\ \vspace{.5mm}
\includegraphics[height=0.18\linewidth]{cat/output/ReconstructedImage12.png}\\ \vspace{.5mm}
\includegraphics[height=0.18\linewidth]{cat/output/ReconstructedImage6.png}\\ \vspace{.5mm}
\includegraphics[height=0.18\linewidth]{cat/output/ReconstructedImage7.png}\\ \vspace{.5mm}
\includegraphics[height=0.18\linewidth]{cat/output/ReconstructedImage9.png}\\ \vspace{.5mm}
\includegraphics[height=0.22\linewidth]{cat/EachChainReconSyn/reconstruct_chain-0003.png}
\includegraphics[height=0.22\linewidth]{cat/EachChainReconSyn/reconstruct_chain-0005.png}
\includegraphics[height=0.22\linewidth]{cat/EachChainReconSyn/reconstruct_chain-0020.png}
\includegraphics[height=0.22\linewidth]{cat/EachChainReconSyn/reconstruct_chain-0034.png}
\end{center}
   \caption{Reconstruction and synthesis  (cats). See caption of Fig. \ref{fig:sparseFRAME}. The number of training images is 12, with 6 of them shown in Fig.   \ref{fig:patterns}. The sizes of images are scaled to 100 $\times$ 100. The number of selected wavelets is 300. }
\label{fig:sparseFRAME2}
\end{figure}


This subsection describes the learning algorithm for training the sparse FRAME model, which consists of two stages. (1) Selecting $\B  = (B_{x_i, s_i, \alpha_i}, i = 1, ..., n)$ by shared sparse coding. (2) Estimating $\lambda = (\lambda_i, i = 1, ..., n)$ given the selected $\B$.  

{\em Stage 1: Deformable shared sparse coding}. For training images $\{\I_m, m = 1, ..., M\}$, we select the basis functions $\{B_{x_i, s_i, \alpha_i}, i = 1, ..., n\}$ by minimizing 
\begin{eqnarray}
\sum_{m=1}^{M} \| \I_m - \sum_{i=1}^{n} c_{m, i} B_{x_i +  \Delta x_{m, i}, s_i, \alpha_i + \Delta  \alpha_{m, i}}\|^2.
\end{eqnarray}
The minimization can be accomplished by the shared matching pursuit algorithm \cite{Mallat_Zhang} \cite{Wu_active_basis_IJCV}  that selects basis functions to encode multiple images simultaneously, while inferring local perturbations by local max pooling \cite{Riesenhuber_Poggio}. The algorithm is presented in Algorithm \ref{alg:MP}.  

We can replace the matching pursuit component in the above algorithm by the orthogonal matching pursuit \cite{OMP}, which is more computationally expensive. We can also replace the matching pursuit component by penalized least squares such as basis pursuit \cite{basis pursuit} or Lasso \cite{lasso}, or more precisely using a penalty such as the $\ell_1/\ell_2$ group norm \cite{Supportunion}. The computation can be much more expensive than shared matching pursuit. 

Simultaneous sparse approximation of multiple signals  has been studied in the harmonic analysis  and machine learning literature \cite{Tropp} \cite{Supportunion}. However, perturbations of the selected basis functions are not considered in these papers. Such a deformable shared matching pursuit algorithm was first proposed by \cite{Wu_active_basis_IJCV}, but it implemented a modified version that enforces approximated orthogonality of the selected basis functions. 


 {\em Stage 2: Sparse FRAME as deformable template}. After selecting $\B = \{B_{x_i, s_i, \alpha_i}, i = 1, ..., n\}$, we can then model $\{\I_m\}$ by the sparse FRAME model (\ref{eq:sparseFRAME}), by estimating $\lambda$ via MLE. $p(\I; \B, \lambda)$ in (\ref{eq:sparseFRAME}) now serves as the deformable template in that the log-likelihood  of  $\I_m$ is 
\begin{eqnarray}
L(\I_m| \B, \lambda) = &&\sum_{i=1}^{n} \lambda_i \max_{\Delta x, \Delta \alpha} |\langle \I_m, B_{x_i+\Delta x, s_i, \alpha_i + \Delta \alpha}\rangle| \nonumber 
\\ &&- \log Z(\lambda), \label{eq:matching}
\end{eqnarray}
which serves as the template matching score. We allow each selected $B_{x_i, s_i, \alpha_i}$ to perturb its location and orientation to account for shape deformation, where the perturbation is inferred by the local max pooling in Algorithm \ref{alg:MP}.

In the learning algorithm, again, let $\lambda^{(t)}$ be the current estimate of $\lambda$, and let $\{\tI_{m}, m = 1, ..., \tilde{M}\}$ be the synthesized images drawn from $p(\I; \B, \lambda^{(t)})$ by $\tilde{M}$ parallel chains. Then we update $\lambda$ by 
\begin{equation}
\begin{aligned}
\lambda^{(t+1)}_{i} = \lambda^{(t)}_{i} +
\gamma_t &\Bigg(\frac{1}{M} \sum_{m=1}^{M} \max_{\Delta x, \Delta \alpha} |\langle \I_m, B_{x_i+\Delta x, s_i, \alpha_i + \Delta \alpha}\rangle|  \label{eq:learning} \\ 
               &- \frac{1}{\tilde{M}} \sum_{m=1}^{\tilde{M}} |\langle \tI_m, B_{x_i, s_i, \alpha_i}\rangle|\Bigg).
\end{aligned}
\end{equation}
The learned $p(\I;  \B, \lambda)$ models the appearance of the undeformed template. There is no local max pooling on the synthesized images, which have not undergone shape deformations or warping. The local max pooling is only applied to the observed images to filter out the shape deformations in the observed images. Thus there is an explicit separation between appearance and shape variations. 

Again the synthesized images can be sampled by the HMC algorithm. For HMC computation, the energy function is 
\begin{eqnarray}
U(\I)=-\sum_{i=1}^n \lambda_{i} | \langle \I, B_{x_i, s_i, \alpha_i} \rangle |+ \frac{1}{2} |\I|^2,
\end{eqnarray}
 and the gradient of this energy function is 
\begin{eqnarray}
\frac{\partial U}{\partial \I}=-\sum_{i=1}^{n} \lambda_i {\rm sign}(\langle \I, B_{x_i, s_i, \alpha_i}\rangle) B_{x_i, s_i, \alpha_i}+\I,
\end{eqnarray}
so HMC  is like a generative process based on linear superpositions of $\B=(B_{x_i, s_i, \alpha_i}, i = 1, ..., n)$. With the separation between appearance and shape, the learned model for appearance may not be very multi-modal, therefore the HMC sampling can be quite fast. 

\begin{algorithm}
\caption{Stage 2: Parameter estimation in sparse model}
\label{code:sparseFRAME}
\begin{algorithmic}[1]
\Require 
\Statex (1) training images $\{\I_m, m=1,...,M\}$, 
\Statex (2) selected basis functions $\B=\{B_{x_i,s_i,\alpha_i},i=1,...,n\}$ 
\Statex (3) inferred perturbations $\{\Delta x_{m,i}, \Delta \alpha_{m,i}, m=1,...,M, i=1,...,n\}$ by local max pooling. 

\Ensure
\Statex $\lambda=\{\lambda_i, i=1,...,n\}$ and $\log Z(\lambda)$ 
\Statex
\State Initialize $\lambda^{(0)}_{i} \leftarrow 0, i=1,...,n$.
\State Calculate observed statistics:

$H^{obs}_{i} \leftarrow \frac{1}{M}\sum_{m=1}^{M} |\langle \I_m, B_{x_i +\Delta x_{m,i},s_i,\alpha_i+ \Delta \alpha_{m,i} } \rangle|$, for $i=1,...,n.$ 
\State Initialize synthesized images $\tI_m$ as Gaussian white noise images
\State Initialize $\log Z(\lambda^{(0)})\leftarrow 0$ 
\State Let $t\leftarrow 0$
\Repeat 
\State Generate $\{\tI_m,m=1,...,\tilde{M}\}$ from $p(\I; \B, \lambda^{(t)})$ by HMC  

\State Calculate synthesized statistics:

$   H^{syn}_{i} \leftarrow \frac{1}{\tilde{M}}\sum_{m=1}^{\tilde{M}} |\langle \tI_m, B_{x_i,s_i,\alpha_i}\rangle|$, for $i=1,...,n.$
\State Update $\lambda^{(t+1)}_{i} \leftarrow \lambda^{(t)}_{i}+ \gamma_t ( H^{obs}_{i} - H^{syn}_{i})$, $i=1,...,n$.

\State Compute $Z$ ratio $\frac{Z(\lambda^{(t+1)})}{Z(\lambda^{(t)})}$ by Eq. (\ref{eq:Z2})
\State Update $\log Z(\lambda^{(t+1)}) \leftarrow \log Z(\lambda^{(t)}) + \log \frac{Z(\lambda^{(t+1)})}{Z(\lambda^{(t)})}$ 
\State Let $t \leftarrow t+1$
\Until $\sum_{i} |H^{obs}_{i} - H^{syn}_{i}| \leq \epsilon$
\end{algorithmic}
\end{algorithm}


The algorithm is presented in Algorithm \ref{code:sparseFRAME}. After we learn $\lambda$ and compute $Z(\lambda)$ as in (\ref{eq:Z2}), we can use the learned model as a deformable template to be matched to the testing image, where the template matching score can be computed according to (\ref{eq:matching}).

Figure \ref{fig:sparseFRAME} illustrates the basic idea of training the sparse FRAME model. The training images are scaled to $100 \times 100$. The number of selected basis functions (Gabor and large DoG wavelets), $n$, is set at 300. In principle it can be automatically determined by criteria like BIC.  In the first stage, by using the  deformable shared  matching pursuit algorithm (see Algorithm \ref{alg:MP}) on the training images, we select $n$ wavelets  $\B = (B_{x_i, s_i, \alpha_i}, i = 1, ..., n)$, which are displayed in the first row, where each  $B_{x_i, s_i, \alpha_i}$  is symbolized by a bar. The first four plots in the first row display the selected  $B_{x_i, s_i, \alpha_i}$  at 4 different scales $s_i$, from the largest to the smallest. The last plot in the first row is a superposition of the 4 scales, with smaller scales appearing darker. The next four rows of the figure display four training images $\I_m$, the symbolic sketches of the deformed templates, $\B_m = (B_{x_i +  \Delta x_{m, i}, s_i, \alpha_i + \Delta  \alpha_{m, i}}, i = 1, ..., n)$, the reconstructed images obtained by the linear superpositions of the perturbed basis functions, $\sum_{i=1}^{n} c_{m, i} B_{x_i +  \Delta x_{m, i}, s_i, \alpha_i + \Delta  \alpha_{m, i}}$, and the residual image $\epsilon_m$. At the second stage, we fit the sparse FRAME model with the $n$ selected wavelets (see Algorithm  \ref{code:sparseFRAME}). The synthesized images $\tilde{\I}_m$ generated from the learned model $p(\I; \B, \lambda)$ are  projected onto the subspace spanned by $\B$. The last row displays projections of four synthesized images. These synthesized images show the appearances before shape deformations. Figure \ref{fig:sparseFRAME2} shows another example. 

\begin{figure}[h]
\begin{center}
\includegraphics[height=0.19\linewidth]{sparseFRAME/cat.png}
\includegraphics[height=0.19\linewidth]{sparseFRAME/lion.png}
\includegraphics[height=0.19\linewidth]{sparseFRAME/cougar.png}
\includegraphics[height=0.19\linewidth]{sparseFRAME/tiger.png}
%\includegraphics[height=0.19\linewidth]{sparseFRAME/snowleopard.png}
\includegraphics[height=0.19\linewidth]{sparseFRAME/wolf.png} \\ \vspace{.5mm}
\includegraphics[height=0.19\linewidth]{sparseFRAME/deerface.png}
\includegraphics[height=0.19\linewidth]{sparseFRAME/horsehead.png}
\includegraphics[height=0.19\linewidth]{sparseFRAME/pony.png}
\includegraphics[height=0.19\linewidth]{sparseFRAME/dog.png}\\ \vspace{.5mm}
\includegraphics[height=0.19\linewidth]{sparseFRAME/deer.png}
\includegraphics[height=0.19\linewidth]{sparseFRAME/camel.png}
\includegraphics[height=0.19\linewidth]{sparseFRAME/elephant.png}
\includegraphics[height=0.19\linewidth]{sparseFRAME/cow.png}\\ \vspace{.5mm}
\includegraphics[height=0.18\linewidth]{sparseFRAME/zebra.png}
\includegraphics[height=0.18\linewidth]{sparseFRAME/leopard.png}
%\includegraphics[height=0.18\linewidth]{sparseFRAME/tigerbody.png}
\includegraphics[height=0.26\linewidth]{sparseFRAME/giraffe.png}
\includegraphics[height=0.26\linewidth]{sparseFRAME/woman.png}\\ \vspace{.5mm}
%\includegraphics[height=0.26\linewidth]{sparseFRAME/women3.png}
\includegraphics[height=0.22\linewidth]{sparseFRAME/egret.png} 
\includegraphics[height=0.22\linewidth]{sparseFRAME/crane.png}
\includegraphics[height=0.22\linewidth]{sparseFRAME/flamingo.png}
\includegraphics[height=0.19\linewidth]{sparseFRAME/seagull.png}
\includegraphics[height=0.19\linewidth]{sparseFRAME/goose.png}\\ \vspace{.5mm}
\includegraphics[height=0.19\linewidth]{sparseFRAME/pigeon.png}
\includegraphics[height=0.19\linewidth]{sparseFRAME/duck.png}
\includegraphics[height=0.19\linewidth]{sparseFRAME/swan.png}
\includegraphics[height=0.19\linewidth]{sparseFRAME/forgetmenot.png}\\ \vspace{.5mm}
\includegraphics[height=0.19\linewidth]{sparseFRAME/springbeauty.png}
\includegraphics[height=0.19\linewidth]{sparseFRAME/eight.png}
\includegraphics[height=0.19\linewidth]{sparseFRAME/bike.png}
\includegraphics[height=0.14\linewidth]{sparseFRAME/car.png}
%\includegraphics[height=0.19\linewidth]{sparseFRAME/teapot.png}
\end{center}
   \caption{Synthesis by sparse FRAME. Images generated by the sparse FRAME models learned from different categories of objects. Typical sizes of the images are $80 \times 80$. Typical number of selected wavelets is 300. The training images are the same as or similar to those for training the dense FRAME models in Fig. \ref{fig:denseFRAME}.  }
\label{fig:sparseFRAME3}
\end{figure}

{\em Experiment 2: Synthesis by Sparse FRAME}. Figure \ref{fig:sparseFRAME3} displays some images generated by the sparse models learned from roughly aligned images. The experiment setting is the same as that in Figure \ref{fig:sparseFRAME} except that the image sizes are typically 80 $\times$ 80, and the allowed displacement of Gabor wavelet is up to 2 pixels. The number of wavelets is 300. We run $\tilde{M} = 36$ parallel chains in the learning algorithm. Even though the number of wavelets are greatly reduced compared to the dense model, the sparse model can still generate realistic object patterns, including  highly textured patterns. Because of the relatively small number of parameters, it is unlikely that the model memorizes the training images. 

As to the actual running time, for the cat example, with 12 training images, the shared matching pursuit in stage 1 takes 95 seconds. For the  algorithm of learning $\lambda$ in stage 2, each iteration takes 2.8 second. The total running time is 6.5 minutes. 


\begin{figure}
\begin{center}
\textbf{(a)}
\includegraphics[height=0.20\linewidth]{comparison/0203.png}
\includegraphics[height=0.20\linewidth]{comparison/cat_dense.png}
\includegraphics[height=0.20\linewidth]{comparison/zebra_dense.png}
\includegraphics[height=0.20\linewidth]{comparison/lion_dense.png}\\ \vspace{.5mm}
\textbf{(b)}
%\includegraphics[height=0.20\linewidth]{comparison/bird_OMP.png}
%\includegraphics[height=0.20\linewidth]{comparison/cat2.png}
%\includegraphics[height=0.20\linewidth]{comparison/zebra_OMP.png}
%\includegraphics[height=0.20\linewidth]{comparison/lion_OMP.png}\\ \vspace{.5mm}
%(c)
%\includegraphics[height=0.20\linewidth]{comparison/bird_twostages.png}
\includegraphics[height=0.20\linewidth]{hummingbird/EachChainReconSyn/reconstruct_chain-0004.png}
\includegraphics[height=0.20\linewidth]{cat/EachChainReconSyn/reconstruct_chain-0003.png}
%\includegraphics[height=0.20\linewidth]{comparison/cat_twostages.png}
\includegraphics[height=0.20\linewidth]{comparison/zebra_twostages.png}
\includegraphics[height=0.20\linewidth]{comparison/lion_twostage.png}
\end{center}
\caption{Comparison of synthesized images generated by 
(a) dense FRAME and (b) sparse FRAME, where the number of selected wavelets is 300. Image sizes are about 100 $\times$ 100. 
}
\label{fig:synComparison}
\end{figure}

For a comparison of different models and learning methods, Figure \ref{fig:synComparison} displays synthesized images generated by the dense FRAME and the sparse FRAME respectively. 

\section{Detection} \label{sect:det}

After learning the sparse FRAME model $p(\I; \B, \lambda)$, where $\B = (B_{x_i, s_i, \alpha_i}, i = 1, ..., n)$ and $\lambda = (\lambda_i, i = 1, ..., n)$, from the roughly aligned training images, we can use the learned model to detect the object in a testing image by deformable template matching. 

Let $\I$ be a testing image defined on the domain ${\cal D}$. We can scan the template over ${\cal D}$, and at each location $X \in \cal D$, we match the template to the image patch of $\I$ within the bounding box centered at $X$ by computing the log-likelihood or the template matching score  based on (\ref{eq:matching}), 
\begin{eqnarray}
L(\I \mid \B_{X}, \lambda) 
= &&\sum_{i=1}^{n}  \lambda_i
\max_{\Delta x, \Delta \alpha} |\langle \I, B_{X+x_i+\Delta x, s_i, \alpha_i + \Delta \alpha}\rangle| \nonumber \\
&&- \log Z(\lambda), \label{eq:A0}
\end{eqnarray}
where we use $\B_X = (B_{X+x_i+\Delta x, s_i, \alpha_i + \Delta \alpha}, i = 1, ..., n)$ to denote the spatially translated and deformed version of the template $\B$. The perturbations of the basis functions are inferred by local max pooling as above. We then choose the location $X$ that achieves  the maximum template matching score as the center of the detected object. In practice, the template can be partially outside the image domain ${\cal D}$  when we scan the template near the boundary of ${\cal D}$. In this case, we only need to set the filter responses outside ${\cal D}$ to be zero. Also, to deal with the scaling issue, we can apply the above algorithm at multiple resolutions of the testing image, and then choose the resolution that achieves the maximum template matching score as the optimal resolution.


\begin{figure}
\begin{center}
\includegraphics[width=1.\linewidth]{transformation3.jpg}
\end{center}
 \caption{Geometric transformation. Flipping: the first row shows an example of left/right flipping transformation, where the first two images are the synthesized image and the symbolic sketch template of the learned model, while the next two images correspond to the flipped model derived from the learned model. Rotation: the other example illustrated in the next two rows displays the rotated models at four different orientations (-90, -45, 45, and 90 degrees) by showing their synthesized images and symbolic sketches. The middle column is from the original learned model.}
\label{fig:transformation}
\end{figure}


\begin{figure*}
\begin{center}
\includegraphics[width=.8\linewidth]{cat_detection5.jpg}
\end{center}
 \caption{Detection. (a) Symbolic sketch template representing 250 selected wavelets.  (b) A synthesized image generated by the learned model. We do not include the large DoG filters in the model, so the synthesized image lacks regional contrast. (c) Testing images with bounding boxes locating the detected objects.}
\label{fig:detection}
\end{figure*}

In addition to spatial translation in scanning, we can also allow geometric transformations such as rotation and left-right flipping of the template. The geometrically transformed versions of the learned model can be obtained by directly applying the operations of dilation, rotation, flipping, or even changing the aspect ratio to $\B=\{B_{x_i, s_i,\alpha_i},i=1,...,n\}$ without changing the values of $\lambda$. This amounts to simple affine transformations of $(x_i, s_i,\alpha_i,i=1,...,n)$. Figure \ref{fig:transformation} shows two examples of geometric transformations of the sparse FRAME model: flipping and rotation. It displays synthesized images generated by the transformed models, which are derived from the learned model, as well as the corresponding symbolic sketches representing the selected wavelets. For better performance in detection, we can first generate a collection of models at different orientations and aspect ratios from the learned model. After that, we use these transformed models to detect the object. We choose the combination of the transformed template and image resolution that gives the best match in terms of the template matching score, and infer the hidden location, orientation, and scale of the detected object in the testing image. 

{\em Experiment 3:  Detection by sparse FRAME.} Figure \ref{fig:detection} shows examples of detection. We learn the model from eight roughly aligned training images, with $\tilde{M} = 36$. The template size is $100 \times 100$. The two images displayed in Figure \ref{fig:detection}(a) and \ref{fig:detection}(b) are symbolic sketches  showing 250 wavelets selected by deformable shared matching pursuit algorithm and a synthesized image generated by the learned model. We transform the learned model into a collection of models at 9 different orientations, and then run the detection algorithm over 17 resolutions of the testing images using these transformed templates. Figure \ref{fig:detection}(c) displays the detection results by drawing bounding boxes on the detected objects. 

This detection algorithm can be combined with the two-stage learning algorithm to learn from training images that are not well aligned by alternating the following two steps. (1) Re-learning the model from the currently aligned training images by the two-stage algorithm. (2) Re-aligning the training images by the detection algorithm. 

\section{Clustering} \label{sect:cluster}

\begin{figure}[t]
\begin{center}
\includegraphics[height=0.17\linewidth, width=0.21\linewidth]{dogclustering/2.png}
\includegraphics[height=0.17\linewidth, width=0.21\linewidth]{dogclustering/12_gray.jpg}
\includegraphics[height=0.17\linewidth, width=0.21\linewidth]{dogclustering/1.png}
\includegraphics[height=0.17\linewidth, width=0.21\linewidth]{dogclustering/w_gray.jpg}\\ \vspace{.5mm}
\includegraphics[height=0.17\linewidth, width=0.21\linewidth]{horseclustering/2.png}
\includegraphics[height=0.17\linewidth, width=0.21\linewidth]{horseclustering/z2_gray.jpg}
\includegraphics[height=0.17\linewidth, width=0.21\linewidth]{horseclustering/1.png}
\includegraphics[height=0.17\linewidth, width=0.21\linewidth]{horseclustering/R_gray.jpg}\\ \vspace{.5mm}
\includegraphics[height=0.17\linewidth, width=0.21\linewidth]{cowclustering/1.png}
%\includegraphics[height=0.1\linewidth]{cowclustering/LearnedTemplate4-1.png}
\includegraphics[height=0.17\linewidth, width=0.21\linewidth]{cowclustering/C11_gray.jpg}
\includegraphics[height=0.17\linewidth, width=0.21\linewidth]{cowclustering/2.png}
%\includegraphics[height=0.1\linewidth]{cowclustering/LearnedTemplate4.png}
\includegraphics[height=0.17\linewidth, width=0.21\linewidth]{cowclustering/C5_gray.jpg}\\ \vspace{.5mm}
\includegraphics[height=0.17\linewidth, width=0.23\linewidth]{bigcat/leopardS.png}
\includegraphics[height=0.17\linewidth, width=0.23\linewidth]{bigcat/leopard.jpg}
\includegraphics[height=0.17\linewidth, width=0.23\linewidth]{bigcat/tigerS.png}
\includegraphics[height=0.17\linewidth, width=0.23\linewidth]{bigcat/tiger.jpg}\\ \vspace{.5mm}
\includegraphics[height=0.11\linewidth, width=0.156\linewidth]{insectclustering/2.png}
%\includegraphics[height=0.15\linewidth]{insectclustering/LearnedTemplate4.png}
\includegraphics[height=0.11\linewidth, width=0.156\linewidth]{insectclustering/B14_gray.jpg}
\includegraphics[height=0.11\linewidth, width=0.156\linewidth]{insectclustering/3.png}
%\includegraphics[height=0.15\linewidth]{insectclustering/LearnedTemplate4-1.png}
\includegraphics[height=0.11\linewidth, width=0.156\linewidth]{insectclustering/10_gray.jpg}
\includegraphics[height=0.11\linewidth, width=0.156\linewidth]{insectclustering/1.png}
%\includegraphics[height=0.15\linewidth]{insectclustering/LearnedTemplate4-2.png}
\includegraphics[height=0.11\linewidth, width=0.156\linewidth]{insectclustering/b_gray.jpg}
\end{center}
   \caption{Clustering. Each row illustrates one clustering experiment by displaying a synthesized image and a training example for each cluster. The number of images within each cluster is around 15 to 20. Typical template sizes are $100 \times 100$. Typical number of wavelets for each template is 300. 
 }
\label{fig:clustering}
\end{figure}

\begin{figure*}
\begin{center}
\includegraphics[width=1.0\linewidth]{clusteringDataset.jpg}
\end{center}
 \caption{The clustering dataset. One example image is shown for each of the 22 clusters distributed in 7 clustering tasks. }
\label{fig:clusteringDataset}
\end{figure*}


\begin {table*}
\begin{center}
{\small
\caption {Comparison of conditional purity (the first two rows) and conditional entropy (the last two rows) between sparse FRAME and  
k-means for clustering.}
\label{table:clustering}
\begin{tabular}{|c|c|c|c|c|c|c|c|}
\hline
  & Exp 1 & Exp 2 & Exp 3 & Exp 4 & Exp 5 & Exp 6 & Exp 7\tabularnewline
\hline
\hline
k-means (purity) & 0.623$\pm$0.016 & 0.870$\pm$0.043 & 0.933$\pm$0.141 &
0.825$\pm$0.121 & 0.911$\pm$0.086 & 0.888$\pm$0.091 & 0.687$\pm$0.110  
\tabularnewline
\hline
FRAME (purity) & 0.943$\pm$0.063 & 0.990$\pm$0.016 & 0.938$\pm$0.131 &
0.895$\pm$0.132 & 1.000$\pm$0.000 & 0.879$\pm$0.141 & 0.741$\pm$0.111  
\tabularnewline
\hline
\hline
k-means (entropy) & 0.652$\pm$0.009 & 0.376$\pm$0.086 & 0.092$\pm$0.195 &
0.243$\pm$0.167 & 0.226$\pm$0.084 & 0.199$\pm$0.126 & 0.639$\pm$0.161  
\tabularnewline
\hline
FRAME (entropy)  & 0.145$\pm$0.157 & 0.037$\pm$0.060 & 0.090$\pm$0.191 &
0.155$\pm$0.189 & 0.000$\pm$0.000 & 0.179$\pm$0.208 & 0.497$\pm$0.192  
\tabularnewline
\hline
\end{tabular}}
\end{center}
\end{table*}


Model-based clustering can be accomplished by the EM-type algorithm \cite{EM} that fits mixtures of sparse FRAME models. Suppose we have $M$ images from $K$ clusters. For each image $\I_m$, we define $(z_m^{(k)},k=1,...,K)$ as a hidden indicator vector, where $z_m^{(k)}=1$ if $\I_m$ comes from cluster $k$, otherwise $z_m^{(k)}=0$. The EM-like clustering algorithm is a greedy scheme that infers $\{z_m^{(k)} \}$ and $\{\B^{(k)}, \lambda^{(k)}, k = 1, ..., K\}$ by maximizing the overall log-likelihood 
\begin{eqnarray}
\sum_{m=1}^M \sum_{k=1}^K  z_m^{(k)} L(\I_m| \B^{(k)}, \lambda^{(k)}),
\end{eqnarray}
where $\B^{(k)}$ are the basis functions selected for cluster $k$, $\lambda^{(k)}$ are the corresponding parameters, and $L(\I_m| \B^{(k)}, \lambda^{(k)})$ is the log-likelihood or template matching score defined by (\ref{eq:matching}). 

The algorithm is initialized by randomly generating $\{z_m^{(k)}\}$, and then iterates the following two steps:

{\em Re-learning:} Given $\{(z_m^{(k)},k=1,...,K), m=1,...,M \}$, learn the sparse FRAME model $p(\I;\B^{(k)},\lambda ^{(k)})$ from images classified into the $k$-th cluster: $\{\I_m, z_m^{k}=1 \}$, for each $k=1,...,K$. 


{\em Classification: } Given the learned models of the $K$ clusters: $\{p(\I;\B^{(k)},\lambda ^{(k)}), k=1,...,K \}$, assign each image $\I_m$ to a cluster $k_*$ that maximizes the template matching score $L(\I_m|\B^{(k)}, \lambda^{(k)})$ over all $k=1,...,K$. Set $z_m^{(k_*)}=1$, and set  $z_m^{(k)}=0$ for $k \neq k_*$.


In the above algorithm, the classification step corresponds to the E-step of the EM algorithm, except that we adopt hard classification instead of computing the expectation of $z_m$ for each image $\I_m$. The re-learning step corresponds to the M-step of the EM algorithm. The algorithm usually converges within a few iterations. 

{\em Experiment 4:  Model-based clustering.} Figure  \ref{fig:clustering} illustrates 5 experiments. The EM-type algorithm usually converges within 3-5 iterations, at which point all the images are correctly separated into their respective clusters. For each cluster, we generate $\tilde{M} = 144$ parallel chains in learning because we need to compute $Z(\lambda)$ accurately for each model, as multiple models compete to explain the images. The same  $\tilde{M} = 144$ is used for the experiments in the remaining part of the paper. 

{\em Experiment 5:  Numerical evaluation on clustering}. To evaluate the clustering accuracies, we use two measures: conditional purity and conditional entropy \cite{purity}.  For a random training image, let $x$ be its true category label and $y$ be the inferred category label.  The conditional purity is defined as $\sum_{y } p(y) \max_{x} p(x|y)$ (the larger the better), and the conditional entropy is  $\sum_{y} p(y) \sum_x p(x|y) \log (1/p(x|y))$ (the smaller the better), where both $p(y)$ and $p(x|y)$ can be estimated from the training data. We also introduce a new dataset for clustering. Figure \ref{fig:clusteringDataset} provides an overview of the dataset. It contains 7 clustering tasks. The numbers of clusters vary from 2 to 5 and are assumed known in these tasks. The number of images in each cluster is typically 15 except in one experiment. We compare the performance of the sparse FRAME with that of the k-means method based on HoG \cite{HOG} features by performing experiments on these 7 clustering tasks. Table \ref{table:clustering} displays the clustering accuracies and standard errors based on 10 repetitions of each experiment. 

\section{Unsupervised learning  from non-aligned images} \label{sect:our}

\begin{figure}[h]
\begin{center}
\includegraphics[height=0.2\linewidth]{codebooks/flyingT_new.png}\\ \vspace{.5mm}
\includegraphics[height=0.3\linewidth]{codebooks/flying_gray.jpg}
\includegraphics[height=0.3\linewidth]{codebooks/flyingR_new.png}\\
\textbf{(a) Seagull flying}
\\ \vspace{2.5mm}
\includegraphics[height=0.2\linewidth]{codebooks/bT1.png}
\includegraphics[height=0.2\linewidth]{codebooks/bT2.png}\\ \vspace{.5mm}
\includegraphics[height=0.4\linewidth]{codebooks/b2_gray.png}
\includegraphics[height=0.4\linewidth]{codebooks/bR.png}\\ 
\textbf{(b) Brick walls and floor tiles}
\vspace{.5mm}
%\includegraphics[height=0.4\linewidth]{codebooks/b3.png}
%\includegraphics[height=0.4\linewidth]{codebooks/b3R.png}
\end{center}
\caption{Unsupervised learning. In each experiment, a codebook of sparse FRAME templates (of size $100 \times 100$ pixels) is learned from the training image. The images on the top row are generated from the learned templates. The image on the left of the second row is the observed training image. The image on the right of the second row is reconstructed using spatially translated, rotated and deformed versions of the learned templates.  }
\label{fig:unsupervised}
\end{figure}

In the previous section, we consider learning a single sparse FRAME model or template from roughly aligned images. The two-stage learning algorithm can serve as a basis for learning a codebook of sparse FRAME templates from non-aligned images without any annotation and labeling, so that the training images can be represented by spatially translated, rotated, scaled and deformed versions of templates selected from the learned codebook. Here we follow the learning scheme in our previous work on compositional sparse coding \cite{ABC}.

\subsection{Learning a codebook of sparse FRAME models}

Figure \ref{fig:unsupervised} shows two experiments. In the first experiment, a single template is learned. In the second experiment, a codebook of two templates (brick and floor tile patterns)  are learned. In each experiment, the images on the top row are generated from the learned models. The image on the left of the second row is the observed image, and the image on the right of the second row is reconstructed by the learned templates.  

{\em Single template with spatial translation}. To fix the notation, we shall assume temporarily that the templates are only allowed spatial translations in encoding the training images. We start from generalizing the representation (\ref{eq:deform}) by assuming that the template may appear at location $X_m$ in image $\I_m$, then we can write the representation as
\begin{eqnarray}
\I_m &=& \sum_{i=1}^{n} c_{m, i} B_{X_m + x_i + \Delta x_{m, i}, s_i, \alpha_i + \Delta \alpha_{m, i}} + \epsilon_m  \\
&=& C_m \B_{X_m}+\epsilon_m, \label{eq:pack1}
\end{eqnarray}
where $\B_{X_m} = (B_{X_m + x_i + \Delta x_{m, i}, s_i, \alpha_i + \Delta \alpha_{m, i}}, i = 1, ..., n)$ is the deformed template spatially translated to $X_m$, $C_m = (c_{m, i}, i = 1, ..., n)$, and by definition 
\begin{eqnarray}
C_m \B_{X_m} = \sum_{i=1}^{n} c_{m, i} B_{X_m + x_i + \Delta x_{m, i}, s_i, \alpha_i + \Delta \alpha_{m, i}}.
\end{eqnarray}
$\B_{X_m}$ explains the part of $\I_m$ that is covered by $\B_{X_m}$.  For each image $\I_m$ and each $X_m$, the log-likelihood  is
\begin{eqnarray}
L(\I_m \mid \B_{X_m}) 
= &&\sum_{i=1}^{n}  \lambda_i
\max_{\Delta x, \Delta \alpha} |\langle \I_m, B_{X_m+x_i+\Delta x, s_i, \alpha_i + \Delta \alpha}\rangle| \nonumber \\
&&- \log Z(\lambda), \label{eq:A}
\end{eqnarray}
which is a slight generalization of (\ref{eq:matching}) and which is the log-likelihood score (\ref{eq:A0}) used for object detection. For notational simplicity, we drop $\lambda$ in $L(\I_m \mid \B_{X_m})$. We always assume that $\lambda$ is estimated by MLE. 


\begin{figure}[h]
\begin{center}
\includegraphics[height=0.15\linewidth]{codebooks/reedT.png}\\ \vspace{.5mm}
\includegraphics[height=0.35\linewidth]{codebooks/reed_gray.jpg}
\includegraphics[height=0.35\linewidth]{codebooks/reedR.png}\\
\textbf{(a) Reed}
\\ \vspace{2.5mm}
\includegraphics[height=0.15\linewidth]{codebooks/wT2.png}
\includegraphics[height=0.15\linewidth]{codebooks/wT1.png}\\ \vspace{.5mm}
\includegraphics[height=0.3\linewidth]{codebooks/w_gray.jpg}
\includegraphics[height=0.3\linewidth]{codebooks/wR.png}\\ \vspace{.5mm}
\includegraphics[height=0.27\linewidth]{codebooks/w2_gray.jpg}
\includegraphics[height=0.27\linewidth]{codebooks/w2R.png}\\
\textbf{(b) Brick walls and ivy leaves}
\\ \vspace{2.5mm}
\includegraphics[height=0.15\linewidth]{codebooks/nT1.png}
\includegraphics[height=0.15\linewidth]{codebooks/nT2.png}\\ \vspace{.5mm}
\includegraphics[height=0.3\linewidth]{codebooks/n_gray.jpg}
\includegraphics[height=0.3\linewidth]{codebooks/nR.png}\\
\textbf{(c) Net}
\end{center}
\caption{Codebook learning. See caption of Fig. \ref{fig:unsupervised}.  In the second experiment (brick walls and ivy leaves), the image on the left of the third row is the testing image. The image on the right of the third row is reconstructed by the templates learned from the training image, which is on the left of the second row.}
\label{fig:codebook0}
\end{figure}



{\em A codebook of templates and objective function}. With the above notation such as that in (\ref{eq:pack1}), now suppose we have a codebook  of $T$ templates, and let us denote them by $\{\B^{(t)}, t = 1, ..., T\}$. Then we can represent the image $\I_m$ by $K_m$ templates that are spatially translated and deformed versions of these $T$ templates in the codebook:
\begin{eqnarray}
\I_m = \sum_{k=1}^{K_m} C_{m, k} \B^{(t_{m, k})}_{X_{m, k}} + \epsilon_m, \label{eq:MODEL}
\end{eqnarray}
where each $\B^{(t_{m, k})}_{X_{m, k}}$ is obtained by translating the template of type $t_{m, k}$, i.e., $\B^{(t_{m, k})}$, to location $X_{m, k}$, and deforming it to match $\I_m$ by local max pooling in (\ref{eq:A}). For now, let us assume that the $K_m$ templates do not overlap with each other, i.e.,  $\B^{(t_{m, k})}$ span orthogonal subspaces for $k = 1, ..., K_m$, such as in the first example of Figure \ref{fig:unsupervised}. Then the dimensions that they explain are independent of each other, and the log-likelihood score is 
\begin{eqnarray}
L(\I_m \mid \B^{(t_{m, k})}_{X_{m, k}} , k = 1, ..., K_m) = \sum_{k=1}^{K_m} L(\I_m \mid \B_{X_{m, k}}^{(t_{m, k})}).\label{eq:ll}
\end{eqnarray}
Our goal is to learn the codebook of $T$ templates from the training images $\{\I_m\}$, while inferring the representation of each $\I_m$, i.e., $\{(t_{m, k}, X_{m, k}), k = 1, ..., K_m\}$, by maximizing the objective function which is defined as the sum of the log-likelihood  (\ref{eq:ll}) over all the training images $\{\I_m\}$, 
\begin{eqnarray}
\sum_{m=1}^{M} \left[\sum_{k=1}^{K_m} L(\I_m \mid \B_{X_{m, k}}^{(t_{m, k})}, k = 1, ..., K_m) \right],  \label{eq:lll}
\end{eqnarray}
subject to the constraint that for each $\I_m$, the encoding templates $\{\B^{(t_{m, k})}_{X_{m, k}}, k = 1, ..., K_m\}$ do not overlap. 


\begin{figure}[h]
\begin{center}
\includegraphics[height=0.15\linewidth]{codebooks/eT1.png}
\includegraphics[height=0.15\linewidth]{codebooks/eT2.png}
\includegraphics[height=0.15\linewidth]{codebooks/eT3.png}
\includegraphics[height=0.15\linewidth]{codebooks/eT4.png}\\ \vspace{.5mm}

\includegraphics[height=0.244\linewidth]{codebooks/e3_gray.jpg}
\includegraphics[height=0.244\linewidth]{codebooks/e5_gray.jpg}
\includegraphics[height=0.244\linewidth]{codebooks/e1_gray.jpg}
\includegraphics[height=0.244\linewidth]{codebooks/e6_gray.jpg}\\ \vspace{.5mm}
\includegraphics[height=0.244\linewidth]{codebooks/eR3.png}
\includegraphics[height=0.244\linewidth]{codebooks/eR5.png} 
\includegraphics[height=0.244\linewidth]{codebooks/eR1.png}
\includegraphics[height=0.244\linewidth]{codebooks/eR6.png}\\ \vspace{.5mm}


\end{center}
\caption{Codebook learning. A codebook of 4 models (each has 250 wavelets) is learned from 20 images. The first row displays the synthesized images ($100 \times 100$) from the 4 models. The second and third rows display 4 training images and their reconstructions by the 4 models.}
\label{fig:codebook01}
\end{figure}


\begin{figure}
\begin{center}
\includegraphics[height=0.18\linewidth]{grape/iter15_projectedTemplate1.png}
\includegraphics[height=0.18\linewidth]{grape/iter15_projectedTemplate2.png} 
\includegraphics[height=0.31\linewidth]{grape/grape.jpg}
\includegraphics[height=0.31\linewidth]{grape/train_colorsketch_image1.png}\\
\textbf{(a) Grapes}
\\ \vspace{2.5mm}
\includegraphics[height=0.18\linewidth]{lotus/iter15_projectedTemplate1.png}
\includegraphics[height=0.18\linewidth]{lotus/iter15_projectedTemplate2.png}\\ \vspace{.5mm}
\includegraphics[height=0.21\linewidth]{lotus/1.jpg}
\includegraphics[height=0.21\linewidth]{lotus/2.jpg}
\includegraphics[height=0.21\linewidth]{lotus/3.jpg}
\includegraphics[height=0.21\linewidth]{lotus/4.jpg}
\\ \vspace{.5mm}
\includegraphics[height=0.21\linewidth]{lotus/train_colorsketch_image2.png}
\includegraphics[height=0.21\linewidth]{lotus/train_colorsketch_image3.png}
\includegraphics[height=0.21\linewidth]{lotus/train_colorsketch_image4.png}
\includegraphics[height=0.21\linewidth]{lotus/train_colorsketch_image5.png}
\\ \vspace{.5mm}
\includegraphics[height=0.21\linewidth]{lotus/train_sketch_image2.png}
\includegraphics[height=0.21\linewidth]{lotus/train_sketch_image3.png}
\includegraphics[height=0.21\linewidth]{lotus/train_sketch_image4.png}
\includegraphics[height=0.21\linewidth]{lotus/train_sketch_image5.png}\\
\textbf{(b) Lotus}
%\\ \vspace{2.5mm}
%\includegraphics[height=0.18\linewidth]{catcode/iter12_projectedTemplate1.png}
%\includegraphics[height=0.18\linewidth]{catcode/iter12_projectedTemplate2.png}
%\includegraphics[height=0.18\linewidth]{catcode/iter12_projectedTemplate8.png}
%\includegraphics[height=0.18\linewidth]{catcode/iter12_projectedTemplate3.png}\\ \vspace{.5mm}
%\includegraphics[height=0.18\linewidth]{catcode/iter12_projectedTemplate5.png}
%\includegraphics[height=0.18\linewidth]{catcode/iter12_projectedTemplate7.png}
%\includegraphics[height=0.18\linewidth]{catcode/iter12_projectedTemplate6.png}
%\includegraphics[height=0.18\linewidth]{catcode/iter12_projectedTemplate4.png}\\ \vspace{.5mm}
%\includegraphics[height=0.2\linewidth]{catcode/cat16.jpg}
%\includegraphics[height=0.2\linewidth]{catcode/cat2.jpg} 
%\includegraphics[height=0.2\linewidth]{catcode/cat19.jpg} 
%\includegraphics[height=0.2\linewidth]{catcode/cat20.jpg} 
%%\includegraphics[height=0.19\linewidth]{catcode/cat18_gray.jpg} 
%\\ \vspace{.5mm}
%\includegraphics[height=0.2\linewidth]{catcode/train_sketch_image16.png}
%\includegraphics[height=0.2\linewidth]{catcode/train_sketch_image2.png} 
%\includegraphics[height=0.2\linewidth]{catcode/train_sketch_image19.png}
%\includegraphics[height=0.2\linewidth]{catcode/train_sketch_image20.png}\\
%(c) Cats
\\ \vspace{2.5mm}
\includegraphics[height=0.18\linewidth]{newcatcode/iter15_projectedTemplate1.png}
\includegraphics[height=0.18\linewidth]{newcatcode/iter15_projectedTemplate2.png}
\includegraphics[height=0.18\linewidth]{newcatcode/iter15_projectedTemplate3.png}
\includegraphics[height=0.18\linewidth]{newcatcode/iter15_projectedTemplate8.png}\\ \vspace{.5mm}
\includegraphics[height=0.18\linewidth]{newcatcode/iter15_projectedTemplate5.png}
\includegraphics[height=0.18\linewidth]{newcatcode/iter15_projectedTemplate7.png}
\includegraphics[height=0.18\linewidth]{newcatcode/iter15_projectedTemplate4.png}
\includegraphics[height=0.18\linewidth]{newcatcode/iter15_projectedTemplate6.png}\\ \vspace{.5mm}
\includegraphics[height=0.19\linewidth]{newcatcode/cat2.jpg}
\includegraphics[height=0.19\linewidth]{newcatcode/cat8.jpg} 
\includegraphics[height=0.19\linewidth]{newcatcode/cat18.jpg} 
\includegraphics[height=0.19\linewidth]{newcatcode/cat14.jpg} 
\includegraphics[height=0.19\linewidth]{newcatcode/cat17.jpg} 
%\includegraphics[height=0.19\linewidth]{catcode/cat18_gray.jpg} 
\\ \vspace{.5mm}
\includegraphics[height=0.19\linewidth]{newcatcode/train_sketch_image2.png}
\includegraphics[height=0.19\linewidth]{newcatcode/train_sketch_image8.png} 
\includegraphics[height=0.19\linewidth]{newcatcode/train_sketch_image18.png}
\includegraphics[height=0.19\linewidth]{newcatcode/train_sketch_image14.png}
\includegraphics[height=0.19\linewidth]{newcatcode/train_sketch_image17.png}\\
\textbf{(c) Cats}
\end{center}
   \caption{Codebook learning. In each of the 3 experiments, synthesized images ($100 \times 100$) from the models of the learned codebook are displayed together with the training images and their sketches by the learned models, where each Gabor wavelet is illustrated by a bar, and the templates appear in different colors (red and green) or with their bounding boxes (in green). (a) Grape experiment: each model has 37 wavelets, learned from 1 image. (b) Lotus experiment: each model has 30 wavelets, learned from 7 images. (c) Cat experiment: each model has 40 wavelets, learned from 20 images.}
\label{fig:codebook}
\end{figure}

{\em Codebook learning algorithm}. To initialize the unsupervised learning algorithm, we first learn the codebook of  templates from randomly cropped image patches. Specifically, for each $\B^{(t)}$, we randomly cropped some image patches from training images, and then we learn $\B^{(t)}$ and the associated parameters $\lambda^{(t)}$ from these image patches using the two-stage algorithm described in the previous sections. Then we iterate the following two steps that seek to maximize the objective function (\ref{eq:lll}):

(1) {\em Image encoding by template matching pursuit}.  This step assumes that the codebook $\{\B^{(t)}\}$ is given, and it seeks to maximize (\ref{eq:lll}) over the encoding of each image $\I_m$ by $\{(t_{m, k}, X_{m, k}), k = 1, ..., K_m\}$. Specifically, for each template in the codebook, we scan it over each image $\I_m$ and compute the log-likelihood score, i.e., compute $\R_m^{(t)}(X) = L(\I_m \mid \B^{(t)}_{X})$ for all $t$ and $X$. Starting from $k=1$, we sequentially select $(X_{m, k}, t_{m, k}) = \arg \max_{X, t} \R_{m}^{(t)}(X)$ subject to the constraint that $\B_{X_{m, k}}^{(t_{m, k})}$ does not overlap with previously selected templates and that the log-likelihood score of the selected template is above a threshold such as 0. 

(2) {\em Template re-learning}. This step assumes that the encoding of each image $\I_m$, i.e.,   $\{(t_{m, k}, X_{m, k})\}$, is given, and it seeks to maximize   (\ref{eq:lll})  by re-learn the codebook of templates $\{(\B^{(t)}, \lambda^{(t)}),  t = 1, ..., T\}$. Specifically, for each template $t$ in the codebook, we re-learn $(\B^{(t)}, \lambda^{(t)})$ from the image patches currently encoded by this template  using the two-stage learning algorithm. 

The above algorithm is a greedy algorithm for maximizing the objective function (\ref{eq:lll}). In fact, it can be considered a combination of detection and clustering studied in the previous sections. Even though the initial templates are random and meaningless, meaningful templates can usually be learned after a small number of iterations. These templates seek to explain different patterns in the observed images. According to our experience, meaningful templates can usually be learned regardless of the starting point of the algorithm. 

In practical implementation of the above learning algorithm, we allow the templates to vary their rotations and scales in addition to spatial translation. We also allow the templates to have limited overlap, that is, after each template is selected, it only inhibits other templates within a limited distance from its center. Our experiences show that prohibiting overlap between the selected templates can result in parts of images left unexplained. Allowing limited overlap can avoid this problem. 

In the template matching pursuit process, when a template is selected, it explains away part of the residual image by least squares projection. So after the template matching pursuit process, the observed images are reconstructed according to (\ref{eq:MODEL}). 

The number of templates in the codebook as well as the numbers of basis functions in the templates can be selected by BIC-like criteria, as suggested by \cite{ABC}. In this paper, we hand picked these parameters. 


{\em Experiment 6:  Unsupervised learning of codebooks}. We can learn a codebook of sparse FRAME models from non-aligned images without annotation. Figure \ref{fig:codebook0} illustrates 3 experiments of codebook learning. In each experiment, the images on the top row are synthetic images generated by the learned models. The input image is shown on the left of the second row. The image on the right of the second row is the reconstructed image using the learned templates. In the second experiment of brick walls and ivy leaves, the templates are learned from the training image in the second row, and they can be used to reconstruct the testing image in the third row.  Figure \ref{fig:codebook01} displays another example of a codebook learned from multiple images. Figure \ref{fig:codebook} displays results from another set of experiments, where for the sake of efficiency, we select $n=40$ Gabor wavelets of a single scale, so the synthesized images mainly capture the edge patterns. Each experiment displays $100 \times 100$ images synthesized by the models in the learned codebook, together with the training image and the sketch of the image by the learned models (in different colors in the first two experiments or with green bounding boxes in the last two experiments). There is one training image in the first experiment, while there are multiple training images in the other two experiments. As to the running time, for the lotus example, each encoding and re-learning iteration takes about 2.6 minutes. We run 15 iterations. 


\subsection{Using learned codebooks for object classification} \label{sect:classification} 

The learned codebook of sparse FRAME models can serve as ``words" in the ``bag-of-word" method for object classification. Suppose we have a codebook of $T$ models $\{\B^{(t)}, t = 1, ..., T\}$ learned from training images. For each image $\I_m$, we denote $\R_m^{(t)}(X, S, A) = L(\I_m\mid\B^{(t)}_{X, S, A})$ as the log-likelihood of $\B^{(t)}$ at location $X$, scale $S$, and orientation $A$. Both $S$ and $A$ are assumed to take values within a finite and properly discretized range. Let 
\begin{eqnarray}
r^{(t)}_m (A) = \max(\max_{X,S} R^{(t)}_m (X, S, A),0)
\end{eqnarray}
 be the maximum score at orientation $A$. Then each image $\I_m$ can be represented by a $T\times N_A$-dimensional feature vector $(r^{(t)}_m (A), t = 1, ..., T, \forall A)$, where $N_A$ is the number of possible values $A$ can take. After extracting features, we can use any discriminative method to train classifiers (e.g. linear logistic regression or SVM \cite{svm}) on such feature vectors for object classification. Spatial pyramid matching (SPM) \cite{SPM} can also be  utilized to further boost the classification performance.

 
{\em Experiment 7: Binary classification}. We evaluate the above ``bag-of-word" representation extracted by a codebook of sparse FRAME templates on a binary classification task. We test it on a collection of 16 categories from Caltech-101 \cite{Caltech101}, all 5 categories from ETHZ Shape \cite{ETHZ} and all 3 categories from Graz-02 \cite{Graz} datasets. The task is to separate each category from a negative background category. We resize all images to $150 \times 150$ pixels without changing their aspect ratios and convert them to grey level images. We randomly choose 30 positive and 30 negative images respectively as training data, and keep the rest as testing data. For Caltech-101 and Graz-02, negative images are chosen from the background category, while for ETHZ, negative examples are chosen from images other than the target category. For each category, we learn a codebook of $T = 10$ sparse FRAME templates. Each template is of the size $100 \times 100$ and has $n = 40$ wavelets. We set scale $S \in \{0.8,1,1.2\}$ and orientation $A \in \{\pm 1,0\} \times \pi/16$. Binary classification is done with linear logistic regression with regularization by $\ell_2$ norm \cite{LIBLINEAR}. We compare our results with those obtained by SIFT \cite{SIFT} features and SVM classifier, where SIFT features are quantized into ``words'' by k-means clustering (K = 50, 100, 500) and fed into linear or kernel SVM. The best result among these six combinations (3 numbers of words $\times$ 2 types of SVM) is then reported. Table \ref{table::binary} shows the comparison results of the binary classification experiments. All experiments are repeated five times with different randomly selected training and testing images, and the average accuracies and the $95\%$ confident intervals are reported. It can be seen that our method generally outperforms the SIFT + SVM method, despite the fact that we use much smaller codebooks (10 ``words'' versus 50, 100, 500 ``words''). 


\begin{table*}
\begin{center}
\caption{Accuracies ($\%$) on binary classification tasks for $24$ categories from Caltech-101, ETHZ Shape and Graz-02 data sets. } \label{table::binary}
\resizebox{.75\textwidth}{!}{
{\begin{tabular}{|lcr|lcr|}
\hline
Datasets   & SIFT+SVM      & Our method     & Datasets       & SIFT+SVM        &   Our method         \\
\hline
Watch      & {\bf 90.1$\pm$1.0}  & 89.1$\pm$1.6   & Sunflower      & 76.0$\pm$2.5    & {\bf 89.6$\pm$3.7}         \\
Laptop     & 73.5$\pm$5.3  & {\bf 89.8$\pm$2.7}   & Chair          & 62.5$\pm$5.0    & {\bf 82.9$\pm$4.7}         \\
Piano      & 84.5$\pm$4.2  & {\bf 93.8$\pm$2.6}   & Lamp           & 61.5$\pm$4.5    & {\bf 86.6$\pm$4.3}         \\
Ketch      & 82.2$\pm$0.8  & {\bf 83.3$\pm$6.5}   & Dragonfly           & 66.0$\pm$4.0    & {\bf 89.9$\pm$5.7}       \\
Motorbike  & {\bf 93.9$\pm$1.2}  & 92.2$\pm$2.9   & Umbrella       & 73.4$\pm$4.4    & {\bf 90.0$\pm$0.7}       \\
Guitar     & 70.0$\pm$2.4  & {\bf 77.3$\pm$6.3}   & Cellphone      & 68.7$\pm$5.1    & {\bf 95.7$\pm$1.8}         \\
Schooner   & 64.3$\pm$2.2  & {\bf 87.7$\pm$2.8}   & Face           & 91.8$\pm$2.3    & {\bf 94.4$\pm$2.3}   \\
Ibis       & 67.8$\pm$6.0  & {\bf 85.3$\pm$2.7}   & Starfish       & 73.1$\pm$6.7    & {\bf 90.0$\pm$2.3} \\
ETHZ-Bottle & 68.6$\pm$3.2  &  {\bf 77.5$\pm$5.6}  & ETHZ-Cup        & {\bf 66.0$\pm$3.3}    &  62.5$\pm$3.0\\
ETHZ-Swans  & 64.2$\pm$1.5  & {\bf 74.0$\pm$7.5}   & ETHZ-Giraffes   & 61.5$\pm$6.4    & {\bf 73.3$\pm$4.8}\\
ETHZ-Apple  & 55.0$\pm$1.8  & {\bf 65.8$\pm$6.1}  &
Graz02-Person& {\bf 70.4$\pm$1.2}  &  68.2$\pm$3.8\\   Graz02-Car       & {\bf 64.0$\pm$6.7}    & 59.6$\pm$5.5 & Graz02-Bike      & 68.5$\pm$2.8    &  {\bf 71.3$\pm$5.1}\\
\hline
\end{tabular}}}
\end{center}

\end{table*}


{\em Experiment 8:  Multi-class classification}. Our second set of experiments is on the LHI-Animal-Faces dataset \cite{HIT}, which consists of around 2200 images of 20 categories of animal or human faces. We randomly select half of the images per class for training and the rest for testing. We learn a codebook of 10 sparse FRAME models for each category in an unsupervised way. We then combine the codebooks of all the categories (in total $20 \times 10 = 200$ codewords). The maps of the template matching scores from the models in the combined codebook are computed for each image, and they are then fed into SPM, which equally divides an image into 1, 4,  16 areas, and concatenates the maximum scores at different image areas into a feature vector. We use multi-class SVM  to  train image classifiers based on the feature vectors, and then evaluate the classification accuracies of these classifiers on the testing data using the one-versus-all rule. Our classification rate is $79.4\%$. For comparison, Table \ref{table::multiclass} lists 4 published results \cite{HIT}  on this dataset obtained by other methods: (a) HoG feature trained with SVM, (b) Hybrid Image Template (HIT) \cite{HIT}, (c) multiple transformation invariant HITs (Mixture of HIT) \cite{HIT}, and (d) part-based HoG feature trained with latent SVM \cite{LSVM}. Our method outperforms the other methods in terms of classification accuracy on this dataset. 


%\begin{table}
%\begin{center}
%\caption{Classification accuracies  ($\%$) on the animal faces dataset.} \label{table::multiclass}
%\resizebox{.48\textwidth}{!}{
%{\begin{tabular}{|c|c|c|c|c|}
%\hline
%HoG+SVM  & HIT & Mixture of & Part-based  &  Our         \\
% &  & HIT  & LSVM   &     method    \\
%\hline
%70.8     & $71.6$  & $75.6$   & $77.6$      & ${\bf 79.4}$         \\
%\hline
%\end{tabular}}}
%\end{center}
%\end{table}

\begin{table}
\begin{center}
\caption{Classification accuracies  ($\%$) on the animal faces dataset.} \label{table::multiclass}
\resizebox{.49\textwidth}{!}{
{\begin{tabular}{|c|c|c|c|c|}
\hline
HoG+SVM  & HIT & Mixture & Part-based  & Our         \\
 &  & of HIT  & LSVM   &    method    \\
\hline
70.8     & $71.6$  & $75.6$   & $77.6$  & ${\bf 79.4}$         \\
\hline
\end{tabular}}}
\end{center}
\end{table}

\begin{table*}[htdp]
\caption{Results on the domain transfer experiment}\label{table::transfer}
\begin{center}
\subfloat[Classification accuracies  ($\%$) on single source four domains benchmark ( C: caltech, A: amazon, D: DSLR, W: webcam)]{
\resizebox{.9\textwidth}{!}{
\begin{tabular}{|c|c|c|c|c|c|c|c|c|}
\hline
Method & C$\rightarrow$A & C$\rightarrow$D & A$\rightarrow$C &  
A$\rightarrow$W & W$\rightarrow$C & W$\rightarrow$A & D$\rightarrow$A  
& D$\rightarrow$W\tabularnewline
\hline
\hline
Metric \cite{1} & 33.7$\pm$0.8 & 35.0$\pm$1.1 & 27.3$\pm$0.7 & 36.0$\pm$1.0 &
21.7$\pm$0.5 & 32.3$\pm$0.8 & 30.3$\pm$0.8 & 55.6$\pm$0.7\tabularnewline
\hline
SGF \cite{2} & 40.2$\pm$0.7 & 36.6$\pm$0.8 & 37.7$\pm$0.5 & 37.9$\pm$0.7 &
29.2$\pm$0.7 & 38.2$\pm$0.6 & 39.2$\pm$0.7 & 69.5$\pm$0.9\tabularnewline
\hline
GFK \cite{3} & 46.1$\pm$0.6 & 55.0$\pm$0.9 & 39.6$\pm$0.4 & 56.9$\pm$1.0 &
32.8$\pm$0.7 & 46.2$\pm$0.7 & 46.2$\pm$0.6 & {\bf 80.2$\pm$0.4}\tabularnewline
\hline
FDDL \cite{4} & 39.3$\pm$2.9 & 55.0$\pm$2.8 & 24.3$\pm$2.2 & 50.4$\pm$3.5 &
22.9$\pm$2.6 & 41.1$\pm$2.6 & 36.7$\pm$2.5 & 65.9$\pm$4.9\tabularnewline

\hline
MMDT \cite{MMDT} & 49.4$\pm$0.8 & 56.5$\pm$0.9 & 36.4$\pm$0.8 & 64.6$\pm$1.2 &
32.2$\pm$0.8 & 47.7$\pm$0.9 & 46.9$\pm$1.0 & 74.1$\pm$0.8\tabularnewline
\hline
SDDL \cite{SDDL} & 49.5$\pm$2.6 & {\bf 76.7$\pm$3.9} & 27.4$\pm$2.4 & {\bf 72.0$\pm$4.8} &
29.7$\pm$1.9 & 49.4$\pm$2.1 & 48.9$\pm$3.8 & 72.6$\pm$2.1\tabularnewline
\hline
Our method & {\bf 62.2$\pm$1.6} & 52.2$\pm$4.0 & {\bf 46.7$\pm$2.5} & 53.2$\pm$4.9
& {\bf 39.1$\pm$3.0} & {\bf 53.2$\pm$4.4} & {\bf 55.3$\pm$2.9} &
72.4$\pm$3.1\tabularnewline
\hline
\end{tabular}
}
}


\subfloat[Classification accuracies  ($\%$) on multiple sources three domains benchmark ]{
\resizebox{.7\textwidth}{!}{
\begin{tabular}{|c|c|c|c|c|c|}
\hline
Source & Target & SGF \cite{2} & RDALR \cite{5} & FDDL \cite{4} & Our method\tabularnewline
\hline
\hline
DLSR, amazon & webcam & 52$\pm$2.5 & 36.9$\pm$1.1 & 41.0$\pm$2.4 &
{\bf 52.2$\pm$1.4}\tabularnewline
\hline
amazon, webcam & DSLR & 39$\pm$1.1 & 31.2$\pm$1.3 & 38.4$\pm$3.4 &
{\bf 54.5$\pm$3.3}\tabularnewline
\hline
webcam, DSLR & amazon & 28$\pm$0.8 & 20.9$\pm$0.9 & 19.0$\pm$1.2 &
{\bf 32.1$\pm$1.6}\tabularnewline
\hline
\end{tabular}}

}
\end{center}

\end{table*}

{\em Experiment 9:  Domain transfer}. Classifiers learned from one domain (the source domain) may perform poorly on other domains (the target domains), because the training and testing data may not come from the same distribution. Learning domain-invariant feature representations can deal with this problem. In this experiment, we test our proposed representation for the task of domain transfer on four domain datasets, and compare with published results \cite{1} \cite{2} \cite{3} \cite{4} \cite{5} \cite{MMDT} \cite{SDDL}. The four datasets are: Amazon, Webcam,  DSLR and Caltech-256  \cite{6}. Each dataset is regarded as a domain. For the experiment with  single source training, 10 classes common to all 4 datasets are extracted: backpack, touring-bike, calculator, head-phones, computer-keyboard, laptop-101, computer-monitor, computer-mouse, coffee-mug, and video-projector. For the experiment with multiple sources training, all 31 classes in Amazon, Webcam and DSLR  are used. We use the evaluation protocol in \cite{3}. We randomly sample labeled data in the source domain as training examples, and unlabeled  data in the target domain as testing examples. We learn a combined codebook (by learning a codebook of 3 templates with $n=40$ wavelets for each category and combining them together), then use it to extract feature vectors and train classifiers by multi-class SVM using the same scheme as in Experiment 8. We evaluate the classification accuracies of these classifiers on the testing domain. For each pair of source and target domains, we report averaged accuracies on  target domains as well as standard errors. Table \ref{table::transfer} shows  the comparisons of recognition accuracies on target domains  for single source training and multiple source training, where the accuracies and standard errors are obtained from 10 repetitions. It can be  seen that our method performs significantly better than previous methods  on 8 out of 11 sub-tasks, and on-par with the best performing method  on the other sub-tasks, even though we do not make use of any domain adaptation techniques. This suggests that the learned codebooks of models capture intrinsically meaningful patterns. 



\section{Conclusion} \label{sect:discuss} 


We propose that the sparse FRAME models form the layer of representational units above the layer of wavelets sparse coding.  A sparse FRAME model makes use of wavelet sparse coding to generate image intensities, while accounting for the distribution of the coefficients of the selected wavelets as well as perturbations of their locations and orientations. 

As a generative model, the sparse FRAME model has the following characteristics. (1) It can reconstruct the training images, and reconstruction is used for selecting the basis functions. (2) It can synthesize new images, and synthesis is required for estimating parameters and calculating the normalizing constant. (3) It separates shape deformations and appearance variations. (4) It gives interpretable sketches. (5) Codebooks of models can be learned in an unsupervised manner.  (6) It combines rich traditions of harmonic analysis and Markov random field models. 

While we have shown that it is possible to learn codebooks of sparse FRAME models, much remains to be understood about learning large codebooks reliably from big training data sets. 

\section*{Reproducibility} 

{\small \url{http://www.stat.ucla.edu/~jxie/sparseFRAME.html} }
The above webpage contains the full data sets and exact parameter settings, and matlab/C code for producing the experimental results presented in this paper. 

\begin{acknowledgements}
The work is supported by NSF DMS 1310391, NSF IIS 1423305, ONR MURI N00014-10-1-0933, DARPA MSEE FA8650-11-1-7149. We thank the three reviewers for their insightful comments and valuable suggestions that have helped us improve the presentation and the content of this paper. We are grateful to one reviewer for sharing the insights on the analysis prior models. Thanks also go to an editor of the special issue for helpful suggestions. We thank Adrian Barbu for discussions. 

\end{acknowledgements}


\section{Appendix} 

\subsection{Simulation by Hamiltonian Monte Carlo} \label{sect:HMC}

To approximate $\E_{p(\I;\lambda^{(t)})}[|\langle \I, B_{x, s, \alpha}\rangle|]$ in equation (\ref{eq:learning0}), we need to draw a synthesized sample set $\{\tI_m\}$ from $p(\I;\lambda^{(t)})$ by HMC \cite{HMC2}. We can write $p(\I; \lambda)$ as $p(\I) \propto \exp(-U(\I))$, where $\I \in R^{|{\cal D}|}$ and 
\begin{eqnarray}
U(\I)=-\sum_{x, s, \alpha} \lambda_{x, s, \alpha} | \langle \I, B_{x, s, \alpha} \rangle |+ \frac{1}{2} |\I|^2
\end{eqnarray}
(assuming $\sigma^2 = 1$). In physics context, $\I$ can be regarded as a position vector and $U(\I)$ the potential energy function. To allow Hamiltonian dynamics to operate, we need to introduce an auxiliary momentum vector $\bm{\phi}\in R^{|{\cal D}|}$ and the corresponding kinetic energy function $K(\bm{\phi})=|\bm{\phi}|^2/2m$, where $m$ represents the mass. After that, a fictitious physical system described by the canonical coordinates $(\I,\bm{\phi})$ is defined, and its total energy is $H(\I,\bm{\phi})=U(\I)+K(\bm{\phi})$. Instead of sampling from $p(\I)$ directly, HMC samples from the joint canonical distribution $p(\I,\bm{\phi}) \propto \exp(-H(\I,\bm{\phi}))$, under which  $\I \sim p(\I)$ marginally and $\bm{\phi}$ follows a Gaussian distribution and is independent of $\I$.  Each time HMC draws a random sample from the marginal Gaussian distribution of $\bm{\phi}$, and then evolves according to the Hamiltonian dynamics that conserves the total energy. 

In practical implementation,  the leapfrog algorithm is used to discretize the continuous Hamiltonian dynamics as follows, with $\epsilon$ being the step-size: 
\begin{eqnarray}
&& \bm{\phi}^{(t+\epsilon/2)}=\bm{\phi}^{(t)}-(\epsilon/2)\frac{\partial U}{\partial \I}(\I^{(t)}), \\
&& \I^{(t+\epsilon)}= \I^{(t)} + \epsilon \frac {\bm{\phi}^{(t+\epsilon /2)}}{m}, \\
&& \bm{\phi}^{(t+\epsilon)}=\bm{\phi}^{(t+\epsilon/2)}-(\epsilon/2)\frac{\partial U}{\partial \I}(\I^{(t+\epsilon)}),
\end{eqnarray} 
that is, a half-step update of $\bm{\phi}$ is performed first and then it is used to compute $\I^{(t + \epsilon)}$ and $\bm{\phi}^{(t + \epsilon)}$. 

A key step in the leapfrog algorithm is the computation of the derivative of the potential energy function
\begin{eqnarray}
\frac{\partial U}{ \partial \I}=-\sum_{x, s, \alpha} \lambda_{x, s, \alpha} \mbox{sign}( \langle \I, B_{x, s, \alpha} \rangle )B_{x, s, \alpha}+ \I, 
\end{eqnarray} 
where the map of responses $r_{x, s, \alpha} = \langle \I, B_{x, s, \alpha} \rangle$ is computed by bottom-up convolution of the filter corresponding to $(s, \alpha)$ with $\I$ for each $(s, \alpha)$. Then the derivative is computed by top-down linear superposition of the basis functions:  $-\sum_{x, s, \alpha} \lambda_{x, s, \alpha} \mbox{sign}( r_{x, s, \alpha} )B_{x, s, \alpha} + \I$, which can again be computed by convolution. Both bottom-up and top-down convolutions can be carried out efficiently by GPUs. 

The discretization of the leapfrog algorithm cannot keep $H(\I, \bm{\phi})$ exactly constant, so a Metropolis acceptance/rejection step is used to correct the discretization error. Starting with the current state, $(\I,\bm{\phi})$, the new state $(\I^ \star,\bm{\phi}^ \star)$, after $L$ leapfrog steps, is accepted as the next state of the Markov chain with probability 
$
\min [1, \exp(-H(\I^ \star,\bm{\phi}^ \star)+H(\I,\bm{\phi})) ].
$
If it is not accepted, the next state is the same as the current state. 

In summary, a complete description of the HMC sampler for inhomogeneous FRAME is as follows:

(i) Generate the momentum vector $\bm{\phi}$ from its marginal distribution $p(\bm{\phi}) \propto \exp(-K(\bm{\phi}))$, which is the zero-mean Gaussian distribution with covariance matrix $m I$ ($I$ is the identity matrix).

(ii) Perform $L$ leapfrog steps to reach the new state $(\I^{\star},\bm{\phi}^{\star}).$

(iii) Perform acceptance/rejection of the proposed state $(\I^{\star},\bm{\phi}^{\star}).$

$L$, $\epsilon$, and $m$ are parameters of the algorithm, which need to be tuned to obtain good performance.

\subsection{Maximum entropy justification} \label{sect:app1}

The inhomogeneous FRAME model can be justified by the maximum entropy principle. Suppose the true distribution that generates the observed images $\{\I_m\}$ is $f(\I)$. Let $\lambda^{\star}$ solve  the population version of the maximum likelihood equation: 
\begin{eqnarray}
\E_{p(\I; \lambda)}[|\langle \I, B_{x, s, \alpha}\rangle|] = \E_{f}[|\langle \I, B_{x, s, \alpha}\rangle|], \;\forall x, s, \alpha.
\end{eqnarray}
 Let $\Omega$ be the set of all the distributions $p(\I)$ such that 
 \begin{eqnarray}
 \E_{p}[|\langle \I, B_{x, s, \alpha}\rangle|] = \E_{f}[|\langle \I, B_{x, s, \alpha}\rangle|], \; \forall x, s, \alpha.
 \end{eqnarray}
  Then $f \in \Omega$. Let $\Lambda$ be the set of all the distributions $\{p_{\lambda}, \forall \lambda\}$, where $p_{\lambda}(\I) = p(\I; \lambda)$. Then $q \in \Lambda$ since $q(\I) = p(\I; \lambda = 0)$. Thus $p_{\lambda^\star}$ is the intersection between $\Lambda$ and $\Omega$. In Figure \ref{fig:ME}, $\Lambda$ and $\Omega$ are illustrated by  blue and green  curves respectively, where each point on the curves is a probability distribution. The two curves $\Lambda$ and $\Omega$ are ``orthogonal'' in the sense that for any $p_{\lambda} \in \Lambda$ and for any $p \in \Omega$, it can be easily proved that the Pythagorean property 
  \begin{eqnarray}
  {\rm KL}(p || p_{\lambda}) = {\rm KL}(p || p_{\lambda^{\star}}) + {\rm KL}(p_{\lambda^{\star}}||p_{\lambda})
  \end{eqnarray}
   holds \cite{Pietra_Pietra}, where ${\rm KL}(p||q)$ is the Kullback-Leibler divergence from $p$ to $q$. This Pythagorean property leads to the following dual properties of $p_{\lambda^{\star}}$:
   
   (1) Maximum likelihood:  Among all $p_{\lambda} \in \Lambda$, $p_{\lambda^{\star}}$ achieves the minimum of ${\rm KL}(f||p_{\lambda})$. 
   
   (2) Maximum entropy or minimum divergence: Among all $p \in {\Omega}$, $p_{\lambda^{\star}}$ achieves the minimum of ${\rm KL}(p||q)$. Thus $p_{\lambda^{\star}}$ can be considered the minimal modification of the reference distribution $q$ to match the statistical properties of the true distribution $f$. 
   
   The above justification is also true for the sparse FRAME model. 

\begin{figure}[h]
\begin{center}
\includegraphics[height=0.4\linewidth]{graph.png}
\end{center}
\caption{Illustration of the maximum entropy principle. Each curve illustrates a set of probability distributions. $\Omega$ is the set of distributions that reproduce statistical properties of filter response of the true distribution $f$. $\Lambda$ is the set of distributions of the model. The two curves are orthogonal to each other in the sense of the Pythagorean property of the Kullback-Leibler divergences. So $p_{\lambda^{\star}}$ can be considered the minimal modification of the reference distribution $q$ to match the statistical properties of $f$. }
\label{fig:ME}
\end{figure}

For sparsification, in principle, we can select $B_{x_i, s_i, \alpha_i}$ sequentially using a procedure like projection pursuit \cite{Friedman_projection_pursuit} or filter pursuit \cite{Minimax}. Suppose we have selected $k$ basis functions $(B_{x_i, s_i, \alpha_i}, i = 1, ..., k)$, and let $p_k$ be the fitted model with the corresponding $\lambda = (\lambda_i, i = 1, ..., k)$ estimated by MLE. Suppose we are to select the next basis function $B_{x_{k+1}, s_{k+1}, \alpha_{k+1}}$.  Let $p_{k+1}$ be the fitted model. Then we want to minimize ${\rm KL}(f||p_{k+1}) = {\rm KL}(f||p_{k}) - {\rm KL}(p_{k+1}||p_k)$, that is, we want to maximize ${\rm KL}(p_{k+1}||p_k)$, which serves as the pursuit index. The problem with such a procedure is that each time we need to fit $p_k$ which involves MCMC computation, and the pursuit index is also difficult to compute. So we choose to pursue a different approach by exploring the connection between sparse FRAME and the shared sparse coding. 


\subsection{Sparse FRAME and shared sparse coding} \label{sect:app2}

{\em From sparse FRAME to shared sparse coding}.  Let us assume that the reference distribution $q(\I)$ in the sparse FRAME model (\ref{eq:sparseFRAME}) is a Gaussian white noise model so that the pixel intensities follow ${\rm N}(0, \sigma^2)$ independently. For sparse FRAME, it is natural to assume that the number of selected basis functions $n$ is much less than the number of pixels in $\I$, i.e., $n \ll |{\cal D}|$, where ${\cal D}$ is the image domain. For notational convenience, we can make $\I$ and $B_i = B_{x_i, s_i, \alpha_i}$, $i = 1, ..., n$ into $|{\cal D}|$-dimensional vectors, and let $\B = (B_1, ..., B_n)$ be the resulting $|{\cal D}| \times n$ matrix. 

The connection between sparse FRAME and shared sparse coding is most evident if we temporarily assume that the selected basis functions $(B_i, i = 1, ..., n)$ are orthogonal (with unit $\ell_2$ norm as assumed before). Extension to non-orthogonal $\B$ is straightforward but requires tedious notation (such as $(\B^{T}\B)^{-1}$). For $\B$, we can construct $\bar{n} = |{\cal D}| - n$ basis vectors of unit norm $\bb_1, ..., \bb_{\bar{n}}$ that are orthogonal to each other and that are also orthogonal to $(B_i, i = 1, ..., n)$. Thus each image $\I = \sum_{i=1}^{n} r_i B_i + \sum_{i=1}^{\bar{n}} \bar{r}_i \bb_i$, where $r_i = \langle \I, B_i \rangle$, and $\bar{r}_i = \langle \I, \bb_i\rangle$. So we have the linear additive model $\I = \sum_{i=1}^{n} r_i B_i + \epsilon$, with $\epsilon = \sum_{i=1}^{\bar{n}} \bar{r}_i \bb_i$ being the least squares residual image. 

Under the Gaussian white noise $q(\I)$, $r_i$ and $\bar{r}_i$ are all independent ${\rm N}(0, \sigma^2)$ random variables because of the orthogonality of $(\B, \bB)$. Let $R$ be the column vector whose elements are $r_i$, and $\bar{R}$ be the column vector whose elements are $\bar{r}_i$. Then under the sparse FRAME model (\ref{eq:sparseFRAME}), only the distribution of $R$ is modified during the change from $q(\I)$ to $p(\I; \B, \lambda)$, which changes the distribution of $R$ from Gaussian white noise $q(R)$ to 
\begin{eqnarray}
p(R; \lambda) = \frac{1}{Z(\lambda)}\exp\left(\sum_{i=1}^{n} \lambda_i |r_i|\right) q(R),
\end{eqnarray}
 while the distribution of the residual coordinates $\bar{R}$ remains Gaussian white noise, and $R$ and $\bar{R}$ remain independent. That is, $p(R, \bar{R}; \lambda) = p(R; \lambda) q(\bar{R}) $.

Thus the sparse FRAME model implies a linear additive model $\I = \sum_{i=1}^{n} r_i B_i + \epsilon$, where $R \sim p(R; \lambda)$ and $\epsilon$ is a Gaussian white noise in the $\bar{n}$-dimensional residual space, and $\epsilon$ is independent of $R$. If we observe independent training images $\{\I_m, m = 1, ..., M\}$ from the model, then $\I_m = \sum_{i=1}^{n} r_{m, i} B_i + \epsilon_m$, i.e., $\{\I_m\}$ share a common set of basis functions $\B = (B_i, i = 1, ..., n)$ that provide sparse coding for multiple images simultaneously. 


{\em From shared sparse coding to sparse FRAME}. Conversely, suppose we are given a shared sparse coding model of the form $\I =\sum_{i=1}^{n} c_i B_i + \epsilon =  \B C + \epsilon$, where $C$ is a column vector whose components are $c_i$. Assume $C \sim p(C)$ and $\epsilon \sim {\rm N}(0, I \sigma^2)$, where $I$ is the $|{\cal D}|$-dimensional identity matrix, and $\epsilon$ and $C$ are independent. Let $\delta = \B^T \epsilon$,  each component of which $\delta_i = \langle \epsilon, B_i\rangle \sim {\rm N}(0, \sigma^2)$ independently. Then we can write $\I = \B R + \bB \bar{R}$, where $R = C + \delta$, and $\bar{\epsilon} = \bB \bar{R}$ is the projection of $\epsilon$ onto the space of $\bB$. Let $\tilde{p}(R)$ be the density of $R = C+ \delta$, which is obtained by convolving $p(C)$ with Gaussian white noise density. Then $p(\I) = \tilde{p}(R) q(\bar{R}) = q(\I) \tilde{p}(R) /q(R)$ since $q(\I) = q(R)q(\bar{R})$ under Gaussian white noise model ($d\I = dR d\bar{R}$ under orthogonality so there is no Jacobian term). If we choose to model $\tilde{p}(R)/q(R) = \exp\left(\sum_{i=1}^{n} \lambda_i |r_i|\right)/Z(\lambda)$, we arrive at the sparse FRAME model. 


{\em Selection of basis functions}. For orthogonal $\B$, as shown above, the probability density $p(\I; \B, \lambda) = q(\bar{R}) p(R; \lambda) = q(\bar{R}) q(R) \exp\left(\sum_{i=1}^{n} \lambda_i |r_i|\right)/Z(\lambda)$. Given a set of training images $\{\I_m, m = 1, ..., M\}$, and for a candidate set of basis functions $\B = (B_i, i= 1, ..., n)$, we can estimate $\lambda = (\lambda_i, i = 1, ..., n)$ by MLE, giving us $\lambda^{\star}$, and the resulting log-likelihood is 
\begin{eqnarray}
   && \sum_{m=1}^{M} \log p(\I_m; \B, \lambda^{\star}) \nonumber \\
   &=&  \sum_{m=1}^{M} \left[\log q(\bar{R}_m) + \log p(R_m; \lambda^{\star})\right] \\
    &=&- \frac{1}{2\sigma^2}  \sum_{m=1}^{M}  ||\I_m - \B R_m||^2  - \frac{M \bar{n}}{2} \log (2\pi\sigma^2) \label{eq:Gauss}\\
      &&+ \sum_{m=1}^{M} \log p(R_m; \lambda^{\star}). \label{eq:R}
\end{eqnarray}
Suppose we are to choose a $\B$ from a collection of candidates. Ideally we should maximize the sum of (\ref{eq:Gauss}) and (\ref{eq:R}). We may interpret (\ref{eq:Gauss}) to be the negative coding length of the residual image $\epsilon$ by the Gaussian white noise model, and interpret (\ref{eq:R}) to be the negative coding length of the coefficients $R_m$ by the fitted model $p(R; \lambda^{\star})$. If $\sigma^2$ is small, (\ref{eq:Gauss}) can be more important, while the coding length of $R_m$ for different $\B$ may not differ too much in comparison. So we choose to seek a $\B$ to maximize only (\ref{eq:Gauss}) or equivalently minimize the overall reconstruction error $\sum_{m=1}^{M}  ||\I_m - \B R_m||^2$. This reflects a two-stage strategy in modeling $\{\I_m\}$. First, we find a set of basis functions $\B$ to reconstruct $\{\I_m\}$ as accurately as possible. Then we fit a statistical model for the reconstruction coefficients. 

 {\em Non-orthogonality}. Even if $\B$ is not orthogonal, which is the case in our work, the connection between the sparse FRAME and shared sparse coding still holds. The responses $R = \B^{T} \I$, but the reconstruction coefficients become $C = (\B^{T}\B)^{-1}R$. The projection of $\I$ onto the subspace spanned by $\B$ is $\B C$. We can continue to assume the implicit $\bB = (\bb_i, i = 1, ..., \bar{n})$ to be orthonormal, and that they are orthogonal to the columns of $\B$. We can also continue to let $\bar{R} = \bB^{T}\I$. In this setting, $R$ and $\bar{R}$ are still independent under the Gaussian white noise model $q(\I)$ because $\B$ and $\bB$ are still orthogonal to each other. Under the sparse FRAME model (\ref{eq:sparseFRAME}), it is still the case that only the distribution of $R$ is modified during the change from $q(\I)$ to $p(\I; \B, \lambda)$, while the distribution of $\bar{R}$ remains white noise and is independent of $R$. The distribution of $R$ implies a distribution of the reconstruction coefficients $C$ because they are linked by a linear transformation.  In fact, the distribution of $C$ is:
 \begin{eqnarray}
 p_C(C; \lambda) = \frac{1}{Z(\lambda)} \exp(\langle \lambda, |\B^{T}\B C|\rangle) q_C(C),
 \end{eqnarray}
  where $q_C(C)$ is the distribution of $C$ under the reference distribution $q(\I)$, and for a vector $u$, $|u|$ means the vector obtained by taking the absolute values of $u$ component-wise.  Now the distributions of $R$ and $C$ involve the Jacobian terms such that $dR d\bar{R}  = |{\rm det}(\B^{T}\B)|^{1/2} d\I = |{\rm det}(\B^{T}\B)| dC d\bar{R}$. In fact $p(\I; \B,  \lambda) = p_C(C; \lambda)q_{\bar{R}}(\bar{R}) |\det(\B^{T}\B)|^{-1/2}$. By the same logic as in (\ref{eq:Gauss}) and (\ref{eq:R}), we still want to find $\B$ to minimize the overall reconstruction error $\sum_{m=1}^{M}\|\I_m - \B C_m\|^2$. 
 
Under the shared sparse coding model, it is tempting to model the coefficients $C$ of the selected basis functions directly. However, $C$ is still a multi-dimensional vector, and direct modeling of $C$ can be difficult. One may assume that the components of $C$ are statistically independent for simplicity, but this assumption is unlikely to be realistic. So after selecting the basis functions, we choose to model the image intensities by the inhomogeneous FRAME model. Even though this model only matches the marginal distributions of filter responses of the selected basis functions, the model does not assume that the responses are independent. 


% BibTeX users please use one of
%\bibliographystyle{spbasic}      % basic style, author-year citations
%\bibliographystyle{spmpsci}      % mathematics and physical sciences
%\bibliographystyle{spphys}       % APS-like style for physics
%\bibliography{}   % name your BibTeX data base

% Non-BibTeX users please use
\begin{thebibliography}{}
%
% and use \bibitem to create references. Consult the Instructions
% for authors for reference list style.
%

\bibitem{BM}
D. H. Ackley, G. E. Hinton, T. J. Sejnowski. A learning algorithm for Boltzmann machines. {\em Cognitive Science}, {\bf 9}, 147--169. 1985. 

\bibitem{subspace}
A. Adler, M. Elad, and Y. Hel-Or. Probabilistic Subspace Clustering via Sparse Representations, {\em IEEE Signal Processing Letters}, {\bf 20},  63--66,  2013.

\bibitem{KSVD}
M. Aharon, M. Elad, and A. M. Bruckstein. The K-SVD: an algorithm for designing of overcomplete dictionaries for sparse representation, {\em IEEE Trans. Signal Process}., {\bf 54}, 4311--4322, 2006.

\bibitem{Bengio}
Y. Bengio, A. C. Courville, and P. Vincent. Representation learning: A review and new perspectives. {\em TPAMI}., {\bf 35}, 1798--1828, 2013.


\bibitem{sparsity}
A. M. Bruckstein, D. L. Donoho, and M. Elad. From sparse solutions of systems of equations to sparse modeling of signals and images. {\em SIAM Review}, {\bf 51},  34--81, 2009.

\bibitem{basis pursuit}
S. S. Chen, D. L. Donoho, and M. A. Saunders. Atomic decomposition by basis pursuit.  {\em SIAM Review}, {\bf 43}, 129--159, 2001.

\bibitem{Chen_Huo}
J. Chen and X. Huo. Sparse representations for multiple measurements vectors (mmv) in an overcomplete dictionary. {\em Proceedings of IEEE International Conference on Acoustics, Speech, and Signal Processing},  {\bf 4}, 257--260, 2005.

\bibitem{HOG}
N. Dalal and B. Triggs. Histograms of oriented gradients for human detection. {\em CVPR}, 2005.

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

\bibitem{HMC2}
S. Duane, A. D. Kennedy, B. J. Pendleton, and D. Roweth. Hybrid Monte Carlo, {\em Physics Letters}, {\bf 195}, 216--222, 1987.

\bibitem{Elad}
M.  Elad. {\em Sparse and Redundant Representations: From Theory to Applications in Signal and Image Processing}, Springer, 2010.

\bibitem{Denoising}
M. Elad and M. Aharon. Image denoising via sparse and redundant representations over learned dictionaries. {\em IEEE Trans. Image Process.}, {\bf 15}, 3736--3745, 2006.

\bibitem{AnalysisSyn}
M. Elad, P. Milanfar, and R. Rubinstein. Analysis versus synthesis in signal priors. {\em Inverse Problems}, {\bf 23}, 2007. 

\bibitem{LIBLINEAR}
R. E. Fan, K. W. Chang, C. J. Hsieh, X. R. Wang, and C. J. Lin. LIBLINEAR: A library for large linear classiÞcation. {\em Journal of Machine Learning Research}, {\bf 9}, 1871--1874, 2008.

\bibitem{Caltech101}
L. Fei-Fei, R. Fergus, and P. Perona. Learning generative visual models from few training examples: an incremental Bayesian approach tested on 101 object categories. {\em CVPR Workshop}, 2004.

\bibitem{LSVM}
P. Felzenszwalb, R. Girshick, D. McAllester, and D. Ramanan. Object detection with discriminatively trained part-based models. {\em IEEE Trans. PAMI}, {\bf 32}, 1627--1645, 2010.

\bibitem{ETHZ}
V. Ferrari, F. Jurie, and C. Schmid. From images to shape models for object detection. {\em IJCV}, {\bf 87}, 284-303, 2010.

\bibitem{Fidler}
S. Fidler, M. Boben, and A. Leonardis. Similarity-based cross-layered hierarchical representation for object categorization.  {\em CVPR}, 2008.

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

\bibitem{bridgeSampling}
A. Gelman and X. L. Meng. Simulating normalizing constants: from importance sampling to bridge sampling to path sampling. {\em Statistical Science}, {\bf 13}, 163--185, 1998.

\bibitem{Gibbs}
S. Geman, and D. Geman. Stochastic relaxation, Gibbs distribution, and the Bayesian restoration of images. {\em IEEE Trans. PAMI}, {\bf 6}, 721--741, 1984.

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

\bibitem{3}
B. Gong, Y. Shi, F. Sha. and  K. Grauman. Geodesic flow kernel for unsupervised domain adaptation. {\em CVPR},  2012. 

\bibitem{2}
R. Gopalan, R. Li, and R. Chellappa. Domain adaptation for  object recognition: an unsupervised approach. {\em ICCV}, 2011. 

\bibitem{6}
G. Griffin, A. Holub, and P. Perona. Caltech-256 object category  dataset. Technical report, Caltech, 2007. 

\bibitem{PoE}
G. E. Hinton. Training products of experts by minimizing contrastive divergence. {\em Neural Computation}, {\bf 14}, 1771--1800, 2002. 

\bibitem{Hinton}
G. E. Hinton, S. Osindero, and Y. Teh.  A fast learning algorithm for deep belief nets. {\em Neural Computation}, {\bf 18}, 1527--1554, 2006.

\bibitem{MMDT}
J. Hoffman, E. Rodner, J. Donahue, K. Saenko, and T. Darrell. Efficient learning of domain-invariant image representations. {\em ICLR}, 2013.

\bibitem{ABC}
Y. Hong, Z. Si, W. Hu, S. C. Zhu, and Y. N. Wu, Unsupervised learning of compositional sparse code for natural image representation. {\em Quarterly of Applied Mathematics}, {\bf 72}, 373--406, 2013. 

\bibitem{5}
I. Jhou, D. Liu, D. T. Lee, and S. Chang. Robust visual  domain adaptation with low-rank reconstruction. {\em CVPR}, 2012. 

\bibitem{SPM}
S. Lazebnik, C. Schmid, and J. Ponce. Beyond bags of features: spatial pyramid matching for recognizing natural scene categories. {\em CVPR}, 2006.

\bibitem{LeeNg}
H. Lee, R. Grosse, R Ranganath, and A. Y. Ng. Convolutional deep belief networks for scalable unsupervised learning of hierarchical representations. {\em ICML}, 2009.

\bibitem{LiuCe}
C.  Liu, S.-C. Zhu, and H.-Y. Shum. Learning Inhomogeneous Gibbs Model of Faces by Minimax Entropy.  {\em ICCV}, 281-287, 2001.

\bibitem{Multitask}
K. Lounici, A. B. Tsybakov, M. Pontil, and S. A. van de Geer. Taking advantage of sparsity in multi-task learning. {\em Proceedings of the 22nd Conference on Learning Theory}, 2009

\bibitem{SIFT}
D. Lowe. Distinctive image features from scale-invariant keypoints. {\em IJCV}, {\bf 60}, 91--110, 2004.

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

\bibitem{Graz}
M. Marszalek and C. Schmid. Accurate object localization with shape masks. {\em CVPR}, 2007.

\bibitem{AnalysisModel}
S. Nama, M. E. Daviesb, M. Eladc, R. Gribonval. The cosparse analysis model and algorithms. {\em Applied and Computational Harmonic Analysis}, {\bf 34}, 30-56, 2013.

\bibitem{ImportanceSampling}
R.  Neal. Annealed importance sampling. {\em Statistics and Computing}, {\bf 11}, 125--139, 2001.

\bibitem{HMC}
R.  Neal. MCMC using Hamiltonian dynamics. {\em Handbook of Markov Chain Monte Carlo}, 2011. 

\bibitem{Supportunion}
G. Obozinski, M. J. Wainwright, and M. I. Jordan, Support union recovery in high-dimensional multivariate regression, {\em Annals of Statistics}, {\bf 39}, 1--47, 2011.

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

\bibitem{OMP}
Y. C. Pati, R. Rezaiifar, and P. S. Krishnaprasad. Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition, {\em  The 27th Asilomar Conference on Signals, Systems and Computers}, 40--44, 1993.

\bibitem{Pietra_Pietra}
S. D. Pietra, V. D. Pietra, and J. Lafferty. Inducing features of random fields.  {\em TPAMI}, {\bf 19}, 380--393, 1997.

\bibitem{Ranzato}
M. Ranzato and G. E. Hinton. Modeling pixel means and covariances using factorized third-order Boltzmann machines. {\em CVPR}, 2010.

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

\bibitem{FoE}
S. Roth and M. Black. Fields of experts. {\em IJCV},  {\bf 82}, 205--229, 2009.


\bibitem{doubleSparsity}
R. Rubinstein, M. Zibulevsky, and M. Elad, Double Sparsity: Learning Sparse Dictionaries for Sparse Signal Approximation, {\em IEEE Transactions on Signal Processing}, {\bf 58}, 1553--1564, 2010.

\bibitem{1}
K. Saenko, B. Kulis, M. Fritz, and T. Darrell. Adapting visual category models to new domains. {\em ECCV}, 2010. 

\bibitem{SDDL}
S. Shekhar, V. M. Patel, H. V. Nguyen, and R. Chellappa. Generalized domain adaptive dictionaries. {\em CVPR}, 2013.

\bibitem{HIT}
Z. Si and S. C. Zhu. Learning Hybrid image Template (HIT) by Information Projection.
{\em IEEE Trans. PAMI}, {\bf 34}, 1354--1367, 2012.


\bibitem{RBM}
P. Smolensky. Information processing in dynamical systems: foundations of harmony theory. In D. E. Rumelhart and J. L. McClelland, editors, {\em Parallel Distributed Processing}, volume 1, chapter 6, pages 194--281. MIT Press, Cambridge, 1986. 

\bibitem{Energy}
Y. W. Teh,  M. Welling, S. Osindero, and G. E. Hinton. Energy-based models for sparse overcomplete representations. {\em Journal of Machine Learning Research}, {\bf 4}, 1235--1260, 2003.

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

\bibitem{Tropp}
J. Tropp, A.  Gilbert, and M. Straus.  Algorithms for simultaneous sparse approximation. part I: Greedy pursuit. {\em Journal of Signal Processing}, {\bf 86}, 572--588, 2006.

\bibitem{purity}
T. Tuytelaars, C. H. Lampert, M. B. Blaschko, and W. Buntine. Unsupervised object discovery:  a comparison. {\em IJCV}, 2009.

\bibitem{svm}
V. N. Vapnik. {\em The Nature of Statistical Learning Theory}. Springer, 2000.

\bibitem{PoT}
M. Welling,  G. E. Hinton, and S. Osindero. Learning sparse topographic representations with products of student-t distributions. {\em NIPS}, 2003.

\bibitem{Wu_active_basis_IJCV}
Y. N. Wu, Z. Si, H. Gong, and S. C. Zhu. Learning active basis model for object detection and recognition. {\em IJCV}, {\bf 90}, 198--235, 2010.

\bibitem{iFRAME}
J. Xie, W. Hu, S. C. Zhu, Y. N. Wu. Learning Inhomogeneous FRAME Models for Object Patterns. {\em CVPR}, 2014.

\bibitem{4}
M. Yang, L. Zhang, X. Feng, and D. Zhang. Fisher discrimination dictionary learning for sparse representation. {\em ICCV}, 2011. 

\bibitem{Younes}
L. Younes. On the convergence of Markovian stochastic algorithms with rapidly decreasing ergodicity rates. {\em Stochastics and Stochastic Reports}, {\bf 65}, 177--228, 1999.

\bibitem{Fergus}
M. Zeiler, G. Taylor, and R. Fergus. Adaptive deconvolutional networks for mid and high level feature learning. {\em ICCV}, 2011.

\bibitem{LeoZhu}
L. Zhu, C. Lin, H. Huang, Y. Chen, and A. Yuille. Unsupervised structure learning: hierarchical recursive composition, suspicious coincidence and competitive exclusion.  {\em ECCV}, 2008.

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

\bibitem{Minimax}
S. C. Zhu, Y. N. Wu,  and D. B. Mumford. Minimax entropy principle and its application to texture modeling. {\em Neural Computation}, {\bf 9}, 1627--1660, 1998.


\end{thebibliography}

\end{document}
% end of file template.tex

