Parameter Learning for Log-supermodular Distributions

We consider log-supermodular models on binary variables, which are probabilistic models with negative log-densities which are submodular. These models provide probabilistic interpretations of common combinatorial optimization tasks such as image segmentation. In this thesis, we focus primarily on parameter estimation in the models from known upper-bounds on the intractable log-partition function. We show that the bound based on separable optimization on the base polytope of the submodular function is always inferior to a bound based on "perturb-and-MAP" ideas. Then, to learn parameters, given that our approximation of the log-partition function is an expectation (over our own randomization), we use a stochastic subgradient technique to maximize a lower-bound on the log-likelihood. This can also be extended to conditional maximum likelihood. We illustrate our new results in a set of experiments in binary image denoising, where we highlight the flexibility of a probabilistic model to learn with missing data. This chapter is based on the work "Parameter Learning for Log-supermodular Distributions", T.

Abstract

Probabilistic graphical models encode hidden dependencies between random variables for data modelling. Parameter estimation is a crucial and necessary part of handling such probabilistic models. These very general models have been used in plenty of fields such as computer vision, signal processing, natural language processing and many more. We mostly focused on log-supermodular models, which is a specific part of exponential family distributions, where the potential function is assumed to be the negative of a submodular function. This property will be very handy for maximum a posteriori and parameter learning estimations. Despite the apparent restriction of the models of interest, they cover a broad part of exponential families, since there are plenty of functions that are submodular, e.g., graph cuts, entropy and others. It is well known that a probabilistic treatment is challenging for most models, however we were able to tackle some of the challenges at least approximately.

In this manuscript, we exploit perturb-and-MAP ideas for partition function approximation and thus efficient parameter learning. Moreover, the problem can be also interpreted as a structure learning task, where each estimated parameter or weight represents the importance of the corresponding term. We propose a way of approximating parameter estimation and inference for models where exact learning and inference is intractable in general case due to the partition function calculation complexity.

The first part of the thesis is dedicated to theoretical guarantees. Given the log-supermodular models, we take advantage of the efficient minimization property related to submodularity. Introducing and comparing two existing upper bounds of the partition function, we are able to demonstrate their relation by proving a theoretical result. We introduce an approach for missing data as a natural subroutine of probabilistic modelling. It appears that we can apply a stochastic technique over the proposed perturb-and-map approximation approach and still maintain convergence while make it faster in practice.

The second main contribution of this thesis is an efficient and scalable generalization of the parameter learning approach. In this section we develop new algorithms to perform parameter estimation for various loss functions, different levels of supervision and we also work on the scalability. In particular, working with mostly graph cuts, we were able to incorporate various acceleration techniques.

As a third contribution we deal with the general problem of learning continuous signals. In this part, we focus on the sparse graphical models representations. We consider common sparsity-inducing regularizers as negative log-densities for the prior distribution. The proposed denoising techniques do not require choosing any precise regularizer in advance. To perform sparse representation learning, the signal processing v community often uses symmetric penalties such as ℓ 1 , but we propose to parameterize the penalty and learn the weight of each loss component from the data. This is feasible via an approach which is similar to what we proposed in the previous sections.

For all aspects of the parameter estimation mentioned above we performed computational experiments to illustrate the ideas or compare with existing benchmarks, and demonstrate its performance in practice.

vi

Résumé

Les modèles graphiques probabilistes codent des dépendances cachées entre des variables aléatoires pour la modélisation des données. L'estimation des paramètres est une partie cruciale et nécessaire du traitement de ces modèles probabilistes. Ces modèles très généraux ont été utilisés dans de nombreux domaines tels que la vision par ordinateur, le traitement du signal, le traitement du langage naturel et bien d'autres. Nous nous sommes surtout concentrés sur les modèles log-supermodulaires, qui constituent une partie spécifique des distributions de la famille exponentielle, où la fonction potentielle est supposée être l'opposée d'une fonction sous-modulaire. Cette propriété sera très pratique pour l'estimation par maximum a posteriori et l'apprentissage des paramètres. Malgré la restriction apparente des modèles d'intérêt, ils couvrent une grande partie des familles exponentielles, puisqu'il y a beaucoup de fonctions qui sont sous-modulaires, par exemple, les coupes dans de graphes, l'entropie et d'autres. Il est bien connu que le traitement probabiliste est un défi pour la plupart des modèles, mais nous avons été en mesure de relever certains de ces défis au moins approximativement.

Dans ce manuscrit, nous exploitons les idées "perturb-and-MAP" pour l'approximation de la fonction de partition et donc un apprentissage efficace des paramètres. De plus, le problème peut être considéré comme une tâche d'apprentissage structuré, où chaque paramètre ou poids estimé représente l'importance du terme correspondant. Nous proposons une méthode d'estimation et d'inférence approchée des paramètres pour les modèles où l'apprentissage et l'inférence exacts sont insolubles dans le cas général en raison de la complexité du calcul des fonctions de partition.

La première partie de la thèse est consacrée aux garanties théoriques. Étant donnés les modèles log-supermodulaires, nous tirons parti de la propriété de minimisation efficace liée à la sous-modularité. En introduisant et en comparant deux bornes supérieures existantes de la fonction de partition, nous sommes en mesure de démontrer leur relation en prouvant un résultat théorique. Nous introduisons une approche pour les données manquantes comme sous-routine naturelle de la modélisation probabiliste. Il semble que nous puissions appliquer une technique stochastique à l'approche d'approximation par perturbation tout en maintenant la convergence et la rendant plus rapide en pratique.

La deuxième contribution principale de cette thèse est une généralisation efficace de l'approche de l'apprentissage paramétrique. Dans cette section, nous développons de nouveaux algorithmes pour effectuer l'estimation des paramètres pour diverses fonctions de perte, différents niveaux de supervision et travaillons sur l'efficacité vii à grande échelle. En particulier, en travaillant principalement avec des coupes de graphes, nous avons pu intégrer différentes techniques d'accélération.

Comme troisième contribution, nous traitons d'un problème général d'apprentissage des signaux continus. Dans cette partie, nous nous concentrons sur les représentations de modèles graphiques parcimonieux. Nous traitons des penalités classiques en estimation parcimonieuse en les considérant comme l'opposée de la log-vraiseullante de la distribution a priori. Les techniques de débruitage proposées ne nécessitent pas de paramêtre de régularisation précis à l'avance. Pour effectuer l'apprentissage de la représentation parcimonieuse la communauté utilise souvent des pénalités symétriques comme la meme ℓ 1 , mais nous proposons de paramétrer la perte et d'apprendre le poids de chaque composante de la pénalité à partir des données. C'est faisable avec l'approche sur laquelle nous avons travaillé auparavant.

Pour tous les aspects de l'estimation des paramètres mentionnés ci-dessus, nous avons effectué des expériences pour illustrer l'idée ou la comparer aux algorithmes existants, et démontrer sa performance en pratique. viii I am grateful that I had a chance to work on my PhD thesis under Francis Bach supervision during these three and a half years. Francis gave me the best thing he could, the enthusiasm and belief in academical research. I am thankful that you introduce the exciting top-tier scientific world to me and I hope I was a good student to accept and succeed. Thank you for making my dream to become a doctor come true and encouraging me to continue my academic career. I will always be grateful for that.

I also acknowledge support the European Union's H2020 Framework Programme (H2020-MSCA-ITN-2014) under grant agreement n o 642685 MacSeNet.

I would like to acknowledge Florence d'Alché-Buc and Matthew Blaschko for reviewing this thesis. I am also thankful to the jury members Ivan Laptev and Umut Şimşekli who agreeing to participate. I deeply thank you all for the correction and remarks that I have got after your careful reading.

I would like also to mention Mike Davies here, who was not my co-advisor, but I often felt that he was. He teached me a different side of view and showed himself in the best qualities during our collaborative time. From Mike I have also learned to be sensitive to the details which is always important.

I am pleased that Anton Osokin and Dmitry Ostrovsky were always there and thank them for collaborations and interesting discussions related to this thesis. I would like to express my deepest gratitude for their continuous help as well.

I would like to highlight that SIERRA team is flourishing space for personal grow and idea sharing. The time I spent in SIERRA/WILLOW environment was the most productive in my life. The concentration of motivation and smartness was the highest ever and I am proud that I was a part of it. I am thankful for everyone who accompanied me on my path to the graduation during all these years.

I would especially like to thank my partner and college from SIERRA Dmitry Babichev for all his unlimited love and support, and thanks for pushing me forward sometimes. This PhD journey would be half less joyful and half harder without you by my side. I am thankful for my family as well, who always supports me, though probably do not have an idea what is this manuscript about.

There are many more people who had influenced me to become a researcher and conduct a PhD. I hope all of you feel my gratitude.

Contributions and thesis outline

Parameter estimation is a big part of probabilistic graphical model handling.

Given the optimal parameter one can perform inference and evaluate marginal and conditional probabilities. Optimal parameters are often considered as a solution of the maximum likelihood problem, which is log-concave for exponential family distributions. However, the partition function term is hard to compute, and thus causes an unsolvable problem for the maximum likelihood learning process in the general case, and thereby we propose to use an effective and tight upper bound of the partition function based on the Gumbel perturbations and apply it for the broad but specific log-supermodular family of distributions.

Thus, this thesis brings new light on parameter learning for challenging probabilistic graphical models. We mostly provide some guarantees on the approach and develop new directions and possible applications in the rest of the thesis. Below we summarize each of the chapter contribution.

Chapter 1: We introduce some basic mathematical concepts that will be essential for this work. We provide an overview of probabilistic graphical models and its properties, submodularity and based on it log-supermodular distributions, and we also cover challenges and suitability of parameter learning for these models of interest.

Chapter 2: As the first contribution of this thesis, we present a new automatic parameter estimation procedure. Based on the idea of efficient maximum likelihood optimization, we propose to incorporate a suitable upper bound on the partition function, and thus make the parameter estimation problem feasible. The only assumption required is the existence of a maximum a posteriori solver. The sufficient condition for this is to consider the log-supermodular models, for which we can show the superiority over existing approaches.

Chapter 3: The next contribution we present is a generalization in terms of scalability, levels of supervision and various performance evaluation functions. To be specific, our intention was to make our approach more practical: we worked on the real losses like Hamming and weighted Hamming losses, we worked on a large-scale image segmentation problem with unlabeled data or only bounding-box given data, as these both cases are more common than the fully-labeled ones.

Chapter 4: This chapter considers the problem of a sparsity-inducing density learning supported on the continuous domain. Basically, we propose to consider the learned density function as a prior that encourages signal sparsity and thereby incorporate the derived prior into decoding, e.g., denoising problem.

Chapter 5: In this chapter we conclude the thesis, list the main contribution Chapter 1

Introduction

In this chapter we introduce in details the topic of probabilistic graphical models and the associated subroutines of the parameter learning problem. We will outline log-supermodular models which are crucial for this thesis, and discuss submodularity along with its properties.

Probabilistic Graphical Models

Probabilistic graphical models are a versatile tool proposed for modeling structured data and capturing the dependencies between variables via concepts from graph theory. Random variables of interest (𝑋 1 , . . . , 𝑋 𝑛 ) = 𝑋, where 𝑋 is distributed according to a distribution 𝑝(𝑥), are dependent and connected in a predefined way. The corresponding graphical model can be described by a notation of a graph 𝐺 = (𝑉, 𝐸) and encode the random variables via the node set 𝑉 and the pairwise dependencies between them via the edge set 𝐸. A probabilistic graphical model takes advantage of probability theory to make a consistent probabilistic interpretation given the values of random variables (Wainwright and Jordan [2008b]). The goal would be to manipulate the model and evaluate various marginal and conditional probabilities. There are mainly two broad families: undirected and directed graphical models. We will briefly cover both, however undirected graphical models (they are also called Markov Random Fields) are more in demand for our needs in this work. In general there are all four types of distributions: that can be described by both UGM and DGM, or by only one of them, or by neither of them. Thus, none of two GM is more powerful or general than the other (see [START_REF] Murphy | Machine learning : a probabilistic perspective[END_REF] for more details). There are also models based on the chain graphs, graphs that include both ideas of undirected and directed graphical model. We present several examples in Figure 1-1.

Markov Random Fields

Undirected graphical models assume a representation where the local dependencies between the nodes are symmetric. This assumption makes sense for some domains, for example for images or other symmetric structured data as relational or spatial data. Another common way to model data is to use directed graphical models which are considered conditional dependencies and based on directed graphs. It is worth noting, that for DGMs is natural to decompose the density function using the chain rule, where each factor is a conditional probability, however is not natural for UGMs anymore. In this case we have to additionally normalize the model via the crafted normalizer term which is called the partition function. Thus, it becomes much harder to manipulate the model.

Conditional Random Fields

The general concept of Conditional Random Fields (CRFs, see [START_REF] Murphy | Machine learning : a probabilistic perspective[END_REF]) is that instead of working with joint probability distribution, we model the a-posteriori distribution, conditioned on the feature data. To make it more concrete, let us introduce the observed variables 𝑋 = (𝑋 1 , . . . , 𝑋 𝑛 ) which is called the evidence and the corresponding unobserved variables 𝑌 = (𝑌 1 , . . . , 𝑌 𝑛 ) which are labels. Thus, using a discriminative model 𝑝(𝑦|𝑥) instead of a generative model 𝑝(𝑥, 𝑦), we can consider the labels-to-predict behaviour only. This problem statement gives us more flexibility and allows us to focus on the target variables only. See Figure 12for the example of a CRF, where the shaded nodes are meant to be the observed variables.

X 1 X 2 X n-1 X n Y 1 Y 2 Y n-1 Y n . . . . . .

Conditional Independence

By now the general idea was defined intuitively. Here we are going to discuss in more details the UGM's parameterization. As the chain rule can not be applied in the undirected case, there is a way to formulate parametrization formally with the help of the Hammersley-Clifford theorem [START_REF] Lauritzen | Graphical models[END_REF]).

Before formulating this, it is necessary to introduce the concept of conditional independence. Conditional independence's goal is to state whether two sets of variables are independent or not conditioned on a third set of variables. This can be decided by means of factoring the distribution into pieces or factors. We have already said that we assume some local dependencies between variables, as, e.g., in Figure 1-1, but now we would like to define it clearly. We illustrate this definition with the following example in Figure 1-3. It is worth noting, that two nodes which are not connected with an edge, are conditionally independent given the rest of the graph. The Markov property makes concrete the intuitive statement that variables interact locally on the graph.

X A X B X C

Model Parameterization

There is a common way for the undirected graphical model parametrization that we present next. We associate a potential function to each region or maximal clique in the graph and can formalize the probability mass function 𝑝(𝑥 1 , . . . , 𝑥 𝑛 ) = 𝑝(𝑥), where 𝑥 𝑖 is a realization of a discrete variable 𝑋 𝑖 , by means of the Hammersley-Clifford theorem [START_REF] Hammersley | On the partition function and random maximum aposteriori perturbations[END_REF]).

Theorem 1.1 (Hammersley-Clifford). Any positive distribution 𝑝(𝑥) > 0 that satisfies the Markov property of an undirected graph 𝐺 = (𝑉, 𝐸) can be expressed via the product of factors, where one factor depends on one maximal clique of the graph:

𝑝(𝑥) = 1 𝑍 ∏ 𝑐∈𝐶 𝜑 𝑐 (𝑥 𝑐 ),
where 𝐶 is a set of all maximal cliques of the graph 𝐺 and 𝑍 =

∑ 𝑥∈X ∏ 𝑐∈𝐶 𝜑 𝑐 (𝑥 𝑐 ) is the normalization constant.
Any non-negative function can serve as a potential function 𝜑 𝑐 . We are going to represent the probability density function via the Gibbs distribution using logpotentials 𝑓 (𝑥) =

∑ 𝑐∈𝐶 𝑓 𝑐 (𝑥 𝑐 ) = ∑ 𝑐∈𝐶 log 𝜑 𝑐 (𝑥 𝑐 ).
Thus, for a MRF it has the following form: 𝑥) .

𝑝(𝑥) = 𝑒 𝑓 (𝑥) 𝑍 = 𝑒 𝑓 (𝑥) ∑ 𝑥∈X 𝑒 𝑓 (
And a CRF model can be written as:

𝑝(𝑦|𝑥) = 𝑒 𝑓 (𝑦|𝑥) 𝑍(𝑥) = 𝑒 𝑓 (𝑦|𝑥) ∑ 𝑦∈Y 𝑒 𝑓 (𝑦|𝑥) .
While the interpretation with maximal cliques is classic but rather heavy, in practice pairwise models are more popular, where all cliques are of size one or two nodes, and thus the probability density function is written as:

𝑝(𝑥) ∝ ∏ 𝑖,𝑗∈𝐸 𝜑 𝑖,𝑗 (𝑥 𝑖 , 𝑥 𝑗 ) ∏ 𝑖∈𝑉 𝜑 𝑖 (𝑥 𝑖 ).
This form is used due to its simplicity and strong coverage of various models. Another common treatment is to represent the log potentials as a linear function of some parameters 𝜃: 𝑓 (𝑥|𝜃) = 𝜓(𝑥) 𝑇 𝜃, where 𝜓(𝑥) is a feature vector of the data vector 𝑥. In many applications these features are created by hand, depending on the domain knowledge and user's needs. We will cover this part in the experimental settings. Thus, our final model has the form:

log 𝑝(𝑥|𝜃) = 𝜓(𝑥) 𝑇 𝜃 -log 𝑍(𝜃),
which is called the log-linear model and the potential vector 𝜓(𝑥) contains only pairwise and unary terms.

The distribution of interest 𝑝(𝑥|𝜃) clearly belongs to the exponential family as it can be rewritten as 𝑝(𝑥|𝜃) = exp[𝜓(𝑥) 𝑇 𝜃 -𝐴(𝜃)], where 𝜓(𝑥) is called a vector of sufficient statistics and 𝐴(𝜃) = log 𝑍(𝜃) is a log-partition function. The distribution enjoys the exponential family's properties, most importantly it attains maximum entropy, and thus makes the weakest assumptions about the data [START_REF] Murphy | Machine learning : a probabilistic perspective[END_REF]). Besides that, using the probability density function allows us to make use of the log-likelihood function concavity. For more details see Wainwright and Jordan [2008b], [START_REF] Murphy | Machine learning : a probabilistic perspective[END_REF], [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF].

Examples

In this part we will present some prominent examples of Markov Random Fields like Ising and Potts Models, and Gaussian MRFs.

Ising Model. This one is a pioneer and one of the classical graphical models created for the statistical physics needs [START_REF] Ising | Beitrag zur theorie des ferromagnetismus[END_REF]). The idea is to model the system of directed magnetic spins, which are directed in one of the two possible directions: 𝑥 𝑖 ∈ {+1, -1}. This can be written down in the following way:

log 𝑝(𝑥|𝑊, 𝑏) = ∑ 𝑖∼𝑗 𝑤 𝑖𝑗 𝑥 𝑖 𝑥 𝑗 + ∑ 𝑖 𝑏 𝑖 𝑥 𝑖 -log 𝑍(𝑊, 𝑏) = 1 2 𝑥 𝑇 𝑊 𝑥 + 𝑏 𝑇 𝑥 -log 𝑍(𝑊, 𝑏).
A sign of the weight 𝑤 𝑖𝑗 will influence the tendency of the two connected nodes to have the same or the opposite sign. It is worth noting, that the normalization constant 𝑍(𝑊, 𝑏) is a sum of potential functions of all possible values of vector 𝑥, i.e., 2 𝐷 terms, and its calculation is #𝑃 hard in general case [START_REF] Jerrum | Polynomial-time approximation algorithms for the Ising model[END_REF]). A variation of this model will be covered in our experiments.

Potts Model.

It is a straightforward generalization of the Ising Model for multiple number of discrete states: 𝑥 𝑖 ∈ {0, 1, . . . , 𝐾} [START_REF] Potts | Some generalized order-disorder transformations[END_REF]). This model will also be covered in the experimental section of Chapter 4. Both Potts and Ising models can be applied for an image segmentation problem, since it is feasible to encourage an images smoothness using these types of models [START_REF] Boykov | An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision[END_REF]). This usually is done via Conditional Random Fields (CRFs discussed before), where the features variables 𝑥 𝑖 are observed and the labels 𝑦 𝑖 are not. We present a grid model common for image segmentation problems in the Figure 1234. Gaussian MRFs. In the case of Gaussian MRFs, the joint probability density function is assumed to be Gaussian and has the following form:

𝑝(𝑥|𝜇, Σ) ∼ N(𝜇, Σ).
If Σ -1 𝑖,𝑗 is zero, then there is no edge between nodes 𝑖 and 𝑗 and as follows no pairwise term of 𝑖 and 𝑗 (this follows from the conditional independence properties).

Probabilistic Inference

Given the probabilistic model we can perform inference with the tools provided by probabilistic theory. Using the axioms of probability, any marginal or conditional probability of interest can be defined, but not all calculations can be done in polynomial time. Let us discuss this problem of computing in more details as this knowledge will be crucial further in this thesis.

Let us consider the conditional probability of a subset 𝑋 𝐴 given the values for a subset 𝑋 𝐵 , while the rest random variables of the graphical model we denote as 𝑋 𝐶 . Our goal would be to evaluate the conditional probability 𝑝(𝑋 𝐴 |𝑋 𝐵 ), which is the classical problem of probabilistic inference.

X A X B X C Figure 1-5:
An example of the conditional probability calculation. The dark nodes correspond to the nodes on which we condition, light-shaded nodes are the nodes we marginalize over, and the white nodes are the nodes we want to calculate conditional probability for.

To make this easy to present, we will use Figure 1-5 with 3 shades of grey for the nodes coloring. The target variables 𝑋 𝐴 are unshaded, the observed variables 𝑋 𝐵 , which are also called the evidence nodes, are dark shaded, while light shaded are the rest nodes values 𝑋 𝐶 and they must be marginalized out to calculate 𝑝(𝑋 𝐴 |𝑋 𝐵 ). We derive the conditional probability in the following way:

𝑝(𝑋 𝐴 |𝑋 𝐵 ) = 𝑝(𝑋 𝐴 , 𝑋 𝐵 ) 𝑝(𝑋 𝐵 ) = ∑ 𝑋 𝐶 𝑃 (𝑋 𝐴 , 𝑋 𝐵 , 𝑋 𝐶 ) ∑ 𝑋 𝐴 ,𝑋 𝐶 𝑃 (𝑋 𝐴 , 𝑋 𝐵 , 𝑋 𝐶 ) .
Thus, the point of discussion in this section is an efficient calculation of this expression, which contains two expensive summations and can be of exponential complexity.

Exact Inference. There are plenty of algorithms for exact inference for the discrete graphical models. The most naive approach would be an elimination algorithm, that attains at a conditional probability of a single node (Wainwright and Jordan [2008b]), which is a quite specific problem to solve. However, a good point is that this approach is applicable for any graphical model. The most well-spread set of techniques are called message passing approaches. For example, the follow-up of the elimination algorithm is the belief propagation or sum-product message passing algorithm that can be applied for any acyclic graphical model, i.e., tree-based ones (Wainwright and Jordan [2008b]), and it can compute all single-node marginals altogether. The next step of generalization would be the junction tree algorithm (Wainwright and Jordan [2008b]), which is mostly the combination of the previous two algorithms.

All of these algorithms are restricted either to a special graph structure, mostly to acyclic graphs or to the low-treewidth cases (the computation time is exponential in the treewidth of the graph, see more details in Wainwright and Jordan [2008b]). Moreover, these approaches do not work for continuous variables, except for jointly Gaussian distributions. It is explainable by the fact that exact inference is NP-hard in general (Wainwright and Jordan [2008b]). To perform general inference, one needs to call for approximated techniques covered in the next paragraph.

Approximate Inference. When exact inference is often intractable or too inefficient, an alternative would be approximate inference. It is interesting to note, that the sum-product message-passing algorithm from the previous paragraph can also be seen as an approximate inference technique applied to a cycled factor graph [START_REF] Szeliski | A comparative study of energy minimization methods for markov random fields with smoothness-based priors[END_REF]). Here we aim to cover the most prominent of concepts like variational and Monte Carlo inference.

Variational Inference. The idea of variational inference is to approximate the intractable distribution of interest 𝑝(𝑥) with another close, but tractable distribution 𝑞 (𝑥). Some examples of 𝑞(𝑥) could be from factored or Gaussian distributions. We usually evaluate the "closeness" via the KL divergence KL(𝑞||𝑝) = ∑ 𝑥 𝑞(𝑥) log 𝑞 (𝑥) 𝑝(𝑥) and, thus, can fit free parameters of 𝑞(𝑥) by minimizing the KL objective.

One of the most classic variational inference approaches is called "mean field" approximation, where the idea is to use a fully factorized posterior as 𝑞(𝑥)

= 𝑛 ∏ 𝑖=1 𝑞(𝑥 𝑖 )
approximation, which is a product of marginals (Wainwright and Jordan [2008b]). Thus, the KL optimization problem is solvable by the coordinate descent method. The mean field approach belongs to the variational Bayes family of approaches, which aim to find the optimal parameters of the distribution 𝑞(𝑥) [START_REF] Murphy | Machine learning : a probabilistic perspective[END_REF]). Moreover, there is also an extension: variational Bayes EM, where we aim to find both the parameters and the latent variables simultaneously [START_REF] Murphy | Machine learning : a probabilistic perspective[END_REF]).

Another milestone is the loopy belief propagation technique, which is a generalization of the message passing (Murphy and Weiss [1999]). There are also more techniques that variational inference paradigm contains, however we will not cover it here (for more details see Wainwright and Jordan [2008b]). It is worth to mention, that variational inference is always approximate even if run to convergence.

Monte Carlo Inference.

To break the bounds caused by variational inference, such as the special form of approximation or necessary properties as conjugate priors, we are going to introduce an alternative class of algorithms based on Monte Carlo approximation. The idea is to approximate any target value with averaged value

obtained from sampling E[𝑓 |D] ≈ 1 𝑁 𝑁 ∑ 𝑛=1 𝑓 (𝑥 𝑛 ).
The most prominent example for sampling in high dimensions is Markov Chain Monte Carlo (MCMC). This approach constructs a Markov chain on the state space X, where the stationary distribution matches the target distribution 𝑝(𝑥). Thus, while performing a random walk on the states and converging to the stationary distribution, we are sampling from 𝑝(𝑥) and then integrate with respect to 𝑝(𝑥) [START_REF] Robert | Monte Carlo statistical methods[END_REF]).

A special case of MCMC is Gibbs sampling, which is very similar to mean field in the sense that we should compute each node's conditional and average out the neighbors. Random walk Monte Carlo methods include such approaches as Metropolis Hastings algorithm, slice sampling and others (for more details see Wainwright and Jordan [2008b], [START_REF] Doucet | An introduction to sequential monte carlo methods[END_REF]).

Variational inference advantages would be its deterministic solution and a bound on loglikelihood it gives, but MCMC can be applied for broader types of distributions and usually is easier to implement, and, given enough time will converge to exact inference.

Submodular Functions and Log-supermodular Models

Moving away from the general theory of graphical models, in this section we would like to present a brand-new family of distributions, which is called family of log-supermodular distributions. They are practically based on the crucial notion of submodularity, and here we present its definition and properties in details.

Submodularity

We first start to consider a set function defined on the ground set 𝑉 in the form

𝐹 (𝑆) : 2 𝑉 → R,
where 𝑆 ⊆ 𝑉 and 2 𝑉 is a set of all possible subsets of the ground set 𝑉 . Optimization of this function is a crucial problem of discrete optimization. The solution could represent the optimal set of pixels for the foreground/background image segmentation problem, or optimal locations for the fire alarms. Without loss of generality we always assume that 𝐹 (∅) = 0.

To make it looks more like a discrete optimization problem, we can biject each subset 𝑆 ⊆ 𝑉 into a binary vector 𝑥 ∈ {0, 1} |𝑉 | , where 𝑥 𝑖 = 1 means that the 𝑖-th element belongs to the subset 𝑆 and |𝑉 | is the ground set cardinality. In this interpretation the target function 𝑓 (𝑥) has its values on the nodes of the hypercube {0, 1} |𝑉 | (see example in Figure 123456).

We first introduce submodular functions, which are set functions with strong properties, kind of analogous to convex functions in the continuous space. They are both very common in practice and both are feasible to be minimized, as we can minimize submodular function in polynomial time [START_REF] Fujishige | Submodular Functions and Optimization[END_REF], [START_REF] Bach | Submodular functions: from discrete to continuous domains[END_REF]).

x 1 x 3 x 2 (0, 0, 0) ∼ { } (1, 0, 0) ∼ {a} (0, 1, 0) ∼ {b} (1, 1, 0) ∼ {a, b} (0, 0, 1) ∼ {c} (1, 0, 1) ∼ {a, c} (0, 1, 1) ∼ {b, c} (
𝐹 (𝑥 + 𝑒 𝑖 ) -𝐹 (𝑥) ≥ 𝐹 (𝑦 + 𝑒 𝑖 ) -𝐹 (𝑦).
We should also introduce a concept of a supermodular function, which is simply the negative of a submodular one. Another concept convenient for our purposes is the base polyhedron, which is a polyhedron that is associated to a submodular function. This polyhedron will play a crucial role for analysis and optimization, such as probabilistic inference in the log-supermodular models.

Submodular Function Minimization

As already been mentioned, there is a polynomial time algorithm for submodular function minimization, while maximization can be tackled only approximately. To substantiate the statement, we introduce the Lovász extension, which is a convex extension that can be defined for any submodular function [START_REF] Lovász | Submodular functions and convexity[END_REF]).

𝑓 (𝑤) = |𝑉 | ∑ 𝑣=1 𝑤 𝑖𝑣 [𝐹 ({𝑖 1 , . . . , 𝑖 𝑣 }) -𝐹 ({𝑖 1 , . . . , 𝑖 𝑣-1 })],
where (𝑖 1 , . . . , 𝑖 |𝑉 | ) is the following permutation of components of 𝑤 :

𝑤 𝑖 1 ≥ 𝑤 𝑖 1 ≥ • • • ≥ 𝑤 𝑖 |𝑉 | .
It is worth noting, that the Lovász extension 𝑓 (𝑤) is convex if and only if the initial function 𝐹 (𝑥) is submodular [START_REF] Lovász | Submodular functions and convexity[END_REF]). Moreover, the solution of the Lovász extension minimization is binary (i.e., can be restricted to belong to {0, 1} |𝑉 | ) and thus appears to be the minimizer for the submodular function 𝐹 (𝑥) as well. This statement gives us a polynomial time guarantee. To sum up, we can link the submodular and convex analysis and profit from some of the convexity guarantees. This is used in machine learning both for regularizers [START_REF] Bach | Structured sparsity-inducing norms through submodular functions[END_REF]) or losses [START_REF] Yu | Learning submodular losses with the lovász hinge[END_REF]).

Another interesting connection to mention is that the Lovász function is a support function of the base polyhedron 𝐵(𝐹 ) itself, and, thus, may be computed in the closed form 𝑓 (𝑤) = min 𝑠∈𝐵(𝐹 ) 𝑤 𝑇 𝑠.

Examples of Submodular Functions

The simplest example of a submodular function would be the cardinality function.

Having a look at Definition 2, we notice the inequality always happens to be equality. It is interesting to notice, that this function is supermodular as well. The functions which are both submodular and supermodular we will call modular function and it is known that they can be respesented as 𝐹 (𝑥) = 𝑤 𝑇 𝑥, where 𝑤 ∈ R |𝑉 | is a vector of weights and 𝑥 is a binary representation of a set.

Another crucial example is a cut function, which is defined for the undirected graph 𝐺 = (𝑉, 𝐸) that builds on the ground set of nodes 𝑉 . Basically the cut function 𝐹 (𝑆) counts the number of edges connecting a node from 𝑆 and a node from 𝑉 ∖𝑆, i.e.

𝐹 (𝑆) = 𝐹 (𝑥 𝑆 ) = ∑ 𝑖∼𝑗 |𝑥 𝑆 𝑖 -𝑥 𝑆 𝑗 |,
where 𝑥 𝑆 is a binary vector representation of the subset 𝑆, and 𝑥 𝑆 𝑖 and 𝑥 𝑆 𝑗 are representations of the nodes 𝑖 and 𝑗 correspondingly. This function is easily generalized for the weighted case.

Other prominent examples of the submodular functions are network flows, concave functions of a set cardinality, mutual entropy and matroids (see more details in [START_REF] Fujishige | Submodular Functions and Optimization[END_REF], [START_REF] Bach | Submodular functions: from discrete to continuous domains[END_REF]).

Log-supermodular Distributions

Now we marry the two big parts introduced before: probabilistic models and submodular functions. Let us consider a probabilistic model defined on the binary space X = {0, 1} 𝐾 in the form of a Gibbs distribution:

𝑝(𝑥) = exp(-𝑓 (𝑥)) 𝑍(𝑓 ) ,
where now we assume that the potential 𝑓 (𝑥) is a submodular function. Hence its negative -𝑓 (𝑥) is a supermodular one and, thus, the model is called log-supermodular probabilistic model or log-supermodular distribution. It was recently introduced by Djolonga and Krause [2014] and will be investigated throughout this thesis.

This common form defines a broad family of distributions on the one side, and allow us to perform efficient parameter learning of the other side. Maximum a posteriori estimation (MAP) for these models is known to be feasible (reduced to the submodular minimization problem), however the partition function calculation is not [START_REF] Kolmogorov | Convergent tree-reweighted message passing for energy minimization[END_REF]).

Examples of Log-supermodular Distributions

One of the examples to discuss is the binary pairwise MRF, which is also known as the Ising model and can be defined via the log-supermodular distribution, if the potentials satisfy the property of submodularity (so called attractive potentials that encourage the same signal sign among pairwise neighbors). Binary models can be naturally applied for computer vision purposes, e.g., for an image segmentation problem.

It is also possible to extend these models in two ways: to consider not only binary values, but discrete and even continuous (part of these work is focused on this), and also consider higher-order potentials instead of pairwise ones. The submodular machinery allow us to perform these extensions naturally. Thus, such models can be incorporated for example as sparsity-inducing priors or applied for semantic segmentation (see Chapter 4).

Parameter Learning and Inference

In this section we discuss a learning procedure for probabilistic graphical models and for the log-supermodular models in particular. There are plenty of approaches for parameter learning. We will cover the most prominent ones and will see that the learning procedure is usually computationally expensive.

Maximum Likelihood Estimation

We again consider a probabilistic model written as

𝑝(𝑥|𝜃) = 1 𝑍(𝜃) exp(𝜃 𝑇 𝜓(𝑥)).
Given the datapoints 𝑥 1 , . . . 𝑥 𝑁 , the averaged log-likelihood function ℓ(𝜃) has the form:

ℓ(𝜃) = 1 𝑁 𝑁 ∑ 𝑛=1 log 𝑝(𝑥 𝑛 |𝜃) = 𝜃 𝑇 ( 1 𝑁 𝑁 ∑ 𝑛=1 𝜓(𝑥 𝑛 ) ) -log 𝑍(𝜃) = 𝜃 𝑇 ψ(𝑥) -log 𝑍(𝜃),
where ψ(𝑥) is an empirical value of the expectation E𝜓(𝑥).

Optimization

Since we work with exponential families only, the given log-likelihood is concave and thereby has one global optimum that can be achieved by an approach from gradient ascent family. A gradient step would be possible since we evaluate the log-likelihood derivative:

𝜕ℓ(𝜃) 𝜕𝜃 = ψ(𝑥) - 𝜕 log 𝑍(𝜃) 𝜕𝜃 = ψ(𝑥) -E𝜓(𝑥|𝜃),
where we use the exponential family property of the log partition function derivative. We can notice that to make a gradient step, we have to perform an inference and evaluate the expectation E𝜓(𝑥|𝜃), which is usually computational expensive and can be done in several ways we discussed in the probabilistic inference section. There is also a moment matching effect, that we can observe while converging to the optimum. For the optimal parameter 𝜃 the expectation E𝜓(𝑥|𝜃) is going to match the empirical expectation ψ(𝑥).

Missing Data Treatment

With the tools given by the probabilistic approach we can also tackle the missing data case:

𝑝(𝑥, 𝑧|𝜃) = 1 𝑍(𝜃) exp(𝜃 𝑇 𝜓(𝑥, 𝑧)),
where we denote 𝑥 as observed variables and 𝑧 as unobserved or hidden. In the given circumstances, we will fit the optimal parameter 𝜃 to the observed data 𝑥 1 , . . . , 𝑥 𝑁 by solving an optimization problem in the form:

max 𝜃∈Θ ℓ(𝜃) = max 𝜃∈Θ 1 𝑁 𝑁 ∑ 𝑛=1 log ∑ 𝑧∈Z 𝑝(𝑥 𝑖 , 𝑧|𝜃) = max 𝜃∈Θ 1 𝑁 𝑁 ∑ 𝑛=1 log ∑ 𝑧∈Z 𝜃 𝑇 𝜓(𝑥 𝑖 , 𝑧) -log 𝑍(𝜃).
We should notice, that the summation term (where we marginalize out the hidden values 𝑧) is of similar complexity as the partition function 𝑍(𝜃) and thus provides another challenge to this optimization problem.

Conditional Maximum Likelihood

In this section we discuss the maximum likelihood formulation for the conditional random field model:

𝑝(𝑦|𝑥, 𝜃) = 1 𝑍(𝑥, 𝜃) exp(𝜃 𝑇 𝜓(𝑥, 𝑦)),
where 𝜓(𝑥, 𝑦) is a potential vector of feature vector 𝑥 and the label vector is 𝑦. This is the model of our particular interest in this thesis and we are going to stick with it in most of the times.

The new log-likelihood optimization problem can be written in the following form:

max 𝜃 1 𝑁 𝑁 ∑ 𝑛=1 log 𝑝(𝑦 𝑛 |𝑥 𝑛 , 𝜃) -𝜆||𝜃|| 2 2 ,
where the last term is a regularization term and can be seen as a log prior 𝑝(𝜃) on the parameters 𝜃. In this common formulation the Gaussian prior was used, however there are plenty of priors that can be incorporated, such as Laplacian, mixed priors or even parameterized customized prior (see Chapter 4).

Partition Function Approximation

As we seen in the previous section, the partition function term restricts the tractability of the parameter estimation via maximum likelihood. In this section we review various partition function calculation techniques. When the challenge of exact partition function computation becomes impossible, we resort to approximation techniques. Some of the approximations are upper bounds.

L-Field bound

This approximation was designed by [START_REF] Djolonga | From MAP to Marginals: Variational Inference in Bayesian Submodular Models[END_REF] directly for the log-supermodular distributions:

𝑝(𝐴) = exp(-𝐹 (𝐴)) 𝑍(-𝐹 ) ,
where 𝐹 (𝐴) is a submodular function. The approach is named L-field as there is a natural relation to mean field method discussed earlier and it is based on the Lovász extension (explains the L in the approach's name).

The main idea is to benefit from the submodular function property, that it can be bounded from the both sides by modular functions 𝑙, 𝑢 : 2 𝑉 → R.

Theorem 1.2. If we consider bounds 𝑙(𝐴) and 𝑢(𝐴) of the submodular function 𝐹 (𝐴):

𝑙(𝐴) ≤ 𝐹 (𝐴) ≤ 𝑢(𝐴) holds for any 𝐴 ⊆ 𝑉 , then the following is correct:

log 𝑍(-𝑢) ≤ log 𝑍(-𝐹 ) ≤ log 𝑍(-𝑙), (1.4.1)
where log 𝑍(-𝐹 ) = log

log 𝑍(𝑠) = log ∑ 𝑥∈{0,1} |𝑉 | exp(-𝑠 𝑇 𝑥) = log ∑ 𝑥∈{0,1} |𝑉 | |𝑉 | ∏ 𝑖=1 exp(-𝑠 𝑖 𝑥 𝑖 ) = = log |𝑉 | ∏ 𝑖=1 ∑ 𝑥 𝑖 ∈{0,1} exp(-𝑠 𝑖 𝑥 𝑖 ) = |𝑉 | ∑ 𝑖=1 log(1 + 𝑒 𝑠 𝑖 ).
Equation (1.4.1) becomes obvious if we imagine the summation and log term before the bound inequalities 𝑙(𝐴) ≤ 𝐹 (𝐴) ≤ 𝑢(𝐴). Via Theorem 1.2 we can bound the partition function from the both sides using some decomposable distributions based on modular functions. But there is no guarantee on the bounds tightness. Our motivation here would be to find the tightest possible modular function, to make our approximation the most accurate. Thus, we got a variational approach, that can be used to approximate the partition function from above:

log 𝑍(-𝐹 ) ≤ min 𝑙 log 𝑍(-𝑙) = min 𝑙 |𝑉 | ∑ 𝑖=1 log(1 + 𝑒 -𝑙 𝑖 ) subject to 𝑙(𝐴) ≤ 𝐹 (𝐴). Theorem 1.3. The optimization problem min 𝑙 |𝑉 | ∑ 𝑖=1 log(1 + 𝑒 -𝑙 𝑖 ) subject to 𝑙(𝐴) ≤ 𝐹 (𝐴)
is equivalent to the following one:

min 𝑙∈𝐵(𝐹 ) |𝑉 | ∑ 𝑖=1 log(1 + 𝑒 -𝑙 𝑖 ),
where 𝐵(𝐹 ) is a base polytope of the submodular function 𝐹 (𝐴).

This optimization problem is an existed problem of convex function optimization over the polytope, and can be solved via Frank-Wolfe algorithm [START_REF] Lacoste-Julien | Block-coordinate Frank-Wolfe optimization for structural SVMs[END_REF]). Finally, we can write down the L-Field upper bound of the partition function in the following way

𝐴 L-field = min 𝑙∈𝐵(𝐹 ) |𝑉 | ∑ 𝑖=1 log(1 + 𝑒 -𝑙 𝑖 ) ≥ log 𝑍(-𝐹 ).

Gumbel bound

In this section we work with a general Gibbs distribution supported on discrete (not binary) set

X = 𝑋 1 × • • • × 𝑋 𝑛 in the form: 𝑝(𝑥) = exp(𝑓 (𝑥)) 𝑍(𝑓 ) .
The Gumbel approximation assumes that MAP (maximum a posteriori) inference problem is feasible, i.e., there is a solver for the following problem:

arg max 𝑥∈X 𝑝(𝑥) = arg max 𝑥∈X 𝑓 (𝑥),
as the partition function is a constant and does not influence the solution. Thus, this makes the MAP inference problem much easier than the partition function calculation.

While the partition function evaluation is #𝑃 -hard problem [START_REF] Jerrum | Polynomial-time approximation algorithms for the Ising model[END_REF]), MAP is NP-hard in general, however it can be solved efficiently for many practical cases, e.g., for log-supermodular models.

Let us introduce the notation of Gumbel distribution as the approach we are going to present is based on it.

Definition 5. We will call a distribution a Gumbel distribution, if it has the following cumulative function

𝐹 (𝑡) = 𝑒 -𝑒 -(𝑡+𝑐) ,
where 𝑐 ≈ 0.577 is the Euler constant. The Gumbel distribution has a zero mean and a mode at the negative Euler constant. Its cumulative distribution and probability density functions are presented in Figure 1-7. The idea of the approximation is to use Gumbel-distributed perturbations to substitute the partition function calculation with a bunch of related MAP inference problems [START_REF] Hammersley | On the partition function and random maximum aposteriori perturbations[END_REF]).

Theorem 1.4. Let {𝛾(𝑥)} 𝑥∈X be a collection of i.i.d. random variables from the Gumbel distribution. Then the random variable max 𝑥∈X {𝑓 (𝑥) + 𝛾(𝑥)} belongs to the shifted Gumbel distribution, i.e., represents as Gumbel variable plus a constant, and the partition function can be represented as:

log 𝑍 = E 𝛾 [max 𝑥∈X {𝑓 (𝑥) + 𝛾(𝑥)}].
Proof. By definition we know that

𝐹 (𝑡) = 𝑃 (𝛾(𝑦) ≤ 𝑡).
Thus, we have, by independence,

∏ 𝑥∈X 𝐹 (𝑡 -𝑓 (𝑥)) = 𝑃 (max 𝑥∈X {𝑓 (𝑥) + 𝛾(𝑥)} ≤ 𝑡).
Since the shifted Gumbel cumulative distribution function is closed under multiplication, the following term is a cumulative distribution function as well:

∏ 𝑥∈X 𝐹 (𝑡-𝑓 (𝑥)) = exp ( - ∑ 𝑥∈X exp(-(𝑡-𝑓 (𝑥)+𝑐)) ) = exp(-exp(𝑡+𝑐)𝑍) = 𝐹 (𝑡-log 𝑍).
Expected value of the random variable max

𝑥∈X {𝑓 (𝑥) + 𝛾(𝑥)} consequently equals to log 𝑍. In this way log 𝑍 = E 𝛾 [max 𝑥∈X {𝑓 (𝑥) + 𝛾(𝑥)}].
However, this theorem just proposes another appearance of the partition function term. This is still a #𝑃 hard problem, since we have exponential many Gumbel perturbations {𝛾(𝑥)} 𝑥∈X . The next step would be to propose an upper bound with lower-dimensional perturbations.

Theorem 1.5. Using the setup from the Theorem 1.4, we can rewrite and bound the partition function in the following way:

log 𝑍 = E 𝛾 1 max 𝑥 1 . . . E 𝛾𝑛 max 𝑥𝑛 {𝑓 (𝑥) + 𝑛 ∑ 𝑖=1 𝛾 𝑖 (𝑥 𝑖 )} ≤ E 𝛾 max 𝑥∈X {𝑓 (𝑥) + 𝑛 ∑ 𝑖=1 𝛾 𝑖 (𝑥 𝑖 )}
See [START_REF] Hammersley | On the partition function and random maximum aposteriori perturbations[END_REF] for the detailed proof. The final upper bound and approximation has the following form:

𝐴 Gumbel = E 𝛾 max 𝑥∈X {𝑓 (𝑥) + 𝑛 ∑ 𝑖=1 𝛾 𝑖 (𝑥 𝑖 )}.
Thus, if we can perform MAP inference max 𝑥∈X {𝑓 (𝑥)}, then it is typically feasible to calculate MAP for single node perturbations max 𝑥∈X {𝑓 (𝑥) + ∑ 𝑛 𝑖=1 𝛾 𝑖 (𝑥 𝑖 )} as well, and also we can evaluate its expected value with Monte Carlo techniques.

Other partition function approximations can be constructed as well. The Mean Field approach gives a bound on the partition function and this bound can be used inside a learning algorithm. "Loopy" belief propagation algorithm also proposes an approximate treatment for the problem of inference and thus, provides a partition function approximation as a subroutine via message-passing technique. Other variations of this idea are Fractional Belief Propagation [START_REF] Wiegerinck | Fractional belief propagation[END_REF]), Tree-Reweighted Belief Propagation (Wainwright and Jordan [2008a]) and Generalized Belief Propagation [START_REF] Yedidia | Constructing free-energy approximations and generalized belief propagation algorithms[END_REF]). 19

Introduction

Submodular functions provide efficient and flexible tools for learning on discrete data. Several common combinatorial optimization tasks, such as clustering, image segmentation, or document summarization, can be achieved by the minimization or the maximization of a submodular function [START_REF] Bach | Learning with submodular functions: a convex optimization perspective[END_REF], [START_REF] Golovin | Adaptive Submodularity: Theory and Applications in Active Learning and Stochastic Optimization[END_REF], [START_REF] Lin | A class of submodular functions for document summarization[END_REF]). The key benefit of submodularity is the ability to model notions of diminishing returns, and the availability of exact minimization algorithms and approximate maximization algorithms with precise approximation guarantees [START_REF] Krause | Submodular function maximization[END_REF]).

In practice, it is not always straightforward to define an appropriate submodular function for a problem at hand. Given fully-labeled data, e.g., images and their foreground/background segmentations in image segmentation, structured-output prediction methods such as the structured-SVM may be used [START_REF] Szummer | Learning CRFs using graph cuts[END_REF]). However, it is common (a) to have missing data, and (b) to embed submodular function minimization within a larger model. These are two situations well tackled by probabilistic modelling.

Log-supermodular models, with negative log-densities equal to a submodular function, are a first important step toward probabilistic modelling on discrete data with submodular functions [START_REF] Djolonga | From MAP to Marginals: Variational Inference in Bayesian Submodular Models[END_REF]). However, it is well known that the log-partition function is intractable in such models. Several bounds have been proposed, that are accompanied with variational approximate inference [START_REF] Djolonga | Scalable Variational Inference in Log-supermodular Models[END_REF]). These bounds are based on the submodularity of the negative log-densities. However, parameter learning (typically by maximum likelihood), which is a key feature of probabilistic modeling, has not been tackled yet.

Contributions

In this chapter we highlight the following contributions:

-In Section 2.4, we review existing variational bounds for the log-partition function and show that the bound of [START_REF] Hammersley | On the partition function and random maximum aposteriori perturbations[END_REF], based on "perturb-and-MAP" ideas, formally dominates the bounds proposed by Djolonga andKrause [2014, 2015].

-In Section 2.5.1, we show that for parameter learning via maximum likelihood the existing bound of Djolonga andKrause [2014, 2015] typically leads to a degenerate solution while the one based on "perturb-and-MAP" ideas and logistic samples of [START_REF] Hammersley | On the partition function and random maximum aposteriori perturbations[END_REF] does not.

-In Section 2.5.2, given that the bound based on "perturb-and-MAP" ideas is an expectation (over our own randomization), we propose to use a stochastic subgradient technique to maximize the lower-bound on the log-likelihood, which can also be extended to conditional maximum likelihood.

-In Section 2.6, we illustrate our new results on a set of experiments in binary image denoising, where we highlight the flexibility of a probabilistic model for learning with missing data.

Submodular functions and log-supermodular models

In this section, we review the relevant theory of submodular functions and recall typical examples of log-supermodular distributions.

Submodular functions

We consider submodular functions on the vertices of the hypercube {0, 1} 𝐷 . This hypercube representation is equivalent to the power set of {1, . . . , 𝐷}. Indeed, we can go from a vertex of the hypercube to a set by looking at the indices of the components equal to one and from set to vertex by taking the indicator vector of the set.

For any two vertices of the hypercube, 𝑥, 𝑦 ∈ {0, 1} 𝐷 , a function 𝑓 : {0, 1} 𝐷 → R is submodular if

𝑓 (𝑥) + 𝑓 (𝑦) 𝑓 (min{𝑥, 𝑦}) + 𝑓 (max{𝑥, 𝑦}),
where the min and max operations are taken component-wise and correspond to the intersection and union of the associated sets. Equivalently, the function 𝑥 ↦ → 𝑓 (𝑥 + 𝑒 𝑖 ) -𝑓 (𝑥), where 𝑒 𝑖 ∈ R 𝐷 is the 𝑖-th canonical basis vector, is non-increasing. Hence, the notion of diminishing returns is often associated with submodular functions. Most widely used submodular functions are cuts, concave functions of subset cardinality, mutual information, set covers, and certain functions of eigenvalues of submatrices [START_REF] Bach | Learning with submodular functions: a convex optimization perspective[END_REF], [START_REF] Fujishige | Submodular Functions and Optimization[END_REF]). Supermodular functions are simply negatives of submodular functions.

In this chapter, we are going to use a few properties of such submodular functions (see [START_REF] Bach | Learning with submodular functions: a convex optimization perspective[END_REF], [START_REF] Fujishige | Submodular Functions and Optimization[END_REF] and references therein). Any submodular function 𝑓 can be extended from {0, 1} 𝐷 to a convex function on R 𝐷 , which is called the Lovász extension. This extension has the same value on {0, 1} 𝐷 , hence we use the same notation 𝑓 . Moreover, this function is convex and piecewise linear, which implies the existence of a polytope 𝐵(𝑓 ) ⊂ R 𝐷 , called the base polytope, such that for all 𝑥 ∈ R 𝐷 , 𝑓 (𝑥) = max

𝑠∈𝐵(𝑓 ) 𝑥 ⊤ 𝑠, that is, 𝑓 is the support function of 𝐵(𝑓 ).
The Lovász extension 𝑓 and the base polytope 𝐵(𝑓 ) have explicit expressions that are presented in the Chapter 1. We will only use the fact that 𝑓 can be efficiently minimized on {0, 1} 𝐷 , by a variety of generic algorithms, or by more efficient dedicated ones for subclasses such as graph-cuts.

Log-supermodular distributions

Log-supermodular models are introduced in Djolonga and Krause [2014] to model probability distributions on a hypercube, 𝑥 ∈ {0, 1} 𝐷 , and are defined as

𝑝(𝑥) = 1 𝑍(𝑓 ) exp(-𝑓 (𝑥)),
where 𝑓 : {0, 1} 𝐷 → R is a submodular function such that 𝑓 (0) = 0 and the partition function is

𝑍(𝑓 ) = ∑ 𝑥∈{0,1} 𝐷 exp(-𝑓 (𝑥))
. It is more convenient to deal with the convex log-partition function

𝐴(𝑓 ) = log 𝑍(𝑓 ) = log ∑ 𝑥∈{0,1} 𝐷 exp(-𝑓 (𝑥)).
In general, the calculation of the partition function 𝑍(𝑓 ) or the log-partition function 𝐴(𝑓 ) is intractable, as it includes simple binary Markov random fields-the exact calculation is known to be #𝑃 -hard [START_REF] Jerrum | Polynomial-time approximation algorithms for the Ising model[END_REF]). In Section 2.4, we review upper-bounds for the log-partition function.

Examples

Essentially, all submodular functions used in the minimization context can be used as negative log-densities (Djolonga andKrause [2014, 2015]). In computer vision, the most common examples are graph-cuts, which are essentially binary Markov random fields with attractive potentials, but higher-order potentials have been considered as well [START_REF] Kohli | Robust higher order potentials for enforcing label consistency[END_REF]). In our experiments, we use graph-cuts, where submodular function minimization may be performed with max-flow techniques and is thus efficient [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF]). Note that there are extensions of submodular functions to continuous domains that could be considered as well [START_REF] Bach | Submodular functions: from discrete to continuous domains[END_REF]).

Upper-bounds on the log-partition function

In this section, we review the main existing upper-bounds on the log-partition function for log-supermodular densities. These upper-bounds use several properties of submodular functions, in particular, the Lovász extension and the base polytope. Note that lower bounds based on submodular maximization aspects and superdifferentials [START_REF] Djolonga | From MAP to Marginals: Variational Inference in Bayesian Submodular Models[END_REF]) can be used to highlight the tightness of various bounds, which we present in Figure 2-1.

Base polytope relaxation with L-Field (Djolonga and Krause [2014])

This method exploits the fact that any submodular function 𝑓 (𝑥) can be lower bounded by a modular function 𝑠(𝑥), i.e., a linear function of 𝑥 ∈ {0, 1} 𝐷 in the hypercube representation. The submodular function and its lower bound are related by 𝑓 (𝑥) = max 𝑠∈𝐵(𝑓 ) 𝑠 ⊤ 𝑥, leading to:

𝐴(𝑓 ) = log ∑ 𝑥∈{0,1} 𝐷 exp (-𝑓 (𝑥)) = log ∑ 𝑥∈{0,1} 𝐷 min 𝑠∈𝐵(𝑓 ) exp (-𝑠 ⊤ 𝑥),
which, by swapping the sum and min, is less than min

𝑠∈𝐵(𝑓 ) log ∑ 𝑥∈{0,1} 𝐷 exp (-𝑠 ⊤ 𝑥) = min 𝑠∈𝐵(𝑓 ) 𝐷 ∑ 𝑑=1 log (1 + 𝑒 -𝑠 𝑑 ) def = 𝐴 L-field (𝑓 ).
(2.4.1)

Since the polytope 𝐵(𝑓 ) is tractable (through its membership oracle or by maximizing linear functions efficiently), the bound 𝐴 L-field (𝑓 ) is tractable, i.e., computable in polynomial time. Moreover, it has a nice interpretation through convex duality as the logistic function log(1 + 𝑒 -𝑠 𝑑 ) may be represented as max

𝜇 𝑑 ∈[0,1] -𝜇 𝑑 𝑠 𝑑 -𝜇 𝑑 log 𝜇 𝑑 -(1 - 𝜇 𝑑 ) log(1 -𝜇 𝑑 ), leading to: 𝐴 L-field (𝑓 ) = min 𝑠∈𝐵(𝑓 ) max 𝜇∈[0,1] 𝐷 -𝜇 ⊤ 𝑠 + 𝐻(𝜇) = max 𝜇∈[0,1] 𝐷 𝐻(𝜇) -𝑓 (𝜇),
where 𝐻(𝜇) = -

∑ 𝐷 𝑑=1 { 𝜇 𝑑 log 𝜇 𝑑 + (1 -𝜇 𝑑 ) log(1 -𝜇 𝑑 ) }
. This shows in particular the convexity of 𝑓 ↦ → 𝐴 L-field (𝑓 ). Finally, [START_REF] Djolonga | Scalable Variational Inference in Log-supermodular Models[END_REF] shows the remarkable result that the minimizer 𝑠 ∈ 𝐵(𝑓 ) may be obtained by minimizing a simpler function on 𝐵(𝑓 ), namely the squared Euclidean norm, thus leading to algorithms such as the minimum-norm-point algorithm [START_REF] Fujishige | Submodular Functions and Optimization[END_REF]).

"Pertub-and-MAP" with logistic distributions

Estimating the log-partition function can be done through optimization using "pertub-and-MAP" ideas. The main idea is to perturb the log-density, find the maximum a-posteriori configuration (i.e., perform optimization), and then average over several random perturbations [START_REF] Hammersley | On the partition function and random maximum aposteriori perturbations[END_REF], [START_REF] Papandreou | Perturb-and-MAP random fields: Using discrete optimization to learn and sample from energy models[END_REF], [START_REF] Tarlow | Randomized optimum models for structured prediction[END_REF]).

The Gumbel distribution on R, whose cumulative distribution function is 𝐹 (𝑧) = exp(-exp(-(𝑧 + 𝑐))), where 𝑐 is the Euler constant, is particularly useful. Indeed, if {𝑔(𝑦)} 𝑦∈{0,1} 𝐷 is a collection of independent random variables 𝑔(𝑦) indexed by 𝑦 ∈ {0, 1} 𝐷 , each following the Gumbel distribution, then the random variable max 𝑦∈{0,1} 𝐷 𝑔(𝑦) -𝑓 (𝑦) is such that we have

log 𝑍(𝑓 ) = E 𝑔 [ max 𝑦∈{0,1} 𝐷 {𝑔(𝑦) -𝑓 (𝑦)} ]
, from the [START_REF] Hammersley | On the partition function and random maximum aposteriori perturbations[END_REF], Lemma 1]. The main problem is that we need 2 𝐷 such variables, and a key contribution of [START_REF] Hammersley | On the partition function and random maximum aposteriori perturbations[END_REF] is to show that if we consider a factored collection {𝑔 𝑑 (𝑦 𝑑 )} 𝑦 𝑑 ∈{0,1},𝑑=1,...,𝐷 of i.i.d. Gumbel variables, then we get an upper-bound on the log partition-function, that is,

log 𝑍(𝑓 ) ≤ E 𝑔 max 𝑦∈{0,1} 𝐷 { 𝐷 ∑ 𝑑=1 𝑔 𝑑 (𝑦 𝑑 ) -𝑓 (𝑦)}.
Writing 𝑔 𝑑 (𝑦 𝑑 ) = [𝑔 𝑑 (1) -𝑔 𝑑 (0)]𝑦 𝑑 + 𝑔 𝑑 (0) and using the fact that (a) 𝑔 𝑑 (0) has zero expectation and (b) the difference between two independent Gumbel distributions has a logistic distribution (with cumulative distribution function 𝑧 ↦ → (1+𝑒 -𝑧 ) -1 ) [START_REF] Nadarajah | A generalized logistic distribution[END_REF]), we get the following upper-bound:

𝐴 Logistic (𝑓 ) = E 𝑧 1 ,...,𝑧 𝐷 ∼logistic [ max 𝑦∈{0,1} 𝐷 {𝑧 ⊤ 𝑦 -𝑓 (𝑦)} ] , (2.4.2)
where the random vector 𝑧 ∈ R 𝐷 consists of independent elements taken from the logistic distribution. This is always an upper-bound on 𝐴(𝑓 ) and it uses only the fact that submodular functions are efficient to optimize. It is convex in 𝑓 as an expectation of a maximum of affine functions of 𝑓 .

Comparison of bounds

In this section, we show that 𝐴 L-field (𝑓 ) is always dominated by 𝐴 Logistic (𝑓 ). This is complemented by another result within the maximum likelihood framework in Section 2.5. Proposition 2.1. For any submodular function 𝑓 : {0, 1} 𝐷 → R, we have:

𝐴(𝑓 ) 𝐴 Logistic (𝑓 ) 𝐴 L-field (𝑓 ).
(2.4.3)

Proof. The first inequality was shown by [START_REF] Hammersley | On the partition function and random maximum aposteriori perturbations[END_REF]. For the second inequality, we have:

𝐴 Logistic (𝑓 ) = E 𝑧 [ max 𝑦∈{0,1} 𝐷 𝑧 ⊤ 𝑦 -𝑓 (𝑦) ] = E 𝑧 [ max 𝑦∈{0,1} 𝐷 𝑧 ⊤ 𝑦 -max 𝑠∈𝐵(𝑓 ) 𝑠 ⊤ 𝑦 ] from properties of the base polytope 𝐵(𝑓 ), = E 𝑧 [ max 𝑦∈{0,1} 𝐷 min 𝑠∈𝐵(𝑓 ) 𝑧 ⊤ 𝑦 -𝑠 ⊤ 𝑦 ] , = E 𝑧 [ min 𝑠∈𝐵(𝑓 ) max 𝑦∈{0,1} 𝐷 𝑧 ⊤ 𝑦 -𝑠 ⊤ 𝑦 ] by convex duality, min 𝑠∈𝐵(𝑓 ) E 𝑧 [ max 𝑦∈{0,1} 𝐷 (𝑧 -𝑠) ⊤ 𝑦 ]
by swapping expectation and minimization, = min

𝑠∈𝐵(𝑓 ) E 𝑧 [ 𝐷 ∑ 𝑑=1 (𝑧 𝑑 -𝑠 𝑑 ) + ] by explicit maximization, = min 𝑠∈𝐵(𝑓 ) [ 𝐷 ∑ 𝑑=1 E 𝑧 𝑑 (𝑧 𝑑 -𝑠 𝑑 ) + ]
by using linearity of expectation, = min

𝑠∈𝐵(𝑓 ) [ 𝐷 ∑ 𝑑=1 ∫ +∞ -∞ (𝑧 𝑑 -𝑠 𝑑 ) + 𝑃 (𝑧 𝑑 )𝑑𝑧 𝑑 ] by definition of expectation, = min 𝑠∈𝐵(𝑓 ) [ 𝐷 ∑ 𝑑=1 ∫ +∞ 𝑠 𝑑 (𝑧 𝑑 -𝑠 𝑑 ) 𝑒 -𝑧 𝑑 (1 + 𝑒 -𝑧 𝑑 ) 2 𝑑𝑧 𝑑 ]
by substituting the density function, = min

𝑠∈𝐵(𝑓 ) 𝐷 ∑ 𝑑=1 log(1 + 𝑒 -𝑠 𝑑 )
, which leads to the desired result.

In the inequality above, since the logistic distribution has full support, there cannot be equality. However, if the base polytope is such that, with high probability ∀𝑑, |𝑠 𝑑 | ≥ |𝑧 𝑑 |, then the two bounds are close. Since the logistic distribution is concentrated around zero, we have equality when |𝑠 𝑑 | is large for all 𝑑 and 𝑠 ∈ 𝐵(𝑓 ).

Running-time complexity of 𝐴 L-field and 𝐴 logistic . The logistic bound 𝐴 logistic can be computed if there is efficient MAP-solver for submodular functions (plus a modular term). In this case, the divide-and-conquer algorithm can be applied for L-Field [START_REF] Djolonga | From MAP to Marginals: Variational Inference in Bayesian Submodular Models[END_REF]). Thus, the complexity is dedicated to the minimization of 𝑂(|𝑉 |) problems. Meanwhile, for the method based on logistic samples, it is necessary to solve 𝑀 optimization problems. In our empirical bound comparison (next paragraph), the running time was the same for both methods. Note however that for parameter learning, we need a single SFM problem per gradient iteration (and not 𝑀 ).

Empirical comparison of 𝐴 L-field and 𝐴 logistic . We compare the upper-bounds on the log-partition function 𝐴 L-field and 𝐴 logistic , with the setup used by [START_REF] Djolonga | From MAP to Marginals: Variational Inference in Bayesian Submodular Models[END_REF]. We thus consider data from a Gaussian mixture model with 2 clusters in R 2 . The centers are sampled from N( [3,3], 𝐼) and N([-3, -3], 𝐼), respectively. Then we sampled 𝑛 = 50 points for each cluster. Further, these 2𝑛 points are used as nodes in a complete weighted graph, where the weight between points 𝑥 and 𝑦 is equal to 𝑒 -𝑐||𝑥-𝑦|| .

We consider the graph cut function associated to this weighted graph, which defines a log-supermodular distribution. We then consider conditional distributions, one for each 𝑘 = 1, . . . , 𝑛, on the events that at least 𝑘 points from the first cluster lie on the one side of the cut and at least 𝑘 points from the second cluster lie on the other side of the cut. For each conditional distribution, we evaluate and compare the two upper bounds. We also add the tree-reweighted belief propagation upper bound (Wainwright and Jordan [2008b]) and the superdifferential-based lower bound [START_REF] Djolonga | From MAP to Marginals: Variational Inference in Bayesian Submodular Models[END_REF]).

In Figure 2-1, we show various bounds on 𝐴(𝑓 ) as functions of the number on conditioned pairs. The logistic upper bound is obtained using 100 logistic samples: the logistic upper-bound 𝐴 logistic is close to the superdifferential lower bound from Djolonga and Krause [2014] and is indeed significantly lower than the bound 𝐴 L-field . However, the tree-reweighted belief propagation bound behaves a bit better in the second case, but its calculation takes more time, and it cannot be applied for general submodular functions.

From bounds to approximate inference

Since linear functions are submodular functions, given any convex upper-bound on the log-partition function, we may derive an approximate marginal probability for each 𝑥 𝑑 ∈ {0, 1}. Indeed, following [START_REF] Hammersley | On the partition function and random maximum aposteriori perturbations[END_REF], we consider an exponential family model 𝑝(𝑥|𝑡) = exp(-𝑓 (𝑥) + 𝑡 ⊤ 𝑥 -𝐴(𝑓 -𝑡)), where 𝑓 -𝑡 is the function 𝑥 ↦ → 𝑓 (𝑥) -𝑡 ⊤ 𝑥. When 𝑓 is assumed to be fixed, this can be seen as an exponential family with the base measure exp(-𝑓 (𝑥)), sufficient statistics 𝑥, and 𝐴(𝑓 -𝑡) is the log-partition function. It is known that the expectation of the sufficient statistics under the exponential family model E 𝑝(𝑥|𝑡) 𝑥 is the gradient of the log-partition function (Wainwright and Jordan [2008b]). Hence, any approximation of this log-partition gives an approximation of this expectation, which in our situation is the vector of marginal probabilities that an element is equal to 1. For the logistic bound, this is the inference mechanism from [START_REF] Hammersley | On the partition function and random maximum aposteriori perturbations[END_REF], with 𝜕 𝑡 𝑑 𝐴 logistic (𝑓 -𝑡) = E 𝑧 𝑦 * (𝑧), where 𝑦 * (𝑧) is the maximizer of max 𝑦∈{0,1} 𝐷 𝑧 ⊤ 𝑦 -𝑓 (𝑦). In practice, in order to perform approximate inference, we only sample 𝑀 logistic variables. We could do the same for parameter learning, but a much more efficient alternative, based on mixing sampling and convex optimization, is presented in the next section.

Parameter learning through maximum likelihood

An advantage of log-supermodular probabilistic models is the opportunity to learn the model parameters from data using the maximum-likelihood principle. In this section, we consider that we are given 𝑁 observations 𝑥 1 , . . . , 𝑥 𝑁 ∈ {0, 1} 𝐷 , e.g., binary images such as shown in Figure 2-2.

We consider a submodular function 𝑓 (𝑥) represented as

𝑓 (𝑥) = 𝐾 ∑ 𝑘=1 𝛼 𝑘 𝑓 𝑘 (𝑥) -𝑡 ⊤ 𝑥.
The modular term 𝑡 ⊤ 𝑥 is explicitly taken into account with 𝑡 ∈ R 𝐷 , and 𝐾 base submodular functions are assumed to be given with 𝛼 ∈ R 𝐾 + so that the function 𝑓 remains submodular. Assuming the data 𝑥 1 , . . . , 𝑥 𝑁 are independent and identically (i.i.d.) distributed, then maximum likelihood is equivalent to minimizing: min

𝛼∈R 𝐾 + , 𝑡∈R 𝐷 - 1 𝑁 𝑁 ∑ 𝑛=1 log 𝑝(𝑥 𝑛 |𝛼, 𝑡) = min 𝛼∈R 𝐾 + , 𝑡∈R 𝐷 1 𝑁 𝑁 ∑ 𝑛=1 { 𝐾 ∑ 𝑘=1 𝛼 𝑘 𝑓 𝑘 (𝑥 𝑛 ) -𝑡 ⊤ 𝑥 𝑛 + 𝐴(𝑓 ) } ,
which takes the particularly simple form min

𝛼∈R 𝐾 + , 𝑡∈R 𝐷 𝐾 ∑ 𝑘=1 𝛼 𝑘 ( 1 𝑁 𝑁 ∑ 𝑛=1 𝑓 𝑘 (𝑥 𝑛 ) ) -𝑡 ⊤ ( 1 𝑁 𝑁 ∑ 𝑛=1 𝑥 𝑛 )
+ 𝐴(𝛼, 𝑡), (2.5.1) where we use the notation 𝐴(𝛼, 𝑡) = 𝐴(𝑓 ). We now consider replacing the intractable log-partition function by its approximations defined in Section 2.4.

Learning with the L-field approximation

In this section, we show that if we replace 𝐴(𝑓 ) by 𝐴 L-field (𝑓 ), we obtain a degenerate solution. Indeed, we have

𝐴 L-field (𝛼, 𝑡) = min 𝑠∈𝐵(𝑓 ) 𝐷 ∑ 𝑑=1 log (1 + 𝑒 -𝑠 𝑑 ) = min 𝑠∈𝐵( ∑ 𝐾 𝑘=1 𝛼 𝐾 𝑓 𝐾 ) 𝐷 ∑ 𝑑=1 log (1 + 𝑒 -𝑠 𝑑 +𝑡 𝑑 ).
This implies that Eq. ( 2.5.1) becomes min

𝛼∈R 𝐾 + , 𝑡∈R 𝐷 min 𝑠∈𝐵( ∑ 𝐾 𝑘=1 𝛼 𝐾 𝑓 𝐾 ) 𝐾 ∑ 𝑘=1 𝛼 𝑘 ( 1 𝑁 𝑁 ∑ 𝑛=1 𝑓 𝑘 (𝑥 𝑛 ) ) -𝑡 ⊤ ( 1 𝑁 𝑁 ∑ 𝑛=1 𝑥 𝑛 ) + 𝐷 ∑ 𝑑=1 log (1 + 𝑒 -𝑠 𝑑 +𝑡 𝑑 ).
The minimum with respect to 𝑡 𝑑 may be performed in closed form with 𝑡 𝑑 -𝑠 𝑑 = log ⟨𝑥⟩ 𝑑 1-⟨𝑥⟩ 𝑑 , where ⟨𝑥⟩ =

∑ 𝐾 𝑘=1 𝛼 𝐾 𝑓 𝐾 ) 𝐾 ∑ 𝑘=1 𝛼 𝑘 ( 1 𝑁 𝑁 ∑ 𝑛=1 𝑓 𝑘 (𝑥 𝑛 ) ) -𝑠 ⊤ ( 1 𝑁 𝑁 ∑ 𝑛=1 𝑥 𝑛 ) + const ,
which is equivalent to, using the representation of 𝑓 as the support function of 𝐵(𝑓 ):

min 𝛼∈R 𝐾 + 𝐾 ∑ 𝑘=1 𝛼 𝑘 [ 1 𝑁 𝑁 ∑ 𝑛=1 𝑓 𝑘 (𝑥 𝑛 ) -𝑓 𝑘 ( 1 𝑁 𝑁 ∑ 𝑛=1 𝑥 𝑛 )]
.

Since 𝑓 𝑘 is convex, by Jensen's inequality, the linear term in 𝛼 𝑘 is non-negative; thus maximum likelihood through L-field will lead to a degenerate solution where all 𝛼's are equal to zero.

Learning with the logistic approximation with stochastic gradients

In this section we consider the problem (2.5.1) and replace 𝐴(𝑓 ) by 𝐴 Logistic (𝑓 ):

min

𝛼∈R 𝐾 + , 𝑡∈R 𝐷 𝐾 ∑ 𝑘=1 𝛼 𝑘 ⟨𝑓 𝑘 (𝑥)⟩ emp. -𝑡 ⊤ ⟨𝑥⟩ emp. + + E 𝑧∼logistic [ max 𝑦∈{0,1} 𝐷 𝑧 ⊤ 𝑦 + 𝑡 ⊤ 𝑦 - 𝐾 ∑ 𝑘=1 𝛼 𝑘 𝑓 (𝑦) ] , (2.5.2)
where ⟨𝑀 (𝑥)⟩ emp. denotes the empirical average of 𝑀 (𝑥) (over the data). Denoting by 𝑦 * (𝑧, 𝑡, 𝛼) ∈ {0, 1} 𝐷 the maximizers of 𝑧 ⊤ 𝑦 + 𝑡 ⊤ 𝑦 -∑ 𝐾 𝑘=1 𝛼 𝑘 𝑓 (𝑦), the objective function may be written:

𝐾 ∑ 𝑘=1 𝛼 𝑘 [ ⟨𝑓 𝑘 (𝑥)⟩ emp. -⟨𝑓 𝑘 (𝑦 * (𝑧, 𝑡, 𝛼))⟩ logistic ] - -𝑡 ⊤ [ ⟨𝑥⟩ emp. -⟨𝑦 * (𝑧, 𝑡, 𝛼)⟩ logistic ] + ⟨𝑧 ⊤ 𝑦 * (𝑧, 𝑡, 𝛼)⟩ logistic .
This implies that at optimum, for 𝛼 𝑘 > 0, then ⟨𝑓 𝑘 (𝑥)⟩ emp. = ⟨𝑓 𝑘 (𝑦 * (𝑧, 𝑡, 𝛼))⟩ logistic , while, ⟨𝑥⟩ emp. = ⟨𝑦 * (𝑧, 𝑡, 𝛼)⟩ logistic , the expected values of the sufficient statistics match between the data and the optimizers used for the logistic approximation [START_REF] Hammersley | On the partition function and random maximum aposteriori perturbations[END_REF]).

In order to minimize the expectation in Eq. ( 2.5.2), we propose to use the projected stochastic gradient method, not on the data as usually done, but on our own internal randomization. The algorithm then becomes, once we add a weighted ℓ 2 -regularization Ω(𝑡, 𝛼):

-Input: functions 𝑓 𝑘 , 𝑘 = 1, . . . , 𝐾, and expected sufficient statistics ⟨𝑓 𝑘 (𝑥)⟩ emp. ∈ R and ⟨𝑥⟩ emp. ∈ [0, 1] 𝐷 , regularizer Ω(𝑡, 𝛼).

-Initialization: 𝛼 = 0, 𝑡 = 0 -Iterations: for ℎ from 1 to 𝐻 -Sample 𝑧 ∈ R 𝐷 as independent logistics -Compute 𝑦 * = 𝑦 * (𝑧, 𝑡, 𝛼) ∈ arg max 𝑦∈{0,1} 𝐷 𝑧 ⊤ 𝑦 + 𝑡 ⊤ 𝑦 - ∑ 𝐾 𝑘=1 𝛼 𝑘 𝑓 (𝑦) -Replace 𝑡 by 𝑡 -𝐶 √ ℎ [ 𝑦 * -⟨𝑥⟩ emp. + 𝜕 𝑡 Ω(𝑡, 𝛼) ] -Replace 𝛼 𝑘 by ( 𝛼 𝑘 -𝐶 √ ℎ [ ⟨𝑓 𝑘 (𝑥)⟩ emp. -𝑓 𝑘 (𝑦 * ) + 𝜕 𝛼 𝑘 Ω(𝑡, 𝛼) ]) + .
-Output: (𝛼, 𝑡).

Since our cost function is convex and Lipschitz-continuous, the averaged iterates are converging to the global optimum [START_REF] Nemirovski | Robust stochastic approximation approach to stochastic programming[END_REF]) at rate 1/ √ 𝐻 (for function values).

Extension to conditional maximum likelihood

In experiments in Section 2.6, we consider a joint model over two binary vectors 𝑥, 𝑧 ∈ R 𝐷 , as follows

𝑝(𝑥, 𝑧|𝛼, 𝑡, 𝜋) = 𝑝(𝑥|𝛼, 𝑡)𝑝(𝑧|𝑥, 𝜋) = exp(-𝑓 (𝑥) -𝐴(𝑓 )) 𝐷 ∏ 𝑑=1 𝜋 𝛿(𝑧 𝑑 ̸ =𝑥 𝑑 ) 𝑑 (1 -𝜋 𝑑 ) 𝛿(𝑧 𝑑 =𝑥 𝑑 ) ,
(2.5.3) which corresponds to sampling 𝑥 from a log-supermodular model and considering 𝑧 that switches the values of 𝑥 with probability 𝜋 𝑑 for each 𝑑, that is, a noisy observation of 𝑥. We have:

log 𝑝(𝑥, 𝑧|𝛼, 𝑡, 𝜋) = -𝑓 (𝑥) -𝐴(𝑓 ) + 𝐷 ∑ 𝑑=1 { -log(1 + 𝑒 𝑢 𝑑 ) + 𝑥 𝑑 𝑢 𝑑 + 𝑧 𝑑 𝑢 𝑑 -2𝑥 𝑑 𝑧 𝑑 𝑢 𝑑 } , with 𝑢 𝑑 = log 𝜋 𝑑 1-𝜋 𝑑 which is equivalent to 𝜋 𝑑 = (1 + 𝑒 -𝑢 𝑑 ) -1 . Using Bayes rule, we have 𝑝(𝑥|𝑧, 𝛼, 𝑡, 𝜋) ∝ exp(-𝑓 (𝑥) -𝐴(𝑓 ) + 𝑥 ⊤ 𝑢 -2𝑥 ⊤ (𝑢 ∘ 𝑧)),
which leads to the log-supermodular model

𝑝(𝑥|𝑧, 𝛼, 𝑡, 𝜋) = exp(-𝑓 (𝑥) + 𝑥 ⊤ (𝑢 -2𝑢 ∘ 𝑧) -𝐴(𝑓 -𝑢 + 2𝑢 ∘ 𝑧)).
Thus, if we observe both 𝑧 and 𝑥, we can consider a conditional maximization of the log-likelihood (still a convex optimization problem), which we do in our experiments for supervised image denoising, where we assume we know both noisy and original images at training time. Stochastic gradient on the logistic samples can then be used. Note that our conditional ML estimation can be seen as a form of approximate conditional random fields [START_REF] Lafferty | Conditional random fields: Probabilistic models for segmenting and labeling sequence data[END_REF]).

While supervised learning can be achieved by other techniques such as structuredoutput-SVMs [START_REF] Szummer | Learning CRFs using graph cuts[END_REF], [START_REF] Taskar | Learning structured prediction models: A large margin approach[END_REF], [START_REF] Tsochantaridis | Large margin methods for structured and interdependent output variables[END_REF]), our approach also applies when we do not observe the original image, which we now consider.

Missing data through maximum likelihood

In the model in Eq. (2.5.3), we now assume we only observed the noisy output 𝑧, and we perform parameter learning for 𝛼, 𝑡, 𝜋. This is a latent variable model for which maximum likelihood can be readily applied. We have: In practice, we will assume that the noise probability 𝜋 (and hence 𝑢) is uniform across all elements. While we could use majorization-minization approaches such as the expectation-minimization algorithm (EM), we consider instead stochastic subgradient descent to learn the model parameters 𝛼, 𝑡 and 𝑢 (now a non-convex optimization problem, for which we still observed good convergence).

log 𝑝(𝑧|𝛼, 𝑡, 𝜋) = log ∑ 𝑥∈{0,1} 𝑝(𝑧, 𝑥|𝛼, 𝑡, 𝜋) = log ∑ 𝑥∈{0,1} 𝐷 exp(-𝑓 (𝑥) -𝐴(𝑓 )) 𝐷 ∏ 𝑑=1 𝜋 𝛿(𝑧 𝑑 ̸ =𝑥 𝑑 ) 𝑑 (1 -𝜋 𝑑 ) 𝛿(𝑧 𝑑 =𝑥 𝑑 )

Experiments

The aim of our experiments is to demonstrate the ability of our approach to remove noise in binary images, following the experimental set-up of [START_REF] Hammersley | On the partition function and random maximum aposteriori perturbations[END_REF]. We consider the training sample of 𝑁 𝑡𝑟𝑎𝑖𝑛 = 100 images of size 𝐷 = 50 × 50, and the test sample of 𝑁 𝑡𝑒𝑠𝑡 = 100 binary images, containing a horse silhouette from the Weizmann horse database [START_REF] Borenstein | Combining top-down and bottom-up segmentation[END_REF]). At first we add some noise by flipping pixels values independently with probability 𝜋. In Figure 2-2, we provide an example from the test sample: the original, the noisy and the denoised image (by our algorithm).

We consider the model from Section 2.5.3, with the two functions 𝑓 1 (𝑥), 𝑓 2 (𝑥) which are horizontal and vertical cut functions with binary weights respectively, together with a modular term of dimension 𝐷. To perform minimization we use graph-cuts [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF]) as we deal with positive or attractive potentials.

Supervised image denoising.

We assume that we observe 𝑁 = 100 pairs (𝑥 𝑖 , 𝑧 𝑖 ) of original-noisy images, 𝑖 = 1, . . . , 𝑁 . We perform parameter inference by maximum likelihood using stochastic subgradient descent (over the logistic samples), with regularization by the squared ℓ 2 -norm, one parameter for 𝑡, one for 𝛼, both learned by cross-validation. Given our estimates, we may denoise a new image by computing the "max-marginal", e.g., the maximum a posteriori max 𝑥 𝑝(𝑥|𝑧, 𝛼, 𝑡) through a single graph-cut, or computing "mean-marginals" with 100 logistic samples. To calculate the error we use the normalized Hamming distance and 100 test images. Results are presented in Table 2.1, where we compare the two types of decoding, as well as a structured output SVM (SVM-Struct [START_REF] Tsochantaridis | Large margin methods for structured and interdependent output variables[END_REF]) applied to the same problem. Results are reported in proportion of correct pixels. We see that the probabilistic models here slightly outperform the max-margin formulation1 and that using mean-marginals (which is optimal given our loss measure) lead to slightly better performance. Unsupervised image denoising. We now only consider 𝑁 = 100 noisy images 𝑧 1 , . . . , 𝑧 𝑁 to learn the model, without the original images, and we use the latent model from Section 2.5.4. We apply stochastic subgradient descent for the difference of the two convex functions 𝐴 logistic to learn the model parameters and use fixed regularization parameters equal to 10 -2 . We consider two situations, with a known noise-level 𝜋 or with learning it together with 𝛼 and 𝑡. The error was calculated using either max-marginals and mean-marginals.

𝜋

Note that here, structured-output SVMs cannot be used because there is no supervision. Results are reported in Tables 2.2 and 2.3. One explanation for a better performance for max-marginals in this case is that the unsupervised approach tends to oversmooth the outcome and max-marginals correct this a bit.

When the noise level is known, the performance compared to supervised learning is not degraded much, showing the ability of the probabilistic models to perform parameter estimation with missing data. When the noise level is unknown and learned as well, results are worse, still better than a trivial answer for moderate levels of noise (5% and 10%) but not better than outputting the noisy image for extreme levels (1% and 20%). In challenging fully unsupervised case the standard deviation is up to 2.2% (which shows that our results are statistically significant).

Conclusion

In this paper, we have presented how approximate inference based on stochastic gradient and "perturb-and-MAP" ideas could be used to learn parameters of log-supermodular models, allowing to benefit from the versatility of probabilistic modelling, in particular in terms of parameter estimation with missing data. While our experiments have focused on simple binary image denoising, exploring larger-scale applications in computer vision (such as done by [START_REF] Zhang | Higher-order inference for multi-class logsupermodular models[END_REF], [START_REF] Tschiatschek | Learning probabilistic submodular diversity models via noise contrastive estimation[END_REF]) should also show the benefits of mixing probabilistic modelling and submodular functions.

Chapter 3

Marginal Weighted Maximum Log-likelihood for Efficient Learning of Perturb-and-Map Models Abstract

In this part we consider the structured-output prediction problem through probabilistic approaches and generalize the "perturb-and-MAP" framework to more challenging weighted Hamming losses, which are crucial in applications. While in principle our approach is a straightforward marginalization, it requires solving many related MAP inference problems. We show that for log-supermodular pairwise models these operations can be performed efficiently using the machinery of dynamic graph cuts. We also propose to use double stochastic gradient descent, both on the data and on the perturbations, for efficient learning. Our framework can naturally take weak supervision (e.g., partial labels) into account. We conduct a set of experiments on medium-scale character recognition and image segmentation, showing the benefits of our algorithms. This chapter is based on the work "Marginal Weighted Maximum Log-likelihood for Efficient Learning of Perturb-and-Map Models", T. Shpakova, F. Bach, A. Osokin, published in Proceedings of the Conference on Uncertainty in Artificial Intelligence (UAI), 2018.

Introduction

Structured-output prediction is an important and challenging problem in the field of machine learning. When outputs have a structure, often in terms of parts or elements (e.g., pixels, sentences or characters), methods that do take it into account typically outperform more naive methods that consider outputs as a set of independent elements. Structured-output methods based on optimization can be broadly separated in two main families: max-margin methods, such as structured support vector machines (SSVM from [START_REF] Taskar | Learning structured prediction models: A large margin approach[END_REF], [START_REF] Tsochantaridis | Large margin methods for structured and interdependent output variables[END_REF]) and probabilistic methods based on maximum likelihoods such as conditional random fields (CRF) [START_REF] Lafferty | Conditional random fields: Probabilistic models for segmenting and labeling sequence data[END_REF]).

Structured-output prediction faces many challenges: (1) on top of large input dimensions, problems also have large outputs, leading to scalability issues, in particular when prediction or learning depends on combinatorial optimization problems (which are often polynomial-time, but still slow given they are run many times); ( 2) it is often necessary to use losses which go beyond the traditional 0-1loss to shape the behavior of the learned models towards the final evaluation metric; (3) having fully labelled data is either rare or expensive and thus, methods should be able to deal with weak supervision.

Max-margin methods can be used with predefined losses, and have been made scalable by several recent contributions (see, e.g., Lacoste-Julien et al. [2013] and references therein), but do not deal naturally with weak supervision. However, a few works [START_REF] Yu | Learning structural SVMs with latent variables[END_REF], [START_REF] Kumar | Self-paced learning for latent variable models[END_REF], [START_REF] Girshick | Object detection with grammar models[END_REF]) incorporate weak supervision into the max-margin approach via the CCCP [START_REF] Yuille | The concave-convex procedure[END_REF]) algorithm.

The flexibility of probabilistic modeling naturally allows (a) taking into consideration weak supervision and (b) characterizing the uncertainty of predictions, but it comes with strong computational challenges as well as a non-natural way of dealing with predefined losses beyond the 0-1 loss. The main goal of this chapter is to provide new tools for structured-output inference with probabilistic models, thus making them more widely applicable, while still being efficient. There are two main techniques to allow for scalable learning in CRFs: stochastic optimization [START_REF] Vishwanathan | Accelerated training of conditional random fields with stochastic gradient methods[END_REF]) and piecewise training (Sutton andMcCallum [2005, 2007], [START_REF] Kolesnikov | Closed-form training of conditional random fields for large scale image segmentation[END_REF]); note that the techniques above can also be used for weak supervision (and we reuse some of them in this work).

Learning and inference in probabilistic structured-output models recently received a lot of attention from the research community (e.g., [START_REF] Bakir | Predicting Structured Data[END_REF], [START_REF] Nowozin | Structured learning and prediction in computer vision[END_REF], [START_REF] Smith | Linguistic structure prediction[END_REF]). In this work we consider only models for which maximum-a posteriori (MAP) inference is feasible (a step often referred to as decoding in max-margin formulations, and which typically makes them tractable). A lot of efforts were spent to explore MAP-solvers algorithms for various problems, leveraging various structures, e.g., graphs of low tree-width [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF], Wainwright and Jordan [2008b], [START_REF] Sontag | Tightening LP relaxations for MAP using message-passing[END_REF], [START_REF] Komodakis | MRF energy minimization and beyond via dual decomposition[END_REF]) and function submodularity [START_REF] Boros | Pseudo-boolean optimization[END_REF], [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts[END_REF], [START_REF] Bach | Learning with submodular functions: a convex optimization perspective[END_REF], [START_REF] Osokin | Submodular relaxation for inference in Markov random fields[END_REF]).

Naturally, the existence of even an exact and efficient MAP-solver does not mean that the partition function (a key tool for probabilistic inference as shown below) is tractable to compute. Indeed, the partition function computation is known to be #𝑃 -hard [START_REF] Jerrum | Polynomial-time approximation algorithms for the Ising model[END_REF]) in general. For example, MAP-inference is efficient for log-supermodular probabilistic models, while computation of their partition function is not [START_REF] Djolonga | From MAP to Marginals: Variational Inference in Bayesian Submodular Models[END_REF]).

For such problems where MAP-inference is efficient, but partition function computation is not, "perturb-and-MAP" ideas such as proposed by [START_REF] Papandreou | Perturb-and-MAP random fields: Using discrete optimization to learn and sample from energy models[END_REF], [START_REF] Hammersley | On the partition function and random maximum aposteriori perturbations[END_REF] are a very suitable treatment. By adding random perturbations, and then performing MAP-inference, they can lead to estimates of the partition function. In Section 3.3, we review the existing approaches to the partition function approximation, parameter learning and inference.

An attempt to learn parameters via "perturb-and-MAP" ideas was made by [START_REF] Hazan | Learning efficient random maximum a-posteriori predictors with non-decomposable loss functions[END_REF], where the authors have developed a PAC-Bayesian-flavoured approach for the non-decomposable loss functions. While the presented algorithm has something in common with ours (gradient descent optimization of the objective upper bound), it differs in the sense of the objective function and the problem setup, which is more general but that requires a different (potentially with higher variance) estimates of the gradients. Such estimates are usual in reinforcement learning, e.g., the log-derivative trick from the REINFORCE algorithm [START_REF] Williams | Simple statistical gradient-following algorithms for connectionist reinforcement learning[END_REF]).

The goal of this chapter is to make the "perturb-and-MAP" technique applicable to practical problems, in terms of (a) scale, by increasing the problem size significantly, and (b) losses, by treating structured losses such as the Hamming loss or its weighted version, which are crucial to obtaining good performance in practice.

Contributions

Overall, in this chapter we make the following contributions:

-In Section 3.4, we generalize the "perturb-and-MAP" framework to more challenging weighted Hamming losses which are commonly used in applications. In principle, this is a straightforward marginalization but this requires solving many related MAP inference problems. We show that for graph cuts (our main inference algorithm for image segmentation), this can be done particularly efficiently. Besides that, we propose to use a double stochastic gradient descent, both on the data and on the perturbations.

-In Section 3.5, we show how weak supervision (e.g., partial labels) can be naturally dealt with. Our method in this case relies on approximating marginal probabilities that can be done almost at the cost of the partition function approximation.

-In Section 3.6, we conduct a set of experiments on medium-scale character recognition and image segmentation, showing the benefits of our new algorithms.

Perturb-and-MAP

In this section, we introduce the notation and review the relevant background. We study the following probabilistic model (a.k.a. a Gibbs distribution) over a discrete product space

𝑌 = 𝑌 1 × • • • × 𝑌 𝐷 , 𝑃 (𝑦) = 1 𝑍(𝑓 ) 𝑒 𝑓 (𝑦) , (3.3.1)
which is defined by a potential function 𝑓 : 𝑌 → R. The constant 𝑍(𝑓 ) = ∑ 𝑦∈𝑌 𝑒 𝑓 (𝑦) is called the partition function and normalizes 𝑃 (𝑦) to be a valid probability function, i.e., to sum to one. 𝑍(𝑓 ) is in general intractable to compute as the direct computation requires summing over exponentially (in 𝐷) many elements. Various partition function approximations methods have been used in parameter learning algorithms [START_REF] Parise | Learning in Markov random fields: an empirical study[END_REF]), e.g., mean-field [START_REF] Jordan | An introduction to variational methods for graphical models[END_REF]), tree-reweighted belief propagation (Wainwright and Jordan [2008b]) or loopy belief propagation [START_REF] Weiss | Comparing the mean field method and belief propagation for approximate inference in MRFs[END_REF]). We will work with the upper bound on the partition function proposed by [START_REF] Hammersley | On the partition function and random maximum aposteriori perturbations[END_REF] as it allows us to approximate the partition function via MAP-inference, calculate gradients efficiently, approximate marginal probabilities and guarantee tightness for some probabilistic models. We introduce this class of techniques below.

Gumbel perturbations

Recently, [START_REF] Hammersley | On the partition function and random maximum aposteriori perturbations[END_REF] provided a general-purpose upper bound on the log-partition function 𝐴(𝑓 ) = log 𝑍(𝑓 ), based on the "perturb-and-MAP" idea [START_REF] Papandreou | Perturb-and-MAP random fields: Using discrete optimization to learn and sample from energy models[END_REF]): maximize the potential function perturbed by Gumbel-distributed noise.1 Proposition 3.1 [START_REF] Hammersley | On the partition function and random maximum aposteriori perturbations[END_REF], Corollary 1). For any function 𝑓 : 𝑌 → R, we have 𝐴(𝑓 ) 𝐴 G (𝑓 ), where

𝐴 G (𝑓 ) = E 𝑧 1 ,...,𝑧 𝐷 ∼Gumbel [ max 𝑦∈𝑌 ( 𝑓 (𝑦) + 𝐷 ∑ 𝑑=1 𝑧 𝑑 (𝑦 𝑑 ) )] . (3.3.2)
Gumbel denotes the Gumbel distribution and {𝑧 𝑑 (𝑦 𝑑 )} 𝑑=1,...,𝐷 𝑦 𝑑 ∈𝑌 𝑑 is a collection of independent Gumbel samples.

The bound is tight when 𝑓 (𝑦) is a separable function (i.e., a sum of functions of single variables), and the tightness of this bound was further studied by [START_REF] Shpakova | Parameter learning for log-supermodular distributions[END_REF] for log-supermodular models (where 𝑓 is supermodular). They have shown that the bound 𝐴 𝐺 is always lower (and thus provide a better bound) than the "L-field" bound proposed by Djolonga andKrause [2014, 2015], which is itself based on separable optimization on the base polytope of the associated supermodular function.

The partition function bound 𝐴 𝐺 can be approximated by replacing the expectation by an empirical average. That is, to approximate it we need to solve a large number (as many as the number of Gumbel samples used to approximate the expectation) of MAP-like problems (i.e., maximizing 𝑓 plus a separable function) which are feasible by our assumption. Strictly speaking, the MAP-inference is NP-hard in general, but firstly, it is much easier than the partition function calculation, secondly, there are solvers for special cases, e.g., for log-supermodular models (which include functions 𝑓 which are negatives of cuts [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts[END_REF], [START_REF] Boykov | An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision[END_REF]) and those solvers are often efficient enough in practice. In this chapter, we will focus primarily on a subcase of supermodular potentials, namely negatives of graph cuts.

Parameter learning and Inference

In the standard supervised setting of structured prediction, we are given 𝑁 pairs of observations D = {(𝑥 𝑛 , 𝑦 𝑛 )} 𝑁 𝑛=1 , where 𝑥 𝑛 is a feature representation of the 𝑛-th object and 𝑦 𝑛 ∈ 𝑌 = 𝑌 1 × • • • × 𝑌 𝐷 𝑛 is a structured vector of interest (e.g., a sequence of tags, a segmentation MAP or a document summarization representation). In the standard linear model, the potential function 𝑓 (𝑦|𝑥) is represented as a linear combination: 𝑓 (𝑦|𝑥) = 𝑤 𝑇 Ψ(𝑥, 𝑦), where 𝑤 is a vector of weights and the structured feature map Ψ(𝑥, 𝑦) contains the relevant information for the feature-label pair (𝑥, 𝑦). To learn the parameters using the predefined probabilistic model, one can use the (regularized) maximum likelihood approach: (3.3.3) where 𝜆 > 0 is a regularization parameter and the likelihood 𝑃 (𝑦|𝑥, 𝑤) is defined as

max 𝑤 1 𝑁 𝑁 ∑ 𝑛=1 log 𝑃 (𝑦 𝑛 |𝑥 𝑛 , 𝑤) - 𝜆 2 ‖𝑤‖ 2 ,
𝑃 (𝑦|𝑥, 𝑤) = exp(𝑓 (𝑦|𝑥)) 𝑍(𝑓, 𝑥) = exp(𝑓 (𝑦|𝑥) -𝐴(𝑓, 𝑥)),
where 𝐴(𝑓, 𝑥) is the log-partition function (that now depends on 𝑥, since we consider conditional models). [START_REF] Hammersley | On the partition function and random maximum aposteriori perturbations[END_REF] proposed to learn parameters based on the Gumbel bound 𝐴 𝐺 (𝑓, 𝑥) instead of the intractable log-partition function:

log 𝑃 (𝑦|𝑥) = 𝑓 (𝑦|𝑥) -𝐴(𝑓, 𝑥) ≤ 𝑓 (𝑦|𝑥) -𝐴 𝐺 (𝑓, 𝑥) = 𝑓 (𝑦|𝑥) -E 𝑧 [ max 𝑦∈𝑌 { 𝐷 ∑ 𝑑=1 𝑧 𝑑 (𝑦 𝑑 ) + 𝑓 (𝑦) }] ≈ 𝑓 (𝑦|𝑥)- 1 𝑀 𝑀 ∑ 𝑚=1 max 𝑦 (𝑚) ∈𝑌 { 𝐷 ∑ 𝑑=1 𝑧 (𝑚) 𝑑 (𝑦 (𝑚) 𝑑 ) + 𝑓 (𝑦 (𝑚) ) } .
Hazan and Jaakkola [2012] considered the fully-supervised setup where labels 𝑦 𝑛 were given for all data points 𝑥 𝑛 . [START_REF] Shpakova | Parameter learning for log-supermodular distributions[END_REF] developed the approach, but also considered a setup with missing data (part of the labels 𝑦 𝑛 are unknown) for the small Weizmann Horse dataset from [START_REF] Borenstein | Combining top-down and bottom-up segmentation[END_REF]. Leveraging the additional stochasticity present in the Gumbel samples, [START_REF] Shpakova | Parameter learning for log-supermodular distributions[END_REF] extend the use of stochastic gradient descent, not on the data as usually done, but on the Gumbel randomization. It is equivalent to the choice of parameter 𝑀 = 1 for every gradient computation (but with a new Gumbel sample at every iteration). In our work, we use the stochastic gradient descent in a regime stochastic w.r.t. both the data and the Gumbel perturbations. This allows us to apply the method to large-scale datasets.

For linear models, we have 𝑓 (𝑦|𝑥) = 𝑤 𝑇 Ψ(𝑥, 𝑦) and Ψ(𝑥, 𝑦) is usually given or takes zero effort to compute. We assume that the gradient calculation ∇ 𝑤 𝑓 (𝑦|𝑥) = Ψ(𝑥, 𝑦) does not add complexity to the optimization algorithm. The gradient of log 𝑃 (𝑦|𝑥) is equal to

∇ 𝑤 𝑓 (𝑦|𝑥) -∇ 𝑤 max 𝑦∈𝑌 { 𝐷 ∑ 𝑑=1 𝑧 𝑑 (𝑦 𝑑 ) + 𝑓 (𝑦|𝑥) } = ∇ 𝑤 𝑓 (𝑦|𝑥) -∇ 𝑤 𝑓 (𝑦 * |𝑥),
where 𝑦 * lies in arg max of the perturbed optimization problem. The gradient of ⟨log 𝑃 (𝑦|𝑥)⟩ (the average over a subsample of data, typically a mini-batch) has the form Sample data mini-batch of small size 𝑇 (that is, 𝑇 pairs of observations)

⟨∇ 𝑤 𝑓 (𝑦|𝑥)⟩-⟨∇ 𝑤 𝑓 (𝑦 * |𝑥)⟩ = ⟨Ψ(𝑥,

4:

Calculate sufficient statistics ⟨Ψ(𝑥, 𝑦)⟩ 𝑒𝑚𝑝. from the mini-batch 5:

for 𝑡=1 to 𝑇 do 6:

Sample 𝑧 𝑑 (𝑦 𝑑 ) as independent Gumbels for all 𝑦 𝑑 ∈ 𝑌 𝑑 and for all 𝑑 Make a gradient step:

𝑤 ℎ+1 → 𝑤 ℎ + 𝛾 ℎ ( ⟨Ψ(𝑥, 𝑦)⟩ 𝑒𝑚𝑝. -⟨Ψ(𝑥, 𝑦 * )⟩ -𝜆𝑤 ℎ )
Note, that the classic log-likelihood formulation (3.3.3) is implicitly considering a "0-1 loss" 𝑙 0-1 (𝑦, ŷ) = [𝑦 ̸ = ŷ] as it takes probability of the entire output object 𝑦 𝑛 conditioned on the observed feature representation 𝑥 𝑛 . However, in many structured-output problems 0-1 loss evaluation is not an adequate performance measure. The Hamming or weighted Hamming losses that sum mistakes across the 𝐷 elements of the outputs, are more in demand as they count misclassification per element.

Marginal probability estimation

Either at testing time (to provide an estimate of the uncertainty of the model) or at training time (in the case of weak supervision, see Section 3.5), we need to compute marginal probabilities for a single variable 𝑦 𝑑 out of the 𝑑 ones, that is,

𝑃 (𝑦 𝑑 |𝑥) = ∑ 𝑦 -𝑑 𝑃 (𝑦 -𝑑 , 𝑦 𝑑 |𝑥),
where 𝑦 -𝑑 is a sub-vector of 𝑦 obtained by elimination of the variable 𝑦 𝑑 . Following [START_REF] Hammersley | On the partition function and random maximum aposteriori perturbations[END_REF] and [START_REF] Shpakova | Parameter learning for log-supermodular distributions[END_REF], this can be obtained by taking 𝑚 Gumbel samples and the associated maximizers

𝑦 𝑚 ∈ 𝑌 = 𝑌 1 × • • • × 𝑌 𝐷 ,
and, for any particular 𝑑, counting the number of occurrences in each possible value in all the 𝑑-th components 𝑦 𝑚 𝑑 of the maximizers 𝑦 𝑚 . While this provides an estimate of the marginal probability, this is not an easy expression to optimize at it depends on several maximizers of potentially complex optimization problems. In the next section, we show how we can compute a different (and new) approximation which is easily differentiable and on which we can apply stochastic gradient descent.

Marginal Likelihood

In this section, we demonstrate the learning procedure for the element-decoupled losses. We consider the regularized empirical risk minimization problem in a general form: (3.4.1) where ℓ(𝑤, 𝑥, 𝑦) can take various forms from Table 3.1 and 𝜆 is the regularization parameter. The choice of the likelihood form is based on the problem setting such as presence of missing data and the considered test-time evaluation function.

max 𝑤 1 𝑁 𝑁 ∑ 𝑛=1 ℓ(𝑤, 𝑥 𝑛 , 𝑦 𝑛 ) - 𝜆 2 ‖𝑤‖ 2 ,

Hamming loss

The Hamming loss is a loss function that counts misclassification per dimension:

𝑙 ℎ (𝑦, ŷ) = 1 𝐷 𝐷 ∑ 𝑑=1 [𝑦 𝑑 ̸ = ŷ𝑑 ].
For this type of loss instead of the classic log-likelihood objective it is more reasonable to consider the following decoupling representation from Table 3.1:

ℓ(𝑤, 𝑥, 𝑦) = 𝐷 ∑ 𝑑=1 log 𝑃 (𝑦 𝑑 |𝑥, 𝑤), (3.4.2)
where

𝑃 (𝑦 𝑑 |𝑥, 𝑤) = ∑ 𝑦 -𝑑 exp(𝑓 (𝑤, 𝑦 -𝑑 |𝑦 𝑑 , 𝑥) -𝐴(𝑓, 𝑥)) = exp(𝐵(𝑓, 𝑦 𝑑 , 𝑥) -𝐴(𝑓, 𝑥))
is the marginal probability of the single element 𝑦 𝑑 given the entire input 𝑥, and

𝐵(𝑓, 𝑦 𝑑 ) = log ∑ 𝑦 -𝑑 exp(𝑓 (𝑤, 𝑦 -𝑑 |𝑦 𝑑 )),
where

𝑦 -𝑑 ∈ 𝑌 1 × . . . 𝑌 𝑑-1 × 𝑌 𝑑+1 × • • • × 𝑌 𝐷 .
Thus, the log-marginal probability may be obtained from the difference of two log-partition functions (which we will approximate below with Gumbel samples). This idea of considering the marginal likelihood was proposed by [START_REF] Kakade | An alternate objective function for markovian fields[END_REF]. Our contribution is to consider the approximation by "perturb-and-MAP" techniques. We thus have a new objective function ℓ(𝑤, 𝑥, 𝑦):

ℓ(𝑤, 𝑥, 𝑦) = 𝐷 ∑ 𝑑=1 [(𝐵(𝑓, 𝑥, 𝑦 𝑑 ) -𝐴(𝑓, 𝑥))] ,
and now the following approximation could be applied:

𝐴(𝑓 ) ≈ 𝐴 𝐺 (𝑓 ) = E 𝑧 { max 𝑦∈𝑌 𝐷 ∑ 𝑑=1 𝑧 𝑑 (𝑦 𝑑 ) + 𝑓 (𝑦) } , 𝐵(𝑓 |𝑦 𝑑 ) ≈ 𝐵 𝐺 (𝑓 |𝑦 𝑑 ) = E 𝑧 { max 𝑦 -𝑑 ∈𝑌 -𝑑 𝐷 ∑ 𝑠=1:𝑠̸ =𝑑 𝑧 𝑠 (𝑦 𝑠 ) + 𝑓 (𝑦 -𝑑 |𝑦 𝑑 ) } .
It is worth noting, that the approximation is not anymore an upper bound of the marginal likelihood; moreover it is a difference of convex functions. Remarkably, the objective function exactly matches the log-likelihood in the case of unary potentials (separable potential function) as the log-likelihood function becomes the sum of marginal likelihoods.

As noted, the objective ℓ(𝑤, 𝑥, 𝑦) is not convex anymore, but it is presented as the difference of two convex functions. We can still try to approximate with stochastic gradient descent (which then only converges to a stationary point, typically a local minimum). Algorithm 2 describes the implementation details. Sample data mini-batch of small size 𝑇 (that is, 𝑇 pairs of observations)

4:
for 𝑡=1 to 𝑇 do 5:

Sample 𝑧 𝑑 (𝑦 𝑑 ) as independent Gumbels for all 𝑦 𝑑 ∈ 𝑌 𝑑 and for all 𝑑 6:

Find 𝑦 * 𝐴 ∈ arg max 𝑦∈𝑌 { 𝐷 ∑ 𝑑=1 𝑧 𝑑 (𝑦 𝑑 ) + 𝑓 (𝑦) } 7:
for 𝑑=1 to 𝐷 do 8:

Find 𝑦 * 𝐵 ∈ arg max 𝑦 -𝑑 ∈𝑌 -𝑑 { 𝐷 ∑ 𝑠=1:𝑠̸ =𝑑 𝑧 𝑠 (𝑦 𝑠 ) + 𝑓 (𝑦 -𝑑 |𝑦 𝑑 )} 9:
Make a gradient step:

𝑤 ℎ+1 → 𝑤 ℎ + 𝛾 ℎ ( ⟨ ⟨Ψ(𝑥, 𝑦 * 𝐵 )⟩ -Ψ(𝑥, 𝑦 * 𝐴 ) ⟩ -𝜆𝑤 ℎ )
Acceleration trick. Interestingly we can use the same Gumbel perturbation realizations for approximating 𝐴 𝐺 (𝑓 ) and 𝐵 𝐺 (𝑓 |𝑦 𝑑 ) through an empirical average. On the one hand, this restriction should not influence on the result as with a sufficient large averaging number 𝑀 , 𝐴 𝐺 (𝑓 ) and 𝐵 𝐺 (𝑓 |𝑦 𝑑 ) converges to their expectations. This is the same for stochastic gradients: on every iteration, we use a different Gumbel perturbation, but we share this one for the estimation of the gradients of 𝐴 𝐺 (𝑓 ) and 𝐵 𝐺 (𝑓 |𝑦 𝑑 ). This allows us to save some computations as shown below, while preserving convergence (the extra correlation added by using the same samples for the two gradients does not change the unbiasedness of our estimates). Moreover, if 𝑦 * 𝐴 has the same label value 𝑦 𝑑 as the ground truth, then the MAP inference problem for 𝑦 * 𝐴 exactly matches the one for 𝑦 * 𝐵 (with the element 𝑦 𝑑 fixed from the ground truth). Then 𝑦 * 𝐴 = 𝑦 * 𝐵 and the corresponding difference of gradients gives zero impact into the gradient step. This fact allows us to reduce the number of MAPinference problems. We should thus calculate 𝑦 * 𝐵 only for those indices 𝑑 that leads to a mismatch between 𝑑-th label of 𝑦 * 𝐴 and the ground truth one. Remarkably, during the convergence to the optimal value, the reduction will occur more often and decrease the execution time with the number of iteration increase. Besides that, in the experiments with graph cuts in Section 3.6 we use dynamic graph cut algorithm for solving several optimization problems of similar structure (here 𝐷 marginal probabilities calculation).

We describe it in more details in Section 3.4.3.

Weighted Hamming loss

The weighted Hamming loss is used for performance evaluation in the models, where each dimension has its own level of importance, e.g., in an image segmentation with superpixels, proportional to the size of superpixels. It differs from the usual Hamming loss in this way:

𝑙 ℎ (𝑦, ŷ) = 1 𝐷 𝐷 ∑ 𝑑=1 𝜃 𝑑 (𝑦 𝑑 )[𝑦 𝑑 ̸ = ŷ𝑑 ].
Thus we consider a dimension-weighted model as it can be adjusted for the problem of interest that gives the model more flexibility. The optimization problem of interest is transformed from the previous case by weighted multiplication: .4.3) To justify this objective function, we notice that in the case of unit weights, the weighted loss and objective function (3.4.3) match the loss and the objective from the previous section. Furthermore, 𝑦 𝑑 with a large weight 𝜃 𝑑 (𝑦 𝑑 ) puts more importance towards making the right prediction for this 𝑦 𝑑 , and that is why we put more weight on the 𝑑-th marginal likelihood. This corresponds to the usual rebalancing used in binary classification (see, e.g., [START_REF] Bach | Considering cost asymmetry in learning classifiers[END_REF] and references therein). Then, the algorithm for this case duplicated the one for the usual Hamming loss and the acceleration trick can be used as well.

ℓ(𝑤, 𝑥, 𝑦) = 𝐷 ∑ 𝑑=1 𝜃 𝑑 (𝑦 𝑑 ) [(𝐵(𝑓, 𝑥, 𝑦 𝑑 ) -𝐴(𝑓, 𝑥))] . ( 3

Scalable algorithms for graph cuts

As a classical efficient MAP-solver for pairwise potentials problem we will use graph cut algorithms from [START_REF] Boykov | An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision[END_REF]. The function 𝑓 (𝑦|𝑥) should then be supermodular, i.e., with pairwise potentials, all pairwise weights of 𝑤 should remain negative.

In both Sections 3.4.1 and 3.4.2 we can apply the dynamic graphcut algorithm proposed by [START_REF] Kohli | Dynamic graph cuts for efficient inference in markov random fields[END_REF], which is a modification of the Boykov-Kolmogorov graphcut algorithms. It is dedicated to situations when a sequence of graphcut problems with slightly different unary potentials need to be solved. Then, instead of solving them separately, we can use the dynamic procedure and find the solutions for slightly different problems with less costs. This makes graphcut scalable for a special class of problems.

It can easily be seen that our sequence of problems

𝑦 * 𝐵 ∈ arg max 𝑦 -𝑑 ∈𝑌 -𝑑 { ∑ 𝑠=1:𝑠̸ =𝑑 𝑧 𝑠 (𝑦 𝑠 ) + 𝑓 (𝑦 -𝑑 |𝑦 𝑑 )}
for 𝑑 = 1, . . . , 𝐷 is a perfect application for the dynamic graph cut algorithm. At each iteration we solve 𝑇 sets of graph cut problems, each of set contains 1 + 𝑎 𝑡 problems solvable by the same dynamic cut, where 𝑎 𝑡 is the number of not matched pixels between 𝑦 * 𝐴 and ground truth 𝑦 𝑛 . Finally, using acceleration trick and dynamic cuts we reduce the gradient descent iteration complexity from 𝑇 ∑ 𝑡=1 (1 + 𝐷 𝑡 ) graphcut problems to 𝑇 dynamic graph cut problems. We make the approach scalable and can apply it for large datasets.

Parameter Learning in the Semisupervised Setup

In this section we assume the presence of objects with unknown labels in the train dataset. We can separate the given data in two parts: fully annotated data

D 1 = {(𝑥 𝑛 , 𝑦 𝑛 )} 𝑁
𝑛=1 as in the supervised case and unlabeled data D 2 = {𝑥 𝑙 } 𝐿 𝑙=1 . Then, the optimal model parameter 𝑤 is a solution of the following optimization problem:

max 𝑤 𝐿 1 (𝑤) + 𝜅𝐿 2 (𝑤) - 𝜆 2 ‖𝑤‖ 2 , (3.5.1)
where

𝐿 1 (𝑤) = ∑ 𝑁 𝑛=1 ℓ 1 (𝑤, 𝑥 𝑛 , 𝑦 𝑛 ), 𝐿 2 (𝑤) = ∑ 𝐿 𝑙=1 ℓ 2 (𝑤, 𝑥 𝑙
) and the parameter 𝜅 governs the importance of the unlabeled data. ℓ 1 (𝑤, 𝑥 𝑛 , 𝑦 𝑛 ) can have a form from the left column of the Table 3.1, and ℓ 2 (𝑤, 𝑥 𝑙 ) from the right one.

Marginal calculations.

It is worth reminding from Section 3.3.3, that we can approximate marginal probabilities 𝑞(𝑦) of holding 𝑦 𝑑 = 𝑘 along with the partition function approximation almost for free. This can be obtained by taking 𝑚 Gumbel samples and the associated maximizers 𝑦 𝑚 ∈ 𝑌 = 𝑌 1 × • • • × 𝑌 𝐷 , and, for any particular 𝑑, counting the number of occurrences in each possible value in all the 𝑑-th components 𝑦 𝑚 𝑑 of the maximizers 𝑦 𝑚 . The approximation accuracy depends on number of samples 𝑀 . To calculate this we already need to have a trained weight vector 𝑤 which we can obtained from the fully annotated dataset D = {(𝑥 𝑛 , 𝑦 𝑛 } 𝑁 𝑛=1 ). We will calculate 𝑞(𝑦) for the unlabelled data D 2 = {𝑥 𝑙 } 𝑙 𝑙=1 . Those marginal probabilities contain much more information than MAP inference for the new data as can be seen on the example in Figure 3-2. We believe that proper use of the marginal probabilities will help to gain better result than using labels from the MAP inference (which we observe in experiments).

It is worth noting that for the inference and learning phases we use a different number of Gumbel samples. During the learning phase, we incorporate the double stochastic procedure and use 1 sample per 1 iteration and 1 label. For the marginal calculation (inference) we should use large number of samples (e.g. 100 samples) to get accurate approximation.

We provide the sketch of the proposed optimization algorithm in Algorithm 3. The optimization of 𝐿 1 is fully supervised and this can be done with tools of the previous section. The optimization of 𝐿 2 requires the specification of ℓ 2 (𝑤, 𝑥), which we take as

ℓ 2 (𝑤, 𝑥) = 𝐷 ∑ 𝑑=1 ∑ 𝑦 𝑑 ∈{0,...,𝐾} 𝑞 𝑑 (𝑦 𝑑 ) log 𝑃 (𝑤, 𝑦 𝑑 |𝑥 𝑑 ) = 𝐷 ∑ 𝑑=1 ∑ 𝑦 𝑑 ∈{0,...,𝐾} 𝑞 𝑑 (𝑦 𝑑 )𝐵(𝑓 |𝑦 𝑑 ) -𝐷𝐴(𝑓 ),
that is, the average of the fully supervised cost function with labels generated from the model 𝑞. The term 𝐿 2 corresponds to the common way of treating unlabeled data through the marginal likelihood. The sub-algorithm for the optimization of ℓ 2 is presented as Algorithm 4. Sample data mini-batch of small size 𝑇 (that is, 𝑇 pairs of observations)

4:

for 𝑡=1 to 𝑇 do 5:

Sample 𝑧 𝑑 (𝑦 𝑑 ) as independent Gumbels for all 𝑦 𝑑 ∈ 𝑌 𝑑 and for all 𝑑 6:

Find 𝑦 * 𝐴 ∈ arg max 𝑦∈𝑌 { 𝐷 ∑ 𝑑=1 𝑧 𝑑 (𝑦 𝑑 ) + 𝑓 (𝑦) } 7:
for 𝑑=1 to 𝐷 and 𝑘=0 to 𝐾 do 8:

Find 𝑦 * 𝐵,𝑑,𝑘 ∈ arg max 𝑦 -𝑑 ∈𝑌 -𝑑 { 𝐷 ∑ 𝑠=1:𝑠̸ =𝑑 𝑧 𝑠 (𝑦 𝑠 ) + 𝑓 (𝑦 -𝑑 |𝑦 𝑑 = 𝑘)} 9:
Make a gradient step:

𝑤 ℎ+1 → 𝑤 ℎ + 𝛾 ℎ ( ⟨ ⟨ 𝐾 ∑ 𝑘=0 𝑞 𝑑 (𝑘)Ψ(𝑥, 𝑦 * 𝐵,𝑑,𝑘 )⟩ -Ψ(𝑥, 𝑦 * 𝐴 ) ⟩ -𝜆𝑤 ℎ )
Acceleration trick. Suppose, that 𝑦 𝑑 can take values in the range {0, . . . , 𝐾}. Again we use the same Gumbel perturbation for estimating 𝐴 𝐺 (𝑓 ) and 𝐵 𝑑𝑘𝐺 (𝑓 |𝑦 𝑑 = 𝑘) for all 𝑘 ∈ {0, . . . , 𝐾}. The consequence of using the same perturbations is that if the 𝑑-th label 𝑦 𝑑 of 𝑦 * 𝐴 takes value 𝑘, than the corresponding 𝑑-th gradient will cancel out with one of the 𝑦 * 𝐵𝑘 . Thus, we will calculate only 𝐾 (instead of 𝐾 + 1 labels) structured labels 𝑦 * 𝐵𝑙 (𝑙 ̸ = 𝑘) and reduce the number of optimization problems to be solved. Dynamic graph cuts are applied here as well.

Finally in Table 3.1 we see the relationships between the proposed objective functions. Firstly, the known labels 𝑦 𝑛 in the supervised case are equivalent to the binary marginal probabilities 𝑞(𝑦 𝑛 ) ∈ {0, 1} 𝐷 𝑛 . Secondly, the unit weights 𝜃 𝑑 (𝑦 𝑑 ) = 1 in the weighted Hamming loss are equivalent to the basic Hamming loss.

Partial labels. Another case that we would like to mention is annotation with partial labels, e.g., in an image segmentation application, the bounding boxes of the images are given. Then denote 𝑦 𝑔𝑖𝑣𝑒𝑛 as the set of given labels. In this setup the marginal probabilities become conditional ones 𝑞(𝑦 𝑑 |𝑦 𝑔𝑖𝑣𝑒𝑛 ) and to approximate this we need to solve several conditional MAP-inference problems. The objective function

ℓ 2 (𝑤) = ∑ 𝑑 ∑ 𝑦 𝑑 ∈{0,...,𝐾} 𝑞 𝑑 (𝑦 𝑑 |𝑦 𝑔𝑖𝑣𝑒𝑛 ) log 𝑃 (𝑤, 𝑦 𝑑 |𝑥 𝑑 , 𝑦 𝑔𝑖𝑣𝑒𝑛 )
remains feasible to optimize.

Experiments

The experimental evaluation consists of two parts: Section 3.6.1 is dedicated to the chain model problem, where we compare the different algorithms for supervised learning; Section 3.6.2 is focused on evaluating our approach for the pairwise model on a weakly-supervised problem.

OCR dataset

The given OCR dataset from [START_REF] Taskar | Learning structured prediction models: A large margin approach[END_REF] consists of handwritten words which are separated in letters in a chain manner, see examples in Figure 3-1. The OCR dataset contains 10 folds of ∼ 6000 words overall. The average length of the word is ∼ 9 characters. Two traditional setups of these datasets are considered: 1) "small" dataset when one fold is considered as a training data and the rest is for test, 2) "large" dataset when 9 folds of 10 compose the train data and the rest is the test data. We perform cross-validation over both setups and present results in Table 3.2.

As the MAP oracle we use the dynamic programming algorithm of [START_REF] Viterbi | Error bounds for convolutional codes and an asymptotically optimum decoding algorithm[END_REF]. The chain structure also allows us to calculate the partition function and marginal probabilities exactly. Thus, the CRF approach can be applied. We compare its performance with the structured SVM from [START_REF] Osokin | Minding the gaps for block Frank-Wolfe optimization of structured SVMs[END_REF], perturb-and-MAP [START_REF] Hammersley | On the partition function and random maximum aposteriori perturbations[END_REF]) and the one we propose for marginal perturb-and-MAP (as Hamming loss is used for evaluation). [START_REF] Taskar | Learning structured prediction models: A large margin approach[END_REF]).

The goal of this experiment is to demonstrate that the CRF approach with exact marginals shows a slightly worse performance as the proposed one with approximated marginals but correct Hamming loss. For the OCR dataset, we performed 10-fold cross-validation and the numbers of Table 3.2 correspond to the averaged loss function (Hamming loss) values over the 10 folds. As we can see from the result in Table 3.2, the approximate probabilistic approaches slightly outperforms the CRF on both datasets. The Gumbel approximation (with or without marginal likelihoods) does lead to a better estimation for the Hamming loss. Note that S-SVM performs better in the case of a larger dataset, which might be explained by stronger effects of model misspecification that hurts probabilistic models more than S-SVM [START_REF] Pletscher | Putting MAP back on the map[END_REF]).

HorseSeg dataset

The problem of interest is foreground/background superpixel segmentation. We consider a training set of images {𝑥 𝑛 } 𝑛=1...𝑁 that contain different numbers of superpixels. A hard segmentation of the image is expressed by an array 𝑦 𝑛 ∈ {0, 1} 𝐷 𝑛 , where 𝐷 𝑛 is the number of superpixels for the 𝑛-th image.

The HorseSeg dataset was created by [START_REF] Kolesnikov | Closed-form training of conditional random fields for large scale image segmentation[END_REF] and contains horse images. The "small" dataset has images with manually annotated labels and contains 147 images. The second "medium" dataset is partially annotated (only bounding Table 3.3: HorseSeg Dataset. Performance Comparison. method "small" "medium" "large" S-SVM+BCFW 12.3 10.9 10.9 perturb&MAP 20.9 21.0 20.9 w.m. perturb&MAP 11.6 10.9 10.9 boxes are given) and contains 5974 images. The remaining "large" one has 19317 images with no annotations at all. A fully annotated hold out dataset was used for the test stage. It consists of 241 images.

The graphical model is a pairwise model with loops. We consider log-supermodular distribution and thus, the max oracle is available as the graph cut algorithm by [START_REF] Boykov | An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision[END_REF]. Note that CRFs with exact inference cannot be used here.

Following [START_REF] Kolesnikov | Closed-form training of conditional random fields for large scale image segmentation[END_REF], for the performance evaluation the weighted Hamming loss is used, where the weight is governed by the superpixel size and foreground/background ratio in the particular image.

That is,

𝑙 ℎ (𝑦, ŷ) = 1 𝐷 𝐷 ∑ 𝑑=1 𝜃 𝑑 (𝑦 𝑑 )[𝑦 𝑑 ̸ = ŷ𝑑 ],
where

𝜃 𝑑 (𝑦 𝑑 ) = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 𝑉 𝑑 2𝑉 𝑓 𝑜𝑟𝑒𝑔𝑟𝑜𝑢𝑛𝑑 , if 𝑦 𝑑 = 1. 𝑉 𝑑 2𝑉 𝑏𝑎𝑐𝑘𝑔𝑟𝑜𝑢𝑛𝑑 , if 𝑦 𝑑 = 0.
𝑉 𝑑 is the size of superpixel 𝑑, 𝑉 𝑏𝑎𝑐𝑘𝑔𝑟𝑜𝑢𝑛𝑑 and 𝑉 𝑓 𝑜𝑟𝑒𝑔𝑟𝑜𝑢𝑛𝑑 are the sizes of the background and the foreground respectively. In this way smaller object sizes have more penalized mistakes.

Since we incorporate 𝜃(𝑦) into the learning process and for its evaluation we need to know the background and foreground sizes of the image, this formulation is only applicable for the supervised case, where 𝑦 𝑑 is given for all superpixels. However, in this dataset we have plenty of images with partial or zero annotation. For these set of images D 2 = {𝑥 𝑙 } 𝐿 𝑙=1 , we handle approximate marginal probabilities 𝑞 𝑙 𝑑 associated to the unknown labels. Using them we can approximate the foreground and background volumes:

𝑉 𝑙 𝑓 𝑜𝑟𝑒𝑔𝑟𝑜𝑢𝑛𝑑 ≈ 𝐷 𝑙 ∑ 𝑑=1 𝑞 𝑙 𝑑 and 𝑉 𝑙 𝑏𝑎𝑐𝑘𝑔𝑟𝑜𝑢𝑛𝑑 ≈ 𝐷 𝑙 ∑ 𝑑=1
(1 -𝑞 𝑙 𝑑 ). We provide an example of the marginal and MAP inference in Figure 3-2. The difference of the information compression between these two approaches is visually comparable. We believe that the smoother and accurate marginal approach should have a positive impact on the result, as the uncertainty about the prediction is well propagated.

As an example of the max-margin approaches we take S-SVM+BCFW from the paper of [START_REF] Osokin | Minding the gaps for block Frank-Wolfe optimization of structured SVMs[END_REF] which is well adapted to large-scale problems. For S-SVM+BCFW and perturb-and-MAP methods we use MAP-inference for labelling unlabelled data using 𝑤 1 (see Section 3.5). For the HorseSeg dataset (Table 3.3), the numbers correspond to the averaged loss function (weighted Hamming loss) values over the hold out test dataset. The results of the experiment in Table 3.3 demonstrate that the approaches taking into account the weights of the loss 𝜃 𝑑 (•) (S-SVM+BCFW and w.m. perturb&MAP) give a much better accuracy than the regular perturb&MAP. S-SVM+BCFW uses loss-augmented inference and thereby augments the weighted loss structure into the learning phase. Weighted marginal perturb-and-MAP plugs the weights of the weighted Hamming loss inside the objective log-likelihood function. Basic perturb-and-MAP does not use the weights 𝜃 𝑑 (•) and loses a lot of accuracy. This shows us that the predefined loss for performance evaluation has a significant influence on the result and should be taken into account.

Small dataset size influence. We now investigate the effect of the reduced "small" dataset. We preserve the setup from the previous section and the only thing that we change is 𝑁 , the size of the "small" fully-annotated dataset D 1 = {(𝑥 𝑛 , 𝑦 𝑛 )} 𝑁 𝑛=1 . The new "small" dataset is 10% the size of the initial one, i.e., only 14 images. By taking a small labelled dataset, we test the limit of supervised learning when few labels are present.

For the HorseSeg dataset (Table 3.4), the numbers correspond to the averaged loss function (weighted Hamming loss) values over the hold out test dataset. The results of this experiments are presented in Table 3.4. In this setup, the probabilistic approach "weighted marginal perturb-and-MAP" gains more than max-margin "S-SVM+BCFW". This could happen because of very limited fully supervised data. The learned parameter 𝑤 1 gives a noisy model and this noisy model produces a lot of noisy labels for the unlabeled data, while weighted perturb-and-MAP is more cautious as it uses probabilities that contain more information (see .

Acceleration trick impact

We now compare the execution time of the algorithm with and without our acceleration techniques (namely Dynamic Cuts [DC] and Gumbel Reduction [ GR]) to get an idea on how helpful they are. Table 3.5 shows the execution time for calculating all 𝑦 * 𝐵 (Algorithm 2) for different numbers of iterations on the HorseSeg small dataset. We conclude that the impact of DC does not depend on the total number of iterations always leading to acceleration around 1.3. For GR, acceleration goes from 3.5 for 100 iterations to 7.6 for one million iterations. Overall, we get acceleration of factor around 10 for one million iterations. method ∖ it 100 10 3 10 4 10 5 10 6 basic 0.9 9.2 89.5 900 8993 DC 0.7 6.9 69.0 696 7171 GR 0.3 2.1 15.5 133.4 1186 DC+GR 0.2 1.5 10.9 83.5 727 Table 3.5: Execution time comparison in seconds. HorseSeg small dataset.

Experiments analysis

The experiments results mainly show that not taking into account the right loss in the learning procedure is detrimental to probabilistic technique such as CRFs, while taking it into account (our novelty) improves results. Also, Tables 3.2 and3 .3 show that the proposed methods achieves (and sometimes surpasses) the level performance of the max-margin approach (with loss-augmented inference).

Further, we observed that the size of the training set influences the SSVM and perturb-and-MAP approaches differently. For smaller datasets, the max-margin approaches tend to lose information due to usage of the hard estimates for the unlabelled data (e.g. in Table 3.4: 16.1 against 14.4 for "medium" dataset without bounding boxes labeling).

Table 3.4 reports an experiment about using weakly-labeled data at the training stage (the results on the partially annotated "medium" dataset). This experiment studied the impact on the final prediction quality of the training set of "medium" size on top of the reduced "small" fully-labelled set. The results of Table 3.4 mean that the usage of our approach adopted to the correct test measure outperforms the default perturb-and-MAP by a large margin. Our approach also significantly outperformed the comparable baseline of SSVM due to reduced size of the "small" fully-labelled set.

Conclusion

In this chapter, we have proposed an approximate learning technique for problems with non-trivial losses. We were able to make marginal weighted log-likelihood for perturb-and-MAP tractable. Moreover, we used it for semi-supervised and weaklysupervised learning. Finally, we have successfully demonstrated good performance of the marginal-based and weighted-marginal-based approaches on the middle-scale experiments. As a future direction, we can go beyond the graph cuts and image segmentation application and consider other combinatorial problems with feasible MAP-inference, e.g., matching.

specifying a certain number of hyper-parameters, for example, for tree-structured sparsity, the weights to be given to each depth of the tree. The goal of this chapter is to propose data-driven estimation procedures for all of these hyper-parameters for a class of priors based on submodular functions. These include tree-structured priors or group-based priors and are commonly used in signal estimation problems [START_REF] Kim | Tree-guided group lasso for multi-task regression with structured sparsity[END_REF], [START_REF] Zhao | Grouped and hierarchical model selection through composite absolute penalties[END_REF]).

Contributions

In this chapter, we make the following contributions:

-We propose in Section 4.3 approximate maximum likelihood estimation for estimating parameters in continuous log-supermodular distributions, whose negative log-densities are commonly used as structured sparse priors in signal processing applications.

-We incorporate sparsity-inducing and challenging 𝑙 𝑝 norms treatment.

-We show in Section 4.4 how "perturb-and-MAP" ideas based on the Gumbel distribution and efficient discretization can be used to approximate the log-partition function for these models, which is the key step for approximate maximum likelihood estimation. Here, the fact that submodular functions can be efficiently minimized is crucial.

-We illustrate our estimation procedure in Section 4.5 on a set of experiments with flow-based priors and signal denoising.

Log-supermodular Distributions

As a probabilistic model, we are going to work with the family of log-supermodular distributions discussed by [START_REF] Djolonga | From MAP to Marginals: Variational Inference in Bayesian Submodular Models[END_REF]. These distributions are a special case of a Gibbs distribution over some variable 𝑥 ∈ X, which could at this point be discrete or continuously valued:

𝑑𝑝(𝑥) = 𝑒 -𝑓 (𝑥) 𝑍(𝑓 ) 𝑑𝜇(𝑥),
where 𝑑𝜇(𝑥) is a base measure, and 𝑓 (𝑥) is a potential function and the normalizer

𝑍(𝑓 ) = ∫ X 𝑒 -𝑓 (𝑥) 𝑑𝜇(𝑥).
It is worth noting that the partition function 𝑍(𝑓 ) is intractable in the general case, continuous or discrete, and its handling constitutes the core computational difficulty of probabilistic inference (Wainwright and Jordan [2008b]). If the potential function 𝑓 (𝑥) is submodular (see definition in the next section), then the distribution above is called log-supermodular (because -𝑓 is supermodular).

We use these models as they cover a broad family of distributions and allow us to perform an efficient gradient descent optimization of the likelihood objective due to submodularity (see below).

Supermodular and Submodular Functions

Submodular functions can be defined on sets X which are products of intervals. These functions have the particular property to have polynomial-time minimization algorithms [START_REF] Krause | Submodular function maximization[END_REF]). In this work, we consider the special cases X = R 𝑛 and X = {0, 1} 𝑛 , which will lead respectively to continuous and discrete submodular functions.

Discrete functions. A function 𝑓 : {0, 1} 𝑛 → R can be uniquely identified to a set-function defined over subsets of {1, . . . , 𝑛}, by defining 𝐹 (𝐴) = 𝑓 (1 𝐴 ) where 1 𝐴 ∈ {0, 1} 𝑛 is the indicator vector of the set 𝐴. It is then said submodular if it satisfies the following diminishing return property ∀𝐴, 𝐵 ⊆ {1, . . . , 𝑛} such that 𝐴 ⊆ 𝐵 and for all 𝑖, then 𝐹 (𝐴 ∪ 𝑖) -𝐹 (𝐴) ≥ 𝐹 (𝐵 ∪ 𝑖) -𝐹 (𝐵). Classical submodular functions include network flows, graph cuts, as well as concave functions of the cardinality and appear in many areas of signal processing and machine learning [START_REF] Fujishige | Submodular Functions and Optimization[END_REF], [START_REF] Bach | Learning with submodular functions: a convex optimization perspective[END_REF]). In particular, non-decreasing submodular functions are commonly used to penalize the support of signal in compressed sensing [START_REF] Obozinski | Convex relaxation for combinatorial penalties[END_REF]). Classical examples include group-based priors (counting the number of active groups) or tree-based priors (cardinality of the smallest rooted tree containing a given set), which are commonly used in signal processing [START_REF] Baraniuk | Model-based compressive sensing[END_REF]).

Continuous functions. In the continuous setting, a twice differentiable function 𝑓 : R 𝑛 → R is submodular if for all 𝑥 ∈ R 𝑛 , the cross second-order derivatives are non-positive, that is, for all 𝑖 ̸ = 𝑗, then 𝜕 2 𝑓 𝜕𝑥 𝑖 𝜕𝑥 𝑗 (𝑥) 0. For non-differentiable functions, see more details in [START_REF] Bach | Submodular functions: from discrete to continuous domains[END_REF]. Important examples include the Lovász extensions of discrete submodular functions defined on {0, 1} 𝑛 [START_REF] Bach | Learning with submodular functions: a convex optimization perspective[END_REF], [START_REF] Fujishige | Submodular Functions and Optimization[END_REF]). These extended functions are typically used for regularization as follows: when a vector 𝑤 whose support 𝐴 = {𝑖, 𝑤 𝑖 ̸ = 0} should have a small value 𝐹 (𝐴) where 𝐹 is the corresponding set-function, then the penalty 𝑓 (|𝑤|), where 𝑓 is the Lovász extension, and the absolute values are taken component-wise, is a natural convex relaxation [START_REF] Bach | Learning with submodular functions: a convex optimization perspective[END_REF], [START_REF] Obozinski | Convex relaxation for combinatorial penalties[END_REF]). Classical example is the Lovász extension of a set cover function and it has the following form:

𝑓 (|𝑤|) = ∑ 𝐴∈A 𝑑 𝐴 max 𝑎∈𝐴 |𝑤 𝑎 |,
where A is a set of subsets of {1, . . . , 𝑛} [START_REF] Bach | Structured sparsity through convex optimization[END_REF]); see more details in Section 4.5.2. Other non-convex but tighter relaxations are of the form 𝑓 (|𝑤| 𝑝 ) where the power is taken component-wise [START_REF] Halabi | Combinatorial penalties: Which structures are preserved by convex relaxations[END_REF]):

𝑓 (|𝑤| 𝑝 ) = ∑ 𝐴∈A 𝑑 𝐴 max 𝑎∈𝐴 |𝑤 𝑎 | 𝑝 .
By changing the free parameter 𝑝 we can manipulate the desired level of sparsity In a sparse set-up we would like to encourage the argument 𝛼 to contain plenty of zero elements and thus we can use the usual ℓ 1 -formulation [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF]

): min 𝛼 1 2 ‖𝑥 -𝐷𝛼‖ 2 + Λ(𝛼), 𝑠.𝑡. Λ(𝛼) = 𝑓 (|𝛼|),
where |𝛼| is meant component-wise. Here 𝑓 (|𝛼|) is encoding the structured sparsity, see an example in the Section 4.5.2. In this paper, for simplicity, we consider only orthonormal dictionaries.

Maximum likelihood for parameter learning. The main goal here is to consider Λ(𝛼) as a negative log-density with the functional parameter 𝑓 (•) and to find its optimal data-dependent parameters by maximum likelihood. To do so, we need to solve an optimization problem that involves the log-partition function in the way described below. We use the notation, for a submodular function 𝑓 defined on

R 𝑘 + , log 𝑍(𝑓 ) = 𝐴(𝑓 ) = log ∫ R 𝑘 + exp(-𝑓 (𝛼))𝑑𝛼.
Note that, by a simple change of variable, we have: 𝑑𝛽.

which changes 𝑔 to 𝑔 plus some separable terms. As a result, 𝑔 is submodular because it is obtained by a monotonic separable change of variable. Later in the discretization, this simply means adding a power 𝑝 in the discretization.

In order to learn the distribution of 𝛼 we need to perform maximum likelihood for the densities defined above. To perform this we need to approximate the log-partition function 𝐴, which is always convex, but usually hard to compute. Two questions arise:

(1) approximation of 𝐴(𝑓 ), and (2) parameterization of 𝑓 in a suitable form so that our optimization problems are easily solved.

Perturb-and-MAP

To perform an effective parameter learning, we would like to approximate the log-partition function 𝐴(𝑓 ) by making use of the "perturb-and-MAP" approach from [START_REF] Hammersley | On the partition function and random maximum aposteriori perturbations[END_REF]. One necessary point is to have an access to an efficient MAP-oracle, e.g., a graphcut solver for our particular example of flow-based priors (see Section 4.5.2).

Our algorithm is based on an approximation result from [START_REF] Hammersley | On the partition function and random maximum aposteriori perturbations[END_REF], which states that for any real-valued function 𝑔 defined on a discrete set T

= ∏ 𝑛 𝑖=1 T 𝑖 , then log ∑ 𝑡∈T 𝑒 𝑔(𝑡) E ℎ 1 ,...,ℎ𝑛∼Gumbel [ max 𝑡∈T ( 𝑔(𝑡) + 𝑛 ∑ 𝑖=1 ℎ 𝑖 (𝑡 𝑖 ) )]
, where Gumbel denotes the Gumbel distribution1 and with a collection {ℎ 𝑖 (𝑡 𝑖 )} 𝑖=1,...,𝑛 𝑡 𝑖 ∈T 𝑖 of independent Gumbel samples. This allows to define an upper-bound on the logpartition function, which is based on perturbing 𝑔 and performing maximization; it is thus efficient only for functions 𝑔 for which adding a separable terms leads to efficient optimization. This is exactly the case for negatives of submodular functions.

Via the proposed approximation we achieve a direct way of parameter learning in the discrete case. We can approximate the expectation over Gumbels using Monte-Carlo ideas by replacing the expectation by sampling. An efficient number of Gumbel samples depends on application, however in our experiments, 𝑀 = 100 seems to be enough for most of the setups. Following [START_REF] Shpakova | Marginal weighted maximum log-likelihood for efficient learning of perturb-and-map models[END_REF], in practice, when embedding the approximation result above in an optimization problem (for maximum likelihood) we can use stochastic gradient descent as a subroutine [START_REF] Shpakova | Marginal weighted maximum log-likelihood for efficient learning of perturb-and-map models[END_REF]), rather than using a fixed set of Gumbel samples.

Extension to the Continuous Case

To work with continuous data, we need to perform a discretization of the partition function to cast its calculation as a discrete optimization problem.

Approximation of 𝐴(𝑓 ) = log

0 = 𝑢 1 -𝑢 0 2 , 𝜋 𝑟-1 = 𝑢 𝑟-1 -𝑢 𝑟-2 2
, with the rest as 𝜋 𝑗 = 𝑢 𝑗+1 -𝑢 𝑗-1

2

. We then discretize the integral in the following way

Â(𝑓 ) = log ∑ 𝑧∈{0,...,𝑟-1} 𝑘 ( 𝑘 ∏ 𝑖=1 𝜋 𝑧 𝑖 ) exp(-𝑓 (𝑢 𝑧 𝑖 , . . . , 𝑢 𝑧 𝑘 )).
This is done by considering 𝑘𝑟 Gumbel variables ℎ 𝑖,𝑗 , 𝑖 ∈ {1, . . . , 𝑘} and 𝑗 ∈ {0, . . . , 𝑟 -1}, and the perturb-and-MAP approximation

Â𝐺 (𝑓 ) = E ℎ ( max 𝑧 { -𝑓 (•) + 𝑘 ∑ 𝑖=1 ℎ 𝑖,𝑧 𝑖 + 𝑘 ∑ 𝑖=1 log 𝜋 𝑧 𝑖 })
.

The separable term ℎ 𝑖,𝑧 𝑖 + log 𝜋 𝑧 𝑖 does not change the nature of the problem and the function remains submodular (but now defined on a finite number of values).

Decoding with MMSE

In this section we discuss several ways of performing inference. Decoding with MAP is straightforward, as we can get the approximation with a discrete solver or an exact solution via divide-and-conquer algorithms [START_REF] Groenevelt | Two algorithms for maximizing a separable concave function over a polymatroid feasible region[END_REF]), if the objective is convex. Instead of using MAP decoder we can use MMSE, which is known to be the Bayesian optimal classifier for the ℓ 2 -loss function (see, e.g., [START_REF]Should penalized least squares regression be interpreted as maximum a posteriori estimation[END_REF]). However, it is a challenge to calculate, and as one of our contributions in this work we propose a way of its approximation via perturb-and-MAP ideas. Indeed, we have the estimator

𝜓 𝑀 𝑀 𝑆𝐸 (𝑥) = E(𝛼|𝑥) 𝑃 (𝛼|𝑥) = 𝑃 (𝑥, 𝛼) ∑ 𝛼 𝑃 (𝑥, 𝛼) = exp(-𝑠(𝛼, 𝑥)) 𝑍(𝑠) , where 𝑠(𝛼, 𝑥) = 1 2𝜎 2 ||𝑥 -𝐷𝛼|| 2 + Λ(𝛼) = 1 2𝜎 2 ||𝐷 𝑇 𝑥 -𝛼|| 2 + Λ(𝛼)
for an orthonormal dictionary 𝐷. Then, the optimal estimator can be computed as:

𝜓 𝑀 𝑀 𝑆𝐸 (𝑥) = E(𝛼|𝑥) ≈ 1 𝑀 𝑀 ∑ 𝑚=1 arg min 𝛼𝑚 { 𝑠(𝛼 𝑚 ) - ∑ 𝑘 ℎ 𝑘 (𝛼 𝑘𝑚 ) } ,
as the expectation of the sufficient statistics is known to be a gradient of the logpartition function (Wainwright and Jordan [2008b]). Thus, we can get an approximate solution via approximate mean marginals and discretization.

Experiments

In the experiment section we illustrate the proposed parameter learning technique and compare the performance of two decoding approaches.

Synthetic Data. Experiments on decoding

We consider a one dimensional distribution 𝑃 (𝛼) = exp -|𝛼| 𝑍 . Closed-form solution for MMSE in discrete setup can then be calculated exactly for this one dimensional setup. We present denoising results for 10,000 randomly sampled data points in

Flow-based priors

We now show examples where our maximization problems can be cast a a max-flow / min-cut. Following Sections 6. We work with patches 8 × 8 of the 'boat' image(see . Half of them are considered as a train dataset to learn hyperparameteters of our priors, and another half as a test dataset. The original image is 512 × 512. After Gaussian noise was added to the test dataset, we try to denoise it with the proposed decoding techniques. We compare two ways of denoising: 1) model the prior directly for the pixels values and 2) model the prior for some wavelet coefficients.

We work under the following assumption

𝑃 (𝛼) = 𝑒 -𝑓 (𝛼) 𝑍(𝑓 ) ,
where 𝑓 (𝛼) = ∑ 𝐴∈A 𝑑 𝐴 max 𝑎∈𝐴 |𝛼 𝑎 |, where A is a set of predefined groups of variables and 𝛼 could be either the pixels values, either the wavelet coefficients. These groups could be based on the image grid or on the wavelet tree structure. Moreover, 𝐷 is a orthogonal wavelet transform (we use Haar wavelets).

Before we report results in Tables 4.2 and 4.3, we first comment on the discretization grid which influences a lot the performance. We use a uniform grid with a step Δ (either 1 or 10). The quality and the running time consumption of the discrete MAP and MMSE depends on the discretization grid and on Δ correspondingly. For a fair comparison we spend the same amount of time for both: MAP and Δ = 1 consume approximately the same amount of time as MMSE and Δ = 10. For MMSE we use 𝑀 = 100 Gumbel perturbations. The method "c. MAP" refers to the exact MAP solution achieved through direct continuous optimization via divide-and-conquer algorithm [START_REF] Groenevelt | Two algorithms for maximizing a separable concave function over a polymatroid feasible region[END_REF]). In Tables 4.2 and4 ||x-𝑥|| 2 and 𝑥 is a test image and x its prediction. We accompany the results with their standard deviation, across 10 different noise samples.

First set of experiments. We compare several models: "baseline", "unary" and "grid" models. "baseline" is a model with no learning and we set parameters 𝑑 𝐴 equal to zeros, "unary" corresponds to groups of size one (independent but not identically distributed variables), and "grid" corresponds to groups of size one and two in a grid manner. Summary: From Tables 4.2 and 4.3, we can see that we are able to learn the structure. In the primal approach (directly on pixels), the Markov random field ("grid") performs best, while for the wavelet approach, unary potentials already work best (with little gains in likelihood and denoising performance for the "tree" approach). Note that MMSE sometimes outperform MAP, however it depends on the discretization grid.

Conclusion

We proposed a new parameter learning approach in the non-trivial structured and continuous setup. We demonstrated its performance for a denoising problem. We also propose a way to perform approximate MMSE estimation which is the optimal decoder for ℓ 2 -loss evaluation. There is still some space for further investigation: non-uniform discretization grids, ℓ 𝑝 -looking norms such as done by [START_REF] Shervashidze | Learning the structure for structured sparsity[END_REF], where more heavy-tailed priors were considered.

Chapter 5 Conclusion and Future Work

Summary of the thesis

In this thesis, we focused on efficient parameter learning methods specific to various submodularity-based probabilistic models. We start our analysis with basic binary models, and expand it to discrete and continuous cases by the end of the thesis.

As our first contribution, we consider the problem of the partition function approximation for the specific log-supermodular distributions and provide some analysis. We presented two existed bounds and have shown that the L-Field bound by Djolonga andKrause [2014, 2015] is always inferior than the Gumbel bound by [START_REF] Hammersley | On the partition function and random maximum aposteriori perturbations[END_REF]. Then, we investigate the ability of these two approaches to perform approximate parameter learning under the maximum likelihood framework, and demonstrated that L-Field bound leads to a degenerate solution under linear parametrization, while the Gumbel approach can be applied successfully.

Another part of the thesis covers the learning rate guarantees. We bring the problem of parameter learning down to the optimization problem of a convex function and tackle it by gradient descent techniques. Our novelty is to incorporate the stochastic techniques over the Gumbel perturbations to accelerate the approach, thus, make it significantly faster. We can do either stochastic over Gumbels either double-stochastic over both data and Gumbels, if necessary.

Taking advantage of the probabilistic modelling, we also tackle the problem of partially missing data and propose an approach for parameter learning treatment. This extension does not conserve the convexity of the objective, but we use a local optimum as a solution.

We highlight that we extended the base approach into such directions as high scale and low levels of supervision. Besides that, we also covered various popular structured losses and reformulated the MLE objective in accordance with the picked loss function.

And the last part of the thesis contributes to the possible structure learning problem and sparsity-inducing solutions of the signal decoding problems. In this part we also move away from standard MAP inference solution and propose to deal with appropriate for squared loss MMSE approximation.

Perspectives

There are several possible directions that can be a prolongation of this work.

-Recently log-supermodular models with higher order potentials are gaining attention [START_REF] Zhang | Higher-order inference for multi-class logsupermodular models[END_REF], [START_REF] Tschiatschek | Learning probabilistic submodular diversity models via noise contrastive estimation[END_REF]). They possibly can be modeled and tackled by our approach.

-Another attractive field is to consider an application of the other submodular functions with feasible MAP solvers besides graphcuts, e.g., matching [START_REF] Taskar | Learning structured prediction models: A large margin approach[END_REF]). And in this framework we also can consider models, where MAP inference can be solved approximately, e.g., solvers based on linear programming relaxations [START_REF] Sontag | Tightening LP relaxations for MAP using message-passing[END_REF]).

-We can also try to tackle a dictionary learning problem altogether with prior learning in a expectation-maximization way [START_REF] Mairal | Online dictionary learning for sparse coding[END_REF]). This can give promising results.

-For the continuous distribution support our treatment is still limited and this can be extended for wiser non-uniform grid approximations, as long as Monte Carlo integration.

-Even larger-scale application can be considered if we see our approach as a last layer of a neural network [START_REF] Goodfellow | Deep learning[END_REF]). We can also possibly achieve better accuracy, if we incorporate neural networks during the feature engineering process.

-We could incorporate non-convex optimization approaches for parameter learning in the unsupervised case.

Résumé

Abstract

Probabilistic graphical models encode hidden dependencies between random variables for data modelling. Parameter estimation is a crucial part of handling such probabilistic models. These very general models have been used in plenty of fields such as computer vision, signal processing, natural language processing. We mostly focused on log-supermodular models, which is a specific part of exponential family distributions, where the potential function is assumed to be the negative of a submodular function. This property is handy for maximum a posteriori and parameter learning estimations. Despite the apparent restriction of the model, is covers a broad part of exponential families, since there are plenty of functions that are submodular, e.g., graph cuts, entropy and others. Probabilistic treatment is challenging for most models, however we were able to tackle some of the challenges at least approximately.

In this manuscript, we exploit perturb-and-MAP ideas for partition function approximation and efficient parameter learning. Moreover, the problem can be also interpreted as a structure learning task, where each estimated parameter or weight represents the importance of the corresponding term. We propose a way of approximate parameter estimation and inference for models where exact learning and inference is intractable in general case due to the partition function calculation complexity. The first part of the thesis is dedicated to theoretical guarantees. Given the log-supermodular models, we take advantage of the efficient minimization property related to submodularity. Introducing and comparing two existing upper bounds of the partition function, we are able to demonstrate their relation by proving a theoretical result.

We introduce an approach for missing data as a natural subroutine of probabilistic modelling. It appears that we can apply a stochastic technique over the proposed perturb-and-map approximation approach and still maintain convergence while make it faster in practice.

The second main contribution is an efficient and scalable generalization of the parameter learning approach. In this section we develop new algorithms to perform parameter estimation for various loss functions, different levels of supervision and we also work on the scalability. In particular, working with mostly graph cuts, we were able to incorporate various acceleration techniques.

As a third contribution we deal with the general problem of learning continuous signals. We focus on the sparse graphical models representations. We consider common sparsity-inducing regularizers as negative log-densities for the prior distribution. The proposed denoising techniques do not require choosing any precise regularizer in advance. To perform sparse representation learning, the signal processing community often uses symmetric losses such as 1 , but we propose to parameterize the loss and learn the weight of each loss component from the data. This is feasible via an approach which is similar to what we proposed in the previous sections. For all aspects of the parameter estimation mentioned above we performed computational experiments to illustrate the idea or compare with existing benchmarks, and demonstrate its performance in practice.

  Example of DGM. (c) Example of chain model.

Figure 1 - 1 :

 11 Figure 1-1: Undirected and directed graphical models on the same set of nodes 𝑉 .

Figure 1 - 2 :

 12 Figure 1-2: Example of a Linear-chain CRF.

Figure 1 - 3 :

 13 Figure 1-3: Conditionally independent subsets 𝑋 𝐴 and 𝑋 𝐵 given the subset 𝑋 𝐶 .

Figure 1 - 4 :

 14 Figure 1-4: Grid-structured CRF for image segmentation task.

Definition 3 .

 3 A base polyhedron 𝐵(𝐹 ) associated with a submodular function 𝐹 (𝑆) is defined in the following form:𝐵(𝐹 ) = {𝑠 ∈ R |𝑉 | , 𝑠(𝑉 ) = 𝐹 (𝑉 ), ∀𝐴 ⊆ 𝑉, 𝑠(𝐴) ≤ 𝐹 (𝐴)},where 𝑠(𝐴) = ∑ 𝑖∈𝐴 𝑠 𝑖 . Thus this polyhedron is defined as the intersection of hyperplanes {𝑠 ∈ R |𝑉 | , 𝑠 𝑇 𝑥 𝐴 = 𝑠(𝐴) ≤ 𝐹 (𝐴)}, where 𝑥 𝐴 are normals of the corresponding hyperplanes.

Definition 4 .

 4 For any submodular function 𝐹 (𝑥) : {0, 1} |𝑉 | → R we can define a Lovász extension 𝑓 (𝑤) : R |𝑉 | → R, which is a convex function and matches the initial function 𝐹 (𝑥) at all vertices of the hyper-cube. It has the following form:

  (𝐴)) is the log-partition function of the target distribution, and a log-partition function of any modular function 𝑠(𝐴) = 𝑠(𝑥) = 𝑠 𝑇 𝑥 has the form :

Figure 1 - 7 :

 17 Figure 1-7: Gumbel distribution cdf and pdf.

  Mean bounds with confidence intervals, 𝑐 = 3.

Figure 2 - 1 :

 21 Figure 2-1: Comparison of log-partition function bounds for different values of 𝑐. See text for details.

Figure 2 - 2 :

 22 Figure 2-2: Denoising of a horse image from the Weizmann horse database[START_REF] Borenstein | Combining top-down and bottom-up segmentation[END_REF]).

Find

  𝑦 * ∈ arg max 𝑦∈𝑌 { 𝐷 ∑ 𝑑=1 𝑧 𝑑 (𝑦 𝑑 ) + 𝑓 (𝑦) } 8:

Algorithm 2

 2 Double SGD for Marginal Likelihood Input: dataset D = {(𝑥 𝑛 , 𝑦 𝑛 )} 𝑁 𝑛=1 , number of iterations 𝐻, size of the mini-batch 𝑇 , stepsize sequence {𝛾 ℎ } 𝐻 ℎ=1 , regularization parameter 𝜆 Output: model parameters 𝑤 1: Initialization: 𝑤 = 0 2: for ℎ = 1 to 𝐻 do 3:

Algorithm 3 1 : 2 : 3 : 4 :

 31234 Sketch for the semisupervised algorithm. Input: fully annotated dataset D 1 = {(𝑥 𝑛 , 𝑦 𝑛 )} 𝑁 𝑛=1 , number of iterations 𝐻, size of the mini-batch 𝑇 , stepsize sequence {𝛾 ℎ } 𝐻 ℎ=1 , regularization param. 𝜆 Output: model parameters 𝑤 1 Initialization: 𝑤 1 = 0 Find 𝑤 1 via Algorithm 2 Input: fully annotated dataset D 1 = {(𝑥 𝑛 , 𝑦 𝑛 )} 𝑁 𝑛=1 , unlabeled dataset D 2 = {𝑥 𝑙 } 𝐿 𝑙=1 , number of iterations 𝐻, size of the mini-batch 𝑇 , stepsize sequence {𝛾 ℎ } 𝐻 ℎ=1 , regularization parameter 𝜆 Output: model parameters 𝑤 1,2 Initialization: 𝑤 1,2 = 𝑤 1 Calculate: 𝑞(𝑦) for unlabeled data via 𝑤 1 5: Find 𝑤 1,2 via mixture of Algorithms 2 and 4

Algorithm 4

 4 Double SGD for Unsupervised Learning Input: unlabeled dataset D 2 = {𝑥 𝑙 } 𝐿 𝑙=1 , parameter estimate 𝑤 1 , number of iterations 𝐻, size of the mini-batch 𝑇 , stepsize sequence {𝛾 ℎ } 𝐻 ℎ=1 , regularization parameter 𝜆 Output: model parameters 𝑤 1,2 1: Initialization: 𝑤 1,2 = 𝑤 1 2: for ℎ = 1 to 𝐻 do 3:

Figure 3 - 1 :

 31 Figure 3-1: Samples from OCR dataset[START_REF] Taskar | Learning structured prediction models: A large margin approach[END_REF]).

Figure 3 - 2 :

 32 Figure 3-2: Example of the marginal and MAP inference for an image from the HorseSeg database[START_REF] Kolesnikov | Closed-form training of conditional random fields for large scale image segmentation[END_REF]).

  log 𝑃 (𝛼) = 𝑓 (𝛼) + 𝐴(𝑓 ) + 𝑘 log 2. The next step is to go from 𝑝 = 1 to 𝑝 < 1. We can done the transformation by changing of variable 𝛽 = |𝛼| 𝑝 (component-wise): log 𝑍(𝑓 ) = log ∫ R 𝑘 exp(-𝑔(|𝛼| 𝑝 ))𝑑𝛼 = 𝑘 log 2

+

  exp(-𝑓 (𝛼))𝑑𝛼. In order to approximate 𝐴(𝑓 ), we are going to discretize each 𝛼 𝑖 into 𝑟 values 0 = 𝑢 0 < 𝑢 1 < • • • < 𝑢 𝑟-1 , and consider the measure on R + define as the weighted sum of Diracs 𝑟-1 ∑ 𝑗=0 𝜋 𝑗 𝛿(𝛽 = 𝑢 𝑗 ). A simple example (based on the trapezoidal rule) for the weights (𝜋 𝑗 ) is 𝜋

Figure 4 - 1 :

 41 Figure 4-1: Classical boat image for image denoising.

  .3 we compute signal-to-noise ratios (SNR), where SNR = 20 log 10 ||𝑥|| 2

  Definition 1 (Markov property). Two subsets of variables 𝑋 𝐴 and 𝑋 𝐵 are conditionally independent given the third subset of variables 𝑋 𝐶 or equivalently 𝑋 𝐴 ⊥ ⊥ 𝑋 𝐵 |𝑋 𝐶 , if every path from a node in 𝑋 𝐴 to a node in 𝑋 𝐵 in the graph 𝐺 intersects the subset 𝑋 𝐶 at least in one node.

  1} |𝑉 | → R is submodular, if for any two vertices of hypercube 𝑥, 𝑦 ∈ {0, 1} |𝑉 | and a canonical basis vector 𝑒 𝑖 ∈ {0, 1} |𝑉 | , such that 𝑥 ≤ 𝑦 component-wise and 𝑦 𝑖 = 0, the following is true:

	1, 1, 1) ∼ {a, b, c}
	Figure 1-6: Hypercube representation for the ground set 𝑉 = {𝑎, 𝑏, 𝑐}.
	Definition 2. The set function 𝐹 (𝑆) : 2 𝑉 → R is called submodular, if for any two
	subset 𝐴, 𝐵 ⊆ 𝑉 and an element 𝑠 ∈ 𝑉 , such that 𝐴 ⊆ 𝐵, 𝑠 / ∈ 𝐵, the diminishing
	marginal returns property holds:
	𝐹 (𝐴 ∪ {𝑠}) -𝐹 (𝐴) ≥ 𝐹 (𝐵 ∪ {𝑠}) -𝐹 (𝐵),
	or equivalently 𝐹 (𝑥) : {0,

Table 2 .

 2

		max-marg.	std	mean-marg.	std
	1%	0.5%	<0.1%	0.5%	<0.1%
	5%	0.9%	0.1%	1.0%	0.1%
	10%	1.9%	0.4%	2.1%	0.4%
	20%	5.3%	2.0%	6.0%	2.0%
	𝜋	max-marg. std mean-marg. std
	1%	1.0%	-	1.0%	-
	5%	3.5%	0.9%	3.6%	0.8%
	10%	6.8%	2.2%	7.0%	2.0%
	20%	20.0%	-	20.0%	-

2: Unsupervised denoising results. Level of noise 𝜋 is fixed.

Table 2 .

 2

3: Unsupervised denoising results. Level of noise 𝜋 is not fixed.

  𝑦)⟩ 𝑒𝑚𝑝. -⟨Ψ(𝑥, 𝑦 * )⟩, where ⟨Ψ(𝑥, 𝑦)⟩ 𝑒𝑚𝑝. denotes the empirical average over the data. Algorithm 1 contains this double stochastic gradient descent (SGD) with stochasticity w.r.t. both sampled data and Gumbel samples. The choice of the stepsize 𝛾 ℎ = 1 𝜆ℎ is standard for strongly-convex problems (Shalev-Shwartz et al. [2011]). Double SGD: stochasticity w.r.t. data and Gumbel samples Input: dataset D = {(𝑥 𝑛 , 𝑦 𝑛 )} 𝑁 𝑛=1 , number of iterations 𝐻, size of the mini-batch 𝑇 , stepsize sequence {𝛾 ℎ } 𝐻

	Algorithm 1

ℎ=1 , regularization parameter 𝜆 Output: model parameters 𝑤 1: Initialization: 𝑤 = 0 2: for ℎ = 1 to 𝐻 do 3:

Table 3 . 1 :

 31 Variants of the Objective Loss ℓ(𝑤, 𝑥, 𝑦) Function. {𝜃 𝑑 (𝑦 𝑑 )} 𝐷 𝑑=1 are the weights of the weighted Hamming loss, {𝑞 𝑑 (𝑦 𝑑 )} 𝐷 𝑑=1 are the marginal probabilities 𝑃 (𝑦 𝑑 |𝑥). 𝑞 𝑑 (𝑦 𝑑 ) log 𝑃 (𝑤, 𝑦 𝑑 |𝑥 𝑑 ) 𝑦 𝑑 ∈𝑌 𝑑 𝑞 𝑑 (𝑦 𝑑 )𝜃 𝑑 (𝑦 𝑑 ) log 𝑃 (𝑤, 𝑦 𝑑 |𝑥 𝑑 )

	Loss	Labelled Data		Unlabelled Data
	0-1	log 𝑃 (𝑤, 𝑦|𝑥)		log	∑	𝑃 (𝑤, 𝑦, 𝑥)
						𝑦∈𝑌
		𝐷		𝐷
	Hamming	∑	log 𝑃 (𝑤, 𝑦 𝑑 |𝑥 𝑑 )	∑	∑
		𝑑=1		𝑑=1	𝑦 𝑑 ∈𝑌 𝑑
	Weighted Hamming

𝐷 ∑ 𝑑=1 𝜃 𝑑 (𝑦 𝑑 ) log 𝑃 (𝑤, 𝑦 𝑑 |𝑥 𝑑 ) 𝐷 ∑ 𝑑=1 ∑

Table 3 .

 3 2: OCR Dataset. Performance Comparison.

	method	small dataset large dataset
	CRF	19.5 ± 0.4	13.1 ± 0.8
	S-SVM+BCFW	19.5 ± 0.4	12.0 ± 1.1
	perturb&MAP	19.1 ± 0.3	12.5 ± 1.1
	marg. perturb&MAP	19.1 ± 0.3	12.8 ± 1.2

Table 3 .

 3 4: Reduced HorseSeg Dataset. Performance Comparison.

	method	10% of "medium" "medium"
		"small" with bbox w/o bbox
	S-SVM+BCFW	17.3	14.0	16.1
	perturb&MAP	23.2	23.4	23.0
	w.m. p.&MAP	18.1	13.7	14.4

Table 4 .

 4 1. We can clearly see that MMSE outperforms MAP in this synthetic setup, in a set-up where the Gumbel approximation of the log-partition function is exact.

	Approach	ℓ 2 loss value
	continuous MAP solution	0.693
	MAP-oracle	0.744
	MMSE exact	0.626
	MMSE approx [𝑀 = 100]	0.632
	MMSE approx [𝑀 = 1000]	0.628

Table 4 .

 4 1: Denoising results. Synthetic data.

  3 and 6.4 of Bach [2013], we consider prior Λ(𝛼) as functions 𝑓 of the form 𝑓 (𝛼) = 𝐴 are parameters that we are interested of. In order to minimize the function 𝑓 (𝑢 𝑧 1 , . . . , 𝑢 𝑧 𝑘 ) + ∑ 𝑘 𝑖=1 𝑣 𝑖,𝑧 𝑖 , which is required inSection 4.4.1, following Ishikawa [2003], we can create a weighted directed graph where the (𝑠𝑡)-minimum cut gives the optimal solution, thus making the optimization efficient.

	∑	𝑑 𝐴 max
	𝐴∈A

𝑎∈𝐴

𝛼 𝑎 ,

where 𝑑

Table 4 .

 4 MAP 22.33 ± 0.02 22.32 ± 0.02 26.16 ± 0.02 MAP Δ10 21.95 ± 0.01 21.95 ± 0.01 25.36 ± 0.02 MAP Δ1 22.32 ± 0.02 22.31 ± 0.02 26.16 ± 0.02 MMSE 22.31 ± 0.02 22.30 ± 0.01 25.99 ± 0.01 2: Primal approach. Second set of experiments. The model "tree" corresponds to groups of size one and two that represent the wavelet quad-tree dependencies of two-dimensional Haar MAP 22.33 ± 0.02 24.96 ± 0.01 24.66 ± 0.01 MAP Δ10 21.98 ± 0.01 24.63 ± 0.08 24.44 ± 0.05 MAP Δ1 22.33 ± 0.02 25.02 ± 0.03 24.72 ± 0.03 MMSE 22.29 ± 0.03 25.17 ± 0.02 25.15 ± 0.02

	method	baseline	unary	grid
	c. Lhood	356.4	355.7	253.5

Table 4 .

 4 3: Wavelets approach.

  Les modèles graphiques probabilistes codent les dépendances entre les variables aléatoires et l'estimation des paramètres fait partie du traitement des modèles probabilistes. Ces modèles ont été utilisés dans des domaines tels que la vision par ordinateur, le traitement du signal, le traitement du langage naturel. Nous nous sommes concentrés sur les modèles log-supermodulaires, qui font partie des distributions familiales exponentielles, où la fonction potentielle est la fonction négative d'une fonction sous-modulaire. Malgré la restriction du modèle, est couvre une grande partie des familles exponentielles, car il ya beaucoup de fonctions qui sont sousmodulaires, par exemple, les coupes graphiques, entropie et autres. Le traitement probabiliste est habituellement difficile, mais nous avons été en mesure de relever certains des défis au moins approximativement. Nous exploitons les idées perturb-and-MAP pour l'approximation des fonctions de partition et l'apprentissage efficace des paramètres. Nous proposons une méthode d'estimation et d'inférence approximative des paramètres pour les modèles où l'apprentissage et l'inférence exacts sont difficiles à gérer dans le cas général. La première partie de la thèse est consacrée aux garanties théoriques. Étant donné les modèles logsupermodulaires, nous tirons parti de la propriété de minimisation efficace liée à la sous-modularité. En introduisant et en comparant deux limites supérieures existantes de la fonction de partition, nous démontrons leur relation en prouvant un résultat théorique. Nous introduisons une approche pour les données manquantes comme sous-routine naturelle de la modélisation probabiliste. Il semble que nous puissions appliquer une technique stochastique à l'approche d'approximation par perturbation et carte proposée tout en maintenant la convergence tout en la rendant plus rapide dans la pratique. Une autre contribution est une généralisation efficace et évolutive de l'approche d'apprentissage des paramètres. Nous développons des algorithmes pour effectuer l'estimation des paramètres pour diverses fonctions de perte, différents niveaux de supervision et nous travaillons sur l'évolutivité. Nous incorporons également certaines techniques d'accélération. Comme troisième contribution, nous abordons le problème général de l'apprentissage des signaux continus. Nous nous concentrons sur les représentations de modèles graphiques clairsemés et nous considérons les régularisateurs à faible densité comme des densités logarithmiques négatives pour la distribution antérieure. Les techniques de débruitage proposées ne nécessitent pas le choix d'un redresseur précis à l'avance. Pour effectuer un apprentissage de représentation clairsemée, la communauté du traitement du signal utilise souvent des pertes symétriques telles que 1 , mais nous proposons de paramétrer la perte et d'apprendre le poids de chaque composante de perte à partir des données. Nous avons effectué des expériences informatiques pour illustrer l'idée générale ou la comparer à des repères existants, et démontrer sa performance dans la pratique.

[START_REF] Hammersley | On the partition function and random maximum aposteriori perturbations[END_REF] shows a stronger difference, which we believe (after consulting with authors) is due to lack of convergence for the iterative algorithm solving the max-margin formulation.

The Gumbel distribution on the real line has cumulative distribution function 𝐹 (𝑧) = exp(-exp(-(𝑧 + 𝑐))), where 𝑐 is the Euler constant.

The Gumbel distribution on the real line has the cumulative distribution function 𝐹 (𝑧) = exp(-exp(-(𝑧 + 𝑐))), where 𝑐 is the Euler constant.

Acknowledgements

and sum up the work that was done. We also discuss the possible future directions.

We note that the content of this thesis is based on the publications/submissions we list below: (a) Chapter 2 is based on the work "Parameter Learning for Log-supermodular Distributions", T. Shpakova, F. Bach, published in Advances of the Conference on Neural Information Processing Systems (NIPS), 2016. (b) Chapter 3 is based on the work "Marginal Weighted Maximum Log-likelihood for Efficient Learning of Perturb-and-Map Models", T. Shpakova, F. Bach, A. Osokin, published in Proceedings of the Conference on Uncertainty in Artificial Intelligence (UAI), 2018. (c) Chapter 4 is based on the work "Hyper-parameter Learning for Sparse Structured Probabilistic Models", T. Shpakova, F. Bach, M. Davies, which is under review for the International Conference on Acoustics, Speech, and Signal Processing (ICASSP), 2019.

Chapter 4

Hyper-parameter Learning for Sparse Structured Probabilistic Models Abstract

In this chapter we consider the estimation of hyperparameters for regularization terms commonly used for obtaining structured sparse parameters in signal estimation problems, such as signal denoising. By considering the convex regularization terms as negative log-densities, we propose approximate maximum likelihood estimation for estimating parameters for continuous log-supermodular distributions, which is a key property that many sparse priors have. We then show how "perturb-and-MAP" ideas based on the Gumbel distribution and efficient discretization can be used to approximate the log-partition function for these models, which is a crucial step for approximate maximum likelihood estimation. We illustrate our estimation procedure on a set of experiments with flow-based priors and signal denoising.

This chapter is based on the work "Hyper-parameter Learning for Sparse Structured Probabilistic Models", T. Shpakova, F. Bach, M. Davies, is under review for the International Conference on Acoustics, Speech, and Signal Processing (ICASSP), 2019.

Introduction

Structured sparsity has emerged a versatile tool to go beyond plain parsimonious models. Indeed, taking into account the potential structure between the signal coefficients to be set to zero, both interpretability and predictive performance can be improved. Structured sparse priors can be handled in the various frameworks that have emerged for sparse methods, within convex optimization [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF], [START_REF] Jacob | Group lasso with overlap and graph lasso[END_REF]) or non-convex optimization [START_REF] Baraniuk | Model-based compressive sensing[END_REF], [START_REF] He | Exploiting structure in wavelet-based bayesian compressive sensing[END_REF], [START_REF] Huang | Learning with structured sparsity[END_REF]).

While encoding structure has several benefits, it comes with the extra task of in some sense (as we know that the Laplacian prior does not really provide sparse samples [START_REF] Davies | Restricted isometry constants where ℓ 𝑝 sparse recovery can fail for 0 < 𝑝 ≤ 1[END_REF])).

Discrete Log-supermodular Distributions

An important example of a submodular minimization problem is a graphcut problem, which is very applicable in a variety of fields. For example it is popular to incorporate graphcut functions as potential functions in Markov random fields (MRFs) and perform an image segmentation as a side effect [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts[END_REF], [START_REF] Boykov | An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision[END_REF]). Other examples can be mutual entropy, certain functions of eigenvalues of submatrices and many others [START_REF] Fujishige | Submodular Functions and Optimization[END_REF]). Seeing them as log-densities thus leads to a probabilistic treatment (see, e.g., [START_REF] Djolonga | From MAP to Marginals: Variational Inference in Bayesian Submodular Models[END_REF]), which similar to what we try to achieve in this paper for continuous distributions.

Continuous Log-supermodular Distributions

Several examples have been considered in other settings and then called "multivariatetotally positive of order 2", and include the multivariate logistic, Gamma distributions, as well as characteristic roots of random Wishart matrices [START_REF] Karlin | Classes of orderings of measures and related correlation inequalities. i. multivariate totally positive distributions[END_REF]). We propose to extend their use in signal processing, in order to learn parameters of the associated negative log-densities (which are submodular functions), in structured sparsity problems.

Log-partition Function for Bayesian Learning

Let us consider the following standard denoising problem, with 𝑥 = 𝐷𝛼 + 𝜀, where 𝜀 is a Gaussian noise ∼ N(0, 𝜎 2 𝐼). Given the noisy signal 𝑥 ∈ R 𝑛 and a dictionary 𝐷 ∈ R 𝑛×𝑘 , we try to recover the initial representation 𝛼 ∈ R 𝑘 . The decoding problem often has the following form:

where Λ(𝛼) serves as a regularizer. This formulation can be considered from three points of view: 1) as penalized least squares regression without any probabilistic meaning, as 2) maximum a posteriori (MAP) or 3) minimum-mean-square-error (MMSE) estimation within some probabilistic model [START_REF]Should penalized least squares regression be interpreted as maximum a posteriori estimation[END_REF]). The probabilistic interpretation allows to learn parameters of Λ(𝛼), which we aim to do here. The joint probability has the form: 𝑃 (𝑥, 𝛼) = 𝑃 (𝛼)𝑃 (𝑥|𝛼) = 𝑃 (𝛼)𝑃 (𝜀 = 𝑥 -𝐷𝛼),

where 𝑃 (𝛼) = exp(-Λ(𝛼))

𝑍 is a log-supermodular prior on 𝛼. Thus, we would like to parameterize and learn the prior 𝑃 (𝛼) from the training dataset before performing the denoising on the test dataset.

List of Figures