Statistical Inference and Probabilistic Modeling in Compositional Vision by Wei Zhang B.A., Fudan University, P.R.China, 2003 M.S., Brown University, RI, 2006, M.A., Brown University, RI, 2007 A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in the Division of Applied Mathematics at Brown University PROVIDENCE, RHODE ISLAND May 2009 Abstract of “Statistical Inference and Probabilistic Modeling in Compositional Vi- sion” by Wei Zhang, Ph.D., Brown University, May 2009 This thesis is a mathematical and computational study of compositional vision. Three topics are covered: (1) ROC performance in a compositional world; (2) the construction of a probabilistic model for compositional structure; and (3) the con- struction of probabilistic model of image gray levels for a given vocabulary of ele- mentary parts. Chapter 1 introduces compositional vision and a probabilistic framework for model- ing hierarchy, reusability, and conditional data models. Chapter 2 focuses on theoretical questions about the ROC performance of various approaches to recognition in hypothetical compositional worlds. The results suggest that even sub-optimal decisions within a hierarchical framework will substantially outperform a decision process that does not explicitly allow for part-based decom- position. Chapter 3 focuses on the first component of the Bayesian approach to composi- tional vision: a prior probability model on hierarchical image interpretations. Non- Markovian (context-sensitive) distributions are investigated, and two theoretical questions are addressed. The existence of a class of non-Markovian distributions is established, and the convergence of an iterative perturbation scheme for achieving these distributions is proven. Chapter 4 focuses on the second component of the Bayesian approach to composi- tional vision: a probability model on pixel intensities conditioned on a given hier- archical structure. In particular, a generative approach to modeling object parts is developed through a probabilistic extension of the idea of fragment-based templates. Chapter 5 makes some conclusions and suggests future directions. c Copyright 2009 by Wei Zhang This dissertation by Wei Zhang is accepted in its present form by the Division of Applied Mathematics as satisfying the dissertation requirement for the degree of Doctor of Philosophy. Date Stuart Geman, Ph.D., Advisor Recommended to the Graduate Council Date Elli Bienenstock, Ph.D., Reader Date Basilis Gidas, Ph.D., Reader Approved by the Graduate Council Date Sheila Bonde, Dean of the Graduate School iii The Vita of Wei Zhang Wei Zhang was born in Baiquan, Heilongjiang province, China on November 15, 1981. She graduated from First High School of Baiquan. In 1999, she started her under- graduate study in Fudan University and received her Bachelor of Science degree in Mathematics four years later. She came to the United States in 2003 and has been attending the Ph.D. program in the Division of Applied Mathematics at Brown University. During the Ph.D. program, she received a Master of Science in Applied Mathematics in 2006 and a Master of Arts in Economics in 2007. She received the Stella Dafermos award in 2009 from the Division Applied Mathematics at Brown University. This dissertation was defended on April 13th, 2009. iv Acknowledgments I would like to thank all the people who have helped me to complete my thesis. Most of all, I want to extend my deepest gratitude to my advisor, Professor Stuart Geman, for his invaluable guidance, support, and encouragement throughout my graduate studies. His inspiring ideas and valuable suggestions shape my thought and guide my career. This thesis would not be possible without his persistent help. I want to thank my other committee members, Professor Elie Bienenstock and Pro- fessor Basilis Gidas, for reading my dissertation and giving me valuable comments. I want to acknowledge Eran Borenstein and Ya Jin, who have been very helpful during my research. Also I want to thank Sergey Kushnarev, Jingmei Qiu, Lo-Bin Chang, Christian Pfrang, Yanchun Wu, Yanqiu Li, and all the people in the Division of Applied Mathematics for helping me in every aspect. Last but certainly not least, I owe my thanks to my family. Thanks to my parents and brother, who support me spiritually. Thanks to my fiance, Wei Guo, for his love, support, and encouragement. Without their encouragement and understanding it would have been impossible for me to finish this work. v Contents Acknowledgments v 1 Introduction to Compositional Vision 1 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Markov Backbone . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.3.1 Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.3.2 Interpretations . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.3.3 Markov Probabilities . . . . . . . . . . . . . . . . . . . . . . . 9 1.4 Content Sensitivity and Non-Markovian Perturbations . . . . . . . . 11 1.5 Images and Image Probabilities . . . . . . . . . . . . . . . . . . . . . 14 1.5.1 Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.5.2 Independence Assumptions . . . . . . . . . . . . . . . . . . . . 15 1.5.3 A Note on Modeling Terminal Bricks . . . . . . . . . . . . . . 16 1.6 Scene Parsing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2 ROC Performance in a Compositional World 19 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.2 ROC Curves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.3 Comparison of Asymptotic ROC Curves of the Optimal Model, the Parts Model and the Universal-Null Model . . . . . . . . . . . . . . . 24 vi 2.3.1 Theoretical Comparison of Asymptotic ROC Curves . . . . . . 25 2.3.2 Demonstration of the Three ROC Curves for Finite n . . . . . 33 2.3.3 Generalization . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.3.4 Note . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.4 Comparison of ROC Curves of the Optimal Model, the Parts Model and the Universal-Null Model within a Hierarchical Setting . . . . . . 37 2.4.1 Theoretical Comparison of Asymptotic ROC Curves . . . . . . 39 2.4.2 Demonstration of the Three ROC Curves for Finite n Within the Hierarchical Setting . . . . . . . . . . . . . . . . . . . . . 44 2.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 2.5.1 Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3 On the Correctness of Compositional Probabilities 54 3.1 Existence of Probability Distribution Satisfying the Conditional Con- straints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.2 How to Achieve the Probability Distribution Satisfying the Condi- tional Constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4 Maximum-Likelihood Templates 77 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 4.2 Generative Probabilistic Model . . . . . . . . . . . . . . . . . . . . . 80 4.2.1 Basic Probabilistic Model with a Single Template . . . . . . . 80 4.2.2 Generalized Probabilistic Model with Multiple Templates. . . 85 4.2.3 Further Generalized Probabilistic Model with Multiple Scales, Rotations and Location Shifts. . . . . . . . . . . . . . . . . . . 90 4.3 Experiments on Learning Facial Part Templates and Applications on Ethnicity Classification . . . . . . . . . . . . . . . . . . . . . . . . . . 97 4.3.1 Facial Part Templates . . . . . . . . . . . . . . . . . . . . . . 97 4.3.2 Ethnicity Classification of Face Images . . . . . . . . . . . . . 103 vii 4.4 Background Template Learning . . . . . . . . . . . . . . . . . . . . . 106 4.4.1 A Basic Probabilistic Model of Background Image Patches . . 109 4.4.2 Experiment for the Basic Background Model . . . . . . . . . . 116 4.4.3 Improvement on the Background Model . . . . . . . . . . . . 117 4.4.4 Experiment for the Improved Background Model . . . . . . . 120 4.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 4.6 Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 5 Conclusion and Future Directions 130 viii List of Figures 1.1 Thought experiment. A simple compositional world. There are only four states to the world. Each state generates a binary random image of independent binary pixel intensities. ‘⊕’ in the figure stands for binary summation. Each “part” (a vertical or a horizontal bar) has only one pose, and the presence of both parts constitutes an ‘L’. . . 5 1.2 Semantic hierarchy for plate-reading application . . . . . . . . . . . . 7 1.3 Architecture. Left. A hierarchy of “bricks,” each representing a dis- junction of conjunctions. Bottom row is the image (pixel) data and the row above it is the set of terminal bricks. The state of a brick signals a chosen set of children. Right. An “interpretation,” which is an assignment of states such that the chosen children of any “on” brick are also on. There can be multiple roots and shared subtrees. Filled circles represent on bricks (non-zero states), and highlighted edges represent chosen children. . . . . . . . . . . . . . . . . . . . . 9 1.4 Samples from Markov backbone (left panel, ’4850’) and compositional distribution (right panel, ’8502’). . . . . . . . . . . . . . . . . . . . . 10 2.1 An example of two ROC curves. The ROC curve α is always strictly higher on the Y axis than the ROC curve β. Hence the decision rule associated with α is strictly better than the one associated with β. . 23 ix 2.2 Graphical model. Panel “A” is a graphical representation for the non-hierarchical setting in Section 2.3. Panel “B” is a graphical rep- resentation for the hierarchical settings in Section 2.4. . . . . . . . . 25 2.3 Comparison of empirical ROC curves of the three models. The upper left panel, the upper right panel, the lower left panel, and the lower right panel show the comparison results associated with n=3, 6, 9, and 12 respectively. The dotted lines represent the ROC curves produced by the optimal ROC curve given by the Neyman-Pearson Lemma; The solid lines represent the ROC curves produced by the parts model; The dashed lines represent the ROC curves produced by the universal-null model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.4 Comparison of empirical ROC curves of the three models within the hierarchical setting, for a non-trivial case. The upper left panel, the upper right panel, the lower left panel, and the lower right panel show the comparison results associated with n=2, 5, 10, and 50 respectively. The dotted lines represent the ROC curves produced by the optimal ROC curve given by the Neyman-Pearson Lemma; The solid lines represent the ROC curves produced by the parts model; The dashed lines represent the ROC curves produced by the universal-null model. The six points {(FPRk , TPRk )}5k=0 were plotted as six star points (*) on the X-Y plane in each panel. . . . . . . . . . . . . . . . . . . . . 48 x 2.5 Comparison of empirical ROC curves of the three models within the hierarchical setting, for a trivial case. The upper left panel, the upper right panel, the lower left panel, and the lower right panel show the comparison results associated with n=2, 5, 10, and 50 respectively. The dotted lines represent the ROC curves produced by the optimal ROC curve given by the Neyman-Pearson Lemma; The solid lines represent the ROC curves produced by the parts model; The dashed lines represent the ROC curves produced by the universal-null model. The three star points (*) in each panel correspond to points (0,0), (P (X α > 0, X β > 0|X γ = 0) = 0.25,1), and (1,1). . . . . . . . . . . . 51 2.6 Comparison of empirical ROC curves of the three models within the hierarchical setting, for a non-trivial case. The upper left panel, the upper right panel, the lower left panel, and the lower right panel show the comparison results associated with n=12, 15, 20, and 25 respectively. The dotted lines represent the ROC curves produced by the optimal ROC curve given by the Neyman-Pearson Lemma; The solid lines represent the ROC curves produced by the parts model; The dashed lines represent the ROC curves produced by the universal-null model. The three star points (*) in each panel correspond to points (0,0), (P (X α > 0, X β > 0|X γ = 0) = 0.25,1), and (1,1). . . . . . . . 52 3.1 Architecture. Left. A hierarchy of “bricks,” each representing a dis- junction of conjunctions. Bottom row is the image (pixel) data and the row above it is the set of terminal bricks. The state of a brick signals a chosen set of children. Right. An “interpretation,” which is an assignment of states such that the chosen children of any “on” brick are also on. There can be multiple roots and shared subtrees. Filled circles represent on bricks (non-zero states), and highlighted edges represent chosen children. . . . . . . . . . . . . . . . . . . . . 55 xi 3.2 The special hierarchy structure in Step 1 and Step 2, when N=4. . . . 58 4.1 12 face images from Feret Face database, each with 17 landmarks labeled manually. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 4.2 The left panel shows 70 training left eye images. The right panel shows the 16 learned templates. . . . . . . . . . . . . . . . . . . . . . 99 4.3 The evolution of the first 8 templates as the EM algorithm ran. . . . 100 4.4 The evolution of the last 8 templates as the EM algorithm ran. . . . . 100 4.5 The left panel shows 70 training left eye images, each with size 15×23). The right panel shows the 16 learned templates, each with size 12×19, under the model with mixtures over multiple templates and spatial shifts. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 4.6 The left panel shows the 16 learned templates from the model not con- sidering spatial shifts. The right panel shows the 16 learned templates from the model considering spatial shifts templates. . . . . . . . . . . 102 4.7 The left panel shows 70 training left eye images, with different sizes ranging from 12by18 to 18by27; eyes are not in the center. The right panel shows the 16 learned templates, each with size 12 by 19, under the model considering spatial shifts and mixtures over two scales 0.9 and 1.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 4.8 The evolution of the 16 templates during 8 runs of the EM algorithm, starting from a random initialization. . . . . . . . . . . . . . . . . . . 104 4.9 The left panel shows 120 training image patches. The right panel shows the 16 learned templates, each with size 15 × 18, from the fully generalized model, with mixtures over multiple templates, spatial shifts, scales, and rotations. . . . . . . . . . . . . . . . . . . . . . . . 105 4.10 The evolution of the 16 templates during 8 runs of the EM algorithm, starting from a random initialization. . . . . . . . . . . . . . . . . . . 106 xii 4.11 The left panel shows 70 training 13× 18 nose or mouth image patches. The right panel shows the 16 learned templates each with size 11 × 16. 107 4.12 The evolution of the first 16 templates as the EM algorithm ran, starting from a random initialization. . . . . . . . . . . . . . . . . . . 108 4.13 The evolution of the last 16 templates as the EM algorithm ran, start- ing from a random initialization. . . . . . . . . . . . . . . . . . . . . . 109 4.14 The left panel shows 70 East Asian eyes, each with height of 10 pixels. The right panel shows the 8 templates, each with size 8 × 20, learned from the model with mixtures over 4 scales: 1.2, 1.1, 1, and 0.9. . . . 110 4.15 The left panel shows 70 Indian eyes, each with height of 10 pixels. The right panel shows the 8 templates, each with size 8 × 20, learned from the model with mixtures over 4 scales: 1.2, 1.1, 1, and 0.9. . . . 111 4.16 59 natural images collected from Internet. . . . . . . . . . . . . . . . 117 4.17 The left panel shows 500 randomly selected training image patches (15 × 15). The right panel shows the 32 learned templates, each with size 10 × 10, learned from the basic background model. . . . . . . . . 118 4.18 The evolution of the 32 templates during 20 runs of the EM algorithm, under the basic background model, where the first row shows the random initialization of the templates. . . . . . . . . . . . . . . . . . 119 4.19 The left panel shows 500 randomly selected training image patches (15 × 15). The middle panel shows the 32 learned templates from the improved background model, each with size 10 × 10. The right panel shows the 32 learned templates from the basic background model, each with size 10 × 10. . . . . . . . . . . . . . . . . . . . . . . . . . . 121 4.20 The evolution of the 32 templates during 20 runs of the EM algorithm under the improved background model, where the first row shows the random initialization of the templates. . . . . . . . . . . . . . . . . . 122 xiii Chapter 1 Introduction to Compositional Vision 1.1 Introduction Most approaches to object recognition or detection in computer vision can be broadly characterized as either generative or discriminative according to whether or not a probability distribution of the image features is modeled. Discriminative approaches are popular in classification related tasks. The idea is to compute a direct mapping from an observed variable Y to a hidden variable X for classification; no direct attempt is made to model the underlying distributions. Various forms of penalty functions, regularization, and kernel functions are used to prevent overfitting. Many discriminative models have been proposed and widely used in computer vision, involving, for example, support vector machines [1], neural networks [2] and boosting [3]. Generative approaches specify a joint probability distribution P (X, Y ), by spec- ifying a prior distribution P (X) and a likelihood function P (Y |X). Once the gener- ative model is built, classifiers can be derived in a straightforward way by exploiting the posterior distribution P (X|Y ), calculated through Bayes’ formula. In contrast 1 2 with generative models, discriminative models focus only on the posterior distribu- tion P (X|Y ). Hence, a generative model features being able to simulate (i.e., to generate, or to sample) Y , but a discriminative model does not. Gaussian mix- ture models, hidden Markov models, and Markov random fields are commonly used generative models. In the recent years, hierarchical generative models have gained more attention due to their rich representation and search efficiency. And a growing number of them are achieving state-of-art performance in a growing number of ap- plications. Some are biologically motivated (e.g. [34]), others are computationally motivated (e.g. [4]). Some involve learned hierarchies (e.g. [9]), others are hand- designed (e.g. [16, 15]). Some are deep hierarchies (e.g. [21, 22] and [23]), some are medium ([29]), others are shallow (e.g. POP model [19, 39], Constellation model [20, 36]). This thesis work falls into the category of hierarchical generative model. It is motivated by both biological evidence in human vision and theoretical evidence for improved ROC performance in compositional (hierarchical) models. 1.2 Motivation • Primates, especially humans, are remarkably good at learning, recognizing, and generalizing objects from a few examples. The performance of their vi- sual system and its robustness surpasses the best state-of-art computer vision systems. Over the last decade, evidence has accumulated about some key fea- tures of the ventral visual pathway. The human vision system depends on a hierarchy of successive layers in the visual cortex. The first layers of the cortex detect an object’s simpler features, such as oriented edges, and higher layers integrate that information to form our perception of the object as a whole. It is a simple-to-complex cell hierarchy, [33, 34]. And there is an apparent increase in both invariance and selectivity in moving from the primary visual cortex 3 to the infero-temporal cortex (see, e.g. [35] and references therein). In addi- tion to the biological evidence from the ventral visual pathway, it is evident that humans have the tendency to represent entities as hierarchies of reusable parts. In a visual world, objects and scenes naturally decompose into hierar- chies of meaningful and generic parts. Natural languages are also hierarchical: hierarchies of words, phrases, and sentences whose composition is governed by grammatical rules. • In a miniature compositional world, better ROC performance of compositional model has been justified theoretically, compared to the models that do not accommodate compositionality. These issues will be introduced and studied in detail in Chapter 2. The basic idea is as follows. Consider a very simple and perfectly composi- tional world. There are two parts—vertical and horizontal bars. Each part, when it appears, appears in a fixed pose. When both parts appear, we declare the presence of an ‘L’. The world, then, has four states: H0 which generates no parts; H1 which generates only a vertical bar; H2 which generates only a horizontal bar; and H3 which generates an L. Figure 1.1 shows the four hy- potheses in this compositional world. For each state, a binary scene (image) is generated by adding independent Bernoulli noise to the binary clean image (‘1’ for foreground, ‘0’ for background), where the addition is binary summation. The problem is to build a classifier that recognizes L’s. According to the Neyman-Pearson lemma, the optimal classifier would base the decision (L or not L) on the likelihood ratio: If Y is the observed array of pixel intensities (the image) then declare L if P (Y |H3 ) ¯ 3) ≥ c P (Y |H (1.1) ¯ 3 is short for the compound (mixture) event and not L otherwise. Here, H 4 H0 ∪ H1 ∪ H2 , so ¯ 3 ) = P0 · P (Y |H0 ) + P1 · P (Y |H1 ) + P2 · P (Y |H2 ) P (Y |H The threshold c, over the range [0, ∞), sweeps out the optimal ROC curve. Specifically, for each value of c there is a detection probability (i.e. true pos- itive rate) pd and a false alarm probability (i.e. false positive rate) pf , and among all classifiers with false alarm probability pf , pd is the highest possible detection probability. The problem with the Neyman-Pearson prescription is that it is impractical. Not actually impractical for the thought experiment, but impractical in anything resembling a real experiment, with multiple parts, multiple objects, variable poses, and so on. An expedient alternative is to devise a “universal-null” model, a serviceable probability on Y under the “not object” condition. The “white noise” model is analyzed here as a universal null (i.e. P (Y |H0 )), and there would be a direct generalization to alternative models for P (Y |H0 ). In particular, concerning the thought experiment, we will consider the performance of the likelihood ratio ¯ 3 ) is replaced by the background distribution test (equation (1.1)) when P (Y |H P (Y |H0 ): declare L if P (Y |H3 ) ≥c (1.2) P (Y |H0 ) and not L otherwise. Another alternative is to accommodate the fact that a ‘L’ is composed of two parts, and to declare L only when there was sufficient evidence for both the vertical and horizontal bars. In and of themselves, these simple building blocks are easy to devise tests for, and in fact in the artificial world of the thought experiment the optimal tests involve no mixtures, as the alternative in each case (the denominator) really is the background, white-noise, model. Consider, 5 then, the “parts” strategy that declares L if P (Y |H1 ) P (Y |H2 ) ≥ c and ≥c (1.3) P (Y |H0 ) P (Y |H0 ) (or, equivalently, if min( PP (Y |H1 ) P (Y |H2 ) , (Y |H0 ) P (Y |H0 ) ) ≥ c), and not L otherwise. Within this problem framework, two theorems – one under a non-hierarchical setting and the other under a hierarchical settings – have been established in Chapter 2 to compare the ROC performance of these three strategies, as the image resolution (i.e. the number of pixels representing a given area of the image) goes to infinity. The ROC curve produced by the parts strategy is exponentially better than the one produced by the universal-null strategy, and is comparable and eventually merges together with the one produced by the optimal strategy given by the Neyman-Pearson lemma. Figure 1.1: Thought experiment. A simple compositional world. There are only four states to the world. Each state generates a binary random image of independent binary pixel intensities. ‘⊕’ in the figure stands for binary summation. Each “part” (a vertical or a horizontal bar) has only one pose, and the presence of both parts constitutes an ‘L’. If the world is indeed compositional, we would do better to take advantage of this fact in devising machines intended to mimic or even just approach human per- 6 formance. The thought experiment of detecting ’L’ above gives us a hint on how to accommodate the structure feature of our world to achieve better performance in computer vision tasks. Certainly, our vision world is way more complicated than the miniature compositional world in the thought experiment. But the philosophy is similar, combining semantic parts into semantic objects. For example, it is obviously crucial to consider the relative coordinates of the parts before combining parts into higher level semantic objects. For example, ‘L’ and ‘T’ are both composed of a hori- zontal bar and a vertical bar, but the relative location of these two parts determines whether to combine them into ‘L’, or ‘T’, or neither. This thesis work is built on the platform of the“compositional machine” – a generative probabilistic model on hierarchies of reusable parts under the Bayesian framework, and has been pioneered by Geman, Bienenstock, and their colleagues [11, 14]. It has been further developed in the Ph.D. work of Huang [12] (preliminary computational experiments), Harrison [13] (some learning related work), and Jin [16, 15] (implementation in license-plate reading). The compositional machine has two components, the prior distribution on the (image) interpretations, and the conditional data distribution (i.e. the likelihood function) on the image given the interpretations. Section 1.3 through section 1.6 will introduce the framework of compositional machine as follows: Section 1.3 will introduce the Markov distribution on the interpretations. We call this the “Markov Backbone”, which serves as a reference distribution. Markov systems like branching processes and probabilistic context-free grammars qualify, but are generally too weak to capture context and content-dependent likelihoods. The Markov property can be saved, but only at the cost of a very large state space. A better approach might be to look for workable non-Markovian distributions. Section 1.4 will extend the “Markov Backbone” to a “compositional system” by introducing a non-Markovian term through a “perturbation” argument. Section 1.5 will describe the conditional data distribution. Section 1.6 will discuss scene parsing via the posterior distribution 7 and Bayesian inference. 1.3 Markov Backbone The application that we have in mind is to Bayesian scene analysis through a prior distribution on scene “parses” (interpretations). Parses are represented in a graphical model. The components of a parse are low-to-high-level abstract variables such as “edge,” “eye,” “face,” “person,” “people,” or “crowd”. To emphasize the reusability of these “parts”, the vertices of the graph are called bricks (as in Lego bricks). The specific assignment is application dependent. For example, in the application of reading license plates [15, 16], the semantic bricks represent different meanings as shown in Figure 1.2. Figure 1.2: Semantic hierarchy for plate-reading application 8 1.3.1 Notation B finite set of bricks T ⊆B terminal bricks xα ∈ {0, 1, . . . , nα }, nα ≥ 1 states of the α brick, α ∈ B ~x = {xβ : β ∈ B} state of all bricks α Pnα {αi}ni=0 , 0 ≤ αi ≤ 1, α i=0 i =1 state probabilities, α ∈ B Ciα ⊆ B, α ∈ B \ T i’th set of children of α, 1 ≤ i ≤ nα , (Ciα 6= Cjα when i 6= j) 1.3.2 Interpretations Consider a directed acyclic graph (DAG) G defined by • A vertex for every brick β ∈ B • A directed edge from α to β if β ∈ Ciα for some i ∈ {1, 2, . . . , nα } An “interpretation” ~x is defined as an assignment of states to {xβ }β∈B such that α ∈ B \ T and xα > 0 ⇒ xβ > 0 ∀β ∈ Cxαα . Let I be the set of interpretations. If we declare a brick α “on” when xα > 0, and if we call Cxαα the chosen children of brick α in state xα > 0, then an interpretation is a state vector ~x in which the chosen children of every non-terminal on brick are themselves on. 9 Figure 1.3: Architecture. Left. A hierarchy of “bricks,” each representing a disjunc- tion of conjunctions. Bottom row is the image (pixel) data and the row above it is the set of terminal bricks. The state of a brick signals a chosen set of children. Right. An “interpretation,” which is an assignment of states such that the chosen children of any “on” brick are also on. There can be multiple roots and shared subtrees. Filled circles represent on bricks (non-zero states), and highlighted edges represent chosen children. 1.3.3 Markov Probabilities For ~x ∈ I, we define the below set B = B(~x) by B = {β ∈ B : β ∈ Cxαα , for some α ∈ B \ T with xα > 0} The Markov (“context-free”) probability of an interpretation ~x ∈ I is defined as Q β β∈B (xβ ) P (~x) = Q β (1.4) x) (1 − 0 ) β∈B(~ Remarks: P 1. ~ ∈I x P (~x) = 1, as can be seen by ordering G by generations, starting with the roots, and then generating a random ~x in the same order, according to α αi, i ∈ {0, 1, . . . , nα }, for any brick not chosen by a parent, and i 1−α , i ∈ 0 {1, 2, . . . , nα}, otherwise. 10 2. P (~x) is a ‘Bayes Net’ with respect to the DAG G, and hence Markov with respect to the undirected ‘moral’ graph derived from G. 3. There is an obvious connection to probabilistic context-free grammars: think of α α → {β : β ∈ Ciα} as a production, chosen with probability i 1−α . But keep in 0 mind that there is no unique “start” symbol, that an interpretation can include many trees, that trees can share parts (instantiations overlap), and that there is a fixed topology (hence no recursion). In the license application [15, 16], Jin and Geman sampled from the Markov backbone, given the semantic assignment of bricks (in Figure 1.2) and the manually hardwired children sets. The left panel of Figure 1.4 shows a 4-digit sample under the Markov backbone. As seen from the figure, although the parts of each digit are present and in roughly the correct locations, neither the parts nor the digits are properly situated. Figure 1.4: Samples from Markov backbone (left panel, ’4850’) and compositional distribution (right panel, ’8502’). 11 1.4 Content Sensitivity and Non-Markovian Per- turbations Most of the proposed generative models in the literature share the Markov property ([20], [27], [17], [7], [47, 48]), due to its computational advantage. But this context- free (Markov) property is problematic. Constituents, in vision and language, are composed with a likelihood that depends not just on their “labels,” (stroke, let- ter, noun phrase, verb phrase, etc.), but also on the details of their instantiations (position, font, gender, tense, etc.). Biological-level ROC performance of an image- analysis system will almost certainly need to be content sensitive. This raises the difficult question of constructing useful non-Markovian probability distributions on hierarchical models. One approach, beginning with coding and description length, was explored in [14]. A different approach, through perturbations, is explored here. Imagine that we have, associated with every brick β ∈ B, an attribute function (scalar valued or vector valued), aβ (~x). A prototypical example is the set of pose coordinates (or relational pose coordinates) of the chosen children of β. Depending on the depth of the instantiation of the children, aβ (~x) may depend on the states of bricks that are several generations removed from β itself (grandchildren, great grandchildren, etc.). Start with the Markov probability P , as defined in (1.4), and fix a particular brick γ ∈ B. Under (1.4), aγ has some distribution, P0γ (aγ |xγ ), for every state, xγ ∈ {0, 1, . . . , nγ }. If, say, aγ (~x) is the vector of poses of the chosen children of γ, then there is no reasonable hope that P0γ corresponds to the empirical (or “real- world”) distribution on the positions of the parts of γ. After all, (1.4) is context free and, in particular, the instantiations of the chosen children of γ are independent (Markov property). Let Pcγ (aγ |xγ ) (as opposed to P0γ (aγ |xγ )) be the correct conditional distribution on the attribute aγ . One way to “perturb” P in (1.4), so as to correct the conditional 12 aγ distributions, is to choose the new distribution, call it P ∗ , that is closest to P , subject to the constraint that P ∗ (aγ |xγ ) = Pcγ (aγ |xγ ) for all aγ and all xγ . If, by “closer,” we mean that D(P ||P ∗ ) (Kullback-Leibler divergence) is minimized, then it is an easy calculation to show that Pcγ (aγ (~x)|xγ ) P ∗ (~x) = P (~x) P0γ (aγ (~x)|xγ ) Remarks: 1. The particular distribution, Pcγ (aγ |xγ = 0), is largely irrelevant to the problem of modeling an object γ, would be very hard to measure, and in any case can be taken as P0γ (aγ |xγ = 0) so that there is no perturbation at all unless the γ brick is on. 2. Bearing in mind the considerations of the previous remark, P ∗ is a “perturba- tion” in the sense that P is only altered in the event of xγ > 0 (γ “on”), which is presumably quite rare for most bricks, as they represent particular parts, objects, or collections of objects. 3. In general P ∗ is no longer Markov, but it is still normalized. 4. In most cases of interest, aγ (~x) would be a function only of xγ and its possible progeny, meaning every brick that could appear in its instantiations. Evidently, the process can be repeated, at other bricks, enforcing a brick-dependent attribute distribution at each step. For any “brick visitation schedule,” γ1 , γ2 , . . . , γ|B| , with {γ1 , γ2 , . . . , γ|B| } = B, we end up with a distribution |B| ∗ Y P γv (aγv (~x)|xγv ) c P (~x) = P (~x) v=1 P˜ γv (aγv (~x)|xγv ) Q β |B| β∈B (xβ ) P γv (aγv (~x)|xγv ) Y c = Q , (1.5) x) (1 − β∈B(~ β0 ) v=1 P˜ γv (aγv (~x)|xγv ) 13 where P˜ γv (aγv (~x)|xγv ) is the distribution on aγv given xγv at the time of the visit to the γv brick, and P˜ γv (aγv (~x)|xγv ) = P0γv (aγv (~x)|xγv ) when v = 1. The result is unsatisfactory in two regards: 1. The distribution turns out to be different for different visitation schedules. 2. Each perturbation, while establishing a desired conditional distribution Pcγ (aγ |xγ ), perturbs the previously established distributions, so that the already-visited bricks no longer have, precisely, the desired attribute distributions. (This ap- plies to the epsilon probabilities as well.) The study of specific examples suggests that the attribute functions {aγ (~x)}γ∈B together with the attribute (conditional) distributions {Pcγ (aγ |xγ )}γ∈B will in general under-determine the distribution on I: there are typically many distributions with the desired constraints. But we do not yet have a satisfactory theory guaranteeing the existence of such a distribution. See Chapter 3 for a guarantee under some restrictions on the conditional constraints. There is a related question of convergence: is it true that an iterative procedure, that visits every site infinitely often, converges to a distribution with the desired attribute probabilities? This turns out to be true under quite general conditions (Chapter 3). The license-plate application explored in [15, 16] used a simple approximation. Each pre-perturbation probability, P˜ γ (aγ |xγ )) in (1.5), was assumed to be close to, and was replaced by, the corresponding conditional probability under the Markov distribution (1.4), which is denoted by pγo (aγ |xγ ). The right panel of Figure 1.4 shows a compositional 4-digit sample generated by Jin and Geman [15, 16] from this non-Markovian model. As we can see, dramatic improvement is achieved, com- pared to the sampling result from the Markov backbone (in the left panel of Figure 1.4). Although the dynamic programming machinery is no longer available for non- Markovian models, certain coarse-to-fine computational engines are available. The Markov distribution is easy to work with and estimates (even exact values) 14 can be derived for the conditional attribute probabilities. Since the target distribu- tions, {Pcγ (aγ |xγ )}γ∈B , are fixed (by hand or inference) and the “null” probabilities {pγo (aγ |xγ )}γ∈B , all derive from the Markov distribution (1.4), there is no depen- dence on order. These considerations lead to the useful (and order-independent) approximation: Q β β∈B (xβ ) Y P β (aβ (~x)|xβ ) ∗ P (~x) ∝ Q x) (1 β∈B(~ − β0 ) β∈B pβo (aβ (~x)|xβ ) A price is paid in that the normalization is no longer exact. On the other hand, the parameters in the Markov “backbone” (1.4), as well as the null probabilities under the Markov distribution, can be estimated by more-or-less standard approaches, and the remaining terms, the brick-conditioned attribute probabilities, are in principle available from examples of the objects of interest. 1.5 Images and Image Probabilities The Bayesian (generative) framework is completed by specifying a “data model”: a probability distribution on images given an interpretation, ~x ∈ I. 15 1.5.1 Notation R index set (pixels) of the “image” Y = {yj : j ∈ R} image (pixel grey levels) YD = {yj : j ∈ D} image values at locations j ∈ D, for any D ⊆ R Rτi ⊆ R, τ ∈ T image locations in the support of terminal brick τ ∈ T when xτ = i > 0 Snτ i=1 Rτi “receptive field” of brick τ ∈ T Given an interpretation ~x ∈ I, define D = D(~x) = {τ : xτ > 0}. The support of an interpretation ~x ∈ I is defined as [ RD = RD (~x) = Rτxτ τ ∈T xτ >0 The support is the set of pixels directly addressed by an interpretation. 1.5.2 Independence Assumptions These are assumptions about the conditional distribution on pixel intensities given an interpretation. They are not unreasonable, as approximations, and they make data modeling much easier. Use xT to indicate the configuration of the terminal bricks, {xτ : τ ∈ T }. A1. P (Y |~x) = P (Y |xT ) the conditional distribution on image data depends only on the states of the terminal bricks 16 Let ~x0 ∈ I be the “zero” interpretation: ~x0 = {xβ0 }β∈B where xβ0 = 0 ∀β ∈ B. P (Y |xT ) P (YRD |xT ) A2. P (Y |xT = P (YRD |xT the (data) likelihood ratio of interpretation 0 ) 0 ) ~x to the “zero” interpretation, ~x0, depends only on the data in the support of ~x Remark: A2 holds if, for example, the image data that is not supported is i.i.d. from a fixed “null” distribution. 1.5.3 A Note on Modeling Terminal Bricks The support of an interpretation ~x is covered by the supports of the active (“on”) terminal bricks. These define disjoint connected components of pixels (connected by overlapping supports), and if the independence assumptions are expanded to connected components, then the task of data modeling is the task of data modeling a set of overlapping supports, conditioned on states of the corresponding terminal bricks. These models can be built from individual templates – one for each support Rτi , i ∈ {1, 2, . . . , nτ }. What is a reasonable model for yRτi , given that xτ = i? Here is an outline of a model that is learnable and has performed well in experiments. 4 Fix a terminal brick τ and a state i (xτ = i). For ease of notation, let G = Rτi . The task is to model YG , under the condition that xτ = i and that there is no other (active) support intersecting G. (The extension to overlapping supports follows similar reasoning.) The idea is that we want to model a distribution P (YG ) in terms of something resembling a “sufficient statistic” S(YG ): P (YG ) = P (S(YG ))/c(S(YG )) 17 where c(s) = #{Y˜G : S(Y˜G ) = s} The probability on the data depends only on the statistic (“feature”) S, and hence all configurations of YG that give the same value of S are equally likely. Notice that, whereas P (S) is the marginal distribution on S (conditioned on xτ = i), it would be a mistake to model the the probability of the actual data, YG , as being proportional to P (S(YG )). The combinatorial factor, c(s), can substantially tilt the data likelihood in favor of configurations that have low “redundancy” (small combinatorial factors, or what are sometimes miss-labeled as low entropies). A simple, and effective, concrete example is S(YG ) = Corr(T, YG ) where T = {tk }k∈G is a template, and Corr(T, YG ) is the normalized correlation between T and YG (and hence between −1 and +1). A good fit to real data (say patches of left eyes, or portions of characters or strokes) can be made by taking P (S = s) to be monotone increasing on s ∈ [−1, 1], as in the backwards exponential P (S = s) ∝ e−λ(1−s) ≈ λe−λ(1−s) . What is needed is an approximation to the combinatorial factor, c(s), since the direct calculation is intractable. Chapter 4 studies image probabilities and terminal bricks in detail. It gives a derivation for an approximation for c(s), and an explicit likelihood model. The parameters of the model (the template T and the scalar λ) can be learned from data through the maximum likelihood equations. In fact, templates can be learned despite pose (e.g. scale, rotation, and translation) variation in the training data, and can be learned at super or sub resolution. Chapter 4 also includes recognition experiments with maximum-likelihood templates. 18 1.6 Scene Parsing In the previous sections, we have specified the prior distribution P (~x) and the data model P (Y |~x). At this point, the generative model has been fully built up. Equipped with this generative model, the computer vision task of scene parsing (i.e. the image interpretation) can be approached through the posterior distribution on an image interpretation ~x given the input image data Y , P (~x|Y ): P (Y |~x)P (~x) P (~x|Y ) = ∝ P (Y |~x)P (~x). (1.6) P (Y ) The image interpretation ~x provides a rich semantic and syntactic parsing of the scene through the brick variables of the compositional machine. With ~x, each image pixel can be identified as either plain background (with no structure), or cluttered background (with structure), or objects at different complexity levels. For example, the license demonstration system [15, 16] could read out from interpretation ~x the detailed locations and identifications of license plate, string, boundary, partial string, L-junction, character, line, and parts, etc., which may occur anywhere across the image. Ideally, the MAP estimator could be computed exactly given the prior model P (~x) and data model P (Y |~x). But the sample space is very high dimensional con- sidering that a compositional machine can easily contain thousands of bricks. Hence, exact inference is intractable and some computationally feasible approximations are needed. The computation engine explored in [15, 16] took advantage of the coarse- to-fine search strategy, motivated by [4]. Jin and Geman [15, 16] focused on a bottom-up pass and touched on a general depth-first search strategy, aiming at fast object detection, for example, early detection of a license boundary. Chapter 2 ROC Performance in a Compositional World 2.1 Introduction Human vision is both highly invariant (invariance refers to the extent to which objects are detected independent of their pose and rendering) and highly selective (selectivity refers to avoiding the misclassification of other structures). Computer vision systems are rarely both. In particular, systems which are highly invariant often suffer from an unacceptable number of false detections. Such an “invariance and selectivity” dilemma is due to a combination of extreme variation in object presentations and the high structure of background clutter. When facing the problem of testing for the presence of an object O in an input image I, we can not do better than to base our decision on the likelihood ratio, i.e., to threshold the ratio P (I|O)/P (I|O), where O stands for “object present” and O stands for “object absent.” This is the Neyman-Pearson Lemma. We refer to the decision model suggested by this lemma as the “optimal model” in this chapter. The problem with Neyman-Pearson Lemma prescription is that it is not generally practical to enumerate O. It is one thing to have a model (a likelihood function) of 19 20 the data given the presence of a particular object of interest (as in P (I|O)), but it is quite another thing to have a model of the data given that the object is not present (as in P (I|O)). “Not present” is a mixture—a very big mixture—and it would be bad news indeed if good classification or object detection required a good data model under this mixture. One way around this is to approximate P (I|O) by P (I), since the set O is very nearly the entire image set. Then P (I) is usually referred as the background or natural-scene model. But of course P (I) is complicated to calculate as well. An expedient alternative is to devise a “universal-null model,” a serviceable prob- ability on I under the “object absent” condition. We will analyze here the “white noise” model as a universal null. Given that the world is highly structured, this naive assumption does not seem proper, although it is widely used (often implicitly) in recognition algorithms. For example, it is common to see in the literature that P (I|O) is thresholded for object O detection. But thresholding P (I|O) is the same as thresholding P (I|O)/P (I|O) while assuming the denominator to be a constant, for example (1/256)|I| , as in the i.i.d. white noise background model. We argue that we can do better than the universal-null model, at least if we accept that the principles of hierarchy and reusability are operating in the real world. Since the clutter of background shares the same reusable parts as the target (object), the idea of compositionality suggests that we devise tests to address the inaccuracy caused by an unnecessary binary choice between target and universal null. The idea is to perform multiple hypothesis tests on object parts, i.e., to detect the parts of the object O, and only claim the existence of O if all the part tests succeed. We call this model the “parts model.” In this chapter, we demonstrate the efficient discrimination available through compositional and hierarchical representation via two theorems establishing a better receiver operating characteristic (ROC) curve. We study the asymptotic (as the res- olution of the data goes to infinity) ROC performance of three models: the optimal 21 model (given by the Neyman-Person Lemma), the parts model, and the universal- null model. This chapter is organized as follows: Section 2.2 will introduce some background materials about ROC curves. Section 2.3 will employ the Large Devia- tion Theory to prove that the parts model is comparable to the optimal model, and is strictly better than the universal-null model (in terms of asymptotic ROC per- formance). Section 2.4 will compare the three models within a hierarchical problem setting. Finally, Section 2.5 will conclude with a discussion and a summary. 2.2 ROC Curves In signal detection theory, a ROC curve is a graphical plot of the true positive rate (TPR) versus the false positive rate( FPR) for a binary classification system as its discrimination threshold varies. ROC analysis is an efficient tool for researchers to select possibly optimal models and to discard suboptimal ones. It is widely used in multiple-disciplines – for example, decision making in medical diagnosis, model selecting in biosciences, machine learning and data mining, etc.. We now introduce some basic concepts in ROC analysis that will be used in this Chapter. Let Y be the observed data. The goal is to make a binary decision between two hypotheses Y ∈ H1 or Y ∈ H0 = H1c . The discrimination rule is Y ∈ H1 if F (Y ) ≥ c and Y ∈ H0 otherwise, where F is a determinant function designed by the user, mapping Y to a real number, and c is a threshold, c ∈ R. Theoretically, if H1 is associated with a probability distribution P1 on domain Ω, and H0 is associated with a distribution P2 on domain Ω, we can calculate the theoretical true positive rate and false positive rate as TPR = P1 ({Y ∈ Ω : F (Y ) ≥ c}), FPR = P2 ({Y ∈ Ω : F (Y ) ≥ c}). In practice, we can only observe samples from H1 and H0 . The observed dataset 22 S can be divided into positive dataset O and negative dataset B, where O ∩ B = ∅. O contains data Y sampled from H1 , while B contains data sampled from H0 . Let D = {Y ∈ S : F (Y ) ≥ c}; then TPR and FPR are calculated as follows: |D ∩ O| |D ∩ B| TPR = , FPR = , |O| |B| where | · | is the measure of a set. For both the theoretical case and the practical case, if we plot true positive rate versus false positive rate in the X-Y plane, ∀c ∈ (−∞, +∞), we get the ROC curve. Each single threshold c corresponds to a point on the ROC curve. The best possible prediction result corresponds to the point (0,1) in the X-Y plane, which stands for 100 percent detection rate and 0 percent false positive rate. The (0,1) point is also called a perfect classification. A random classification, say by flipping a biased coin, would give a point along the diagonal line from the left bottom to the top right corners (the so-called line of no-discrimination). In general, the more closely the ROC curve approaches the (0,1) point, the better the classification result. If one ROC curve is always strictly higher on the Y axis, then we refer to this case as “the first ROC curve is strictly better than the other.” See the example in Figure 2.1. In statistics, the Neyman-Pearson Lemma states that when performing a hy- pothesis test between two point hypotheses H0 : θ = θ0 and H1 : θ = θ1, among the likelihood-ratio test which rejects H0 in favor of H1 , the following testing rule is optimal: P (Y |θ1) F (Y ) = ≥ c, ∀ c ∈ R. P (Y |θ0) i.e., for any other test Fˆ with equal or smaller false positive rate than F , its true positive rate must be no greater than the one given by F . In other words, the ROC curve given by F is at least as good as the one given by Fˆ . Another concept that we will need later is the area under the ROC curve, denoted by AUROCC. The discussion and notation are simplified by assuming that Y is a 23 Figure 2.1: An example of two ROC curves. The ROC curve α is always strictly higher on the Y axis than the ROC curve β. Hence the decision rule associated with α is strictly better than the one associated with β. continuous random variable under both hypotheses H0 and H1 . AUROCC is a good measurement to assess the classifier, especially when evaluating two ROC curves which cross each other. AUROCC reaches its maximum 1, when the classifier makes error-free decisions. Green and Swets gave a concise formula to evaluate AUROCC theoretically, [24]: AUROCC = P ({(Y1 , Y2 ) ∈ H1 × H0 : F (Y1 ) ≥ F (Y2)}). Accordingly, the area above the ROC curve, denoted as AAROCC, is equal to one minus AUROCC, i.e., AAROCC = P ({(Y1 , Y2 ) ∈ H1 × H0 : F (Y1) ≤ F (Y2)}). (2.1) 24 2.3 Comparison of Asymptotic ROC Curves of the Optimal Model, the Parts Model and the Universal-Null Model Consider a thought experiment. Imagine a target (object) composed of two parts. The generative model is Yn = fn (X) + Zn , where X stands for the presence of one target with two parts. X takes values in the set {(0, 0), (1, 0), (0, 1), (1, 1)} (standing for four cases: nothing, only part 1 present, only part 2 present, the whole target present) with probability p0 , p1 , p2 , p3 respectively (p0 + p1 + p2 + p3 = 1, pi > 0, ∀i). See Figure 2.2A for a graphical representation. fn (X) is obtained by repeating the left part of X n times, followed by repeating the right part of X n times. For example, if n = 3 and X = (1, 0), we will have fn (X) = (1, 1, 1, 0, 0, 0). Zn represents random noise. Zn is a vector with the same length as fn (X), containing 2n i.i.d Bernoulli() random variables, 1 <  < 0.5. So Zn = (z1, . . . , z2n ), and P (zi = 1) = , P (zi = 0) = 1 − . The addition between fn (X) and Zn follows the rule of binary summation. For example, if n = 3, fn (X) = (1, 1, 1, 0, 0, 0), Zn = (1, 0, 0, 1, 1, 0), Yn will be (0, 1, 1, 1, 1, 0). Our task is this: Observing Yn , make a decision whether Yn is derived from X=(1,1) or not. Based on the general philosophy of the optimal model, the parts model and the universal-null model that we stated earlier, we define these three concepts in detail as follows, according to the current problem setting. optimal model: We conclude X = (1, 1) if Yn passes the following test: P (Yn |X = (1, 1)) ≥ c. (2.2) P (Yn |X = (0, 0), (1, 0), or (0, 1)) parts model: Let Yn1 be the first n elements of Yn , and Yn2 be the second n elements of Yn . Let X1 stand for the first element of X, and X2 be the second element of X. 25 Figure 2.2: Graphical model. Panel “A” is a graphical representation for the non- hierarchical setting in Section 2.3. Panel “B” is a graphical representation for the hierarchical settings in Section 2.4. We conclude X = (1, 1), if Yn passes both of the following tests:  P (Yn1 |X1 =1)   P (Yn1 |X1 =0) ≥ c (2.3)  P (Yn2 |X2 =1) ≥ c  P (Yn2 |X2 =0) universal-null model: We conclude X = (1, 1) if Yn passes the following test: P (Yn |X = (1, 1)) P (Yn |X1 = 1, X2 = 1) = ≥ c. (2.4) P (Yn |X = (0, 0)) P (Yn |X1 = 0, X2 = 0) 2.3.1 Theoretical Comparison of Asymptotic ROC Curves We are interested in the performance of the ROC curves associated with the three models as n goes to infinity. We will study this topic through analyzing the area above the three ROC curves. To ease the notations, we denote the area above the ROC curve associated with the optimal model by AG, the one associated with the 26 parts model by AP , and the one associated with the universal-null model by AW . It is obvious that in the limit, all three ROC curves go to the best point (0,1). Hence, there is no difference between asymptotic ROC curves given by these three models, and AG, AP, AW → 0, as n → ∞. However, what we are interested in is the convergence rate of the three ROC curves. We argue that the ROC curve of the parts model converges to (0,1) exponentially faster than the ROC curve of the universal-null model, and its convergence rate is comparable with that of the optimal model. Theorem 1. AP (1) limn→∞ AW = 0, and the ratio converges to zero in an exponential order. AP (2) AG ≤ C when n is large, where C is a positive constant. Proof. We will first prove the first part of the theorem, (1). Let I(Yn1 ) and I(Yn2 ) denote the numbers of 1’s in Yn1 and Yn2 separately. Then,  (1−)I(Yn1 ) n−I(Yn1 )   (1−)n−I(Yn1 ) I(Yn1 ) ≥ c (2.3) ⇐⇒  (1−)I(Yn2 ) n−I(Yn2 )  (1−)n−I(Yn2 ) I(Yn2 ) ≥ c (1 − )I(Yn1 )+I(Yn2 ) · 2n−I(Yn1 )−I(Yn2 ) (2.4) ⇐⇒ ≥c (1 − )2n−I(Yn1 )−I(Yn2 ) · I(Yn1 )+I(Yn2 ) Since  < 0.5, (2.3) ⇐⇒ min{I(Yn1 ), I(Yn2 )} ≥ cˆ; (2.4) ⇐⇒ (I(Yn1 ) + I(Yn2 )) ≥ cˆ. 27 Let S = {Yn : Yn is derived from X = (1, 1)}, S = {Yn : Yn was derived from X = (1, 0) or (0, 1) or (0, 0)}, S10 = {Yn : Yn is derived from X = (1, 0)}, S01 = {Yn : Yn is derived from X = (0, 1)}, S00 = {Yn : Yn is derived from X = (0, 0)}. We have, S = S10 ∪ S01 ∪ S00 . By (2.1), AP = P ( min{I(Y˜n1), I(Y˜n2 )} ≥ min{I(Yn1), I(Yn2 )} | Yn ∈ S, Y˜n ∈ S). AW = P ( (I(Y˜n1 ) + I(Y˜n2 )) ≥ (I(Yn1 ) + I(Yn2 )) | Yn ∈ S, Y˜n ∈ S). We make the following notations for later use: AP10 = P ( min{I(Y˜n1), I(Y˜n2 )} ≥ min{I(Yn1), I(Yn2 )} | Yn ∈ S, Y˜n ∈ S10). AP01 = P ( min{I(Y˜n1), I(Y˜n2 )} ≥ min{I(Yn1), I(Yn2 )} | Yn ∈ S, Y˜n ∈ S01). AP00 = P ( min{I(Y˜n1), I(Y˜n2 )} ≥ min{I(Yn1), I(Yn2 )} | Yn ∈ S, Y˜n ∈ S00). AW10 = P ( (I(Y˜n1 ) + I(Y˜n2 )) ≥ (I(Yn1 ) + I(Yn2 )) | Yn ∈ S, Y˜n ∈ S10). AW01 = P ( (I(Y˜n1 ) + I(Y˜n2 )) ≥ (I(Yn1 ) + I(Yn2 )) | Yn ∈ S, Y˜n ∈ S01). AW00 = P ( (I(Y˜n1 ) + I(Y˜n2 )) ≥ (I(Yn1 ) + I(Yn2 )) | Yn ∈ S, Y˜n ∈ S00). 28 Note, when n is large, there exist positive constants C1 , C2, C3, C4 , and C such that AP10 = AP01 ≤ C1 · P (I(Y˜n2 ) ≥ I(Yn2) | Yn ∈ S, Y˜n ∈ S10), AP00 ≤ C2 · P (I(Y˜n2 ) ≥ I(Yn2 ) | Yn ∈ S, Y˜n ∈ S00) ≤ C3 · P (I(Y˜n2 ) ≥ I(Yn2 ) | Yn ∈ S, Y˜n ∈ S10), X AP ≤ C4 · APij (i,j)∈{(1,0),(0,1),(0,0)} ≤ C · P (I(Y˜n2 ) ≥ I(Yn2 ) |Yn ∈ S, Y˜n ∈ S10 ), (2.5) where the last inequality above is due to the fact that Y˜n ∈ S10 and Y˜n ∈ S00 influence I(Y˜n2 ) the same way. Since AW ≤ AW10, from (2.5) we have AP P ( min{I(Y˜n1 ), I(Y˜n2)} ≥ min{I(Yn1), I(Yn2 )} | Yn ∈ S, Y˜n ∈ S) = AW P ( (I(Y˜n1 ) + I(Y˜n2 )) ≥ (I(Yn1 ) + I(Yn2 )) | Yn ∈ S, Y˜n ∈ S) P (I(Y˜n2 ) ≥ I(Yn2 ) |Yn ∈ S, Y˜n ∈ S10 ) ≤ C· . AW10 Let {Ui }ni=1 , {Vi }ni=1 , {Gi }ni=1 and {Qi }ni=1 be i.i.d Bernoulli() r.v.’s, then ∀ Yn ∈ S, Y˜n ∈ S10 we have X X X X I(Y˜n2 ) = Vi , I(Yn2 ) = n − Qi , I(Y˜n1 ) = n − Gi , I(Yn1 ) = n − Ui . i i i i Hence, n P (I(Y˜n2 ) ≥ I(Yn2) |Yn ∈ S, Y˜n ∈ S10 ) = P ( X (Vi + Qi ) ≥ n ), i=1 n X AW10 = P ( [(Vi + Qi ) + (Ui − Gi )] ≥ n ). i=1 29 Therefore, P ( ni=1 (Vi + Qi) ≥ n } P AP ≤ C· P { ni=1 [(Vi + Qi ) + (Ui − Gi )] ≥ n ) P AW Pn P ( n1 + Qi ) ≥ 1 } i=1 (Vi = C· 1 Pn P{ n i=1 (Vi + Qi + Ui − Gi ) ≥ 1 ) Let J = V1 + Q1 ; K = V1 + Q1 + U1 − G1 ; HJ (α) = log EeαJ ; HK (α) = log EeαK ; LJ (β) = sup(αβ − HJ (α)); α LK (β) = sup(αβ − HK (α)). α Note LJ (β) and LK (β) are convex functions. Since E[J ] = E[K] = 2 , both LJ (β) and LK (β) achieve their minimum value 0 at β = 2, and are positive for any β 6= 2. Since  < 0.5, 2 < 1, we have inf LJ (β) = LJ (1) > 0; β≥1 inf LK (β) = LK (1) > 0. β≥1 30 For α 6= 0, by Jensen Inequality, E[eαK ] = E[E[eαK |J ]] = E[E[eα(J+(U1−G1 ))|J ]] > E[eαE[J+(U1 −G1 )|J] ] = E[eα(J+E(U1−G1 ))] = E[eαJ ] Thus, HK (α) > HJ (α), ∀α 6= 0. Assume, LK (1) = sup(α − HK (α)) α = α0 − HK (α0). Note that LK (1) > 0 and 0 − HK (0) = 0, thus α0 6= 0. Hence, LJ (1) = sup(α − HJ (α)) α ≥ α0 − HJ (α0 ) > α0 − HK (α0 ) = LK (1). (2.6) 31 By Cramer’s Theorem, n 1 1X lim log P ( (Vi + Qi ) ≥ 1 ) = − inf LJ (β) = −LJ (1), n→∞ n n i=1 β≥1 n 1 1X lim log P ( (Vi + Qi + Ui − Gi ) ≥ 1 ) = − inf LK (β) = −LK (1). n→∞ n n i=1 β≥1 Therefore, Pn AP P ( n1 + Qi ) ≥ 1 ) i=1 (Vi ≤ C· 1 Pn AW P( n i=1 (Vi + Qi + Ui − Gi ) ≥ 1 ) e−(LJ (1)+o(1))n = C· e−(LK (1)+o(1))n We have concluded above that LJ (1) < LK (1) in (2.6), hence AP goes to zero exponentially, as n goes to ∞. AW The first part of the theorem has been proved. Now we will prove the second part. Remember I(Yn1 ) and I(Yn2 ) stand for the numbers of 1’s in Yn1 and Yn2 respectively. The decision rule of the optimal model in (2.2) is equivalent to (1 − )I(Yn1 )+I(Yn2 ) · 2n−I(Yn1 )−I(Yn2 ) · p3  ≥ c (1 − )I(Yn1 )+n−I(Yn2 ) · n−I(Yn1 )+I(Yn2 ) · p0    + (1 − )n−I(Yn1 )+I(Yn2 ) · I(Yn1 )+n−I(Yn2 ) · p1      + (1 − )2n−I(Yn1 )−I(Yn2 ) · I(Yn1 )+I(Yn2 ) · p2 Considering the denominator above is a summation over three items, we flip the numerator and the denominator of LHS above, and get an equivalent decision rule 32 as 1− 2n−2I(Yn1 )−2I(Yn2 )   p0  p3 ·    1− n−2I(Yn2 )   p1 ≤ c   +  p3 ·   (2.7)   p2 1− n−2I(Yn1 )  + p3 ·  Let Z(Yn1 ) = n − 2I(Yn1 ), Z(Yn2 ) = n − 2I(Yn2 ), 1− p0 a= > 1, b0 = ≥ 0,  p3 p1 p2 b1 = ≥ 0, b2 = ≥ 0, p3 p3 and (2.7) becomes b0 · aZ(Yn1 )+Z(Yn2 ) + b1 · aZ(Yn1 ) + b2 · aZ(Yn2 ) ≤ c. (2.8) By (2.1), the Area above ROC curve of the optimal model is   Z(Yn1 )+Z(Yn2 ) Z(Yn1 ) Z(Yn2 )  b0 · a + b1 · a + b2 · a AG = P  Yn ∈ S, Y˜n ∈ S  .  ˜ ˜ ˜ ˜ ≥ b0 · aZ1 (Yn )+Z2 (Yn ) + b1 · aZ1 (Yn ) + b2 · aZ2 (Yn ) 33 Since S = S10 ∪ S01 ∪ S00, and Yn1 and Yn2 are independent,   Z(Yn1 )+Z(Yn2 ) Z(Yn1 ) Z(Yn2 )  b0 · a + b1 · a + b2 · a AG ≥ P  Yn ∈ S, Y˜n ∈ S10   ˜ ˜ ˜ ˜ ≥ b0 · aZ1 (Yn )+Z2 (Yn ) + b1 · aZ1 (Yn ) + b2 · aZ2 (Yn ) ≥ P (Z(Yn1 ) ≥ Z1 (Y˜n ), Z(Yn2 ) ≥ Z2 (Y˜n )| Yn ∈ S, Y˜n ∈ S10 ) = P (Z(Yn1 ) ≥ Z1 (Y˜n )| Yn ∈ S, Y˜n ∈ S10) · P (Z(Yn2 ) ≥ Z2 (Y˜n )| Yn ∈ S, Y˜n ∈ S10 ) = P (I(Y˜n1 ) ≥ I(Yn1 ) |Yn ∈ S, Y˜n ∈ S10 ) · P (I(Y˜n2 ) ≥ I(Yn2 ) |Yn ∈ S, Y˜n ∈ S10 ) 1 ≥ · P (I(Y˜n2 ) ≥ I(Yn2 ) |Yn ∈ S, Y˜n ∈ S10). 2 From (2.5), AP ≤ C · P (I(Y˜n2 ) ≥ I(Yn2 ) |Yn ∈ S, Y˜n ∈ S10). Therefore, AP C · P (I(Y˜n2 ) ≥ I(Yn2 ) |Yn ∈ S, Y˜n ∈ S10) ≤ 1 AG 2 · P (I(Y˜n2 ) ≥ I(Yn2 ) |Yn ∈ S, Y˜n ∈ S10) = 2 · C. 2.3.2 Demonstration of the Three ROC Curves for Finite n Within the problem setting in Section 2.3.1, we did a simple experiment to compare the empirical ROC curves for the three models: the optimal model, the parts model, and the universal-null model. It also shows how ROC curves vary with respect to the different image resolutions, i.e., n. We gave a prior distribution of X: P (X = (1, 1)) = 0.09, P (X = (1, 0)) = P (X = (0, 1)) = 0.21, P (X = (0, 0)) = 0.49. We took  = 0.3 for the superposed noise. 34 We generated N = 100, 000 random samples, and calculated empirical ROC curves, for n = 3, 6, 10, 18. In Figure 2.3, we plotted four figures comparing the three ROC curves, where each figure corresponds to a particular n. Solid line represents the ROC curve for the parts model, dashed line represents the one for the universal-null model, while dotted line represents the one for the optimal model. As shown in the experiment, when n is small (n=3 or 6), the universal-null ROC curve can be better than or cross the parts ROC curve. But when n is large (n ≥ 10) and as n gets larger, the parts model gives strictly better ROC curve than the universal-null model, and it merges to the optimal ROC curve. 2.3.3 Generalization In Section 2.3.1, Theorem 1 was established based on the case that X is composed of two parts. Actually, Theorem 1 remains true when X is composed of more than two parts. Without loss of generality, assume that X contains m parts. The problem setting is similar: X = (X1 , X2 , · · · , Xm ), where {Xi }m i=1 takes binary value 0 or 1. Yn = fn (X) + Zn = (fn (X1 ), fn (X2 ), · · · , fn (Xm )) + Zn , where fn (X) is obtained by repeating X1 n times, followed by repeating X2 n times, . . ., followed by repeating Xm n times. For example, if m = 3, n = 2, X = (1, 0, 1), we will have fn (X) = (1, 1, 0, 0, 1, 1). Zn represents the superposed random noise. Zn is a vector with the same length as fn (X), containing m · n i.i.d Bernoulli() random variables, 1 <  < 0.5. So Zn = (z1 , . . . , zmn ), and P (zi = 1) = , P (zi = 0) = 1 − . The addition between fn (X) and Zn follows the rule of binary summation. For example, if m = 3, n = 2, fn (X) = (1, 1, 0, 0, 1, 0), Zn = (1, 0, 0, 1, 1, 0), Yn will be (0, 1, 0, 1, 0, 0). Our task is this: Given a fixed m, observing Yn , make a decision whether Yn is derived from X=(1,1,. . . ,1) or not. To solve the problem, we can formalize the parts model and the universal-null model in almost the same manner as before. And again by Large Deviation Theory, 35 n=3 n=6 1 1 0.9 0.9 0.8 0.8 0.7 0.7 0.6 0.6 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 n = 10 n = 18 1 1 0.9 0.9 0.8 0.8 0.7 0.7 0.6 0.6 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 Figure 2.3: Comparison of empirical ROC curves of the three models. The upper left panel, the upper right panel, the lower left panel, and the lower right panel show the comparison results associated with n=3, 6, 9, and 12 respectively. The dotted lines represent the ROC curves produced by the optimal ROC curve given by the Neyman-Pearson Lemma; The solid lines represent the ROC curves produced by the parts model; The dashed lines represent the ROC curves produced by the universal-null model. 36 we can achieve the same asymptotic result that the ratio between two AAROCCs of the parts model and the universal-null model goes to zero exponentially, as n goes to infinity. 2.3.4 Note Theorem 1 does not hold if the limit condition is  → 0 instead of n → ∞. Since as  → 0, AP10 = AP01 = C1 · P (I(Y˜n2 ) ≥ I(Yn2 ) | Yn ∈ S, Y˜n ∈ S10) Xn = C1 · P ( (Vi + Qi) ≥ n ) i=1 = C1 · O( (1 − )n ) n Xn AW10 = AW01 = P( [(Vi + Qi) + (Ui − Gi )] ≥ n ) i=1 = O( (1 − )3n ) n (See the Appendix for detailed derivation of the big O approximations above.) Since AP00 ≤ AP10, AW00 ≤ AW10, there exists b > a > 0 such that AP AP10 AP10 ∈ [a · ,b· ]. AW AW10 AW10 Note that AP10 → C, as  → 0. AW10 AP Hence AW is bounded by positive numbers, as  → 0. 37 2.4 Comparison of ROC Curves of the Optimal Model, the Parts Model and the Universal- Null Model within a Hierarchical Setting Consider a simple hierarchical tree structure with only three nodes, γ, α, β. The node γ is the parent node, pointing to two child nodes α and β. The states of brick γ, α, and β are donated as xγ , xα, xβ , where xγ ∈ {0, 1}, xα ∈ {0, 1, . . . , nα}, and xβ ∈ {0, 1, . . . , nβ }. From here through out this chapter, we will use X γ , X α , X β to represent the random variables, while using xγ , xα, xβ to represent the realization of these three random variables. We generate image data Y in the following three steps. First, we randomly select X γ according to a Bernoulli(p) distribution, where p ∈ (0, 1). Second, we generate X α and X β given X γ , as follows: • If X γ = 1, we generate X α and X β according to a joint conditional distribution P (X α , X β |X γ = 1). We require that P (X α = 0 or X β = 0|X γ = 1) = 0, i.e. both X α and X β have to be positive when X γ = 1. • If X γ = 0, we generate X α and X β independently, i.e. P (X α , X β |X γ = 0) = P (X α |X γ = 0) · P (X β |X γ = 0). We require that P (X α = xα|X γ = 0) > 0 and P (X β = xβ |X γ = 0) > 0, ∀xα, xβ . Finally, we generate image Y given xα and xβ , according to P (Y |xα , xβ ). See Figure 2.2B for a graphical representation. The task is to decide whether X γ = 1 or not, given image data Y . To ease the notations, we define the testing rules associated with three models (the optimal model, the parts model, and the universal-null model) as follows: Given a threshold c, X γ = 1 if Fo(Y ) ≥ c, ∀o ∈ {g, p, w}, 38 where “g,” “p,” and “w” stands for the optimal model, the parts model, and the universal-null model respectively. With this hierarchical problem setting, we formal- ize Fg , Fp, and Fw in a slightly different way from Section 2.3: optimal model: P (Y |xα, xβ )P (xα, xβ |X γ = 1) P P (Y |X γ = 1) α x >0,x β >0 Fg (Y ) = = P P (Y |X γ = 0) xα ,xβ P (Y |x , x )P (x , x |X = 0) α β α β γ parts model:  P (Y |xα )P (xα |X α>0) P xα >0    P (Y |X α =0) , denoted by Fpα(Y )     P (Y |xβ )P (xβ |X β >0) P xβ >0 Fp (Y ) = min P (Y |X β =0) , denoted by Fpβ (Y )     P (Y |xα ,xβ )P (xα ,xβ |X γ =1) P denoted by Fˆp (Y )  α >0,X β >0 Px ,   xα >0,xβ >0 P (Y |xα,xβ )P (xα,xβ |X γ =0) Or, alternatively,  maxxα >0 P (Y |xα )P (xα |X α>0)    P (Y |X α =0) , denoted by Fpα (Y )     maxxβ >0 P (Y |xβ )P (xβ |X β >0) Fp(Y ) = min P (Y |X β =0) , denoted by Fpβ (Y )     maxxα >0,X β >0 P (Y |xα,xβ )P (xα,xβ |X γ =1) denoted by Fˆp (Y )  ,   maxxα >0,xβ >0 P (Y |xα ,xβ )P (xα ,xβ |X γ =0) universal-null model: P (Y |xα , xβ )P (xα, xβ |X γ = 1) P P (Y |X γ = 1) xα >0,xβ >0 Fw (Y ) = = P (Y |X α = 0, X β = 0) P (Y |X α = 0, X β = 0) 39 Equivalently (since the Fw (Y ) below is an increasing function of the Fw (Y ) above), α β α β γ P xα >0,xβ >0 P (Y |x , x )P (x , x |X = 1) Fw (Y ) = P xα >0,xβ >0 P (Y |xα, xβ )P (xα , xβ |X γ = 1) + P (Y |X α = 0, X β = 0) 2.4.1 Theoretical Comparison of Asymptotic ROC Curves We want to study mathematically the ROC performance of the three models within this hierarchical setting. A few notations need to be defined for later use. Let {(ik , jk )} be a re-ordered sequence of {1, 2, . . . , nα } × {1, 2, . . . , nβ }, such that ∀k ∈ {1, 2, . . . , (nα · nβ − 1)}, θk ≥ θk+1 , where θk is defined as P (X α = ik , X β = jk |X γ = 1) θk = , ∀k ∈ {1, 2, . . . , nα · nβ }. (2.9) P (X α = ik , X β = jk |X γ = 0) The binary decision problem within this hierarchical setting is called “trivial” if θ1 = θ2 = . . . = θk = . . . = θnα ·nβ −1 = θnα ·nβ . By this definition, nα · nβ = 1 automatically falls into the category of “trivial.” This decision problem is called “trivial” if there exist k1 , k2 ∈ {1, 2, . . . , nα · nβ } such that θk1 6= θk2 . The following theorem addresses the non-trivial case. The trivial case essentially reduces to Theorem 1, as discussed in Section 2.5. Theorem 2. As resolution of data Y goes to ∞, the ROC curves produced by the optimal model and the parts model merge together, and strictly better than the one produced by the universal-null model for non-trivial cases. 40 Proof. For all c ∈ [0, ∞], let O = {Y : generated by X γ = 1}, B = {Y : generated by X γ = 0}, D = {Y : F (Y ) ≥ c}. With the definitions above, TPR and FPR associated with threshold c become |D ∩ O| TPR(c) = , |O| |D ∩ B| FPR(c) = . |B| For the time being, we assume that the decision problem is non-trivial. Let S 00 = {Y : generated by X α = X β = 0}, S 01 = {Y : generated by X α = 0, X β > 0}, S 10 = {Y : generated by X α > 0, X β = 0}, Sk = {Y : generated by X α = ik > 0, X β = jk > 0}, for k ∈ {1, 2, . . . , (nα · nβ )}. Claims without proof: As resolution n goes to infinity, (1) Fp(Y ) −→ 0, ∀ Y ∈ S 00 ∪ S 01 ∪ S 10. (2) Fp(Y ) −→ θk , ∀ Y ∈ Sk , ∀k ∈ {1, 2, . . . , (nα · nβ )}. This indicates that, as n goes to infinity, all the Fp(Y )s will concentrate on finitely many values: {0} ∪ {θk : k ∈ {1, 2, . . . , (nα · nβ )}. Note that θk ≥ θk+1 , hence as the threshold c decreases from ∞ to 0, we stay at the origin first, then all the Y s in S1 pass the test when c = θ1 , followed by all the Y s in S2 passing the test when c = θ2, . . . . . . , followed by all the Y s in S( nα · nβ ) pass test when c = θ( nα · nβ ), and finally 41 all the Y s pass the test when c = 0. Therefore, 4 when c > θ1 , TPR0 = TPR(c) = 0, 4 FPR0 = FPR(c) = 0, 4 |S1 ∩ O| when θ2 < c ≤ θ1 , TPR1 = TPR(c) = = P (X α = i1, X β = j1 |X γ = 1), |O| 4 |S1 ∩ B| FPR1 = FPR(c) = = P (X α = ik , X β = jk |X γ = 0), |B| ············ 4 |(∪kl=1 Sl ) ∩ O| when θk+1 < c ≤ θk , TPRk = TPR(c) = |O| k X = P (X α = il, X β = jl |X γ = 1), l=1 4 |(∪kl=1 Sl ) ∩ B| FPRk = FPR(c) = |B| k X = P (X α = il, X β = jl |X γ = 0), l=1 ············ α β 4 |(∪nl=1·n Sl ) ∩ O| when 0 < c ≤ θnα ·nβ , TPRnα ·nβ = TPR(c) = = 1, |O| α β 4 |(∪nl=1·n Sl ) ∩ B| FPRnα ·nβ = FPR(c) = |B| = P (X > 0, X > 0|X γ = 0), α β 4 when c = 0, TPRnα ·nβ +1 = TPR(c) = 1, 4 FPRnα ·nβ +1 = FPR(c) = 1. α β +1 Hence, as c decreases from ∞ to 0, we get (nα · nβ + 2) points {(FPRk , TPRk )}nk=0·n in the X-Y plane, if we connect these points one by one, we get the asymptotic ROC curve for the parts model. Note that this piecewise linear asymptotic ROC is convex, due to the fact that θk ≥ θk+1 . The reason is that ∀k ∈ {1, 2, . . . , nα · nβ }, θk is just the slope of the k th linear piece of the asymptotic ROC curve, where “k th ” is counted from the point (0,0) to the point (1,1) along the asymptotic ROC curve. 42 All the arguments above work the same for the optimal model, hence in the limit, the parts model and the optimal model share the same ROC curve. For the universal-null model: Denote S 00 = {Y : generated by X α = X β = 0}, S 01 = {Y : generated by X α = 0, X β > 0}, S 10 = {Y : generated by X α > 0, X β = 0}, S = {Y : generated by X α > 0, X β > 0}. Claims without proof: As resolution n goes to infinity, (1) Fw (Y ) −→ 0, ∀ Y ∈ S 00. (2) Fw (Y ) −→ ∞, ∀ Y ∈ S. (3) Fw (Y1 ) < Fw (Y2 ) < Fw (Y3 ), ∀ Y1 ∈ S 00 , Y2 ∈ S 01 ∪ S 10, Y3 ∈ S. For Y ∈ S 01 ∪S 10, we do not have a clear statement about the convergence of Fw (Y ), but the three claims above are enough to portrait the asymptotic ROC curve for the universal-null model. As threshold c decreases from ∞ to 0, we stay at the origin first, then when c < ∞ and c keeps decreasing, all the Y s in S will pass the test, earlier than any Y from the set S 01 ∪ S 10 passes the test. Let M = maxY ∈S 01 ∪S 10 Fw (Y ). 43 Asymptotically we have, when c = ∞, TPR(c) = 0, FPR(c) = 0, |S ∩ O| when M < c < ∞, TPR(c) = = P (X α > 0, X β > 0|X γ = 1) = 1, |O| |S ∩ B| FPR(c) = = P (X α > 0, X β > 0|X γ = 0), |B| |S ∩ O| when 0 < c ≤ M, TPR(c) ≥ = P (X α > 0, X β > 0|X γ = 1) = 1, |O| (hence, here, TPR(c) = 1) |S ∩ B| FPR(c) ≥ = P (X α > 0, X β > 0|X γ = 0), |B| when c = 0, TPR(c) = 1, FPR(c) = 1. Thus, as c decreases from ∞ to 0, only three points are distinguishable and matter for the ROC curve: (0,0), (P (X α > 0, X β > 0|X γ = 0),1), and (1,1). When we connect them one by one, we get the asymptotic ROC curve for the universal-null model, as the resolution n goes to infinity. From the computation above, we can conclude that (1) As the resolution n goes to infinity, the ROC curves of the parts model and the optimal model merge together. Their asymptotic ROC curves are con- vex, piecewise linear, and uniquely determined by the (nα · nβ + 2) points α β +1 {(FPRk , TPRk )}nk=0·n . (2) As the resolution n goes to infinity, the ROC curve of the universal-null model converges to a convex and piecewise linear curve, uniquely determined by the three points (0,0), (P (X α > 0, X β > 0|X γ = 0),1), and (1,1). (3) The asymptotic ROC curve of the parts model is no worse than the asymptotic 44 ROC curve of the universal-null model. For the non-trivial case that nα ·nβ > 1 and ∃ k1 , k2 ∈ {1, 2, . . . , nα ·nβ } such that θk1 6= θk2 , the asymptotic ROC curve of the parts model is strictly better than that of the universal-null model. For the trivial case that nα = nβ = 1 or θ1 = θ2 = . . . = θnα ·nβ , the asymptotic ROC curves of the optimal model, the parts model, and the universal-null model converge to the same convex and piecewise linear curve, that is uniquely determined by the three points (0,0), (P (X α > 0, X β > 0|X γ = 0),1), and (1,1). 2.4.2 Demonstration of the Three ROC Curves for Finite n Within the Hierarchical Setting Consider a simple thought experiment where nα = 2, nβ = 2. The probability parameters were set as follows. P (X γ = 1) = 0.09, P (X γ = 0) = 0.91. P (X α = 1, X β = 1|X γ = 1) = 0.4, P (X α = 1, X β = 2|X γ = 1) = 0.1, P (X α = 2, X β = 1|X γ = 1) = 0.1, P (X α = 2, X β = 2|X γ = 1) = 0.4. P (X α = 0|X γ = 0) = 0.2, P (X α = 1|X γ = 0) = 0.3, P (X α = 2|X γ = 0) = 0.5, P (X β = 0|X γ = 0) = 0.3, P (X β = 1|X γ = 0) = 0.2, P (X β = 2|X γ = 0) = 0.5. Given X α and X β , we gave the data model a similar form as in Section 2.3.2, Yn = (fn (X α ) + Znα, fn (X β ) + Znβ ), where fn (X α ) is obtained by repeating X α n times, while fn (X β ) is obtained by repeating X β n times. Znα and Znβ represent the added random noise. Znα is a vector with length n, containing n i.i.d random variable z, P (z = 0) = 1 − , P (z = k) = /nα , ∀k ∈ {1, . . . , nα }. Znβ is a vector with length n, containing n i.i.d random variable z 0, P (z 0 = 0) = 1 − , P (z 0 = k) = /nβ , ∀k ∈ 45 {1, . . . , nβ }. The addition of the ith element of fn (X α ) (call it “a”) and the ith element of Znα (call it “b”) is defined as the remainder of (a+b) divided by (nα +1). Note that the value of this remainder is in the set {0, 1, . . . , nα }. The addition of the ith element of fn (X β ) and the ith element of Znβ is defined similarly, and takes values in the set {0, 1, . . . , nβ }. For example, if n = 2, X α = 1, X β = 2, Znα = (0, 2), Znβ = (1, 0), we will have Yn = ((1, 1) + (0, 2), (2, 2) + (1, 0)) = (1, 0, 0, 2). Our task is this: Given Yn , make a decision whether X γ = 1 or not. In this experiment, we picked  = 0.3. Let (i1 , i2, i3, i4) = (1, 2, 2, 1), (j1 , j2, j3 , j4 ) = (1, 2, 1, 2), and define P (X α = i1 , X β = j1 |X γ = 1) θ1 = = 6.6667, P (X α = i1 , X β = j1 |X γ = 0) P (X α = i2 , X β = j2 |X γ = 1) θ2 = = 1.6, P (X α = i2 , X β = j2 |X γ = 0) P (X α = i3 , X β = j3 |X γ = 1) θ3 = = 1, P (X α = i3 , X β = j3 |X γ = 0) P (X α = i4 , X β = j4 |X γ = 1) θ4 = = 0.6667. P (X α = i4 , X β = j4 |X γ = 0) We can see that θk ≥ θk+1 , ∀k ∈ {1, 2, 3}. The asymptotic theoretical ROC curve of the optimal model and the parts model is determined by the following six points 46 {(FPRk , TPRk )}5k=0 : when c > θ1, TPR0 = 0, FPR0 = 0, when θ2 < c ≤ θ1, TPR1 = P (X α = 1, X β = 1|X γ = 1) = 0.4, FPR1 = P (X α = 1, X β = 1|X γ = 0) = 0.06, X2 when θ3 < c ≤ θ2, TPR2 = P (X α = il , X β = jl |X γ = 1) = 0.8, l=1 2 X FPR2 = P (X α = il , X β = jl |X γ = 0) = 0.31, l=1 3 X when θ4 < c ≤ θ3, TPR3 = P (X α = il , X β = jl |X γ = 1) = 0.9, l=1 X3 FPR3 = P (X α = il , X β = jl |X γ = 0) = 0.41, l=1 X4 when 0 < c ≤ θ4, TPR4 = P (X α = il , X β = jl |X γ = 1) = 1, l=1 X4 FPR4 = P (X α = il , X β = jl |X γ = 0) = 0.56, l=1 when c = 0, TPR5 = 1, FPR5 = 1. We want to compare the empirical ROC curves of three models for different ns, and to examine whether the empirical ROC curves produced by the optimal model and the parts model will converge to the asymptotic theoretical ROC curve determined by the six points above. To get the empirical ROC curves, we pick four different ns: n = 2, 5, 10, 50. For each n, we randomly generated 100000 random Yn according the generating procedure described earlier, and calculated the empirical ROC curves. For each n, we plotted the three ROC curves in one single figure. Within the same figure we also 47 plotted the six points {(FPRk , TPRk )}5k=0 for reference. Figure 2.4 shows the four cases associated with four different ns. As we can see, as n is large enough, the ROC curves of the optimal model and the parts model merge together, and are strictly better than the ROC curve produced by the universal-null model. Also, when n is large (=50), the empirical ROC curves of the optimal model and the parts model merge to the asymptotic theoretical ROC curve determined by {(FPRk , TPRk )}5k=0 . In addition, when n is large, the ROC curve of the universal-null model also converges in the same way as we concluded in Theorem 2. 2.5 Discussion The Neyman-Pearson Lemma gives us the optimal model for hypothesis testing in object detection or recognition. However, this optimal model is impractical. Not actually impractical for the thought experiment, but impractical in anything resem- bling a real experiment, with multiple parts, multiple objects, variable poses, and so on. Hence we need to seek alternative models. We studied and compared three models in this chapter – the optimal model, the parts model, and the universal-null model. We have shown for both the non- hierarchical setting and the hierarchical setting that the parts model has better ROC performance than the universal-null model and has comparable performance with the optimal model, at least when n is large. From the testing rule associated with the parts model, it seems that, at least in the thought experiments, the parts model does not save us a dramatic amount of computation, compared to the optimal model. However, it does in practice, if we follow a sequence of testings for the parts model. For example, within the setting of Theorem 2, we can examine brick α first, and only move forward to examine brick β if Fpα(Y ) ≥ c, and only move forward further to examine brick γ if Fpβ (Y ) ≥ c. By this way, it decreases the computation to a great extent, compared to the optimal model. This amounts to a simple example 48 n=2 n=5 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 n = 10 n = 50 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 Figure 2.4: Comparison of empirical ROC curves of the three models within the hierarchical setting, for a non-trivial case. The upper left panel, the upper right panel, the lower left panel, and the lower right panel show the comparison results associated with n=2, 5, 10, and 50 respectively. The dotted lines represent the ROC curves produced by the optimal ROC curve given by the Neyman-Pearson Lemma; The solid lines represent the ROC curves produced by the parts model; The dashed lines represent the ROC curves produced by the universal-null model. The six points {(FPRk , TPRk )}5k=0 were plotted as six star points (*) on the X-Y plane in each panel. 49 of the coarse-to-fine algorithm developed in [6, 4, 8]. We introduced two theorems on the asymptotic ROC performance of these three models. Now, the question is how the two theorems, Theorem 1 and Theorem 2, are related with each other. As we concluded earlier in Theorem 2, the asymptotic ROC curve of the parts model is no worse than the asymptotic ROC curve of the universal-null model. Espe- cially, for the non-trivial cases that ∃ k1 , k2 ∈ {1, 2, . . . , nα · nβ } such that θk1 6= θk2 , the asymptotic ROC curve of the parts model is strictly better than the asymptotic ROC curve of the universal-null model. This has been demonstrated in the thought experiment earlier. In Theorem 2, for the trivial cases that θ1 = θ2 = . . . = θnα ·nβ , the asymptotic ROC curves of the optimal model, the parts model, and the universal- null model converge to the same convex and piecewise linear curve (call it η), that is determined by the three points (0,0), (P (X α > 0, X β > 0|X γ = 0),1), and (1,1). In this case, how to compare the ROC performance of the three models? This is where Theorem 1 comes in. Theorem 1 can be adapted to prove that even though the three ROC curves merge together to η, the ROC curve given by the parts model converges to η exponentially faster than the universal-null model. The only difference from the proof of Theorem 1 is that η here is determined by the three points (0,0), (P (X α > 0, X β > 0|X γ = 0),1), and (1,1), while the η in Theorem 1 is determined by the three points (0,0), (0,1), and (1,1). The difference between two ηs from two theorems is due to their different problem settings. In Theorem 2, the confusion resulted by the case when X α > 0, X β > 0 can not be eliminated, due to the fact that P (X α > 0, X β > 0|X γ = 0) > 0. Luckily, this difference does not prevent us to adapt the proof of Theorem 1 to the trivial case of Theorem 2. To demonstrate the connection between the two theorems, we did a simple ex- periment for one of the trivial cases, nα = nβ = 1, in the hierarchical setting of 50 Theorem 2: P (X γ = 1) = 0.09, P (X γ = 0) = 0.91. P (X α = 1, X γ = 1|X γ = 1) = 1, P (X α = 0|X γ = 0) = 0.5, P (X α = 1|X γ = 0) = 0.5, P (X β = 0|X γ = 0) = 0.5, P (X β = 1|X γ = 0) = 0.5. We generated data Yn in the same way as in Section 2.4.2 with  = 0.3. Our task is this: Observing Yn , make a decision whether X γ = 1 or not. First we plotted ROC curves for n=2, 5, 10, 50 in Figure 2.5. We can see that as n gets larger and larger, the three empirical ROC curves converge to the same curve η determined by the three points (0,0), (P (X α > 0, X β > 0|X γ = 0) = 0.25,1), and (1,1). To illustrate better that the ROC curve of the parts model converges to η exponentially faster than that of the universal-null model, we plotted the ROC curves for more ns: n=12, 15, 20, 25, in Figure 2.6. 2.5.1 Appendix The large O derivation in section 2.3.4: AP10 = AP01 = C1 · P (I(Y˜n2 ) ≥ I(Yn2 ) | Yn ∈ S, Y˜n ∈ S10) n X = C1 · P ( (Vi + Qi ) ≥ n ) i=1 n X n X = C1 · P( (Vi + Qi ) = n + k ) k=0 i=1 Xn = C1 · n+k (1 − )n−k . k=0 51 n=2 n=5 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 n = 10 n = 50 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 Figure 2.5: Comparison of empirical ROC curves of the three models within the hierarchical setting, for a trivial case. The upper left panel, the upper right panel, the lower left panel, and the lower right panel show the comparison results associated with n=2, 5, 10, and 50 respectively. The dotted lines represent the ROC curves produced by the optimal ROC curve given by the Neyman-Pearson Lemma; The solid lines represent the ROC curves produced by the parts model; The dashed lines represent the ROC curves produced by the universal-null model. The three star points (*) in each panel correspond to points (0,0), (P (X α > 0, X β > 0|X γ = 0) = 0.25,1), and (1,1). 52 n = 12 n = 15 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 n = 20 n = 25 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 Figure 2.6: Comparison of empirical ROC curves of the three models within the hierarchical setting, for a non-trivial case. The upper left panel, the upper right panel, the lower left panel, and the lower right panel show the comparison results associated with n=12, 15, 20, and 25 respectively. The dotted lines represent the ROC curves produced by the optimal ROC curve given by the Neyman-Pearson Lemma; The solid lines represent the ROC curves produced by the parts model; The dashed lines represent the ROC curves produced by the universal-null model. The three star points (*) in each panel correspond to points (0,0), (P (X α > 0, X β > 0|X γ = 0) = 0.25,1), and (1,1). 53 Since  < 0.5, we have AP10 = AP01 = C1 · O(n (1 − )n ). n X AW10 = AW01 = P ( [(Vi + Qi ) + (Ui − Gi )] ≥ n ) i=1 3n X n X = P( [(Vi + Qi ) + (Ui − Gi )] = n + k ) k=0 i=1 3n XX n n X = P( [(Vi + Qi ) + (Ui − Gi )] = n + k, Gi = j ) k=0 j=0 i=1 3n X X n n X = P( Vi + Qi + Ui = n + k + j, Gi = j ) k=0 j=0 i=1 3n X X n = n+2k+j (1 − )3n−2k−j k=0 j=0 Since  < 0.5, we have AW10 = AW01 = O(n (1 − )3n ). Chapter 3 On the Correctness of Compositional Probabilities Chapter 1 introduced a probabilistic framework for hierarchical generative models, which is composed of a prior distribution on image interpretations and a conditional data distribution on an image given an interpretation. This chapter will focus on the study of the prior distribution, especially to answer the question proposed in Chapter 1, how to perturb a Markovian distribution of interpretations to obtain a non-Markovian distribution that satisfies a set of conditional (attribute) constraints. This question is actually two-fold: First, given the conditional distributions of a set of attribute functions, does there exist any distribution on interpretations that satisfies all of the conditional (attribute) constraints? Second, if there exists such a distribution, how would we achieve it starting with a Markov distribution on interpre- tations? This Chapter will answer these two questions through two theorems. The first theorem will prove that there exits a distribution on interpretations satisfying a “special” set of conditional constraints. The second theorem will give a complete solution to the second question above. Its proof is general and is independent of the proof for the first theorem. We will first review the problem setting and some important concepts defined in 54 55 Chapter 1. Consider a directed acyclic graph G defined by • A vertex for every brick α ∈ B • A directed edge from α to β if β ∈ Ciα for some i ∈ {1, 2, . . . , nα } An “interpretation” ~x is defined as an assignment of states to {xα }α∈B such that α ∈ B \ T and xα > 0 ⇒ xβ > 0 ∀β ∈ Cxαα . (T is the set of terminal bricks defined in Section 1.3.1.) Let I be the set of interpretations. If we declare a brick α “on” when xα > 0, and if we call Cxαα the chosen children of xα > 0, then an interpretation is a state vector ~x in which the chosen children of every non-terminal on brick are themselves on. (See Figure 3.1.) Figure 3.1: Architecture. Left. A hierarchy of “bricks,” each representing a disjunc- tion of conjunctions. Bottom row is the image (pixel) data and the row above it is the set of terminal bricks. The state of a brick signals a chosen set of children. Right. An “interpretation,” which is an assignment of states such that the chosen children of any “on” brick are also on. There can be multiple roots and shared subtrees. Filled circles represent on bricks (non-zero states), and highlighted edges represent chosen children. The Markov (“context-free”) probability of an interpretation ~x ∈ I is defined as Q β β∈B (xβ ) P0 (~x) = Q β , (3.1) x) (1 − 0 ) β∈B(~ where B = B(~x) is the below set defined in Section 1.3.3. 56 α ax (~x) is an attribute function (scaler valued or vector valued) associated with non-terminal brick α ∈ B \ T , when α takes value xα and xα > 0. To ease the notation, from here through the end of this Chapter, we will use aα (·) to stand for α ax (·). We require that aα (~x) depends only on the state of the descendants (i.e., children, grandchildren, great grandchildren, etc.) of α, by the graph structure G. Let Dα be the set of the descendant bricks of α, i.e. for any brick β ∈ Dα , there either exists a directed edge pointing from α to β (‘α → β’), or there exist bricks η1 ,. . . ,ηm s.t. ‘α → η1 → · · · → ηm → β.’ Let xDα be the state of Dα . Then under this notation, aα(~x) = aα(xDα ). A prototypical example is the set of pose coordinates (or relational pose coordi- nates) of the chosen children of α. Depending on the depth of the instantiation of the α children, ax (~x) may depend on the states of bricks that are several generations re- α moved from α itself (grandchildren, great grandchildren, etc.). Pcx (aα |xα) (shorted as Pc (aα|xα)) is the desired conditional distribution on the attribute aα, which will not, in general, agree with P0 defined in (3.1). The questions we address here are: 1. Given {Pc (aα |xα) : xα > 0}α∈B\T , does there exist a distribution P on ~x ∈ I such that P (aα |xα) = Pc (aα|xα), ∀xα > 0? 2. If there exists such a P , how can it be constructed? Theorem 3 (Section 3.1) and Theorem 4 (Section 3.2) will answer these two questions, respectively. 57 3.1 Existence of Probability Distribution Satisfy- ing the Conditional Constraints In this section, we will prove the existence of a distribution on interpretations that satisfies all the conditional constraints under a “positive constraint” assumption: Pc (aα(S)|xα) = 1, α ∈ B \ T , ∀xα > 0, (3.2) where S is a subset of I, defined as S = {~x ∈ I : xβ > 0, ∀β ∈ Dα }. We call {Pc (aα |xα) : xα > 0}α∈B\T a set of “positive conditional constraints,” if it satisfies (3.2). Theorem 3. For the directed acyclic graph G, given any set of positive conditional constraints {Pcα : xα > 0}α∈B\T , there exists at least one distribution P on I such that P (aα |xα) = Pc (aα|xα), ∀α ∈ B \ T , ∀xα > 0. Proof. We will start from a special case and make this case more and more general – We will first prove the existence of P for a special hierarchy structure and a special set of conditional constraints in Step 1. Then we will move one step forward, to prove the existence of P in Step 2 for the same special hierarchy structure (as in Step 1) and the original conditional constraints {Pc (aα|xα ) : xα > 0}α∈B\T . Finally, in Step 3, we will prove the existence of P for the original hierarchy structure as show in Figure 3.1 and the original conditional constraints {Pc (aα|xα) : xα > 0}α∈B\T . To ease the notation, let N be the number of bricks in B. Also denote B as B = {α1 , α2, . . . , αN }. Step 1. We define a special hierarchy structure Gˆ (one example is shown in Figure 3.2) as follows, • A directed edge from brick αi to brick αj if and only if j < i, ∀i, j ∈ {1, . . . , N}, and a special set of given conditional constraints as follows, 58 • {Pc ({xαj }j 0}N i=2 satisfying ∀x αi > 0, Pc ({xαj }j 0, ∀xα1 , . . . , xαi−1 > 0, (3.3) Pc ({xαj }j 1 and xαi > 0. And this will be done through four successive minor steps: Step 1.1, Step 1.2, Step 1.3, and Step 1.4. Figure 3.2: The special hierarchy structure in Step 1 and Step 2, when N=4. Step 1.1. In this step, we only deal with the ancestor brick αN , the brick on top of the hierarchy as shown in Figure 3.2. Let P ({xαj }j 0. Step 1.2. In this step, we will focus on the second to the lowest brick, α2, in addition to brick αN when xαN = 0 (See Figure 3.2). ∀s ∈ {2, 3, . . . , nα1 }, ∀t ∈ {1, 2, . . . , nα2 }, 59 the following equation always holds,  P  nαN P (xα1 = 1, xα2 = t|xαN = k) · P (xαN = k)  k=1  α1 P (x = 1|x α2 = t) + P (xα1 = 1, xα2 = t|xαN = 0) · P (xαN = 0) =  P α  .(3.5) P (xα1 = s|xα2 = t) n N P (xα1 = s, xα2 = t|xαN = k) · P (xαN = k)  k=1  + P (xα1 = s, xα2 = t|xαN = 0) · P (xαN = 0) We require P (x) to satisfy the following constraints: ∀s ∈ {2, 3, . . . , nα1 }, ∀t ∈ {1, 2, . . . , nα2 }, ∀k ∈ {1, 2, . . . , nαN },     P (xα1 = 1|xα2 = t) = Pc (xα1 = 1|xα2 = t),     P (xα1 = s|xα2 = t) = Pc (xα1 = s|xα2 = t),   (3.6) P (xα1 = 1, xα2 = t|xαN = 0) = 2·n1α2 · Pc (xα1 = 1|xα2 = t),         P (xαN = k) = δ, where δ is a tiny positive number.  From (3.3), P (xα1 = 1, xα2 = t|xαN = 0) constructed above is positive. Since from Step 1.1, we have already had P (xα1 = 1, xα2 = t|xαN = k) = Pc (xα1 = 1, xα2 = t|xαN = k), (3.7) P (xα1 = s, xα2 = t|xαN = k) = Pc (xα1 = s, xα2 = t|xαN = k). (3.8) After plugging (3.6), (3.7), and (3.8) into (3.5), we solve out P (xα1 = s, xα2 = t|xαN = 0), ∀s ∈ {2, 3, . . . , nα1 }, ∀t ∈ {1, 2, . . . , nα2 }. Before moving forward to other bricks, we first need to show that P (xα1 = s, xα2 = t|xαN = 0) defined above is proper, i.e. the summation of P (xα1 = s, xα2 = t|xαN = 0) over positive s and positive t is in [0,1]. Since P (xαN = k) = δ and δ is very small, 60 from equation (3.5) we have the following approximation, P (xα1 = 1, xα2 = t|xαN = 0) P (xα1 = 1|xα2 = t) ≈ . P (xα1 = s, xα2 = t|xαN = 0) P (xα1 = s|xα2 = t) It implies that ∀s ∈ {1, . . . , nα1 }, ∀t ∈ {1, . . . , nα2 }, P (xα1 = s, xα2 = t|xαN = 0) > 0 and P (xα1 = s, xα2 = t|xαN = 0) P (xα1 = s|xα2 = t) ≤ 2 · P (xα1 = 1, xα2 = t|xαN = 0) · P (xα1 = 1|xα2 = t) 1 α1 α2 P (xα1 = s|xα2 = t) = P (x = 1|x = t) · nα2 P (xα1 = 1|xα2 = t) 1 = α Pc (xα1 = s|xα2 = t). n 2 Hence, α n 2 X n 1α α n 2 n 1α X α1 α2 αN 1 X X P (x = s, x = t|x = 0) ≤ α Pc (xα1 = s|xα2 = t) ≤ 1. t=1 s=1 n 2 t=1 s=1 Step 1.3. In this step, we will focus on the third to the lowest brick, α3, in ad- dition to brick αN when xαN = 0. ∀s ∈ {1, 2, . . . , nα1 }, ∀t ∈ {1, 2, . . . , nα2 }, ∀r ∈ 61 {1, 2, . . . , nα3 } and (s, t) 6= (1, 1), the following equality holds P (xα1 = 1, xα2 = 1|xα3 = r) P (xα1 = s, xα2 = t|xα3 = r)   P (xα1 = 1, xα2 = 1, xα3 = r|xαN = k) · P (x αN = k) PnαN    +  k=1     P (xα1 = 1, xα2 = 1, xα3 = r|xαN = 0) · P (xαN = 0) =   . (3.9) P (xα1 = s, xα2 = t, xα3 = r|xαN = k) · P (x αN = k) PnαN    +  k=1     P (xα1 = s, xα2 = t, xα3 = r|xαN = 0) · P (xαN = 0) From Step 1.2, we have already defined P (xαN = k) = δ, ∀k ∈ {1, 2, . . . , nαN }. Now we require P (x) to satisfy more constraints as follows: ∀s ∈ {2, . . . , nα1 }, ∀t ∈ {2, . . . , nα2 }, ∀r ∈ {1, . . . , nα3 }, ∀k ∈ {1, 2, . . . , nαN },     P (xα1 = 1, xα2 = 1|xα3 = r) = Pc (xα1 = 1, xα2 = 1|xα3 = r),    P (xα1 = s, xα2 = t|xα3 = r) = Pc (xα1 = s, xα2 = t|xα3 = r),    1 P (xα1 = 1, xα2 = 1, xα3 = r|xαN = 0) = · Pc (xα1 = 1, xα2 = 1|xα3 = r)   2nα3       · mins,t>0 {P (xα1 = s, xα2 = t|xαN = 0)}.  (3.10) From Step 1.2, ∀s ∈ {1, . . . , nα1 }, ∀t ∈ {1, . . . , nα2 }, P (xα1 = s, xα2 = t|xαN = 0) > 0, hence P (xα1 = 1, xα2 = 1, xα3 = r|xαN = 0) defined above is positive. Since from Step 1.1, we have already had P (xα1 = 1, xα2 = 1, xα3 = r|xαN = k) = Pc (xα1 = 1, xα2 = 1, xα3 = r|xαN = k) (3.11) P (xα1 = s, xα2 = t, xα3 = r|xαN = k) = Pc (xα1 = s, xα2 = t, xα3 = r|xαN = k). (3.12) 62 After plugging (3.10), (3.11), and (3.12) into (3.9), we solve out P (xα1 = s, xα2 = t, xα3 = r|xαN = 0). ∀s ∈ {1, . . . , nα1 }, ∀t ∈ {1, . . . , nα2 }, ∀r ∈ {1, . . . , nα3 }. Similar as in Step 1.2, we have to justify that P (xα1 = s, xα2 = t, xα3 = r|xαN = 0) defined above is proper. It P α3 is easy to see that we only need to justify nr=1 P (xα1 = s, xα2 = t, xα3 = r|xαN = 0) ≤ P (xα1 = s, xα2 = t|xαN = 0). Due to the fact that P (xαN = k) = δ and δ is a small positive number, we have the following approximation from equation (3.9), P (xα1 = 1, xα2 = 1, xα3 = r|xαN = 0) P (xα1 = 1, xα2 = 1|xα3 = r) ≈ . P (xα1 = s, xα2 = t, xα3 = r|xαN = 0) P (xα1 = s, xα2 = t|xα3 = r) It implies α n 3 X P (xα1 = s, xα2 = t, xα3 = r|xαN = 0) r=1 α n 3 X P (xα1 = s, xα2 = t|xα3 = r) ≤ 2 · P (xα1 = 1, xα2 = 1, xα3 = r|xαN = 0) · r=1 P (xα1 = 1, xα2 = 1|xα3 = r) n 3α X 1 = α2 P (xα1 = 1, xα2 = 1|xα3 = r) · min{P (xα1 = s, xα2 = t|xαN = 0)} r=1 n s,t>0 P (xα1 = s, xα2 = t|xα3 = r) · P (xα1 = 1, xα2 = 1|xα3 = r) α n 3 α1 α2 αN 1 X = min{P (x = s, x = t|x = 0)} · α3 Pc (xα1 = s, xα2 = t|xα3 = r) s,t>0 n r=1 ≤ P (xα1 = s, xα2 = t|xαN = 0). Step 1.4. For brick α4, α5 , . . . , αN −1 , by going through the similar procedure as in Step 1.3, we construct P ({xαl }jl=1 |xαN = 0), ∀j ∈ {4, . . . , N −1}. And similarly as in Step 1.3, it is not hard to be justified that these constructed conditional probabilities 63 P ({xαl }jl=1 |xαN = 0) are proper and do not contradict with each other. So far, we have constructed N −1 αN N −1 αN P ({xαl }l=1 |x = k) = Pc ({xαl }l=1 |x = k), ∀k > 0, (3.13) and we have defined N −1 αN N −1 P (xαN = k) & P ({xαl }l=1 |x = 0), ∀k > 0, ∀{xαl > 0}l=1 . (3.14) In order to specify the full probability P ({xαl }N l=1 ), we only need to fill in the condi- N −1 tional probability P ({xαl }l=1 | xαN = 0) for the case when there exists at least one xαl equal to zero, and this filling procedure is trivial. According to the construction process of P from Step 1.2 through Step 1.4, automatically, P satisfies P ({xαj }j 0. At this point, we have proved the existence of P for the special hierarchy structure (as shown in Figure 3.2) and the special set of conditional constraints. Now let us move forward, to prove the existence of P for the same hierarchy structure as in Step 1 and the original conditional constraints {Pc (aα|xα) : xα > 0}α∈B\T . Step 2. Consider the same compositional machine structure Gˆ as in Step 1, but the original given conditional constraints {Pc (aα|xα) : xα > 0}α∈B\T . In order to exploit the result from Step 1, we need to make a connection between the special set of con- ditional constraints and the original conditional constraints. The connection is that there necessarily exists at least one conditional distribution P˜c (xα1 , . . . , xαi−1 |xαi ), satisfying (3.3) and (3.4), and ∀i ∈ {2, . . . , N}, ∀xαi > 0, ∀a in the domain of aαi (xα1 , xα2 , . . . , xαi−1 ), Z P˜c (xα1 , . . . , xαi−1 |xαi ) = Pc (aαi = a|xαi ). α {x j }i−1 j=1 , s.t. aαi (xα1 ,xα2 ,...,xαi−1 )=a 64 This existence of P˜c (xα1 , . . . , xαi−1 |xαi ) is guaranteed by the assumption that {Pcα : xα > 0}α∈B\T is a set of “positive conditional constraints.” Hence, as long as we require the distribution P (~x) to satisfy the conditional constraint posed by P˜c as in Step 1, P will satisfy the conditional constraint posed by {Pc (aα|xα) : xα > 0}α∈B\T automatically. At this point, we have proved the existence of P for the special hierarchy structure Gˆ (as shown in Figure 3.2) and the original conditional constraints {Pc (aα|xα) : xα > 0}α∈B\T . Now let us move to the final stage to prove the existence of P for the original hierarchy structure G (as shown in Figure 3.1 ) and the original conditional constraints. Step 3. Consider the original hierarchy structure G (shown in Figure 3.1) and the original conditional constraints {Pc (aα |xα) : xα > 0}α∈B\T . We construct P through the following procedure: (1). Transform the original hierarchy structure G (shown in Figure 3.1) to the spe- cial hierarchy structure Gˆ (as shown in Figure 3.2), by adding more edges into the graph. (2). Do Step 2, and get P (~x). It is easy to see that the resulting distribution P (~x) satisfies ∀α ∈ B \ T , ∀xα > 0, P (aα|xα ) = Pc (aα|xα ). But is it a proper distribution on set of interpretations I defined by the hierarchy structure G? The answer is yes. We reason as follows: Let Iˆ represent the set of interpretations defined by the special hierarchy structure Gˆ in Step 3.(1) above. By the definition of “an interpretation” (refer to the beginning of this Chapter), we have Iˆ ⊆ I, considering that Gˆ was obtained after more directed ˆ hence P (~x) is edged were added to G. From Step 3.(1), this P (~x) is supported on I, supported on I as well. 65 3.2 How to Achieve the Probability Distribution Satisfying the Conditional Constraints We have proved in the previous section the existence of a distribution P on I that satisfies all the conditional constraints {Pc (aα|xα) : xα > 0}α∈B\T , if this set of conditional constraints is a set of “positive” conditional constraints. In point of fact, the attribute functions {aα(~x)}α∈B\T together with the attribute (conditional) distributions {Pc (aα|xα)}α∈B\T will in general under-determine the distribution on I. Hence, given the existence of one distribution P satisfying the desired constraints, there are many such distributions satisfying the desired constraints. In this section, we will present a scheme to iteratively perturb the Markovian distribution P0 defined in (3.1) such that, under a hypothesis of “non triviality,” it will converge to an asymptotic distribution P ∗ on I, and P ∗ will satisfy the conditional constraints as well. “Non trivial” means that ∀α ∈ B \ T , ∀xα > 0, Pc (aα(~x)|xα) has the same support as P0 (aα(~x)|xα), (3.15) where P0 (aα(~x)|xα) is the marginal conditional distribution of aα(~x) under P0 defined in (3.1). Let N be the number of bricks in B, and denote B as B = {α1 , α2 , . . . , αN }. Let 4 Mi = nαi , i.e., xαi ∈ {0, 1, . . . , Mi }, ∀i ∈ {1, . . . , N}. Let Pk be the distribution after 66 k steps of perturbations. We define an infinite sequence of perturbations as follows: 1{xα1 =1} Pc (aα1 (~x)|xα1 = 1)  P1 (~x) = P0 (~x) · P0 (aα1 (~x)|xα1 = 1) 1{xα1 =2} Pc (aα1 (~x)|xα1 = 2)  P2 (~x) = P1 (~x) · P1 (aα1 (~x)|xα1 = 2) .. . 1{xα2 =1} Pc (aα2 (~x)|xα2 = 1)  Pnα1 +1 (~x) = Pnα1 (~x) · Pnα1 (aα2 (~x)|xα2 = 1) .. . PN In general, let Ms = i Mi . With this notation, ∀m ∈ {0, 1, 2, . . .}, ∀t ∈ {1, 2, . . . , Ms }, Pl−1 Pl ∀l ∈ {1, 2, . . . , N}, if t ∈ [ i=1 Mi + 1, i=1 Mi ], we have a general perturbation formula as follows: !1 {xαl =t− l−1 Mi } Pl−1 Pc (aαl (~x)|xαl = t − P i=1 Mi ) i=1 PmMs +t (~x) = PmMs +t−1 (~x) · (3.16) = t − l−1 P PmMs +t−1 (aαl (~x)|xαl i=1 Mi ) with the exception that PmMs +t (~x) = PmMs +t−1 (~x), if αl ∈ T . We stop whenever Pk (aα(~x)|xα) = Pc (aα(~x)|xα), ∀α ∈ B \ T , ∀xα > 0, where k = mMs + t. We continue otherwise. (The “non-triviality” condition defined in (3.15) guarantees that the denominator of the ratio in each perturbation defined above is non-zero if its corresponding numerator is non-zero.) To ease the notation, we will simply use the form in (3.16) for all the brick αl ∈ B, while assuming that !1 {xαl =t− l−1 Pl−1 Pc (aαl (~x)|xαl = t − P i=1 Mi } i=1 Mi ) Pl−1 = 1, if αl ∈ T . PmMs +t−1 (aαl (~x)|xαl = t − i=1 Mi ) Theorem 4. Under the non-triviality condition, if there exists a distribution P (~x) on the interpretations I s.t. ∀α ∈ B \ T , ∀xα > 0, P (aα (~x)|xα ) = Pc (aα(~x)|xα), then the sequence of perturbations {Pk }k defined above gives us a pointwise convergent 67 k→∞ distribution P ∗ on I, i.e. Pk (~x)) −→ P ∗ (~x), ∀~x ∈ I. The asymptotic distribution P ∗ satisfies: ∀α ∈ B \ T , ∀xα > 0, P ∗ (aα(~x)|xα) = Pc (aα(~x)|xα). Proof. We will prove the convergence of {Pk } and the property with P ∗ in three steps, Step 1, Step 2, and Step 3. Step 1. ∀k ∈ {0, 1, . . .}, w.o.l.g., we perturb brick α at xα = j, j > 0, at (k + 1)th step: 1{xα =j} P (aα (~x)|xα = j)  Pk+1 (~x) = Pk (~x) · (3.17) Pk (aα(~x)|xα = j) We want to show D(P kPk+1 ) ≤ D(P kPk ) in Step 1. Let G = {~x : xα = j}. From (3.17) we have four observations: Pk+1 (~x) = Pk (~x), ∀~x ∈ Gc , (3.18) Pk (G) = Pk+1 (G), since Pk (Gc ) = Pk+1 (Gc ), (3.19) Pk+1 (aα(~x)|G) = P (aα(~x)|G), (3.20) Pk (~x|aα(~x), G) = Pk+1 (~x|aα(~x), G), (3.21) where (3.18) through (3.20) are straightforward, and (3.21) can be derived as follows: Pk (~x|G) Pk (~x|aα (~x), G) = . Pk (aα (~x)|G) Pk+1 (~x|G) Pk+1 (~x|aα (~x), G) = Pk+1 (aα (~x)|G) P (aα (~ x)|xα=j) Pk (~x|G) ·Pk (aα(~x)|xα=j) = α P (a (~x)|G) Pk (~x|G) = . Pk (aα (~x)|G) 68 Now we split D(P kPk ) into two parts, P (~x) Z D(P kPk ) = P (~x) log Pk (~x) Zx~ P (~x) P (~x) Z = P (~x) log + P (~x) log Pk (~x) Gc Pk (~x) ZG P (G) · P (a (~x)|G) · P (~x|aα(~x), G) α P (~x) Z = P (~x) log + P (~x) log . G Pk (G) · Pk (a (~x)|G) · P (~x|a (~x), G) α α Gc Pk (~x) With the observations (3.18) to (3.21), we get D(P kPk ) − D(P kPk+1 ) P (G) · P (~x|aα(~x), G) P (G) · P (~x|aα(~x), G) Z Z  = P (~x) log − P (~x) log G Pk (G) · Pk (~x|aα (~x), G) G Pk+1 (G) · Pk+1 (~x|aα(~x), G) P (aα (~x)|G) P (aα(~x)|G) Z Z  + P (~x) log − P (~x) log G Pk (aα(~x)|G) G Pk+1 (aα(~x)|G) Z  P (~x) P (~x) Z + P (~x) log − P (~x) log Gc Pk (~x) Gc Pk+1 (~x) P (aα(~x)|G) Z = P (~x) log G Pk (aα (~x)|G) P (aα (~x)|G) Z = P (G) P (~x|G) log G Pk (aα(~x)|G) P (aα (~x)|G) Z Z = P (G) P (~x|G) log α (~ d~x daα a α {~ α x∈G: a (~ α x)=a } P k (a x )|G) P (aα|G) α Z Z  = P (G) P (~x|G)d~x log da aα x∈G: aα (~ {~ x)=aα } Pk (aα |G) P (aα |G) α Z = P (G) P (aα |G) log da aα Pk (aα|G) = P (G) · D(P (aα |G)kPk (aα |G)) ≥ 0. 69 At this moment, we have proved ∀k ∈ {0, 1, . . .}, D(P kPk ) ≥ D(P kPk+1 ), i.e., {D(P kPk )}∞ k=0 is monotonically decreasing. (3.22) Since D(P kPk ) ≥ 0, ∀k, we conclude that {D(P kPk )}∞ k=0 is converging. Hence {D(P kPk ) − D(P kPk+1 )}∞ k=0 is non-negative and converges to 0. Step 2. ∀l ∈ {1, 2, . . . , N}, ∀j ∈ {1, 2, . . . , Mj }, (i.e., considering brick αl , when xαl = j,) ∀m ∈ {1, 2, . . .}, let klj (m) = mMs + l−1 j x : xαl = j}. P i=1 Mi + j, Gl = {~ D(P kPkj (m)−1 ) − D(P kPkj (m) ) = P (Gjl ) · D(P (aαl |Gjl )kPkj (m)−1 (aαl |Gjl )). l l l Since LHS above is converging to 0, as m → ∞, and P (Gjl ) 6= 0, we conclude D(P (aαl |Gjl )kPkj (m)−1(aαl |Gjl )) is converging to 0, as m → ∞. (3.23) l From here until the end of Step 2, our goal is to find a universal bound for Pkj (m)(aαr |Gsr ) l , Pkj (m)−1 (aαr |Gsr ) l ∀l, r ∈ {1, 2, . . . , N}, ∀j, s ∈ {1, 2, . . . , Mr }, ∀aαr ∈ {aαr (~x) : ~x ∈ I}, where Gsr = {~x : xαr = s}. First, we need to derive a few inequalities. ∀l ∈ {1, 2, . . . , N}, ∀j ∈ {1, 2, . . . , Ml }, ∀aαl ∈ {aαl (~x)}, from the non-trivial assumption of Pc , we have P (aαl |xαl = j) = Pc (aαl |xαl = j) > 0. In addition, ~x has a finite domain, hence the domain of aα(~x) is also finite. Therefore, the number MIN defined below is positive, 4 MIN = min min min P (aαl |xαl = j) > 0. {l : αl ∈ B\T } j∈{1,...,Ml} aαl ∈{aαl (~ x)} 70 By (3.23), ∀ > 0 s.t. (MIN)2 −2Ms < ·2 , (3.24) 128 ln 2 ∃ m∗ s.t. ∀m > m∗ , ∀l, ∀j, D(P (aαl |Gjl )kPkj (m)−1 (aαl |Gjl )) < . l Lemma 12.6.1 from Elements of Information Theory by Cover and Thomas, [25], gives a lower bound for the K-L divergence between any two density functions, 1 D(P1 kP2 ) ≥ kP1 − P2 k21 . 2 ln 2 Here it implies 1 D(P (aαl |Gjl )kPkj (m)−1 (aαl |Gjl )) ≥ kP (aαl |Gjl ) − Pkj (m)−1(aαl |Gjl )k21 l 2 ln 2 l " #2 1 X αl j αl j = |P (a |Gl ) − Pkj (m)−1 (a |Gl )| 2 ln 2 αl l a 1  2 ≥ P (aαl |Gjl ) − Pkj (m)−1(aαl |Gjl ) , ∀aαl . 2 ln 2 l Since ~x has a finite domain, the inequality above implies ∀aαl , √ √ |P (aαl |Gjl ) − Pkj (m)−1 (aαl |Gjl )| ≤ 2 ln 2 · . l 71 √ √ √ √ Hence Pkj (m)−1 (aαl |Gjl ) ∈ (P (aαl |Gjl ) − 2 ln 2 · , P (aαl |Gjl ) + 2 ln 2 · ), and l √ √ √ √ P (aαl |Gjl ) − 2 ln 2 ·  2 ln 2 ·  −1 ∈ ( √ √ , √ √ ). Pkj (m)−1 (aαl |Gjl ) P (aαl |Gjl ) − 2 ln 2 ·  P (aαl |Gjl ) − 2 ln 2 ·  l P (aαl |Gj ) √ √ l 2 ln 2 ·  − 1 ≤ √ √ Pkj (m)−1 (aαl |Gjl ) P (aαl |Gjl ) − 2 ln 2 ·  l √ √ 2 ln 2 ·  ≤ (due to (3.24)) P (aαl |Gjl ) − 12 P (aαl |Gjl ) √ 2 2 ln 2 √ = ·  P (aαl |Gjl ) √ 2 2 ln 2 √ ≤ ·  MIN 4 √ 1 = C·  (< · 2−Ms ), (3.25) 4 where the last inequality in the parentheses is due to the upper bound of  in (3.24). Now we consider all the bricks at klj (m) step, after perturbing the brick l. ∀r ∈ 72 {1, 2, . . . , N}, ∀s ∈ {1, 2, . . . , Mr }, ∀aαr , from the definition of Pkj (m) (~x), l Pkj (m) (Gsr ) − Pkj (m)−1 (Gsr ) l l Z Z = Pkj (m) (~x) − Pkj (m)−1(~x) l l x∈Gsr ~ ~ ∈Gsr x "Z Z # = Pkj (m)(~x) − Pkj (m)−1 (~x) l l x∈Gsr ∩Gjl ~ ~ ∈Gsr ∩Gjl x "Z Z # + Pkj (m) (~x) − Pkj (m)−1(~x) l l ~ ∈Gsr ∩(Gjl )c x ~ ∈Gsr ∩(Gjl )c x "Z # P (aαl (~x)|Gjl ) Z = Pkj (m)−1(~x) · − Pkj (m)−1(~x) x∈Gsr ∩Gjl ~ l Pkj (m)−1 (aαl (~x)|Gjl ) x∈Gsr ∩Gjl ~ l l "Z Z # + Pkj (m)−1 (~x) − Pkj (m)−1 (~x) j l j l ~ ∈Gsr ∩(Gl )c x ~ ∈Gsr ∩(Gl )c x ! P (aαl (~x)|Gjl ) Z = Pkj (m)−1 (~x) · −1 . x∈Gsr ∩Gjl ~ l Pkj (m)−1 (aαl (~x)|Gjl ) l Hence, Z P (aαl (~x)|Gj ) |Pkj (m)(Gsr ) − Pkj (m)−1 (Gsr )| ≤ Pkj (m)−1 (~x) · l − 1 α j l l x s j ~ ∈Gr ∩Gl l Pkj (m)−1 (a (~x)|Gl ) l l √ ≤ C ·  · Pkj (m)−1 (Gsr ∩ Gjl ) l √ ≤ C·  · Pkj (m)−1 (Gsr ), l where C ·  is from (3.25). We can rewrite the inequality above as √ √ Pkj (m)−1 (Gsr ) · (1 − C · ) ≤ Pkj (m) (Gsr ) ≤ Pkj (m)−1 (Gsr ) · (1 + C · ). (3.26) l l l With the inequalities derived by far, now we are ready to look at the ratio between 73 Pkj (m)(aαr |Gsr ) and Pkj (m)−1 (aαr |Gsr ). l l Pkj (m)(aαr |Gsr ) l Z = Pkj (m) (~x|Gsr ) l x∈Gsr : {~ aαr (~ x)=aαr } 1 Z = Pkj (m) (~x) Pkj (m)(Gsr ) x∈Gsr : : aαr (~ {~ x)=aαr } l l 1 Z = Pkj (m) (~x) Pkj (m)(Gsr ) j x∈Gsr ∩Gl : aαr (~ {~ x)=aαr } l l 1 Z + Pkj (m)(~x) Pkj (m) (Gsr ) x∈Gsr ∩(Gjl )c : aαr (~ {~ x)=aαr } l l 1 P (aαl |Gjl ) Z = Pkj (m)−1 (~x) · Pkj (m)(Gsr ) l j x∈Gsr ∩Gl : aαr (~ {~ x)=aαr } l Pkj (m)−1 (aαl (~x)|Gjl ) l 1 Z + Pkj (m)−1(~x). Pkj (m) (Gsr ) x∈Gsr ∩(Gjl )c : aαr (~ {~ x)=aαr } l l By (3.26), Pkj (m) (aαr |Gsr ) l 1 √ Z ≤ s √ · Pkj (m)−1 (~x) · (1 + C · ) Pkj (m)−1 (Gr ) · (1 − C · ) x∈Gsr ∩Gjl : aαr (~ {~ x)=aαr } l l 1 √ Z + √ · Pkj (m)−1 (~x) · (1 + C · ) · (1 − C · ) Pkj (m)−1 (Gsr ) {~ j x∈Gsr ∩(Gl )c : aαr (~ x)=aαr } l l √ αr s 1+C ·  = Pkj (m)−1 (a |Gr ) · √ . l 1−C ·  Similarly, by (3.26), √ αr 1−C ·  Pkj (m) (a |Gsr ) ≥ Pkj (m)−1 (a αr |Gsr ) · √ . l l 1+C·  74 Hence, ∀l, r ∈ {1, 2, . . . , N}, ∀j, s ∈ {1, 2, . . . , Mr }, ∀aαr , Pkj (m) (aαr |Gsr ) √ √ 1−C ·  1+C ·  l ∈ [ √ , √ ] Pkj (m)−1(aαr |Gsr ) 1+C ·  1−C ·  l √ √ ⊂ [1 − 4C · , 1 + 4C · ]. (3.27) Step 3. From the general definition of the perturbation (3.16), we can see that ∀r ∈ {1, 2, . . . , N}, ∀s ∈ {1, 2, . . . , Mr }, there exists krs ∈ {1, 2, . . . , Ms } s.t. Pkj (m)−krs (aαr |Gsr ) = P (aαr |Gsr ), ∀αr . l Hence, by (3.27), krs −1 Pkj (m) (aαr |Gsr ) Y Pkj (m)−i (aαr |Gsr ) l = l Pkj (m)−krs (aαr |Gsr ) i=0 Pkj (m)−(i+1)(aαr |Gsr ) l l √ s √ s ∈ [(1 − 4C · )kr , (1 + 4C · )kr ] √ √ Ms ⊆ [(1 − 4C · )Ms , (1 + 4C · ) ] √ √ ⊆ [1 − 2Ms · 4C · , 1 + 2Ms · 4C · ]. (3.28) (See the Appendix for the derivation of (3.28).) Since Pkj (m)−krs (aαr |Gsr ) = P (aαr |Gsr ), l we have Pkj (m)(aαr |Gsr ) √ √ l ∈ [1 − 2Ms · 4C · , 1 + 2Ms · 4C · ]. P (aαr |Gsr ) (M IN )2 By far, we have proved: ∀ > 0 s.t.  < 128 ln 2 · 2−2Ms , ∃m∗ s.t. ∀m > m∗, ∀l, r ∈ {1, 2, . . . , N}, ∀j, s ∈ {1, 2, . . . , Mr }, ∀aαr , Pkj (m) (aαr |Gsr ) √ √ l ∈ [1 − 2Ms · 4C · , 1 + 2Ms · 4C · ], P (aαr |Gsr ) 75 where √ 2 2 ln 2 C = > 0 MIN MIN = min min min P (aαl |xαl = j) > 0. {l : αl ∈ B\T } j∈{1,...,Ml } aαl ∈{aαl (~ x)} √ The upperbound of  guarantees 2Ms · 4C ·  < 1. For any 0 < δ << 1, take (M IN )2  = 128 ln 2 · 2−2Ms · δ 2, and there exists m∗ s.t. ∀m > m∗ , ∀l, r ∈ {1, 2, . . . , N}, ∀j, s ∈ {1, 2, . . . , Mr }, ∀aαr , Pkj (m)(aαr |Gsr ) √ √ l ∈ [1 − 2Ms · 4C · , 1 + 2Ms · 4C · ] P (aαr |Gsr ) = [1 − δ, 1 + δ]. Since this is true ∀l, ∀j, it is straightforward to rephrase it as: ∀δ > 0 small enough, ∃k ∗, ∀k > k ∗ , ∀r ∈ {1, 2, . . . , N}, ∀s ∈ {1, 2, . . . , Nr }, ∀aαr , Pk (aαr |Gsr ) ∈ [1 − δ, 1 + δ]. P (aαr |Gsr ) This implies Pk (aαr |Gsr ) αr Z αr D(Pk (a |Gsr )kP (aαr |Gsr )) = Pk (aαr |Gsr ) · log da P (aαr |Gsr ) Z ≤ Pk (aαr |Gsr ) log(1 + δ) daαr ≤ δ. Note that D(Pk (aαr |Gsr )kP (aαr |Gsr )) is non-negative, therefore D(Pk (aαr |Gsr )kP (aαr |Gsr )) → 0, uniformly for r, s, aαr . (3.29) On the other hand, ~x has a finite domain, hence there exists a subsequence {Pkl }∞ l=1 76 and a limit probability distribution P ∗ s.t. ∀~x, Pkl (~x) → P ∗ (~x), as l → ∞. This also indicates D(Pkl (aαr |Gsr )kP ∗ (aαr |Gsr )) → 0, as l → ∞. (3.30) Combining (3.30) with (3.29), we have D(P ∗ (aαr |Gsr )kP (aαr |Gsr )) = 0. ∀r, s, aαr . With this equality, we can replace P by P ∗ in the arguments for (3.22) and get D(Pk kP ∗ ) monotonically decreases w.r.t k. But on the other hand we also have D(Pkl kP ∗) → 0, as l → ∞, therefore, D(Pk kP ∗ ) → 0, as k → ∞. Chapter 4 Maximum-Likelihood Templates 4.1 Introduction In Bayesian image analysis, a “prior” probability distribution expresses common or learned knowledge about likely and unlikely interpretations, and a (conditional) data model places a probability on image-based observations given any particular inter- pretation. In conjunction, a prior distribution with a conditional data distribution constitute the “forward” or “generative” model. To the extent that the generative model generates features, as opposed to pixel intensities, the “inverse” or “posterior distribution” on interpretations given images is based on incomplete information; feature vectors are generally insufficient to recover the original intensities. There- fore, it is of great interest to develop a class of data models that generate pixel intensities rather than image features. The models covered in this chapter are based on image-fragment templates, sim- ilar to that introduced by Ullman and his collaborators [26, 27, 28]. Ullman et al. selected as templates the informative image patches that have the highest mu- tual information with the object class they represent, and use these templates for object classification and segmentation. Other researchers applied interest point de- tectors on training images, and then extracted image patches around the interest 77 78 points obtained, [20], [36], [37] and [38]. They used these cropped image patches as templates or template candidates. Heisele et al. designed a SVM algorithm to select good image patches around manually selected seed points – that minimized the foreground-versus-background classification error – as their templates for facial parts, [30, 31]. In all these methods the templates were selected from image patches cropped out of natural images. However, if there actually exists a model that governs the natural image generation from the templates, it would be improper to use the generated image patches as the hidden templates. Allassonniere et al. accordingly developed a Bayesian framework to learn deformable templates and proposed a de- tailed parameter estimation by the Expectation-Maximization (EM) algorithm, [39]. Its application on recognition of handwritten digits showed promising results. How- ever, due to the additive noise model, this generative probabilistic model still suffers from not allowing linear change of the gray-scale map. In addition, it is not possible to use this data model to formulate a credible background model for observations outside the object. The deformable templates may handle small variation on scales, but perform badly on spatial shifts, rotations and big scale variations. Sabuncu et al. adapted the method ([39]) on clustering and registering MRI images, [40]. From a training data set of whole brain MR volumns of subjects, they used the EM algo- rithm to learn two (or three) templates, which corresponded to two (or three) age groups, and used them to partition new datasets. Their model still suffers from the same drawbacks as [39]’s. Similar to [39], we learn templates through a generative probabilistic model. Good performance in most image analysis applications requires some degree of in- variance. In the context of a probability model, the notion of invariance is closely connected to the statistical notion of sufficiency. The normalized correlation between an image patch and an image template is a convenient statistic that is invariant to linear transformations. Only gray-scale images are considered in this chapter. We model pixel grey levels by assuming that their distribution depends only on the 79 normalized correlation. Rather than using a correlation between a template and a given image patch as an extracted feature, we define a conditional-data model on pixel intensities under which the correlation is assumed to be a sufficient statistic. This produces a tractable forward model, and furthermore a fully specified likeli- hood function that can be used to learn templates from image data. A training set of eyes, for example, yields an ensemble of templates of left and right eyes, of familiar and natural character, but not actually coming from any particular individ- uals in the training set. By including mixtures over spatial shifts, scale and rotation variations, our model is able to accommodate invariance with respect to the spatial transformations of objects, in addition to linear transformations of pixel intensities. The same idea can be adapted to model background image patches. When a small window slides across a natural image, most of patches seen through the window are (or are close to) uniform-colored. This indicates that a uniform-color template needs to be included into the mixture model for background image patches. However, the normalized correlation between an image patch and a uniform-color template is not well defined. Certain adaptation and modification needs to be made, and that will discussed later in this Chapter. Chapter 1 has presented the compositional machinery, that is composed of two parts: a prior model on image interpretations, and a data model given the terminal bricks of the hierarchy structure. This chapter will focus on the second part of the compositional machinery, i.e. the data model, and is organized as follows. Section 4.2 will develop a probabilistic framework for modeling image patches based on the concepts of templates and sufficient statistics. A basic probabilistic model with a single template will be introduced, and followed by a generalization to multiple tem- plates. A full generalization will be presented after that, to embrace the invariance to scales, rotations and spatial shifts. All the unknown parameters, including the templates themselves, will be learned through an Expectation-Maximization (EM) 80 algorithm. Section 4.3 will implement this method to model facial parts, and demon- strates its direct application on ethnicity classification of face images. Section 4.4 will focus on the adapted model for background image patches. Finally, Section 4.5 will conclude with a discussion about this model and suggestions about future work. 4.2 Generative Probabilistic Model Key notations: T represents an image template. Y represents an observed gray- scale image patch with n pixels. Each pixel of Y takes its value from the finite set {0, 1, . . . , 255}, unless otherwise specified. S = S(Y, T ) is the normalized correlation between T and Y . (Unless otherwise specified, we will use S(Y, T ), S(Y ), and S interchangeably.) 4.2.1 Basic Probabilistic Model with a Single Template We want to model the probability distribution of Y conditioned on it representing a certain object category. For the purpose of discussion, let this object category be an “eye.” (Unless otherwise specified, P (Y |eye) and P (Y ) are going to be used interchangeably through this chapter, to represent the conditional probability of Y .) Our plan is to model P (Y ) based on S being a sufficient statistics of Y , by as- suming that all the Y s are equally likely given the same S. We will first describe the sampling procedure of Y from P (Y ), and this will uniquely determine the ex- pression of P (Y ). As specified earlier, Y is a discrete variable from the domain {0, 1, . . . , 255}n , but at this moment, let us suppose for a moment that Y is con- tinuous – assuming each pixel of Y takes its value from interval [0, 256) – and this temporary modification will ease the understanding of the sampling procedure of Y . Assume T is an eye template, the same size as Y . S is defined as the normalized correlation between Y and T , hence, it takes its value from the interval [−1, 1]. The sampling procedure of Y is composed of two steps: step 1, S = s is sampled from 81 [−1, +1] according to a certain probability distribution of S, denoted as PS ; step 2, given S = s, Y is sampled uniformly from the set {Y ∈ [0, 256)n : S(Y, T ) = s}, i.e., the set of all possible image patches satisfying S(Y, T ) = s. (Now it becomes clear that why Y ’s continuity is needed for describing its sampling procedure here, since otherwise the set {Y ∈ [0, 256)n : S(Y, T ) = s} might be empty.) The generating procedure of Y described above uniquely determines P (Y ), P (S(Y )) P (Y ) = P (Y, S(Y )) = P (S(Y )) · P (Y |S(Y )) = , ∀Y ∈ [0, 256)n , c(S(Y )) where c(S(Y )) is the combinatorial factor, the measure of the set {Y˜ ∈ [0, 256)n : S(Y˜ ) = S(Y )}., Notice. An essential assumption is made here, and will be in effect through out this chapter: Given S(Y, T ) = s, all the Y s are equally likely. S can be seen as a sufficient statistics of Y . Now let us come back to the discrete case, i.e. Y ∈ {0, 1, . . . , 255}n . From this point until the end of this Chapter, Y will be a discrete variable taking its value from {0, 1, . . . , 255}n : S(Y, T ) = s}. We will model the distribution of a discrete Y similarly as the continuous case, i.e. P (S(Y )) P (Y ) = , ∀Y ∈ {0, 1, . . . , 255}n , (4.1) c(S(Y )) where c(S(Y )) is the combinatorial factor, counting the number of Y˜ s that satisfy S(Y˜ ) = S(Y ), i.e., c(S(Y )) = #{Y˜ ∈ {0, 1, . . . , 255}n : S(Y˜ ) = S(Y )}. But we need to be careful, since the set {Y ∈ {0, 1, . . . , 255}n : S(Y, T ) = s} might be empty for a random s ∈ [−1, +1]. (See the Appendix for the importance sampling method we proposed for sampling discrete Y uniformly from the set {Y ∈ {0, 1, . . . , 255}n : S(Y, T ) = s}.) Some rigorous mathematical derivation is needed here, in order to evaluate P (Y ). Two things need to be specified in (4.1): the marginal distribution of S, P (S(Y )), and the combinatorial factor, c(S(Y )). 82 Let us first look at c(S(Y )). It is not hard to realize that a direct calculation of c(S(Y ) = s) is intractable, considering the high dimensionality of Y . What is needed is an approximation to c(s). Here, we approximate it through the Central Limit Theorem (CLT). Consider a thought experiment: Let Y˜ be another image patch with n pixels. Each pixel of Y˜ is distributed according to i.i.d. uniform distribution (denoted as P˜ ) on {0, 1, . . . , 255}. Since P˜ (Y˜ ) = (1/256)n ∀ Y˜ , the “sufficiency assumption” is exactly true under P˜ - all configurations of Y˜ that give the same value of S(Y˜ ) = s are equally likely. The marginal distribution of S under P˜ is X 1 n P˜ (S = s) = P˜ (Y˜ ) = ( ) · c(s). 256 Y˜ :S(Y˜ )=s It implies, c(s) = (256)n · P˜ (S = s). Let F˜s (S) be the cumulative distribution function of S under P˜ . By the Central Limit Theorem (CLT), as n goes to infinity, F˜s (S) converges to the cumulative function of a normal distribution with mean µ and variance σ 2, where µ = 0, σ = √1 . (See the n Appendix for the detailed derivation.) Hence, when n is large, s √ n Z nx2 F˜s (S = s) ≈ √ e− 2 dx. (4.2) −∞ 2π Since the asymptotic distribution N(µ, σ 2 ) is continuous while the distribution of S is discrete for any finite n, we have to be more careful in order to take advantage of the density function of N(µ, σ 2). Let δ be a small positive number. δk = k · δ, k = 0, ±1, ±2, . . .. Let sδ (Y ) be S(Y ) truncated to the largest δk s.t. δk ≤ S(Y ), i.e., 4 sδ (Y ) = δk∗ , where k ∗ = argmax{k∈{0,±1,±2,...}: δk ≤S(Y )} δk . 83 A few notations are defined as follows for later use: 4 Pδ (sδ ) = P ( S(Y ) ∈ [sδ , sδ + δ) ), sδ ∈ {iδ}∞ i=−∞ 4 P˜δ (sδ ) = P˜ ( S(Y˜ ) ∈ [sδ , sδ + δ) ), 4 cδ (sδ ) = #{ Y : S(Y ) ∈ [sδ , sδ + δ) } = #{ Y˜ : S(Y˜ ) ∈ [sδ , sδ + δ) }. Simply by the definition of the notations above, we have X Pδ (sδ ) = P (Y ) ≈ P (Y ) · cδ (sδ ). (4.3) Y : S(Y )∈[sδ ,sδ +δ) 1 n P˜δ (sδ ) = P˜ (Y˜ ) = ( X ) · cδ (sδ ). (4.4) 256 Y˜ : S(Y˜ )∈[sδ ,sδ +δ) Eqn. (4.4) gives cδ (sδ ) = (256)n · P˜δ (sδ (Y )). After plugging it into (4.3), we get Pδ (sδ ) Pδ (sδ ) P (Y ) ≈ = . (4.5) cδ (sδ ) (256)n · P˜δ (sδ ) In order to get the final expression for P (Y ), we need to specify Pδ (sδ ) and P˜δ (sδ ). From (4.2), when n is large, P˜δ (sδ ) = F˜s (sδ + δ) − F˜s (sδ ) Z sδ +δ √ n − nx2 ≈ √ e 2 dx sδ 2π √ n n(sδ )2 ≈ δ · √ e− 2 2π √ n n(S(Y ))2 ≈ δ · √ e− 2 , (4.6) 2π where the last approximation in (4.6) is due to sδ ≈ S(Y ), by the definition of sδ . 84 Concerning Pδ (sδ ), theoretically, this is a discrete probability function of {S(Y ) : Y ∈ {0, 1, . . . , 255}n }. Given that Fˆs (S = s) converges to a cumulative distribu- tion function of a normal distribution, the discrete domain of S(Y ) is dense in [-1,1] asymptotically, as n goes to infinity. Therefore, when n is large, the cumulative distribution function of S(Y ), namely F (S = s), can be approximated by a cumu- lative distribution function Fˆs (S = s) of a “continuous” distribution with a density function Pˆs (S = s). (The hat “ ˆ ” stands for “continuous.”) Hence, Pδ (sδ ) = F (sδ + δ) − F (sδ ) ≈ Fˆs (sδ + δ) − Fˆs (sδ ) Z sδ +δ = Pˆs (x) dx sδ ≈ δ · Pˆs (sδ ) ≈ δ · Pˆs (S(Y )). (4.7) Combining (4.5), (4.6), and (4.7), we have, Pδ (sδ ) P (Y ) ≈ (4.8) (256)n · P˜δ (sδ ) δ · Pˆs (S(Y )) ≈ √ n(S(Y ))2 (4.9) (256)n · δ · √n e− 2 2π Pˆs (S(Y )) = √ n(S(Y ))2 . (4.10) (256)n · √n e− 2 2π Now we need to specify Pˆs (S), a density function of a continuous random variable on [−1, 1], approximating P (S). Given that Y is generated from an eye template T , it is reasonable to model Pˆs (S = s) to be monotone increasing on s ∈ [−1, 1]. We model Pˆs (S = s) as a backwards truncated exponential distribution. (Actually, this distribution is close to the empirical histogram of S(Y, T ) as well. For example, 85 simply pick T to be one of the training eye image patches, and plot the histogram of the normalized correlation S between T and all the rest of eye image patches. With enough eye image patches, say 100, we will observe that the histogram of S is roughly proportional to a backwards truncated exponential distribution.) Pˆs (S = s) ∝ λ e−λ (1−s) , ∀s ∈ [−1, 1], where λ is the distribution parameter, controlling the steepness. The greater λ is, the more S concentrates at +1. After normalization, Pˆs (S = s) becomes 1 Pˆs (S = s) = −2λ λ e−λ (1−s) , ∀s ∈ [−1, 1]. (4.11) 1−e After plugging (4.11) into (4.10), we get the final expression for P (Y ), 1 1−e−2λ λ e−λ (1−S(Y )) P (Y ) = √ n(S(Y ))2 , (4.12) Q · (256)n · √n e− 2 2π where Q is a number that changes the approximations from (4.8, 4.9) to the equality from (4.12). Q is very close to 1 when n is large. 4.2.2 Generalized Probabilistic Model with Multiple Tem- plates. In general, objects from a particular category can not be well represented by only one single template. Object classes are often defined as a combination of structures with distinct characteristics. Taking the “eye” category for example, there is a big variation of its appearance among humans from different ethnic groups. It is there- fore natural to generalize the model to embrace multiple templates, and learn the associated weights of each mixture together with all the other unknown parameters. Suppose there exist Nt eye templates {T1, . . . , Tt, . . . , Tn }, each with the same 86 P size as Y and each associated with a weight t , where t t = 1. We model the distribution of Y as a mixture over Nt templates, Nt X P (Y |eye) = t · P (Y |Tt), (4.13) t=1 where each P (Y |Tt ) is modeled the same as P (Y ) with a single template from Section 4.2.1. After (4.12) is plugged in, the mixture distribution (4.13) becomes Nt X P (Y |eye) = t · P (Y |Tt) t=1 Nt 1 X 1−e−2λt λt e−λt (1−St (Y )) = t · √ 2 , (4.14) n − nSt (Y ) t=1 Qt · (256)n √2π e 2 where St(Y ) = S(Y, Tt) is the normalized correlation between Y and the tth template Tt; Qt is a constant, and very close to 1 when n is large. Parameter Learning. There are a large amount of unknown parameters in this probabilistic model (4.14), that include all the pixel values of Nt image templates, (t) (t) (t) {Tt = (τ1 , τ2 , . . . , τn )0 }N t=1 , in addition to 2Nt scalar parameters, {λt } and {t }. t As mentioned earlier in this chapter, we learn all the unknown parameters (including the templates themselves) by Expectation-Maximization algorithm. Suppose there (i) (i) (i) are N eye image patches {Yi = (y1 , y2 , . . . , yn )0 }N i=1 for training. Assuming that they are i.i.d. samples from the model P (Y |eye), the likelihood function is: Y P (Y~ |eye) = P (Yi |eye) N X Nt 1 Y 1−e−2λt λt e−λt (1−St (Yi )) = t · √ nSt (Yi )2 . (4.15) n − i=1 t=1 Qt (256)n √2π e 2 Note that S is invariant to linear transformation of templates, hence in order to uniquely identify the templates, we require that all the templates are normalized, 87 with mean 0 and variance 1. To simplify the computation, all the training eye image P (i) P (i)2 patches are pre-normalized as well, such that k yk = 0, k yk = 1. In this way, P (i) (t) St (Yi ) becomes the inner product of Y and Tt , i.e., St (Yi ) = k yk · τk . Since in (4.14) Qt is very close to 1 when n is large, we replace Qt with 1 to simplify our computation. The vector of all the unknown parameters is denoted as ~θ. Expectation Step. ∀i, = P (Xi = t|Yi , ~θ(c) ) (i) Pˆt (c) (c) Pt (Yi ) · t = PN (c) (c) , t t=1 Pt (Yi ) · t where Xi = t means that the tth template Tt generated Yi ; the form θ(c) stands for the “current” guess of θ and 1 (c) (c) (λt )(c) e−(λt) (1−St (Yi )) (c) 1−e−2(λt ) Pt (Yi ) = √ 2 . n − nSt (Yi ) Qt (256)n √2π e 2 Maximization Step. Maximize h i Pˆt · log t · P (Yi |Xi = t, ~θ) (i) XX B = t i   P (i) (t) 1 −λt (1− k yk ·τk ) XX (i) XX (i) 1−e−2λt λt e = Pˆt · log(t ) + Pˆt · log  ,   √ (i) (t) 2 ·τ ) P n( k y n − k k t i t i (256)n √2π e 2 over {t}t , {λt }t , {Tt}t subject to: Nt n n (t) (t) X X X t = 1; τk = 0; (τk )2 = 1. t=1 k=1 k=1 88 It is straightforward to solve out {λt } and {t}, N 1 X ˆ (i) t = P N i=1 t ˆ (i) PN 2λt e−2λt + e−2λt − 1 i=1 Pt λt = · PN P (i) (t) , e−2λt − 1 ˆ (i) i=1 Pt · (1 − yk τk ) where λt can be identified by a simple numerical method, for example, Newton’s method or binary search. The non-trivial part is to solve out the Nt unknown (t) (t) (t) templates, {Tt = (τ1 , τ2 , . . . , τn )0 }N t=1 . For any t ∈ {1, 2, . . . , Nt }, maximizing B t over Tt is equivalent to n n ˆ (t) (t) X X max B, s.t. τk = 0, (τk )2 = 1, (4.16) (t) {τk }nk=1 k=1 k=1 where " # ˆ= X (i) X (i) (t) n X (i) (t) 2 B Pˆt λt · yk τk + ·( yk τk ) i k 2 k This is a constrained quadratic maximization problem, and can be solved through traditional methods, for example, the gradient ascent method. However, there is a better way to attack this problem directly, through matrix decomposition and changing variables. This direct method will yield a closed form solution in one singe step. Let X (i) Xn (i) V = λt Pˆt Yi A= Pˆt Yi (Yi )> . i i 2 V is a column vector; A is a n × n semi-positive definite matrix. Let {(am , em )}nm=1 be the pairs of eigenvalues and corresponding eigen vectors of A. Let {vm }nm=1 and τm }nm=1 be the decomposition coefficients of V and Tt with respect to the new base {ˆ 89 vectors {em }m . We have the following decomposition of A, V , and Tt : X X X A= am em e> m, V = vm em , Tt = τˆm em. m m m Denote e = (1, . . . , 1)> , a column vector with n 1s. Let sm = e> · em , then the constrained maximization problem defined in (4.16) becomes ˆ X X 2 max n B, s.t. τˆm sm = 0, τˆm = 1, (4.17) {ˆ τm }m=1 m m where ˆ = V > TT + T > ATt = X X 2 B t vm τˆm + amτˆm . m m Now the question is how to solve this transformed constrained optimization problem (4.17). We will take the traditional Lagrange Multipliers method. Define X X X X ˆ= L vm τˆm + 2 am τˆm − ζ( τˆm sm ) − η( 2 τˆm − 1). m m m m ˆ w.r.t. τˆm , ζ, and η, and we have Take derivatives of L ∂L ˆ vm − ζsm =0 =⇒ τˆm = . (4.18) ∂ τˆm 2(η − am ) P vm sm X η−am τˆm sm = 0 =⇒ ζ=P s2m . (4.19) m η−am 2 X 2 X vm τˆm =1 =⇒ = 1. (4.20) m m 4(η − am )2 Equation (4.19) implies ζ = 0, due to the fact that vm · sm = 0, ∀m ∈ {1, 2, . . . , n}. (4.21) 90 (See the Appendix for the proof of (4.21)). Now consider the equation w.r.t. η in n−2 (4.20). Generally, there exist (n− 2) roots, {ηr }r=1 , that solve this equation. Among ˆ i.e., these (n − 2) roots of η, we pick η ∗ that maximizes B, ˆ η ∗ = argmax B {ηr }n−2 r=1 X X 2 = argmax vm τˆm + am τˆm {ηr }n−2 r=1 m m 2 2 vm am vm = argmax + . {ηr }n−2 2(η − am ) 4(η − am )2 r=1 After plugging ζ = 0 and η = η ∗ into (4.18), we get ∗ vm τˆm = . 2(η ∗ − am ) And this uniquely determines Tt, X ∗ Tt = τˆm em . m 4.2.3 Further Generalized Probabilistic Model with Multi- ple Scales, Rotations and Location Shifts. This section studies the case when an eye in an image patch is neither horizontally located in the center nor with a fixed scale. To capture the invariance of the object to different scales, rotations and spatial shifts, we generalize the previous model (4.14) to include mixtures over discrete scales, rotations and spatial shifts of templates within the image. Besides shifts by pixels, the discrete location shifts can be fraction of a pixel, that enables the templates to find better locations to fit the image data. With this generalized model, to generate an eye image patch Y , we first randomly select a template T , a scale s, a rotation r, and a spatial shift l, then scale and rotate T according to s and r. Second, we use this scaled and rotated template to generate 91 Y˜ – a sub-region of Y that is determined uniquely by (s, r, l) – within Y through sufficient statistics as before. Finally, we fill the rest of Y , i.e. the complement of Y˜ , with i.i.d. random noise from {0, 1, . . . , 255}. Notice. In this further generalized model, T can have a different size from Y due to the fact that T can scale, rotate, and shift within Y . For example, if Y is 10 × 15, then T can be 3 × 2 or 30 × 40, as long as all the templates have the same size. This indicates that T and Y˜ (mentioned in the sampling procedure above) will have different resolutions. Hence, in order to generate Y˜ from T , we will have to “stretch” T (thinking T as a rubber sheet) to match the pixel coordinate of Y first. We call this process a “projection” of T . This projection – at rotation r, scale s, and spatial shift l – will be obtained through a projection matrix Mr,s,l , whose derivation will be described in the Appendix. Suppose each template Tt is associated with Ns scales and Nr rotations. Let Qs,r be the set of possible discrete spatial shifts of a template under scale s and rotation r within the image patch Y . And ∀l ∈ Qs,r , let gs,r,l (Tt ) be the projection of Tt down to the pixel coordinate of Y , under scale s, rotation r, and spatial shift l. Let t be the prior weight associated with each template Tt , δst be the chance that scale s is selected for Tt , and ηrt be the chance that rotation r is selected for Tt . Hence, PNt PNs t P Nr t t=1 t = 1, and ∀t, s=1 δs = 1, r=1 ηr = 1. Considering an object is almost equally likely to appear anywhere in an image patch, we model the distribution of l 1 is uniform from the set Qs,r , i.e. P (l) = kQs,r k , ∀l ∈ Qs,r , where kQs,r k is the counting measure of set Qs,r . Let Y s,r,l (i.e. Y˜ in the previous paragraph) be the sub-region of Y covered by gs,r,l (Tt). Let n be the total number of pixels in Y , while ns,r,l is the number of pixels in Y s,r,l . Nt X Ns X Nr X X 1 P (Y |eye) = t δst ηrt · P (Y s,r,l |gs,r,l (Tt ))· Po ((Y s,r,l )c ), (4.22) t=1 s=1 r=1 l∈Qs,r kQs,r k where Po ((Y s,r,l )c ) stands for the probability density function of the area of Y that 92 is not covered by gs,r,l (Tt ). Since (Y s,r,l )c is filled with i.i.d. uniform noise, 1 n−ns,r,l Po ((Y s,r,l )c ) = ( ) . (4.23) 256 Now consider P (Y s,r,l |gs,r,l (Tt)). Note that Y s,r,l and gs,r,l (Tt ) contain the same amount of pixels, ns,r,l . It is no different to model the distribution of Y s,r,l given gs,r,l (Tt), from modeling the distribution of an image patch given a single template in Section 4.2.1. Therefore, 1 (1−St,s,r,l (Y s,r,l )) s,r,l 1−e−2λt λt e−λt P (Y |gs,r,l (Tt)) = √ ns,r,l ·(St,s,r,l (Y s,r,l ))2 , (4.24) n Qt,s,r,l · (256)ns,r,l · √s,r,l e− 2 2π where St,s,r,l (Y ) is the normalized correlation between Y s,r,l and gs,r,l (Tt)), λt is the parameter associated with the backward truncated exponential distribution of St,s,r,l (Y ), and Qt,s,r,l is a constant, that is close to 1 when ns,r,l is large enough. Plugging (4.23) and (4.24) back into (4.22), we get 1 s,r,l X 1 n 1−e−2λt λt e−λt (1−St,s,r,l (Y )) P (Y |eye) = t δst ηrt × ( ) √ s,r,l ))2 . (4.25) 256 ns,r,l − ns,r,l ·(St,s,r,l (Y t,s,r,l Qt,s,r,l · √2π e 2 Parameter Learning. All the unknown parameters in (4.25) are learned by the EM algorithm, similarly as in Section 4.2.2. However, the matrix decomposition method in Section 4.2.2 does not work here for updating Tt in M-step. We will describe later an alternative way to update Tt, the direct gradient ascent method. Define gs,r,l (Tt) = Ms,r,l · Tt, where Ms,r,l stands for the projection matrix of Tt down to the coordinate of image patch Y , under scale s, rotation r, and location shift l. Ms,r,l can be adjusted such 93 that the mean of Ms,r,l · Tt is equal to zero. For example, Ms,r,l (i, j) can be simply re- placed with (Ms,r,l (i, j)− the mean of the j th column of Ms,r,l ). (See the Appendix for (i) (i) (i) the derivation of Ms,r,l .) Suppose there are N eye images {Yi = (y1 , y2 , . . . , yn )0 }N i=1 (can be of different sizes) for training. Assuming that they are i.i.d. samples from the density function (4.25), the likelihood function of these N training images is: N Y P (Y~ |eye) = P (Yi |eye) i=1 s,r,l 1 Y X 1 n 1−e−2λt λt e−λt (1−St,s,r,l (Yi )) = t δst ηrt × ( ) √ s,r,l 2 . 256 n ns,r,l ·St,s,r,l (Y i ) i t,s,r,l Qt,s,r,l · √s,r,l e− 2 2π To simply the computation and notation, we pre-normalized all Yis,r,l s.t. it has mean 0 and variance 1, and that leads to (Y s,r,l )> Ms,r,l Tt St,s,r,l (Yis,r,l ) = q i . > > Tt Ms,r,l Ms,r,l Tt Again we approximated Qt,s,r,l by 1, if we assume ns,r,l is large enough for all s, r, l. For example, mins,r,l ns,r,l ≥ 100. Expectation Step. ∀i, Pˆ(t,s,r,l) = P (Xi = (t, s, r, l)|Yi, ~θ(c)) (i) (4.26) 1 (c) (c) 1l∈Qs,r ·  (δst )(c) (ηrt )(c) · P(t,s,r,l)(Yi ) kQs,r k t = P (c) t (c) t (c) (c) , 1 t,s,r,l kQs,r k  t (δ s ) (ηr ) · P (t,s,r,l) (Yi ) where ~θ stands for all the unknown parameters, and (c) 1 (c) (1−St,s,r,l (Yis,r,l )) (c) λt e−λt 1 n 1−e−2λt P(t,s,r,l)(Yi ) = P (Yi |Xi = (t, s, r, l), ~θ(c ) = ( (c) ) √ s,r,l 2 . 256 n ns,r,l ·(St,s,r,l (Yi )) √s,r,l e− 2 2π 1 n Actually, ( 256 ) get canceled in both the numerator and the denominator in (4.27). 94 Maximization Step. Maximize   1 t δst ηrt · P (Yi |Xi = (t, s, r, l), ~θ ) (i) XX B = Pˆ(t,s,r,l) · log i kQs,r k t,s,r,l (i) 1 Pˆ(t,s,r,l) · log( XX = t δst ηrt ) i kQs,r k t,s,r,l   s,r,l 1 −λt (1−St,s,r,l (Yi )) XX (i)  1 n 1−e−2λt λt e + Pˆ(t,s,r,l) · log ( ) ,  √ s,r,l 2 256 ns,r,l − ns,r,l ·(St,s,r,l (Y i )) t,s,r,l i √ 2π e 2 over {Tt}N Nt Nt  t t t=1 , {t }t=1 , {λt }t=1 , {δs }s,t , {ηr }r,t subject to: t Nt X Ns X Nr X t = 1; δst = 1, ηrt = 1, ∀t. t=1 s=1 r=1 It is straightforward to solve out {λt }, {t }, {δst } and {ηrt }, (i) Pˆ(t,s,r,l) P i,s,r,l t = P ˆ (i) i,t,s,r,l P(t,s,r,l) 1 X ˆ (i) = P(t,s,r,l), N i,s,r,l ˆ (i) P i,r,l P(t,s,r,l) δst = P , P ˆ (i) i,s,r,l (t,s,r,l) ˆ (i) P i,s,l P(t,s,r,l) ηrt = P , P ˆ (i) i,s,r,l (t,s,r,l) ˆ (i) P 2λt e−2λt + e−2λt − 1 i,s,r,l P(t,s,r,l) λt = ·P , e−2λt − 1 (i) Pˆ(t,s,r,l) · (1 − St,s,r,l (Yis,r,l )) i,s,r,l where λt can be identified by a numerical searching method, e.g., Newton’s method or binary search. The non-trivial part is to solve out the unknown templates {Tt}N t=1 . t 95 ∀l ∈ {1, 2, . . . , L}, maximizing B over Tt is equivalent to ˆ max B, Nt {Tt }t=1 where, ˆ= (i) h ns,r,l i Pˆ(t,s,r,l) λt · St,s,r,l (Yis,r,l ) + X B · (St,s,r,l (Yis,r,l ))2 2 i,s,r,l   2  s,r,l > s,r,l > (Y ) Ms,r,l Tt ns,r,l  (Yi ) Ms,r,l Tt   Pˆ(t,s,r,l) λt · q i X (i) = + · q T >M > M T 2 T >M > M T i,s,r,l t s,r,l s,r,l t t s,r,l s,r,l t At this moment, its becomes clear that the matrix decomposition in Section 4.2.2 does not work for the maximization problem above. That is due to the non-linear q denominators Tt> Ms,r,l > Ms,r,l Tt , resulted from the projection of Tt . As we mentioned earlier, the gradient ascent method is taken to attack this problem instead. First ˆ w.r.t Tt, take the derivative of B ˆ ∂B X s,r,l  ∂S t,s,r,l (Yi s,r,l ) = λt + ns,r,l · St,s,r,l (Yi ) · , ∂Tt ∂Tt i,s,r,l where q > M Ms,r,l s,r,l Tt > Ms,r,l Yis,r,l Tt> Ms,r,l > Ms,r,l Tt − (Yis,r,l )> Ms,r,l Tt √ ∂St,s,r,l(Yis,r,l ) Tt> Ms,r,l > M s,r,l Tt = ∂Tt Tt> Ms,r,l > Ms,r,l Tt > Ms,r,l Yis,r,l (Y s,r,l )> Ms,r,l Tt > Ms,r,l Ms,r,l Tt = q −q i · > > Tt> Ms,r,l > Ms,r,l Tt Tt> Ms,r,l > Ms,r,l Tt Tt Ms,r,l Ms,r,l Tt > > Ms,r,l Ms,r,l Ms,r,l Tt = q · Yis,r,l − > > · St,s,r,l (Yis,r,l ). Tt> Ms,r,l > Ms,r,l Tt Tt Ms,r,l Ms,r,l Tt 96 Combining two derivatives above, ∂B ˆ X > Ms,r,l (i)   = q · Pˆ(t,s,r,l) λt + ns,r,l · St,s,r,l (Yis,r,l ) Yis,r,l ∂Tt Tt> Ms,r,l > Ms,r,l Tt i,s,r,l X Ms,r,l > Ms,r,l Tt   − · ˆ (i) P λ t + n s,r,l · S (Y t,s,r,l i s,r,l ) St,s,r,l (Yis,r,l ) (4.27) > > (t,s,r,l) i,s,r,l T t M s,r,l M T s,r,l t A note on Model Simplification. The model given by (4.25) accommodates the case of shifting a template by a fraction of a pixel, i.e., Qs,r can includes a fraction, besides an integer. Compared to pure integer shifts, the fraction shifts give the model more flexibility to better accommodate the image data Y . However, on the other hand, this results in a high computation and memory cost in parameter learning, especially in updating Nt templates through gradient ascent given in (4.27). If we only consider integer pixel shifts, the computation and memory cost can be reduced to a great extent, since Ms,r,l and ns,r,l will be independent of l. And (4.27) will become ˆ ∂B X τ Ms,r X (i)   = q ˆ P(t,s,r,l) λt + ns,r · St,s,r,l (Yi ) Yis,r,l s,r,l ∂Tt i,s,r Tt> Ms,r >M T s,r t l∈Qs,r X Ms,r > Ms,r Tt X (i)   − Pˆ λt + ns,r · S (Y t,s,r,l i s,r,l ) St,s,r,l (Yis,r,l ) > > (t,s,r,l) T Ms,r Ms,r Tt l∈Q i,s,r t s,r 97 4.3 Experiments on Learning Facial Part Tem- plates and Applications on Ethnicity Classi- fication 4.3.1 Facial Part Templates We implemented different versions of the maximum-likelihood template model on the Feret Face database. This database is composed of 499 gray-scale face images, each 215 × 214 pixels. Each face image has 15 facial landmarks manually labeled in advance. Figure 4.1 shows twelve face images from this dataset and the corresponding landmarks. With the help of the landmarks, we cropped out different groups of facial parts (left eyes, right eyes, noses and mouths) and scaled them down for training. Two sets of experiments will be described in this section. The first set sought to learn templates from training image patches that belong to the same facial part. We used left eye image patches as the training data for a sequence of experiments (Experiment 1 through Experiment 4). We will later show how the learned template changed with respect to different models. The second set of experiments sought to learn templates from training image patches that belong to different facial parts. The purpose of the second set of experiments was to check whether different types of facial templates will pop up automatically and be associated with the right weights. Here we used a mixture of mouth image patches and nose image patches as the training data, and presented the implementation result in Experiment 5. 1. Experiment for the model with mixtures only over multiple templates, (4.14). The training data was composed of 499 left eye image patches, each with size 12×19. Starting from random initialization, the EM algorithm learned 16 templates each with size 12×19. By random initialization, we mean each pixel of the initial template 98 Figure 4.1: 12 face images from Feret Face database, each with 17 landmarks labeled manually. is i.i.d. random noise. In the parameter initialization for the EM, all the λs were set to 3, and all the s were set to be the same (i.e. each initial template was given the same initial weight). Figure 4.2 shows 70 training left eye images. Figure 4.3 shows the evolution of the first 8 templates as the EM algorithm ran. The first row shows the 8 random initial templates, and the (i + 1)th row shows the updated 8 templates after i runs of the EM algorithm. As shown in Figure 4.3, the EM algorithm converged quickly. Figure 4.4 shows the evolution of the other 8 templates from random noise as the EM algorithm ran. 2. Experiments for the model with mixtures only over multiple templates and spatial shifts. The training dataset was composed of 499 left eye images where the left eye was not necessarily centered in the middle, each with size 15 × 23. The EM algorithm learned 16 templates each with size 12 × 19 from random initialization. λs and s were initialized in the same way as in the previous experiment. Figure 4.5 shows 70 training image patches and the 16 learned templates. To demonstrate that the model with mixtures over multiple templates and spatial shifts yields a better result than the basic model with mixtures only over multiple 99 Figure 4.2: The left panel shows 70 training left eye images. The right panel shows the 16 learned templates. templates (4.14), we did another experiment under the basic model, to learn 16 templates each with the same size, 12 × 19. To cut down the training image patches to size 12 × 19 for the basic model, we trimmed each of the original training images (15 × 23) by a few pixels on the margin. Based on the trimmed training images, the EM algorithm learned 16 templates each with size 12 × 19 under the basic model (4.14). Figure 4.6 compares the two sets of 16 learned templates side by side. It is obvious that the model with mixtures over both multiple templates and spatial shifts yielded sharper templates with greater divergence, than the basic model with mixtures only over multiple templates (4.14). 3. Experiment for the model with mixtures over multiple templates, spatial shifts and scales. The dataset was composed of 499 left eye images, each with a random size ranging 100 Figure 4.3: The evolution of the first 8 templates as the EM algorithm ran. Figure 4.4: The evolution of the last 8 templates as the EM algorithm ran. from 12 × 18 to 18 × 27. The eye was not necessarily located in the center of the training image patch. 16 templates each with size 12 × 19 were learned after 8 runs of the EM algorithm. Only two discrete scales (0.9 and 1.1) were considered for each template, hence, Ns = 2, Nr = 1 in the model (4.25). Due to the fact that two discrete scales are not sufficient to cover a continuous size range of training images, the learned templates were not as sharp as in previous experiments. The quality of learned templates can be improved simply by adding more discrete scales for each template, i.e., by increasing Ns . In Figure 4.8 , we show how the 16 templates were updated from random noise after each run of the EM algorithm. The first row is the initialization, and the (i + 1)th row shows the updated templates after i runs of the 101 Figure 4.5: The left panel shows 70 training left eye images, each with size 15 × 23). The right panel shows the 16 learned templates, each with size 12 × 19, under the model with mixtures over multiple templates and spatial shifts. EM algorithm. As it shows in the figures, the EM algorithm converged quickly for our case. 4. Experiment for the full model with mixtures over multiple templates, spatial shifts, scales and rotations (4.25). The dataset was composed of 499 nose images cropped from 499 face images that had been randomly rotated (the rotation angel ∈ [−10o , 10o ]) and scaled (the scaling factor ∈ [0.3, 0.5]). Hence the nose in each image patch was tilted and not always in the center; the size of each image patch ranged from 16 × 18 to 30 × 33. We implemented the fully generalized model, given in (4.25), on this training dataset. The model parameters were set as follows: 16 templates, each with size 15×18, three discrete scales {1.17, 1, 0.83}, and three discrete spatial shifts {−6.7o , 0o , 6.7o }, 102 Figure 4.6: The left panel shows the 16 learned templates from the model not consid- ering spatial shifts. The right panel shows the 16 learned templates from the model considering spatial shifts templates. hence Nt = 16, Ns = Nr = 3 in (4.25). Figure 4.9 shows 120 training data and the 16 learned templates through the EM algorithm. Figure 4.10 shows the evolution of the 16 templates as the EM algorithm ran. 5. Experiment on a training dataset composed of two types of facial parts, under the model with mixtures over multiple templates and spatial shifts. The training dataset was composed of 499 nose images and 499 mouth images, each with size 13 × 18. The EM algorithm learned 32 templates after 15 runs, each with size 11 × 16. Figure 4.11 shows 120 training image patches and the 32 learned templates. Besides the apparent difference between a nose and a mouth, there was a big variation of facial features and expressions among the training image patches – for example, with or without moustache, smiling or not smiling. Hence, as expected, the learned templates from the EM algorithm revealed this variation of facial features and expressions, besides distinguishing noses from mouths. In addition, both the summation of the weights associated with nose templates and the summation of the weights associated with mouth templates were very close to 0.5, which indicated that our model was properly weighted. Figure 4.12 and 4.13 shows the evolution of the first 16 templates and the last 16 templates as the EM algorithm ran. Again, 103 Figure 4.7: The left panel shows 70 training left eye images, with different sizes ranging from 12by18 to 18by27; eyes are not in the center. The right panel shows the 16 learned templates, each with size 12 by 19, under the model considering spatial shifts and mixtures over two scales 0.9 and 1.1. the first row shows the 16 random initial templates, and the (i + 1)th row shows the updated templates after i runs of the EM algorithm. The EM algorithm converged quickly. 4.3.2 Ethnicity Classification of Face Images The learned templates of facial parts in the previous section can be useful building blocks for a complete face model. For example, one way to achieve this complete model would be to combine these learned templates with information on the rela- tive positions of facial parts. However, these templates are also directly useful by themselves in fulfilling certain computer vision tasks. In this section, we will take the example of ethnicity classification of face images to show the direct usefulness of these learned templates of facial parts. We had a dataset composed of 352 East Asian male faces and 272 South Asian male faces. The task was to distinguish East Asian faces from South Asian faces. Certainly a complete face model could be built for each ethnic group and used for classification. But, we fulfilled this face classifi- cation task exclusively by examining the region around the eyes – i.e. East Asian eyes indicated an East Asian face while South Asian eyes indicated an South Asian 104 Figure 4.8: The evolution of the 16 templates during 8 runs of the EM algorithm, starting from a random initialization. face. In other words, we classified the eye image patches cropped from the face im- ages first, and then used the eye classification result to make an decision about the original face images. The classification of eye image patches was done based on our eye model involving templates. We designed the experiment as follows. First the region of the pair of eyes was cropped out of each face and scaled to have the same height, 10 pixels. Now we had two sets of image patches, 352 East Asian eyes (call it set Ae ) and 272 South Asian eyes (call it set Ai ). We used half (selected randomly) of the image patches from Ae as training data for the East Asian group, and half (selected randomly) of the image patches from Ai as training data for the South Asian group. The other half from Ae and the other half from Ai were merged together, and played the role of testing data. We implemented the model with mixtures over 8 templates (each with size 8 × 20), spatial shifts, and 4 discrete scales (1.2, 1.1, 1, and 0.9). The EM algorithm was performed on both the East Asian training data and South Asian training data, and two models were learned: P (Y |East Asian eyes) and P (Y |South Asian eyes), where Y represented an image patch. Figure 4.14 shows 70 East Asian eyes from set Ae and the corresponding 8 learned templates. As for testing the performance of 105 Figure 4.9: The left panel shows 120 training image patches. The right panel shows the 16 learned templates, each with size 15 × 18, from the fully generalized model, with mixtures over multiple templates, spatial shifts, scales, and rotations. image classification, for Y from the testing data, we classified Y and the associated face image as East Asian if P (Y |East Asian eyes) ≥ P (Y |South Asian eyes), and as South Asian if P (Y |East Asian eyes) < P (Y |South Asian eyes). Figure 4.15 shows 70 South Asian eyes from set Ai and the corresponding 8 learned templates. As shown in the figures, East Asian eye templates and South Asian eye templates looked very different from each other, and captured well the facial features associated with the corresponding ethnic group. Since we only had a few hundred face images, in order to achieve a less biased result, we performed 50 cross-validations. Within each cross validation, we repeated the training and testing procedure described above. Each cross validation gave a correct classification rate for the East Asian group, and a correct classification rate for the South Asian group. These two rates were averaged and recorded as Ri , i ∈ {1, . . . , 50}. Finally, after 50 cross-validations were finished, we averaged all the Rs (R1 through R50 ), giving us a classification rate of 97 percent. 106 Figure 4.10: The evolution of the 16 templates during 8 runs of the EM algorithm, starting from a random initialization. 4.4 Background Template Learning The need for prior models of natural image structures occurs in many machine vi- sion problems including denoising, optical flow, super-resolution, stereo etc., even in object detection. In their attempts at object detection, many researchers have considered only object-class models for object detection while neglecting background models, but mathematically this is essentially equivalent to assuming that all im- ages under the background model have the same probability, i.e. are made of i.i.d. uniform random noise (white noise). It is obvious that the probability of generating anything resembling a natural scene from images with random pixel intensities is ex- tremely low. This suggests that in this state-space of all possible scenes, the region of the space occupied by natural scenes is also extremely low, hence the white noise background model is improper. Far from white noise, the world is delicately structured. When we slide a small window across a natural image, most of what we see through the window is edges and uniform colors. These structured patches are the bricks building up the world. 107 Figure 4.11: The left panel shows 70 training 13 × 18 nose or mouth image patches. The right panel shows the 16 learned templates each with size 11 × 16. Hence, the statistics of small image patches have received extensive attention in the literature. In particular, sparse coding approaches attempt to model structural properties of images in terms of a set of linear filter responses. Olshausen and Field [44] have represented an image patch in terms of a liner combination of learned filters. Welling et al. [46] built a probabilistic model for natural image patches based on the Product-of-Experts (POE) framework [45], taking the product of student- t distributions on linear filter responses, and learned all the unknown parameters including the filter themselves by the principle of maximum likelihood. One of the influential statistical models for natural images is the Markov Random Field (MRF). MRFs define a distribution over images that is based on simple and local interactions between pixels. Roth and Black [47, 48] introduced the Fields of Experts (FOE) model under the MRF framework, taking as the potentials the POE model from [46]. They used the principle of maximum likelihood, as well, to find the optimum filters. However, as in previous models of images based on filter outputs 108 Figure 4.12: The evolution of the first 16 templates as the EM algorithm ran, starting from a random initialization. [49], the probability of an image given the model involves an intractable partition function. This makes learning extremely slow, since MCMC sampling has to be performed before every step of gradient ascent of the log likelihood. (Note: From this point until the end of this chapter, we will use “background” and “natural” interchangeably. By “background” image patches, we mean image patches cropped randomly from random natural images.) As in our modeling of the image patches of object parts, we model background image patches based on templates. Our world is structured, hence we expect the background templates to be structured with similar patterns. We will again only consider gray-scale background images. And we will follow the framework of learning 109 Figure 4.13: The evolution of the last 16 templates as the EM algorithm ran, starting from a random initialization. object templates to learn background templates. However, the correlation between a uniform-color template and an image patch is not defined due to the fact that the variance of a uniform-color template is 0. Hence, certain modification needs to be made in order to accommodate the characteristics of the background image patches. 4.4.1 A Basic Probabilistic Model of Background Image Patches Let Y = (y1 , . . . , yn ) be a background image patch with n pixels. Let P (Y |bg) represent the distribution of Y conditioned on it is a background image patch. We start by giving P (Y |bg) the same expression as in 4.14, which is a mixture over 110 Figure 4.14: The left panel shows 70 East Asian eyes, each with height of 10 pixels. The right panel shows the 8 templates, each with size 8 × 20, learned from the model with mixtures over 4 scales: 1.2, 1.1, 1, and 0.9. normalized (mean 0, variance 1) Nt templates each with the same size as Y : Nt X P (Y |bg) = t · P (Y |Tt) t=1 Nt 1 X 1−e−2λt λt e−λt (1−St (Y )) = t · √ 2 , (4.28) n − nSt (Y ) t=1 Qt · (256)n √2π e 2 where St is the normalized correlation between Y and Tt, ∀t ∈ {1, . . . , Nt }. To capture the big portion of uniform-color background image patches mentioned earlier, we add one more mixture into (4.28), and this new mixture involves a new template, T0, that is uniformly colored. We will define this (Nt + 1)th mixture in the following paragraph and afterwards we will see that T0 is actually an “abstract” template. The (Nt + 1)th mixture is also based on a sufficient statistics S0 , a function of Y . Different from the definition of {St }N t=1 , S0 is defined to be uniquely determined by t the variance of Y : X S0 (Y ) = (yk − y)2 /(2552 · n), k P where y = (1/n) · k yk . From the definition, S0 does not involve T0 explicitly. But it is connected to T0 in the following sense: The smaller S0 is, the more Y looks like 111 Figure 4.15: The left panel shows 70 Indian eyes, each with height of 10 pixels. The right panel shows the 8 templates, each with size 8 × 20, learned from the model with mixtures over 4 scales: 1.2, 1.1, 1, and 0.9. a uniform-color image patch generated from T0, due to the fact that the variance of T0 is equal to zero. This special template T0 is our trick to handle the specialty of background image patches. After adding the (Nt + 1)th mixture into (4.28), we have Nt X P (Y |bg) = t · P (Y |Tt) t=0 Nt 1 X 1−e−2λt λt e−λt (1−St (Y )) = t · √ 2 + P (Y |T0) n − nSt (Y ) t=1 Qt · (256)n √2π e 2 Nt 1 X 1−e−2λt λt e−λt (1−St (Y )) P0 (S0 (Y )) = t · √ + , (4.29) #{Y : S0 (Yˆ ) = S0(Y )} ˆ 2 n − nSt (Y ) t=1 Qt · (256)n √2π e 2 where #{Yˆ : S0 (Yˆ ) = S0 (Y )} is the counting measure of the set {Yˆ : S0 (Yˆ ) = S0(Y )}. (4.29) also indicates the sampling procedure of a background image patch Y : we first randomly select t∗ from the set {0, 1, . . . , Nt } according to their weight {t}N t=0 , followed by sampling the corresponding sufficient statistics St∗ = s from its t marginal distribution, and finally uniformly sample Y from the set {Y ∈ {0, 1, . . . , 255}n : St∗ (Y ) = s}. If Y is generated from the uniform-color template T0, S0(Y ) will be more likely to be close to zero. Hence it is reasonable to model P0 (S0(Y ) = s) as a 112 decreasing function of s. Here we model it as a discretized exponential distribution with parameter λ0 . As for #{Yˆ : S0 (Yˆ ) = S0 (Y )}, we again approximate it through CLT similarly as in Section 4.2.1. Hence, P0 (S0 (Y )) λ0 e−λ0 S0 (Y ) P (Y |T0) = = √ − n2 (S0 −µ0 )2 , (4.30) #{Yˆ : S0(Yˆ ) = S0 (Y )} n Q0 · √2πσ e 2σc 0 where the numerator is from the exponential marginal distribution of S0(Y ), and the denominator is due to the fact that if Yˆ is a random vector with length n and each of it element is a i.i.d. random sample from a uniform distribution on {0, 1, . . . , 255}, then S0(Yˆ ) converges in distribution to a normal distribution N(µ0 , σ02/n) as n goes to infinity, where µ0 = 5461.25/(2552 ) and σ0 = 4884.6/(2552 ). (See the Appendix for the detailed derivation of the asymptotic distribution of S0 (Yˆ ) through CLT.) After plugging (4.30) into (4.29), we get Nt 1 X 1−e−2λt λt e−λt (1−St (Y )) λ0 e−λ0 S0 (Y ) P (Y |bg) = t · √ 2 + 0 · √ − n2 (S0 (Y )−µ0 )2 , (4.31) n − nSt (Y ) n t=1 Qt · (256)n √2π e 2 Q0 · √ 2πσ0 e 2σc where Qt is close to 1 when n is large, ∀t ∈ {0, 1, . . . , Nt }. It is straightforward to generalize the model given in (4.31) one step further, by adding in one layer of mixtures over spatial shifts. We only consider here spatial shifts in the unit of integer pixels. Since scales and rotations are not considered for this moment, {Tt}N t=0 have to be no bigger than Y in order to shift around within t the image patch Y . Let nT be the pixel number of each template Tt , and n be the number of pixels in Y . nT ≤ n. Let Q be the set of all the discretized spatial shifts. Note here Q only contains integers. We assume each spatial shift is equally likely, i.e. 1 P (l) = kQk , ∀l ∈ Q, where kQk is the counting measure of set Q. Let Y l represent the sub-region of Y covered by a template, and (Y l )c represent the rest area of Y that is not covered by a template. Each pixel of (Y l )c is modeled as i.i.d. random 1 n−nT noise from {0, 1, . . . , 255}, hence P ((Y l )c ) = ( 256 ) . The generalized background 113 model becomes: P (Y |bg) X 1 X Nt = t · P ((Y l )c ) · P (Y l |Tt) kQk t=0 l∈Q Nt 1 l X 1 1 n−nT X 1−e−2λt λt e−λt (1−St (Y )) = · ( ) · t · √ nT St (Y l )2 kQk 256 t=1 n Qt · (256)nT · √2πT e− 2 l∈Q l X 1 1 n−nT λ0 e−λ0 S0 (Y ) + · ( ) · 0 · √ n kQk 256 n − T (S0 (Y l )−µ0 )2 Q0 · (256)nT · √2πσT 0 e 2σc2 l∈Q X 1 1 n = · ( ) · l∈Q kQk 256   Nt 1 −λt (1−St (Y l )) −λ0 S0 (Y l ) X 1−e−2λt λt e λ0 e  t · √ nT St (Y l )2 + 0 · √ n , (4.32) n nT − 2σT2 (S0 (Y l )−µ0 )2 t=1 Qt · √ T e− 2 Q0 · √ e c 2π 2πσ0 where Qt is close to 1 when nT is large, ∀t ∈ {0, 1, . . . , Nt }. Let Y l = (y1l , . . . , ynl T ), S0(Y l ) is defined as  2 nT X ykl − y l S0 (Y l ) = , 2552 · nT k=1 PnT where y l = (1/nT ) · k=1 ykl . And St (Y l ) is defined to be the normalized correlation between Y l and Tt , for all t ∈ {1, . . . , Nt }. Parameter Learning. We again apply the EM algorithm to learn all the unknown parameters. As we indicated earlier, there is no T0, just a parameter λ0 governing the distribution on the variance of Y l . Hence only Nt templates {Tt = (τ1 , . . . , τnT )}N t t=1 P Nt need to be learned, besides the other parameters {λt }N Nt t=0 and {t }t=1 , with t t=0  = 1 and t , λt ≥ 0. We also require that {Tt}N t=1 are normalized, with mean 0 and variance t (i) (i) (i) 1. Suppose that we have N background image patches {Yi = (y1 , y2 , . . . , yn )0 }N i=1 for training. Assume that they are i.i.d. samples from the model P (Y |bg) in (4.32), 114 then the likelihood function is: Y P (Y~ |bg) = P (Yi |bg) i Y X 1 1 n = · ( ) · i l∈Q kQk 256   Nt 1 −λt (1−St (Yil )) −λ0 S0 (Yil ) X 1−e−2λt λt e λ0 e t · + 0 · .   √ nT St (Yil )2 √ − nT (S (Y l )−µ )2  n Q˜t · n 0 i 0 t=1 Qt · √ T e− 2 √ T e 2σc2 2π 2πσ0 Since Qt is close to 1 ∀t ∈ {0, . . . , Nt } when nT is large, we replace it with 1 in the likelihood above to simplify our computation. Let θ~ represent the vector of all unknown parameters. Expectation Step. ∀i, ∀t ∈ {0, . . . , Nt }, we define Pˆ(t,l) = P (Xi = (t, l)|Yi, ~θ(c) ) (i) (c) 1 (c) t · kQk · P(t,l)(Yi ) = PN (c) (c) , 1 P t t=1  l∈Q t · kQk · P(t,l) (Yi ) (c) where θ(c) stands for the “current” guess of θ and P(t,l) (Yi ) is defined as (c) (c) 1 (1−St (Yil )) (c) λt e−λt (c) 1 n 1−e−2λt P(t,l) (Yi ) = ( ) · √ nT St (Y l )2 , ∀t ∈ {1, . . . , Nt}; 256 n √ T e− 2 i 2π (c) (c) S0 (Yil ) (c) 1 n λ0 e−λ0 P(0,l) (Yi ) = ( ) · √ n − T (S0 (Yil )−µ0 )2 256 n √ T e 2σc2 2πσ0 115 Maximization Step. Define Nt X X   X ˆ (i) 1 ~ B = P(t,l) · log t · · P (Yi |Xi = (t, l), θ) t=0 l∈Q i kQk   N N 1 l −λt (1−St (Yi ))  1 n 1−e−2λt λt e t X t X (i) (i) X X = Pˆ(t,l) · log(t ) + Pˆ(t,l) · log ( ) ·  √ nT St (Y l )2 256  n i t=0 l,i t=1 l,i √ T e− 2 2π   l X (i) ˆ 1 n λ0 e−λ0S0 (Yi ) + P(0,l) · log ( ) · √ n l . 256 n0 − 2σ02 (S0 (Yi )−µ0 )2 l,i √ e c2πσ0 We want to maximize B w.r.t. {t}t , {λt }t , {Tt}t subject to: Nt n n (t) (t) X X X λt ≥ 0; t = 1; τk = 0; (τk )2 = 1. t=0 k=1 k=1 It is straightforward to solve out {λt } and {t} as follows 1 X ˆ (i) t = Pt,l , ∀t ∈ {0, 1, . . . , Nt }, N i,l ˆ (i) P i,l Pt λ0 = P (i) , Pˆt · S0 (Y l ) i,l i ˆ (i) P 2λt e−2λt + e−2λt − 1 i,l P(t,l) λt = ·P , ∀t ∈ {1, . . . , Nt }, e−2λt − 1 ˆ (i) l i,l P(t,l) · (1 − St (Yi )) where λt can be identified by a simple numerical searching method, e.g. New- ton’s method or Binary search. As for updating the unknown templates {Tt = (τ1 , . . . , τnT )}N t=1 , ∀t ∈ {1, 2, . . . , Nt }, maximizing B over Tt is equivalent to t nT nT (t) (t) X X max ˆ B, s.t. τk = 0, (τk )2 = 1, n (t) T {τk }k=1 k=1 k=1 116 where ˆ= X ˆ (i) h l > n l > 2 i B P(t,l) λt · (Yi ) Tt + · ((Yi ) Tt) . i,l 2 This is a constrained quadratic maximization problem, and can be solved through matrix decomposition and changing variables similarly as in the M-step in Section 4.2.2. Let X (i) Xn (i) V = λt · Pˆ(t,l) · Yil A= · Pˆ(t,l) · Yil · (Yil )> . i,l i,l 2 V is a column vector with length nT , and A is a nT ×nT semi-positive definite matrix. ˆ becomes The original maximization problem of B nT nT (t) (t) X X ˆ max B, s.t. τk = 0, (τk )2 = 1, (t) nT {τk }k=1 k=1 k=1 where ˆ = V > · Tt + T > · A · Tt . B t The following steps to solve this maximization problem above are exactly the same as the M-step for learning eye templates in Section 4.2.2, where matrix eigen decom- position was applied. 4.4.2 Experiment for the Basic Background Model We collected 59 natural images from Internet for training, as shown in Figure 4.16. These images have the same width, 240 pixels, and different heights (ranging from 160 pixels to 360 pixels). We tried to cover a good diversity of background images, from man-made objects to animals to human, from indoor offices to outdoor urban scenes. We cropped out 15×15 non-overlapping image patches from these 59 images, and that yielded 12753 training image patches, i.e., N = 12753, n = 15 · 15. Under the model given by (4.32), the EM algorithm learned from this dataset 32 templates 117 {Tt}32 t=1 each with size 10 × 10, starting from a random initialization. By the end of the EM algorithm, 0 was equal to 0.9092, and λ0 was equal to 174.8316, that means the mean of the exponential distribution of S0 was equal to 0.0057. This indicates that most of the natural image patches look very like uniform-color patches. Figure 4.17 shows 600 training image patches and the 32 learned templates. Figure 4.16: 59 natural images collected from Internet. 4.4.3 Improvement on the Background Model From the learned the background templates in Section 4.4.2, we observe that most of the time, templates came in pairs. For example, both the left-white/right-black 118 Figure 4.17: The left panel shows 500 randomly selected training image patches (15 × 15). The right panel shows the 32 learned templates, each with size 10 × 10, learned from the basic background model. template and the left-black/right-white template showed up, and either of the two can be seen as a flip of the other. Since all the templates are normalized, flipping a template is equivalent to multiplying the original by ’-1’. Another observation is that some templates flipped during the EM iterations, from the 2nd iteration to the 3rd iteration, or from the 3rd iteration to the 4th iteration, as shown in Figure 4.18. This indicates that templates tried to flip themselves back and forth, in order to better fit the training data. Both “Templates come in pairs” and “Templates flipped while the EM algorithm ran” are not surprising, since most of our world is symmetric and well structured. Upon realizing the automatic template flipping, we want to improve the background model by encoding the flipping pattern into the model, to gain better representation of our world. We re-formalize the generating procedure of Y from {Tt}N t=0 . To sample a background image patch Y , we first select t t∗ randomly from the set {0, 1, . . . , Nt } and that determines which template to use. 119 Figure 4.18: The evolution of the 32 templates during 20 runs of the EM algorithm, under the basic background model, where the first row shows the random initializa- tion of the templates. Independently we also select a random spatial shift l and that determines where to put down this template in Y . Now, if t∗ 6= 0, we do the following: With probability 0.5 the template Tt∗ generates Y , and with probability 0.5 the flip of Tt∗ (i.e. −Tt∗ ) generates Y . If t∗ = 0, then no flipping is involved, and Y is generated from T0 directly. Besides this modification involving flipping of templates, we make another im- provement by generalizing the distribution of S0 (Y ) from an exponential distribu- tion to a more flexible distribution, a gamma distribution with parameter λ0 and k. Therefore, the new probability density function of a background image patch Y becomes (Note that St (Y, −Tt) = −St(Y, Tt), ∀t ∈ {1, 2, . . . , Nt }) P (Y |bg)   Nt 1 −λt (1−St (Y l )) 1 −λt (1+St (Y l )) 1 1 n t λt e λt e  1−e−2λt 1−e−2λt X X = ·( ) · · √ nT St (Y l )2 + √ nT St (Y l )2  l∈Q kQk 256 t=1 2 Qt · n √ T e− 2 ˜t · Q n √ T e− 2 2π 2π l X 1 1 n S0 (Y l )k−1 · λk0 · e−λ0 S0 (Y ) / Γ(k) + ·( ) · 0 · √ n . (4.33) kQk 256 nT − 2σT2 (S0 (Y l )−µ0 )2 Q0 · 2πσ0 e √ c l∈Q The set of unknown parameters is composed of {Tt = (τ1 , . . . , τnT )}N Nt t=0 , {λt }t=0 , k, t 120 P Nt and {t}N 2 P P t=0 , subject to t t=0  = 1, t, λt ≥ 0, k τk = 0, and k τk = 1. The EM algorithm is applied again to learn the unknown parameters. Its updating procedure is very similar as that of the model (4.32) in Section 4.4.1. In M-step, t , λt , and k are solved out directly, and matrix eigen decomposition is applied to get the closed form solution for Tt. 4.4.4 Experiment for the Improved Background Model Compared with the experiment section 4.4.2, here we used the same dataset and learned the same amount of templates (32 templates), each with the same size 10×10. The only difference is that we used the improved background model given by (4.33). The learned 0 was equal to 0.9057, and λ0 = 81.0450, k = 0.4831. This implies that the mean of the gamma distribution of S0 (Y ) was equal to 0.0060 (compared to 0.0057 in Section 4.4.2). Figure 4.19 shows 600 training image patches and the 32 learned templates. For the purpose of comparison between the basic background model and the improved background model, Figure 4.19 also shows in the third panel the 32 templates learned from the basic background model given in Section 4.4.2). By comparing two sets of learned templates in Figure 4.19, it is obvious that the improved model produced a bigger diversity among the same amount of templates. For example, under the improved model, templates with texture showed up as well. It is also worthy to notice that the flipping phenomena that appeared from the 2nd to 4th row in Figure 4.17 (the template evolution under the basic model) did not appear in Figure 4.20) (the template evolution under the improved model). 4.5 Discussion In this chapter, we developed a generative probabilistic framework to model image patches, based on deformable representative templates. It can be used to model image patches from both foreground objects and background scenes. The templates 121 Figure 4.19: The left panel shows 500 randomly selected training image patches (15 × 15). The middle panel shows the 32 learned templates from the improved background model, each with size 10 × 10. The right panel shows the 32 learned templates from the basic background model, each with size 10 × 10. learned from the model are useful building blocks in developing models for object detection and recognition. There are a few aspects of this model that can be improved. First, the marginal distribution of the sufficient statistics was modeled as backwards truncated exponen- tial distribution, each with a single parameter λ. This distribution approximated the empirical histogram of the sufficient statistics of the real image data relatively well, but not perfectly. The empirical histogram looked more like a gamma distribution. Its peak was close to +1, but not exactly at +1. Hence it will better fit the training data if the marginal distribution is modeled in a more flexible form – for example, a gamma distribution with two parameters instead of exponential distribution with a single parameter. And we can learn its unknown parameters through the EM al- gorithm, much as what we did with the current model. Second, the shape of the 122 Figure 4.20: The evolution of the 32 templates during 20 runs of the EM algo- rithm under the improved background model, where the first row shows the random initialization of the templates. templates was modeled as a rectangle. An important generalization would be to use non-rectangular templates. For example, we might use oval-shaped templates to model a face and almond-shaped templates to model an eye. One way to achieve this is to manually design the shape of the templates, i.e., to give binary weights, either 0 or 1, to each single pixel of the rectangular templates. Alternatively, a better approach would be to associate a weight with each pixel of the rectangular templates. These weights could then be learned automatically along with the pixel values of templates. Third, the number of templates in our model was manually selected. It has not been studied how to pick an optimal number of templates and how this number influences the model performance. There are a few questions left to answer, e.g., does it always improve performance to have more templates? At this point one possible practical solution would be to compare the performance of several models with different numbers of templates and to pick the best one. Certainly, this type of model selection would be very time consuming. It would be better to do a theoretical analysis on the relation between the number of the templates and the model’s performance, and to select the optimal number of templates accordingly. Another aspect worth attention is data collection. In this chapter, all the facial 123 image patches used in our experiments were cropped from Feret face images with the help of their pre-labeled landmarks. Without these landmarks, we would have to manually collect the facial image patches one by one, and that would be labor intense. The bad news is that in reality useful landmarks are not available in general, let alone good training image patches that are ready to use. However, the good news is that we may be able to take advantage of existing methods to ease the data collection for certain applications. For example, face detection has been thoroughly studied in computer vision; one successful method of face detection has been the Viola-Jones boosting method, [43]. The existing methods work well for face detection, but not face recognition. These face detectors can be used to collect datasets for our model as follows: First we can apply the face detector to target the position of faces in images, then crop the facial parts according to the holistic spatial constraints of the parts layout. These cropped image patches will be a dataset for our model. Certainly, most of the collected image patches will not be calibrated – the interesting region will be off center; objects will have different sizes and rotations. But they are still good enough to provide training and serve the goal of face recognition, due to the fact that our model is invariant to spatial shifts, scales and rotations. 4.6 Appendix 1. How to sample random image uniformly given its normalized correla- tion with a fixed template, in Section 4.2.1. Question: Given gray-scale template T with n pixels, T = {t1, t2, . . . , tn }, and a real number ρ0 ∈ [−1, 1], how would we sample a gray-scale image Y = {y1, y2, . . . , yn } uniformly from the set {Y ∈ {0, . . . , 255}n : S(Y, T ) = ρ0 }? Here S(Y, T ) is the normalized correlation between Y and T , and is defined as Pn 1 Pn Pn 1 Pn i=1 (ti − n j=1 tj ) i=1 (yi − n j=1 yj ) S(Y, T ) = qP qP n 1 Pn 2 n 1 Pn i=1 (ti − n j=1 tj ) i=1 (yi − n j=1 y j )2 124 Solution: Let T = ~t, Y = ~y. Without loss of generality, the original question is equivalent to: given T = ~t ∈ Rn , with ni=1 ti = 0 and ni=1 t2i = 1, we want to P P sample Y = ~y ∈ Rn with ni=1 yi = 0 and ni=1 yi2 = 1 such that ni=1 ti yi = ρ0 . P P P Observe that ~t, ~y ∈ (n−1) dimension unit sphere and cos ∠(~t, ~y ) =< ~t, ~y >= S(T, Y ), ~ denoted by ~z , then a p hence kOAk = ρ0 , kABk = 1 − ρ20 . If we can sample AB, desirable ~y will be ~y = ρ0 · ~t + ~z. Now the question comes down to how to sample ~z . Pn Pn To guarantee that i=1 yi = 0 and i=1 yi2 = 1, ~z should satisfy  Pn i=1 zi = 0      ~z⊥~t    p  kzk = 1 − ρ2 .  0 Let ~e = (1, 1, . . . , 1) ∈ Rn , observe the first two conditions above is equivalent to ~z ⊥α(~t, ~e), where α(~t, ~e) is the space spanned by ~t and ~e. From this observation, if we can sample ~x uniformly from the (n − 1) dimension unit sphere, and project it onto √ 2 c 1−ρ ~ ~ α(t, ~e) , the complement space of α(t, ~e), we will get a desirable ~z as ~z = x~p · kx~pk 0 , where x~p denotes the projected vector from ~x. Detail: (1). How to get a random sample from the (n − 1) dimensional unit sphere: Let ~0 = (0, . . . , 0) ∈ Rn and I be the unit matrix with size n. Sample a n dimensional vector ~v = (v1 , v2, . . . , vn ) ∼ N(~0, I), the n-D Gaussian distribution. Let ~x = ~ v kvk , then ~x is a random sample from (n − 1) dimension unit sphere. 125 c (2). How to project ~x onto α(~t, ~e) : Define a matrix   ~t   ~e       A=  ~0   . ..   .       ~0 n×n Let µ~1 , . . . , µ~k be all the eigenvectors corresponding to the ”0” eigenvalue of A. Define a size n-by-k matrix B = (µ~1 , . . . , µ~k ) and a size n-by-n matrix P = B ·(B ∗ ·B)−1 ·B ∗, then P is our projection matrix, i.e., P ~x = x~p . 2. Approximate P (S(Y )) by Central Limit Theorem, in Section 4.2.1. Let Y = (y1, y2 , . . . , yn ), T = (t1 , t2, . . . , tn ), where, {yi } are i.i.d. random variable P with uniform distribution on {0, 1, . . . , 255}; {ti} are constants with i ti = 0 and P 2 2 P i ti = 1. Denote m = E(yi ), σ = V ar(yi), y = (1/n) yi , and P P (yi − y)ti 1 (ti yi /σ) S(Y ) = corrcoef(Y, T ) = pP =√ ·q . (yi − y)2 n 1 · P (y − y) 2 nσ2 i By Strong Law of Large numbers, (y− m)2 P P (yi − m) → E[(y1 − m)2 ] = σ 2; = → E[y1] − m = 0, n n hence 1 X 2 1 X 2 1 X 2 · (y i − y) = 2 · (y i − m) − ( · (yi − m))2 −→ 1 a.s., as n → ∞. nσ nσ nσ 126 By generalized Central Limit Theorem, X t i yi D −→ N(0, 1), as n → ∞, i σ if ∀ > 0, n X E (ti yi − tim)2 · 1{|tiyi −ti m|>} −→ 0, as n → ∞.   (4.34) i=1 Condition (1) is satisfied if max |ti | −→ 0, as n → ∞. (4.35) i Under our settings, T stands for the normalization for an meaningful image template, hence (2) is true except for the trivial templates, for example T only has a fixed number of pixels with non-zero values as n → ∞. But we don’t consider such trivial templates in our model. Therefore, when n is large, we can approximate the distribution of S(Y ) by N(0, 1/n). 3. Proof for “vm · sm = 0, ∀m = 1, 2, . . . , n,” Eqn. (4.21). Recall that X (i) Xn (i) X V = λl Pˆl Yi A= Pˆl Yi (Yi )> = am em e> m, i i 2 m n-D column vector e = (1, 1, . . . , 1)> , and sm = e> em . By the definition of A and that Yi has been normalized, i.e. e> Yi = 0, we have X e> Ae = 0 =⇒ am e> em e> me = 0 X =⇒ am (e> em )2 = 0. Since am ≥ 0, (e> em )2 ≥ 0, we have ∀m, either e> em = 0 or am = 0. If e> em = 0, 127 then sm = 0, hence vm sm = 0 is proved. If am = 0, since am is the eigen value of A, (i) X Aem = 0 =⇒ Pˆl Yi (Yi )> em = 0 i (i) Pˆl Yi (Yi )> em = 0 X =⇒ e> m i (i) X =⇒ Pˆl (e> 2 m Yi ) = 0. i (i) (i) Since Pˆl ≥ 0, it must be that ∀i, either Pˆl = 0 or e> m Yi = 0, hence (i) X vm = e> mV = Pˆl e> m Yi = 0, i and vmsm = 0 is true automatically. 4. Approximate P0 (S0 (Y )) by Central Limit Theorem, in Section 4.4.1. Let Y = (y1, y2 , . . . , yn ), where {yk } are i.i.d. random variable with uniform dis- tribution on {0, 1, . . . , 255}. Denote m = E(yk ), σ 2 = V ar(yk ), y = (1/n) yk P and 1X 1X 1 X S0 (Y ) = (yk − y)2 = (yk − m)2 − 2 ( (yk − m))2. n n n k k k Observe that " # 1X E (yk − m)2 = Var(y1), n k " # " # 1 X 1 X X E ( (yk − m))2 = E (yk − m)2 + (yk − m)(yj − m) n2 k n2 k k6=j 1 1 = 2 [ nVar(y1) + 0 ] = Var(y1), n n hence, E[S0(Y )] ≈ Var(y1 ) = 5461.25. 128 As for the variance of S0 , from the derivation of E[S0] above we can see that as n → +∞, 1X 1X S0 = (yk − y)2 ≈ (yk − m)2. n n k k Let zk = (yk − m)2, then {z1 , . . . , zn } are i.i.d. uniformly distributed from the set {0.52 , 1.52 , . . . , 126.52 , 127.52 }. Then 1X 1X 1X S0 = (yk − y)2 ≈ (yk − m)2 = zk . n n n k k By Central limit theorem, 1 X D √ (zk − E[z1]) −→ N(0, Var(z1)), as n → +∞. n Since   1X 1 1 X S0 ≈ zk = √ √ (zk − E[z1]) + E[z1], n n n when n is big enough, we can approximate the distribution of S0 by a Gaussian 4884.62 distribution with mean E[z1] = Var(y1 ) = 5461.25 and variance Var(z1)/n = n . 5. Derivation of Ms,r,l , the projection matrix of Tt onto the coordinate of image patch Y , under scale s, rotation r, and location shift l, in Section 4.2.3. gs,r,l (Tt) = Ms,r,l ·Tt , is the projected template on the coordinate of Y . The projection from Tt to gs,r,l (Tt) is composed of two steps. Step 1 scales, rotates, and shifts the pixel grid of Tt under s, r, and l. It moves the (i, j) pixel of Tt to a new location (ˆi, ˆj) w.r.t. the coordinate of Y . Note that i, j ∈ Z, ˆi, ˆj ∈ R. For each pixel (k1 , k2 ) of gs,r,l (Tt), step 2 computes its value by averaging the values of pixels of Tt whose transformed location (ˆi, ˆj) falls in a small neighborhood of (k1 , k2 ). (s,r,l) Step 1: (i, j) −→ (ˆi, ˆj). Step 2: We take the Blackman window method, [42]. 129 For each pixel (k1 , k2 ) of gs,r,l (Tt), its pixel value is equal to X 1{dij ≤1} · [0.42 + 0.5 cos(π · dij ) + 0.08 cos(2π · dij )] · T (i, j), (ˆi,ˆ j) where  k(k1 ,k2 )−(ˆi,ˆ   R·bs j)k2 , if k(k1 , k2) − (ˆi, ˆj)k2 ≤ R · bs , dij = ,  0,  otherwise. where k · k2 above stands for L-2 norm on the Y coordinate; bs is the sth scaling factor for Tt; R is a constant, in charge of the size of the neighborhood of (k1 , k2 ), and we pick R = 1.6. Ms,r,l is uniquely determined by the two steps above. Chapter 5 Conclusion and Future Directions Conclusion We believe the “hierarchy of reusabe parts” moves us closer to bridging the ROC gap which persistently occurs in comparison of human and machine performance in vision. As we know, with quite a bit of training, current machines can achieve a reasonable detection rate at an allowable false positive rate. Nevertheless, every percentage of improvement towards perfect detection (i.e. no missing target) is overwhelmed by massive false positive targets. The underlying reason is that there was never a proper background model. In some sense, people all try to separate background from foreground. They model the generative background as simple as white noise or as complicated as MRF depending on the assumption, while modeling the foreground object in a totally different manner, one not compatible with the background model. We all know the mistakes (false positive) occur at the most ambiguous region (the cluttered region). In this region, the image patch looks like everything including the target to the machine if you have to set one threshold for a test of object versus non-object. It is easy to see that those cluttered regions are made of the same parts as the foreground. This prevailing sharing phenomena paired with predominant false detection in this region suggests an object equipped with its own background model, that is, a model that accommodate both objects 130 131 and background which consist of the same reusable parts, in a uniform way. The composition system embodies this idea, and we have shown that theoretically the ROC gap can be narrowed more than ever. The maximum likelihood template gives a solution to how to model image data given the upper level image interpretations, and also provides a way to model objects and background on a same platform. Future Directions 1. Continuous study of maximum likelihood templates. The joining of constituents to make a composition (e.g. a left eye, a right eye, and a forehead to make a top half of a face) will in general depend on the poses (e.g. position, scale, rotation) of the constituents. This dependency can be conveniently formulated in terms of a likelihood ratio: the likelihood of an interpretation involving the constituents composed into a single entity divided by the likelihood of the same interpretation in which the constituents are related only by chance. The entire distri- bution on image interpretations can be stitched together from the specification of all such likelihood ratios. These likelihood ratios, and their dependence on the poses of constituent parts, raise the important question of coordinate systems. A compelling argument can be made for using only relational coordinate systems, so as to define an object independent of it’s particular presentation in three-dimensional space. But the formulation of a suitable relational coordinate system is not straightforward, due mostly to the nonlinearities introduced by scale. 2. Applications in medical imaging. In medical imaging, the term atlas usually refers to a (probabilistic) model of a population of images, with the parameters learned from a training dataset. In its simplest form, an atlas is a mean intensity image. Atlases are used for various pur- poses including normalization of new subjects for structure and function locations, segmentation or parcellation of certain structures of interest and group analyses to identify pathology related changes or developmental trends. We will explore the applications of our maximum likelihood templates in these areas, especially in the 132 study of MRI images. The learned templates will play the role of atlases. Since MRI images are 3D, we need to generalize our model from 2D to 3D, and the unit concept from “pixel” to “voxel”. Since our model is not constrained by the dimensionality of the image data, this generalization is straightforward. Another issue that we need to address to is the specialty of medical images compared to natural images. Minor local differences between two MRI images can lead to opposite diagnosis results. Currently our models only consider the rigid global transformation of templates, e.g. scales and rotations, and this is not enough for modeling medical images. We have to consider local deformation in addition to the current global transformation, e.g. 3D nonlinear transformation model parameterized via B-splines [50]. 3. Video Analysis. One advantage for using generative (Bayesian) models is the availability of a mathematically coherent extension to multi-source data. Thus, there is nothing in principle that confines a compositional model to a single frame. This genera- tive framework can be adapted to solve, for example, video-related tasks. Given a compositional hierarchy for a face or body, and a “state equation” for a trajec- tory of parts, the generative form of the model generates image sequences (video) rather than single images. The “inverse” problem, the problem of identifying high- likelihood interpretations under the posterior distribution, becomes a version of the tracking problem. Bibliography [1] V. N. Vapnik. The nature of Statistical Learning Theory Springer-Verlag, Berlin, 1995. [2] C. M. Bishop. Neural networks for pattern recognition. Oxford University Press, Oxford, UK, 1995. [3] Y. Freund and R. E. Schapire. A decision-theoretic generalization of on-line learning and an application to boosting. Journal of Computer and System Sci- ences, 55(1): 119-139, 1997. [4] F. Fleuret and D. Geman. Coarse-to-fine face detection. Inter. Journal of Com- puter Vision, 41, 85-107, 2001. [5] U. Grenander. General pattern theory: a study of regular structures. Oxford University Press, 1993. [6] S. Geman, K. Manbeck, and E. McClure. Coarse-to-fine search and rank-sum statistics in object recognition. Technical report, Division of Applied Mathemat- ics, Brown University, 1995. [7] S. Geman and M. Johnson. Article title Probability and statistics in compu- tational linguistics, a brief review . Mathematical foundations of speech and language processing. M. Johnson; S. P. Khudanpur; M. Ostendorf; R. Rosenfeld (Eds.), X, ISBN: 0-387-20326-5, pp. 1-26, 2004. 133 134 [8] Y. Amit, D. Geman, and X. Fan. A coarse-to-fine strategy for multi-class shape detection. IEEE Trans. PAMI, 2004. [9] B. Ommer and J. M. Buhmann. Learning the compositional nature of visual objects. In CVPR’07, IEEE, 2007. [10] R. Fergus, P. Perona, and A. Zisserman. Object Class Recognition by Unsuper- vised Scale-Invariant Learning. CVPR’03, IEEE, Vol. 2 (2003), pp. 264-271, 2003. [11] E. Bienenstock, S. Geman, and D. Potter. Compositionality, MDL priors, and object recognition. M.C. Mozer, M.I. Jordan, and T. Petsche, editors, Advances in Neural Information Processing Systems, Vol. 9, pp. 838-844, MIT Press, 1997 [12] S. H. Huang. Compositional approach to recognition using multi-scale computa- tions. Ph.D. thesis, Division of Applied Mathematics, Brown University, 2001. [13] M. Harrison Discovering compositional structure. Ph.D. thesis, Division of Ap- plied Mathematics, Brown University, 2005. [14] S. Geman, D. F. Potter, and Z. Chi. Composition Systems. Quarterly of Applied mathematics, LX: 707-736, 2002. [15] Ya Jin. Non-Markovian Hierarchy Vision System. PhD thesis, Brown University, Division of Applied Mathematics, 2006. [16] Y. Jin and S. Geman. Context and hierarchy in a probabilistic image model. In CVPR’06, volume (2), pp. 2145-2152. IEEE, 2006. [17] W. T. Freeman, J. Yedidia, and Y. Weiss. Understanding belief propagation and its generalizations. International Joint Conference on Artificial Intelligence, 2001. 135 [18] M. Riesenhuber and T. Poggio. Models of object recognition. Nature Neuro- science, supplement, volume 3: 1199-1204, 2000. [19] Y. Amit and A. Trouve. Pop: Patchwork of parts models for object recognition. Technical report, University of Chicago, 2004. [20] M. Weber, M. Welling, and P. Perona. Unsupervised learning of models for recognition. Proc. Sixth European Conf. computer Vision, pp. 18-32, 2000. [21] Z.J. Xu, H. Chen and S.C. Zhu. A high resolution grammatical model for face representation and sketching. Proc. IEEE Conf. on Computer Vision and Pat- tern Recognition (CVPR), San Diego, June 2005.[pdf] [22] F. Min, J.L. Suo, and S.C. Zhu. An and-or graph model for face representa- tion, sketching and aging. Chapter in Encyclopedia of Biometric Recognition , Springer, 2009. [23] L. Zhu, Y. Chen, X. Ye, and A. Yuille. Structure-Perceptron Learning of a Hierarchical Log-Linear Model. In CVPR, 2008. [24] David M. Green & John A. Swets (1966). Signal Detection Theory and Psy- chophysics. New York: Wiley. [25] Thomas M. Cover and Joy A. Thomas (1991). Elements of Information Theory, pp. 300-301. New York: Wiley. [26] E. Sali and S. Ullman. Combining class-specific fragments for object classifi- cation. In proc. 10th British Machine Vision Conference, volume 1, 203-213, 1999. [27] S. Ullman, E. Sali, and M. Vidal-Niquet. A fragment-based approach to object representation and classification. In IWVF, pp. 85-100, 2001. 136 [28] E. Borenstein and S. Ullman. Class-specific, top-down segmentation, In ECCV02, LNCS 2353, pp. 109-122, 2002. [29] B. Epshtein and S. Ullman. Feature hierarchies for object classification. In ICCV, 2005. [30] B. Heisele, T. Serre, M. Pontil, T. Vetter, and T. Poggio. Categorization by learning and combining object parts. In NIPS, 2001. [31] B. Heisele, T. Serre, and T. Poggio. A component-based framework for face detection and identification. Int. J. of Comp. Vision 74(2), pp. 167-181, 2007. [32] T. Serre, A. Oliva, and T. Poggio. A feedforward architecture accounts for rapid categorization. Proceedings of the National Academy of Science, 104(15), pp. 6424-6429, April 2007. [33] D. H. Hubel and T. N. Wiesel. Receptive fields and functional architecture of monkey striate cortex. J. Physiol., 195, pp. 215-243, 1968. [34] T. Serre, A. Oliva, T. Poggio. A feedforward architecture accounts for rapid categorization. Proceedings of the National Academy of Sciences, vol. 104, no. 15, 2007. [35] S. Geman. Invariance and selectivity in the ventral visual pathway. Journal of Physiology-Paris, Vol. 100, No. 4., pp. 212-224, October 2006. [36] R. Fergus, P. Perona, and A. Zisserman. Object class recognition by unsupervised scale-invariant learning. In CVPR 2003. [37] B. Leibe and B. Schiele. Interleaved object categorization and segmentation. Pro- ceedings of British Machine Vision Conference (BMVC), 2003. [38] S. Agarwal, A. Awan, and D. Roth. Learning to detect objects in images via a sparse part-based representation. IEEE TPAMI, 26(11): 1475-1490, 2004. 137 [39] S. Allassonni`ere, Y. Amit, and A. Trouv´e. Towards a coherent statistical frame- work for dense deformable template estimation. Journal Of The Royal Statistical Society Series B, 2007, vol. 69, issue 1, pp. 3-29. [40] M. R. Sabuncu, S. K. Balci, and P. Golland. Discovering Modes of an Image Population through Mixture Modeling. Proceedings of the International Confer- ence on Medical Image Computing and Computer Assisted Intervention (MIC- CAI), LNCS 5242, p. 381-389, 2008. [41] Z. Si, H. Gong, S.C. Zhu, and Y.N. Wu. Learning Active Basis Models by EM- Type Algorithms. Statistical Science, 2009. [42] R. B. Blackman and J. W. Tukey. The measurement of power spectra. Dover Publications, New York, 1958. [43] P. Viola and M. Jones. Rapid object detection using a boosted cascade of sim- ple features. In IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2001. [44] B. A. Olshausen and D. J. Field. Sparse coding with an overcomplete basis set: A strategy employed by V1?. Vision Research, 37(23): 3311-3325, Dec. 1997. [45] G.E. Hinton. Products of experts. In Int. Conf. on Art. Neur. Netw. (ICANN), vol. 1, pp. 1-6, Sept. 1999. [46] M. Welling, G.E. Hinton, and S. Osindero. Learning Sparse Topographic Repre- sentationswith Products of Student-t Distributions. In Adv. in Neur. Inf. Proc. Sys. (NIPS), vol 15, pp. 1359-1366, 2003. [47] S. Roth and M. J. Black. Fields of experts: A framework for learning image priors. In IEEE Conf. on Computer Vision and Pattern Recognition, vol. II, pp. 860-867, June 2005. 138 [48] S. Roth and M. J. Black. Fields of experts. In International Journal of Computer Vision, 82(2): 205-229, Apr. 2009. [49] S. C. Zhu and D. Mumford. Prior learning and Gibbs reaction-diffusion. IEEE Trans. Patter Anal. Mach. Intell., 19(11): 1236-1250, Nov. 1997. [50] D. Rueckert, et al.. Nonrigid registration using free-form deformations: appli- cation to breast MR images. IEEE Transaction of Medical Imaging 18, 712-721, 1999.