Piotr Bojanowski

Rémi Lajugie

Francis Bach

Ivan Laptev

Jean Ponce

Cordelia Schmid

Josef Sivic

Weakly Supervised Action Labeling in Videos Under Ordering Constraints

We are given a set of video clips, each one annotated with an ordered list of actions, such as "walk" then "sit" then "answer phone" extracted from, for example, the associated text script. We seek to temporally localize the individual actions in each clip as well as to learn a discriminative classifier for each action. We formulate the problem as a weakly supervised temporal assignment with ordering constraints. Each video clip is divided into small time intervals and each time interval of each video clip is assigned one action label, while respecting the order in which the action labels appear in the given annotations. We show that the action label assignment can be determined together with learning a classifier for each action in a discriminative manner. We evaluate the proposed model on a new and challenging dataset of 937 video clips with a total of 787720 frames containing sequences of 16 different actions from 69 Hollywood movies.

Introduction

Significant progress towards action recognition in realistic video settings has been achieved in the past few years [START_REF] Laptev | Learning realistic human actions from movies[END_REF][START_REF] Liu | Recognizing human actions by attributes[END_REF][START_REF] Niebles | Modeling Temporal Structure of Decomposable Motion Segments for Activity Classification[END_REF][START_REF] Sadanand | Action bank: A high-level representation of activity in video[END_REF][START_REF] Wang | Action recognition by dense trajectories[END_REF]. However action recognition is often cast as a classification or detection problem using fully annotated data, where the temporal boundaries of individual actions, e.g. in the form of presegmented video clips, are given during training. The goal of this paper is to exploit the supervisory power of the temporal ordering of actions in a video stream, as illustrated in figure 1.

Gathering fully annotated videos with accurately time-stamped action labels is quite time consuming in practice. This limits the utility of fully supervised machine learning techniques on large-scale data. Using data redundancy, weakly and semi-supervised methods are a promising alternative in this case. On the other hand, it is easy to gather videos with some level of textual annotation but poor temporal localization, from movie scripts for example. This type of weak supervisory signal has been used before in classification [START_REF] Laptev | Learning realistic human actions from movies[END_REF] and temporal localization [START_REF] Duchenne | Automatic annotation of human actions in video[END_REF] tasks. However, the crucial information on the ordering of actions Fig. 1. Examples of video clips with associated actions sequence annotations such as provided in our dataset. Both examples contain the same set of actions but occurring in a different order. In this work we use the type and order of events as a supervisory signal to learn a classifier of each action and temporally localize each action in the video.

has, to the best of our knowledge, been ignored so far in the weakly supervised setting. Following recent work on discriminative clustering [START_REF] Bach | DIFFRAC: a discriminative and flexible framework for clustering[END_REF][START_REF] Xu | Maximum Margin Clustering[END_REF], image [START_REF] Joulin | Discriminative Clustering for Image Cosegmentation[END_REF] and video [START_REF] Bojanowski | Finding Actors and Actions in Movies[END_REF] cosegmentation, we propose to exploit this information in a discriminative framework where both the action model and the optimal assignments under temporal constraints are learned together.

Related Work

The temporal ordering of actions, e.g. in the form of Markov models or action grammars, have been used to constrain action prediction in videos [START_REF] Hongeng | Large-scale event detection using semi-hidden markov models[END_REF][START_REF] Ivanov | Recognition of visual activities and interactions by stochastic parsing[END_REF][START_REF] Kwak | Scenario-based video event recognition by constraint flow[END_REF][START_REF] Laxton | Leveraging temporal, contextual and ordering constraints for recognizing complex activities in video[END_REF][START_REF] Ryoo | Recognition of composite human activities through context-free grammar based representation[END_REF][START_REF] Vu | Automatic video interpretation: A novel algorithm for temporal scenario recognition[END_REF]. These kinds of spatial and temporal constraints have been also used in the context of group activity recognition [START_REF] Amer | Monte carlo tree search for scheduling activity recognition[END_REF][START_REF] Khamis | Combining per-frame and per-track cues for multi-person action recognition[END_REF]. Similar to us, these papers exploit the temporal structure of videos, but focus on inferring action sequences from noisy but pre-defined action detectors, often in constrained surveillance and laboratory settings with a limited number of actions and static cameras. In contrast, in this work we explore the temporal structure of actions for learning action classifiers in a weakly supervised set-up and show results on challenging videos from feature length movies.

Related is also work on recognition of composite activities [START_REF] Rohrbach | Script Data for Attribute-Based Recognition of Composite Activities[END_REF], where atomic action models ("cut", "open") are learned given full supervision on a cooking video dataset. Composite activity models ("prepare pizza") are learned on top of the atomic actions, using the prediction scores for the atomic actions as features. Annotations are, however, used without taking into account the ordering of actions.

Temporal models for recognition of individual actions have been explored in e.g. [START_REF] Laptev | Learning realistic human actions from movies[END_REF][START_REF] Niebles | Modeling Temporal Structure of Decomposable Motion Segments for Activity Classification[END_REF][START_REF] Tang | Learning latent temporal structure for complex event detection[END_REF]. Implicit models in the form of temporal pyramids have been used with bag-of-features representations [START_REF] Laptev | Learning realistic human actions from movies[END_REF]. Others have used more explicit temporal models in the form of, e.g. latent action parts [START_REF] Niebles | Modeling Temporal Structure of Decomposable Motion Segments for Activity Classification[END_REF] or hidden Markov models [START_REF] Tang | Learning latent temporal structure for complex event detection[END_REF]. Contrary to these methods, we do not use an a priori model of the temporal structure of individual actions, but instead exploit the given ordering constraints between actions to learn better individual actions models.

Weak supervision for learning actions has been explored in [START_REF] Bojanowski | Finding Actors and Actions in Movies[END_REF][START_REF] Duchenne | Automatic annotation of human actions in video[END_REF][START_REF] Laptev | Learning realistic human actions from movies[END_REF]. These methods use uncertain temporal annotations of actions provided by movie scripts. Contrary to these works our method learns multiple actions simultaneously and incorporates temporal ordering constraints on action labels obtained, e.g. from the movie scripts.

Dynamic time warping algorithms (DTW) can be used to match temporal sequences, and are extensively used in speech recognition, e.g. [START_REF] Gold | Speech and Audio Signal Processing -Processing and Perception of Speech and Music, Second Edition[END_REF][START_REF] Rabiner | Fundamentals of speech recognition[END_REF]. In computer vision, the temporal order of events has been exploited in [START_REF] Nguyen | Joint segmentation and classification of human actions in video[END_REF], where a DTW-like algorithm is used at test time to improve the performance of nonmaximum suppression on the output of pre-trained action detectors.

Discriminative clustering is an unsupervised method that partitions data by minimizing a discriminative objective, optimizing over both classifiers and labels [START_REF] Bach | DIFFRAC: a discriminative and flexible framework for clustering[END_REF][START_REF] Xu | Maximum Margin Clustering[END_REF]. Convex formulations of discriminative clustering have been explored in [START_REF] Bach | DIFFRAC: a discriminative and flexible framework for clustering[END_REF][START_REF] Guo | Convex Relaxations of Latent Variable Training[END_REF]. In computer vision these methods have been successfully applied to co-segmentation [START_REF] Joulin | Multi-class cosegmentation[END_REF]. The approach presented in this paper is inspired by this framework, but adds to it the use of ordering constraints.

In this work, we make use of the Frank-Wolfe algorithm (a.k.a conditional gradient) to minimize our cost function. The Frank-Wolfe algorithm [START_REF] Frank | An algorithm for quadratic programming[END_REF][START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-free sparse convex optimization[END_REF] is a classical convex optimization procedure that permits optimizing a continuously differentiable convex function over a convex compact domain only by optimizing linear functions over the domain. In particular, it does not require any projection steps. It has recently received increased attention in the context of large-scale optimization [START_REF] Harchaoui | Conditional gradient algorithms for machine learning[END_REF][START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-free sparse convex optimization[END_REF].

Problem Statement and Contributions

The temporal assignment problem addressed in the rest of this paper and illustrated by Fig. 1 can be stated as follows: We are given a set of N video clips (or clips for short in what follows). A clip is defined as a contiguous video segment consisting of F frames, and may correspond, for example, to a scene (as defined in a movie script) or a collection of subsequent shots. Each clip is divided into T small time intervals (chunks of videos consisting of F/T = 10 frames in our case), and annotated by an ordered list of K elements taken from some action set A of size A = |A| (that may consist of labels such as "open door", "stand up", "answer phone", etc., as in Fig. 1 for example). Note that clips are not of the same length but for the sake of simplicity, we assume they are. We address the problem of assigning to each time interval of each clip one action in A, respecting the order in which the actions appear in the original annotation list (Fig. 2). Fig. 2. Right: The goal is to find assignment of video intervals 1 to T (x-axis) to the ordered list of action annotations a(1) to a(K) indexed by integer k from 1 to K (yaxis). Left: The ordered annotation index k is mapped, through mapping a to action labels from the set A, in that example a(3) ="Eat". To preserve the given ordering of annotations we only consider assignments M that are non-decreasing. One such assignment m is shown in red.

Contributions. We make the following contributions: (i) we propose a discriminative clustering model (section 2) that handles weak supervision in the form of temporal ordering constraints and recovers a classifier for each action together with the temporal localization of each action in each video clip; (ii) we design a convex relaxation of the proposed model and show it can be efficiently solved using the conditional gradient (Frank-Wolfe) algorithm (section 3); and finally (iii) we demonstrate improved performance of our model on a new action dataset for the tasks of temporal localization (section 6) and action classification (section 7). All the data and code are publicly available at http://www.di.ens.fr/willow/research/ordering.

Discriminative Clustering with Ordering Constraints

In this section we describe the proposed discriminative clustering model that incorporates label ordering constraints. The input is a set of video clips, each annotated with an ordered list of action labels specifying the sequence of actions present in the clip. The output is the temporal assignment of actions to individual time intervals in each clip respecting the ordering constraint provided by the annotations together with a learnt classifier for each action, common for all clips. In the following, we first formulate the temporal assignment of actions to individual frames as discriminative clustering (section 2.1), then introduce a parametrization of temporal assignments using indicator variables (section 2.2), and finally we describe the choice of a loss function for the discriminative clustering that leads to a convex cost (section 2.3).

Problem Formulation

Let us now formalize the temporal assignment problem. We denote by x n (t) in R d some local descriptor of video clip number n during time interval number t. For every k in {1, . . . , K}, we also define a n (k) as the element of A corresponding to annotation number k (Fig. 2). Note that the set of actions A itself is not ordered: even if we represent A by a table for convenience, the elements of this table are action labels and have no natural order. The annotations, on the other hand, are ordered, for example according to where they occur in a movie script, and are represented by some integer between 1 and K. Thus a n maps (ordered) annotation indices onto (unordered) actions, and depends of course on the video clip under annotation. Parts of any video clip may belong to the background. To account for this fact, a dummy label ∅ is inserted in the annotation list between every consecutive pair of actual labels.

Let us denote by M the set of admissible assignments on {1, . . . , T }, that is, the set of sequences m = (m 1 , . . . , m T ) with elements in {1, . . . , K} such that m 1 = 1, m T = K, and m t+1 = m t or m t+1 = m t + 1 for all t in {1, . . . , T -1} Such an assignment is illustrated in Fig. 2.

Let us also denote by F the space of classifiers of interest, by Ω : F → R some regularizer on this space and by ℓ : A × R A → R + an appropriate loss function. For a given clip n and a fixed classifier f , the problem of assigning the clip intervals to the annotation sequence can be written as the minimization of the cost function:

E(m, f, n) = 1 T T t=1 ℓ (a n (m t ), f (x n (t))) (1)
with respect to assignment m in M. The regularizer Ω prevents overfitting and we therefore define a scalar parameter λ to control this effect. Jointly learning the classifiers and solving the assignment problem corresponds to the following optimization problem:

min f ∈F N n=1 min m∈M E(m, f, n) + λΩ(f ). (2)

Parameterization Using an Assignment Matrix

As will be shown in the following sections, it is convenient to reformulate our problem in terms of indicator variables. The corresponding multi-class loss is ℓ : {0, 1} A × R A → R + , and the classifiers are functions f : R d → R A . For a clip n, let us define the assignment matrix Z n ∈ R T ×A which is composed of entries z n ta such that z n ta = 1 if the interval t of clip n is assigned to class a. Let Z n t denote the row vector of dimension A corresponding to the t-th row of Z n . The cost function E(m, f, n), defined in Eq. (1) can be rewritten as

1 T T t=1 ℓ (Z n t , f (x n (t))
). Note: To avoid cumbersome double summations, we suppose from now that we work with a single clip. This allows us to drop the superscript notation, we replace Z n by Z and skip the sum over clips. We also replace the descriptor notation x n (t) by x t and the row extraction notation Z n t by Z t . This is without Fig. 3. Illustration of the correspondence between temporal assignments (left) and associated valid assignment matrices that map action labels a to time intervals t (right).

m = [1,1,1, 2,2,2,3,3,4,4,4,4,5,5,6,6,6] a =
Left: a valid assignment non-decreasing mt = k. Right: the corresponding assignment matrix Z. One can build the assignment matrix Z given the assignment m and the annotation sequence a by putting a 1 at index (t, a(mt)) in Z for every t. One obtains m given Z by iteratively constructing a sequence of integers of length T such that mt+1 = mt if the t-th and (t + 1)-th row of Z are identical, and mt+1 = mt + 1 otherwise.

loss of generality, and our method as described in the sequel handles multiple clips with some simple bookkeeping. Because of temporal constraints, we want the assignment matrices Z to correspond to valid assignments m. This amounts to imposing some constraints on Z. Let us therefore define Z, the set of all valid assignment matrices as:

Z = Z ∈ {0, 1} T ×A | ∃ m ∈ M, s.t., ∀ t, Z ta = 1 ⇐⇒ a(m t ) = a . (3)
There is a bijection between the sets Z and M. For each m in M there exists a unique corresponding Z in Z and vice versa. Figure 3 gives an intuitive illustration of this bijection.

The set Z is a subset of the set of stochastic matrices (positive matrices whose rows sum up to 1), formed by the matrices whose columns consist of exactly K blocks of contiguous ones occurring in a predefined order (K = 6 in Fig. 3). There are as many elements in Z as ways of choosing (K -1) transitions among (T -1) possibilities, thus |Z| = T -1 K-1 , which can be extremely large in our setting (in our setting T ≈ 100 and K ≈ 10). Furthermore, it is very difficult to describe explicitly the algebraic constraints on stochastic matrices that define Z. This point will prove important in Sec. 3 when we propose an optimization algorithm for learning our model. Using these notations, Eq. ( 2) is equivalent to:

min f ∈F ,Z∈Z 1 T T t=1 ℓ (Z t , f (x t )) + λΩ(f ).
(4)

Quadratic Cost Functions

We now choose specific functions ℓ and f that will lead to a quadratic cost function. This choice leads, to a convex relaxation of our problem. We use multi-class linear classifiers of the form f (x) = x T W + b, where W ∈ R d×A and b ∈ R 1×A . We choose the square loss function, regularized with the Frobenius norm of W , because in that case the optimal parameters W and b can be computed in closed form through matrix inversion. Let X be the matrix in R T ×d formed by the concatenation of all 1 × d matrices x t . For this choice of loss and regularizer, our objective function can be rewritten using the matrices defined above as:

1 T T t=1 ℓ(Z t , f (x t )) + λΩ(f ) = 1 T Z -XW -b 2 F + λ 2 W 2 F . (5)
This is exactly a ridge regression cost. Minimizing this cost with respect to W and b for fixed Z can be done in closed form [START_REF] Bach | DIFFRAC: a discriminative and flexible framework for clustering[END_REF][START_REF] Hastie | The elements of statistical learning: data mining, inference and prediction[END_REF]. Setting the partial derivatives with respect to W and b to zero and plugging the solution back yields the following equivalent problem:

min Z∈Z Tr ZZ T B , where B = 1 T Π T (I T -X X T Π T X + T λI d -1 X T )Π T , (6)
and the matrix Π p is the p × p centering matrix I p -1 p 1 p 1 T p . This corresponds to implicitly learning the classifier while finding the optimal Z by solving a quadratic optimisation problem in Z. The implicit classifier parameters W and b are shared among all video clips and can be recovered in closed-form as:

W = (X T Π d X + λI) -1 X T Π T ZD 1/2 , b = 1 T 1 T (Z -Xw)D 1/2 . ( 7
)
3 Convex Relaxation and the Frank-Wolfe Algorithm

In Sec. 2, we have seen that our model can be interpreted as the minimization of a convex quadratic function (B is positive semidefinite) over a very large but discrete domain. As is usual for this type of hard combinatorial optimization problem, we replace the discrete set Z by its convex hull Z. This allows us to find a continuous solution of the relaxed problem using an appropriate and efficient algorithm for convex optimization.

The Frank-Wolfe Algorithm

We want to carry out the minimization of a convex function over a complex polytope Z, defined as the convex hull of a large but finite set of integer points defined by the constraints associated with admissible assignments. When it is possible to optimize a linear function over a constraint set of this kind, but other usual operations (like projections) are not tractable, a good way to optimize a convex objective function is to use the iterative Frank-Wolfe algorithm (a.k.a. conditional gradient method) [START_REF] Bertsekas | Nonlinear Programming[END_REF][START_REF] Frank | An algorithm for quadratic programming[END_REF]. We show in Sec. 3.2 that we can minimize linear functions over Z, so this is an appropriate choice in our case.

The idea behind the Frank-Wolfe algorithm is rather simple. An affine approximation of the objective function is minimized yielding a point Z * on the Fig. 4. Illustration of a Frank-Wolfe step (see [START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-free sparse convex optimization[END_REF] for more details). Left: the domain Z interest, objective function (red), and its linearization at current point (blue). Right: top view of Z. Note that, in the algorithm, we actually minimize a linear function at each step. Adding a constant to it does not affect the solution of the minimization problem. That is why, we depicted an hyperplane that seems shifted from the origin. edge of Z. Then a convex combination of Z * and the current point Z k is computed. This is repeated until convergence (see Alg. 1). The interpolation parameter γ can be chosen either by using the universal step size 2 p+1 , where p is the iteration counter (see [START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-free sparse convex optimization[END_REF] and references therein) or, in the case of quadratic functions, by solving a univariate quadratic equation. In our implementation, we use the latter. A good feature of the Frank-Wolfe algorithm is that it provides a duality gap (referred to as the linearization duality gap [START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-free sparse convex optimization[END_REF]) that can be used as a certificate of sub-optimality and stopping criterion. The procedure is described in the special case of our relaxed problem in Algorithm 1. Figure 4 illustrates one step of the optimization.

Minimizing Linear Functions over Z by Dynamic Programming

It is possible to minimize linear functions over the binary set Z. Simple arguments (see for instance Prop B.21 of [START_REF] Bertsekas | Nonlinear Programming[END_REF]) show that the solution over Z is also a solution over Z. We will therefore focus on the minimization problem on Z and keep in mind that it also gives a solution over Z as required by the Frank-Wolfe algorithm. Minimizing a linear function on Z amounts to solving the problem: min Z∈Z Tr C T Z = T t=1 A a=1 Z ta C ta , where C is a cost matrix in R T ×A . Using the equivalence between the assignment matrix (Z) and

k ← 0 while Tr(∇ f (Z k )(Z k -Z * )) ≥ ǫ do Compute the current gradient in Z, ∇ f (Z k ) = Z T k B. Choose Z * in arg min Z∈Z Tr(Z∇ f (Z k )) using dynamic programming.
Compute the optimal Frank-Wolfe step size γ.

Z k+1 = Z k + γ(Z * -Z k ) k ← k + 1. end
Algorithm 1: The Frank-Wolfe optimization procedure.

the plain assignment (m) representations (Fig. 3), this is equivalent to solving min m∈M T t=1 A a=1 1 a(mt)=a C ta . To better deal with the temporal structure of the assignment, let us denote by D ∈ R T ×K the matrix with entries D tk = C ta(k) . The minimization problem then becomes min m∈M T t=1 D tmt , which can be solved using dynamic time warping. Indeed, let us define for all t ∈ {1, . . . , T } and k ∈ {1, . . . , K}: P * t (k) = min m∈M t s=1 D sms . We can think of P * t (k) as the cost of the optimal path from (1, 1) to (t, k) in the graph defined by admissible assignments, and we have the following dynamic programming recursion: P * t (k) = D tk + min(P * t-1 (k -1), P * t-1 (k)). The optimal value P * T (K) can be computed in O(T K) using dynamic programming, by precomputing the matrix D, incrementally computing the corresponding P * t (k) values, and maintaining at each node (t, k) back pointers to the appropriate neighbors.

Rounding

At convergence, the Frank-Wolfe algorithm finds the (non-integer) global optimum Ẑ of Eq. ( 6) over Z. Given Ẑ, we want to find an appropriate nearby point Z in Z. The simplest geometric rounding scheme consists in finding the closest point of Z according to the Frobenius distance : min Z∈Z Ẑ -Z 2 F . Expanding the norm yields: Ẑ -Z 2 F = Tr( ẐT Ẑ) + Tr(Z T Z) -2Tr( ẐT Z). Since Ẑ is fixed, its norm is a constant. Moreover, since Z is an element of Z, its squared norm is constant and equal to T . The rounding problem is therefore equivalent to: min Z∈Z -2Tr( ẐT Z), that is to the minimization of a linear function over Z. This can be done, as in Sec. 3.2, using dynamic programming.

Practical Concerns

In this section, we detail some refinements of our model. First, we show how to tackle a semi-supervised setting where some time-stamped annotations are available. Second, we discuss how to avoid the trivial solutions, a common issue in discriminative clustering methods [START_REF] Bach | DIFFRAC: a discriminative and flexible framework for clustering[END_REF][START_REF] Guo | Convex Relaxations of Latent Variable Training[END_REF][START_REF] Joulin | Discriminative Clustering for Image Cosegmentation[END_REF].

Semi-supervised Setting

Let us suppose that we are given some fully annotated clips (in the sense that they are labeled with time-stamped annotations), corresponding to a total of L time intervals. For every interval l we have a descriptor X l in R d and a class label a l in A. We can incorporate this data by modifying the optimization problem as follows:

min f ∈F min m∈M E(m, f, n) + 1 L L l=1 ℓ(a l , f (X l )) + λΩ(f ). ( 8
)
The first term is the weakly supervised assignment cost and the second one is the loss on annotated data. This supervised model does not change the optimization procedure, which remains valid.

Minimum Size Constraints

There are two inherent problems with discriminative clustering. First, the constant assignment matrix is typically a trivial optimum because, as explained in [START_REF] Guo | Convex Relaxations of Latent Variable Training[END_REF], the optimization domain is symmetruc over the permutations of labels. Thanks to our temporal constraints, the set Z is not symmetric and thus we are not subject to this effect. The second difficulty is linked to the use of the centering matrix Π T in Eq. ( 6). Indeed, we notice that the constant vector of length T is an eigen vector of Π T and thus the column-wise constant matrices are trivial solutions to our problem. Due to the temporal ordering constraints these solutions are not admissible, but in practice we have noticed that the algorithm returned an assignment with the label ∅ being dominant. We cope with these issues by adding a class penalizing linear term and by weighting the different classes. We refer the reader to [START_REF] Bojanowski | Weakly Supervised Action Labeling in Videos Under Ordering Constraints[END_REF] for more details.

Dataset and Features

Dataset. Our input data consists of challenging video clips annotated with sequences of actions. One possible source for such data is movies with their associated scripts [START_REF] Bojanowski | Finding Actors and Actions in Movies[END_REF][START_REF] Duchenne | Automatic annotation of human actions in video[END_REF][START_REF] Laptev | Learning realistic human actions from movies[END_REF][START_REF] Sivic | Who are you?" -Learning person specific classifiers from video[END_REF]. The annotations provided by this kind of data are noisy and do not provide ground-truth time-stamps for evaluation. To address this issue, we have constructed a new action dataset, containing clips annotated by sequences of actions. We have taken the 69 movies from which the clips of the Hollywood2 dataset were extracted [START_REF] Laptev | Learning realistic human actions from movies[END_REF], and manually added full time-stamped annotation for 16 classes (12 of these classes are already present in Hollywood2).

To build clips that form our input data, we search in the annotations for action chains containing at least two actions. To do so, we pad the temporal action annotations by 250 frames and search for overlapping intervals. A chain of such overlapping annotations forms one video clip with associated action sequence in our dataset. In the end we obtain 937 clips, with the number of actions ranging from 2 to 11. We subdivide each clip into temporal intervals of length 10 frames. Clips contain on average 84 intervals, the shortest containing 11, the longest 289.

Feature representation. We have to define a feature vector for every interval of a clip. We build a bag-of-words vector x t per interval t. Recall that intervals are of length 10 frames. To aggregate enough features, we decided to pool features from the 30-frame-long window centered on the interval. We compute video descriptors following [START_REF] Wang | Action Recognition with Improved Trajectories[END_REF]. We generate a vocabulary of size 2000 for HOF features. We restricted ourselves to one channel to improve the running time, while being aware that by doing so we sacrifice some performance. In our informal experiments, we also tried the MBH channels yielding very close performance. We use the Hellinger kernel to obtain the explicit feature map by square-rooting the l 1 normalized histograms. Every data point is associated with a vector x t in R 2000 .

Fig. 5. Splitting of the data described in Sec. 6.

Action Labeling Experiments

Experimental Setup. To carry out the action labeling experiments, we split 90% of the dataset into three parts (Fig. 5) that we denote Sup (for supervised), Eval (for evaluation) and Val (for validation). Sup is the part of data that has time-stamped annotations, and it is used only in the semi-supervised setting described in Sec. 4.1. Val is the set of examples on which we automatically adjust the hyper-parameters for our method (λ, κ, D). In practice we fix the Val set to contain 5% of the dataset. This set is provided with fully time-stamped annotations, but these are not used during the cost optimization. None of the reported results are computed on this set. We evaluate the quality of the assignment on the Eval set. Note that we carry out the Frank-Wolfe optimization on the union of all three sets. The annotations from the Sup set are used to constrain Z in the semi-supervised setup while those from the Val set are only used for choosing our hyper parameters. The supervisory information used over the rest of the data are the ordered annotations without time stamps. Please also keep in mind that there are no "training" and "testing" phases per se in this primary assignment task. All our experiments are conducted over five random splits of the data. This allows us to present results with error bars.

Performance Measure. Several measures may be used to evaluate the performance of discriminative clustering algorithms. Some authors propose to use the output classifier to perform a classification task [START_REF] Duchenne | Automatic annotation of human actions in video[END_REF][START_REF] Xu | Maximum Margin Clustering[END_REF] or use the output partition of the data as a solution of the segmentation task [START_REF] Joulin | Discriminative Clustering for Image Cosegmentation[END_REF]. Yet another way to evaluate is to use a loss between partitions [START_REF] Hubert | Comparing partitions[END_REF] as in [START_REF] Bach | DIFFRAC: a discriminative and flexible framework for clustering[END_REF]. Note that because of temporal constraints, for every clip we have a set of corresponding (prediction, ground-truth) pairs. We have thus chosen to measure the assignment quality for every ground-truth action interval I * and prediction I as |I ∩ I * |/|I|. This measure is similar to the standard Jaccard measure used for comparing ensembles [START_REF] Jaccard | The distribution of the flora in the alpine zone[END_REF]. Therefore, with a slight abuse of notation, we refer to this measure as the Jaccard measure. This performance measure is well suited for our problem since it respects the following properties: (1) it is high if the action predicted is included in the ground-truth annotation, (2) it is low if the prediction is bigger than the annotation, (3) it is lowest if the prediction is out of the annotation, (4) it does not take into account the prediction of the background class. The score is averaged across all ground-truth intervals. The perfect score of 1 is achieved when all actions are aligned to the correct annotations, but accurate temporal segmentation is not required as long as the predicted labels are within the ground truth interval. Baselines. We compare our method to the three following baselines. All these are trained using the same features as the ones used for our method. For all baselines, we round the obtained solution Z using the scheme described in Sec. 3.3. Normalized Cuts (NCUT). We compare our method to normalized cuts (or spectral clustering) [START_REF] Shi | Normalized Cuts and Image Segmentation[END_REF]. Let us define B as the symmetric Laplacian of the matrix

E: B = I -D -1 2 ED -1 2
, where D = Diag (E1). E measures both the proximity and appearance similarity of intervals. For all (i, j) in {1, . . . , T } 2 , we compute:

E ij = e -α|i-j|-βd χ 2 (Xi,Xj ) 1 |i-j|<dmin , where d χ 2 is the Chi-squared distance.
More precisely, we minimize over all cuts Z the cost g(Z) = Tr ZZ T B . g is convex (B is positive semidefinite) and we can use the Frank-Wolfe optimization scheme developed for our model. Intuitively, this baseline is searching for a partition of the video such that time intervals falling into the same segments have close-by features according to the Chi-squared distance. Bojanowski et al. [START_REF] Bojanowski | Finding Actors and Actions in Movies[END_REF]. We also consider our own implementation of the weaklysupervised approach proposed in [START_REF] Bojanowski | Finding Actors and Actions in Movies[END_REF]. We replace our ordering constraints by the corresponding "at least one" constraints. When an action is mentioned in the sequence, we require it appears at least once in the clip. This corresponds to a set of linear constraints on Z. We adapt this technique in order to work on our dataset. Indeed, the available implementation requires storing a square matrix of the size of the problem. Instead, we choose to minimize the convex objective of [START_REF] Bojanowski | Finding Actors and Actions in Movies[END_REF] using the Frank-Wolfe algorithm which is more scalable. Supervised Square Loss (SL). For completeness, we also compare our method to a fully supervised approach. We train a classifier using the square loss over the annotated Sup set and score all time intervals in Eval. We use the square loss since it is used in our method and all other baselines. Weakly Supervised Setup. In this setup, all baselines except (SL) have only access to weak supervision in the form of ordering constraints. Figure 6 (left) illustrates the quality of the predicted asignmentss and compares our method to baselines. Our method performs better than all other weakly-supervised methods. Both the Bojanowski et al. and NCUT baselines have low scores in the weakly-supervised setting. This shows the advantage of exploiting temporal constraints as weak supervisory signal. The fully supervised baseline (blue) eventually recovers a better alignment than our method as the fraction of fully annotated data increases. This occurs (when the red line crosses the blue line) at the 25% mark, as the supervised data makes up for the lack of ordering constraints. Fully time-stamped annotated data are expensive to produce whereas movies scripts are often easy to get. It appears thus that manually annotated videos are not always necessary since good performance is reached simply by using weak supervision. Figure 7 shows the results for all weakly-supervised methods for all classes. We notice that we outperform the baselines on the most frequent classes (such as "Open Door", "Sit Down" and "Stand Up"). Semi-supervised Setup. Figure 6 (right) illustrates the performance of our model when some supervised data is available. The fraction of the supervised data is given on the x-axis. First, note that our semi-supervised method (red) is always and consistently (Cf error bars) above the square loss baseline (blue). Of course, during the optimization, our method has access to weak annotations over the whole dataset, and to full annotations on the Sup set whereas the SL baseline has access only to the latter. This demonstrates the benefits of exploiting temporal constraints during learning. The semi-supervised Bojanowski et al. baseline (orange) has low performance, but it improves with the amount of full supervision provided.

Classification Experiments

The experiments in the previous section evaluate the quality of the recovered assignment matrix Z. Here we evaluate instead the quality of the recovered classifiers on a held-out test set of data for an action classification task. We recover these classifiers as explained later in this section. We can treat them as K independent, one-versus-rest classifiers and use them to score the samples from the test set. We evaluate this performance by computing per-class precision and recall and report the corresponding average precision for each class. Experimental setup. The models are trained following the procedure described in the previous section. To test the performance of our classifiers, we use the held out set of clips. This set is made of 10% of the clips from the original data. The clips from this set are identical in nature to the ones used to train the models. We also perform multiple random splits to report results with error bars.

Recovering the classifiers. One of the nice features of our method is that we can estimate the implicit classifiers corresponding to our solution Z * . We do so using the expression from Eq. 7.

Baselines. We compare the classifiers obtained by our method to those obtained by the Bojanowski et al. baseline [START_REF] Bojanowski | Finding Actors and Actions in Movies[END_REF]. We also compare them to the classifiers learned using the (SL) baseline.

Weakly Supervised Setup. Classification results are presented in Fig. 8 (left).

We observe a behavior similar to the action labeling experiment. But the supervised classifier (SL) trained on the Sup set using the square loss (blue) always performs worse than our model (red). This can be explained by the fact that the proposed model makes use of mode data. Even though our model has only access to weak annotation, it can prove sufficient to train good classifiers. The weaklysupervised method from Bojanowski et al. (orange) is performing worst, exactly as in the previous task. This can be explained by the fact that this method does not have access to full supervision or ordering constraints.

Semi-supervised Setup. In the semi-supervised setting (Fig. 8 (left)), our method (red) performs better than the supervised SL baseline (blue). The action model we recover is consistently better than the one obtained using only fully supervised data. Thus, our method is able to perform well semi-supervised learning. The Bojanowski et al. baseline (orange) improves when the fraction of annotated examples increases. Nonetheless, we see that making use of ordering constraints as used by our method signigicantly improves over simple linear inequalities ("at least one" constraints as formulated in [START_REF] Bojanowski | Finding Actors and Actions in Movies[END_REF]).

Acknowledgements. This work was supported by the European integrated project AXES, the MSR-INRIA laboratory, EIT-ICT labs, a Google Research Award, a PhD fellowship from the EADS Foundation, the Institut Universitaire de France and ERC grants ALLEGRO, VideoWorld, Activia and Sierra.

Fig. 6 .

 6 Fig.6. Alignment evaluation for all considered models. Left: weakly-supervised methods. This graph is shown for various fractions of fully supervised data only to compare to the SL baseline. "Weak" methods do not make use of this supervision. Right: semisupervised methods. See project webpage [1] for qualitative results.

Fig. 7 .

 7 Fig. 7. Alignment performance for various weakly-supervised methods for all classes.

Fig. 8 .

 8 Fig. 8. Classification performance for various models. Left: weakly-supervised methods. Right: semi-supervised methods. See project webpage [1] for qualitative results.

⋆ WILLOW project-team, DI/ENS, ENS/INRIA/CNRS UMR 8548, Paris, France. ⋆⋆ SIERRA project-team, DI/ENS, ENS/INRIA/CNRS UMR 8548, Paris, France. ⋆ ⋆ ⋆ LEAR team, INRIA Grenoble Rhône-Alpes, France.


Alexandre Le Boité

Giuliano Orso

Cristiano Ciuti

Bose-Hubbard Model: Relation Between Driven-Dissipative Steady-States and Equilibrium Quantum Phases

Keywords: numbers: 42.50.Ar, 03.75.Lm, 42.50.Pq, 71.36.+c

We present analytical solutions for the mean-field master equation of the driven-dissipative Bose-Hubbard model for cavity photons, in the limit of both weak pumping and weak dissipation. Instead of pure Mott insulator states, we find statistical mixtures with the same second-order coherence g (2) (0) as a Fock state with n photons, but a mean photon number of n/2. These mixed states occur when n pump photons have the same energy as n interacting photons inside the nonlinear cavity and survive up to a critical tunneling coupling strength, above which a crossover to classical coherent state takes place. We also explain the origin of both antibunching and superbunching predicted by P-representation mean-field theory at higher pumping and dissipation. In particular, we show that the strongly correlated region of the associated phase diagram cannot be described within the semiclassical Gross-Pitaevski approach.

I. INTRODUCTION

Since the seminal work by Fisher et al. in 1989 [1], the Bose-Hubbard model and its celebrated Mott insulator to superfluid quantum phase transition has attracted a great deal of attention. In spite of being the simplest model of interacting bosons on a lattice, its relevance for the study of many-body phenomena in bosonic systems has become even more prominent since the experimental observation of the predicted phase transition with ultracold atoms in optical lattices [2]. At the time of the first theoretical investigations, the most relevant candidates for experimental verification were undoubtedly atoms or Cooper pairs. But with the tremendous experimental progress in quantum optics and quantum electrodynamics of the past twenty years [3,[START_REF]The Physics of Semiconductor Microcavities[END_REF], the exciting field of many-body physics has ceased to belong exclusively to the realm of genuine particles: the exploration of Bose-Hubbard physics in optical systems is now within reach. It has indeed been demonstrated that in a nonlinear optical medium, light behave as a quantum fluid of interacting photons [START_REF] Carusotto | [END_REF] and some of the most remarkable features of quantum fluids, such as superfluid propagation [6,7] or generation of topological excitations [8][9][10][11] have been observed in experiments with solid-state microcavities. It has also become possible to design arrays of nonlinear cavity resonators, such that the effective on-site photonphoton interactions are large enough to enter the strongly correlated regime [12][13][14][15].

The question of whether key features of equilibrium physics are still present when the bosons have a finite life time has been of crucial importance ever since the first theoretical proposals for implementing the Bose-Hubbard model with photons or polaritons [16][17][18]. Exploiting the analogies between the two models, early works were focused on phenomena close to the equilibrium Mott insulator-Superfluid transition [19,20]. More recent studies were devoted to the interesting driven-dissipative regime, where the cavity resonators are excited by a coherent pump which competes with cavity losses [21][22][23][24][25][26][27][28][29][30]. In such conditions, the role of dissipation is crucial and similarities with equilibrium phenomena are more difficult to uncover. In particular, our recent study of the 2D mean-field phase diagram in the thermodynamical limit showed that the system is driven into steady-state phases whose general properties and collective excitations may be, at first sight, very different from equilibrium [31].

In this context, it is important to clarify the relation between driven-dissipative and equilibrium models and gain more insight into the nature of the stationary states. To tackle this question, we focus in this article on the limit of weak dissipation and weak pumping. Analytical expressions for the density matrix and the relevant observables enable us to identify the non-classical states of light which are the most closely related to a Mott insulator. Such states may be prepared when multiphotonic absorption processes become resonant. We show that beyond a critical value of the hopping rate between neighboring cavities, these non-classical states disappear and the system is driven into a classical coherent state. These results cast light on the photon statistics observed for stronger pumping and stronger dissipation. In particular, they give a clear picture of the two phases presented in the bistability diagram of ref. [31]. We also show that the 'weakly-interacting' sector of such diagram is well understood by means of a Gross-Pitaevski approximation. As expected, this simplified approach fails when on-site repulsion is much larger than the coupling between sites.

The paper is organized as follows: The model is introduced in section II. Section III is devoted to the limit of weak pumping and weak dissipation. The density matrix along with the relevant observables are first computed for single cavity (III A), and then extended to coupled cavities (III B). In section IV we explore the Gross-Pitaevski regime and we conclude in section V.

II. THE MODEL

We consider a driven-dissipative Bose-Hubbard model under homogeneous coherent pumping describing a bidimensional square lattice of cavity resonators. The system is described by the following Hamiltonian [START_REF] Carusotto | [END_REF]:

H = - J z <i,j> b † i b j + N i ω c b † i b i + U 2 b † i b † i b i b i + F e -iωpt b † i + F * e iωpt b i , (1)
where b † i creates a boson on site i, J > 0 is the tunneling strength, and z = 4 is the coordination number. < i, j > indicates that tunneling is possible only between first-neighbors. U > 0 represents the effective on-site repulsion, F is the amplitude of the incident laser field, ω c is the bare cavity frequency and ω p the frequency of the pump. The dynamics of the many-body density matrix ρ(t) is described in terms of a Lindblad master equation:

i∂ t ρ = [H, ρ] + iγ 2 N i 2b i ρb † i -b † i b i ρ -ρb † i b i , ( 2
)
where γ is the dissipation rate. It is convenient to eliminate the time dependency in Eq.( 1) by performing a unitary transformation on the density matrix:

ρ → U ρU † , (3)
where U = e iωpt i b † i bi . This amounts to writing the density matrix in a frame rotating at the pump frequency ω p . In this rotating frame, the Hamiltonian governing the dynamics is now time independent:

H rf = - J z <i,j> b † i b j - N i ∆ωb † i b i + U 2 b † i b † i b i b i +F b † i +F * b i .
(4) We have introduced the detuning between the pump the bare cavity frequency ∆ω = ω p -ω c . While for equilibrium quantum gases, the chemical potential µ is a key quantity, in this non-equilibrium model the steady-state phases depend instead on the pump parameters F and ∆ω, which compete with γ. It is worth pointing out that this is a model which well describes a lattice of cavities whose extra-cavity environment is the electromagnetic vacuum (apart from the applied driving field).

Apart from being conceptually simple and rich in phenomenology, the Bose-Hubbard model has another advantage: mean-field theory gives very good qualitative results at equilibrium. This approach has thus been extended to the driven-dissipative model under homogeneous pumping [31]. It consists of a decoupling approximation b

† i b j → b † i b j + b j b † i
, that reduces the initial N -site Hamiltonian to an effective single-site problem:

H mf = -∆ωb † b+ U 2 b † b † bb+(F -J b )b † +(F * -J b * )b, (5)
the effective external coherent field being:

F → F ′ = F -J b , (6)
where the bosonic coherence b has to be determined self-consistently. An exact solution of the single-cavity problem may be obtained by using the complex Prepresentation for the density matrix [31][32][33]. In this phase space approach, the matrix ρ is expressed in a coherent state basis as:

ρ = C β Cα P (α, β) |α β * | β * |α dαdβ, (7)
where C β and C α are paths of integration on individual complex planes (α, β). This representation allows a mapping of the master equation into a Fokker-Planck equation for the function P (α, β). As a result, all one-time correlation functions in the steady-state can be computed exactly and are given by:

(b † ) j (b) j = 2F U 2j × Γ(c)Γ(c * ) Γ(c + j)Γ(c * + j) × F (j + c, j + c * , 8|F/U | 2 ) F (c, c * , 8|F/U | 2 ) , (8)
with c = 2(-∆ω -iγ/2)/U and the hypergeometric function

F (c, d, z) = ∞ n Γ(c)Γ(d) Γ(c+n)Γ(d+n) z n
n! , Γ being the gamma special function. These results for the single cavity can be directly applied to mean-field theory by performing the substitution of Eq.( 6) and solving the following self-consistent equation for b :

b = (F -J b ) ∆ω + iγ/2 × F (1 + c, c * , 8| F -J b U | 2 ) F (c, c * , 8| F -J b U | 2 ) . (9)
In order to characterize the state of the intracavity electromagnetic field, we will focus mainly on two observables: the mean photon density b † b , and the zero-delay second order autocorrelation function g (2) (0). The latter is expressed by:

g (2) (0) = b † b † bb b † b 2 , (10)
and gives information on the photon statistics. It was shown that the self-consistent parameter in Eq.( 5) is responsible for the appearance of tunneling-induced bistability. That is, for a wide range of parameters, the selfconsistent equation for b has multiple solutions, giving rise to two stable steady-states (for the same values of all the parameters). In particular, we have identified a 'low-density' phase, where the average number of photons per site is very low ( b † b ≪ 1) but fluctuations in the photon statistics are high (g (2) (0) ≫ 1). In the other stable phase, called 'high-density' phase, b † b ∼ 1 and the statistics is sub-poissonian (g (2) (0)) < 1). Interestingly, a related bistable behavior has been predicted for systems of driven-dissipative Rydberg atoms [34].

As mentioned in the introduction, the next question that comes to mind is how to connect these results with the better known physics of the equilibrium model and its quantum phase transition. Are the stationary states related in any way to a Mott insulator or a superfluid? The best way to answer this question is to explore the limit of weak pumping and weak dissipation, since in this regime the Hamiltonian and the dynamics tend to resemble more and more their equilibrium counterpart. It is the subject of the following section.

III. THE LIMIT OF WEAK PUMPING AND WEAK EXCITATION

A. Single Cavity: Exact Solution Multiphotonic resonances. At equilibrium, the ground state of an isolated site is a pure Fock state |n , with n fixed by the value of the chemical potential. Photons, on the other hand, have always to be injected inside the cavity. When the coupling to the external field is very weak, the only way to have b † b ≥ 1 is to be at resonance with multiphotonic absorption processes.

Absorption of photons is favored when n incident laser photons have the same energy than n photons inside the cavity, that is: nω p = nω c + U n(n -1)/2. Expressed in terms of the pump-cavity detuning, this relation reads:

U ∆ω = 2 n -1 . (11)
This resonance condition can be satisfied for n > 1 only if ∆ω > 0. If the pump is resonant with the bare cavity frequency, i.e ∆ω = 0, only single photons can be absorbed resonantly. There is no relation between U and ∆ω in this case and the density matrix is found by expanding the master equation in powers of F/U and γ/U . In the following, we will focus on the more interesting case of multiphotonic resonances (n > 1) and assume ∆ω > 0.

For simplicity, let us first look at H mf for J = 0 and U = 2∆ω (two-photon resonance). It appears that in the rotating frame, the vacuum |0 has the same energy that the two photon-state |2 in the absence of driving. This degeneracy is lifted by the coupling to the external field whose effect on the dynamics can be understood qualitatively in the following way. Suppose that at time t = 0 the intra-cavity field is in the vacuum state |0 . Since the vacuum is no longer an eigenstate of the Hamiltonian, but a linear combination of eigenstates

1 √ 2 (|0 + |2 ) and 1 √ 2 (|0 -|2
), the cavity field will start to oscillate between |0 and |2 . These Rabi oscillations will take place until the occurrence of a quantum jump, resulting from spontaneous emission processes inside the cavity. If the frequency splitting between the two eigenstates is very small when compared with the dissipation rate, the field will not have time to oscillate and will stay mostly in the vacuum state. However, the photons being still absorbed by pairs, the fluctuations in the photon distribution may be very high. On the contrary, Rabi oscillation will take place if the dissipation rate is sufficiently small. Besides, the time between quantum jumps will be longer, thereby reducing the probability of emitting two photons at once.

More quantitative results can be obtained by solving the master equation explicitly. To fully grasp the effect of the resonance, we will assume that F/∆ω ≫ γ/∆ω. We will work is the basis formed by the eigenstates of the total Hamiltonian (up to the lowest order in F/∆ω), given by:

|a = 1 √ 2 (|0 + |2 ), (12)
|b = 1 √ 2 (|0 -|2 ), (13)
|c = |1 .
Since the driving term only couples Fock states |m to |m ± 1 , the coupling of |0 and |2 is of order 2. The energies given by second-order perturbation theory are:

E a ≃ F 2 ∆ω (1 + √ 2), (14)
E b ≃ F 2 ∆ω (1 - √ 2), (15)
E c ≃ -∆ω. (16)
As expected, the energy splitting between states |a and |b is proportional to F 2 /∆ω. In the "dressed-states" basis, the dissipative term of the master equation couples populations and coherences of the density matrix. However, in the lowest order in F/∆ω and γ/∆ω, the coefficients ρ ac and ρ bc vanish. The master equation reads:

∂ t ρ aa = γ[ 1 2 ρ cc -ρ aa + 1 2 (ρ ab + ρ ba )], (17)
∂ t ρ bb = γ[ 1 2 ρ cc -ρ bb + 1 2 (ρ ab + ρ ba )], (18)
∂ t ρ ab = ∆ω[(-i2 √ 2 F 2 ∆ω 2 - γ ∆ω )ρ ab + γ 2∆ω ], (19)
∂ t ρ cc = γ[ρ aa + ρ bb -ρ cc -(ρ ab + ρ ba )]. (20)
The stationary value for ρ ab is then:

ρ ab = 1 2 1 1 + i2 √ 2F 2 γ∆ω , (21)
giving:

ρ ab + ρ ba = 1 1 + 8F 4 γ 2 ∆ω 2 = ξ. ( 22
)
All the other coefficients along with the observables can be expressed as functions of the parameter ξ. Namely, we have:

ρ 11 ≃ ρ cc = 1 2 (1 -ξ), (23)
ρ aa = ρ bb = 1 4 (1 + ξ), (24)
ρ 22 = 1 2 (ρ aa + ρ bb ) - 1 2 (ρ ab + ρ ba ) = 1 4 (1 -ξ). ( 25
)
For this we can compute the mean photon density and g (2) (0):

b † b = 1 -ξ, (26)
g (2) (0) = 1 2(1 -ξ) . (27)
A comparison between theses two expressions and the exact P -representation formula is shown on Fig. (1). The mean photon density and g (2) (0) are plotted as a function of ξ for F/∆ω = 10 -2 . In order to stay in the domain of validity of Eq.( 26) and Eq.( 27), γ/∆ω is ranging from

F 2 10∆ω 2 to F 5∆ω .
In this conditions, the approximations underlying the derivation are justified and the above expressions are very accurate.

The two limits F 2 γ∆ω ≪ 1 and F 2 γ∆ω ≫ 1 correspond to ξ → 1 and ξ → 0 respectively. In the first case, the photon density goes to zero as expected, and g (2) (0) diverges. In the other limit, the field is in a statistical mixture of 3 states and the density matrix in Fock space is:

ρ = 1 4 |0 0| + 1 2 |1 1| + 1 4 |2 2|. (28)
The photon density in then equal to one and g (2) (0) = 0.5. We see that in the particular case of multiphotonic resonances, one must be careful in discussing the limit F ∆ω , γ ∆ω ≪ 1 as the system behavior varies qualitatively depending on the ratio F 2 γ∆ω . The state obtained for ξ → 0 is highly nonclassical and the closest to a Fock state that one can hope for in this context.

Besides, this result is not limited to n = 2. In the general case of n-photon resonance, the coupling between |0 and |n is of order n in F and the energy splitting proportional to (F/∆ω) n . Therefore, a state similar to Eq.( 28) may be obtain in the limit

F n γ∆ω n-1 ≫ 1.
As shown in the Appendix, the corresponding density matrix is:

ρ (n) = 1 2 n n k=0 n k |k k|. ( 29
)
For n > 2, not only the two states |0 and |n are degenerate in the absence of driving, but so are all the states |k , |q with k + q = n. This degeneracy is reflected in the n → n-k symmetry of Eq.( 29). The state of Eq.( 29) is characterized by: Note that the value of g (2) (0) is the same as in the n thlobe of the equilibrium model (pure Fock state with n photons). We emphasize that Eq.(29-31) also apply to n = 1 (∆ω = 0), in the limit F/U → 0 and γ/F → 0.

b † b = n 2 , (30)
g (2) (0) = 1 - 1 n . (31)
Driving out of multiphotonic resonances. In the case where multiphotonic absorption processes are non resonant, the mean photon number is expected to be very low. By performing an expansion of the master equation in powers of F/∆ω and γ/∆ω one can show that the coefficients of the stationary density matrix (in Fock space) obey the following hierarchy: Re(ρ nm ) ∼ ( F ∆ω ) n+m and Im(ρ nm ) ∼ γ ∆ω ( F ∆ω ) n+m . F/∆ω and γ/∆ω are assumed to be both much smaller than one, but it is not necessary to impose F ≫ γ. The following treatment is still valid if F/∆ω and γ/∆ω are of the same order of magnitude.

Neglecting the probability of having 3 or more photons inside the cavity and keeping the lowest order in F/∆ω and γ/∆ω for the remaining coefficients, we obtain the following set of equations (the real part of a complex number z is denoted by z, its imaginary part by z):

ρ 10 = F ∆ω ρ 00 , (32) ρ10
= - γ 2∆ω ρ 10 , (33)
ρ 11 = - 2F γ ρ10 , (34)
ρ 20 = √ 2F/∆ω 2 -U/∆ω ρ 10 , (35) ρ20
= 1 2 -U/∆ω (- γ ∆ω ρ 20 + F ∆ω √ 2ρ 10 ), (36)
ρ 21 = F/∆ω 1 -U/∆ω ( √ 2ρ 11 -ρ 20 ), ( 37
) ρ21 = -1 1 -U/∆ω ( 3γ 2∆ω ρ 21 + F ∆ω ρ20 ), (38)
ρ 22 = - √ 2F γ ρ21 . ( 39
)
Setting ǫ = F/∆ω; η = γ/∆ω and u = U/∆ω, the stationary density matrix reads:

ρ =    1 ǫ(1 + iη 2 ) √ 2ǫ 2 2-u (1 + iη 2 4-u 2-u ) ǫ(1 -iη 2 ) ǫ 2 √ 2ǫ 3 2-u (1 + iη 2-u ) √ 2ǫ 2 2-u (1 -iη 2 4-u 2-u ) √ 2ǫ 3 2-u (1 -iη 2-u ) 2ǫ 4 (2-u) 2    , (40)
From this we can compute the mean photon density and the second-order autocorrelation function:

b † b = F ∆ω 2 , (41)
g (2) (0) = 4 (2 -U/∆ω) 2 . ( 42
)
As expected, when the system becomes linear, i.e. U → 0, the cavity is driven into a coherent state (g (2) (0) = 1). However, the on-site interaction induces large fluctuations in the photon statistics when the two-photon absorption process becomes resonant (U/∆ω = 2).

B. Coupled Cavities: Mean-Field Solution

Let us first go back to the two-photon resonance. The analytical expression for the bosonic coherence in this regime is:

b = F ∆ω (2ξ -1) + i γ 2F (ξ -1), (43)
and the mean-field self-consistent equation is obtained by replacing F with F -J b . Since ξ is also a function of F , this equation is difficult to solve analytically in its general form. For ξ → 0, however, γ/F ≪ 1 and the imaginary part can be neglected. We are left with the simple expression:

b = - F ∆ω . ( 44
)
The substitution F → F -J b then gives:

b = -F/∆ω 1 -J/∆ω . ( 45
)
This shows that the coupling between sites amounts to replacing F with F ′ = F 1-J/∆ω . In other words, the effective pump is enhanced by the coupling between cavities. As a result, the system is driven into the ξ = 0 state and will stay there as long as the approximation F ′ /∆ω ≪ 1 holds. Results for different values of ξ are presented on Fig. (2). When J/∆ω ∼ 1 the above treatment ceases to be valid because F ′ ∼ ∆ω, and the system enters another regime. Exact P -representation calculations show that the mean photon density starts to increase with J, while g (2) (0) goes to 1, thus indicating a crossover from a quantum state to a classical coherent one (see Fig. (2)). As we shall see in section IV , this idea is confirmed by the fact that the linear asymptotic behavior of b † b as a function of J visible on Fig. (2), corresponds to Gross-Pitaevski semi-classical predictions.

Once again, this can be extended to larger values of n. For the n-photon resonance in the limit F n γ∆ω n-1 ≫ 1, Eq.( 44) for a single cavity becomes (see Appendix):

b = -(n -1) F ∆ω . (46)
When the coupling between cavities in switched on, b is given by:

b = -(n -1)F/∆ω 1 -(n -1)J/∆ω , (47)
which means that the system will stay in the state Eq.( 29) until J/∆ω ∼ 1/(n -1). As in the case of two-photon resonance, exact P -representation calculations presented on Fig. (3) show a crossover to a classical coherent state. This crossover is the closest equivalent, in this driven dissipative system, of the equilibrium Mott insulator to superfluid phase transition.

Out of multiphotonic resonances, the coupling between cavities has a different effect. In this regime, the system is in a state described by Eq.(40) at J = 0. The bosonic coherence is then:

b = F ∆ω . ( 48
)
At finite J, it becomes:

b = F/∆ω 1 + J/∆ω , (49)
and the effective pump is given by F → F ′ = F 1+J/∆ω . Contrary to Eq.(47), the intensity of F ′ decreases with J. As a consequence, the system will remain in a state qualitatively similar to Eq.( 40) and no crossover occurs. (continuous blue line) and γ/∆ω = F 2 ∆ω 2 (red dotted-dashed line). The lines show P -representation calculations and the markers the results of Eq.( 45). The effect of coupling J is to drive the system into the state of Eq.( 28), until the critical coupling Jc = ∆ω is reached. (green dotted-dashed line). The mean photon density and g (2) (0) are plotted as a function of the tunneling amplitude J/∆ω for F/∆ω = 10 -2 and F n /(∆ω n-1 γ) = 10. The system stays in the state of Eq.( 29) until a critical coupling Jc = ∆ω/(n -1), after which it is driven towards a coherent state.

C. Relation to Steady-States at Higher Pumping and Dissipation

The results presented above shed light on the nature of the steady-states obtained at higher pumping and dissipation. First, the effect of multiphotonic resonances is visible on the bistability diagram: they are responsible for its peculiar lobe structure (see Fig. [START_REF] Carusotto | [END_REF] or ref. [31]). Moreover, the properties of the states of Eq.( 29) and the behavior showed on Fig. (2) and Fig. (3) is very similar to that of the 'high-density phase mentioned in section II. In this phase the light is antibunched and the photon density increases with increasing J.

As for the 'low-density' phase, it shows photon superbunching, near the two-photon resonance, which is well described by Eq.(40). Besides, the photon density in this phase is decreasing with J, as suggested by Eq.( 49).

Comparing the expression for the bosonic coherence in Eq.( 49) and Eq.( 47) we see that in both regimes, the effect of tunneling is directly related to the sign of its real part. Interestingly, this remains true for higher pumping and dissipation, although the mean-field self-consistent equation takes a more complicated form. The real part of b is negative in the 'high-density' phase and positive in the 'low-density' phase.

IV. THE GROSS-PITAEVSKI REGIME

We have seen in section III, in the case of the multiphotonic resonances, that as J/∆ω becomes large, the system is driven into an almost coherent steady-state with b † b ≫ 1. This indicates that it enters a semi-classical regime where correlation functions can be approximated by :

b †n b m ≃ b † n b m . ( 50
)
As a consequence, all these functions are determined by a single complex number, namely the bosonic coherence β = b . Besides, a general differential equation for correlation functions can be readily obtained from Eq.( 2). Its most general expression in the context of mean-field therory is the following:

∂ t b †n b m = [b †n b m , H mf ] - iγ 2 (n + m) b †n b m . (51)
In the particular case of β = b and under the assumption of Eq.(50), the previous equation yields:

i∂ t β = (-∆ω -J - iγ 2 + U |β| 2 )β + F. ( 52
)
This equation is a single-mode version of the Gross-Pitaevski equation. Note that in this regime, the decoupling of neighboring sites amounts to a shift in the cavity frequency, ∆ω → ∆ω + J. The steady-state value for β is:

β = F ∆ω + J -U |β| 2 + iγ 2 , (53)
which gives a third order polynomial equation for the mean photon density n = |β| 2 :

n((∆ω + J -nU ) 2 + γ 2 4 ) = F 2 . ( 54
)
This equation explains the linear assyptotic behavior of n as a function of J/∆ω visible on Fig. (2). Indeed, when F, γ, → 0 and J → ∞, we find: Gross-Pitaevski approximation is also relevant at higher pumping and dissipation, especially when the coupling between sites and the number of photons are very high.

n ∼ J U , (55)
As it as been widely use in the theory of quantum fluids, whether with cold atoms or polaritons, it is fruitful to compare Gross-Pitaevski results with P -representation calculations presented in [31]. For example, Fig. [START_REF]The Physics of Semiconductor Microcavities[END_REF] shows that for large coupling between sites and in the 'high density' phase, Gross-Pitaevski approximation is sufficient to capture the behavior of the mean photon density as a function of the on-site interaction U .

A. Gross-Pitaevski Criterium for Bistability

Eq.( 54) was introduced in quantum optics as part of a semi-classical theory of optical bistability in a single nonlinear cavity [32]. In the present context, as mentioned in section II, mean-field theory predicts tunneling-induced bistability within a fully quantum framework. For some values of the parameters, Eq.( 54) has three real and positive roots. One of them corresponds to the 'low-density' phase (n ∼ 10 -2 ) and the two others to 'high-density phases'. Although only one 'high-density' phase was mentioned in our previous description of the mean-field phase diagram, a third solution was indeed found using generalized P -representation, but the corresponding phase proved to be always unstable.

Fig. [START_REF] Carusotto | [END_REF] shows the two bistability diagrams obtained respectively from Eq.( 54) and generalized Prepresentation. As expected, Gross-Pitaevski approximation is very good for small values of U , and predicts accurately the appearance of bistability in the lowerright corner of the diagram. It is less accurate when U becomes large and on the whole, bistability is "overestimated" by the Gross-Pitaevski criterium: monostable regions according to Eq.(54) (in orange on Fig.( 5)), are much smaller that the exact ones (in light blue). In particular it fails to predict the lobe structure that is visible on the P -representation diagram. These lobes stem from the n-photon resonances discussed in the previous section. Since a semi-classical approach does not take into account the quantized nature of the field, these resonances are washed-out in the Gross-Pitaevski diagram.

In the framework of Gross-Pitaevski approximation, the number of solutions is given by the sign of the discriminant of Eq.(54). A very good approximation for the critical value of U can be found by noticing that in the high-density phase, the photon density decreases with U . The critical value is then approximately the one for which the density is maximal. This yields:

U c1 ∆ω = γ 2 4F 2 (1 + J ∆ω ). (56)
In fact, this approximate expression corresponds to the first term in the expansion in powers of F/∆ω of the exact solution. A similar expansion for the other frontier in the diagram gives:

U c2 ∆ω = 4∆ω 2 27F 2 (1 + J/∆ω) 3 . ( 57
)
The expansion up to the next term is:

U c1 ∆ω = γ 2 4F 2 (1 + J ∆ω ) - γ 4 64∆ω 2 F 2 (1 + J/∆ω) , (58)
U c2 ∆ω = 4∆ω 2 27F 2 (1 + J/∆ω) 3 + γ 2 12F 2 (1 + J/∆ω). (59)

B. Bogoliubov Theory

As shown in our previous work, one can study fluctuations of the density matrix around mean-field by means of a extended Bogoliubov theory. The fluctuations are defined as follow:

ρ = i (ρ mf + δρ i ). ( 60
)
We also introduce the Fourier transform of the matrices

δρ i , δρ k = 1 √ N N i=1 e -ik•r i δρ i
In their most general formulations, the equations of evolution that stem from linearization around mean-field are:

i∂ t δρ k = L mf [δρ k ] + L k [δρ k ], (61)
where

L mf [δρ k ] = [H mf , δρ k ] - iγ 2 (2bδρ k b † -b † bδρ k -δρ k b † b).
(62) This operator is the usual Liouvillian for the effective single cavity problem. This term in Eq.( 61) is thus independent of k. Propagation effects arise from the second term:

L k [δρ k ] = -t k (Tr(bδρ k )[b † , ρ mf ] + Tr(b † δρ k )[b, ρ mf ]), ( 63
) with t k = J/z(cos k x a + cos k y a).
The situation is greatly simplified in the Gross-Pitaevski regime where the system is described by classical complex field. Fluctuations around the mean-field value β then obey the following equation:

i∂ t δβ k δβ * -k = -∆ω -t k + 2U |β| 2 -iγ/2 U β 2 U β * 2 ∆ω + t k -2U |β| 2 -iγ/2 δβ k δβ * -k . ( 64
)
This leads to a complex Bogoliubov spectrum:

ω ± (k) = ± (-∆ω -t k + 2U |β| 2 ) 2 -U 2 |β| 4 - iγ 2 . (65)
Dispersion relations extracted from Eq.(61) and Eq.( 65) are shown on Fig. (6). For small on-site repulsion and large tunneling amplitude ( U/∆ω = 0.5 and J/∆ω = 3, upper panels), the Gross-Pitaevski approximation give good quantitative results and the corresponding spectrum is included in the more general approached outlined in Eq.(61). As expected, it fails in the regime of strong correlations. Lower panels of Fig. (6) show the dispersion relations for U/∆ω = 2 and J/∆ω = 1. For these parameters, g (2) (0) = 0.69, proving that the hypothesis of a quasi coherent state underlying Gross-Pitaevski approximation scheme is not justified.

V. CONCLUSION

In this paper we have explored the driven-dissipative Bose-Hubbard model in the limit of weak pumping and weak dissipation and provided analytical results for the mean-field density matrix. In this regime, the drivendissipative model bear a formal resemblance to its equilibrium counterpart. However, since photons, unlike, atoms have to be injected by the pump laser, the mean photon number inside the cavity can be large only when multiphotonic absorption processes become resonant. In the case of n-photon resonance we have shown that if the intensity of the pump is sufficiently large when compared with the dissipation rate, a single cavity can be driven into a statistical mixture that has the same second-order correlation function g (2) (0) as a pure Fock state with n photons, and a mean photon number of n/2. At resonance, the effect of the coherent pump is enhanced by the coupling between sites, eventually leading to a crossover from these quantum states to classical coherent ones.This behavior is characteristic of the 'high-density' phase observed at higher pumping and dissipation, in the regime of tunneling-induced bistability. Outside of these multiphotonic resonance processes, the mean-photon density is much smaller than 1 and the effect of the pump is reduced by tunneling. Nevertheless, on-site interactions induce photon superbunching close to the two-photon resonance. This peculiar photon statistics is recovered at higher pumping and dissipation in the 'low-density' phase of the bistable region.

In addition, we have shown that the structure of the bistability diagram cannot be explained without a full quantum treatment of the single cavity. Indeed, a Gross-Piaevski semiclassical approach gives satisfactory results in the 'weakly-interacting' sector of the diagram, but fails in the strongly correlated regime. In particular, the size of the bistable region predicted by Gross-Pitaevski equation is considerably larger when compared with Prepresentation mean-field calculations.

VI. ACKNOWLEDGEMENT

We acknowledge support by the French ANR project QUANDYDE and by ERC grant 'CORPHO'. C.C. is member of Institut Universitaire de France. 61), the two red lines are the branches derived from Gross-Pitaevski equations. For these parameters, Gross-Pitaevski approximation is accurate. Lower panels: γ/∆ω = 0.2 ; F/∆ω = 0.4, U/∆ω = 2 and J/∆ω = 1 (monostable phase). Gross-Pitaevski approximation fails in the regime of strong correlations.

VII. APPENDIX

In this appendix we use the exact solution [32] of the single-cavity problem to prove Eq.( 29) and Eq.(46), holding for multiphotonic resonances in the limit γ ≪ F ≪ ∆ω.

For an isolated cavity, the density matrix of the stationary state is known analytically [32]:

ρ n,m = 1 √ n!m! -2F U n -2F * U m Γ(c)Γ(c * ) Γ(c + n)Γ(c * + m) F (c + n, c * + m, 4|F/U | 2 ) F (c, c * , 8|F/U | 2 ) , ( 66
) with c = 2(-∆ω -iγ/2) U . (67)
In Eq.(66) Γ(z) is the gamma special function which has poles at negative and zero integer values, whereas F is an hypergeometric series given by:

F (c, d, z) = ∞ k=0 Γ(c)Γ(d) Γ(c + k)Γ(d + k) z k k! . ( 68
)
In correspondence of a q-photon resonance, with q > 1, the constant c in Eq.( 67) is given by c = -(q -1)(1 + i γ 2∆ω

),

implying that c ≈ -(q -1) for γ ≪ ∆ω. As a result the quantities Γ(c + k), Γ(c * + k) are diverging for 0 ≤ k < q, implying that certain coefficients of the hypergeometric series (68) will actually diverge in the limit γ ≪ F ≪ ∆ω.

From the above consideration, the leading contribution in the two hypergeometric functions in Eq.(66) are given by F (c, c * , 2z) ≃ Γ(c)Γ(c * ) Γ(c + q)Γ(c * + q) (2z) q q! (70) and

F (c + k, c * + k, z) ≃ Γ(c + k)Γ(c * + k) Γ(c + q)Γ(c * + q) z q-k (q -k)! , (71)
respectively, where z = 4F 2 /U 2 . Substituting Eqs.( 70) and (71) into the general expression for the density matrix, Eq.(66), we find

ρ kk = 1 2 q k! q! (q -k)! , (72)
which corresponds to Eq.( 29) by replacing q with n.

It is also easy to see that off-diagonal terms will instead vanish in the same limit γ ≪ F ≪ ∆ω. Indeed, a similar analysis gives for m < n (the opposite case can be treated in the same way):

ρ nm = 1 √ n!m! -2F U n-m 1 2 q q! (q -m)!(n -m)! , (73)
which indeed vanishes for vanishing pump amplitude, F → 0. The expression for the bosonic coherence is obtained in a similar way. The general formula for b is:

b = F ∆ω + iγ/2 × F (1 + c, c * , 8| F U | 2 ) F (c, c * , 8| F U | 2 ) , (74)
and the leading term in the geometric function appearing in the numerator is:

F (c + 1, c * , 2z) ≃ Γ(c + 1)Γ(c * ) Γ(c + 1 + q)Γ(c * + q) (2z) q q! . (75)
Using Eq.(70), we find:

b = -(q -1) F ∆ω . (76)

FIG. 1 :

 1 photon density, g 2 (0)

F 2 10∆ω 2

 22 to F 5∆ω . Lines are the result of exact P -representation calculations while markers correspond to the simplified expression of Eq.(43). Right panel: mean photon density b † b (continuous blue line) and g(2) (0) (red dashed-dotted line), vs. ξ. Same conditions and conventions as in the other panel.

FIG. 2 :

 2 FIG.2:(Color online) Two-photon resonance. Mean photon density and g(2) (0) as a function of the tunneling amplitude J/∆ω for F/∆ω = 10 -2 . The values of γ/∆ω are F 2 10∆ω 2

4 FIG. 3 :

 43 FIG.3:(Color online) Quantum to classical crossover for n = 2 (continuous blue line), n = 3 (red dashed line), and n = 4 (green dotted-dashed line). The mean photon density and g(2) (0) are plotted as a function of the tunneling amplitude J/∆ω for F/∆ω = 10 -2 and F n /(∆ω n-1 γ) = 10. The system stays in the state of Eq.(29) until a critical coupling Jc = ∆ω/(n -1), after which it is driven towards a coherent state.

  which agrees with the results of Fig.(2) and Fig.(3).

FIG. 4 :FIG. 5 :

 45 FIG. 4: (Color online) Mean photon density versus on-site repulsion U in the high density phase, for J/∆ω = 3 ; F/∆ω = 0.4 ; γ/∆ω = 0.2. Red dots: P -representation calculations ; Blue crosses: Gross-Pitaevski approximation.

FIG. 6 :

 6 FIG.6:(Color online) Energy-momentum dispersion of elementary excitations. Upper panels: γ/∆ω = 0.2 ; F/∆ω = 0.4, U/∆ω = 0.5 and J/∆ω = 3 (high-density phase). Real and imaginary part of the low-energy branches (in units of γ) are plotted vs k. Blue and black lines depict branches obtained with Eq.(61), the two red lines are the branches derived from Gross-Pitaevski equations. For these parameters, Gross-Pitaevski approximation is accurate. Lower panels: γ/∆ω = 0.2 ; F/∆ω = 0.4, U/∆ω = 2 and J/∆ω = 1 (monostable phase). Gross-Pitaevski approximation fails in the regime of strong correlations.


Matthias Seibert

Julian Wörmann

Rémi Gribonval
email: gribonval@inria.fr

Martin Kleinsteuber
email: kleinsteuber@tum.de

Separable Cosparse Analysis Operator Learning

Keywords: Cosparse Analysis Model, Analysis Operator Learning, Sparse Coding, Separable Filters

published or not. The documents may come   L'archive ouverte pluridisciplinaire

INTRODUCTION

Sparse signal models enjoy great popularity for handling various signal processing tasks such as image denoising, super-resolution, compressive sensing, and others. In the last decade, most of the research activity has focused on the wellknown synthesis model, where the signal is synthesized via a linear combination of a few signal atoms taken from a dictionary. The synthesis model has a closely related counterpart called the analysis model where signals are characterized by their sparsity in a transformed domain. This transformation is achieved through an operator Ω ∈ R m×n , with m ≥ n. When applied to a signal x ∈ R n it yields the coefficient This work was partially supported by the German Research Foundation (DFG) under grant KL 2189/8-1, the DFG funded Cluster of Excellence CoTeSys, and by the European Research Council, PLEASE project (ERC- StG-2011-277906).

* These authors contributed equally to this work.

vector α ∈ R m defined as

α = Ωx. (1)
In the cosparse analysis model the signal information is encoded in the zero-entries of α. The number of zero elements in the coefficient vector, i.e. p = mα 0 , is called the cosparsity of the signal, and we refer to Ω as the cosparse analysis operator.

A prominent example for an analytic analysis operator is the finite difference operator in image processing. However, the advantage of the low computational complexity of such an analytically given transformation comes at the cost of a poor adaptation to specific signal classes of interest. It is now well-known that for a particular class, sparser signal representations and thus better reconstruction accuracies can be achieved if the transformation matrix is learned from a representative training set. We follow this approach in our work. Various learning algorithms that aim at learning a cosparse analysis operator have been proposed, cf. [START_REF] Yaghoobi | Analysis operator learning for overcomplete cosparse representations[END_REF][START_REF] Nam | Cosparse analysis modeling-uniqueness and algorithms[END_REF][START_REF] Rubinstein | Analysis K-SVD: A dictionary-learning algorithm for the analysis sparse model[END_REF][START_REF] Hawe | Analysis operator learning and its application to image reconstruction[END_REF][START_REF] Chen | Insights into analysis operator learning: From patch-based sparse models to higher order MRFs[END_REF]. In the following, we give a short summary of the algorithms proposed in [START_REF] Rubinstein | Analysis K-SVD: A dictionary-learning algorithm for the analysis sparse model[END_REF] and [START_REF] Hawe | Analysis operator learning and its application to image reconstruction[END_REF]. [START_REF] Rubinstein | Analysis K-SVD: A dictionary-learning algorithm for the analysis sparse model[END_REF] is of importance as, to our knowledge, the only available separable cosparse analysis operator learning scheme is based on the algorithm therein, while [START_REF] Hawe | Analysis operator learning and its application to image reconstruction[END_REF] is relevant to our cause since we extend its framework to the separable case. The algorithm proposed in [START_REF] Rubinstein | Analysis K-SVD: A dictionary-learning algorithm for the analysis sparse model[END_REF] is based on the famous K-SVD algorithm for dictionary learning. It consists of two alternating steps. First, a sparse coding step is performed, during which the current cosparse analysis operator is fixed and the best approximation of the noisy training signal with regard to the given cosparsity is computed. Second, each row of the cosparse analysis operator is updated using only information from the training signals orthogonal to the respective row. The general idea in [START_REF] Hawe | Analysis operator learning and its application to image reconstruction[END_REF] is to use a geometric optimization approach in order to fulfill the requirements demanded of the cosparse analysis operator. A smooth cost function is designed that, in addition to the sparsity promoting function, includes penalty functions that ensure full rank and incoherence of the learned operator. This cost function is minimized by using a geometric conjugate gradient method.

The cosparse analysis operators produced by these learning algorithms are structured only in the sense that they have unit row norm, while no further structure is enforced. These matrices allow for a high cosparsity in the analyzed signal. However, the size of the learned operators is ultimately restricted due to limitations in memory and computational power. Separable learning schemes tackle this problem by enforcing additional structure on the learned operators.

For the synthesis model the viability of separable approaches has been examined in [START_REF] Rigamonti | Learning separable filters[END_REF][START_REF] Hawe | Separable dictionary learning[END_REF]. To the best of the authors' knowledge, the only work combining a separable structure with the cosparse analysis model can be found in [START_REF] Qi | Two dimensional analysis sparse model[END_REF]. The authors suggest a learning scheme for two-dimensional signals, which is based on the Analysis K-SVD algorithm in [START_REF] Rubinstein | Analysis K-SVD: A dictionary-learning algorithm for the analysis sparse model[END_REF].

We propose a separable learning scheme based on the work in [START_REF] Hawe | Analysis operator learning and its application to image reconstruction[END_REF]. By introducing separability we are able to deal with data of arbitrarily high order while maintaining its inherent structure. We evaluate our proposed method on volumetric Magnetic Resonance Imaging (MRI) data.

Notation. Matrices are denoted as boldface capital letters such as X. Vectors are boldface small letters, i.e. x. Tensors are denoted as calligraphic letters X , whereas the (i 1 , i 2 , . . . , i N ) entry of a tensor with order N is denoted by x i1i2...i N . Entries of vectors and matrices are indexed accordingly. The i th row of a matrix is denoted by ω i,: . Some additional notation necessary for working with tensors will be introduced in the following section.

STRUCTURED COSPARSE ANALYSIS OPERATOR LEARNING

Tensor Operations

In order to be able to deal with multidimensional data, we have to introduce some techniques from the field of multilinear algebra. First, the so called n-mode product, introduced in [START_REF] Lathauwer | A multilinear singular value decomposition[END_REF], provides a way to apply transformations in the form of matrices to the separate modes of the tensor.

Definition 1. Given the tensor S ∈ R I1×I2×...×I N and the matrix U ∈ R Jn×In , the n-mode product is denoted by

S × n U
and results in a

I 1 × I 2 × . . . × I n-1 × J n × I n+1 × . . . × I N tensor.
The entries of this tensor are defined as

(S × n U) i1i2...in-1jnin+1...i N = In in=1 s i1i2...i N • u jnin (2)
for all j n = 1, . . . , J n .

To offer a better understanding of this concept, a visualization for a tensor of order 3 is provided in Figure 1. In addition to this product, we will also make use of the so-called n-mode matrix unfolding that enables us to one-to-one map a given tensor to a matrix of appropriate size. Since we will exclusively use unfoldings along the N th mode, we will only give the definition for this specific case.

= S I 1 I 2 I 3 A J 1 J 2 J 3 Ω (1) I 1 J 1 Ω (2) I 2 J 2 Ω (3) I 3 J 3 Fig. 1: The n-mode product A = S × 1 Ω (1) × 2 Ω (2) × 3 Ω (3)
with tensors A, S and matrices Ω (i) , i = 1, 2, 3.

I 1 I 2 I 3 I 2 I 3 I 1 → Fig. 2:
The 3-mode matrix unfolding of a 3-tensor.

Definition 2. The N -mode matrix unfolding of a tensor S ∈ R I1×...×I N is denoted by S (N ) ∈ R I N ×( j =N Ij ) . Given the tensor S ∈ R I1×...×I N , the unfolding matrix S (N ) contains the element s i1i2...i N in row i N and column

(i 1 -1)I 2 I 3 . . . I N -1 + (i 2 -1)I 3 I 4 . . . I N -1 + . . . + i N -1 .
For a general definition of this mapping we refer the interested reader to [START_REF] Lathauwer | A multilinear singular value decomposition[END_REF]. Figure 2 illustrates the unfolding of a third order tensor.

The unfolding makes it possible to rewrite the n-mode product of a tensor equivalently as a matrix-vector product.

A = S × 1 Ω (1) . . . × N Ω (N ) ⇔ A (N ) = Ω (N ) • S (N ) • Ω (1) ⊗ . . . ⊗ Ω (N -1) ⊤ ⇔ vec(A (N ) ) = Ω (1) ⊗ . . . ⊗ Ω (N ) • vec(S (N ) ) (3)
Here, vec denotes the operator that stacks the columns of a matrix below one another and ⊗ denotes the Kronecker product of matrices. The first reshaping is obtained by unfolding A and S along the N th mode, while the second follows from the definition of the Kronecker product.

We refer to a large, structured cosparse analysis operator Ω that can be represented as a concatenation of smaller operators {Ω (1) , . . . , Ω (N ) } as a separable operator.

Derivation of the Cost Function

With the introduced tensor operations we are able to define a way to apply a cosparse analysis operator to a multidimensional signal S recorded as an N -tensor. Let {S i }, i = 1, . . . , T be a set of training signals where each signal is a tensor of order N . The optimal separable cosparse analysis operator for this set is learned by solving the learning problem arg min

Ω (i) ∈Ci 1 T T j=1 g S j × 1 Ω (1) . . . × N Ω (N ) 2 (4)
with the sparsity promoting function

g(A) := k log 1 + να 2 k . (5)
The parameter ν > 0 serves as a weighting factor to control the sparsity, and the sum is taken over all entries α k of the tensor A. There are two things to point out in the design of Equation ( 4). First, in the proposed optimization problem we minimize the square of g(•). As shown in [START_REF] Hawe | Analysis operator learning and its application to image reconstruction[END_REF] this leads to a balanced minimization of expectation and variance of the samples' sparsity. Second, the minimization is performed over a constraint set C i which prevents the trivial solution of the given problem. The role of this constraint set will be discussed in the following.

As shown in (3), the n-mode product which is used in the cost function can be rewritten as a matrix-vector-product. The Kronecker product of the small operators that appears in (3) can be interpreted as a large cosparse analysis operator with additional structure which is applied to the vectorized signal. The algorithm we propose is based on the work by the authors in [START_REF] Hawe | Analysis operator learning and its application to image reconstruction[END_REF] and extended to a separable case. In [START_REF] Hawe | Analysis operator learning and its application to image reconstruction[END_REF], the authors enforce certain properties on their unstructured operator Ω via regularizers. They demand that (i) The rows of Ω have unit Euclidean norm.

(ii) The operator Ω has full rank, i.e. it has the maximal number of linear independent rows.

(iii) The rows of the operator Ω are not trivially linearly dependent, i.e. ω k,: = ±ω l,: for k = l.

As stated in [START_REF] Hawe | Analysis operator learning and its application to image reconstruction[END_REF] condition (i) is realized by the oblique manifold

OB(m, n) = Ω ∈ R m×n : (ΩΩ ⊤ ) ii = 1, i = 1, . . . , m .
In the following, we show that for a structured cosparse analysis operator enforcing these constraints on the operator components is equivalent to enforcing them on the large operator.

It should be pointed out that the cosparse analysis operators for the separate modes are actually the transposed of an element of the oblique manifold of appropriate size. It follows from the definition of the Kronecker product that if the rows of Ω (1) and Ω (2) have unit norm, then the rows of Ω (1) ⊗ Ω (2) have unit norm as well. Furthermore, the rank of the matrix resulting from Kronecker product Ω (1) ⊗ Ω (2) is the product of the rank of Ω (1) and Ω (2) . Finally, if the operators Ω (1) , Ω (2) do not exhibit trivially linearly dependent rows, then neither does Ω (1) ⊗ Ω (2) . By induction the above mentioned properties can be extended to the Kronecker product of finitely many operators. Hence, the properties demanded in [START_REF] Hawe | Analysis operator learning and its application to image reconstruction[END_REF] are fulfilled for Ω (1) ⊗ . . . ⊗ Ω (N ) if the conditions (i), (ii), and (iii) hold for each operator component Ω (i) , i = 1, . . . , N . To realize the full rank and the incoherence constraints, we adopt the regularizers

h Ω (i) = - 1 ni log(ni) log det 1 mi (Ω (i) ) ⊤ Ω (i) ,
r

Ω (i) = - k<l log 1 -(ω (i) k,: ) ⊤ (ω (i) l,: ) 2 .
The sum of these log-barrier functions is itself a logarithmic barrier functions which ensures that Ω (i) lies within the intersection of the two proposed constraint sets. With the corresponding weights κ and µ which control the importance of the corresponding constraint, we are able to formulate the optimization problem arg min

Ω (i) 1 T T j=1 g S j × 1 Ω (1) . . . × N Ω (N ) 2 + κ N j=1 h Ω (j) + µ N j=1
r Ω (j) subject to:

Ω (i) ∈ OB(m i , n i ), i = 1, . . . , N. (6)
This optimization task is tackled using a geometric conjugate gradient method based on the one introduced in [START_REF] Hawe | Analysis operator learning and its application to image reconstruction[END_REF]. Our algorithm is an extension in the sense that it is designed to work on the product of oblique manifolds. The line search required in the CG algorithm is performed using the nonmonotone line search strategy proposed in [START_REF] Zhang | A nonmonotone line search technique and its application to unconstrained optimization[END_REF].

MULTIDIMENSIONAL SIGNAL RECONSTRUCTION

Multidimensional signals can be found in many applications ranging from video, over hyper-spectral imaging to 3D models, medical images, and others. However, even within one signal class each dimension can represent a completely different physical meaning. While in hyper-spectral data the third dimension encodes the spectral component of a scene, volumetric data typically encodes spatial information. In this work we focus on three-dimensional data with a homogeneous physical meaning in all dimensions. An example for data of this kind is volumetric MRI data on which we evaluate our proposed scheme by conducting several signal reconstruction experiments.

For this purpose, we first extract a training set of T = 20 000 non-flat and normalized tensors of size S i ∈ R 5×5×5 from a synthetic, unbiased MRI brain atlas [START_REF] Fonov | Unbiased nonlinear average ageappropriate brain templates from birth to adulthood[END_REF][START_REF] Fonov | Unbiased average ageappropriate atlases for pediatric studies[END_REF], which is available online 1 . The parameter ν in ( 5) is set to ν = 1 000, while the weightings of the constraints read κ = 500 and µ = 0.5. After learning, we use the operators each of size Ω (i) ∈ R 6×5 , i = 1, 2, 3 to regularize the reconstruction for- mulated as an inverse problem. However, since the operators only act on local signal cubes, we first introduce the operator Π (r,c,k) (•) that extracts a small tensor of appropriate size centered at location (r, c, k) from the entire volume. Given the measurements y = Φ(V) + η ∈ R J1•J2•J3 arranged as a vector, where Φ(V) models the measurement process and η constitutes some noise, we finally state the reconstruction of the original volumetric signal V ∈ R I1×I2×I3 as the inverse problem

V ⋆ ∈ arg min V 1 2 Φ(V) -y 2 2 + λ M (r,c,k) g Π (r,c,k) (V) × 1 Ω (1) × 2 Ω (2) × 3 Ω (3) . (7)
In all our experiments, we employ the same sparsity promoting function g(•) as given in [START_REF] Chen | Insights into analysis operator learning: From patch-based sparse models to higher order MRFs[END_REF]. The parameter λ ∈ R + balances between the data fidelity and the regularization and M serves as a normalization factor. Optimization is performed via a conjugate gradient approach. Our algorithm is named MAOL (Multidimensional Analysis Operator Learning).

In our first experiment we evaluate the influence of preserving the multidimensional structure compared to the standard approach of vectorizing each tensor. For this purpose we conducted a simple denoising experiment. The second experiment illustrates the flexibility of our method which is applicable to a wide range of image processing tasks. A prominent linear inverse problem in medical imaging is the reconstruction of MRI images from partial Fourier measurements. Reducing the amount of measurements and thus reducing the scan time for MRI is highly desired for clinical throughput and convenience for patients. Based on the theory of Compressive Sensing (CS), many algorithms have been proposed in the literature so far that exploit the sparsity of this type of data. Typically, sparsity is measured in a fixed analytic transform domain like in the wavelet basis, while additional priors like total-variation further enhance the reconstruction quality. However, recently learning approaches have made their way into this field due to the fact that they adapt the structure of the signal more accurately. The original MRI data has been scaled such that its intensity values lie in the range [0, 255]. The performance is measured in terms of PSNR and the mean structural similarity (MSSIM) averaged over all slices in the whole volume. All experiments are implemented in MATLAB and executed on a standard desktop PC. [START_REF] Cocosco | Brainweb: Online interface to a 3D MRI simulated brain database[END_REF][START_REF] -S. Kwan | MRI simulation-based evaluation of image-processing and classification methods[END_REF] which has been artificially corrupted by AWGN. We chose the Analysis K-SVD (AKSVD) algorithm proposed in [START_REF] Rubinstein | Analysis K-SVD: A dictionary-learning algorithm for the analysis sparse model[END_REF] for comparison and learned a single operator Ω KSVD ∈ R 216×125 , i.e. the analyzed vector has the same dimension as ours. The operator Ω KSVD is directly learned from 20 000 noise contaminated tensors from the data that has to be reconstructed. Table 1 summarizes the results for two different choices of σ noise as well as the computation time for the operator learning t learn and the reconstruction t rec . We assume σ noise to be known and set λ = 200σ noise . While the performance of the two algorithms in the range of moderate AWGN is comparable in terms of PSNR and MSSIM, the separable structure of our operator enables a significant saving of computation time compared to the vectorizing approach.

Regarding our second experiment, the underlying imaging process in MRI consists in measuring Fourier coefficients of a typically transversal slice rather than directly acquiring pixel values. Thus, the measurement process is modeled via Φ(V) = Ψ u (V), where Ψ u (V) performs a multi-slice 2D undersampled Fourier transform resulting in a vector containing the stacked measurements of each slice. The right choice of the chosen coefficients has a direct impact on the reconstruction quality. Here, we use radial sampling as illustrated in Fig. 3a with an undersampling rate of 20%. The weighting of the regularizer reads λ = 1 500. We compare our results to a standard approach (RecPF) [START_REF] Yang | A fast alternating direction method for TVL1-L2 signal reconstruction from partial fourier data[END_REF] that minimizes an error term consisting of a least squares data fitting term combined with total variation and wavelet sparsity regularization. Furthermore, we give the results of the algorithm proposed in [START_REF] Ravishankar | MR image reconstruction from highly undersampled k-space data by dictionary learning[END_REF] (DLMRI) which utilizes an adaptively trained synthesis dictionary on each slice in the reconstruction process. Results on the same noiseless volume as mentioned above are summarized in Table 2. The RecPF algorithm has the fastest execution time, however its accuracy suffers from the fixed an- alytical transformations. Both learning approaches DLMRI and MAOL achieve significantly better results, while our approach has a runtime eight times faster than DLMRI. Fig. 3 shows 2D cross sections of the reconstructed volumes.

CONCLUSION

We provide a separable cosparse analysis operator learning framework for multidimensional signals called MAOL. Due to the separable nature, the computational complexity of both learning and reconstruction is reduced which makes it possible to train the cosparse analysis operator on larger signal patches compared to methods that generate unstructured operators. MAOL is designed using techniques from multilinear algebra and geometric optimization. The viability of our approach is demonstrated at the hand of experiments performed on volumetric MRI data.

Fig. 3 :

 3 Fig. 3: Sampling mask and reconstruction results for (b) inverse FFT after zero-filling missing Fourier coefficients, (c) the RecPF algorithm, and (d) our proposed separable approach MAOL.

Table 1 :

 1 MRI volume data reconstruction from measurements corrupted by additive white Gaussian noise with standard deviation σ noise . Achieved PSNR in decibels, t in minutes.

	Method σnoise	PSNR MSSIM	tlearn	trec
	AKSVD	5	38.41	0.968	980 840
		15	28.83	0.798
	MAOL	5	38.55	0.971	18	12
		15	30.64	0.851

Table 2 :

 2 Compressive Sensing MRI volume reconstruction from noiseless radial sampled measurements with an undersampling rate of 20%. Achieved PSNR in decibels.

	Method	PSNR MSSIM
	RecPF	26.06	0.788
	DLMRI	27.35	0.818
	MAOL	27.90	0.847
	4. NUMERICAL EXPERIMENTS
	Denoising aims at reconstructing a clean signal version from
	noise corrupted data, i.e. we have Φ(V) = vec(V) and η is ad-
	ditive white Gaussian noise (AWGN) with standard deviation
	σ noise . The volumetric MRI test data is a representative region
	of size 70 × 70 × 70 voxels from the T1 weighted BrainWeb
	phantom available online 2

http://www.bic.mni.mcgill.ca/ServicesAtlases/ICBM152NLin2009

http://brainweb.bic.mni.mcgill.ca/brainweb/