Gibbs field approach and stochastic optimization in image processing and bioinformatics

1. Gibbs field approach and stochastic optimization in image processing and bioinformatics Elena Zhizhina (Dobrushin’ Laboratory, Institute for Information Transmission Problems, Moscow) St.Peterburg May 11-15, 2010

2. Probabilistic approach in image analysis The basic idea of probabilistic approach in image analysis ( J. Besag, Spatial interaction and the statistical analysis of lattice systems (with discussion). J. Roy. Statist. Soc. Ser.B 36, pp. 192-236, 1974, Geman S., Geman D., Stochastic Relaxation, Gibbs Distribution, and the Bayesian Restoration of Images // IEEE Trans. Pattern Anal. Machine Intelligence. V. 6. no 6. P. 721-741, 1984) was to rewrite an image processing procedure in the language of statistical physics using concepts of statistical ensembles, equilibrium and non-equilibrium dynamics. Under this view, images are considered as configurations of a Gibbs field. The implicit assumption behind the probabilistic approach in image analysis is that, for a given problem, there exists a Gibbs field such that its ground states represent regularized solutions of the problem.

3. A model in the framework of Gibbs field approach is defined by • The space of images = Configuration space = The space of realization of a random field • Energy function is a sum of two terms: the first term represents a priori knowledge on a general structure and attributes of images (the energy of interaction), the second term depends on the data (the energy of an external non-homogeneous attractive field) The choice of the energy function depends on a concrete problem, and it is generally guided by experience rather than by formal methods for model-fitting. • Stochastic dynamics (evolution of configurations) as a basis for the stochastic optimization scheme

5. A space of images can be described as a space of realizations of a random field. Following Bayesian setting we search for realizations with maximum of the posterior distribution under given data . If the posterior distribution has the Gibbs form the model is called a Gibbs filed model. Gibbs field approach. I. General setting.

6. Gibbs field approach. II. Global optimization problem. Gibbs field models are defined by energy functions The goal is to find ground states of the model because Stochastic algorithms have been proposed to solve the problem of global optimization.

7. The very high dimensionality of images, as well as the non-convexity of the energy , usually excludes any direct and deterministic method for the optimization. At the same time, the local interaction allows to use stochastic iterative algorithms involving local changes at each step. In these algorithms, the resulting image is constructed as the limit configuration of a stochastic iterative procedure. At each iteration, the new configuration is obtained according to a transition distribution which depends on the current configuration. Using the local interaction property, computations of the transition probabilities become also local. In this connection, a choice of stochastic dynamics maximum adapted to a specific problem under consideration is a crucial step in the construction of the algorithm.

8. Gibbs field approach.III. Stochastic algorithms. For a Gibbs measure with density there exist stationary stochastic processes on such that the evolution of any initial state for a long time tends to a set of typical configurations of the Gibbs measure and

9. Configurations from are typical for The annealing (or cooling) procedure describes how to manage the convergence to the required configurations Gibbs field approach.IV. Simulated annealing.

10. (spin is defined at each pixel) The configuration space S is the spin space at each pixel (grey level intensity) The energy function The Metropolis-Hastings dynamics (or any spin flip type dynamics) with a slowly decreasing cooling parameter Pixel-wise models

11. PALSAR

12. The energy function. I. Denoising problem data term restoration term The data term contains information from the observation. The restoration term models a priori information on a general structure and attributes of images. For example, we say that the “true configuration” is locally constant or smooth, and it doesn’t have frequent sharp discontinuities and indented boundaries. Usually constraints are realized in a soft form: no prohibition but some penalty (any term of the energy is not equal to ). The solution of the problem is a configuration meeting general prior conditions and at the same time well-matched to the data.

13. II. Segmentation (classification) problem Spin K is the number of classes, and each class is characterized by a label Classification term controls a length of a boundary between partitions. data term classification term

14. Marked point models(random set of macro-objects) Recently, there has been again growing interest in the application of Gibbs point fields model to feature extraction, object detection, surface reconstruction, stereo matching problems. These problems become critical in remote sensing with the development of high resolution sensors for which the object geometry is well defined. All these problems related with consideration of strong geometrical constraints in a priori potential. Any geometrical properties can easily be introduced into the MP model through the object geometry. Different types of objects (trees, roads, buildings, etc.) can be considered within the same model with appropriate interactions. Moreover, interactions between points (objects) allows to model some prior information on the object configuration, and the data are taken into account at the object level, thus improving robustness of the algorithms.

15. OrbView-3

16. Marked point models. Configurations If we denote by the set of all point configurations from a finite volume G, by S a space of marks (a spin space) and by the Poisson measures with activity z, then the marked configuration space of the model is

17. Marked point models. Reference measure. and a reference measure on can be written as is the conditional (under given configuration for positions of marks) free marks measure equals to the product of the free mark measures over all points from the configuration .

18. Marked point models. Interactions. The probability distribution on the configuration space is defined then as a Gibbs reconstruction of the reference measure with the energy function involving both objects positions and their marks. To find global minimizers of the energy function, one can consider various stochastic dynamics with a given stationary Gibbs measure under the annealing procedure.

19. Geometrical constraints: the configuration space is a random set of grains (discs) with centers at The energy function The birth-and-death process with non-homogeneous intensities for birth and death. Approximation in time. Marked point models

20. A new stochastic algorithm for object detection. The main idea behind our algorithm is to use the continuous time stochastic dynamics (a stationary, even reversible, process) and then to take the transition operator of the discrete time approximation of the process as a base of stochastic iterative steps of the algorithm (a non-stationary Markov chain). Descombes X., Minlos R.A., Zhizhina E., Object extraction using stochastic birth-and-death dynamics in continuum, Journal of Mathematical Imaging and Vision, Vol. 33(3), p. 347, 2009. Descombes X., Minlos R.A., Zhizhina E., Object extraction using stochastic birth-and-death dynamics in continuum, INRIA research report RR-6135, March 2007; https://hal.inria.fr/inria-00133726 Descombes X., Zhizhina E., The Gibbs fields approach and related dynamics in image processing, Condensed Matter Physics, 2008, Vol.11, No.2, p. 1-20.

21. Generators of the processes A continuous time equilibrium dynamics with intensities Transition operator of the approximation process

22. A new Multiple Birth and Death (MBD) algorithm. The algorithm is defined as follows: Computation of the birth map: To speed up the process, we consider a non homogeneous birth rate B(s) to favor birth where the data term is strong. Main program: Initialise the inverse temperature parameter and the discretization step and alternate birth and death steps. • Birth step is taken with density w.r.t. the Lebesgue measure on V. • Death step: for each point from the configuration, the death probability is defined as follows: • Decrease the temperature and the discretization step by a given factor and go back to the birth step.

23. Convergence of the approximation process to the continuous time process with the generator for all uniformly on bounded intervals of time Here Convergence of the approximation processes (double annealing procedure) is a measure concentrated on global minima of the energy functional

24. Plantation recognitionleft: initial image (provided by French National Forest Inventory); right: detected trees

25. The energy function. Example 1. A priori term. To optimize the overlap between objects and to get a flexible model w.r.t. the data we have only to penalize overlapping discs. The prior term is defined by a pair interaction The total energy is

26. The energy function for object detection. The data term. The data term gives penalty for a bad fitting of the current disc configuration onto the data. The detection is based on a contrast between crown and boundary areas on an image. Let be empirical means and variances associated with internal area and external boundary area of a disc superimposed on the data . We introduce a distance between two distributions associated with above parameters: Then the data term is a decreasing function of the distance .

27. The data term

28. Forest trees recognitionleft: initial image (provided by French National Forest Inventory); right: detected trees

29. Flamingo populationStation Biologique Tour du Valat

30. Flamingo population: detected birds

34. Road detection using multiple birth and death (MBD) optimization algorithm • Configuration space = space of discs with small radius • Energy function is the sum of three terms: - a priori term to model connectivity (by pair potential with non-zero attractive part on fitting distances), and a prior term for managing curvature and junctions of the road network (by multi-particle interactions); - the data driven term is defined through the gradient field of the data digital image (a local road detector) • Optimization using MBD algorithm

38. Gibbs field approach for evolutionary analysis of regulatory signal of gene expression under constraints on secondary structure, V.A. Lyubetsky, E.A. Zhizhina, L.I. Rubanov, Problems of Information Transmission, 2008, Vol. 44, pp. 333-351. We propose a new approach to modeling of a nucleotide sequence evolution in mRNA, subject to constraints on secondary structure. Our approach is based on the optimization problem for a functional that involves both standard evolution of primary structure and a condition of the secondary structure conservatism. The secondary structure here is specified by a sophisticated potential of non-local interaction We discuss results of simulations by the example of evolution for classical attenuation regulation case.

39. The model • The graph is assumed to be known: a finite phylogenetic tree with given phylogenetic lengths proportional to the time of evolution. • A spin space is complicated. Each spin is a long sequence (100-200 symbols) in 4-letter (nucleotides) alphabet {A,C,T,G}, and in addition they should have a specific secondary structure. • The sequences on the tree leaves (current sequences) are fixed. • The problem is to reconstruct ancestral sequences at all inner nodes of the tree in such way that under theevolution the secondary structure is conserved (as a structure!) along the way from leaves to the root of the tree. • The optimization scheme is a MH type dynamics with 3 types of spin modification in one iteration step: a letter substitution at one position, insertion and deletion.

40. Example (classical attenuation regulation of threonine biosynthesis in gamma-proteobacteria). The standard tree of species has 27 nodes, among them 14 leaves, and every edge is assigned a phylogenetic length in conventional units. The leaves are marked with abbreviated names of species as follows: EC -- Escherichia coli, TY -- Salmonella typhi, KP -- Klebsiella pneumoniae, EO -- Erwinia carotovora, YP -- Yersinia pestis, HI -- Haemophylus influenzae, VK -- Pasterella multocida, AB -- Actinobacillus actinomycetemcomitans, PQ -- Mannheimia haemolytica, VC -- Vibrio cholerae, VV -- Vibrio vulnificus, VP -- Vibrio parahaemolyticus, SON -- Shewanella oneidensis, XCA -- Xanthomonas campestris.

41. Example of a spin: a sequence composed of letter in the 4-letter alphabet {A, C, T, G} (the alphabet of nucleotides). This sequence is the regulatory signal for classical attenuation regulation: а) sequence with two shoulders of the anti-terminator helix marked, state A of the regulation site; b) same sequence with two shoulders of the terminator helix marked, state T of the regulation site. T-run is shown in capital letters; c) sequence of the leader peptide gene that contains start and stop codons shown in underlined capital letters, and 12 regulatory codons shown in capital letters. The sequences a) and b) continue the sequence с) starting from the position marked by arrow.

42. Attenuation regulation. Primary and secondary structure. Attenuation regulation is based on the possibility of forming alternative secondary structure such that one structure allows protein synthesis and the other prohibits it. We can see Terminator helix and Anti-terminator helix. Each helix is two helix shoulders which formed a sequence of complementary pairs of nucleotides (G-C, T-A, and G-T with some constraints). The right A-shoulder and the left T-shoulder have non-zero intersection, that means that only one helix (T or A) could be closed at the moment. The sequence itself is called the primary structure and the sequence along with helices is called secondary structure.

43. reflects the energy of pair interaction in the system of spins (for the align sequences) reflects the dependency on data given at leaves reflects the requirement of conserved secondary structure along each edge.

45. )

46. Conclusions Our analysis of the composition of minimal configurations (ground states of our model) under different shows that in the domain of “moderate” values a strong regulatory structure of one type is kept along the whole tree of evolution. In this case, we can observe paths from almost every leaf to the root that conserves a secondary structure. When or small enough, i.e. only a primary structure of the evolving sequence remains significant, all tests show lack of paths with a secondary structure conserved from leaves to the root. Finally, for large when the primary structure interaction becomes less significant in our functional, the composition of minimal configurations changes again. It appears as a sharp decreasing number of long paths with conserved secondary structure from leaves to the root: we can see only pieces of these long paths, more energy-favorable for the corresponding energy term .