Image Compression and Data Clustering: New Takes on Some Old Problems by Wei-Ying Wang B.A., National Taiwan University; Taiwan, 2004 M.Sc., National Taiwan University; Taiwan, 2006 M.Sc., Brown University; Providence, RI, 2012 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 2017 c Copyright 2017 by Wei-Ying Wang This dissertation by Wei-Ying Wang 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 Matthew Harrison, Ph.D., Reader Date Elie Bienenstock, Ph.D., Reader Approved by the Graduate Council Date Andrew G. Campbell, Dean of the Graduate School iii Vita In 2000, Wei-Ying Wang started his undergraduate study at National Taiwan Univer- sity, and he received his Bachelor of Arts degree in Economics in 2004. After that he entered the graduate school at National Taiwan University, and received his Master of Science degree in Mathematics in 2006. Later he joined the military until 2008. He then worked for image denoising and signal processing in Academia Sinica for two years. Wei-Ying Wang came to the United States in 2010, and has been attending the Ph.D. program in the Division of Applied Mathematics at Brown University. He received a Master of Science degree in Applied Mathematics in 2012. This dissertation was defended on May 11th, 2017. iv Acknowledgments I would like to thank all the people who have helped me to complete my thesis. Most of all, I want to express my deepest gratitude to Prof. Geman, who helps me on so many aspects. This thesis is only possible because of him. When he saw my compression algorithm, he related it to a classical theorem so that we have an image compression algorithm that has analytic guarantee. He also came up with the idea of summarizing data with a concave distance function, rendering a very robust estimator. There are a lot of other examples and I can’t tell it all here. I am very fortunate to be his student. I would also like to thank my other committee members, Prof. Harrison and Prof. Bienenstock, for carefully reading my dissertation and giving precious comments. I wish to acknowledge Prof. Chii-Ruey Hwang, Prof. Ting-Li Chen, Prof. Lo-Bin Change, Dr. Hong Zhang, Dr. Jackson Looper, and Wilson Mckerrow, who are extremely helpful during my research. I need to mention about my wife, Yi-Wen Lai, who faithfully walks all the way with me, and see me through highs and lows of my days. I am thankful to have my first son, Chris, for providing so much fun and meaning in my life to keep me going. Also, I owe my thanks to my parents in Taiwan, who have been very supportive and encouraging. Last but not least, thank God for putting me here, experiencing this much, and carrying me out during countless stay-up night. v Contents Acknowledgments v 1 Overview 1 2 Entropy Rate on Images 4 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 The Entropy of Images . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2.1 Coding Rate on a Stationary and Ergodic Source . . . . . . . 13 2.2.2 Entropy Rate on a Stationary Source . . . . . . . . . . . . . . 17 2.2.3 Ergodic Theory on Images . . . . . . . . . . . . . . . . . . . . 20 2.2.4 Markov Approximation on Images . . . . . . . . . . . . . . . . 47 2.2.5 Coding Rate on a Stationary Source . . . . . . . . . . . . . . 57 2.3 The Ideal Compression Algorithm . . . . . . . . . . . . . . . . . . . . 64 2.3.1 The L1 (P ) Convergence . . . . . . . . . . . . . . . . . . . . . 65 2.3.2 The P -a.e. Convergence . . . . . . . . . . . . . . . . . . . . . 67 2.3.3 How to Build an Ideal Compression Algorithm . . . . . . . . . 74 2.4 Application . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 2.4.1 Comparison-Based Image Compression Algorithm . . . . . . . 76 vi 2.4.2 The Performance . . . . . . . . . . . . . . . . . . . . . . . . . 85 3 Robust Generalized Clustering 91 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 3.2 Loss Function for Multiple Instances . . . . . . . . . . . . . . . . . . 100 3.2.1 Reduction to m = 1 . . . . . . . . . . . . . . . . . . . . . . . . 101 3.2.2 Lines in R2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 3.2.3 Hyperplanes in Rd . . . . . . . . . . . . . . . . . . . . . . . . 109 3.2.4 Non-interpolating Optima: Some Partial Results . . . . . . . . 118 3.3 Observations about Sparse-distance Clustering in R1 . . . . . . . . . 131 3.3.1 Asymptotic Properties . . . . . . . . . . . . . . . . . . . . . . 132 3.3.2 Breakdown Point . . . . . . . . . . . . . . . . . . . . . . . . . 136 3.4 Minimizing the Loss Function . . . . . . . . . . . . . . . . . . . . . . 143 3.4.1 Remove-or-Replace Algorithm . . . . . . . . . . . . . . . . . . 143 3.4.2 Robust Nature for Sparse Distance . . . . . . . . . . . . . . . 147 3.4.3 Instances from a Single Model . . . . . . . . . . . . . . . . . 150 3.4.4 Instances from Multiple Models . . . . . . . . . . . . . . . . . 161 3.5 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171 4 Conclusions 173 4.1 On Entropy Rate of Images . . . . . . . . . . . . . . . . . . . . . . . 174 4.2 On Robust Generalized Clustering . . . . . . . . . . . . . . . . . . . . 174 vii List of Figures 2.1.1 A typical context for 3 × 5 black-and-white image patch, where y ∈ {0, 1, · · · , 255} is the value of a target pixel, and its context value is represented as (x1 , x2 , · · · , x12 ) ∈ {0, 1, · · · , 255}12 . . . . . . . . . . . 5 2.2.1 Λn on Z2 plane, note that the coordinator we use is different from Cartesian coordinator system. . . . . . . . . . . . . . . . . . . . . . . 11 0 2.2.2 Bm on Z2 plane. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2.3 The plot of Znm on Z2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.2.4 Example of rectangle sets. . . . . . . . . . . . . . . . . . . . . . . . . . 48 2.4.1 The kernel we use to approximate empirical distribution of P (X0 |XBm ). 77 2.4.2 A example of α-weight corresponding to different tier and sizes of patch. Note that the weight is not yet been normalized and is shown propor- tionally. For example: (a) For a 3 × 5 patch, T1 has a weight proportion to 3, T2 has a weight proportion to 2; T3 has a weight proportion to 1, and B = 3 in this case. Note that a notation × indicates the target pixel location for each patch. . . . . . . . . . . . . . . . . . . . . . . . 78 2.4.3 Demonstration of CBIC algorithm. . . . . . . . . . . . . . . . . . . . . 81 viii 2.4.4 Raster scan as decoding the image from an algorithm built on context model with maximum patchsize 2 × 3. The procedure start from (a): Given the first 3 top-left corner pixel values, using a context model on 1 × 2 patch to get the first row of image (b). Next, apply a context model on 1 × 2 patch to get left and right column of image. Then use that of 2 × 3 to get the second row, and continue to obtain the whole image (e). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 2.4.5 4 test images: Lena, Jetplane, Goldhill, and Barbara. Goldhill has size 720 × 576, and others are 512 × 512. . . . . . . . . . . . . . . . . . . . 86 2.4.6 The bpp result of CBIC algorithm for “Lena” image with different size of patches and different size of library. In these experiments, we choose essential sample amount K = 200 to build the empirical distribution. 87 2.4.7 Comparison between CBIC and other famous lossless compression algo- rithms. The following information in the parenthesis shows the number of library and the patchsize used to generate the result here: Lena: (80M, 5 × 9); Jetplane: (80M, 4 × 7); Goldhill: (40M, 4 × 7); Barbara: (20M, 4 × 7). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 2.4.8 Top: Reconstruct the image with MAP estimators. Bottom: The orig- inal picture. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 3.1.1 “Triangle data” . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 Pn 3.1.2 Minimizers of f (`) = i=1 d(zi , `)β when β = 1 and β = 2, where ` is a straight line in R2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 3.1.3 Left: the data includes 70 points generated from the uniform distribu- Pn tion. Right: The effect of β on the minimizer of f (`) = i=1 d(zi , `)β , over all lines, `, in R2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 ix 3.2.1 Counter example for the optimal line in R3 , which explaining the op- timal line might not necessary go through two data points. The data is z1 , · · · , z4 , when b is high enough, the green line will be the best line (that has smaller or equal minimum target function value) among all the interpolating lines, i.e. the lines that connect two data points. However, when β < 1 and is close to 1, the red line will produce smaller target function value than that of the green one, and the red line only touches one data point. . . . . . . . . . . . . . . . . . . . . . . . . . . 110 3.4.1 We examine the minimizer of Function 3.4.4 for L ⊂ M2 , a set of in- terpolating lines in R2 on X-shaped data with outliers. (a) The original data consists of 100 data points that corrupted by 40 outliers. (b) EM algorithm for 2 instances. (c) RRA result for β = 0.8 and λ = 30 with random initial. In this case, the exact solution coincides RRA result. . 148 3.4.2 Different λ setting to RRA when minimizing Function 3.4.4 for L ⊂ M2 , a set of interpolating lines in R2 , on X-shaped data with outliers. One can see that for λ = 10 ∼ 50, the algorithm gives two lines. Even when other λ is chosen, the instances that RRA found is still informative. Note that here we use random initial for RRA. . . . . . . . . . . . . . 149 3.4.3 The "S-curve" data consists of 50 data points (n = 50). . . . . . . . . 151 3.4.4 The exact solution (plotted as set of lines) of minimizing the loss func- tion (Equation 3.4.5) while we use β = 0.8. (a) When the amount of instances m = 2 and (b) m = 3. . . . . . . . . . . . . . . . . . . . . . 152 x 3.4.5 The result of RRA with locally best initial on S-curve data (n = 50) for different λ, where β = 0.8. When λ = 20 ∼ 110 (the second and the third plots), RRA gives 3 line instances. Note that the sum of error, h, for m = 2 is 18.711 (the first plots) and for m = 3 is 9.094 (the second and the third plot) coincide with that of exact solution, means that RRA approximates the exact solution well in this data. . . . . . . 152 3.4.6 The result of RRA (β = 0.8, λ = 110) with random initial on S-shaped data with different random seed. From left to right, the random seed are set as 0,1,2,3,4,5, respectively. One can see that initialization does not cause huge different on RRA result in this sample. . . . . . . . . 153 3.4.7 Some results of various data in R2 . RRA is initialized with locally best initial. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 3.4.8 “Line 3D” data for the experiment on the interpolating line model on R3 . (a) and (b) are different viewing angle. . . . . . . . . . . . . . . . 156 3.4.9 The exact solution for “Line 3D” data. (a) The exact solution for m = 2 with different viewing angles on top and bottom. (b) The exact solution for m = 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156 3.4.10 Result of “Line 3D” on R3 data. Top row and bottom row are different viewing angle of the result. From left to right is the result of different λ setting. A desired results are occurred when λ ∈ (30, 250). Note that RRA coincides with the exact solution. . . . . . . . . . . . . . . . . . 157 3.4.11 The“Plane 3D” data set for the experiment on the plane model on R3 . (a) and (b) are different viewing angle. . . . . . . . . . . . . . . . . . 159 3.4.12 The exact solution for m = 2 for “Plane 3D” data. (a) and (b) are different viewing angle. . . . . . . . . . . . . . . . . . . . . . . . . . . 159 xi 3.4.13 RRA results of “Plane 3D” on R3 data. Top row and bottom row are different viewing angle of the results. From left to right are the results of different λ setting. A desired results are occurred when λ ∈ (10, 50). Note RRA coincides with the exact solution in this case. . . . . . . . . 160 3.4.14 The 2D data “percent” and “Corrupted Percent” for experiment in this section. (a) “Percent” data consists of n = 90 data points. (b) “Cor- rupted Percent” data has 40 uniformly distributed noises as outliers in “percent” data (30.8% of the data). . . . . . . . . . . . . . . . . . . . 163 3.4.15 Results of RRA with the model penalty parameter [ω1 , ω2 ] on the “Per- cent” data. (n = 90). While setting ω1 = 1, a satisfactory result is appeared when ω2 = 2 ∼ 7 (Plots (b) and (c)). The parameters are denoted on top of each plots. One can see the penalty term acts on RRA result. The big dot markers are the instances in M1 found by RRA; the line is the instance in M2 found by RRA. . . . . . . . . . . 164 3.4.16 Results of RRA with the model penalty parameter [ω1 , ω2 ] on “Cor- rupted Percent” data (n = 130). While setting ω1 = 1, a satisfactory result is appeared when ω2 = 1.5 ∼ 6 (Plots (b) and (c)). The parame- ters are denoted on top of each plots. One can see the penalty term acts on RRA result. The big dot markers are the instances in M1 found by RRA; the line is the instance in M2 found by RRA. . . . . . . . . . . 164 3.4.17 (a) “Corrupted percent” data, where 40 out of 130 data points are out- liers (30.8% of the data). (b) Candidates for MSV method on percent data, required from RRA at step 1 of MSV method. There are total 6 candidates, represented as bigger dots and lines, appeared in this data 168 xii 3.4.18 Results of MSV method on “Corrupted percent” data. Assuming u = 70, that is, roughly speaking, we use 70% of the data in a cluster to calculate the volume. (a) Assume the amount of instances m = 2. (b) m = 3. (c) m = 4. Note that the value of M SV (m) is denoted at the top of each plot. When m = 4, there is no combination that provides the sufficient amount (> ) of data in a cluster, so M M SV (4) = ∞ (on the top of the plot, the MSV value is 3053.3, which is just a big value that indicate this happened). The transperant line is indicated as the “volume” of each instances. . . . . . . . . . . . . . . . . . . . . . . . . 168 3.4.19 The “Corrupted 3D mixture” data, where 15 out of 65 data points are outliers (23.1% of the data). (a) The original data from different viewing angle (top and bottom). (b) Candidates for MSV method on percentile data, required from RRA, from different viewing angle (top and bottom). There are total 9 candidates, represented as bigger dots, lines, and planes, appeared in this data. . . . . . . . . . . . . . . . . 169 3.4.20 Results of MSV method on corrupted 3D mixture data. We provide different viewing angles by separating them at top row and bottom row. Assuming u = 70, i.e. 70% of data in a cluster is used to compute the volume, roughly speaking. (a) Assume the amount of instances m = 2. (b) m = 3. (c) m = 4. Note that the value of M SV (m) is denoted at the top of each plot. . . . . . . . . . . . . . . . . . . . . . 170 3.4.21 m = 3 result for “corrupted 3D mixture“ data. The transparent area is the volume we calculated for each instance. . . . . . . . . . . . . . . . 170 xiii Chapter One Overview 2 An important job of a statistician is to summarize data. To this end, many statistical methods were developed in the 18th century, when mathematicians like Laplace and Gaussian introduced the least square method and the least absolute value regression. These “classical” statistical tools become the foundation of modern statistics, even though people was only able to deal with a very small amount of data. Then, in the early 20th century, computers were invented and these once considered “impractical” methods finally shined. However, as computer technology advanced, it brought more challenges. For example, as data being digitized, we started to consider the most efficient way to store and transmit them. In 1940s, Claude Shannon published topics about data compression and entropy. The idea of compression is to use probabilities to “filter out” the unnecessary information in the data. Just like many statistical tools—it summarizes the data. A related field is called the “lossless compression,” which compresses the data in a way that people can even reconstruct the original data from the summarized information. Another example is a field called data clustering, which emerged with the increasing amount of data. In 1950s, K-means algorithm was invented to cluster, or say, to summarize, data points into several representative points. About the same time, people started to realize that some problems cannot be solved with polynomial computational time, like the problem of minimizing the loss function of K-means. These are called NP-hard problems, where it might need another technology revolution to deal with. We look into these two old topics, data compression and data clustering, in this thesis, and gives some new ideas upon it. In Chapter 2, we will examine theories behind lossless image compression algo- rithms. It is surprising that people didn’t explore the theories of image compression a lot. There were many literature in information theory, but only a few of them were applied for describing information on images. In this chapter we will try to be as 3 self-contained as possible so that people who have some knowledge in probability the- ory should be able to appreciate it. We will expand existing theories for information on images, and come up with an optimal compression algorithm from it. Just like those statisticians invented the classical methods far before computers, the algorithm is very impractical in view of modern computers. However, during experiments we see that it triumphs a state of the art image compression scheme and we can further push our algorithm more to achieve the theoretical optimal compression rate of an image. In Chapter 3, we generalize the loss function of K-means to hyperplanes, where it was defined on points. The problem is still NP-hard and the original K-means algorithm is inapt for this generalization, since it gets stuck on local optima easily. The data concerned require rotation equivariant methods to summarize it, rendering most regression methods useless. We utilize a concave distance function, called the sparse distance, who helps us gain huge robustness. One will see the effect of using such a distance, where it successfully summarizes the data with high amount of noise and outliers. Moreover, the sparse distance makes the solution of our loss function “interpolating”—an important properties to allow us to solve the minimizer. While it is still difficult to solve if there are many inherent clusters, we designed a well-suited algorithm that can approximate the exact solution very well. Chapter Two Entropy Rate on Images 5 Figure 2.1.1: A typical context for 3 × 5 black-and-white image patch, where y ∈ {0, 1, · · · , 255} is the value of a target pixel, and its context value is represented as (x1 , x2 , · · · , x12 ) ∈ {0, 1, · · · , 255}12 . 2.1 Introduction After the arithmetic coding was invented, most research on lossless compression al- gorithms focus on building models. Most of the image compression algorithms, and most successful ones, fall into the category of predictive context modeling, which build a statistical model on the value of a target pixel given the surrounding pixels, or say, the context. A better prediction of target pixels from their context rendering a better compression rate. One of the hardest problem on image compression is the unman- ageable number of contexts. For example, in a typical context in a small 3 × 5 image patches as shown in Figure 2.1.1 , where each pixel has 256 depth of intensities will result in 25612 of possible contexts. This size grow fast when one consider a larger sized context. Of course, most of them are not seen in images, and lots of them are similar to each other. To find suitable and representative contexts is a generally done by data clustering. Here are some successful algorithms that use predictive context modeling. A Famous algorithm LOCO-I/Jpeg-LS[25] uses a very simple predictor with different residual models on several clusters of contexts. The advantage of it is the minimum 6 memory and computation resources requirement. CALIC [41] algorithm uses a simple predictor and ideas such as the error feedback to better classifying contexts. It’s a 1 successful image compression and usually has 0.5 to 1 less bits per pixel than that of LOCO. Practically, most of the time CALIC is among the best when comparing different compression schemes. The reason it is not as widely-used as LOCO-I is the requirement of relatively expensive computation resources. There are some compres- sion scheme that is comparable to the results of CALIC. For example, TMW [29] introduces multiple linear predictors and utilize the idea of blending on different class of context. Another example is Golchin et. al. [15] who build several linear predictors and assigns blocks of images into its best predictor. While these algorithm uses local information, there are also attempts of using global attributes of image, [35] using segmentation method to reduce the number of the contexts to area like edges and smooth textures in the image. These algorithms achieve some success on practical end, and often time the best and the second best have marginal differences of compression rates. However, theo- retical supports are rarely found. One of the algorithm that has theoretical backup is called the universal modeling, which is a variation of the Ziv-Lemple coding. The core of this algorithm is to exploit repetitions of the context inside of the source, and build an empirical distribution on that. The application of the universal modeling on image compression can be found in [26]. Due to the nature of huge amount of contexts in images, its not easy to find similar patches on sources with small size. The result is often poor and cannot compete to some of the best compression schemes, like CALIC. 1 the average number of bits required to code an image 7 The motivation of this article is to get some idea about how much room was left for improvement in lossless compression. Thus, we look for the notion of “minimum coding rate” in the literature, the most related article is written by Ornstein and Weiss [34], where it extends the Shannon-McMillian-Breimen theorem on higher dimensions. It indicates, under stationary and ergodic condition, the minimum coding rate (or 2 say, bits per pixel) for an image will converge to a constant, which happens to be the entropy rate. However, it is hard to believe that all images have the same theoretical minimum bits per pixel. We often see the compression results differ from one another. For example, the image with higher amount of noise will require higher amount of bits to code it, since there is more “surprises” in the image. So what is going wrong with the theory and the real world? Apparently, the theory is valid and the stationary condition is very nature. The inconsistency, we believe, is due to the ergodic assumption. In Section 2.2 we will eventually drop the ergodic condition, and the minimum coding rate, under the stationary condition, converged3 to a random variable depends on the image it self. Furthermore, the random variable can be characterized with an ergodic components of the underlying stationary probability. So, we have the idea of the theoretical bound for every image; the next step to to understand how the state of the art lossless image compression algorithm performed. Specifically, we aim at predicative context modeling. In Section 2.3 we give a theoret- ical result for predictive context modeling method. With only stationary condition and some small assumptions, the coding rate, when the context size grows, will con- verge to the minimum bits per pixel for every image, as we mentioned in the previous 2 The “image” refers to an infinite sized image. Mathematically speaking, an image is considered as a finite alphabets sequence with Z2 indices. 3 The convergence is in the sense of “almost everywhere,” means it converges for almost all images except those have probability 0. 8 paragraph. In the end of this section, we will mention about an idea to utilize the theorem we got, to build an algorithm that can be “push to the limit,” i.e. achieve the minimum coding rate. The algorithm we built to test the idea is described in Section 2.4, where we call it the comparison-based image compression algorithm (CBIC). The key idea is to use an empirical distribution from external source. With the property that a empirical distribution will converge to the true one with larger amount external library4 and larger context, we can have the proposed algorithm converged to the true entropy, or say, the minimum bits per pixel. Through experiments we see that it is able to defeat most of the state of the art algorithms, as we would expect. Note that CBIC is not a practical compression algorithm, since it requires lots of memory usage and expensive computation resources. However, the ability to scale up and provably optimal make it possible to examine the best algorithm we have so far. We think that our result suggests that there is very little room left for compression algorithms, since we only do a tiny bit better than theirs and the fact that the gain becomes marginal when we scale up the context size. This article is arranged as follows: Section 2.2 aims on the theory behind the min- imum coding rate of an image. In the end of the section, we proved that, under the stationary condition, the minimum coding rate of an image converge to the entropy rate of the ergodic component with respect to the image. Section 2.3 gives theoretical 4 The use of the external library for compression can be found in the lossy compression schemes. Algorithms like VQ, the vector quantization [30], codes blocks of images into indices of library. Their work includes reducing the redundancy and dependency of the library so as to minimize the size of the library to gain a better compression ratio. However, in our algorithm, the redundancy of the library is one of the important elements, which forms a nice statistical property for empirical distributions. 9 background of image compression algorithms with predictive context modeling. Sec- tion 2.4 introduce CBIC, which is able to scale up the context size so as to achieve the optimal coding rate. It helps examine the performance of some existing algorithms like JPEG-LS and CALIC. 10 2.2 The Entropy of Images In this section, we would like to give a theoretical background of the optimal com- pression rate. There are abundant theories on the entropy rate on 1D sequence space, see [17, 16]. To extend them to images, i.e. 2D lattice, is rarely seen in the literature. Although they have many similarities, it is not clear that we can directly extend the 1D result. We aimed to fill this void and give a better understanding of the the- ory behind image compression algorithms. In this section, most of theories will be extended from that of 1D, follows the idea of Gray’s book [16, 17]. Before we describe the theoretical result, we will give several definitions. We first define a probability space on Z2 -lattice, say (Ω, F, P ). Let X be the finite space of pixel values, for example, X = {0, 1, · · · , 255}, which are possible intensities in 2 ` a pixel. Define Ω , X Z and F , σ({χΛn }n∈N∪0,`∈Z2 ), a σ-algebra generated by ` generating set C , {χΛn }n∈N∪0,`∈Z2 , where Λn , {(i, j) ∈ Z2 | − n ≤ i, j ≤ n} and Λ`n , Λn + `, which is an (2n + 1) × (2n + 1) square lattices of Z2 centered at `, for ` ∈ Z2 . Λn is shown in Figure 2.2.1. Note that Λn ⊂ Λn+1 with ∪n Λn = Z2 . We then define P to be a probability measure on (Ω, F). For any ω ∈ Ω, define the the transformation, T ` : Ω → Ω for ` ∈ Z2 as (T ` ω)i = ω`+i for ω ∈ Ω, `, i ∈ Z2 . Also, we define X0 (ω) = X(ω) , ω0 , X` (ω) , X(T ` ω) = ω` , and 11 XA , ωA for any ` ∈ Z2 and A ⊂ Z2 . We call that (Ω, F, P, T ) a dynamical system, where T = {T ` }`∈Z2 . We will use the notation (< b) , {a ∈ Z2 |a < b} for b ∈ Z2 , where “<” is the lexicographic order on Z2 . That is, a < b if and only if either a1 < b1 or a1 = b1 , a2 < b2 . The entropy of a random variable Y is defined as HP (Y ) , EP (− log P (Y )), and the entropy rate is defined as: 1 hP , lim HP (XΛn ) (2.2.1) n→∞ |Λn | with |A| being the number of elements in A for A ⊂ Ω. Here we have |Λn | = (2n+1)2 . We sometimes omit subscript of h and H if the underlying probability is obvious. The logarithms throughout the article are based on 2. Figure 2.2.1: Λn on Z2 plane, note that the coordinator we use is different from Cartesian coordinator system. 12 Another similar definition is the ideal coding rate, which describes the coding rate when image size Λn is large: − log P (XΛn (ω)) lim , (2.2.2) n→∞ |Λn | given ω ∈ Ω. Note that we haven’t established whether the entropy rate or the coding rate are existed yet. The necessary condition for their existence will be proved in Subsection 2.2.1. Subsection 2.2.1 describes the theorem saying that the coding rate is actually the entropy rate under stationary and ergodic conditions. Subsection 2.2.2 depicts the entropy rate under stationary condition, in particular, we can view the entropy rate as the conditional entropy given infinitely many history, which will be described in Corollary 2.2.4. Subsection 2.2.3 gives a solid background (which is hard to find in literature) of the ergodic decomposition for 2D lattice, which is essential after we drop the ergodic condition. Subsection 2.2.4 defines the Markov approximation on 2D, and discusses some related properties of it. The Markov approximation, P m for m ∈ N, will have the following property P m (XΛn ) = Y P (X` |XBm ` ∩Λ ) n `∈Λn ` for Bm ⊂ Ω, is the mth order history for location ` ∈ Z2 , depicted in Figure 2.2.2. The Markov approximation is an important idea in practical compression scheme where it uses local information approaching the global one. Note that the coding rate (Equation 2.2.2) is defined on the global information. Subsection 2.2.5 proves the the Theorem 2.2.34, which generalizes Shannon-McMillian-Breimen theorem on 13 0 Figure 2.2.2: Bm on Z2 plane. 2D non-ergodic source. Briefly speaking, given an ω ∈ Ω the coding rate with proba- bility P will be the entropy rate of P ’s ergodic components with respect to ω, which corresponds to the lowest possible compression rate of an image. 2.2.1 Coding Rate on a Stationary and Ergodic Source First let’s introduce the following assumptions: 1. (Stationary condition) A dynamical system (Ω0 , F 0 , µ, S) is stationary if the underlying probability measure, µ, is stationary. i.e. µ(F ) = µ(S −1 F ) for F ∈ F 0 and S ∈ S . 2. (Ergodic condition) A dynamical system (Ω0 , F 0 , µ, S) is ergodic if µ(A) ∈ {0, 1} for A ∈ IS , the class of S-invariant sets5 . i.e. IS = {A ∈ F|S −1 A = A} for S ∈ S. Remark. The reader might wondering why we use S −1 in the definition instead of simply S. This is because, for example, in space (Ω0 , F 0 ), given F ∈ F 0 and S is a transformation, SF = {ω : ω ∈ SF } = {ω : S −1 ω ∈ F }, and S −1 ω might not 5 The class of S invariance sets is actually a σ-algebra. See Ch. 6 of [13]. 14 be well-defined when Ω0 = χ{0,1,2,··· } and S is the left shift: (Sω)i = ωi+1 . However, if Ω0 = χZ then S −1 ω is well-defined and the stationary condition can be written 2 as P (F ) = P (SF ) in this case. Of course, in our setting Ω = χZ we can use P (F ) = P (T F ) for all T ∈ T and F ∈ F as the definition of P being stationary without trouble. Note that we sometimes say a sequence of random variable is stationary or ergodic means the underlying dynamical system is. If the dynamical system is 1D (i.e. when Ω0 = X Z ), and has these two conditions, both entropy rate and coding rate are existed and they agree to each other. This result is well-known in 1D setup, which was originally proved by Breiman [6]. Fortunately, 2 with similar arguments, we can establish them in 2D (i.e. Ω = X Z as our original definition of Ω), which will be described in the rest of this subsection. For the existence of coding rate in 2D, an important result is done by Ornstein and Weiss [34] and Lindenstrauss [24], who proved the Shannon-McMillan-Breiman theorem (SMB) for d-dimensional lattice for stationary ergodic sequences (another more comprehensive and easier to understand reference is [37]). More exactly, if we assume {Xi }i∈Z2 is an ergodic and stationary sequence on (Ω, F, P, T ), then for almost every ω, we have the coding rate − log P (XΛn ) lim exist and constant a.e. (2.2.3) n→∞ |Λn | 15 In other words, 1 X lim − log P (X` |X(<`)∩Λn ) exist and constant a.e. (2.2.4) n→∞ (2n + 1)2 `∈Λn One can see that the coding rate is actually the compression rate for infinite-sized image using the correct underlying probability, which is assumed to be stationary and ergodic. In other words, the coding rate is the best compression rate we can achieve for an infinite-sized image from a stationary ergodic source. Note that we will later show that the convergence is actually the limit value of the entropy rate, hP . This conclusion can be shown with uniformly integrability. The next lemma, which is modified from Billingsley [3] who proved it in 1D, demonstrates this in our 2D setup. − log P (XΛn ) Lemma 2.2.1. Let gn = |Λn | , then {gn }n∈N is uniformly integrable. Proof. Let an , |Λn |. Given ω ∈ Ω, if gn (ω) ∈ [k, k + 1), we have P (ωΛn ) ∈ (2(−k−1)an , 2−kan ], so Z gn dP ≤ (k + 1)P (gn ∈ [k, k + 1)) gn ∈[k,k+1) ≤ (k + 1)P (ω : P (ωΛn ) ∈ (2(−k−1)an , 2−kan ]) ≤ (k + 1)2−kan · |χ|an ≤ (k + 1)2−kan +an log |χ| = (k + 1)2−an (k−log |χ|) . 16 Thus Z ∞ Z X gn dP = gn dP gn >k i=1 gn ∈[k+i,k+i+1) ∞ −(k+i−log |χ|)an X ≤ (k + i + 1)e , i=1 which converges to 0 as k → ∞. This leads to the following corollary. Corollary 2.2.2. Under the stationary and ergodic conditions, the entropy rate hP ,exist, and is the same as the coding rate a.e., i.e. 1 − log P (XΛn ) lim HP (XΛn ) = n→∞ lim n→∞ |Λn | |Λn | a.e. Proof. Let gn be defined as in Lemma 2.2.1. Since we have gn converges a.e. to a constant, say h (Equation 2.2.3) and gn ∈ L1 (P ) is uniformly integrable, the theory about uniform integrability tell us: h = E(h) = E( lim gn ) n→∞ = lim E(gn ) n→∞ 1 = lim E(− log P (XΛn )) n→∞ |Λ | n = lim H(XΛn ) n→∞ = hP . 17 2.2.2 Entropy Rate on a Stationary Source 1 Here we want to show that if the dynamic system is stationary, h , limn→∞ |Λn | H(XΛn ) exists and equals to H(X0 |X(<0) ), where H(X0 |X(<0) ) , n→∞ lim H(X0 |X(<0)∩Λn ). Note that it always exists since {H(X0 |X(<0)∩Λn )}n∈N is a decreasing non-negative sequence in R. We assume that {Xi }i∈Z2 is a stationary but not necessary a ergodic sequence. First, note that we can calculate H(XΛn ) as follows: X H(XΛn ) = H(X` |X(<`)∩Λn ) ). `∈Λn The following lemma prepare the way to prove hP = H(X0 |X(<0) ). Lemma 2.2.3. Under stationary condition, we have H(X` |X(<`)∩Λn ) ↓ H(X0 |X(<0) ) as n → ∞. Proof. Since{H(X0 |X(<0)∩Λn )}n∈N is a decreasing non-negative sequence in R, and H(X` |X(<`)∩Λn ) ↓ H(X` |X(<`) ) = H(X0 |X(<0) ) 18 for all ` ∈ Z2 , where the last equation is because of the stationary condition. Figure 2.2.3: The plot of Znm on Z2 . Corollary 2.2.4. Under the stationary assumption, the entropy rate, defined as 1 lim H(XΛn ) n→∞ |Λn | exists and equals H(X0 |X(<0) ). ˆ , H(X0 |X(<0) ). First note that it is easy to see that limn→∞ Proof. Let h 1 H(XΛn ) ≥ |Λn | ˆ since for all n ∈ N: h, 1 1 X H(XΛn ) = H(X` |X(<`)∩Λn ) ) |Λn | |Λn | `∈Λn 1 X ≥ H(X` |X(<`) ) |Λn | `∈Λn 1 X stationary H(X0 |X(<0) ) = |Λn | `∈Λn = ˆ h. So the rest is to prove limn→∞ 1 H(XΛn ) ˆ Because of Lemma 2.2.3, given an ≤ h. |Λn | 19 ˆ < . Set arbitrary  > 0, there is an M ∈ N s.t. for all m > M , H(X0 |X(<0)∩Λm ) − h Znm , {` ∈ Λn |Bm ` ` ⊂ Λn }, where Bm 0 = Bm 0 + `, Bm = (< 0) ∩ Λm for ` ∈ Z2 (Figure 0 2.2.2 and Figure 2.2.3 depict Bm and Znm , respectively). Calculate 1 1 X H(XΛn ) = H(X` |X(<`)∩Λn ) ) |Λn | |Λn | `∈Λn 1 X X = [ H(X` |X(<`)∩Λn ) ) + H(X` |X(<`)∩Λn ) )] |Λn | `∈Znm `∈Λn \Z m n 1 X X ≤ [ H(X` |XBm ` ) + H(X` |X(<`)∩Λn ) )] |Λn | `∈Znm `∈Λn \Z m n 1 X X = [ H(X0 |XBm 0 ) + H(X` |X(<`)∩Λn ) )] |Λn | `∈Znm `∈Λn \Z m n because of stationary condition. Note there |Znm | = (2n + 1 − m)(2n + 1 − 2m) and all the term in the summation are uniformly bounded. So 1 lim H(XΛn ) ≤ H(X0 |XBm 0 ) n→∞ |Λn | ˆ + . < h Since  is arbitrary, we have limn→∞ 1 H(XΛn ) ˆ ≤ h. |Λn | Thus, if we allow ergodic assumption, we can combine this with Corollary 2.2.2 to get: Corollary 2.2.5. Under the stationary and ergodic assumptions, the entropy rate, h, 1 − log P (XΛn ) exist and can be view as limn→∞ |Λn | H(XΛn ), H(X0 |X(<0) ) ,or limn→∞ |Λn | . Remark. Note that this is the Shannon-McMillian-Breiman theorem on 2D-lattice version. 20 2.2.3 Ergodic Theory on Images When looking at Corollary 2.2.5, it is nature to ask a question: What happens if one drops the ergodic assumption? To answer this, one has to understand the er- godic decomposition of a stationary probability. Briefly speaking, the entropy rate 1 limn→∞ |Λn | H(XΛn ) or H(X0 |X(<0) ) will still converge, due to Corollary 2.2.4, but the convergence will be the mixture of entropy of ergodic components. We will see that in − log P (XΛn (ω)) Theorem 2.2.20. And the coding rate limn→∞ |Λn | will still converge a.s., but the convergence will be a random variable depending on which ergodic component the ω belongs to. This part will be proved in Theorem 2.2.34. Before proving these the- orems, one needs to make sure that some famous ergodic theorems for 1D sequences is valid for generalizing to 2D images. Subsection 2.2.3.1 introduces Birkhoff ergodic theorem on 2D, which says 1 X lim f (T ` ω) = EP (f |I) P -a.e, n→∞ |Λ | n `∈Λn for I = IT is the T -invariant σ-algebra6 . In Subsection 2.2.3.2, we generalize the ergodic decomposition theorem to 2D, which basically says that every stationary distribution can be decomposed as mixture of ergodic stationary distributions. In the end of this section, Subsection 2.2.3.3 states that the entropy rate can be thought to be a mixture of entropies of ergodic components, which is relatively easy to generalized from 1D setup. 6 Remember that T = {T ` }`∈Z2 . 21 2.2.3.1 Birkhoff Ergodic Theorem on Images Let’s first introduce the pointwise ergodic theorem on Z2 lattice. It is well-known that in 1D sequence: For f ∈ L1 (µ) and µ is a stationary probability on 1D sequence with respect to a transformation S, we have 1 n−1 f (S i ω) = Eµ (f |IS ) µ − a.e., X lim n→∞ n i=0 where IS is an S-invariant σ-algebra (see [13]). A similar statement is valid in 2D, which can be found in [31]: for any f ∈ L1 (P ), 1 X lim f (T ` ω) = EP (f |I) P -a.e, (2.2.5) n→∞ |Λn | `∈Λn where I = IT is the T -invariant σ-algebra. i.e. the class of T invariant events, I , {B ∈ F|T −` B = B} for any ` ∈ Z2 . However, the explanation of this statement is hard to find (in [31], it is quoted from [33], where it might need more mathematical background to see it). The reader can assume this is right and proceed to the next subsection. For those who are interested about this, we would like to give a complete explanation to it. In the following content we assume a stationary dynamical system (Ω, F, P, T ) as we had on Z2 lattice defined in the beginning Section 2.2. First let’s introduce the special averaging sequence defined in [32]. Definition 2.2.6. (Special averaging sequence) {An }n∈N is a special averaging se- quence on group G = (Ω0 , ◦) if: (a) lim |(gAn )4An |/|An | = 0 for all g ∈ G, 22 (b) A1 ⊂ A2 ⊂ · · · , and (c) |A−1 n An | ≤ M · |An | for some constant M and for all n. The symbol “4” is the symmetric difference of two sets and AB , {a ◦ b : a ∈ A, b ∈ B}. Apparently, {Λn }n∈N is a special averaging sequence on group (Z2 , +). 0 Another example is {Bm }m∈N , defined as Λm ∩ (< 0), (we sometimes denoted as Bm ) which is shown in Figure 2.2.2. Note that the sequence that just run on one direction, e.g.. {(0, i) ∈ Z2 : i ≤ n}n∈N , is not a special averaging sequence on (Z2 , +) as it fails on (a) with g = (1, 0). With the fact that (Z2 , +) is an amenable group, we have the following theorem as shown in [32]. (Note that we assume P is a stationary distribution.) Theorem 2.2.7. If {An } is a special averaging sequence on (Z2 , +), then for any f ∈ L1 (P ), 1 X lim n→∞ |A | f (T ` ω) = f ? (ω) n `∈An exist for P-a.e. ω ∈ Ω and is a T -invariant function, i.e. for any T ∈ T , f ? (T ω) = f ? (ω) for P -a.e. ω ∈ Ω. We would like to characterize f ? by saying f ? = E(f |I), where I is the σ-field of T -invariant events. We will prove it with uniform integrable property in Theorem 2.2.10. Let’s first prepares it with following two lemmas. Lemma 2.2.8. For any f ∈ L1 (P ) and any transformation S = T ` , ` ∈ Z2 , we have E(f (S)) = E(f ) 23 Proof. We know that if f is an indicator function on Ω, say 1A , A ∈ F, we have E(1A (S)) = P (S −1 A) = P (A) = E(1A (ω)). So the equation is true for indicator function. Thus, the result is true when f is a simple function (finite linear combination of indicator functions). Assume f ≥ 0. Since there is a sequence of simple functions {fn }n∈N , s.t. fn ↑ f a.e. as n → ∞, we have Z Z f (Sω) = lim fn (Sω) n→∞ Z MCT lim fn (Sω). = n→∞ (MCT stands for the monotone converge theorem.) Thus, the equation is true for f ≥ 0. We can now finish the proof by writing f = f + − f − , where f + and f + are both non-negative functions. Lemma 2.2.9. If {An } is a special averaging sequence on Z2 , then for any f ∈ L1 (Ω, F, P ), the sequence {gn }n∈N defined as 1 X gn (ω) , f (T ` ω) |An | `∈An is uniformly integrable. 24 |f |dP = M < ∞, then calculate R Proof. First, set Z 1 X Z |gn |dP ≤ |f (T ` ω)|dP (ω) |An | `∈An 1 X Z = |f (ω)|dP (ω) |An | `∈An = M. Thus gn is uniformly L1 -bounded. Second, for any E ∈ F, calculate Z 1 X Z |gn | ≤ |f (T ` ω)| E |An | `∈An E 1 X Z = |f (ω)|. |An | `∈An T −` E We know that since f ∈ L1 , for any  > 0, there is a δ > 0, s.t. for any event E 0 ∈ F with P (E 0 ) < δ, we have Z |f | < . E0 Note that because of stationary, P (E) = P (T −` E) for any ` ∈ An . Thus, for any  > 0, there is a δ > 0, s.t. for any event E ∈ F with P (E) < δ, we have Z |gn | < , E and is independent of n. Now we have prepared to characterize the convergence, f ? , in Theorem 2.2.7. 25 Theorem 2.2.10. For f ∈ L1 (P ), we have f ? (ω) = E(f |I)(ω) P -a.e, where I is the σ-field of T -invariant events. Proof. We need to show that: (a) f ? ∈ I. (b) For any B ∈ I, f? = R R B B f. Note that (a) is satisfied because f ? is T -invariant function. As for the property (b), given any B ∈ I, calculate 1 X E(f ? 1B ) = lim E(n→∞ f (T ` )1B ) |An | `∈An 1 X u.i. lim E( f (T ` )1B ) = n→∞ |An | `∈An 1 X = n→∞lim E(f (T ` )1B ). |An | `∈An (u.i means uniformly integrable.) And calculate Z E(f (T ` )1B ) = f (T ` ω)dP (ω) B Z = f (ω)dP (ω) −` B ZT B∈I f (ω)dP (ω) = B = E(f 1B ). So, we have proved (b). 26 Thus we have proved Equation 2.2.5 on more general averaging sequence case. Note that E(f |I) doesn’t depend on the choice of the averaging set{An }, means that 1 f (T ` ) is the same for any averaging sequence {An }. P the limit of |An | `∈An 2.2.3.2 Ergodic Decomposition on Images The main purpose of this subsection is to provide a theoretical background on ergodic decomposition for Z2 lattice. Again, we assume our dynamic system (Ω, F, T, P ) on Z2 is stationary. To have a little understand of the ergodic decomposition theorem, the most simple form of the theorem is Z P (F ) = Pδ (F )dQψ (δ) for F ∈ F, where Pδ is a (stationary and ergodic) measure on (Ω, F), Qψ (·) = P (ψ −1 (·)) is a measure on the measures of (Ω, F), and ψ is the ergodic component function that assigns ω ∈ Ω to its ergodic component ψω , a stationary ergodic prob- ability on (Ω, F). As one can see, the stationary probability P is decomposed to ergodic components {Pδ }δ∈ψ(Ω) . Furthermore, the probability distribution ψω (·) can be characterized as P (·|ψ)(ω). Let’s give an simple example to describe this: 27 Example. Set Ω0 = {ω01 , ω10 , ω001 , ω010 , ω100 }, where ω01 = {0, 1, 0, 1, · · · }; ω10 = {1, 0, 1, 0, · · · }; ω001 = {0, 0, 1, 0, 0, 1, · · · }; ω010 = {0, 0, 1, 0, 0, 1, · · · }; ω100 = {1, 0, 0, 1, 0, 0, · · · }. 0 Define F 0 = 2Ω and µ(ω) , 21 P1 (ω) + 21 P2 (ω), where 1 P1 (ω) = I{ω01 ,ω10 } (ω), and 2 1 P2 (ω) = I{ω001 ,ω010 ,ω001 } (ω), 3 where I is an indicator function. Define the Transformation S as (Sω)i = ωi+1 . In this example we have a non-ergodic stationary dynamical system (Ω0 , F 0 , µ, S). One can verify it is not ergodic by the following fact: Set A , {ω10 , ω01 }, then A ∈ IS , the T 0 -invariant σ-algebra., and we have µ(A) = 1 2 ∈ / {0, 1}. On the other hand, (Ω0 , F 0 , µ, S) is an ergodic and stationary system for either i = 1, 2. In this case, µ can be decomposed to P1 , P2 , the ergodic components. One can also see the ergodic component function ψµ for µ can be defined as ψµ (ω) = P1 if ω ∈ {ω01 , ω10 } and ψµ (ω) = P2 if ω ∈ {ω001 , ω010 , ω001 }. The verification of the ergodic decomposition theorem for 2D is similar to that in 1D. However, it is hard to find literature for this. To prove the ergodic decomposi- tion theorem on images, we will follow the idea of Ch.8 of Gray’s book [16], where 28 everything is proved for 1D and asymptotic mean stationary (AMS) case, which is a weaker form of stationary property. Note that as we study through Gray’s book, there are some proofs that we find difficult to understand, so we will try to complete them as much as we can. Before introducing the ergodic decomposition, let’s first introduce the concept of standard measurable space. Some terminology will be used: First, a field under a sample space Ω0 is a collection of subsets of Ω0 that (1) has Ω0 in it, (2) is closed under complement, and (3) is closed under finite unions7 . Second, a field generated by a collection of subsets, say A, means it is the smallest field contains A. Last, a countable generating field means the field is generated by countable subsets. Definition 2.2.11. A measurable space (Ω0 , F 0 ) is standard if it has the countable extension property. That is, there is a countable generating field, G 0 s.t. σ(G 0 ) = F 0 and any set function µ on G can be uniquely extended to a probability measure on F 0 if µ has the following properties: (a) (normalization) µ(Ω0 ) = 1, (b) (nonnegativity) µ(G) ≥ 0 for all G ∈ G 0 , and Pn (c) (finite additivity) µ(∪ni=1 Gi ) = i=1 µ(Gi ) for Gi ∈ G 0 , Gi ∩ Gj = φ with i 6= j = 1, 2, · · · , n. There are many different definitions of the standard measure space, see Ch.2 of [16]. We define it like this because it is ready to use. In our setting, as we have seen 7 Note that σ-algebra is closed under countable unions 29 in the beginning of this section, the measurable space (Ω, F) has a nature countable ` generating set C , {χΛn }n∈N∪0,`∈Z2 . Thus if (Ω, F) is a standard space then we can define a set function on the field that generated by C, say G, then extend to F. The verification of (Ω, F) being a standard space involves the Carathéodory extension theorem and the Kolmogorov extension theorem, which can be found in Ch2 of [16]. We will simply put the theorem here. Theorem 2.2.12. (Ω, F) is a standard space. In the following, we will go through the formation of the ergodic components. Let’s define the set of the ergodic points 1 X E , ∩G∈G {ω : lim 1G (T ` ω) exists}. (2.2.6) n→∞ |A | n `∈An Note that G is countable, so we know that E is a measurable (because 1G (T ` ω) is measurable function for any ` ∈ Z2 ) and P (E) = 1. Consider a set function for G ∈ G   EP (1G |I)(ω) , ω∈E   µω (G) ,  , ω ∈ Ec  η(G)  for every ω ∈ Ω and some stationary ergodic distributions η on (Ω, F). Since µω is normalized, nonnegative, finitely additive set function on G, it can be extend to a probability measure on (Ω, F) (by Definition 2.2.11). So we can now write µω (F ) = Eµω (1F ) for any F ∈ F and ω ∈ Ω. Also, it is easy to verify that µω is a stationary ergodic probability distribution on (Ω, F). In fact,{µω }ω∈Ω is referred as the ergodic decomposition of the stationary distribution P . 30 It will be worth mentioning that it will be difficult to claim measurability if we define E as the intersection of sets on F, which is an uncountable intersection. Also, the reader might be tempted (like the author did!) to define µω (G) = EP (1G |I)(ω) 1 1G (T ` ω) converges for every G ∈ G, and µω (G) = η(G) P whenever limn→∞ |An | `∈An otherwise. However, it will be hard to claim the finite additivity of the set function, since there might be an ω (although in the measure 0 set) that the limit exists for set F1 but not set F2 . An important property of µω is that it is invariant. Lemma 2.2.13. µω , as a function of ω ∈ Ω, is T -invariant. i.e. µω (F ) = µT k ω (F ) for all ω ∈ Ω, F ∈ F, and k ∈ Z. Proof. Because of unique extendability of the standard space, we only need to consider events in G. For ω ∈ E and G ∈ G, we first need to show that T k ω ∈ E, for any 1 k ∈ Z2 . Define Zn , 1G (T ` ω) and Z , limn→∞ Zn (ω), then calculate for P |An | `∈An any k ∈ Z2 : |Zn T k − Z| ≤ |Zn T k − Zn | + |Zn − Z|. 31 The limit of the second term is 0. And the first term 1 X 1 X |Zn T k − Zn | = | 1G (T ` T k ω) − 1G (T ` ω)| |An | `∈An |An | `∈An 1 1G (T ` ω) − 1G (T ` ω)| X X = | |An | `∈T −k An `∈An 1 1G (T ` ω) X ≤ |An | `∈T −k An 4An |T −k An 4An | ≤ |An | → 0 as n → ∞ from the definition of special averaging set. Thus |Zn T k − Z| → 0. So we showed that if ω ∈ E, then T k ω ∈ E, and moreover, 1 X µT ω (G) = lim Zn T k = Z = lim 1G (T ` ω) = µω (G). n→∞ n→∞ |A | n `∈An In the last, for those ω ∈ / E for any T = T k , k ∈ Z2 . / E, it is easy to imply that T ω ∈ So µω = η = µT ω . With the definition of µω , we can further characterize Theorem 2.2.10 by the following Lemma. Lemma 2.2.14. For f ∈ L1 (P ), we have E(f |I)(ω) = Eµω (f ) for P -a.e. ω. 32 Proof. We only need to consider ω ∈ E. To prove the equation, it is easy to see that the result is true if f is a simple function. Next, suppose f ≥ 0, then we know there is a sequence of simple function {fk }, s.t. fk ↑ f P -a.e., Thus E(f |I)(ω) = E(lim fk |I)(ω) k M CT lim E(fk |I) = k = lim Eµω (fk ) k MCT Eµω (f ). = Thus the equality is true if f ≥ 0 and is integrable. For integrable f , we can then write f = f + − f − to complete the proof. Next, we introduce the ergodic component function. Let P = P((Ω, F)) be the space of probabilities on (Ω, F). And the ergodic component function is defined as ψ : Ω → P with ψ(ω) = ψω , µω , which assigns ω to its own stationary ergodic distribution. It will be convenient to have the following characterization: Eµω (f ) = EP (f |ψ(ω)) for f ∈ L1 (P ), which will require a σ-algebra on P. Also we would like to define a measure Qψ (·) = P (ψ −1 (·)) to measure the size of the partition of ω that have the same ergodic component, which means we need to have a measure space for P. The σ-algebra for the measure space can be induced by defining a metric on P, where we use the distributional distance here: ∞ 2−i |µ(Gi ) − ν(Gi )|, X dG (µ, ν) , i=1 33 where {Gi }∞ i=1 = G. (Remember that G is the field countably generated by C = ` {χΛn }n∈N∪0,`∈Z2 .) With this metric we can define the σ-algebra, say L, to be the Borel σ-algebra: L , σ({µ : dG (µ, ν) < r}v∈P,r≥0 ). Thus we have the following lemma, which is from Lemma 8.5.1 of [16]. Lemma 2.2.15. There is a measurable space (P, L) s.t. the map ψ : Ω → P defined as ψ(ω) , µω is measurable. Proof. We need to show that ψ −1 (A) ∈ F for any A ∈ L. Define S , {A ∈ L : ψ −1 (A) ∈ F} to be the collection of sets in L that make ψ measurable. It’s easy to see that S is a σ-algebra because ψ −1 preserve unions and complements (a general property for a mapping). We first show that for any ν ∈ P, and r ≥ 0, we have Br (v) ∈ S, where Br (v) , {µ ∈ P : dG (µ, ν) < r}. i.e. ψ −1 (Br (v)) ∈ F. To see this, one can write ψ −1 (Br (v)) = {ω ∈ Ω : ψ(ω) ∈ Br (v)} = {ω : dG (ψ(ω), ν) < r} ∞ 2−i |µω (Gi ) − ν(Gi )| < r}. X = {ω : i=1 Then we have ψ −1 (Br (v)) ∈ F because µω (Gi ) ∈ {E(1Gi |I)(ω), η(Gi )}, as a function of ω, is measurable. So Br (v) ∈ S, which means σ(Br (v)) ⊂ S and we have ψ −1 (A) ∈ F for any A ∈ L. 34 Since ψ actually maps Ω into stationary and ergodic distributions, say Pe , a subspace of P. We can further restrict the measure space (P, L) to (Pe , Le ), where Le is the Borel σ-algebra of Pe . However, it will require that Pe is a measurable set in (P, L), since we would like ψ to be measurable in Pe , too. The following lemma (slightly modified from Lemma 6.7.4 in [16], where it was proved in 1D and AMS case) is necessary for us to demonstrate this fact in Theorem 2.2.18. Lemma 2.2.16. A stationary dynamical system (Ω, F, Q, T ) is ergodic if and only if 1 X lim Q(T −` G ∩ G) = Q2 (G) (2.2.7) n→∞ |A | n `∈An for any G ∈ G, and a special averaging sequence {An }. Note that G is the (countable) generating field defined as the smallest field contains C. Proof. Note that for any f ∈ L1 , EQ (f |I) = EQ (f ) if Q is ergodic, since Q(A) ∈ {0, 1} for any A ∈ I from the definition of ergodicity. 1. (⇒) For any G ∈ G, calculate 35 1 X 1 X lim Q(T −` G ∩ G) = lim EQ (1T −` G∩G ) n→∞ |An | `∈An n→∞ |An | `∈An 1 X Z = lim 1T −` G (ω)dQ(ω) n→∞ |A | n `∈An G Z 1 X u.i. lim 1T −` G (ω)dQ(ω) = G n→∞ |An | `∈A n Z ergodic E(1G )dQ(ω) = G = Q2 (G). Note that the uniformly integrability (u.i.) is from Lemma 2.2.1, which requires Q to be a stationary measure. 2. (⇐) First, understand that for any A ∈ I, limn→∞ 1 Q(T −` A ∩ A) = P |An | `∈An Q(A). So we have Q(A) = Q2 (A), and this implies Q(A) ∈ {0, 1}, which proves the system is ergodic. The hard part is that the Equation 2.2.7 is valid only for events on the generating field G, and I is not necessary contained in G. So we have to argue that Equation 2.2.7 holds for any event F ∈ F. We know that for any  ∈ (0, 1) and F ∈ F there is an F0 ∈ G s.t. Q(F 4F0 ) <  (see Corollary 1.5.3 in [16]). Set 1 X (?) = | Q(T −` F ∩ G) − Q2 (F )|, |An | `∈An 36 and from the triangular inequality we can calculate 1 X 1 X (?) ≤ | Q(T −` F ∩ F ) − Q(T −` F0 ∩ F0 )| |An | `∈An |An | `∈An 1 X +| Q(T −` F0 ∩ F0 ) − Q2 (F0 )| |An | `∈An +|Q2 (F0 ) − Q2 (F )|. The second term will converge to 0 when n → ∞ by the condition we have. The third term will be less than 2 , because |Q(C) − Q(D)| ≤ Q(C4D) for any event C and D. And the first term will be less than 1 X |Q(T −` F ∩ F ) − Q(T −` F0 ∩ F0 )|. |An | `∈An Note that we have 1 |Q(T −` F ∩ F ) − Q(T −` F0 ∩ F0 )| ≤ Q(T −` F ∩ F 4T −` F0 ∩ F0 ) 2 = () We can do a further estimate of () by Q(C4D) ≤ Q(C4E) + Q(D4E), so () ≤ Q(T −` F ∩ F 4T −` F0 ∩ F ) + Q(T −` F0 ∩ F0 4T −` F0 ∩ F ) ≤ Q(T −` F 4T −` F0 ) + Q(F0 4F ) = Q(T −` (F 4F0 )) + Q(F0 4F ) = 2Q(F0 4F ). 37 Thus the first term is bounded by 2, this implies 1 X lim | Q(T −` F ∩ G) − Q2 (F )| ≤ 2. n→∞ |An | `∈An We can then let  → 0 and concludes that Equation 2.2.7 holds for any event F ∈ F. We are now prepared to claim that Pe is a measurable set. Define Ps to be the collection of stationary probability distributions in P. The next lemma is modified from Lemma 8.5.2 in [16]. Lemma 2.2.17. We have (a) Ps , the collection of stationary distributions in P, is a close set in (P, L), and (b) Pe , the collection of stationary ergodic distributions in P, is a measurable set in (P, L). Proof. (a) Given any sequence {mn }n∈N ⊂ Ps , s.t. mn (F ) → m(F ) as n → ∞ for any F ∈ F. Given any G ∈ G, from the definition of G (G is a field generated by C) we know that T ` G ∈ G ⊂ F, for any ` ∈ Z2 . So m(G) = m(T ` G) for G ∈ G since m(G) = lim mn (T ` G) = lim mn (G) = m(T ` G). n→∞ n→∞ 38 By the uniqueness of the measure extension property (i.e. Definition 2.2.11 on standard space (⊗, F)), we conclude that m(F ) = m(T ` F ) for any F ∈ F, i.e. m ∈ Ps . Hence it is closed. (b) From Lemma 2.2.16, we know that 1 X Pe = Ps ∩ {m ∈ P : n→∞ lim m(T −` G ∩ G) = m2 (G) for any G ∈ G} |An | `∈An (note that Pe ⊂ Ps by definition). From (a) we know that Ps is closed, hence measurable. So we only need to show the second set is measurable. Furthermore, since the second set equals to 1 X ∩G∈G {m : n→∞ lim m(T −` G ∩ G) = m2 (G)} |An | `∈An and the fact that G is countable and T −` G ∈ G, we only need to show that for any G ∈ G, m(G), as a function of m, is a measurable function on (P, L). We can show this by proving CG , {m0 ∈ P : m0 (G) < r} is an open set on the metric space the (P, dG ) for any r > 0. The proof follows: Given any m0 ∈ CG , there is an  > 0 s.t. m0 (G) < r − . Given δ < 2−k , for any v ∈ {v 0 ∈ P : dG (m0 , v 0 ) < δ}, we have 2−i |m0 (Gi ) − v(Gi )| < δ, X Gi ∈G which implies |m0 (G) − v(G)| < 2k δ for some k ∈ N. If v(G) < m0 (G) then v ∈ CG . If v(G) > m0 (G), we have v(G) < m0 (G) + 2k δ < r −  + 2k δ < r. 39 Thus we have proved that for any m0 ∈ CG , there is a δ s.t. for any v ∈ {v 0 ∈ P : dG (m0 , v 0 ) < δ}, v ∈ CG , i.e. CG is open. This implies CG ∈ L. Hence m(G) is a measurable function on (P, L) for any G ∈ G. We have prepared ourselves all the background that needed for the ergodic de- composition theorem on images, i.e. (Ω, F), with respect to ergodic component ψ. The following is modified from Theorem 8.5.1 in [16]. Theorem 2.2.18. (Ergodic Decomposition) Under (Pe , Le ), the measurable space defined on stationary and ergodic distributions on (Ω, F), a function ψ : Ω → Pe is defined as   E(1F |I)(ω) ,ω∈E   ψ(ω)(F ) , ψω (F ) ,  , ω ∈ Ec  η(F )  for any F ∈ F and for some η ∈ Pe , where E is defined in Equation 2.2.6 with P (E) = 1. Then we have the following properties: (a) ψ is measurable and T -invariant, i.e. for any T = T k , k ∈ Z2 , ψ(T ω) = ψ(ω), or say ψT ω (F ) = ψω (F ) for any F ∈ F. (b) On (Pe , Le ), there is a well-defined probability distribution Qψ (C) = P (ψ −1 (C)) for C ∈ Le . And for any f ∈ L1 (P ), we have Z Z EP f = Eψω f dP (ω) = EPδ (f )dQψ (δ), ω∈Ω ψω (F )dP (ω) = Pδ (F )dQψ (δ) for F ∈ F. R R which implies P (F ) = ω∈Ω 40 (c) For any f ∈ L1 (P ), we have Eψ (f ) = EP (f |ψ), which implies ψ(F ) = P (F |ψ) for F ∈ F. Note that we wrote Pδ = δ for δ ∈ Pe in this theorem. Proof. We will prove these statements on generating fields G, since we can extend the result to F because (Ω, F) is standard space. (a) (Measurabilility) From Lemma 2.2.15, we know that for any A ∈ L we have ψ −1 (A) ∈ F. Since Pe is a measurable subset of P from Lemma 2.2.17, any open set in Pe is measurable in (P, L). This implies Le ⊂ L. (Note that both of them are Borel σ-algebra.) So given any A ∈ Le , we have A ∈ L. Thus we can conclude ψ is measurable. (Invariance) The invariance simply follows from Lemma 2.2.13. (b) (Well-definedness) First, by (a) ψ is measurable, so P (ψ −1 (C)) has meaning for any C ∈ Le . It remains to show that it is a probability measure. Obviously it is non-negative. Next let’s check if it have normalization property: Qψ (Pe ) = P (ψ −1 (Pe )) = 1 since for every ω, ψ(ω) ∈ Pe . The countable additivity can be easily proved because ψ −1 preserves unions, and complements. Thus we conclude Qψ is a well-defined probability measure on (Pe , Le ). (First equality) For any f ∈ L1 (P ) calculate EP (f ) = EP (EP (f |I)) Lem 2.2.14 EP (Eψω (f )) = Z = Eψω (f )dP (ω). ω∈Ω 41 (Second equality) Calculate Z Z Z Eψω (f )dP (ω) = f (v)dψω (v)dP (ω) ω∈Ω ω∈Ω v∈Ω Z Z δ = ψ(ω) f (v)dPδ (v)dP (ψ −1 (δ)) = δ∈ψ(Ω) v∈Ω Z Z (?) f (v)dPδ (v)dP (ψ −1 (δ)) = δ∈Pe v∈Ω Z = EPδ (f )dQψ (δ), where in (?), we used Qψ (Pe ) = Qψ (ψ(Ω)) = 1, so Qψ (Pe \ψ(Ω)) = 0. (c) We are going to show: (1) Eψ (f ) ∈ σ(ψ), and (2) f for B ∈ σ(ψ). R R B Eψ (f ) = B (1) For any f ∈ L1 , it suffices to show that a function qf : Pe → R defined as q(µ) , Eµ (f ) is a continuous function. Since for any a ∈ R if qf is continuous then qf−1 ((−∞, a)) is an open set in (Pe , Le ), thus {ω : Eψω f < a} = ψ −1 ({µ ∈ Pe : Eµ (f ) < a}) = ψ −1 qf−1 ((−∞, a)) ∈ ψ −1 (Le ) , σ(ψ) Let’s first show that when f = 1G , G ∈ G, the generating field, q1G is continuous. So, we have q(µ) = µ(G) for µ ∈ Pe . Set δ > 0, for any µ, ν ∈ Pe with dG (µ, v) < δ. Thus, by definition of dG , we have for G = {Gi }ni=1 ∞ 2−i |µ(Gi ) − ν(Gi )| < δ, X i=1 42 which implies |µ(G) − ν(G)| < 2k δ for some k ∈ N. So, for any µ ∈ Pe and for any  > 0, we can choose δ < 2−k , s.t. for any ν ∈ Pe with dG (µ, v) < δ, we have |q1G (µ) − q1G (ν)| = |µ(G) − ν(G)| < , i.e. q1G is continuous for G ∈ G. Next, to see that q1F is continuous for F ∈ F, define S , {F ∈ F : q1F is continuous.}. It is easy to see that S is a σ-algebra. (Since, for any µ ∈ Pe , first, given any A ∈ S, q1Ac (µ) = 1−µ(A) = 1−q1A (µ) ∈ S ; second, given any {Ai }∞ i=1 ⊂ S ∈ S.) P P that is pairwise disjointed, we have q1∪Ai (µ) = i µ(Ai ) = i q1Ai And, by previous statement, G ∈ S for G ∈ G, so F = σ(G) ⊂ S. So we have proved q1F is continuous for F ∈ F. Thus it is easy to see that q1f is continuous if f is any simple function. For f ≥ 0, we can then define {fk } to be the sequence of simple functions s.t. fk ↑ f P -a.e., and using monotone convergence theorem to conclude q1f is continuous for f ≥ 0. In the last, for general f ∈ L1 we have q1f is continuous since we can write f = f + − f − . (2) Since ψ is invariant, we have ψ −1 (A) ∈ I, for any A ∈ Le . This implies 43 σ(ψ) ⊂ I. So for any B ∈ σ(ψ), we have B ∈ I, and Z f (ω)dP (ω) = EP (EP (f 1B |I)) B B∈I EP (1B EP (f |I)) = Z Lem 2.2.14 Eµω (f )dP (ω) = ZB = Eψω (f )dP (ω). B So we have the result that Eψ (f ) = EP (f |ψ). The theorem gives many useful conclusions. First, a stationary probability distri- bution can be viewed as weighted sum of stationary ergodic distributions, where the weights is indexed by elements in Pe and not by that in Ω. Second, We can divide Ω into several ergodic components, the ergodic components about ω can be defined as the set ψ −1 (ψ(ω)) which tells you about which ergodic stationary distribution that ω belong to, with size Qψ (ψ(ω)). Third, ψ is actually the probability of P condition on the ergodic component ψ. Also note that ψ and Pe doesn’t depend on averaging sequences, so we can say that ψ and Pe are induced by the stationary dynamical system (Ω, F, P, T ). Let’s write down the useful corollary we got so far, for a future use. Corollary 2.2.19. Assume (Ω, F, P, T ) we defined on 2D-lattice is stationary with induced ergodic component function ψ. For any f ∈ L1 (P ) and ant special averaging 44 sequence {An } 1 X lim f (T ` ω) = E(f |I)(ω) = Eψ (f ) n→∞ |An | `∈An for P -a.e. ω ∈ Ω. 2.2.3.3 Ergodic Decomposition of the Entropy Rate In this subsection, we will introduce one of the important results regarding the ergodic decomposition. Since, from Theorem 2.2.18, we saw that a stationary P can be viewed as the combination of stationary ergodic distributions {Pδ }δ∈Pe , a nature question is that: Can we view the entropy rate hP as a combination of {hPδ }δ∈Pe ? The answer is yes. The following fundamental theorem helps us understand the relationship between entropy rate and ergodic decomposition, which is due to Jacobs [21], and is well-described in Section 2.4 in [17]. We will not give a complete proof here, since most of them is the same as in 1D setting. Instead, we will briefly describe the proof and give a complementary argument for properties that are not clearly extendable to 2D-lattice. Theorem 2.2.20. (The Ergodic Decomposition of the Entropy Rate) We have Z hP = hψω dP (ω) ω∈Ω To prove this, the reader might be tempted to use property (b) on Theorem 2.2.18 and Lebesgue dominate convergence theorem to exchange integral and limit, however, 1 since H (XΛn (ω)) |Λn | ψω converge almost surely in ψ (not in P ), we cannot apply the 45 theorem here. It’s actually a non-trivial result. The basic idea is based on Theorem 8.9.1 of [16], which says if a function D : Ps → R ≥ 0, where Ps is the set of stationary distributions on (Ω, F), satisfies (a) D(ψω ) ∈ L1 (P ) for any ω ∈ Ω, (b) D(λµ1 + (1 − λ)µ2 ) = λDµ1 + (1 − λ)D(µ2 ) for λ ∈ (0, 1), µ1 , µ2 ∈ Ps , and (c) D is upper semi-continuous, then Z D(P ) = D(ψω )dP (ω). 1 So we can prove Theorem 2.2.20 by choosing D(P ) = hP , limn→∞ H (XΛn ). |Λn | P Note that D satisfies (b) (this is from Lemma 2.3.4 of [17]): First, note that we have X I = Hλµ1 +(1−λ)µ2 (X) = −λ µ1 (x) log(λµ1 (x) + (1 − λ)µ2 (x)) x X −(1 − λ) µ2 (x) log(λµ1 (x) + (1 − λ)µ2 (x)) x = λEµ1 (− log(λµ1 (X) + (1 − λ)µ2 (X))) +(1 − λ)Eµ2 (− log(λµ1 (X) + (1 − λ)µ2 (X))). 46 So the the expectation in the first term of I can be estimated: Eµ1 (− log(λµ1 (X) + (1 − λ)µ2 (X))) µ2 (X) =Eµ1 (− log(µ1 (X)) − log(λ + (1 − λ) )) µ1 (X) µ2 (X) =Hµ1 (X) + Eµ1 (− log(λ + (1 − λ) )) µ1 (X) (?) µ2 (X) ≥Hµ1 (X) + Eµ1 (1 − λ − (1 − λ) )) µ1 (X) =Hµ1 (X) + 1 − λ − (1 − λ) =Hµ1 (X) Note that at (?) we use the inequality− log x ≥ 1 − x. With the same reason we have the expectation in the second term of I large or equal to Hµ2 (X). Thus we have the lower bound of I as I ≥ λHµ1 (X) + (1 − λ)Hµ2 (X). Second, the upper bound of I can be estimated with I ≤ λEµ1 (− log(λµ1 (x)) (1 − λ)Eµ2 (− log((1 − λ)µ2 (x)) = λHµ1 (X) + (1 − λ)Hµ2 (X) + g(λ), where g(λ) = −λ log(λ) + (1 − λ) log(1 − λ). Thus, by dividing |Λn | and let n → ∞ on both bounds of I, we have the property (b). Since most of the proof is the same as in Section 2.4 of [17], we will not repeat 47 here. In particular, when proving (c), the upper semi-continuity, it requires hP = 1 inf n∈N H (XΛn ), |Λn | P which can be proved by 1 1 X HP (XΛn ) = HP (X` |X(<`)∩Λn ) |Λn | |Λn | `∈Λn 1 X ≥ HP (X` |X(<`) ) |Λn | `∈Λn 1 X stationary Hp (X0 |X(<0) ) = |Λn | `∈Λn = H(X0 |X(<0) ) Cor.=2.2.4 hP for all n ∈ N. 2.2.4 Markov Approximation on Images Here we introduce the Markov approximation, for stochastic process on Z2 . Markov approximation will play a crucial role as proving Theorem 2.2.34. Note that in 1D, the kth-order Markov approximation is defined as n−1 pk (Y0:n−1 ) = p(Y0:k−1 ) Y p(Yi |Yi−k:i−1 ) i=k (see ch16.8 of [9], where the indices are defined on {0, 1, · · · , ∞} and the notation a : b means {a, a + 1, · · · , b}). This is an vital part when proving the original Shannon- McMillian-Breiman (SMB) theorem. And we will extend the concept to 2D-lattice. Specifically, the properties we will develop in this section will help Lemma 2.2.32, the upper bound of the non-ergodic SMB theorem in the next section. The reader 48 Figure 2.2.4: Example of rectangle sets. can first proceed to the next subsection, where we will prove the non-ergodic SMB theorem on 2D, and come back to this section when one needs supports on Lemma 2.2.32. First, we define the set that can be approximated. Definition 2.2.21. Set A ⊂ Z2 is a rectangle set if for some `, i ∈ Z2 and n ∈ N∪{0}, A = (< i) ∩ Λ`n . Also, we define A to be the collection of all possible rectangle sets. A quick examples are Λ`n and Bm ` ` for ` ∈ Z2 and n, m ∈ N ∪ {0}, where Bm , 0 ` + Bm ⊂ Z2 and Λ`n = ` + Λ0n . Remember that Bm 0 = Bm , Λ0m ∩ (< 0) and Λ0n = Λn . (as depicted as in Figure 2.2.4). Definition 2.2.22. (Markov approximation) A set function P m : X A → R is the mth-order Markov approximation P m to the probability distribution P s.t. P m (xA ) = Y P (X` = x` |XBm ` ∩A = xB ` ∩A ) m `∈A 49 for all A ∈ A. Moreover, define  P m (xA∪B ) , if P m (xB ) > 0     P m (x m B) P (xA |xB ) ,   0  , o.w. if A ∪ B ∈ A and B ∈ A. We will use the abbreviation P m (XA ) as a random variable P m (XA (ω)). Same abbreviation rule applies to P . In the next lemma, we introduce some property of the Markov approximation. Lemma 2.2.23. (P m properties) P m has the following properties: P m (xA ) = 1 for all A ∈ A. P 1. xA ∈XA 2. P m (xA ) = P m (xA+` ) for A ∈ A, ` ∈ Z2 . 3. The conditional probability on mth-order area is the same as the original one, i.e. P m (x0 |xBm 0 ) = P (x0 |xBm 0 ). 4. P m (x0 |xA ) = P m (x0 |xBm ) if for some k ∈ Z2 , A = (< 0) ∩ Λkn ⊃ Bm 0 . 5. P  P m on X A . That is, P m (xA ) = 0 implies P (xA ) = 0 for all A ∈ A. 6. Pδ  P m on X A for all Pδ ∈ Pe , where Pe is the collection of all the stationary ergodic distribution induced by the dynamical system (Ω, F, P, T ). 50 Proof. All the proofs are straightforward: 1. Given any rectangle set A, for some ` ∈ Z2 , and denote XA as the set of all the configuration on A. Set {`1 > `2 > · · · > `|A| } = A Then |A| m X X Y P (xA ) = P (x`i |xBm `i ∩A ) xA ∈XA xA ∈XA i=1 |A| X X Y = P (x |x `1 `1 ∩A Bm ) P (x`i |xBm `i ∩A ) x`1 ∈X xA\{`1 } ∈XA\{`1 } i=2 |A| X Y = 1· P (x`i |xBm `i ∩A ). xA\{`1 } ∈XA\{`1 } i=2 P m (xA ) = 1. P Then iteratively we can get xA ∈XA 2. Trivial result since P is stationary. 3. This is because P m (xBm0 ∪{0} ) P m (x0 |xBm 0 ) = m P (xBm 0 ) P (x0 |xBm 0 P (x` |xB ` ∩B 0 ) Q 0 ) `∈Bm m m = |x Q `∈Bm0 P (x ` Bm` ∩B 0 m ) = P (x0 |xBm 0 ). 4. Since P m (xA∪{0} ) P m (x0 |xA ) = P m (xA ) P (x0 |xBm `∈A P (x` |xBm Q 0 ) ` ∩A ) = `∈A P (x` |xBm Q ` ∩A ) = P (x0 |xBm 0 ) = P m (x0 |xBm 0 ). 51 5. Given any rectangle set A ⊂ Z2 and xA ∈ B, s.t. P m (xA ) = 0, this will imply ` ∩A ) = 0 for some ` ∈ A. Thus P (x` , xB ` ∩A ) = 0 for some `. And that P (x` |xBm m P (xA ) = P (xA |x` , xBm ` ∩A )P (x` , xB ` ∩A ). m 6. Because of Pδ  P. The reader must be aware that although P m looks like a probability distribution in a lot of ways, it is not a probability distribution, since the definition of rectangles would not admit a σ-algebra. We will just treat P m as a special set function on rectangle sets, which is enough for the use here. The following lemma give the coding rate for Markov approximation, which is modified to 2D from Lemma 3.1.1 of [17]. Lemma 2.2.24. For any m ∈ N, there is a T -invariant function hm : Ω → R, s.t. −1 lim log P m (XΛn ) = hm P -a.e. n→∞ |Λn | Furthermore, hm can be specified as  −1 log P m (XΛn )(ω) converges.  Eψ(ω) (− log P (X0 |XBm )) , if limn→∞   m |Λn | h (ω) ,  log |X |   , o.w. Proof. Calculate −1 −1 X log P m (XΛn ) = log P (X` |XBm` ∩Λ ) n |Λn | |Λn | `∈Λn −1 X −1 X = log P (X` |XBm` ∩Λ ) + n log P (X` |XBm ` ∩Λ ), n |Λn | `∈Znm |Λn | `∈Λn \Z m n 52 where Znm , {` ∈ Λn |Bm ` ⊂ Λn }, as depicted in Figure 2.2.3. Since E(− log P (X` |XBm ` ∩Λ ) n is finite (because X is of finite alphabet), we know that − log P (X` |XBm ` ∩Λ ) is also n finite a.e. Thus the second term goes to 0 when n get sufficiently larger. With the same reason, the limit won’t change if we replace the second term with −1 X ` log P (X0 |XBm 0 )T (ω) |Λn | `∈Λn \Z m n which is also goes to 0 when n goes to ∞. And the first term: −1 X −1 X ` log P (X` |XBm ` ∩Λ )(ω) n = log P (X0 |XBm 0 )T (ω) |Λn | `∈Znm |Λn | `∈Znm because of stationary condition. Thus −1 −1 X lim log P m (XΛn ) = lim log P (X0 |XBm ` 0 )T (ω) n→∞ |Λ | n→∞ |Λ | n n `∈Λn a.e.Eψω (− log P (X0 |XBm 0 )) = which is from the result of ergodic theorem, i.e. Corollary 2.2.19. Let’s give the definitions of relative entropy and relative entropy rate about P m in the following definitions, which will be used in the next lemma. Definition 2.2.25. Given a stationary distribution P 0 on (Ω, F) such that P 0  P , the relative entropy of P 0 and P m on Λn is defined as P m (XΛn ) HP 0 ||P m (XΛn ) , −EP 0 (log ) P 0 (XΛn ) 53 for n ∈ N. Also, the relative entropy rate of P 0 and P m on Λn is defined as 1 lim hP 0 ||P m , n→∞ HP 0 ||P m (XΛn ) |Λn | P m (XΛn ) Note that the definition of relative entropy, −EP 0 (log P 0 (XΛn ) ), is well-defined because of the absolute continuous assumption (Lemma 2.2.23 (5) and (6)). Most of the time P 0 will be used as P or ψ. The following lemmas (Lemma 2.2.26, 2.2.28 and 2.2.29) will be used in proving SMB theorem for stationary non-ergodic source (Theorem 2.2.34). They next lemma is modified to 2D version from Lemma 2.4.3 [17]. Lemma 2.2.26. Given a stationary distribution P 0 on (Ω, F) such that P 0  P we have hP 0 ||P m = −hP 0 + EP 0 (− log P (X0 |XBm 0 )), 1 where hP 0 , limn→∞ H 0 (XΛn ). |Λn | P Proof. Since 1 1 P 0 (xΛn ) log P m (xΛn )}. X HP 0 ||P m (XΛn ) = {−HP 0 (XΛn ) − |Λn | |Λn | x ∈X Λn Λn 54 The first term will be −hP 0 as n → ∞. Calculate the second term: P 0 (xΛn ) log P m (xΛn ) = P 0 (xΛn ) X X X log P (x` |xBm ` ∩Λ ) n xΛn xΛn `∈Λn P 0 (xΛn ) log P (x` |xBm X X = ` ∩Λ ) n `∈Λn xΛn P 0 (xΛn ) log P (x` |xBm X X = ` ∩Λ ) n m xΛ `∈Zn n P 0 (xΛn ) log P (x` |xBm X X + ` ∩Λ ) n m xΛ `∈Λn \Zn n = (1) + (2) ` where Znm , {` ∈ Λn |Bm ⊂ Λn }. Calculate term (1): P 0 (xBm X X (1) = ` ∪{`} ) log P (x` |xB ` ∩Λ ) m n m x ` `∈Zn B m ∪{`} P 0 (xBm X X = ` ∪{`} ) log P (x` |xB ` ) m m `∈Zn xB ` m ∪{`} P 0 (xBm X X stationary 0 ∪{0} ) log P (x0 |xB 0 ) m = m `∈Zn xB 0 m ∪{0} = |Znm |EP 0 (log P (x0 |xBm 0 )). We know that |Znm |/|Λn | → 1 as n → ∞, so 1 X X 0 lim P (xΛn ) log P (x` |xBm ` ∩Λ ) n = EP 0 (log P (x0 |xBm 0 )) n→∞ |Λ | n `∈Z m xΛn n As for term (2), we know that there are only finitely many values (with order O(m2 )) P 0 (xΛn ) log P (x` |xBm m ` ∩Λ ) for different ` ∈ Λn \Zn , and each value is finite. P of xΛn n Along with the fact that |Λn \Znm |/|Λn | → 0 as n → ∞, we know that 1 P 0 (xΛn ) log P (x` |xBm X X lim ` ∩Λ ) n = 0. n→∞ |Λn | `∈Λn \Z m xΛn n 55 Thus we have 1 lim HP 0 ||P m (XΛn ) = −hP 0 + EP 0 (− log P (X0 |XBm 0 )). n→∞ |Λn | Lemma 2.2.27. We have hP ||P m ≥ 0 for all m ∈ N. Proof. Since for any n ∈ N, P m (XΛn ) HP ||P m , −EP (log ) P (XΛn ) 1 1 = EP (log m − log ) P (XΛn ) P (XΛn ) X 1 X 1 = EP ( log − log ) `∈Λn P (X` |XBm` ∩Λ ) n `∈Λn P (X` |X(<`)∩Λn ) X 1 1 = EP (log − log ). `∈Λn P (X` |XBm` ∩Λ ) n P (X` |X(<`)∩Λn ) And we have 1 1 EP (log ) ≥ EP (log ), P (X` |XBm ` ∩Λ ) n P (X` |X(<`)∩Λn ) ` which is because Bm ∩ Λn ⊂ (< `) ∩ Λn . So HP ||P m (XΛn ) > 0 for any n ∈ N, which 1 implies hP ||P m , limn→∞ H m (XΛn ) |Λn | P ||P > 0. 56 The following result is intuitive, says that when m gets larger, then P m will be closer to P . The original 1D version of this can be found in Lemma 2.6.1 of [17]. Lemma 2.2.28. limm→∞ hP ||P m ↓ 0. Proof. From Lemma 2.2.26, since hP ||P m (XΛn ) = −hP + EP (− log P (X0 |XBm 0 )) = −hP + HP (X0 |XBm 0 ) ↓ 0 when m → ∞, which is from Lemma 2.2.3. The following lemma gives the relation between ψ and P m . It is extracted and modified to 2D from the prove from Lemma 3.3.2 of [17]. This will be used in proving Lemma 2.2.32 on the next subsection. Lemma 2.2.29. inf m∈N hψω ||P m = 0 P − a.e. Proof. The proof follows the idea of Lemma 2.2.28. First, calculate Z Z inf hψω ||P m dP (ω) ≤ inf hψω ||P m dP (ω) = (?). m∈N m∈N 57 Then Z (?) Lem =2.2.26 inf dP (ω){−hψω + Eψω (− log P m (XΛn (ω))} m∈N Z Thm 2.2.20 −hP + inf dP (ω)Eψω (− log P m (XΛn (ω))) = m∈N Z = −hP + inf dP (ω)Eψω (− log P (X0 |XBm 0 )) m∈N Thm 2.2.18 (b) −hP + inf EP (− log P (X0 |XBm 0 )) = m∈N Lem 2.2.26 inf hP ||P m = m∈N Lem 2.2.28 0 = Along with the fact that the integrand, inf m∈N hψω ||P m , is non-negative, so inf m∈N hψω ||P m must be 0 P − a.e. 2.2.5 Coding Rate on a Stationary Source In this subsection we extend the Shannon-McMillan-Breimen (SMB) theorem to a 2D source which is stationary but not necessary ergodic. The dynamical system (Ω, F, P, T ), as we defined it for 2D-lattice in the beginning of this section, is only assumed to be stationary. In this case, the coding rate with P will no longer be a constant but a random variable. In fact, we will see, in Theorem 2.2.34, the coding rate converges to the entropy rate of the source’s ergodic component, i.e. log P (XΛn )(ω) lim − = hψω , n→∞ |Λn | where the ergodic component function ψ is defined in Theorem 2.2.18. Note that 58 ψ is induced from (Ω, F, P, T ). We will prove this by the idea from Gray’s book [17], where it was proved in 1D, in the end of this subsection (Theorem 2.2.34) after we prepared several lemmas. The following Lemma is Lemma 3.3.1 of [17] for 1D version. Lemma 2.2.30. (Coding rate for the ergodic components) Assume (Ω, F, P, T ) is stationary. We have 1 P ({v ∈ Ω : n→∞ lim − log ψω (XΛn (ω)) = hψω }) = 1, |Λn | 1 1 where hψω , limn→∞ H (XΛn ) |Λn | ψω = limn→∞ E (− log ψω (XΛn )) |Λn | ψω Remark. That this is not the result of Corollary 2.2.2, where, in terms of the notation here, it says 1 log ψω (XΛn (v)) → hψω |Λn | for ψω − a.e. v ∈ Ω, instead of P -a.e. v ∈ Ω. Proof. Define 1 G = {ω : lim − log ψω (XΛn (ω)) = hψω } |Λn | 1 Gδ = {ω : lim − log Pδ (XΛn (ω)) = hPδ } |Λn | for all Pδ ∈ Pe . Because of the Theorem 2.2.18(b), we know that Z P (G) = Pδ (G)dQψ (δ), 59 and Pδ (G) = P (G|ψ = δ) = P (G ∩ {ψ = δ}|ψ = δ) = P (Gδ |ψ = δ) = Pδ (Gδ ) = 1. The last equality is because of the Corollary 2.2.2 and the fact that Pδ is a stationary ergodic distribution for XΛn (ω), for all ω with ψ(ω) = δ. Thus P (G) = 1. We can now prove the first part of the Theorem 2.2.34. It is similar to the first part of the proof of Lemma 3.3.2 in Gray’s book [17] for 1D version. Note that Gray’s argument in the proof has some defects inside, and we manage to correct that here. Lemma 2.2.31. (Non-ergodic SMB lower bound) Assume (Ω, F, P, T ) is stationary, then for P -a.e. ω 1 ψω (XΛn (ω)) lim inf log ≥ 0. |Λn | P (XΛn (ω)) Proof. Note that for any event F ∈ F, if P (F ) = 0, we have ψ(F ) = 0 P -a.e. from the ergodic decomposition. In the following, we will omit ω for simplicity. Given any  > 0, we have 1 P (XΛn ) P (XΛn ) P( log > ) = P( > e|Λn | ) |Λn | ψ(XΛn ) ψ(XΛn ) Z P (XΛn ) Thm 2.2.18(b) P( > e|Λn | |ψ = δ)dQψ (δ) = ψ(XΛn ) Z P (XΛn ) = Pδ ( > e|Λn | )dQψ (δ) Pδ (XΛn ) Z P (XΛn ) −|Λn | ≤ EPδ ( )e dQψ (δ), Pδ (XΛn ) by Markov inequality. Note that we use the notation δ = Pδ in equations. Next, 60 define Aδn , {an ∈ X Λn |Pδ (an ) > 0}, we can evaluate the expectation by the following observation: P (XΛn ) X P (an ) EPδ ( ) = Pδ (an ) Pδ (XΛn ) an ∈Aδn Pδ (an ) ≤ 1, which concludes 1 P (XΛn ) P( log > ) ≤ e−|Λn | . |Λn | ψ(XΛn ) So we have ∞ X P (XΛn ) P ( |Λ1n | log ψ(XΛn ) > ) < ∞. n=1 P (XΛn ) Then by Borel-Cantelli Lemma, P ( |Λ1n | log ψ(XΛn ) > , i.o.) = 0, which means 1 P (XΛn ) lim sup log ≤ . |Λn | ψ(XΛn ) P (XΛn ) Since  is arbitrary, we have lim sup |Λ1n | log ψ(XΛn ) ≤ 0, i.e. 1 ψ(XΛn ) lim inf log ≥ 0. |Λn | P (XΛn ) Next Lemma gives the other direction, which is the second part of the proof of Lemma 3.3.2 in [17] for 1D version. Note that it is tempted to use the same method as 1 ψ(X ) in the Lemma 2.2.31, and using the Markov inequality on |Λn | log P (XΛΛn ) ,which involve n 61 ψ(X ) the calculation of EP ( P (XΛΛn ) ). This expectation is hard to bound since it is difficult n R to evaluate the integral ψω (XΛn (ω) )dω. Instead, we use the Markov approximation to help estimate. Lemma 2.2.32. (non-ergodic SMB upper bound) Assume (Ω, F, P, T ) is stationary, then for P -a.e. ω 1 ψω (XΛn (ω)) lim sup log ≤ 0. |Λn | P (XΛn (ω)) Proof. We are going to use the Markov approximation, P m , to estimate the upper bound. First, let’s establish the following claim. P m (XΛn (ω)) Claim 2.2.33. For any m ∈ N, lim sup |Λ1n | log P (XΛn (ω)) ≤ 0 P -a.e. ω Proof. We will using the method similar to the way we used to prove the lower bound. For any  > 0, the Markov inequality shows 1 P m (XΛn ) P m (XΛn ) P( log ≥ ) = P ( ≥ e|Λn | ) |Λn | P (XΛn ) P (XΛn ) P m (XΛn ) −|Λn | ≤ EP ( )e P (XΛn ) = 1 · e−|Λn | . So by Borel-Cantelli Lemma we have 1 P m (XΛn ) P( log ≥ , i.o) = 0. |Λn | P (XΛn ) Thus, we have 1 P m (XΛn ) lim sup log ≤ 0 P − a.e., |Λn | P (XΛn ) 62 since  > 0 is arbitrary. Because the this claim, we can calculate 1 1 1 1 lim sup log ≤ lim sup log m |Λn | P (XΛn (ω)) |Λn | P (XΛn (ω)) = Eψω (− log P (X0 |XBm 0 )), which is the result of Lemma 2.2.24 (SMB for P m ). Subtract lim |Λ1n | log ψω (XΛ1 , n (ω)) i.e. hψω (Lemma 2.2.30), on both side, we have 1 ψω (XΛn (ω)) lim sup log ≤ −hψω + Eψω (− log P (X0 |XBm 0 )) |Λn | P (XΛn (ω)) = −hψω + Eψω (− log P m (X0 |XBm 0 )) Lemma 2.2.26 hψω ||P m = for any m ∈ N. So 1 ψω (XΛn (ω)) lim sup log ≤ inf hψω ||P m |Λn | P (XΛn (ω)) m∈N ≤ 0, which is from Lemma 2.2.29. We are now ready to prove the non-ergodic SMB theorem. The following theorem is proved in 1D fashion in section 3.3 of [17]. Theorem 2.2.34. (Coding Rate on Stationary Source, the generalized SMB theorem) Assume the dynamical system (Ω, F, P, T ) is stationary, we have for almost every 63 ω∈Ω log P (XΛn )(ω) lim − = hψω , n→∞ |Λn | for ψ is the induced ergodic component function defined in Theorem 2.2.18. Proof. (Theorem 2.2.34) Combine Lemma 2.2.31 and Lemma 2.2.32, we have 1 ψω (XΛn (ω)) lim log = 0 P -a.e. n→∞ |Λn | P (XΛn (ω)) Also we have, from Lemma 2.2.30 we have 1 1 lim − log ψω (XΛn (ω)) = lim Hψω (XΛn ) , hψω . n→∞ |Λn | n→∞ |Λ | n So we have the result 1 lim − log P (XΛn (ω)) = hψω . n→∞ |Λn | 64 2.3 The Ideal Compression Algorithm Like many compression algorithms, we predict the pixel value by looking at its former ` decoded neighboring pixels. The set Bm = (< `) ∩ Λ`m as we defined before can be considered as the prediction window of size m for pixel location ` ∈ Λn (see Figure 2.2.2). In the end of this section, Corollary 2.3.7 will tell us that a successful algorithm can be built by estimating, P (X` |XBm ` ∩Λ ). In other words, it uses the local n information, Bm , to achieve the theoretical lower bound for an image. More exactly, we would like to show −1 X lim lim log P (X` |XBm ` ∩Λ )(ω) n = hψ(ω) , m→∞ n→∞ |Λn | `∈Λn or, in terms of the Markov approximation (see Definition 2.2.22), −1 lim lim log P m (XΛn )(ω) = hψ(ω) , m→∞ n→∞ |Λn | under the stationary condition in the sense of 1. L1 (P ) and 2. P -a.e. when (a) ψ(Ω) is countable, or (b) there is a ξ ∈ (0, 1) s.t. for any m ∈ N, we have P (X0 |XBm )(ω) > ξ for P -a.e. ω ∈ Ω. 2 About the condition 2(a), even thought X has finite alphabets, Ω = X Z is not 65 a countable set (since {0, 1}N has uncountable many elements). Also, note that a easy situation to achieve condition 2(b) is when one assumes there is always certain amount of noise in the image, where every value has chance to appear in any location of image. 2.3.1 The L1 (P ) Convergence Remember that hψ (ω) is the ideal entropy rate for every ω ∈ Ω defined as in Lemma 2.2.30: 1 hψ (ω) , n→∞ lim Hψ(ω) (XΛn ) |Λn | And from Lemma 2.2.24, we defined the entropy rate for the Markov approxima- tion as  −1 log P m (XΛn )(ω) converges.  Eψω (− log P (X0 |XBm )) , if limn→∞   m |Λn | h (ω) ,  log |X |   , o.w. −1 −1 such that limn→∞ |Λn | log P m (XΛn )(ω) = hm (ω) P -a.e. Since limn→∞ |Λn | log P m (XΛn )(ω) is convergent P -a.e. (Lemma 2.2.24), our objective is equivalent to show lim hm (ω) = hψ(ω) m→∞ in the sense of L1 (P ). 66 Theorem 2.3.1. Under the stationary condition, we have hm → hψ in L1 (P ) as m → ∞. Proof. Since hm (ω) ≥ hψ (ω) for any ω ∈ Ω, we can calculate EP |hm − hψ | = EP (hm ) − EP (hψ ) Thm 2.2.20 EP (Eψ (− log P (X0 |XBm 0 )) − hP = Thm 2.2.18(b) EP (EP (− log P (X0 |XBm 0 )|ψ) − hP = = EP (− log P (X0 |XBm 0 ) − hP → 0 as m → ∞, which is from Lemma 2.2.3. One the other hand, the P -a.e. convergence of hm is not as easy to see, we will discuss that in the next subsection. 67 2.3.2 The P -a.e. Convergence Define   − log P (X0 |XBm )(ω) , if P (XBm )(ω) > 0   fm (ω) ,   0  , o.w. for ω ∈ Ω, m ∈ N. Note that {ω ∈ Ω|P (XBm )(ω) = 0} is a measure 0 set with respect of P . We will see that, in Theorem 2.3.4, the P -a.e. convergence of hm to hψ can be achieved when lim Eψ(ω) (fm ) = Eψ(ω) ( lim fm ) P -a.e. (2.3.1) m→∞ m→∞ The next two lemmas prepares Theorem 2.3.4. And we will discuss two conditions that achieve Equation 2.3.1. Lemma 2.3.2. We have lim fm = − log P (X0 |X(<0) ) m→∞ P -a.e. and Pδ a.e. for δ ∈ ψ(Ω). Proof. We will show this with martingale theory. Define Fn , σ(∪nm=1 XBm ), thus Fn ⊂ Fn+1 for all n ∈ N, i.e. {Fn } is a filtration on(Ω, F, P ). Combine this with the fact that 1X0 =x0 is an integrable random variable for any x ∈ X on Ω, we have 0 lim E(1X0 =x0 |Fn ) = E(1X0 =x0 |F∞ ) a.e. n→∞ 68 by martingale theorem, where F∞ , σ(∪∞ n=1 Fn ) = σ(X(<0) ). That is, lim P (X0 |XBm ) = P (X0 |X(<0) ) a.e. m→∞ So for any x0 ∈ X s.t. P (X0 = x0 |X(<`) ) > 0, we have 1 1 lim log = log a.e. m→∞ P (X0 |XBm ) P (X0 |X(<0) ) Thus 1 1 P ( lim log = log ) ≥ P (ω : P (X0 |X(<0) )(ω) > 0) n→∞ P (X0 |XBm ) P (X0 |X(<0) ) = 1 − P (ω : P (X0 |X(<0) )(ω) = 0) = 1, since P (ω : P (X0 |X(<0) )(ω) = 0) = P (ω : P (X(≤0) )(ω) = 0) = 0. So we have proved the convergence in P -a.e. sense. Then we can use the fact Pδ  P to imply the convergence is also Pδ -a.e. The following lemma is modified to 2D from lemma 8.6.2 of Gray’s book [16]. The basic idea is that: When you conditioning on infinite past, you are actually conditioning on the ergodic component. 69 Lemma 2.3.3. The mapping ψ is measurable with respect to σ(X(<0) ). In particular P (X0 |X(<0) ) = P (X0 |X(<0) , ψ) = ψ(X0 |X(<0) ). Proof. Note that the second equality is just from the definition of ergodic decompo- sition. And the first equality is a direct result of ψ ∈ σ(X(<0) ). To see this, first note that from the Birkhoff ergodic theorem on images (Theorem 2.2.19), for any F ∈ F we have for P -a.e. ω ∈ Ω, 1 X lim 1F (T −` ω) = EP (1F |I), m→∞ |Bm | `∈Bm since {Bm }m∈N is a special averaging sequence. Thus we have 1 X lim 1F (T −` ω) = ψω (F ) m→∞ |B | m `∈Bm P -a.e. Next, using the fact that the probability ψ can be determined by generating set C, we can show ψ is σ(X(<0) )-measurable by choosing an event G ∈ G, say G is σ(XΛkn )-measurable, where k ∈ Z2 and n ∈ N. Without loss of generality, we set k = (0, 0). Set an , (n + 1, n + 1) ∈ Z2 , we have for P -a.e. ω: ψω (G) = ψω (T −an G) 1 X = lim 1T −an G (T −` ω), m→∞ |B | m `∈Bm which is a limit of σ(X(<0) ) measurable functions since T −an F0 is σ(X(<0) )-measurable and 1A (T −` ω) is σ(X(<0) )-measurable if A is σ(X(<0) )-measurable and ` < 0. Thus we have proved that for any event in generating set ψω (G) is σ(X(<0) )-measurable. 70 This fact implies that ψ is σ(X(<0) )-measurable. Theorem 2.3.4. If limm→∞ Eψ(ω) (fm ) = Eψ(ω) (limm→∞ fm ) P -a.e., we have −1 lim lim log P m (XΛn )(ω) = hψ(ω) m→∞ n→∞ |Λn | P − a.e. −1 Proof. Note that limn→∞ |Λn | log P m (XΛn )(ω) = Eψ(ω) (fm ) for P -a.e. ω ∈ Ω, so we only need to prove limm→∞ Eψ(ω) (fm ) = hψ(ω) P -a.e. Given an ω ∈ Ω with ψ(ω) = δ and P (ψ −1 (δ)) > 0, we have lim Eδ (fm ) m→∞ = Eδ (m→∞ lim fm ) Lem. 2.3.2 Eδ (− log P (X0 |X(<0) )) = Thm. 2.2.18(b) E(− log P (X0 |X(<0) )|ψ = δ) = Lem.=2.3.3 E(− log ψ(X0 |X(<0) )|ψ = δ) = Eδ (− log Pδ (X0 |X(<0) )) = hδ . Next, for the rest of the ω, say A = {ω ∈ Ω|P (ψ −1 (ψ(ω))) = 0}, we will prove it is a 2 P measure 0 set. Since Ω = X Z is countable, X P (A) = P (ω) ω∈A 71 and because of ω ∈ ψ −1 ψ(ω), we have P (ω) ≤ P (ψ −1 ψ(ω)) = 0 for ω ∈ A. Hence P (A) = 0. The next two lemmas give two cases that are sufficient for Theorem 2.3.4. Lemma 2.3.5. If ψ(Ω) has only countably many elements, then lim Eψ(ω) (fm ) = Eψ(ω) (m→∞ m→∞ lim fm ) for P -a.e. ω ∈ Ω. Proof. Since ψ(Ω) has only countably many elements, we can write X P (F ) = qη Pη (F ) η∈ψ(Ω) for any F ∈ F, where qη , P (ψ −1 (η)). For any ω s.t. P (ω) > 0, define δ = ψ(ω). Note that qδ > P (ω) > 0. We will prove this by stating that {fm }m∈N is uniformly integrable with respect to probability Pδ . This implies the desired property for any ω s.t. P (ω) > 0, which means the property 72 holds for P -a.e. ω ∈ Ω. Given k ∈ N, calculate Z fm dPδ ≤ (k + 1)Pδ (fm ∈ [k, k + 1)) fm ∈[k,k+1) ≤ (k + 1)Pδ (P (X0 |XBm ) ∈ (2−k−1 , 2−k ]) ≤ (k + 1)Pδ (P (X0 |XBm ) ≤ 2−k ) {Pδ (P (X0 |XBm ) ≤ 2−k |XBm = xBm )Pδ (XBm = xBm )} X = (k + 1) xB m ∈χBm {Pδ (P (X0 |xBm ) ≤ 2−k |xBm )Pδ (xBm )} X = (k + 1) xBm :Pδ (xBm )>0 Note that we can ignore those xBm with P (xBm ) = 0, since P (F |xBm ) = 0 for any F ∈ F and those term will vanish. Also observe that for those xBm s.t. Pδ (xBm ) > 0, we have P (xBm ) > 0. Let’s calculate the the first term inside the summation: P (XBm ∪0 ) Pδ (P (X0 |xBm ) ≤ 2−k |xBm ) = Pδ ( ≤ 2−k |xBm ) P (xBm ) = Pδ (P (XBm ∪0 ) ≤ P (xBm ) · 2−k |xBm ) qη Pη (XBm ∪0 ) ≤ P (xBm ) · 2−k |xBm ) X = Pδ ( η∈ψ(Ω) ≤ Pδ (qδ · Pδ (XBm ∪0 ) ≤ P (xBm ) · 2−k |xBm ) P (xBm ) −k Pδ (xBm ) > 0 Pδ (qδ · Pδ (X0 |xBm ) ≤ · 2 |xBm ) = Pδ (xBm ) P (xBm ) −k = Pδ (Pδ (X0 |xBm ) ≤ qδ−1 · 2 |xBm ) Pδ (xBm ) X = Pδ (x0 |xBm ) P (x ) x0 ∈X :Pδ (x0 |xBm )≤qδ−1 P (xBm ) ·2−k δ Bm P (xBm ) −k ≤ qδ−1 · 2 |X |. Pδ (xBm ) 73 So, Z P (xBm ) −k {qδ−1 X fm dPδ ≤ (k + 1) · 2 |X |Pδ (xBm )} fm ∈[k,k+1) xBm :Pδ (xBm )>0 Pδ (xBm ) {qδ−1 · P (xBm ) · 2−k |X |} X = (k + 1) xBm :Pδ (xBm )>0 ≤ (k + 1)qδ−1 2−k · |X |. Thus we can check the uniformly integrability of fm by Z ∞ Z X fm dPδ = fn dPδ fm >k i=1 fm ∈[k+i,k+i+1) ∞ (k + i + 1)qδ−1 2−k+i+log |X | , X ≤ i=1 which converges to 0 as k → ∞ uniformly in m. Lemma 2.3.6. If there is a ξ ∈ (0, 1) s.t. for any m ∈ N, P (X0 |XBm )(ω) > ξ for P -a.e. ω ∈ Ω , then lim Eψ(ω) (fm ) = Eψ(ω) (m→∞ m→∞ lim fm ) for P -a.e. ω ∈ Ω. Proof. Let Ω0 be the collection of ω ∈ Ω s.t. limm→∞ fm (ω) is convergent and P (X0 |XBm )(ω) > ξ. Note that P (Ω0 ) = 1 by Lemma 2.3.2 and the assumption here. For ω ∈ Ω0 , since P (X0 |XBm )(ω) > ξ > 0, we have fm ≤ − log ξ (note that the condition also implies P (XBm (ω)) > 0). Thus, we have fm ≤ − log ξ P -a.e. Thus, let δ = ψ(ω), by the Lebesgue dominate convergence theorem, we have lim Eδ (fm ) = Eδ ( lim fm ), m→∞ m→∞ 74 which proves the lemma. 2.3.3 How to Build an Ideal Compression Algorithm Let’s summarize the previous results in this section by the following corollary. Corollary 2.3.7. (The Ideal Algorithm) Under the stationary condition, we have −1 X lim lim log P (X` |XBm ` ∩Λ )(ω) n = hψ(ω) m→∞ n→∞ |Λn | `∈Λn in the sense of 1. L1 (P ), and 2. P -a.e. when ψ(Ω) is countable, or there is a ξ ∈ (0, 1) s.t. for any m ∈ N, P (X0 |XBm )(ω) > ξ for P -a.e. ω ∈ Ω Roughly speaking, this corollary tells us that, under light conditions, with large context, the compression algorithm that build upon P (X0 |XBm ∩Λn ) will approach the optimal entropy rate. Thus, one can focus on building a good estimator of 2 P (X` |XBm ` ∩Λ ). One idea is to estimate this empirically as follows: For any A ⊂ Z , n assume {(y u , yAu )}ru=1 to be an i.i.d. sample of random variable (X, XA ), and we define Pˆ r , the conditional empirical distribution of P with the library of size r , through 75 indicator functions as follows: P r  u}  u=1 I{x=yu ,xA =yA Pr    P r I u , if u=1 1{xA =yAu } 6= 0 Pˆ r (X = x|XA = xA ) = u=1 {xA =yA }   0 , o.w.,   where I is an indicator function. Then for any m, n, we define P r (X` = x|XBm ` ∩Λ ), n ` ∩Λ ) to be: the estimator of P (X` = x|XBm n  Pˆ r (X  = x|XBm ) , if Bm ∩ Λn = Bm   P r (X` = x|XBm ` ∩Λ ) n = (2.3.2)   1   |X | , o.w. By estimating with empirical distribution, we can expect that the larger amount of samples, r, we have, the better our result is. That is, “the more you know the better you get.” The obvious con for doing this is that we will need enormous amount of r to get a reasonable result. We will device a feasible way of utilizing the idea described here in the next section. 76 2.4 Application In this section we will introduce a compression scheme, called the comparison-based image compression algorithm (CBIC), which utilize the ideal algorithm idea described in subsection 2.3 to compress a real world image. The algorithm is introduced in Subsection 2.4.1, and we examine its performance in Subsection 2.4.2. Note that it is a valid image compression algorithm; however it is very computationally expensive in view of nowadays computer. The goal here is to build an algorithm that is able to scale up the context size, so as to estimate the minimum bits per pixel for an image. And compare the result with the best image compression algorithm today, CALIC [41]. A (black and white) image is set to have 256 depth of illumination, that is, X = {0, 1, · · · , 255}. The image we want to code is called the original image, and the training image set are called the library images. 2.4.1 Comparison-Based Image Compression Algorithm Instead of estimating P (X` |XBm ) with indicators, like we mentioned in Subsection 2.3.3, we use a kernel estimator to build the probability. The kernel estimator we use is a kind of a Laplace distribution 1 Kerp (z) , P|X | p|z| z˜=−|X | p|˜z| 77 Figure 2.4.1: The kernel we use to approximate empirical distribution of P (X0 |XBm ). on z ∈ Z s.t. |z| ≤ 255. for some p ∈ (0, 1). Figure 2.4.1 gives the plot of this kernel. 1−p |z| Note that Kerp (z) ≈ 1+p p . When p → 0, the kernel becomes the indicator function like 1{z=0} , and when p → 1, it becomes the uniform distribution. For easier demonstration, set the target pixel of original patch locates at (0, 0) 0 and its context, or say, neighbors, locates at T = Bm ⊂ Z2 (as shown in Figure 2.2.2). For m = 2, the corresponding patch is 3 × 5; for m = 3 is 4 × 7; m = 4 is 5 × 9. Note that m = 4 is the biggest context we will set in this experiments. To emphasize the importance of the context that near the target pixel, we will partition the context into several tiers. Set the tier of context Tb ⊂ T , b = 1, · · · , B, such that {Tb }B b=1 partitions T . For a 3 × 5 patch, we set B = 3; for a 4 × 7 patch, B = 4; for a 5 × 9 patch, B = 6. For Tb with lower value of b ∈ {1, 2, · · · , B} means it is the tier that closer to the target pixel location. For each tier we apply different α-weight 78 (a) A 3 × 5 patch (b) A 4 × 7 patch (c) A 5 × 9 patch Figure 2.4.2: A example of α-weight corresponding to different tier and sizes of patch. Note that the weight is not yet been normalized and is shown proportionally. For example: (a) For a 3 × 5 patch, T1 has a weight proportion to 3, T2 has a weight proportion to 2; T3 has a weight proportion to 1, and B = 3 in this case. Note that a notation × indicates the target pixel location for each patch. α = {αb }B b=1 , where B X αb |Tb | = 1 b=1 See Figure 2.4.2 for the information about the α-weight. The value of the original patch will be represented as a pair (·, ·) ∈ X × X T , where the first entry is called the target pixel value, and the second entry is called the context value. Before coding the patch, we first remove the illumination of both original patch and library patches. The mean of a patch is calculated with respect to α weight, then rounded to the closest integer. Note that every patch has its own mean. In the end we get the (value) centered original patch (y, x) and centered library patches {(yi , xi )}N i=1 , where N is the size of library (patches). Also we record the α-weighted 79 mean, µy , for the original target pixel. Set the following parameters: 1. r ∈ [0, ∞): A parameter for calculating distance between two context values x and xˆ, where a distance is defined as B X b−r ||xTb − xˆTb ||L1 dr (x, xˆ) , PB . b=1 b=1b−r |T | b 2. T ∈ (0, ∞): This parameter gives us more control of the distance between two context. 3. K ∈ N: Essential sample size determines the amount of library patches which are similar to the original one. Given parameters (p, r, T, K), we estimate P (y|x) out of library {(yi , xi )}N i=1 as the following (simplified) procedures: 1. Find the suitable library set: (a) Compare the α-weighted L1 -norm between x and xi for i = 1, · · · , N. Say dir = dr (x, xi ) for i = 1, · · · , N . (b) Choose K smallest values among{dr (x, xi )}N (k) K i=1 , say {dr }k=1 , and find the corresponding target pixel value of the library patches, say {y (k) }K k=1 . We will use these samples corresponds to {d(k) }K k=1 to build the empirical dis- tribution in the next step. 80 2. Calculate the the following density for T > 0 : (k) K ˜N X e−dr /T Kerp (y − y˜(k) ) P (y|x; α, p, r, T, K) = PK (j) −d /T e r k=1 j=1 where y˜(k) , max(0, min(255, y (k) + µy )). Then P˜ N (y|x) is the approximated conditional probability to P (y|x). Remember that µy is the α-weighted mean we recorded while standardizing the original path (y, x). Note that the density is consisted of kernel estimators weighted by context similarities. One would expect that when N → ∞, d(k) become 0 for all k = 1, · · · , K. And when the essential sample number K grows, the optimal p would goes to 0, which indicates that P˜ N is an approximate of the empirical conditional distribution. The whole procedure is illustrated in Figure 2.4.3. Assume there are L interior pixels (that is, pixels that have complete Bm neigh- bor), the total bits that need to be transmitted for interior pixels is L log2 P˜`N (y` |x` ; α, p, r, T, K), X f (α, p, r, T, K) = − `=1 where P˜`N is P˜ N specified by pixel `’s neighbors x` . Ideally, one would like to find the parameters α, p, r, T, and K, that minimize the total bits. Since sending these parameter won’t add many bits. Realistically, 81 Figure 2.4.3: Demonstration of CBIC algorithm. 82 we can reduce the problem to just optimize the total bits with p, r, and T . First, the parameter α is used to centering the data, so that we can use the library more efficiently. That is, without the illumination, we leave only the texture of the patches, which makes the matching procedure easier. (We have tried about standardize the patch with its standard deviation but this adds more burden computationally and didn’t gain more success.) If α is chosen badly, we can still overcome it with enormous library. So the role of α should have less impact on the result. In the experiment, we will simply choose α-weight as in Figure 2.4.2. Second, as for essential sample size K, the bigger K means more samples on estimating the empirical distribution, but it also means the quality of sample will drop since total library size N is limited. However, the quality of the sample can also be controlled by r and T , so higher number of K won’t affect much if we can choose r and T freely. The task left is how we optimize the rest of the parameters: p, r, and T . Note that, for different locations of target pixels, the algorithm use different set of K samples from library to determine the empirical probability. The size of total, N , is often enormous. For example, the highest number we use during the experiment is N = 80 million of 5 × 9 patches. It’s impractical when we try to optimize the total bit since every time we choose a triple (p, r, T ), we have to pull K nearest patch out of N library for each (interior) pixel, which requires expensive computation resources. The solution is that we choose κ nearest library patches for each pixel, in the sense that the patches have smaller α-weighted L1 norms between neighbors, where κ > K; often we choose κ = 1000. Note that, again, every different locations of target pixels have its different κ library patch set. To choose κ library for a pixel, we simply use the same α-weight as we did when calculating the mean. Then it is possible to optimize the parameter r, T , and p. In practical, we iteratively optimize between r, T , and 83 p until they converge. Theoretically, function f is not convex with respect to any of the three parameters. However, we find no trouble as we just ask the computer to minimize with respect of either one with brute force. In fact, f shows convexity with respect to p in every our experiments. Of course, one can always use larger κ to downgrade the effect of different α. To have an idea of the parameter, for “Lena” image that we will use in the next subsection, we get the estimates of p ≈ 0.544, r ≈ 1.7, and T ≈ 1.2. A heap based partial sorting algorithm is used in our implementation, as we extract K nearest contexts out of N library. The procedure follows by sorting first K distances out of N distances, which forms a sorted partition of size K, then insert the unsorted N − K distances into the sorted partition and remove the largest element in the partition. It takes O(N log K) times to do so. For users to be able to decode the image, we need to transfer the information about the edge pixels of the original image. The decoding procedure for 2 × 3 patch is demonstrated in Figure 2.4.4. It is similar for a context model that has maximum patchsize 3 × 5, 4 × 7, 5 × 9, and so on. For example, as we code with context model with maximum 3 × 5 patches, which needs 2 columns on each of the left and right, and 2 rows on top of the original image (these are called edge pixels). We use the same method on 1 × 3, 3 × 1,2 × 3, and 3 × 2, patchsizes of context model to code it. To code edge pixels, with the same κ size of library for each edge pixel, we optimize T and p corresponds to the edge pixel while using parameter r acquired during the optimization for the interior pixel. One should understand that when the size of image is large, the bits that generate from edge pixels are marginal and will not affect overall bpp much. 84 (a) (b) (c) (d) (e) Figure 2.4.4: Raster scan as decoding the image from an algorithm built on context model with maximum patchsize 2 × 3. The procedure start from (a): Given the first 3 top-left corner pixel values, using a context model on 1 × 2 patch to get the first row of image (b). Next, apply a context model on 1 × 2 patch to get left and right column of image. Then use that of 2 × 3 to get the second row, and continue to obtain the whole image (e). 85 Finally, when P˜`N fail to give a good estimate, say P˜`N < c for some original pixel `, we replace it with location information and record the pixel value directly. For ¯ + log2 256, where L example, when the illumination depth is 256, c = log2 L ¯ is the total number of pixel of an image. That is, when recording P˜`N costs more than the c, we simply replace it with c. So the total bits should be L log2 P˜`N IP˜ N ≤c + cIP˜ N >c . X − ` ` `=1 When sending the code, we need to add the following bits: • 3 pixel values for upper left corner, 8 bits each. • 3 parameters (T , p, r) for interior bpp, 32 bits for T and r, 8 bits for p. • 2 parameters (T , p) for edge bpp, 32 bits for T , 8 bit for p. • 2 parameters for image dimension, 16 bits each 2.4.2 The Performance The performance of an image compression algorithm can be measured by calculate the bits per pixel (bpp, the total file size divided by number of pixels). We prepared 4 test images as shown in Figure 2.4.5. These images are 8-bit image (28 possible pixel values), so the worst possible bpp result are 8 bits. Modern lossless compression algorithm often has ability to compress them to 3-5 bpp. The “big library” are 647 portraits images from internet (converted to 8 bit images), which comprise about 300 86 Lena Jetplane Goldhill Barbara Figure 2.4.5: 4 test images: Lena, Jetplane, Goldhill, and Barbara. Goldhill has size 720 × 576, and others are 512 × 512. million 4 × 7 image patches, and we will randomly choose several amount of patches in it to build the library for our experiments. Due to the limitation of the modern computer, for CBIC algorithm, we choose the essential sample number K = 200 in general (we tried use higher amount of K, however, the result only improve a little). Figure 2.4.6 compares different patchsize and library sizes while compressing “Lena” image. The library for each experiment is randomly selected from the big library. Roughly speaking, when library size is small, 5 × 9 patchsize has slightly worse result than that from 4 × 7. But when the library size is high, like 40 millions, the compression rate improves from 5 × 9. This is because it is harder for 5 × 9 patches to find quality samples from smaller library. When library size is up to 40 million, a larger patch size gives a better result. In this case, the gain between 5 × 9 and 4 × 7 are marginal, compare to that between 4 × 7 and 3 × 5, which indicates the importance of neighbors diminishes along patchsize. Figure 2.4.7 compares the result of CBIC, the proposed method, with CALIC (see [41]) and Jpeg-LS (see [25]). Note that even though we use the portrait library, the image with different subjects (e.g. Jetplane or Goldhill) still gives compatible result with CALIC. This suggests that a successful compression scheme can rely only 87 Figure 2.4.6: The bpp result of CBIC algorithm for “Lena” image with different size of patches and different size of library. In these experiments, we choose essential sample amount K = 200 to build the empirical distribution. 88 Figure 2.4.7: Comparison between CBIC and other famous lossless compression al- gorithms. The following information in the parenthesis shows the number of library and the patchsize used to generate the result here: Lena: (80M, 5 × 9); Jetplane: (80M, 4 × 7); Goldhill: (40M, 4 × 7); Barbara: (20M, 4 × 7). on the local information of an image, rather than the global one. Also notice that our algorithm shows good result on compressing portrait images “Barbara” with only 20 million library (also, from Table 2.4.6, 20 million library patches also gives good result on “Lena”) , and this might indicate it will be efficient to find patches within the same type of photos. For reader who followed through the article, this might be an evidence that these portrait images have similar ergodic components. In Figure 2.4.8, we reconstruct the Lena image with pixel value that has highest probability the estimator gives (i.e. the MAP, maximum a posterior probability, esti- mators). That is, each pixel intensity is the one that has highest value of P˜`N (y` |x` ), note that x` is the actual pixel values of the neighborhood of `th pixel8 . The estima- tor is obtained from our method with 4 × 7 patchsize, with 5 million library patches. Note that 5 million library patches is considered as the smallest amount of library in 8 Do not confused this with lossy compression. This used the original image context, hence it is not a lossy compression result. 89 our experiment. One can see that a lot of detail is captured by our models, includes some of the textures of the hat and the hair. Despite all the great performance of CBIC, the main drawback of the algorithm is the maximum memory usage and expensive computation resource requirements. While we implement the algorithm with Python combine with C that handles the most intensive sorting algorithm, it would take several hours to encode, and even more time to decode. For example, we used Amazon EC2 to enhance the computation, which gives 16 or more cores for parallel computing, yet it still take 8-10 hours to encode a small 512 × 512 image with 40 million library. Why do we built such an algorithm? The main purpose is to have an image compression algorithm that is able to scale up the context size, and to achieve the optimal coding rate, so that we can see how a state of the art algorithms performed. As one can see from Table 2.4.7, we are only a little better (about 0.05 to 0.15 bpp) than one of the best scheme. Also, when we scale up the context size (see Table 2.4.6), the gain in performance is diminished. This suggests that the state of the art lossless image compression algorithm is already near optimal and it might be hopeless to find another scheme that is overwhelmingly better than the one we have today. 90 Figure 2.4.8: Top: Reconstruct the image with MAP estimators. Bottom: The original picture. Chapter Three Robust Generalized Clustering 92 3.1 Introduction A traditional way to summarize a data is through regression. Over years of research, many generalizations have been developed, some of which search for multiple struc- tures associated with different parts of the data. For example, Breiman proposed [8] the regression trees algorithm (CART, classifications and regression trees) that finds histogram-liked structures to explain the data in each of many cells that are derived by a recursive splitting of the variable associated with each dimension. Freid- man introduced MARS, the multivariate regression splines [14], using hinge function as a basis to learn piece-wise linear functions. The first stage of MARS is to over- fit the data with hinge functions, where hinge functions with the smallest sum of squared residuals are added iteratively. Then a backward deletion method is applied to obtain the final result. Later, Breiman combined the regression trees idea with hinge functions, and proposed an algorithm called the hinge hyperplane method, [7], where a relatively computationally inexpensive hinge-finding algorithm is introduced for learning optimal regression trees. Another example is SVR, or “support vector regression,” proposed by Vapnik, et. al. [40]. SVR seeks the flattest nonlinear struc- ture that fits most of the data, by solving a quadratic programming problem. A clear development of SVR can be found in [39]. The above-mentioned methods accommodate what we call “multiple instance data,” in which there are multiple models, each designed to fit a different subset of the data. In these cited examples, there is a dependent variable and the ultimate goal is regression. We will also assume that multiple models are needed to explain the data, but in something more akin to the clustering problem than the regression 93 problem—e.g. see Figure 3.1.1. The best-known approach is probably the K-means algorithm: Given data {zi }ni=1 ∈ Rd , K-means (see [27] for an early reference) seeks m points {cj }m d j=1 ⊂ R that minimize the loss function, f : n min kzi − cj k2 X f (c1 , · · · , ck ) = j=1:m i=1 A heuristic algorithm is applied: Firstly, arbitrarily choose m points in Rd and as- sociate (cluster) the data that is closest to each point. Second, calculate the mean of each cluster and repeat the first step with these means. Continue recursively. Apparently, this would not be a fruitful approach to the data in Figure 3.1.1. A generalization of this idea, replacing points with hyperplanes (“K-planes”, [28]), was proposed by Manwani at 2012. However, the direct generalization quickly gets stuck stuck in local minimums and is highly sensitive to noise, as the author himself indi- cates in his paper. This might be the reason for why this simple generalization has found little favor in literature. Manwani’s remedy is to add a local penalty function, kzi − µj k2 , to the loss function to ensure that the points in each cluster are close together. This approach can also be viewed as an instance of the EM algorithm [10], in which residuals are modeled as a Gaussian mixture with prior. We will compare our method with an orthogonal version of EM in Section 3.4. Not surprisingly, we will find that EM with Gaussian mixtures is sensitive to noise and outliers. Our approach, which will be introduced shortly, offers a way to learn multiple structures using a very simple loss function, similar to that of the K-means algorithm, but with substantial robustness. The robustness issue already appears in the context of ordinary linear regression. Some noise or outliers can quickly undermine the result. This is the inevitable nature 94 Figure 3.1.1: “Triangle data” of a quadratic loss function, which nevertheless enjoys many excellent theoretical (and even practical properties) under a Gauss-Markov type condition. But most real- world data is not well described by these assumptions. Many methods and resources have been dedicated to remedying this. A good treatment of many of these can be found int the book by Sen and Srivastava [38]. One fruitful line of attack is called, quite generally, “robust regression.” Whereas the objective is still to minimize the sum of residuals, the focus is on using different functions of the residuals. The first approaches in this direction date back to the 18th century, with the idea to replace squared residuals by the absolute value of residuals, “least absolute value (LAV) regression,” proposed by Boscovich [5] and later studied by Laplace [22]. (A comprehensive review of LAV regression can be found in [11].) To demonstrate the difference between LAV and traditional regression, consider the data set displayed in Figure 3.1.2. We generated 100 points from a distribution concentrating near a straight line, and then added 50 outliers uniformly distributed in the rectangle. How does β affect the minimizer of n d(zi , `)β , X f (`) = i=1 95 where d(z, `) is the Euclidean distance from a point z to a line `? When β = 2, the minimizer is the line that passes through the mean of the data and has the largest variance of the projected data (i.e. the eigenvector of the sample covariance matrix corresponding to the largest eigenvalue). Here, we are using “orthogonal distance regression” [4], as defined in the formulation of principle component analysis (PCA). Notice that the optimal line when β = 1 is less affected by outliers than the corresponding line when β = 2. (Traditional LAV applies to the vertical residuals, as in traditional regression, rather than the orthogonal distance used in PCA.) Pn Figure 3.1.2: Minimizers of f (`) = i=1 d(zi , `)β when β = 1 and β = 2, where ` is a straight line in R2 . A unified approach is Huber’s M -estimator [20], which determines regression pa- rameters by minimizing the sum of error function n X ρ(i ), i=1 where i , i = 1, 2, · · · , n, are the residuals, andρ is a convex function. Note that when ρ is the square function the estimator is the ordinary regression estimator; LAV regression is recovered by taking ρ to be the absolute value. As for non-convex ρ, 96 Huber wrote “If ρ is not convex, the estimators... will, in general, no longer converge toward some constant. Apparently, one has to impose not only local but also some global conditions on ρ in order to have consistency.” In contrast to Huber’s ideas, we propose to explore non-convex error functions. We forgo any reasonable notion of consistency, but we gain a strong sense of robustness, and in some cases computational advantages as well. Specifically, we consider residual Pn functions of the form i=1 d(zi , `)β , where d(zi , `) is Euclidean distance and β ∈ (0, 1). When β ∈ (0, 1) we call d(z, `)β a sparse distance. For sparse distances, exact minimums can sometimes be computed, and oftentimes approximated (3.2.2, 3.4). As we shall see, sparser distances (smaller values of β) are more robust. As an illustration, in the data shown in Figure 3.1.3 there 70 outliers out of 100 data points. The 30 “good” data points were generated from a noisy line and the 70 outliers from the uniform distribution on the rectangle. As β is made smaller, the minimizing line “migrates” to a line that is nearly indistinguishable from the one that generated the 30 data points. We are interested in extensions to “multiple instances”: problems for which two or more shapes (e.g. lines, or lines and planes) might be the most appropriate summary. This calls for an extension of the loss function. Given data ~z = {zi }ni=1 ⊂ Rd1 and a class of models M (e.g. all one and two-dimensional linear manifolds in Rd ), 1 We use ~z, ambiguously, to represent the vector (z1 , z2 , . . . , zn ) and the (possibly multi) set {z1 , z2 , . . . , zn }. 97 Figure 3.1.3: Left: the data includes 70 points generated from the uniform distribu- tion. Right: The effect of β on the minimizer of f (`) = ni=1 d(zi , `)β , over all lines, P `, in R2 . we propose an unsupervised clustering algorithm through minimization of the loss function n λX f (L) , γ(L) + min d(zi , `)β , (3.1.1) n i=1 `∈L where d is, again, Euclidean distance and L ∈ 2M is a set of instances (or “structures” or “models”). Obviously, there needs to be a penalty on the number of instances: γ(L) is a monotone increasing function of the number of models in L. The minimization of the loss functions with sparse distance involves the “inter- polated solutions.” By interpolating we means it is fully determined by the data points that it passes through, e.g. a two-dimensional plane passing through three linearly independent points. In LAV regression, it is known that for d dimensional data (xi , yi ) ∈ Rd−1 × R1 , i = 1, 2, · · · , n (with full-rank matrix), the regression hy- perplane will pass through at least d data points, see [1]. (Some other properties of 98 the estimator can be found in Chapter 1 of [23].) Formally, the minimizer, ˆb, of n X d−1 X f (b) = |yi − xij bj |, i=1 j=0 where xi0 = 1 ∀i, defines a hyperplane, d−1 {(x, y) ∈ Rd−1 × R1 : y = ˆb0 + xj ˆbj } X j=1 that includes at least d of the data points. It is then said that the solution is inter- polated. As we shall see, the phenomenon is also present in the minimization of n d(zi , `)β , X i=1 over all ` ∈ M, where M is the set of (d − 1 dimensional) hyperplanes and zi ∈ Rd ∀i. The solution passes through at least d data points. While this result for LAV regression (or for that matter the case here, when β = 1) can be proven through the dual representation, the technique does not work for β < 1. We will treat this more general case in 3.2.2 and 3.2.3, where we will describe a “interpolated solution.” For more general planes, of dimension k < d − 1, the minimum-loss solution does not necessarily go through any data points. In this case, the k-dimensional planes that pass through k + 1 data points can only be considered as local minima. We discuss this in Section 3.2.4, and for special cases k = 0 and k = 1 we will give bounds on β which guarantee that the optimizer is interpolating. It is instructive to take a close look at the simplest example, which is the clustering of points on the line. In Section 3.3 we will examine a limited type of consistency and we will show that the estimator has 50% breakdown point, which is the highest 99 possible for a location estimator. In Section 3.4 we will develop an approximate minimizer of the (multiple-instance) loss function, which we will call the “Remove or Replace Algorithm” (RRA). We will perform numerous experiments, with different amounts of outliers and a variety of model classes, M. For “mixed” model classes—those with two or more parametrized shapes (e.g. lines and d − 1 dimensional planes), care must be taken to not unduly favor model classes with more degrees of freedom, such as a hyperplane versus a line. We will offer an extension of RRA for choosing multiple (mixed) instances from mixed model classes. In particular, we will adapt a method called minimum sum of volumes, which provides a way to select a final subset from a candidate list of multiple mixed instances. Candidates are chosen using the RRA algorithm. In our experiments, we obtain natural summaries of complex data, even in the presence of a high percentage of outliers. The approach, using sums of volumes, is similar, but not the same as, a method known as the minimum volume covering ellipsoid. The idea of clustering with minimum volume ellipsoids was first proposed by Rosen [36]; see also Barnes [2] for further developments of the algorithm. 100 3.2 Loss Function for Multiple Instances The main object of study, in this section, will be a special case of Equation 3.1.1: given β ∈ (0, 1), let n λX f (L) = |L|2 + min d(zi , `)β . (3.2.1) n i=1 `∈L where L ⊂ M and d(zi , `) is the (orthogonal) distance2 from zi ∈ Rd to the model instance ` ∈ M. Before minimizing Function 3.2.1, we standardize the data to achieve scale and zi }ni=1 , translation invariance. That is, for a raw data {ˆ zˆi − z¯ˆ zi = q P , 1 n kˆ z i − ¯ˆk2 z n i=1 where z¯ˆ = 1 Pn Pn n ˆi i=1 z and k · k is the usual L2 norm. Thus we have i=1 zi = 0, and 1 Pn n i=1 kzi k2 = 1. Our main goal will be to prove that when M is the family of d − 1 dimensional planes in Rd then, generically, the minimizer of 3.2.1 passes through at least d linearly independent points in ~z (i.e. “interpolates”). Subsections 3.2.2 and 3.2.3 are devoted to the proof. Section 3.2.4 gives some partial results for the more difficult case of planes of dimension less than d − 1. Precisely, here is what we mean by : 2 We will only discuss linear model classes, but we have in mind more general shapes. In general, “orthogonal” would be replaced by minimum distance. 101 Definition 3.2.1. Fix d ≥ 2 and k ≤ d − 1. Let Mk be the set of all k dimen- sional planes in Rd , let z1 , z2 , . . . , zn be n elements of Rd , and write ~z for the set {z1 , z2 , . . . , zn }. A plane ` ∈ Mk is interpolating (or, sometimes, interpolates) rela- tive to ~z if it contains at least k + 1 linearly independent points from ~z. Before turning to these topics, we will first establish, in Section 3.2.1, a straight- forward result that will greatly simplify our subsequent analyses. 3.2.1 Reduction to m = 1 Given a minimizer of 3.2.1, and a particular model instance ` ∈ M that is a compo- nent of the minimizer, consider the subset of data points, A ⊂ ~z, which are closer to ` than to any other model instance in the solution set (for now, ignore ties). Does d(zi , ·)β ? The obvious answer, ‘yes’, is correct, but to show this we P ` minimize z∈A can not simply assume the contrary and then improve `, since changing ` will, in general, also change A. Here we give a proof, and, later, we will use this result to help delineate some key properties of an optimal solution. Definition 3.2.2. Let L∗ = {`∗1 , `∗2 , . . . , `∗m } be a minimizer of Function 3.2.1. An optimal assignment, given L∗ , is any partition of {z1 , z2 , . . . , zn }, A(L∗ ) = {A1 (L∗ ), A2 (L∗ ), . . . , Am (L∗ )}, such that z ∈ Aj (L∗ ) ⇒ d(z, `∗j )β ≤ d(z, `∗k )β ∀k = 1 : m 102 To ease notation we will sometimes write ρz (`) in place of d(z, `)β . With this conven- tion: Theorem 3.2.3. Let L∗ = {`∗1 , `∗2 , . . . , `∗m } be a minimizer of Function 3.2.1, and let A(L∗ ) be an optimal assignment for L∗ . Then for all j = 1, · · · , m `∗j ∈ arg min X ρz (`), `∈M z∈Aj (L∗ ) Proof. Suppose not. Then without loss of generality we may assume that `∗1 ∈ X / arg min ρz (`) `∈M z∈A1 (L∗ ) In which case there exists `˜1 ∈ M such that ρz (`∗1 ) > ρz (`˜1 ) X X z∈A1 (L∗ ) z∈A1 (L∗ ) ˜ = {`˜l , `∗2 , `∗3 , . . . , `∗m }, and `˜j = `∗ for j = 2, 3, · · · , m. For any L ⊂ M, Define L j Pn define h(L) = i=1 minl∈L ρzi (`). Then n h(L∗ ) = X min∗ ρzi (`) `∈L i=1 m ρz (`∗j ) X X = j=1 z∈Aj (L∗ ) m ρz (`˜1 ) + ρz (`∗j ) X X X > z∈A1 (L∗ ) j=2 z∈Aj (L∗ ) m ρz (`˜j ) X X = j=1 z∈Aj (L∗ ) n ˜ X ≥ min ρzi (`) = h(L) ˜ i=1 l∈L 103 which contradicts the optimality of L∗ . Basically, Theorem 3.2.3 says that the any assignment of ambiguous points (points that have the same distance to two or more instances in a minimizer) would not affect the result of the theorem, i.e. `∗ will still minimize its cluster. An interesting fact can be derived: When β = 2 and point model is used, there is no ambiguous point for the minimizer. This is because, if there were ambiguous points, any different assignment will result in different minimizer, which violate the property of unique solution of a (strictly) convex optimization problem. However, when β < 1, it is possible to see ambiguous points. 3.2.2 Lines in R2 A natural choice of model M is the set of straight lines in R2 , which is called the line model. The main purpose here is to show that the minimizer of n d(zi , `)β , X f (`) = (3.2.2) i=1 over the line model, given data ~z = {zi }ni=1 ⊂ R2 , is interpolating—i.e. passes through at least two points (Theorem 3.2.6). To have a better understanding of this, first consider both data point zi ∈ R1 , M = R1 . It’s easy to see that the minimizer exists since f is continuous and the minimizer must be within a bounded set [mini=1:n zi , maxi=1:n zi ]. Also, because of f 00 (x) < 0 for x ∈ / ~z, the only possible minimizer can only occur when x ∈ ~z. That is, 104 the minimizer interpolates. Our proof of that on R2 is similar to this, as one will see. Here are some notations we will use throughout the section. An element in Rd is considered as a d-dimensional column vector, and the superscript T on a vector means the transpose of the vector. The identity matrix is denoted as I. The superscript ⊥ on a space is the orthogonal complement of the space, i.e. B ⊥ = {y ∈ Rd : x · y = 0 for any x ∈ B} for any B is a subspace in Rd . Let S , R mod π, and v(θ) , (cos(θ), sin(θ))T for θ ∈ S1 . Define a function L that maps R × S to a closed set on R2 , s.t. L(a, θ) = {y ∈ R2 : y T v(θ) = a}. Thus, we define a line L(a, θ) s.t. it has normal vector v(θ) and d(0, L(a, θ)) = |a|. Note that this defines an unique line on R2 , so we can define the space of lines on R2 as L , R × S. Note that the space S has a natural metric dS (θ1 , θ2 ) = min{|θ1 − θ2 |, π − |θ1 − θ2 |}. So we can define the metric of space L as the product metric dL (L(a1 , θ1 ), L(a2 , θ2 )) , |a1 − a2 | + dS (θ1 , θ2 ) for ai ∈ R, θ ∈ S. Thus we have a metric space for lines (L, dL ). To prepare Theorem 3.2.6, we need to show that the minimizer exists (Lemma 3.2.4 and 3.2.5). Lemma 3.2.4. Function 3.2.2 is continuous on (L, dL ). Proof. For any z ∈ R2 , we only need to show that d(z, ·) is continuous on (L, dL ). 105 Pn Then i=1 d(z, ·)β is continuous on (L, dL ). Set a line `a,θ = L(a, θ) for a ∈ [0, ∞) and θ ∈ S1 , since d(·, ·) is the Euclidean distance. Note that `a,θ = av(θ) + `0,θ . since the we have d(z, `a,θ ) = d(z − av(θ), `0,θ ) = k(z − av(θ))T v(θ)k, where v(θ) = (cos(θ), sin(θ)). Thus, we can see that d(z, `a,θ ), as a function of (a, θ), is a continuous function of a, cos θ, and sin θ. This implies d(z, `x,θ ) is continuous on (L, dL ). Lemma 3.2.5. The minimizer of Function 3.2.2 exists. 1 Pn Proof. Define a d-dimensional ball Bv (¯ z ), where z¯ = n i=1 zi , and v = max kzi − z¯k. Define B , {` ∈ L|` passes through cl(B2v (¯ z ))}, where cl(·) means the closure of a set. Thus, for ` ∈ / B then f (`) > nv β . Next, if `z¯ is a line starting at z¯, we have f (`z¯) ≤ nv β for every possible direction. Thus we have min f (`) ≤ nv β ≤ min f (`). `∈B `∈B / 106 This implies min f (`) = min`∈B f (`), `∈L which means we can restrict our search of lines on B, a compact set in (L, dL ). This, along with the fact that the Function 3.2.2 is continuous, the minimizer is always existed. Note that the minimizer of Function 3.2.2 might not be unique: Consider d = 2 and 4 points forms a exact square on a plane, then there are two lines admits the same minimum value of the target function. Theorem 3.2.6. Given data points {zi }ni=1 ⊂ R2 , assume β ∈ (0, 1), the minimizer of Function 3.2.2 interpolates. Proof. (Proof by contradiction) Assume the optimal line doesn’t pass through any data points. Given any point z ? ∈ R2 on the optimal line we have zi 6= z ? for any i = 1, · · · , n. We may further assume z ? = 0 since we can shift the whole data to achieve it. Furthermore, set the optimal line, say `(θ? ), parallel to the vector (cos θ? , sin θ? ) for some θ? ∈ [0, 2π). Define a line, `(θ), that passes the origin and parallel to the vector (cos θ, sin θ) for θ ∈ [0, 2π). Then θ? should be the argument of minimum of the following function: n d(zi , `(θ))β . X g(θ) = i=1 Note that, from our setting, the optimal line `(θ? ) doesn’t pass any data point, i.e. θ? 6= θi + 2πk, for any k ∈ N and i = 1, · · · , n. For any data point zi , i = 1, · · · , n, 107 we can also write zi = ||zi ||(cos θi , sin θi ) for some θi ∈ [0, 2π). First, observe that d(zi , `(θ)) = ||zi || · | sin(θ − θi )|. For i = 1, · · · , n, define gi (θ) , d(zi , `(θ))β , thus we have gi (θ) = ||zi ||β · | sin(θ − θi )|β ∂gi (θ) = ||zi ||β β| sin(θ − θi )|β−1 cos(θ − θi ) · sgn(sin(θ − θi )), ∂θ for any θ 6= θi + 2πk, where sgn(a) , 1a≥0 − 1a≤0 . The Hessian of gi can be calculated whenever θ 6= θi + 2πk,: Hgi (θ) = ||zi ||β β{(β − 1)| sin(θ − θi )|β−2 cos2 (θ − θi ) · [sgn(sin(θ − θi ))]2 +| sin(θ − θi )|β−1 · (− sin(θ − θi ) · sgn(sin(θ − θi )) + 0)} = ||zi ||β β{(β − 1)| sin(θ − θi )|β−2 cos2 (θ − θi ) − | sin(θ − θi )|β−1 · | sin(θ − θi )|} = ||zi ||β β{(β − 1)| sin(θ − θi )|β−2 cos2 (θ − θi ) − | sin(θ − θi )|β } < 0. Thus for θ? 6= θi + 2πk, n Hg (θ? ) = Hgi (θ? ) < 0, X i=1 108 which means θ? is not the argument of minimum of g, we have the contradiction. Thus, the optimal line has to pass at least one data point, say z1 . Next, find the optimal direction θ by minimizing (d(zi , `z1 (θ))β , X gz1 (θ) = i:zi 6=z1 ,i=1:n where `z1 (θ) is a line passes z1 with direction (cos θ, sin θ). Using the same proce- dure as above: start from assuming the optimal line doesn’t touch any data point other than those equal to z1 , and get the contradiction. One can conclude that the minimizer must admits a line that touches at least two data points. One might expect to have the same result for general data in Rd —the optimal line is interpolating. However, unfortunately it is not true. The following example shows it is not necessary true. Example 3.2.7. Consider a set of data that plotted in Figure 3.2.1 z1 = (1, 0, 0) z2 = (−1, 0, 0) √ z3 = (0, 3, 0) √ 3 z4 = (0, , b). 3 Let m = 1, when b is sufficiently high, if we confine the search of the optimal line on only the line that connects two data points, the result line, say `g , will be one of 109 z1 z4 ,z2 z4 , and z3 z4 , which has the target function value of q 1 2 b2 + 3 β f (`g ) = 2 · ( q ) , 4 b2 + 3 which goes to 2 · 2β when b → ∞. Next, let’s consider another line that go through z4 and perpendicular to the x − y plane, say `r , depicted as red line in Figure 3.2.1, will have target function value 2 f (`r ) = 3 · ( √ )β . 3 So, when β ↑ 1, we have √ lim f (`r ) = 2 3 ≈ 3.46 β↑1,b→∞ lim f (`g ) = 4, β↑1,b→∞ That is, we found a line `g that has smaller value then any line that connects 2 data points. Hence, in Rd ,d > 2, the optimal line is not necessary interpolating. 3.2.3 Hyperplanes in Rd In this subsection, we will study the property of the minimizer while M is the family of Rd−1 planes in Rd , d > 2, which is called the plane model. Before we define space of planes, note that a vector r ∈ Rd1 , {r ∈ Rd |rT r = 1} can be represented as a vector of angles η = (θ1 , θ2 , · · · , θd−1 ) ∈ Θd , [0, π]d−2 × [0, 2π) using spherical coordinates. And there is a bijection v : Θ → Rd to convert from spherical coordinates to Cartesian 110 Figure 3.2.1: Counter example for the optimal line in R3 , which explaining the optimal line might not necessary go through two data points. The data is z1 , · · · , z4 , when b is high enough, the green line will be the best line (that has smaller or equal minimum target function value) among all the interpolating lines, i.e. the lines that connect two data points. However, when β < 1 and is close to 1, the red line will produce smaller target function value than that of the green one, and the red line only touches one data point. 111 coordinates. Note that the kth- elements of vector v(η) ∈ Rd , k = 1, 2, · · · , d, can be represented as Qd−k v(η)k = ( j=1 sin θj ) cos θd−k+1 , Q0 where θd = 0 and j=1 sin θj = 1 (see [19] page 645). Thus, an Rd−1 plane in Rd can be represented as follows. Given a ∈ R, and η ∈ Θ , S d−1 , where S = R mod π, a plane can be defined as a function P maps R × Θ to a close set on Rd , s.t. P(a, η) = {y ∈ Rd : y T v(η) = a}. That is, P(a, η) is a plane that has normal vector v(η) and has distance |a| to the origin. As we did in Subsection 3.2.2, we define a product metric on the plane space. The distance between two planes, D(`1 , `2 ), `j = P(aj , ηj ) ∈ L, j = 1, 2, is defined as d−1 X |a1 − a2 | + dS (θ1k , θ2k ), k=1 where θjk = (ηj )k , j = 1, 2. Thus we have a metric space (Pd , D). Given data {zi }ni=1 ⊂ Rd , the target function is n d(zi , `)β X f (`) = (3.2.3) i=1 for ` ∈ P. Note that Pd can also be defined as the collection of planes that passing through d linearly independent points. We would like to show that the minimizer of Function 3.2.3 will pass through at least d data points (Theorem 3.2.10). As before, we first need to show the minimizer of exists (Lemma 3.2.8 and 3.2.9). 112 Lemma 3.2.8. Function 3.2.3 is continuous on (Pd , D). Proof. For any z ∈ Rd , we only need to show that d(z, ·) is continuous on the product space (P, D). Given a plane `a,η = P(a, η) for a ∈ R and η = (θ1 , · · · , θd−1 ) ∈ Θd , since the distance d(·, ·) is considered to be the Euclidean distance, we have d(z, `a,η ) = d(z − av(θ), `0,η ) = k(z − av(θ))T · v(η)k, which is a sum of continuous functions of a and θ1 , · · · , θd−1 . Thus, d(z, ·) is clearly continuous on (P, D). Lemma 3.2.9. The minimizer of Function 3.2.3 exists. Proof. The proof is similar to Lemma 3.2.5. Define a d-dimensional ball Bv (¯ z ), where 1 Pn z¯ = n i=1 zi , and v = max kzi − z¯k. Define B , {` ∈ P|` passes through cl(B2v (¯ z ))}. Thus, if ` ∈ / B then f (`) > nv β . Next, if `z¯ is a plane passing through z¯, we have f (`z¯) ≤ nv β for every possible direction. Thus we have min f (`) ≤ nv β ≤ min f (`). p∈B p∈B / 113 This implies min f (`) = minp∈B f (`), p∈P which means we can restrict our search of planes that passes in B, i.e. a compact set in (P, D). This, along with the fact that the Function 3.2.2 is continuous, the minimizer is always existed. Theorem 3.2.10. Assume β ∈ (0, 1) and span(z1:n ) = Rd , the minimizer of Function 3.2.3 interpolates. Proof. Assuming there is no data points that lies on the optimal plane `. We may assume 0 ∈ ` and have normal vector (1, 0, · · · , 0) since we can shift and rotate the whole data to achieve it. Next, rewrite our target function by restricting the domain of planes that passing through 0 and allowing only the first two coordinates of their normal vectors to move. That is, for some θ ∈ [0, 2π), let `0,θ to be the plane that passes 0 and has normal vector rθ , (cos θ, sin θ, 0, · · · , 0), we search the argument of the minimum of the following function n (d(zi , `0,θ )β . X g(θ) = i=1 And we know that θ = π2 +kπ, k ∈ N, which corresponds to normal vector (1, 0, · · · , 0), minimizes the function. The following argument will show that if the optimal plane doesn’t touch any data point, the Hessian of g, say Hg , will be negative when θ = π 2 + kπ, and we will get the contradiction. π First, assume that when θ = 2 + kπ, the plane `0,θ doesn’t touch any data point. 114 For i = 1, · · · , n, write zi as (xi , yi , zi,3 · · · , zi,d ) ∈ Rd , and calculate d(zi , `0,θ ) = |zi · rθ | = |xi cos θ + yi sin θ| q xi yi = x2i + yi2 | q cos θ + q sin θ| x2i + yi2 x2i + yi2 = αi | cos(θ − θi )|, q where αi , x2i + yi2 > 0, and θi satisfies cos θi = √ x2i and sin θi = √ y2i . xi +yi2 xi +yi2 π Note that d(zi , `0,θ ) is twice differentiable only when θ − θi 6= 2 + kπ. Let’s establish this by the following claim: π Claim 3.2.11. For i = 1, · · · , n, d(zi , `0,θ ) is twice differentiable when θ = 2 + kπ. Proof. Given any i = 1, · · · , n, since we assume the optimal plane doesn’t touch any of the data point, the vector zi = zi − 0 will have the property that zi is not perpendicular to the optimal plane’s normal vector (1, 0, · · · , 0), i.e. zi · (1, 0, · · · , 0) 6= 0, π which implies xi 6= 0, and this means cos θi 6= 0, i.e. θi 6= 2 + kπ. So we conclude π that d(zi , `0,θ ) is twice differentiable when θ is in the neighborhood of 2 + kπ. Let gi (θ) , (d(zi , `0,θ )β for i = 1, · · · , n, we can calculate gi0 s derivative on the 115 π neighborhood of θ = 2 + kπ. ∂gi (θ) ∂ β = α | cos(θ − θi )|β ∂(θ) ∂(θ) i = αiβ β| cos(θ − θi )|β−1 · [− sin(θ − θi ))sgn(cos(θ − θi )]. π And we can calculate the Hessian of gi when θ is in the neighborhood of 2 + kπ : Hgi (θ) = αiβ β{(β − 1) · | cos(θ − θi )|β−2 · [− sin(θ − θi )sgn(cos(θ − θi ))]2 +| cos(θ − θi )|β−1 · [(− cos(θ − θi ))sgn(cos(θ − θi )) + 0]} = αiβ β{(β − 1) · | cos(θ − θi )|β−2 sin2 (θ − θi ) + | cos(θ − θi )|β−1 · [−| cos(θ − θi )|]} = αiβ β{(β − 1) · | cos(θ − θi )|β−2 sin2 (θ − θi ) − | cos(θ − θi )|β ]} < 0 π when θ = 2 + kπ. So we have Pn Hg (θ) = i=1 Hgi (θ) < 0 π when θ = 2 + kπ. This means the plane that passes 0 with the normal vector (1, 0, · · · , 0) is not the optimal plane. Thus we have the contradiction. There is at least one data point touches the optimal plane. Next, say the data point that touches the optimal plane is z1 , we may assume that z1 = 0, if not, we can shift the data to achieve that. As before, we may assume the π optimal plane touches 0 and has normal vector (1, 0, · · · , 0). Then, again θ = 2 + kπ should minimize g (1) (θ) = (d(zi , `0,θ )β . X i∈{1,··· ,n|zi 6=z1 } 116 Note that we excludes the data points that is equal to 0 so that g is twice differentiable π when θ = 2 + kπ. Using the same argument, we conclude that the optimal plane have to pass through another data point that is not equal to 0, say z2 . Finally, assume the optimal plane only go through q distinct data points, 1 < q < d, say z1:q and z1 = 0. And the normal vector of the optimal plane in the space of span(z2 , · · · , zq )⊥ . For j = 1, · · · , d, define ej ∈ Rd1 , such that all the components of ej are 0 except the j th component is 1. And we may assume span(z2 , · · · , zq )⊥ = span(ed−q+2 , ed−q+3 , · · · , ed )⊥ , Aq and the optimal plane’s normal vector is e1 , since we can rotate the whole data to achieve that. Thus, if we restrict the plane that has normal vector in Aq , we should find the plane with normal vector e1 minimize the target function. In another words, π if we restrict the plane that has normal vector in {`0,θ }θ∈[0,2π) , θ = 2 + kπ should minimize g (q) (θ) = (d(zi , p0,θ )β . X i∈{q+1,··· ,n|zi 6=z1:q } Note that `0,θ has normal vector (cos θ, sin θ, 0, · · · , 0) ∈ Aq . 117 Then, again, the same argument will give us the contradiction. One must understand the argument is not valid when q = d, since span(e2 , ed−q+3 , · · · , ed )⊥ is the space that span by only one vector, e1 . So there is no freedom to move the normal vector to get smaller value on the target function. Thus we have proved that the optimal plane must pass at least d data points. Combine with Theorem 3.2.3 and Theorem 3.2.6, let’s sum up the important result with the following corollary: Corollary 3.2.12. Given data {zi }ni=1 ⊂ Rd with span(z1:n ) = Rd , d ≥ 1 and assume the model class M = Pd , the space of Rd−1 plane in Rd (which can be determined by d linear independent points in Rd ). If β ∈ (0, 1), then there is a minimizer L∗ = {`∗1 , · · · , `∗m } ∈ Pd of function n λX f (L) = |L|2 + min d(zi , `)β , (3.2.4) n i=1 `∈L for L ⊂ Pd , such that `∗j passes through at least d linear independent data points of {zi }ni=1 for all j = 1, 2, · · · , m. i.e. the minimizer is interpolating. 118 3.2.4 Non-interpolating Optima: Some Partial Results What can be said about the minimizers of 3.2.4 when M = Pk (k < d, d is the dimension of the data) is the set of planes of dimension k − 1 (or say, the planes that can be determined by k linear independent points in Rd )? In general, the solutions will be non-interpolating. But we can get a sense of what will happen for β small by observing that when β = 0, d(zi , `)β becomes binary3 , zero or one, and equal to zero if and only if zi lies on the plane `. Since, generically, at most k members of ~z will lie on a k − 1 dimensional plane, the minimum cost over Pk is n − k, achieved by every interpolating plane in Pk . With this observation in mind, we introduce the notion of an “-interpolating” plane: Definition 3.2.13. Fix  > 0. A plane ` ∈ Pk is -interpolating if {zi ∈ ~z : d(zi , `) < } contains at least k linearly independent elements of ~z. From a computational viewpoint, knowing that the optimal solution is - interpolating suggests building a search algorithm based on the polynomial (versus exponential) problem of finding an optimal interpolating plane. Although we have no guarantees, the following theorem further supports this approach: Theorem 3.2.14. For any d ≥ 2, k ≤ d, and data {zi }ni=1 ⊂ Rd with span(z1:n ) = Rd , let L = {`1 , `2 , . . . , `m } be a minimizer of Function 3.2.1 over Pk . Fix  > 0 and 3 Here we set 00 = 0 119 choose   n−k+1 log n−k β<  √  log 2  n Then `j is -interpolating for every j = 1, . . . , m. Proof. According to the result from section 3.2.1, it is enough to consider the target function n d(zi , `)β X f (`) = i=1 over ` ∈ Pk . Recall that z1 , . . . , zn is standardized: kzi k2 = n. P P i=1:n zi = 0 and i=1:n Consequently, kzi k2 ≤ n, ∀i and it follows that kzi − zj k2 ≤ 4n for all 1 ≤ i, j ≤ n. Let `˜ be any interpolating plane in Pk . The distance from any zi to `˜ can not be ˜ But `˜ is interpolating and includes at least bigger than its distance to any point on `. ˜ ≤ 2√n. one data point, zj . Since kzi − zj k2 ≤ 4n, it follows that d(zi , `) Now, since `˜ contains at least k data points, n √ ˜ β ≤ (n − k)(2 n)β X d(zi , `) i=1 In contradiction to the claim in the theorem, assume that ` is not -interpolating. Then there are at least n − k + 1 data points that have distance from ` at least as great as : n d(zi , l)β ≥ (n − k + 1)β X i=1 120 Now check that   n−k+1 log n−k β<  √  log 2  n √ implies that (n − k)(2 n)β is smaller than (n − k + 1)β , which means that n n ˜β < d(zi , `)β X X d(zi , `) i=1 i=1 and ` can not be optimal. In the following paragraphs we explore two further special cases: clustering to points in d ≥ 2 dimensions and clustering to lines in d ≥ 3 dimensions. In both cases we will identify a range of values for β within which optimal solutions are guaranteed to contain at least one data point. Clustering to points in Rd , d ≥ 2. When M = Rd , d > 1, i.e. the point model, and m = 1, given data {zi }ni=1 ⊂ Rd , the target function can be written as n d(zi , x)β , X f (x) = (3.2.5) i=1 where x ∈ Rd . Note that the optimal point is not necessary to be some of the data points. We say the optimizer is non-interpolating. One can consider three points in R2 that forms an equilateral triangle, the optimal point is the center of the triangle when β close to 1. The next lemma will state the fact that every data point is a local 121 minimum of f . Lemma 3.2.15. Given data {zi }ni=1 ⊂ Rd . For i = 1, 2, · · · , n, zi is the strictly minimum point of Function 3.2.5 on the d dimensional ball center at zi with radius 1 1 ki (1 − β) 1 min{ min ||zj − zi ||, ( P 1 ) 2−β }, 2 j:zj 6=zi 2 j:zj 6=zi ||zj −zi ||2−β Pn where ki = j=1 I{zj =zi } . Proof. It is sufficient to just consider z1 is a local minimum. May assume z1 = 0. For any radius r > 0, set a point z (r) = ra, a ∈ Rd1 , such that kz (r) k = r. Then n q β f (z (r) ) = X (zi − z (r) )T (zi − z (r) ) . i=1 For abbreviation, fix an a ∈ Rd1 let g(r) , f (z (r) ) q gi (r) , (zi − z (r) )T (zi − z (r) ). We will show that when r is within the radius, we have negative Hessian of g. So for any such a point z (r) , there is a direction (toward 0) such that it produce a negative value of the directional second derivative, i.e. z (r) is not a local minimum. 122 When z (r) = ra 6= zi , calculate 1 gi0 (r) = {−2 · (zi − z (r) )T a} 2gi 1 = {(z (r) − zi )T a} gi 1 = {r − zi T a} gi since z (r)T a = raT a = r. Thus ∂ β g (r) = βgiβ−2 {r − zi T a} ∂ i whenever ra 6= z. Thus, we can calculate the Hessian of giβ 1 Hgβ (r) = β{(β − 2)giβ−3 (r − zi T a)2 + giβ−2 } i gi −1 1 = β{(2 − β) 4−β (r − zi T a)2 + 2−β }. gi gi Note that g1 (r) = r. So for those i such that zi = z1 −1 1 Hgβ (r) = β{(2 − β) 4−β (r)2 + 2−β } i r r −β(1 − β) = . r2−β And for i such that zi 6= z1 , let ri , kzi k, we have β β Hgβ (r) ≤ ≤ , i gi2−β |ri − r|2−β 123 where the second inequality is because of gi (r) = d(zi , ra) ≥ |d(zi , 0) − d(0, ra)| = |ri − r|. ri ri Furthermore, since r < 2 , for any i with zi 6= z1 , ri − r > 2 , which means 2 Hgβ (r) < β( )2−β . i ri Put these estimations together, we get 2 k1 (1 − β) ( )2−β − X Hg (r) < β{ }. i:zi 6=z1 ri r2−β 1 k1 (1−β) So, if r < 12 ( P 1 2−β ) 2−β , we have i:zi 6=z1 ri k1 (1 − β) 1 2−β k1 (1 − β) r2−β < P 2 2−β = ( ) 1 2−β 2 P i:zi 6=z1 ( ri ) j:zj 6=zi ( rj ) Thus we have Hg (r) < 0 when 1 1 k1 (1 − β) 1 2−β }. r < min{ min rj , ( P 1 ) 2 j:zj 6=z1 2 j:zj 6=zi ( rj )2−β This completes the proof. Lemma 3.2.16. Given data {zi }ni=1 ⊂ Rd , if d(zi , zj ) > q for any i 6= j s.t. zi 6= zj , then zi , i = 1, 2, · · · , n is the strictly minimum point on the d dimensional ball center 124 at zi with radius  1  2q ( 1−β  ) 2−β , for n > 2   n−1 α(n, q, β) =  . q   2 , o.w. Proof. Note that n ≤ 2 is a trivial case. Set n > 2, we can estimate the following for every 1 = 1, 2, · · · , n: ki (1 − β) 1 )−1 X 1 = ki (1 − β)( ||zj − zi ||2−β P j:zj 6=zi ||zj −z1 ||2−β j:zj 6=zi 1 )−1 X ≥ ki (1 − β)( j:zj 6=zi q 2−β n − ki −1 = ki (1 − β)( ) q 2−β ki (1 − β) 2−β = q n − ki 1 − β 2−β ≥ q n−1 for any i = 1, 2, · · · , n. Thus, we have for n ≥ 2, q 1 − β 2−β 1 1 ki (1 − β) 1 ( ) ≤ (P 1 ) 2−β . 2 n−1 2 j:zj 6=zi ||zj −zi ||2−β 1 Also, since ( 1−β n−1 ) 2−β < 1, we have q 1 − β 2−β 1 q 1 ( ) < ≤ min ||zj − zi || 2 n−1 2 2 j:zj 6=zi for any i = 1, 2, · · · , n. So 1 1 − β 2−β 1 1 1 ki (1 − β) 1 ( q) ≤ min{ min ||zj − zi ||, ( P 1 ) 2−β } 2 n−1 2 j:zj 6=zi 2 j:zj 6=zi ||zj −zi ||2−β 125 for any i = 1, 2, · · · , n. This completes the proof. One quick observation of the function α(n, q, β) is that it is increasing when β gets smaller. So that the radius will get bigger when we lower the value of β, until 1 the radius hits 2 minj:zj 6=zi kzj − zi k. Even though in Rd , d ≥ 2, the best point that minimize Function 3.2.5 might not be one of the data, we can “force” the optimal point to be one of the data by choosing a more sparse distance, i.e. smaller β. One can say that a sparser distance tends to creat an interpolating optimizer. The next theorem states this fact. Theorem 3.2.17. For data {zi }ni=1 ⊂ Rd , with q, s ∈ R s.t. q < kzi − zj k < s for all i 6= j = 1, 2, · · · , n. Given any β0 ∈ (0, 1), if β ≤ min{β0 , β1 (β0 , n, s, q)}, where  n )/ log( 2s(n−1)  log( n−1 , for n > 2   q(1−β ) 0 β1 (β0 , n, s, q) ,  , for n ≤ 2,  β0  then the global minimum points of Function 3.2.5 are among the data point. Furthermore, for n ≥ 2 if β¯ = β(n, ¯ s, q) satisfies the equation n 2s(n − 1) β¯ = log( )/ log( ¯ ), n−1 q(1 − β) then n 2s(n − 1) β¯ ≥ min{β0 , log( )/ log( )} n−1 q(1 − β0 ) 126 for every β0 ∈ (0, 1). Proof. Note that n ≤ 2 is a trivial case. Set n > 2, given any x ∈ Rd s.t. x ∈ / {zi }ni=1 and is a local minimum point of f , the Function 3.2.5. Note that for α = α(n, q, β) that defined in the previous lemma, i.e. q 1 − β 2−β 1 α(n, q, β) , ( ) . 2 n−1 So x will not be located inside any d dimensional ball Bα (zi ) for any i = 1, 2, · · · , n because of the previous lemma. Thus we have f (x) ≥ nαβ . Furthermore, given β0 ∈ (0, 1), for any β < β0 , we have 1 − β > 1 − β0 , which implies q 1 − β0 2−β 1 α(n, q, β) ≥ ( ) . 2 n−1 We can further estimate it by q 1 − β0 2−β β αβ ≥ ( )β ( ) 2 n−1 q β 1 − β0 β > ( ) ( ) 2 n−1 q(1 − β0 ) β = ( ) , 2(n − 1) β 1−β0 where the second inequality is because of 2−β < β and n−1 < 1. 127 On the other hand, we have min f (zi ) ≤ (n − 1)sβ i for any i = 1, 2, · · · , n because s is the maximum length of two data points. Thus, if n 2s(n − 1) β ≤ log( )/ log( n−1 q(1 − β0 ) n−1 q(1 − β0 ) = log( )/ log( ) n 2s(n − 1) we have q(1 − β0 ) β n−1 ( ) ≥ , 2s(n − 1) n q(1−β0 ) note that log( 2s(n−1) ) < 0. This means q(1 − β0 ) β f (x) ≥ nα ≥ n( ) ≥ (n − 1)sβ ≥ min f (zi ). 2(n − 1) i So we finish the proof of the first part. The second part of the theorem is valid since n 2s(n − 1) log( )/ log( ) n−1 q(1 − β0 ) is increasing when β0 goes to 0. 128 The next lemma explains a property of the bound. ¯ s, q) is Lemma 3.2.18. Given any β0 ∈ (0, 1), the function β1 (β0 , n, s, q) or β(n, non-increasing when n increases. Proof. The scale invariance is easy to see, since any scale imposed on the data will not affect the value qs . For the first property, we only need to look at n > 2 case. Defining η(x) = β1 (β0 , x, s, q) log x − log(x − 1) = , log(x − 1) + c 2s where c = log q(1−β 0) > 0. And calculate for x > 2 1 1 1 η 0 (x) ∝ (log(x − 1) + c)( − )− (log x − log(x − 1)) x x−1 x−1 1 c 1 = (log(x − 1) + c) − − log x x x−1 x−1 log(x − 1) + c log x + c = − x x−1 < 0 So, for n > 2, β1 (β0 , n, s, q) is decreasing when n increases. Clustering to lines in Rd , d ≥ 3. So, when β is small enough, the optimizer of the target function on the point model will be interpolating. An immediate generalization can be derived on the line model: 129 when β is small enough, the optimal line will pass though at least one data point. This will be stated in the next theorem. Theorem 3.2.19. For Rd data, d ≥ 3, with the line model. Given β0 ∈ (0, 1) and data {zi }ni=1 , the the global minimum line of the function n d(zi , `)β X f (`) = i=1 will pass at least one data point when   β1 (β0 , n, s, q) , n>3   β≤   β0  , o.w. Proof. The n ≤ 3 is a trivial case. For n > 3 and given β0 ∈ (0, 1), assume the optimal line, say `? , doesn’t pass any data point. So the target function has value n ? d(zi , `? )β . X f (` ) = i=1 May assume `? = L(0, r), i.e. it passes 0 with direction r ∈ Rd1 . We can project all zi (r)}ni=1 . the data point to the plane Pr that is perpendicular to r and passes 0, say {ˆ Then n f (`? ) = zi (r), 0)β , X d(ˆ i=1 which is equivalent to the point model in Rd−1 . Lemma 3.2.16 says that when β ≤ β1 (β0 , n, s, q), 130 a data point (note that no data point is 0), say z1 , will be the minimum point of the function n zi (r), z)β . X h(z) = d(ˆ i=1 So, we can move `? parallelly (so that it is still the normal vector of plane Pr ) such that it touches z1 . We thus find a better line L(z1 , r) that have lower target function value. So we have the first contradiction, the optimal line must pass at least one data point. 131 3.3 Observations about Sparse-distance Cluster- ing in R1 In this Section we examine the minimizer of the function n d(zi , `)β X f (`) = i=1 when ` ∈ R. The intention is to understand the simplest case so that one could have some general ideas of sparse distance. Suppose our data z = {zi }ni=1 ⊂ R is generated from some distribution F (· − ξ), where ξ ∈ R is a location parameter. To estimate the location parameter. We can define the set [Tn (z)] as n [Tn (z)] = {θ ∈ R1 : |zi − θ|β reaches minimum}, X (3.3.1) i=1 and then define T = Tn = Tn (z) as mean of the set [Tn (z)]. The definition is similar to the M-estimator from Huber [20]. The existence of the estimator for any given data set is easy to verify, since Tn must be within [min zi , max zi ], i=1:n i=1:n then we are solving a minimizing problem of a continuous function on a compact set. Thus the Tn exists. Another easy but important property is that Tn is location invariant, i.e. Tn (z + c) = Tn (z) + c for some constant c. 132 3.3.1 Asymptotic Properties However, the estimator Tn does not necessary converge to a point when β ≤ 1, unlike β > 1 case, which is easy to verify under some regular conditions on F . To 1 Pn quickly see this, note that one can obtain the limit of the function n i=1 |zi − θ|β asymptotically under some condition of strong law of large number (SLLN), this function will converge to φ(θ) = E(|Z1 − θ|β ) a.e. It is easy to verify that when β > 1, φ(θ) is a strictly convex function under some regular condition on distribution F , so Tn exist and unique in this case. For the reader that when β = 1, some choice of F may give multiple solutions, and the estimator might not be convergent. For example, think about the distribution with density 1 (I{x=−1} + I{x=1} ), 2 then Tn will take values −1, 0, and 1 infinite often and won’t converge. Thus, to have a unique asymptotic solution, one must impose some conditions on the distribution F. When β ∈ (0, 1), it’s not hard to see that [Tn ] can only contain points in {zi }ni=1 . Since the Hessian of |zi − θ|β (as a function of θ) is β(β − 1)|zi − θ|β−2 < 0, when θ 6= zi , i = 1, 2, · · · , n, which results in negative Hessian of target function for / {zi }ni=1 . any θ ∈ 133 The following lemma shows that if density of F is symmetric and unimodal, then Tn is a consistent estimator of location parameter under some mild conditions. Lemma 3.3.1. When β ∈ (0, 1), given {xi }ni=1 i.i.d. from a continuous distribution F (· − ξ). Assume the density g of distribution F satisfies 1. g is symmetric and unimodal at 0. 2. g 0 exist and continuous, 3. Eg |X| < ∞, Eg |X|β−1 < ∞ then Tn (x) converge to ξ a.e. Pn Proof. May assume ξ = 0. Since the elements of [Tn ] satisfies minθ i=1 |zi − θ|β , it satisfies n 1X min |zi − θ|β θ n i=1 also. To use SLLN, calculate Z ∞ E(|Z1 − θ|β ) = |x − θ|β g(x)dx −∞ Z θ+1 Z β = |x − θ| g(x)dx + |x − θ|β g(x)dx θ−1 R\(θ−1,θ+1) Z θ+1 Z ≤ g(x)dx + |x − θ|g(x)dx θ−1 R\(θ−1,θ+1) < ∞ because of Assumption 3. So we can use SLLN to conclude n 1X |zi − θ|β → E(|Z1 − θ|β ), n i=1 134 a.e. Thus the elements in the limit of [Tn ] are arguments of minimum of ψ(θ) , E(|Z1 − θ|β ). Observe that ψ is a even function since Z1 and −Z1 have the same distribution. Next, we will show that ψ has only one global minimum at 0, which implies Tn → 0 a.e. Note that Z ∞ ψ(θ) = |x − θ|β g(x)dx. −∞ First we have to argue that ψ 0 (θ) is existent. Since Z 0 ψ (θ) = β |x − θ|β−1 · sgn(x − θ) · g(x)dx Z θ Z ∞ β−1 = −β (θ − x) g(x)dx + β (x − θ)β−1 g(x)dx, −∞ θ R ∞ β−1 and when θ > 0, the second term is bounded from above by β 0 x g(x)dx; when R0 β−1 θ < 0, the first term is bounded from below by −β −∞ (−x) g(x)dx. And the R∞ assumption 3 states that Eg |X|β−1 exist and finite, which implies 0 (−x)β−1 g(x)dx R ∞ β−1 and 0 x g(x)dx are both finite. This indicates that ψ 0 (θ) exist for all θ ∈ R. Calculate Z ∞ ψ(θ) = |x − θ|β g(x)dx Z−∞ ∞ = |y|β g(y + θ)dy. −∞ 135 Thus Z ∞ ψ 0 (θ) = |y|β g 0 (y + θ)dy Z−∞ ∞ = |x − θ|β g 0 (x)dx. −∞ Observe that ψ 0 (0) = 0 since g 0 is an odd function (because g is an even function from assumption 1). Also, we know that ψ 0 is an odd function. Continuing calculating ψ0: Z 0 Z ∞ ψ 0 (θ) = |x − θ|β g 0 (x)dx + |x − θ|β g 0 (x)dx. −∞ 0 The first part can be estimated with the change of variable y = −x, then Z 0 Z 0 β 0 |x − θ| g (x)dx = − | − y − θ|β g 0 (−y)dy −∞ Z∞∞ = − |y + θ|β g 0 (y)dy. 0 Thus Z ∞ Z ∞ ψ 0 (θ) = − |x + θ|β g 0 (x)dx + |x − θ|β g 0 (x)dx Z ∞0 0 β β 0 = (|x − θ| − |x + θ| )g (x)dx. 0 When θ > 0, because g 0 (x) < 0 when x > 0 (Assumption 1) and |x − θ| < |x + θ|, we 136 have ψ 0 (θ) > 0. Combine the observation of ψ 0 being an odd unction, we have   > 0 , when θ > 0   ψ 0 (θ) is   <  0 , when θ < 0. Thus we have the result that ψ attend its global minimum at 0. 3.3.2 Breakdown Point Breakdown point is a way to measure the robustness of an estimator. First concept is proposed by Hampel [18]. Donoho and Huber [12] gives a simple version of it called the finite sample breakdown point, which we will adopt in this subsection. Given a sample x = {xi }ni=1 , the finite sample breakdown points at sample x for T is defined as k BP (T, x) = min{ : sup |T (x) − T (xk )| = ∞}, k n xk ∈X k (x) where X k (x) , {y ∈ Rn : yi = xi except first k elements.} for k < n. The highest breakdown point for a location estimator is 50%. We would like to show that when β ∈ (0, 1), the estimator T has 50% breakdown 137 points. Let’s first prepare some notations: n |xi − θ|β . X Q(x,θ) , i=1 Note that all the local minimizer of Q(x, ·) are data points xi , i = 1, 2, · · · , n. Let’s introduce a lemma. Lemma 3.3.2. If T (x) = c ∈ R, and x ˆ = x except one of the element being changed to c, then T (ˆ x) = c Proof. Note that since T is defined as midpoint of [T ], c may not necessary to be a data point. Assume the changed element is x1 , that is, xˆ1 = c. Define Q(·) = Q(x, ·) Q0 (·) = Q(ˆ x, ·). Note that we have Q(x) ≥ Q(c) for any x ∈ x ˆ since we assume T (x) = c. We will split the proof into two claims: T (ˆ x) = c when (1) x1 < c. (2) x1 > c. Claim. When x1 < c, T (ˆ x) = c. Proof. We will to show that Q0 (x) > Q0 (c) for any x ∈ x ˆ s.t. x 6= c, so that the value of T is not changed. We will split the problem into 3 situations: (a) x ≤ x1 , (b) x ∈ (x1 , c) (c) x > c, where x ∈ x ˆ , which covers all the the possible minimum points ˆ. in x 138 (a) When x ≤ x1 , calculate Q0 (x) = Q(x) − |x1 − x|β + |c − x|β > Q(x) ≥ Q(x) − |c − x1 |β ≥ Q(c) − |c − x1 |β = Q0 (c). Note that the first inequality is because x1 < c. (b) When x ∈ (x1 , c), calculate Q0 (x) = Q(x) − |x1 − x|β + |c − x|β ≥ Q(x) − (x − x1 )β > Q(x) − (c − x1 )β ≥ Q(c) − (c − x1 )β = Q0 (c). Note that the second inequality is because of x1 < x < c. 139 (c) When x > c, calculate Q0 (x) = Q(x) − |x1 − x|β + |c − x|β = Q(x) − [(x − x1 )β − (x − c)β ] > Q(x) − (c − x1 )β ≥ Q(c) − (c − x1 )β = Q0 (c), where the first inequality is because of aβ − bβ ≤ (a − b)β if a > b. Claim. When x1 > c, T (ˆ x) = c. Proof. Again, we split the problem into 3 cases: (a) x < c, (b) x ∈ (c, x1 ), (c) x > x1 , where x ∈ x ˆ. (a) When x < c, calculate Q0 (x) = Q(x) − |x1 − x|β + |c − x|β = Q(x) − [(x1 − x)β − (c − x)β ] > Q(x) − (x1 − c)β ≥ Q(c) − (xkk − c)β = Q0 (c). 140 (b) x ∈ (c, x1 ), calculate Q0 (x) = Q(x) − |x1 − x|β + |c − x|β > Q(x) − (x1 − x)β > Q(x) − (x1 − c)β ≥ Q(c) − (x1 − c)β = Q0 (c). (c) x ≥ x1 , calculate Q0 (x) = Q(x) − |x1 − x|β + |c − x|β > Q(x) > Q(x) − (x1 − c)β ≥ Q(c) − (x1 − c)β = Q0 (c). Theorem 3.3.3. When β ∈ (0, 1), the finite sample breakdown point at data x esti- mator T is 50%. ˆ k ∈ X k (x) with xˆki = xˆkj for Proof. Given x = {xi }ni=1 and xk ∈ X k (x), there is an x 141 i, j = 1, 2, · · · , k s.t. xk ). T (xk ) = T (ˆ This is because of Lemma 3.3.2, i.e. we can change xki , i = 1, 2, · · · , k, to T (xk ) one point at a time and doesn’t affect the estimator T . Define Xˆ k (x) = ∪c∈R x ˆ k,c ˆ k,c ∈ X k (x) s.t. xˆki = xˆkj = c for i, j = 1, 2, · · · , k. So and x |T (x) − T (xk )| ≤ sup xk )| |T (x) − T (ˆ ˆ ∈Xˆ k (x) x k xk,c )|. = sup |T (x) − T (ˆ c∈R ˆ k,c ∈ X k (x), we have Since x sup xk,c )|, |T (x) − T (xk )| = sup |T (x) − T (ˆ xk ∈X k (x) c∈R ˆ k,c ∈ X k (x) and xˆki = c for i = 1, 2, · · · , k. Set k < n/2 and define where x xk,c , ·). Q(·) = Q(ˆ We may assume that for all i = 1, · · · , n, xi ∈ (−, ) for some  > 0, and T (x) = 0. Calculate n |c − xi |β > (n − k)|c − |β . X Q(c) = i=k+1 142 And n Q(0) = k · |c|β + |xi |β X i=k+1 < k · |c| + (n − k)β . β Since n − k > k, there exist a c0 = c0 (x, k) > 0 s.t. Q(0) < Q(c) whenever |c| > c0 . xk,c )| So the minimum point of Q will not be c when |c| > c0 , which means |T (x) − T (ˆ xk,c )| < c0 +  < ∞, can only be among (−,). And when |c| ≤ c0 we have |T (x) − T (ˆ This implies sup |T (x) − T (xk )| < ∞ xk ∈X k (x) when k < n/2. So k BP (T, x) = min{ : sup |T (x) − T (xk )| = ∞} = 50%. k n xk 143 3.4 Minimizing the Loss Function Given data points ~z = {zi }ni=1 ⊂ Rd , a model M, a penalty function γ : 2M → R, and an error function ρ : Rd × M, the target function we try to minimize is n λX f (L) , γ(L) + min ρ(zi , `) (3.4.1) n i=1 `∈L for L ∈ 2M , or say L ⊂ M. When |M| < ∞, solving arg minL⊂M f (L) with brute force involves O(2|M| ) calculations. Here we propose Remove-or-Replace algorithm (RRA) on approximating a solution of arg minL⊂M f (L) when |M| < ∞. 3.4.1 Remove-or-Replace Algorithm The procedure of remove-or-replace algorithm, RRA, is as follows. 1. (Initialization) Choose an initial set L ⊂ M, say L = {`1 , · · · , `|L| }. 0 2. For i = 1, · · · , |L|, set L(i) = L\`i , calculate `i with 0 `i = arg min f (` ∪ L(i) ). `∈φ∪M\L Note that we also include φ, called the removal element, such that φ∪L(i) = L(i) . That is, if the minimizer is φ, we know that the target function is smaller without any choice of `i (among M\L). If there are multiple arguments of minimum, pick φ if φ is among them, otherwise, pick the first occurrence of ` ∈ M\L. 144 Record the value of target function 0 ci = f (`i ∪ L(i) ), 0 and the corresponding set Li = `i ∪ L(i) . |L| |L| 3. After previous step we have {ci }i=1 and {Li }i=1 . Calculate r = arg min ci , i=1:|L| and replace L with Lr if f (L) < f (Lr ). Again, if there are multiple arguments of minimum, pick the first occurrence of i in c1 , · · · , c|L| . 4. Repeat step 2 and 3 until there is no further decreasing on target function. The algorithm found L as an approximate minimizer of f . The complexity of the algorithm is O(q 2 |M|), where q is the number of element on the initial set. It is clear that every step of RRA decreases the target function f , which make it a greed algorithm. The reason we don’t want to specify the number of instance, m, instead of introducing parameter λ in our loss function (Function 3.4.1) is the ability to perform RRA. Note that RRA is applicable for any function f (L) where L ⊂ M with |M| < ∞. In the rest of this subsection, we will go through the setting of experiments on RRA. The error function we choose is denoted as the modified sparse distance, which depends on the data and the model chosen: d(·, ·)β ρ(·, ·) = , ζ~z,M 145 where β ∈ (0, 1) and |M| n 1 d(zi , `j )β XX ζ~z,M = n · |M| j=1 i=1 is the modified parameter ζ = ζ~z,M . The main reason we include it is that we would like to have similar loss function value on the sum of error term between different models. For example, given the same set of data ~z, the point model, will have bigger value of ζ than the line model. This small modification make parameter λ less sensitive to model changes. Note that since we assume β ∈ (0, 1) throughout this article, Corollary 3.2.12 greatly reduces the problem on minimize Equation 3.4.1 under certain situations. When the line model L is considered, the corollary tells us that we only need to con- sider C2n interpolating lines, i.e. C2n < ∞, which enables us to use RRA algorithm. More generally, when using the plane model Pd on Rd , we only need to consider inter- polating Cdn planes. We also found that even when the optimizer is not interpolating, if we restrict the model to those instances that interpolate, the result is still very informative. Thus, given the data ~z = {zi }ni=1 ⊂ Rd , the models we choose will be restricted to interpolating models, Mp , where p ∈ N, p ≤ d, means the model is the collection of instance that passing p linearly independent data points z1 , z2 , · · · zp . For example, M1 is the interpolating point model and M2 is the interpolating line model, etc. Also, we will choose a penalty function γ(L) = |L|2 . 146 Thus, for the most of the time the target function we use is n λ f (L) , |L|2 + min d(zi , `)β X (3.4.2) n · ζ~z,M i=1 `∈L for L ⊂ Mp , p ≤ d. Note that generally speaking we have ζ~z,Mp2 < ζ~z,Mp1 if d ≥ p2 > p1 ∈ N. The error function is modified so that the parameter λ is less susceptible to p. To evaluate how good the approximation of RRA is, we will also calculate the exact solution of the minimizer of f (L) for small amount of instances. Since the penalty function only depends on the number of instances, we can compare RRA and the exact solution by the sum of error, h. That is, for each m, the number of instances, we examine the instances acquired by RRA or by the exact solution with the sum of error n X h(L) , min ρ(zi , `), (3.4.3) `∈L i=1 with |L| = m and L ⊂ Mp . The h value will be shown on figures that we provided in the next few subsections. An easy observation is that when m = 1, the exact solution will be the same as RRA gives, so we will omit displaying it. Note that when m is N known, to look for exact solution will cause Cm , where N = Cpn , computations, and we usually cannot find the exact solution for larger p and m. About the initial set, since we can easily choose n instances of model that make data term vanished, there is no point to have more then n elements in the initial set L. (Actually, for the line model case, we don’t need more than [ n2 ] elements, because any possible line must connect two data points.) When point model M1 is considered, the initial set is the same as M1 , i.e. the data points. For other models, there are 147 two strategies to choose the initial set: 1. (Locally best initial) Given data {zi }ni=1 , for each data point zi , i = 1, · · · , n, we can choose the best instance, say `i , of the given model that fit lo- ρ(zj , `i ), where Nk (zi ) , {zj ∈ P cally such that `i that minimizes zj ∈Nk (zi ) z|zj is zi0 s k nearest neighbor.}. So we have n instances. Then, remove the redundant instances, and we get the initial set L = {`i }qi=1 to work on the aforementioned algorithm. The nearest neighbor parameter k is [ n5 ] through the experiments if this initialization method is chosen. 2. (Random initial) Randomly choose q instances, where q ≤ n, whichever make data term vanished. 3.4.2 Robust Nature for Sparse Distance We mentioned that we deliberately choose β ∈ (0, 1) due to the robustness issue. For multiple line instances that appeared in the data, our following experiment shows the robust nature is also retained on minimizing the function n λ f (L) , |L|2 + min d(zi , `)β , X (3.4.4) n · ζ~z,M2 i=1 `∈L for L ⊂ M2 in R2 . Figure 3.4.1 plots both the result of EM algorithm (with number of instances m = 2 is known) and RRA result (with the line model and random initial) to the corrupted X-shaped data. The X-shaped data has originally 60 data point, and is corrupted by adding 40 outliers. The EM algorithm assumes two line instances with Gaussian error, and is initialized with the exact 2 lines that generate 148 the 60 un-corrupted data. We can see that the result of EM algorithm gives a poor result in this example, compare to that from RRA. In this case, the exact solution of minL∈M2 f (L) for m = 2 is the same as we got from RRA, which indicates that the loss function f is indeed a great choice for this corrupted data. (a) (b) (c) Figure 3.4.1: We examine the minimizer of Function 3.4.4 for L ⊂ M2 , a set of interpolating lines in R2 on X-shaped data with outliers. (a) The original data consists of 100 data points that corrupted by 40 outliers. (b) EM algorithm for 2 instances. (c) RRA result for β = 0.8 and λ = 30 with random initial. In this case, the exact solution coincides RRA result. Another things to point out is about the choice of λ. For data with this many outliers, one can decrease the value of λ to reduce the number of lines that RRA will find, so we set λ = 30 in this case. Actually, when one choose λ between 10 to 50, RRA will give two lines, see Figure 3.4.2. Note even when λ is miss-specified, the instances that RRA gives is still very informative. 149 Figure 3.4.2: Different λ setting to RRA when minimizing Function 3.4.4 for L ⊂ M2 , a set of interpolating lines in R2 , on X-shaped data with outliers. One can see that for λ = 10 ∼ 50, the algorithm gives two lines. Even when other λ is chosen, the instances that RRA found is still informative. Note that here we use random initial for RRA. 150 3.4.3 Instances from a Single Model The intention of this subsection is to see how RRA work on various of data from a single model. For reader’s reference, the target function we used here is n λ f (L) , |L|2 + min d(zi , `)β , X (3.4.5) n · ζ~z,Mp i=1 `∈L where L ⊂ Mp . (Mp is the model that collects all the interpolating Rp−1 linear structure, i.e. the linear structure that determined by p linear independence data points.) We already see the loss function produces a robust result in Subsection 3.4.2. In particular, in this subsection we would like to know: 1. How general can we choose the parameter λ? Also, does initialization of RRA matters a lot? 2. How is the result from RRA compare to the exact solution? 3. How much do we sacrifice if we only look for interpolating solution? For exam- ple, for data in R2 , if we look for instances from the interpolating point model M1 , i.e. the data point. Another example is for data in R3 , while we look for instances from M2 , the collection of the lines that connect two data points. 3.4.3.1 Line Model M2 on R2 Data For line model on R2 data, to minimize Function 3.4.5 on all the lines in R2 is equivalent to minimize on M2 , since the minimizer is interpolating. The data that 151 we will examine line model on R2 is shown in Figure 3.4.3, which is called the "S-curve" data generated from from (xi , yi ) = (ui , −0.02u2i − 0.005ui + ui + 10 + i ), where ui ∼ U (−10, 10) and i ∼ N (0, 1). The exact solution for the loss function we proposed (i.e. Function 3.4.5) is shown in Figure 3.4.4 for m = 2 and 3. The sum of error, h, as we mentioned in Equation 3.4.3 is shown on the top of the plots which will help examine the result from RRA. Figure 3.4.3: The "S-curve" data consists of 50 data points (n = 50). The result of RRA with locally best initial (described in Subsection 3.4.1) is shown in Figure 3.4.5. We can see that the algorithm obtains the exact solution in this case, since the sum of error is the same as in the exact solution. Observe that the parameter λ controls the weight on fitting the data and the number of lines it could find. Most of the time the algorithm found a reasonable solutions. We also tried RRA with random initial (described in Subsection 3.4.1) on the S-shaped data. The result is plot in Figure 3.4.6 where the different random seed is used. One of six initials didn’t find 3 desired line instances, which gives an idea that 152 (a) Exact Solution for m = 2 (b) Exact Solution for m = 3 Figure 3.4.4: The exact solution (plotted as set of lines) of minimizing the loss function (Equation 3.4.5) while we use β = 0.8. (a) When the amount of instances m = 2 and (b) m = 3. Figure 3.4.5: The result of RRA with locally best initial on S-curve data (n = 50) for different λ, where β = 0.8. When λ = 20 ∼ 110 (the second and the third plots), RRA gives 3 line instances. Note that the sum of error, h, for m = 2 is 18.711 (the first plots) and for m = 3 is 9.094 (the second and the third plot) coincide with that of exact solution, means that RRA approximates the exact solution well in this data. 153 RRA is not very sensitive to initials. Generally speaking, random initial still gives a very acceptable result in this example Figure 3.4.6: The result of RRA (β = 0.8, λ = 110) with random initial on S-shaped data with different random seed. From left to right, the random seed are set as 0,1,2,3,4,5, respectively. One can see that initialization does not cause huge different on RRA result in this sample. Some other example is shown in Figure 3.4.7. 154 Original data Exact Solution RRA Figure 3.4.7: Some results of various data in R2 . RRA is initialized with locally best initial. 155 3.4.3.2 Line Model M2 on R3 Data Here we examine RRA on the interpolating line model on R3 . Note that the optimizer will, in general, not interpolate. Given r = (r1 , r2 , r3 ) ∈ R31 with r1 > 0, σ > 0 and a point z0 = (x1 , x2 , x3 ) ∈ R3 , a set of data {zi }ni=1 = {(xi1 , xi2 , xi3 )}ni=1 ⊂ R3 forms roughly a line in 3D where xm ≤ xi1 ≤ xM can be generated by θ = (r, σ, z0 , xm , xM ) xm − x0 xM − x0 ti ∼iid U ( , ) r1 r1 i ∼iid N (0, σ 2 ) x1i ∼iid x1 + ti · r1 x2i ∼iid x2 + ti · r2 + i x3i ∼iid x3 + ti · r3 + i , for i = 1, 2, · · · , n. The test data is generate by three different sets of θ: θ1 = (0, 0, 3, 0.3, 5, −4, 1, −3, 5) θ2 = (1, 2, 5, 0.3, −2, 3, 4, −1, 4) θ3 = (1, 1, 1, 0.3, 1, 10, 3, 0, 4). The original data is called “Line 3D” data, plotted in Figure 3.4.8. Figure 3.4.10 compares RRA result with exact solution. RRA is performed with locally best initial and β = 0.8, λ = 50. One can see that RRA gives a good approximation on the exact solution. We also test it with random initial with various random seeds, the result is very consistent and similar to the result from locally best initial. 156 (a) (b) Figure 3.4.8: “Line 3D” data for the experiment on the interpolating line model on R3 . (a) and (b) are different viewing angle. (a) Exact solution m = 2 (b) Exact solution m = 3 Figure 3.4.9: The exact solution for “Line 3D” data. (a) The exact solution for m = 2 with different viewing angles on top and bottom. (b) The exact solution for m = 3. 157 Figure 3.4.10: Result of “Line 3D” on R3 data. Top row and bottom row are different viewing angle of the result. From left to right is the result of different λ setting. A desired results are occurred when λ ∈ (30, 250). Note that RRA coincides with the exact solution. 158 3.4.3.3 Plane Model M3 on R3 Data In this subsection we examine RRA when line model is considered in R3 . The data is generated from: x = (x1 , x2 , x3 ) ∈ R3 K 1 X x ∼ g(x|θk ), K k=1 where gk (·) is a density function defined as: 1 (ak x1 + bk x2 + dk ) 2 g(x|θk ) ∝ I[m1k ,M1k ] (x1 ) · I[m2k ,M2k ] (x1 ) · exp(− 2 (x3 − )) 2 · σk ck θk = (m1k , M1k , m2k , M2k , ak , bk , ck , dk , σk ) The planes data has 40 data points generated from 2 set of parameters (K = 2): θ1 = (0, 5, −5, 5, 0, 0, 1, 0, 0.2) θ2 = (0, 5, −5, 5, 0, 1, 1, 0, 0.2), which is call “Plane 3D” data set. The original data is plotted in Figure 3.4.11 and the exact solution is shown in Figure 3.4.12 which uses β = 0.8, λ = 50. To represent the result, we plot it with two different perspective. RRA result is shown in Figure 3.4.13 (initial with the locally best initial). We can see that the result is identical to the exact solution. 159 (a) (b) Figure 3.4.11: The“Plane 3D” data set for the experiment on the plane model on R3 . (a) and (b) are different viewing angle. (a) Exact solution (b) Anther viewing angle Figure 3.4.12: The exact solution for m = 2 for “Plane 3D” data. (a) and (b) are different viewing angle. 160 Figure 3.4.13: RRA results of “Plane 3D” on R3 data. Top row and bottom row are different viewing angle of the results. From left to right are the results of different λ setting. A desired results are occurred when λ ∈ (10, 50). Note RRA coincides with the exact solution in this case. 161 3.4.4 Instances from Multiple Models So far we have examined RRA for multiple instances from a single model. We would like to consider the case when multiple instances of different models appears in data in this subsection. In general, we want to find structures from a set of data ~z = {zi }ni=1 ⊂ Rd that have instances from multiple (linear) interpolating models M1 , M2 , · · · , Md . Note that, for p ≤ d, p ∈ N, Mp is the model that contains linear structures in Rd that passing through p linear independent data points. One can still use RRA by setting M = ∪dp=1 Mp , and look for an interpolating minimizer. However, if there is no penalty on using more sophisticated model, the sum of the error function that Pn we preferred, i=1 min`∈L d(zi , `)β , will favor instances from Md , the highest level of models. Due to this problem, we provide two strategies 1. Penalizing higher level models, and 2. Inspect with the sum of minimum volume, as the following two subsections. 3.4.4.1 Penalizing Higher Level Models The first strategy is to penalize the usage of the more complicated model by adding heavier weights on the error function. That is, one considers the following target function n λX ωψ(`) · d(zi , `)β f (L) , |L|2 + min , n i=1 `∈L ζ~z,Mψ(`) 162 for L ⊂ M = ∪dp=1 Mp , and ψ : M → {1, 2, · · · , d} s.t. ψ(`) = p for ` ∈ Mp , p ∈ {1, 2, · · · , d}, and ω : {1, 2, · · · d} → R is the model penalty parameter s.t. ωp1 ≤ ωp2 for p1 < p2 ∈ {1, 2, · · · , d}. Note that the modified parameter ζ~z,Mp also penalize the more advanced model since ζ~z,Mp1 > ζ~z,Mp2 for p1 < p2 . However, in practice, one often need to add model penalty term ωp to acquire a satisfactory result. As our experiments shown, see Figure 3.4.15 and Figure 3.4.16, the results are sensitive to the choice of the model penalty term ωp . A caveat that worth noting is that we use the error function as ω(ψ(`)) · d(zi , `)β ρ(z, `) = , ζ~z,Mψ(`) which helps on finding instances, but does not possess the trait of Euclidean distance, i.e. given a data point z and two linear instances `1 ∈ Mp , `2 ∈ Mq , p 6= q ∈ {1, 2, · · · , d}, ρ(z, `1 ) < ρ(z, `2 ), does not necessary imply d(z, `1 ) < d(z, `2 ). This might cause some problems when we try to cluster with the error function ρ(·, ·). Due to this fact, when we do clustering, we use the usual Euclidean distance d(·, ·). We generate two sets of 2D data, “percent” and “Corrupted Percent” to examine the effect of the model penalty parameter ω, which is shown in Figure 3.4.14. One of the data set is called clean percent shape data (see Figure 3.4.14(a)) that consists of 3 instances, two point instances and a line. For another data set (see Figure 3.4.14(b)) we add 40 uniformly distributed noises onto the clean percent shape data, which help us gain the idea of how robust the method is. The result for the clean percent shape data is displayed in Figure 3.4.15. One can 163 see that when ω1 = 1, and ω2 = 2 ∼ 7, RRA successfully find the 3 desired structure. And when [ω1 , ω2 ] = [1, 1.5], RRA favors line structure, a higher level model, as we mentioned. On the contrary, when [ω1 , ω2 ] = [1, 8], RRA favors point instances, which indicates over compensation happened. The result for the corrupted percenti shape data is in Figure 3.4.16. When ω1 = 1, ω2 = 1.5 ∼ 6, we see that RRA find the structures, even under the presence of 30.8% noises. (a) (b) Figure 3.4.14: The 2D data “percent” and “Corrupted Percent” for experiment in this section. (a) “Percent” data consists of n = 90 data points. (b) “Corrupted Percent” data has 40 uniformly distributed noises as outliers in “percent” data (30.8% of the data). 3.4.4.2 Inspecting with Minimum Sum of Volumes Another strategy is less sensitive to the choice of parameters. We called it the mini- mum sum of volumes (MSV ) method. The idea is to utilize the result instances from RRA with different models as candidates, and pick the candidates by minimize the 164 (a) (b) (c) (d) Figure 3.4.15: Results of RRA with the model penalty parameter [ω1 , ω2 ] on the “Percent” data. (n = 90). While setting ω1 = 1, a satisfactory result is appeared when ω2 = 2 ∼ 7 (Plots (b) and (c)). The parameters are denoted on top of each plots. One can see the penalty term acts on RRA result. The big dot markers are the instances in M1 found by RRA; the line is the instance in M2 found by RRA. (a) (b) (c) (d) Figure 3.4.16: Results of RRA with the model penalty parameter [ω1 , ω2 ] on “Cor- rupted Percent” data (n = 130). While setting ω1 = 1, a satisfactory result is appeared when ω2 = 1.5 ∼ 6 (Plots (b) and (c)). The parameters are denoted on top of each plots. One can see the penalty term acts on RRA result. The big dot markers are the instances in M1 found by RRA; the line is the instance in M2 found by RRA. 165 volume of the structures. That is, the candidates are picked from RRA with Mp for p = 1, 2, · · · , P , by minimizing n λX f (L) , |L|2 + min ρ(zi , `) n i=1 `∈L for L ⊂ Mp , where d(zi , `)β ρ(zi , `) = . ζ~z,Mp Given an instance `p ∈ Mp , a cluster of data A ⊂ ~z, and u ∈ (0, 100], the volume is defined as vu (`p , A) , C(rinclass u , p − 1) · C(rdist , d − (p − 1)), where π k/2 C(r, k) , rk Γ( k2 + 1) is the volume of a k-dimensional ball, and u u rdist = rdist (`p , A) , Qu ({d(z, `p )}z∈A ), where Qu (C) is the uth-percentile of the 1-dimensional data C, u ∈ (0, 100]. rinclass = rinclass (`p , A) , max d(x, x¯), x∈P roj(A,`p ) where P roj(A, `p ) = {P roj(z, `p )}z∈A is the set of points that is projected to `p from data in A. Note that rinclass (`, A) = 0 for ` ∈ M1 , a instance from the point model, and C(0, 0) = 1. For example, when the dimension of data d = 3 is considered, the volume of a point instances along with the data it clustered is a ball in R3 with radius 166 rdist ; the volume of a line instances along with the data it clustered is a volume of a cylinder in R3 , which has a height rinclass and a circle bottom with radius rdist ; the volume of a plane instances along with the data it clustered is another volume of a (flatter) cylinder in R3 , which has a height rdist and a circle bottom with radius rinclass . A example of volume in R3 is shown in Figure 3.4.21. Given u ∈ (0, 100], the MSV method is a procedure for finding m instances that have minimum sum of volumes, M SVu (m), is: 1. For each p ∈ {1, 2, · · · , P }, obtain candidate of instances Lp from model Mp by minimizing n λ f (L) , |L|2 + min d(zi , `)β X n · ζ~z,Mp i=1 `∈L for L ⊂ Mp . (This is done by RRA.) In the end we have Lcand = ∪Pp=1 Lp is the collection of candidate instances from different models. Note that one should expect k > m, otherwise, raise the value of λ so one can have more candidates found. 2. Set k = |Lcand |. Thus we have N = Cm k combinations of m instances out of k candidates. The set of every combination is represented as {M1 , M2 , · · · , MN }. 3. For any combination , say M = {`1 , `2 , · · · , `m }, obtain the cluster Aj = {z ∈ ~z|d(z, `j ) < d(z, `), ∀` ∈ L} for j = 1, 2, · · · , m. If two or more instances minimize the same data, randomly choose one. Calculate the sum of volumes, Pm SVu (M ) = j=1 vu (`j , Aj ), for this combination of candidates. If there is a cluster Aj , j = 1, 2, · · · , m, which has less than  amount of data, i.e. |Aj | < , set SVu (M ) = ∞. 167 ˆ , that has minimum sum of volumes. i.e. M 4. Choose the combination, M ˆ = arg mini=1:N SVu (Mi ) Note that we use d(·, ·), the Euclidean distance, to cluster and calculate volumes, instead of using the sparse distance d(·, ·)β . This is because Euclidean distance is more natural to use when clustering (compare to ρ) and to calculate volumes, where both of them relies on geometrical properties of the data. For a noisy data, we use the parameter u to keep the method robust. Also, in step 3 when calculating M SVu (m), we exclude the combination that contains a cluster that less than  amount of data. This is to prevent some complications that might occur, especially when some outliers are present in the data that makes the computation of volumes difficult. We set n  = [ 15 ]. We apply MSV method on a 2D data .The original data and candidates are shown in Figure 3.4.17, and result of MSV is shown is Figure 3.4.18. In this data set, 30.8% of the data are outliers. The application of MSV method on a 3D data is shown with data “Corrupted 3D mixture” in Figure 3.4.19 and Figure 3.4.20, where we have 15 outliers out of 65 data points. The uncorrupted version of the data consists of a point, a line, and a plane instance (total 3 instances). Also, in Figure 3.4.21 we show the volume derived for each instance in m = 3 case. 168 (a) (b) Figure 3.4.17: (a) “Corrupted percent” data, where 40 out of 130 data points are out- liers (30.8% of the data). (b) Candidates for MSV method on percent data, required from RRA at step 1 of MSV method. There are total 6 candidates, represented as bigger dots and lines, appeared in this data (a) (b) (c) Figure 3.4.18: Results of MSV method on “Corrupted percent” data. Assuming u = 70, that is, roughly speaking, we use 70% of the data in a cluster to calculate the volume. (a) Assume the amount of instances m = 2. (b) m = 3. (c) m = 4. Note that the value of M SV (m) is denoted at the top of each plot. When m = 4, there is no combination that provides the sufficient amount (> ) of data in a cluster, so M M SV (4) = ∞ (on the top of the plot, the MSV value is 3053.3, which is just a big value that indicate this happened). The transperant line is indicated as the “volume” of each instances. 169 (a) (b) Figure 3.4.19: The “Corrupted 3D mixture” data, where 15 out of 65 data points are outliers (23.1% of the data). (a) The original data from different viewing angle (top and bottom). (b) Candidates for MSV method on percentile data, required from RRA, from different viewing angle (top and bottom). There are total 9 candidates, represented as bigger dots, lines, and planes, appeared in this data. 170 (a) (b) (c) Figure 3.4.20: Results of MSV method on corrupted 3D mixture data. We provide different viewing angles by separating them at top row and bottom row. Assuming u = 70, i.e. 70% of data in a cluster is used to compute the volume, roughly speaking. (a) Assume the amount of instances m = 2. (b) m = 3. (c) m = 4. Note that the value of M SV (m) is denoted at the top of each plot. Figure 3.4.21: m = 3 result for “corrupted 3D mixture“ data. The transparent area is the volume we calculated for each instance. 171 3.5 Future Work There are many variations on the theme and as of yet few analytic guidelines. Future investigations will include: 1. Performance Bounds. We seek exact and tight bounds on the performance of the RRA algorithm. Our experiments, as well as the −interpolation theorem of section 3.2.4 suggest that RRA, or perhaps another coarse-to-fine type search strategy, finds an excellent approximate solution. Can we bound the difference in loss between the RRA and the global optimum? The problem is reminiscent of the covering problem, where such bounds do exist. 2. The Role of |L|2 . We have suggested various ways of choosing the number of instances. The performance of any algorithm is likely to be sensitive to the complexity term (complexity “penalty”) in the loss function. We have experi- mented exclusively with the quadratic loss. How will performance differ using other loss functions, e.g. other powers of |L|? More theory is needed. 3. −interpolation: can  be set to 0? Perhaps the most glaring hole in our analytic treatment concerns the theory on −interpolation: In section 3.2.4 we constructed a function βo (n, k) such that β < βo guarantees that all instances of an optimal solution would be −interpolating. Does there exist a similar result for  = 0? In other words, does there exist a function βo (n, k) (or, perhaps, βo (n, k, d)) such that β < βo guarantees that all instances of an optimal solution are interpolating? 4. Extending the library of shapes. We can generalize to a far larger library 172 of surfaces, e.g. the quadratic surfaces in Rd , perhaps defined by a system of implicit equations, M ↔ F (z, θ) = 0, z ∈ Rd , one equation for each θ. If M is closed under Euclidean transformations, and perhaps to scale as well, then many of the same analytic and computational methods developed here for reducing search to the set of interpolating instances are still appropriate, albeit with greater theoretical and computational challenges. We plan to work towards the attractive idea of a having a rich library of manifolds that can be used to create a summary of a complex high-dimensional data set through a collection of multiple, heterogeneous instances. Chapter Four Conclusions 174 4.1 On Entropy Rate of Images We examine how modern lossless image compression algorithms performed. Specifi- cally, we aim at the dominated method in this field, the predictive context modeling method. To this end, the extended Shannon-McMillian-Breimen theorem is proved. First, under stationary condition, the coding rate will converge to the entropy rate of its own ergodic component for almost every image. That is, given an image, the entropy rate of the ergodic component with respect to the image is the theoretical minimum bits per pixel (bpp) of any compression algorithm. Second, with light as- sumptions, the predictive context model with a stationary probability will converge to the theoretical minimum when the size of context grows. In order to find out the minimum bpp for an image, we proposed CBIC, the comparison-based image com- pression algorithm. The compression scheme is a type of “the more you know the better you get” algorithm, which means we can put as many as resources (in our case, the library) we want to get a better result. This is done by utilizing the idea of empir- ical distribution from an external source. We see the bpp result converged when size of context grows, like the theory suggests. In the experiment, we get slightly better (about 0.05 to 0.15 bpp) result compare to one of the best compression scheme we have so far. This indicates that these schemes may already be close to the optimal. 4.2 On Robust Generalized Clustering We propose a new loss function for clustering data points to shapes. Thus far we have focused on “linear shapes,” meaning points, lines, and in general, high-dimensional 175 planes, but there are natural extensions to parametric classes of surfaces, such as quadratic and cubic surfaces. The idea is to fit multiple surfaces to the data (e.g. a heterogeneous collection of planes, including points, lines, and planes of different dimensions), assigning each point to the closest surface. Importantly, the loss function is a concave function of the distance from data points to their nearest surface. This induces “sparse” distances, meaning that optimal surfaces tend to be near or even on top of data points—sparse in the sense of many zero distances. Furthermore, concave loss is, by design, robust to outliers. In the case of planes, depending on the dimension of the data and the dimension of the planes, we can sometimes prove that every member of the minimizing collection is “interpolating,” i.e. completely determined by a set of data points at zero distance. Related results indicate that the optimal planes will be close to or containing data points, provided that the concavity of the loss respects an explicit lower bound. Such results are useful for computation: whereas the minimization of the loss function is in general NP hard, a polynomial complexity “Remove or Replace” algorithm gives excellent approximations in cases where the actual minimum can be (laboriously) computed. 176 Bibliography [1] Nabih N. Abdelmalek. On the discrete linear L1 approximation and L1 solutions of overdetermined linear equations. J. Approximation Theory, 11:38–53, 1974. ISSN 0021-9045. [2] E. R. Barnes. An algorithm for separating patterns by ellipsoids. IBM J. Res. Develop., 26(6):759–764, 1982. ISSN 0018-8646. doi: 10.1147/rd.266.0759. URL http://dx.doi.org/10.1147/rd.266.0759. [3] P. Billingsley. Ergodic Theory and Information. Wiley, New York, 1965. [4] Paul T. Boggs and Janet E. Rogers. Orthogonal distance regression, 1990. URL http://dx.doi.org/10.1090/conm/112/1087109. [5] R. J. Boscovich. De literraria expeditione per pontificiam synopsis amplioris operis, ac habentur plura eius ex exempla ria etiam sensorum impressa. Bonon- iensi Scientiarum et Artium Inst. Atque Acad. Commen., 4:353–396., 1757. [6] L. Breiman. The individual ergodic theorem of information theory. Ann. of Math. Stat, 28:809–811, 1957. [7] Leo Breiman. Hinging hyperplanes for regression, classification, and function 177 approximation. IEEE Trans. Inform. Theory, 39(3):999–1013, 1993. ISSN 0018- 9448. doi: 10.1109/18.256506. URL http://dx.doi.org/10.1109/18.256506. [8] Leo Breiman, Jerome H. Friedman, Richard A. Olshen, and Charles J. Stone. Classification and regression trees. Wadsworth Statistics/Probability Series. Wadsworth Advanced Books and Software, Belmont, CA, 1984. ISBN 0-534- 98053-8; 0-534-98054-6. [9] Thomas M. Cover and Joy A. Thomas. Elements of information theory. Wiley- Interscience [John Wiley & Sons], Hoboken, NJ, second edition, 2006. ISBN 978-0-471-24195-9; 0-471-24195-4. [10] A. P. Dempster, N. M. Laird, and D. B. Rubin. Maximum likelihood from incomplete data via the EM algorithm. J. Roy. Statist. Soc. Ser. B, 39(1): 1–38, 1977. ISSN 0035-9246. URL http://links.jstor.org/sici?sici= 0035-9246(1977)39:1<1:MLFIDV>2.0.CO;2-Z&origin=MSN. With discussion. [11] Terry E. Dielman. Least absolute value rregression recent contributions. Journal of statistical computation and simulation, 75(4):263–286, apr . [12] David Donoho and Peter J. Huber. The notion of breakdown point. In A Festschrift for Erich L. Lehmann, Wadsworth Statist./Probab. Ser., pages 157– 184. Wadsworth, Belmont, CA, 1983. [13] R. Durrett. Probability: Theory and Examples. Duxbury advanced series. Thom- son Brooks/Cole, 2005. ISBN 9780534424411. URL https://books.google. com/books?id=NPYYAQAAIAAJ. [14] Jerome H. Friedman. Multivariate adaptive regression splines. Ann. Statist., 19 (1):1–141, 1991. ISSN 0090-5364. doi: 10.1214/aos/1176347963. URL http: 178 //dx.doi.org/10.1214/aos/1176347963. With discussion and a rejoinder by the author. [15] F. Golchin and K. K. Paliwal. A lossless image coder with context classifica- tion, adaptive prediction and adaptive entropy coding. Proceedings of IEEE International Conference on Acoustics, Speech and Signal Processing, Seattle, Washington, USA, pages 2545– 2548., May 1998. [16] Robert M. Gray. Probability, Random Processes, and Ergodic Properties. Springer, 2010. [17] Robert M. Gray. Entropy and Information Theory. Springer, 2013. [18] Frank R. Hampel. A general qualitative definition of robustness. Ann. Math. Statist., 42:1887–1896, 1971. ISSN 0003-4851. doi: 10.1214/aoms/1177693054. URL http://dx.doi.org/10.1214/aoms/1177693054. [19] Sadri Hassani. Mathematical physics. Springer, Cham, second edition, 2013. ISBN 978-3-319-01194-3; 978-3-319-01195-0. doi: 10.1007/978-3-319-01195-0. URL http://dx.doi.org/10.1007/978-3-319-01195-0. A modern introduc- tion to its foundations. [20] Peter J. Huber. Robust estimation of a location parameter. Annals of Mathe- matical Statistics, 35(1):73–101, 1964. [21] K. Jacobs. The ergodic decomposition of the Kolmogorov-Sinai invariant. Aca- demic Press, New York, 1963. [22] Pierre-Simon Laplace. Théorie analytique des probabilités. Mme Courcies, Paris, 1812. 179 [23] Kenneth D. Lawrence and Jeffrey L. Arthur. Robust Regression: Analysis and Applications. Marcel Dekker, 1990. [24] E. Lindenstrauss. Pointwise theorems for amenable groups. Invent. math, 146, 2001. [25] G. Seroussi M. J. Weinberger and G. Sapiro. LOCO-I: A low complexity, context- based, lossless image compression algorithm. Data Compression Conference (Snowbird, UT), pages 140–149, 1996. [26] J. J. Rissanen M. J. Weinberger and R. B. Arps. Applications of universal context modeling to lossless compression of gray-scale images. IEEE Trans. on Image Processing, 5(4), April 1996. [27] J. MacQueen. Some methods for classification and analysis of multivariate ob- servations. pages Vol. I: Statistics, pp. 281–297, 1967. [28] Naresh Manwani and P. S. Sastry. K-plane regression. CoRR, abs/1211.1513, 2012. URL http://arxiv.org/abs/1211.1513. [29] Bernd Meyer and Peter Tischer. TMW - a new method for lossless image com- pression. pages 533–538, 1997. [30] N. M. Nasrabadi and R. A. King. Image coding using vector quantization: A review. IEEE Trans. on communications, 36(8), August 1988. [31] A. Nevo and E.M. Stein. A generalization of birkhoff’s pointwise ergodic theorem. Acta. Math, 173:135–154, 1994. [32] D. Ornstein and B. Weiss. The shannon-mcmillan-breiman theorem for a class of amenable groups. Isreal Journal of Mathematics, 44(1), 1983. 180 [33] D. Ornstein and B. Weiss. Entropy and isomorphism theorem for actions of amenable group. J. D’Analyse Mathematique, 48, 1987. [34] D. Ornstein and B. Weiss. Entropy and recurrence rate for stationary random fields. IEEE Trans. Inf. Theory, 48(6):1694–1697, June 2002. [35] K. Ratakonda and N. Ahuja. Lossless image compression with multiscale seg- mentation. IEEE Transactions Image Processing, 11(11):1228–1237, 2002. [36] J. B. Rosen. Pattern separation by convex programming. J. Math. Anal. Appl., 10:123–134, 1965. ISSN 0022-247x. doi: 10.1016/0022-247X(65)90150-2. URL http://dx.doi.org/10.1016/0022-247X(65)90150-2. [37] D. J. Rudolph. Fundamentals of measurable dynamics: ergodic theory on Lebesgue spaces. Oxford University Press, 1990. [38] Ashish Sen and Muni Srivastava. Regression Analysis:Theory, Methods, and Applications. Springer, 1990. [39] Alex J. Smola and Bernhard Schölkopf. A tutorial on support vector regres- sion. Stat. Comput., 14(3):199–222, 2004. ISSN 0960-3174. doi: 10.1023/ B:STCO.0000035301.49549.88. URL http://dx.doi.org/10.1023/B:STCO. 0000035301.49549.88. [40] Vladimir Vapnik, Steven E. Golowich, and Alex Smola. Support vector method for function approximation, regression estimation, and signal processing. In Ad- vances in Neural Information Processing Systems 9, pages 281–287. MIT Press, 1996. [41] X. Wu. Context-based, adaptive, lossless image coding. IEEE Trans. Commu- nications, 45(4), April 1997.