Angélique Drémeau

Patrick Héas

Cédric Herzet

SPARSE REPRESENTATIONS IN NESTED NON-LINEAR MODELS

Keywords: Non-linear sparse representation, dynamic programming, ℓ 0 -norm relaxation

Following recent contributions in non-linear sparse representations, this work focuses on a particular non-linear model, defined as the nested composition of functions. Recalling that most linear sparse representation algorithms can be straightforwardly extended to non-linear models, we emphasize that their performance highly relies on an efficient computation of the gradient of the objective function. In the particular case of interest, we propose to resort to a well-known technique from the theory of optimal control to estimate the gradient. This computation is then implemented into the optimization procedure proposed by Candès et al., leading to a non-linear extension of it.

INTRODUCTION

The sparse model assumes that a signal can be represented, exactly or approximatively, by a number of elements much smaller than its dimension. Exploited for more than twenty years, this model has proved to be a good prior for many types of signals in a variety of domains including audio [START_REF] Févotte | Sparse regression with structured priors: application to audio denoising[END_REF] and image [START_REF] Figueras I Ventura | Low rate and scalable image coding with redundant representations[END_REF] processing and is at the heart of the recent compressive-sensing paradigm [START_REF] Donoho | Compressed sensing[END_REF].

Standard sparse representation procedures usually focus on linear observation models, that is

y = Hx + n, (1)
where y is the observed data, say of dimension N , x is assumed to be sparse (i.e. contain very few non-zero elements) in a larger space of dimension M ≥ N and n stands for some observation noise. Recovering x from y then requires to solve a problem of the form (or some variants thereof):

x = argmin x x 0 subject to y -Hx 2 2 ≤ ǫ, (2)
where . 0 denotes the ℓ 0 pseudo-norm which counts the number of non-zero elements in its argument and ǫ > 0 controls the reconstruction error.

In practice, the linear observation model may be poorly adapted to many situations. As an example, we can mention the compressive phase retrieval problem, which aims at recovering a sparse signal from the knowledge of the amplitudes of some complex linear measurements (see e.g., [START_REF] Ohlsson | Compressive phase retrieval from squared output measurements via semidefinite programming[END_REF][START_REF] Schechtman | Gespar: Efficient phase retrieval of sparse signals[END_REF]). Hence, recent contributions have addressed the problem of exploiting sparse priors with non-linear observation models, that is

y = h(x) + n, (3)
for some non-linear observation operator h :

X M → Y N ,
where X, Y = R or C. Extending the approach followed in the linear case, these contributions propose also a generalization of the penalty function used in (2) of the form

x = argmin x x 0 subject to J(x) ≤ ǫ, (4)
where J(x) is some scalar function (e.g., J(x) = yh(x) 2 2 ) accounting for model [START_REF] Donoho | Compressed sensing[END_REF]. In this paper, we are interested in a particular case of nonlinearity, where the penalty function J(x) is defined as the nested composition of some functions. Formally, we write

J(x) = L l=1 J l • f l • . . . • f 1 (x), (5)
where {f l } L l=1 are some differentiable functions and • stands for the function-composition operator. This type of model is for instance of interest in the ubiquitous situations where one collects partial information on the state of a dynamical system whose initial condition admits a sparse decomposition (see section 3). In particular, we emphasize that results from optimal control [START_REF] Cannon | Theory of Optimal Control and Mathematical Programming[END_REF] can be exploited to provide a fast implementation of any gradient-based algorithm by taking benefit of the special structure of the non-linear model [START_REF] Schechtman | Gespar: Efficient phase retrieval of sparse signals[END_REF]. We propose then a practical implementation of this computation into the optimization procedure proposed in [START_REF] Candes | Enhancing sparsity by reweighted l1 minimization[END_REF].

For a sake of conciseness, we restrict our exposition to the real case (X = Y = R), but similar reasoning can be conducted in the complex case.

SPARSE REPRESENTATION AND GRADIENT EVALUATION

Sparsity-constrained linear models

In the literature addressing the standard sparse representation problem (2), many computationally-efficient procedures rely (often implicitly) on the fact that the evaluation of the gradient of J(x) = y -Hx 2 2 involves a low computational burden. More specifically, letting ∇ x [ ∂ ∂x1 , . . . , ∂ ∂x M ] T , we have that the gradient of J(x) evaluated at some point x ⋆ can be written as

∇ x J(x ⋆ ) = -2H T (y -Hx ⋆ ). (6)
We note that the evaluation of ∇ x J(x ⋆ ) only involves multiplications by the dictionary H and its transpose H T . In the case of general dictionaries, the complexity associated with the evaluation of the gradient thus scales as O(M N ). This constitutes one of the key ingredients of the success of several procedures efficiently tackling high-dimensional problems. Among others, we can mention the procedures based on a relaxation of the ℓ 0 pseudo-norm (e.g., FISTA [START_REF] Beck | A fast iterative shrinkagethresholding algorithm for linear inverse problems[END_REF], reweighted ℓ 1 norm [START_REF] Candes | Enhancing sparsity by reweighted l1 minimization[END_REF]), the family of thresholding algorithms (e.g., IHT [START_REF] Blumensath | Iterative thresholding for sparse approximations[END_REF]), or the greedy procedures (e.g., MP [START_REF] Mallat | Matching pursuits with timefrequency dictionaries[END_REF], OMP [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF], CoSaMP [START_REF] Needell | Cosamp: iterative signal recovery from incomplete and inaccurate samples[END_REF]) which sequentially update a support estimate by including the element x j leading to the highest local descent of J(x) (that is the element with the

largest partial derivative | ∂J(x (k) ) ∂xj |, where x (k) is the current estimate).

Sparsity-constrained non-linear models

It is noticeable that many algorithms mentioned above, when expressed in terms of the gradient of J(x), can straightforwardly be applied to non-linear sparse representation problems. Following this idea, the principles underlying IHT, MP, OMP and CoSaMP have for example been extended to the non-linear setting in [START_REF] Blumensath | Compressed sensing with nonlinear observations and related nonlinear optimization problems[END_REF], [START_REF] Beck | Sparsity constrained nonlinear optimization: optimality conditions and algorithms[END_REF], [START_REF] Blumensath | Gradient pursuit for non-linear sparse signal modelling[END_REF] and [START_REF] Bahmani | Greedy sparsity-constrained optimization[END_REF] respectively. Similarly, relaxed problems, more specifically based on a ℓ 1 relaxation of the ℓ 0 pseudo-norm, have been devised for constraints of particular form J(x) = ||y -h(x)|| 2 2 . In [START_REF] Ohlsson | Nonlinear basis pursuit[END_REF], the authors consider the noiseless case (ǫ = 0 in ( 4)) and propose to approximate h by its Taylor expansion reducing the non-linear term to a quadratic expression and allowing then the use of lifting techniques. A similar idea is applied in [START_REF] Ebtehaj | Variational data assimilation via sparse regularization[END_REF] to a non-linear operator h defined as a nested composition of functions. The initial optimization problem is thus reformulate as a quadratic programming problem through a first-order linearization.

The tractability of these extensions is however highly dependent on the efficient evaluation of the gradient ∇ x J(x). In this paper, we elaborate on this problem for the particular family of cost functions J(x) defined in [START_REF] Schechtman | Gespar: Efficient phase retrieval of sparse signals[END_REF]. Our exposition is based on the well-known theory of optimal control which traces back to the 70's (see e.g., [START_REF] Cannon | Theory of Optimal Control and Mathematical Programming[END_REF]).

Efficient gradient computation

Considering [START_REF] Schechtman | Gespar: Efficient phase retrieval of sparse signals[END_REF], the gradient of J(x) evaluated at some point x ⋆ is written as

∇ x J(x ⋆ ) = L l=1 ∇ x (J l • f l • . . . • f 1 ) (x ⋆ ), (7)
by virtue of the linearity of the operator ∇ x . Let us make two remarks. First, the composed function J l • f l • . . . • f 1 does not have any simple analytical expression in many situations; in such cases, we have therefore to resort to the chain rule of derivative for composed functions to evaluate its gradient. Second, the functions {f l } L l=1 appear in each term of J(x) in a structured manner and this fact should be taken into account in any efficient evaluation of ∇ x J. This is the goal of the procedure described hereafter.

Let us define, ∀l ∈ {1, . . . , L}, ∀x ∈ X M ,

s l f l • . . . • f 1 (x), (8)
and at the particular point of interest

x ⋆ : s ⋆ l f l •. . .•f 1 (x ⋆ ).
Clearly, with this definition, ∀l ∈ {1, . . . , L}, s l = f l (s l-1 ), and (5) evaluated at x ⋆ can be rewritten as

J(x ⋆ ) = L l=1 J l (s ⋆ l ). (9)
Therefore, using the chain rule of derivative, we obtain 1

∇ x J(x ⋆ ) = L l=1 ∇ x f l (s ⋆ l-1 ) T ∇ s l J l (s ⋆ l ),
with the convention s 0 = x 0 (resp. s ⋆ 0 = x ⋆ 0 ), and from the dependence between s l and s l-1 ,

∇ x f l (s ⋆ l-1 ) T = ∇ x f l (f l-1 (s ⋆ l-2 )) T , = ∇ x f l-1 (s ⋆ l-2 ) T ∇ s l-1 f l (s ⋆ l-1 ) T . ( 10
)
Applying this expression recursively, we finally have

∇ x J(x ⋆ ) = L l=1   l j=1 ∇ sj-1 f j (s ⋆ j-1 ) T   ∇ s l J l (s ⋆ l ). ( 11
)
On the one hand, we note that the latter expression does no longer involve the derivation of some composed functions but is exclusively based on the derivative of each component function. On the other hand, some care has to be taken in order to avoid unnecessary computational burden in the evaluation of [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF]. This gives rise to the following backwardforward procedure: 1 If f l (s ⋆ l-1 ) [f l,1 , ..., f l,N ] T , the operator ∇x applied to f l (s ⋆ l-1 ) T results in the M × N matrix whose (i, j)-th element is ∂f l,j ∂x i .

• The sequence {s ⋆ l } L l=1 is evaluated via the forward recursion

s ⋆ l = f l (s ⋆ l-1 ), (12) s
⋆ 0 = x ⋆ . (13)
• All multiplications by a same matrix ∇ s l-1 f l (s ⋆ l-1 ) T are gathered in one single operation. This is done through the backward recursion

p L = ∇ L-1 f L (s ⋆ L-1 ) T ∇ s L J L (s ⋆ L ), (14)
p l = ∇ l-1 f l (s ⋆ l-1 ) T (∇ s l J l (s ⋆ l ) + p l+1 ), (15)
leading finally to p 0 = ∇ x J(x ⋆ ). In that way, the multiplication by each matrix ∇ s l-1 f l (s ⋆ l-1 ) T is only performed once during the whole recursion. This backward-forward procedure is widely used in geophysical applications (e.g., [START_REF] Ba | Variational data assimilation for missing data interpolation in sst images[END_REF]). However, to the best of our knowledge, the explicit (and motivated) use of this technique into contexts of sparsity-constrained problems has never been considered. In particular, in [START_REF] Ebtehaj | Variational data assimilation via sparse regularization[END_REF] which focuses on a similar non-linear model, this efficient computation of the gradient is not proposed.

SPARSE REPRESENTATIONS IN DYNAMICAL MODELS

We emphasize that the structure of the cost function in ( 5) is well-suited to the characterization of dynamical systems with partial state information. Let us indeed consider a dynamical system characterized by a generic state evolution equation

s l = f l (s l-1 ) ∀l ∈ {2, . . . , L}. (16)
Assume moreover, that noisy partial observations of the states are collected at each time, that is

y l = g l (s l ) + n. (17)
A typical problem encountered in many domains of applications consists in recovering the sequence of {s l } L l=1 from the collected observation {y l } L l=1 . In the quite common case where the dimension of the collected data is inferior to the number of unknowns, one has to include an extra constraint on the sought vector in order to hope achieving a consistent estimation. Hereafter, we will assume that the initial state is sparse in some redundant dictionary H, that is

s 1 = Hx, (18)
for some sparse vector x. One possible formulation of the state estimation problem is therefore as follows

min x x 0 subject to L l=1 y l -g l (s l ) 2 2 ≤ ǫ, s l = f l (s l-1 ), s 1 = Hx. ( 19
)
Obviously, this problem can be reformulated as ( 4) with a cost function satisfying (5) by setting J l (z) y l -g l (z) 2 2 . The methodology described in this paper is therefore well-suited for gradient evaluation in this type of setup.

It is noticeable that many dynamical models typically evolve in high-dimensional spaces, leading in turn to sparserepresentation problems of very high dimensions. In such settings, an efficient evaluation of the gradient of J(x) turns out to be crucial for the tractable search of a solution of (4). In particular, any attempt to evaluate ∇ x J(x) by any standard numerical means (e.g., finite differences) is computationally intractable.

RESULTS: APPLICATION TO SQG MODEL

As an example of the methodology presented in this paper, we consider a non-linear sparse-representation problem in a particular geophysical application. More specifically, we focus on the characterization of the state of the ocean by exploiting the Surface Quasi-Geostropic (SQG) dynamical model [START_REF] Isern-Fontanet | Three-dimensional reconstruction of oceanic mesoscale currents from surface information[END_REF].

The SQG model assumes that some geophysical quantity s(u, t) (the so-called "buoyancy") obeys the following partial differential equation

∂ ∂t s(u, t) + (∇ ⊥ ∇ -1/2 s(u, t)) T ∇ u s(u, t) = 0, (20)
in which u ∈ R 2 and t ∈ R play the respective roles of spatial and temporal variables, and ∇ ⊥ ∇ -1/2 is a vectorial differential operator whose definition can be found in [START_REF] Isern-Fontanet | Three-dimensional reconstruction of oceanic mesoscale currents from surface information[END_REF]. In the sequel, we will consider a discretized version of (20), of the form of the state equation ( 16). This discretized model is built by means of a standard 4th-order Runge-Kutta numerical integration scheme [START_REF] Atkinson | An Introduction to Numerical Analysis[END_REF]. Satellites collect partial information {y l } L l=1 about the buoyancy at different time instants. We assume hereafter that each observation y l is related to the state of the system s l (but not directly to x) by a noisy linear observation model:

y l = G l s l + n, (21)
where

G l ∈ R N ×N .
The goal is then to recover the buoyancy from the lowdimensional information provided by the satellite by exploiting: i) the geophysical model [START_REF] Isern-Fontanet | Three-dimensional reconstruction of oceanic mesoscale currents from surface information[END_REF], nesting the buoyancy at different time instants; ii) the sparse decomposition of the initial condition s 1 in some redundant dictionary H [START_REF] Ebtehaj | Variational data assimilation via sparse regularization[END_REF].

At this point, the question is posed about the choice of the optimization procedure. As emphasized in section 2.2, providing an efficient evaluation of the gradient of the cost function, different well-known sparse optimization algorithms can be considered. Here, we propose to formulate the optimization problem as where λ > 0 and r(x) is some sparsity-enforcing regularizing function. In our simulation, we chose r(x) = m log(x m + ǫ), with ǫ = 10 -1 . Then our optimization procedure follows the majorization-minimization technique exposed in [START_REF] Candes | Enhancing sparsity by reweighted l1 minimization[END_REF]; at each iteration an upper bound on the goal function is constructed by majorizing r(x) by a weighted ℓ 1 norm. We look for the minimum of each of these majorizing functions by means of descent procedures involving the gradient of J(x) evaluated as presented in section 2.3.

min x L l=1 y l -G l (f l • . . . f 2 (Hx)) 2 2 + λ r(x), (22)
Particularized to the SQG model, the evaluation of the forward and backward recursions ( 12)-( 15) have a complexity of order O(M L). By comparison, using a finite-difference scheme to evaluate the gradient requires to run (at least) two forward recursions by element of x, leading to an overall complexity of O(M 2 L). This order of complexity thus precludes us from using this type of approach in moderate-tohigh dimensional problems.

The simulation setup considered in this paper is as follows. The state vectors s l are assumed to live in 256dimensional space. The initial condition is supposed to have a sparse decomposition in a dictionary H ∈ R 256×512 made up of sine and cosine functions. These dimensions have been chosen for the sake of running extensive simulation. We note that in practical SQG setups, the dimension of x is of the order of 512 2 or 1024 2 . The observations y l ∈ R 32 are collected at four different time instants and the observation matrices G l correspond to random subsampling operators. The ratio between the number of observations and the dimension of x is therefore equal to (32 × 4)/512 = 1/4. In Fig. 1, we represent the relative mean-square error (MSE) ||x -x|| 2 2 /||x|| 2 2 achieved by the minimization of ( 22) via the majorization-minimization procedure described above. As a point of comparison, we run the same algorithm on a linear sparse representation problem having the same problem di- mensions. For each data point, we average the performance over 50 trials. We can notice that the considered procedure can achieve an acceptable relative mean square error over a wide range of sparsity levels. We note also that the non-linear setup suffers from a reasonable degradation with respect to the linear setup. This tendency is confirmed in Fig. 2 which illustrates the probability of making at least one error on the support of the sought sparse vector.

CONCLUSION

In this paper, we address the problem of sparse representations in a non-linear setting. We emphasize that, unlike in the linear setting, the computation of the gradient of the cost function may be a bottleneck for the extension of standard estimation procedures. We show that this computation may be handled efficiently, by applying principles from the theory of optimal control, as long as the cost function satisfies some desirable structural property. Our derivations are illustrated on a particular example dealing with the estimation of the state of a geophysical system from partial observations.

Fig. 1 .

 1 Fig. 1. Relative MSE versus the number of non-zero coefficients K in the sparse vector.

Fig. 2 .

 2 Fig.2. Error probability versus the number of non-zero coefficients K in the sparse vector.

ACKNOWLEDGEMENTS

This work has been supported in part by the ERC under the European Union's 7th Framework Programme Grant Agreement 307087-SPARCS and by the CNRS/INSU through the LEFE funding program.


Angélique Drémeau

Patrick Héas

Cédric Herzet

Combining sparsity and dynamics: an efficient way

Most linear sparse representation algorithms can be straightforwardly extended to non-linear models. Their performance however, relies on an efficient computation of the gradient of the objective function. In this paper, we focus on a particular non-linear model, defined as the nested composition of functions and propose to resort to a well-known technique from the theory of optimal control to compute the gradient. As a proof of concept, this computation is then implemented into the optimization procedure proposed by Candès et al., and applied to a geophysical dynamical model.

Introduction

Recent contributions have addressed the problem of exploiting sparse priors with non-linear observation models, that is

y = h(x) + n, (1)
where h : R M → R N (with M ≥ N ) is a non-linear observation operator and n stands for an observation noise. Extending the approach followed in the linear case, these contributions propose also a generalization of the penalty function, leading to an optimization problem of the form (or some variants thereof)

x = argmin x x 0 subject to J(x) ≤ ǫ, (2)
where J(x) is some scalar function (e.g., J(x) = y -h(x) 2 2 ) accounting for discrepancies from model [START_REF] Blumensath | Iterative thresholding for sparse approximations[END_REF].

Noticing that many sparse representation algorithms dealing with linear observation models rely -implicitly or explicitlyon the computation of the gradient of the function J(x), nonlinear versions of them can be straightforwardly derived. Following this idea, the extensions of the well-known algorithms IHT [START_REF] Blumensath | Iterative thresholding for sparse approximations[END_REF], MP [START_REF] Mallat | Matching pursuits with timefrequency dictionaries[END_REF], OMP [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF] and CoSaMP [START_REF] Needell | Cosamp: iterative signal recovery from incomplete and inaccurate samples[END_REF] have thus been proposed, respectively in [START_REF] Blumensath | Compressed sensing with nonlinear observations and related nonlinear optimization problems[END_REF], [START_REF] Beck | Sparsity constrained nonlinear optimization: optimality conditions and algorithms[END_REF], [START_REF] Blumensath | Gradient pursuit for non-linear sparse signal modelling[END_REF] and [START_REF] Bahmani | Greedy sparsityconstrained optimization[END_REF] .

However, whereas in the linear case, the evaluation of the gradient of J(x) only involves multiplications by the dictionary and its transpose, its computational cost can be prohibitive in some non-linear cases. In this paper, we elaborate on this problem for the particular family of cost functions J(x) defined as the nested composition of some functions. Formally, we write

J(x) = L l=1 J l • f l • . . . • f 1 (x), (3)
where {f l } L l=1 are some differentiable functions and • stands for the function-composition operator. This type of model is for instance of interest in the ubiquitous situations where one collects partial information on the state of a dynamical system whose initial condition admits a sparse decomposition (see section 2.2). In particular, we emphasize that results from optimal control [START_REF] Cannon | Theory of Optimal Control and Mathematical Programming[END_REF] can be exploited to provide a fast implementation of any gradient-based algorithm by taking benefit of the special structure of the non-linear model [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF]. We propose then a practical implementation of this computation into the optimization procedure proposed in [START_REF] Candes | Enhancing sparsity by reweighted l1 minimization[END_REF].

Sparse Representations in Nested

Non-Linear Models

In this section, we elaborate on the efficient evaluation of the gradient of J(x) when structured as in (3). The methodology is then applied to a particular geophysical problem.

Efficient gradient computation

We use the following definitions and notations. Considering model (3), we set, ∀l ∈ {1, . . . , L}, ∀x ∈ R M ,

s l f l • . . . • f 1 (x). (4)
We thus have ∀l ∈ {1, . . . , L},

s l = f l (s l-1 ), (5)
with the convention s 0 = x. We also define the gradient operator as

∇ x [ ∂ ∂x 1 , . . . , ∂ ∂x M ] T , (6)
so that ∇ x applied to a vector z = [z 1 , ..., z N ] T results in the M × N matrix whose (i, j)-th element is ∂zj ∂xi . With these notations in mind, (3) evaluated at x ⋆ can be rewritten as

J(x ⋆ ) = L l=1 J l (s ⋆ l ), (7)
where s ⋆ l is defined as in (4) with x = x ⋆ . Therefore, using the chain rule of derivative, we obtain

∇ x J(x ⋆ ) = L l=1 ∇ x f l (s ⋆ l-1 ) T ∇ s l J l (s ⋆ l ),
and from the dependence between s l and s l-1 ,

∇ x f l (s ⋆ l-1 ) T = ∇ x f l (f l-1 (s ⋆ l-2 )) T , = ∇ x f l-1 (s ⋆ l-2 ) T ∇ s l-1 f l (s ⋆ l-1 ) T . ( 8
) Finally, applying this expression recursively, we have

∇ x J(x ⋆ ) = L l=1   l j=1 ∇ sj-1 f j (s ⋆ j-1 ) T   ∇ s l J l (s ⋆ l ). (9)
The latter expression is exclusively based on the derivative of each function component. Its evaluation can then be performed through the following forward-backward procedure:

• The sequence {s ⋆ l } L l=1 is evaluated via the forward recursion

s ⋆ l = f l (s ⋆ l-1 ), (10) s
⋆ 0 = x ⋆ . (11)
• All multiplications by a same matrix ∇ s l-1 f l (s ⋆ l-1 ) T are gathered in one single operation. This is done through the backward recursion

p L = ∇ L-1 f L (s ⋆ L-1 ) T ∇ s L J L (s ⋆ L ), (12)
p l = ∇ l-1 f l (s ⋆ l-1 ) T (∇ s l J l (s ⋆ l ) + p l+1 ), (13)
leading finally to p 0 = ∇ x J(x ⋆ ). In that way, the multiplication by each matrix ∇ s l-1 f l (s ⋆ l-1 ) T is only performed once during the whole recursion. This forward-backward procedure is widely used in geophysical applications (e.g., [START_REF] Ba | Variational data assimilation for missing data interpolation in sst images[END_REF]). However, to the best of our knowledge, the explicit (and motivated) use of this technique into contexts of sparsity-constrained problems has never been considered. In particular, in [START_REF] Ebtehaj | Variational data assimilation via sparse regularization[END_REF] which focuses on a similar nonlinear model, this efficient computation of the gradient is not proposed.

Application to super-resolution in SQG dynamical model

As a practical example of the proposed methodology, we focus on a super-resolution problem in a geophysical context, namely the high-resolution characterization of the state of the ocean by exploiting: i) the Surface Quasi-Geostropic (SQG) dynamical model [START_REF] Isern-Fontanet | Three-dimensional reconstruction of oceanic mesoscale currents from surface information[END_REF]; ii) a sparse prior on the initial condition of the dynamical model; iii) low-resolution satellite images. We assume that the SQG model, of the form of (5), rules the evolution of some state variable s l . The definition of f l depends on the considered numerical integration scheme (here a 4th-order Runge-Kunta method) but is not specified hereafter for conciseness. Moreover, as a prior knowledge, the initial state is supposed to be sparse in some redundant dictionary H,

s 1 = Hx, (14)
for some sparse vector x. We are interested in recovering the value of {s l } L l=1 from the observation of low-dimensional images {y l } L l=1 (e.g., collected by satellites), with

y l = G l s l + n, (15)
where G l is some known observation matrix and n is an unknown corrupting noise.

In order to solve this inverse problem, we address the following optimization problem

min x L l=1 y l -G l (f l • . . . f 2 (Hx)) 2 2 + λ r(x), (16)
where λ > 0 and r(x) = m log(x m + ǫ), ǫ = 10 -1 , is some sparsity-enforcing regularizing function. A solution to ( 16) is searched by using the majorization-minimization optimization technique exposed in [START_REF] Candes | Enhancing sparsity by reweighted l1 minimization[END_REF]; at each iteration an upper bound on the goal function is constructed by majorizing r(x) by a weighted ℓ 1 norm. We look for the minimum of each of these majorizing functions by means of descent procedures involving the gradient of J(x) (corresponding here to the first term in ( 16)) evaluated as presented in section 2.1.

Particularized to the SQG model, the evaluation of the forward-backward recursions (10)-( 13) have a complexity of order O(M L). By comparison, using a finite-difference scheme to evaluate the gradient requires to run (at least) two forward recursions by element of x, leading to an overall complexity of O(M 2 L). This order of complexity thus precludes us from using this type of approach in moderate-to-high dimensional problems.

The simulation setup considered in this paper is as follows. The state vectors s l are assumed to live in 256-dimensional space. The initial condition is supposed to have a sparse decomposition in a dictionary H ∈ R 256×512 made up of sine and cosine functions. The observations y l ∈ R 32 are collected at four different time instants and the observation matrices G l correspond to random subsampling operators. The ratio between the number of observations and the dimension of x is therefore equal to (32 × 4)/512 = 1/4. In Fig. 1, we represent the relative mean-square error (MSE) ||x -x|| 2 2 /||x|| 2 2 achieved by the minimization of (16) via the majorization-minimization procedure described above. As a point of comparison, we run the same algorithm on a linear sparse representation problem having the same problem dimensions (namely y = GHx where G is a rate-1/2 random subsampling matrix). For each data point, we average the performance over 50 trials.

We can notice that the considered procedure can achieve an acceptable relative mean square error over a wide range of sparsity levels. We note also that the non-linear setup suffers from a reasonable degradation with respect to the linear setup.

Conclusion

In this paper, we address the problem of sparse representations in a non-linear setting. While a high computational cost of the gradient of the goal function may prevent the use of standard estimation procedures, we show that it can be overcome by applying principles from the theory of optimal control, as long as the cost function satisfies some desirable structural property. Our derivations are illustrated on a particular example dealing with the estimation of the state of a geophysical system from partial observations.

Figure 1 :

 1 Figure 1: Relative MSE versus the number of non-zero coefficients K in the sparse vector.

This work has been supported in part by the ERC under the European Union's 7th Framework Programme Grant Agreement 307087-SPARCS and by the CNRS/INSU through the LEFE funding program.


Cédric Herzet

Charles Soussen

Enhanced Recovery Conditions for OMP/OLS by Exploiting both Coherence and Decay

Introduction

In this paper, we focus on two popular instances of greedy algorithms, namely orthogonal matching pursuit (OMP) [START_REF] Pati | Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition[END_REF] and orthogonal least squares 1 (OLS) [START_REF] Chen | Orthogonal least squares methods and their application to non-linear system identification[END_REF][START_REF] Natarajan | Sparse approximate solutions to linear systems[END_REF].

The suboptimal nature of OMP and OLS has led many researchers to study conditions under which these procedures succeed in recovering the true sparse vector. This question has been widely addressed for OMP in the recent years, including worst-case uniform [START_REF] Tropp | Greed is good: algorithmic results for sparse approximation[END_REF][START_REF] Davenport | Analysis of orthogonal matching pursuit using the restricted isometry property[END_REF] and probabilistic analyses [START_REF] Tropp | Signal recovery from random measurements via orthogonal matching pursuit[END_REF]. Although OLS has been known in the literature for a few decades (under different names [START_REF] Blumensath | On the difference between orthogonal matching pursuit and orthogonal least squares[END_REF]), exact recovery analysis for OLS have only appeared very recently, see [START_REF] Foucart | Stability and robustness of weak orthogonal matching pursuits[END_REF][START_REF] Soussen | Joint kstep analysis of orthogonal matching pursuit and orthogonal least squares[END_REF][START_REF] Herzet | Exact recovery conditions for sparse representations with partial support information[END_REF].

Most of the existing works deal with uniform guarantees: these conditions ensure the success of OMP/OLS for a given sparsity level (or a given support) irrespective of the magnitude of the non-zero coefficients. In contrast with these works, we derive new guarantees of success accounting for the decay of the non-zero elements of the sparse vector. Our conditions are expressed in terms of the mutual coherence of the dictionary µ and encompass, as particular cases, some well-known results of the literature. The proofs of the results are reported in our technical report [START_REF] Herzet | Relaxed recovery conditions for OMP/OLS by exploiting both coherence and decay[END_REF].

Context and Main Results

Let us assume that y ∈ R m is a (noisy) linear combination of k columns of A ∈ R m×n indexed by Q ⋆ , that is

y = Ax + w with x i = 0 ⇔ i ∈ Q ⋆ Card{Q ⋆ } = k (1)
where w ∈ R m denotes some noise vector and Card{•} stands for the cardinality operator. We assume that the columns a i of the dictionary are normalized:

a i 2 = 1 ∀i.
Let us first consider the noiseless case (w = 0). Theorem 1 provides sufficient conditions of success for OMP/OLS accounting for the possible decay of non-zero coefficients in x.

In our statement, we assume without loss of generality that

Q ⋆ = {1, 2, . . . , k}, (2)
and

|x 1 | ≥ |x 2 | ≥ . . . ≥ |x k | > 0. ( 3
) Theorem 1. If µ < 1 k , (4)
1 The OLS algorithm is also known as forward selection [START_REF] Miller | Subset Selection in Regression, Second Editon[END_REF], Order Recursive Matching Pursuit (ORMP) [START_REF] Cotter | Forward sequential algorithms for best basis selection[END_REF] and Optimized Orthogonal Matching Pursuit (OOMP) [START_REF] Rebollo-Neira | Optimized orthogonal matching pursuit approach[END_REF] in the literature. and

|x i | > 2µ(k -i) 1 -iµ |x i+1 | ∀i ∈ {1, . . . , k}, (5)
then Oxx selects atoms in Q ⋆ from noiseless data during the first k iterations.

In our technical report [START_REF] Herzet | Relaxed recovery conditions for OMP/OLS by exploiting both coherence and decay[END_REF], we show that this kind of result can also be extended to the characterization of the success of OMP/OLS in the noisy setting ( w 2 ≤ ǫ) or the partial recovery of the support by the latter algorithms. In particular, we show that a well-known result by Donoho et al. [START_REF] Donoho | Stable recovery of sparse overcomplete representations in the presence of noise[END_REF]Th. 5.1] can be further relaxed as

Theorem 2. If µ < 1 2k -1 , (6)
and

|x i | > 2ǫ 1 -(2k -i)µ ∀i ∈ {1, . . . , k}, (7)
then Oxx selects atoms in Q ⋆ during the first k iterations.

More specifically, we see that the condition by Donoho et al. in [START_REF] Donoho | Stable recovery of sparse overcomplete representations in the presence of noise[END_REF]Th. 5.1] is sufficient for [START_REF] Tropp | Greed is good: algorithmic results for sparse approximation[END_REF] to be satisfied and is therefore stronger than the conditions mentioned in Theorem 2. On the other hand, the two types of conditions become equivalent as soon as the nonzero elements of x have the same magnitude.


Chun-Hao Hao Huang

Edmond Boyer
email: edmond.boyer@inria.fr

Nassir Navab
email: navab@in.tum.de

Slobodan Ilic

Human Shape and Pose Tracking Using Keyframes

This paper considers human tracking in multi-view setups and investigates a robust strategy that learns online key poses to drive a shape tracking method. The interest arises in realistic dynamic scenes where occlusions or segmentation errors occur. The corrupted observations present missing data and outliers that deteriorate tracking results. We propose to use key poses of the tracked person as multiple reference models. In contrast to many existing approaches that rely on a single reference model, multiple templates represent a larger variability of human poses. They provide therefore better initial hypotheses when tracking with noisy data. Our approach identifies these reference models online as distinctive keyframes during tracking. The most suitable one is then chosen as the reference at each frame. In addition, taking advantage of the proximity between successive frames, an efficient outlier handling technique is proposed to prevent from associating the model to irrelevant outliers. The two strategies are successfully experimented with a surface deformation framework that recovers both the pose and the shape. Evaluations on existing datasets also demonstrate their benefits with respect to the state of the art.

Introduction

Marker-less human motion capture consists in tracking human shape and pose using visual information. This has become an important research area with many applications in motion analysis or digital content production. Perhaps the most widespread approach to solve this problem is to deform a pre-defined reference surface so as to fit data derived from image observations, e.g. silhouettes or 3D points. This model-based strategy has demonstrated a good success over the past few years [START_REF] Cagniart | Probabilistic deformable surface tracking from multiple videos[END_REF][START_REF] Gall | Motion capture using joint skeleton tracking and surface estimation[END_REF][START_REF] Huang | Robust human body shape and pose tracking[END_REF][START_REF] Liu | Markerless motion capture of interacting characters using multiview image segmentation[END_REF][START_REF] Straka | Simultaneous shape and pose adaption of articulated models using linear optimization[END_REF][START_REF] Vlasic | Articulated mesh animation from multi-view silhouettes[END_REF], because of its ability to enforce strong consistencies over time through the prior models of shape and deformation.

However, this strategy still relies on the assumption that image silhouette point cloud others ours

Figure 1. Our approach recovers the shape and the pose of the subject despite missing data (top row) and outliers (bottom row) whereas other approaches (top: [START_REF] Cagniart | Probabilistic deformable surface tracking from multiple videos[END_REF] and bottom: [START_REF] Straka | Skeletal graph based human pose estimation in real-time[END_REF] + [START_REF] Straka | Simultaneous shape and pose adaption of articulated models using linear optimization[END_REF]) fail.

image observations are complete and relevant, i.e. they do not describe another shape. In practice, it appears to be difficult to maintain such assumption when considering more realistic dynamic scenarios and with fewer constraints on the capture environment. As shown in Fig. 1, top row, background subtraction is often erroneous with the consequence of missing data, e.g. the missing arm. Another type of errors occurs when image observations describe a shape that is not in consideration, and can therefore mislead the surface as well, e.g. the chair in Fig. 1 bottom row. Our aim is to propose an alternative strategy that better handles such situations and hence contributes to the general objective of unconstrained human motion tracking in real environments.

Our framework is inspired by keyframe-based methods, e.g. [START_REF] Vacchetti | Fusing online and offline information for stable 3d tracking in real-time[END_REF] in camera/object tracking, and by non-sequential surface registration methods [START_REF] Budd | Global nonrigid alignment of surface sequences[END_REF][START_REF] Klaudiny | Towards optimal nonrigid surface tracking[END_REF]. In these works the tracking task is eased by reducing the discrepancy between the model and the input data to be matched. We exploit a similar idea that consists in having multiple reference shape models that can be fit to the observations. Numerous existing approaches rely on a one-fit-all strategy where a single reference model is deformed to fit to all observations. This strategy is likely to fail when observations describe a shape significantly different from the model, due to the presence of missing data or outliers. Instead, we propose to build a set of reference models called keyframes, which correspond to several representative shapes and poses that have been explored during tracking. They are identified online using mean-shift clustering and without the need for offline preprocessing. At each frame, the best keyframe is chosen as the reference model. We combine this strategy with a robust surface deformation method. Comparisons with the stateof-the-art confirm the advantages of this approach with real and inaccurate data.

This paper has several contributions. First we introduce the notion of keyframes in 3D human motion tracking, and further propose a keyframe-based tracking framework that updates the keyframe pool incrementally. Second, a new outlier rejection method is presented with the benefit of high integrability into previous probabilistic surface deformation framework. Both contributions increase the robustness and significantly limit the impact of missing data and outliers. To evaluate our method, we recorded new sequences that include static outliers. To the best of our knowledge, none of the current public dataset presents such feature.

Related work

Existing methods that track human poses and shapes generally express the problem as maximum a posteriori (MAP) estimation which involves a data term modeling the likelihood of the estimation and a regularization term modeling the adequacy to the prior information. Methods differ then by the input data and the assumed prior knowledge.

Data term

Data terms measure how well the model explains the observations. In general, silhouettes, point clouds, and photometric information are considered for this purpose.

Silhouettes. Many approaches deform the model such that the contour of the projected surface coincides with the contour of the observed silhouettes, e.g. [START_REF] Aguiar | Performance capture from sparse multi-view video[END_REF][START_REF] Gall | Motion capture using joint skeleton tracking and surface estimation[END_REF][START_REF] Straka | Simultaneous shape and pose adaption of articulated models using linear optimization[END_REF][START_REF] Vlasic | Articulated mesh animation from multi-view silhouettes[END_REF]. In [START_REF] Gall | Motion capture using joint skeleton tracking and surface estimation[END_REF], Gall et al. consider also photometric information to establish 3D-2D model-data correspondences. Later in [START_REF] Liu | Markerless motion capture of interacting characters using multiview image segmentation[END_REF], additional image segmentation information are used to differentiate multiple interacting subjects. In these works, silhouette overlap error is often regarded as a standard error measure, which is sensible only when silhouettes are accurate and fully describe shapes. Also note that distances in 2D images do not necessarily reflect distances in 3D and small errors along the silhouette contour can correspond to large distances in 3D. As a result, some authors, e.g. [START_REF] Cagniart | Probabilistic deformable surface tracking from multiple videos[END_REF][START_REF] Huang | Robust human body shape and pose tracking[END_REF] advocate for considering 3D points as input data.

Point clouds. Given a set of points reconstructed from multiple silhouettes, i.e. [START_REF] Franco | Exact polyhedral visual hulls[END_REF][START_REF] Furukawa | Accurate, dense, and robust multiview stereopsis[END_REF], some authors first estimate correspondences between the model and the 3D observations, and then deform the model accordingly. Although the reconstruction suppresses artifacts resulting from 2D noise, it also introduces new errors, such as missing body parts or fake geometry elements like ghost limbs. To robustly track in the presence of outliers, Huang et al. [START_REF] Huang | Robust human body shape and pose tracking[END_REF] train a linear support vector machine (SVM) that classifies the input data into different body parts. Outliers are then rejected based on the posteriors given by the SVM classifier. This approach depends heavily on the classifier, and is time-consuming since the SVM must be trained at each frame. In [START_REF] Cagniart | Probabilistic deformable surface tracking from multiple videos[END_REF], Cagniart et al. model outliers as an additional component of a Gaussian Mixture Model (GMM) equipped with a uniform distribution defined a priori. The adjustment of this distribution is however difficult and has a strong influence on the results. We propose a more robust outlier rejection that does not depend on user defined parameters.

Regularization term

Evolving a surface with discrete observations is ambiguous by nature and some prior information on the model is usually required. This information varies from generative spatial shape models to discriminative models that are learned from already known shapes and poses. Spatial shape models. Several works employ Laplacian coordinates [START_REF] Sorkine | Laplacian surface editing[END_REF] to preserve local shapes, e.g. [START_REF] Gall | Motion capture using joint skeleton tracking and surface estimation[END_REF][START_REF] Straka | Simultaneous shape and pose adaption of articulated models using linear optimization[END_REF], while others define a rigidity term that serves similar purposes, e.g. [START_REF] Cagniart | Free-form mesh tracking: a patch-based approach[END_REF]. Note that all these methods refer to a single static reference shape to constrain local deformations. This reference shape model is usually in rest pose and built prior to the tracking [START_REF] Cagniart | Probabilistic deformable surface tracking from multiple videos[END_REF][START_REF] Aguiar | Performance capture from sparse multi-view video[END_REF][START_REF] Vlasic | Articulated mesh animation from multi-view silhouettes[END_REF]. However, the observations can significantly deviate from the reference model along time. The shapes and poses that were already recovered during tracking can help in that respect, which motivates the multiplekeyframe tracking framework presented in this paper.

Learned deformation models.

A few works also make benefit from pre-collected information to help the tracking. They seek to learn the possible deformations in advance to regularize the results. In a non-sequential strategy, Budd et al. [START_REF] Budd | Global nonrigid alignment of surface sequences[END_REF] and Klaudiny et al. [START_REF] Klaudiny | Towards optimal nonrigid surface tracking[END_REF] assume that the complete input sequence is available beforehand and they find the best order to traverse it using a minimum spanning tree algorithm. Duveau et al. [START_REF] Duveau | Cage-based motion recovery using manifold learning[END_REF] propose a supervised learning strategy that regularizes the results based on the learned distribution in a latent parameter space. These methods require a pre-processing step either to build a shape-similarity tree from the input sequence [START_REF] Budd | Global nonrigid alignment of surface sequences[END_REF][START_REF] Klaudiny | Towards optimal nonrigid surface tracking[END_REF], or to learn a lowdimensional representation from the gathered motion training data [START_REF] Duveau | Cage-based motion recovery using manifold learning[END_REF][START_REF] Urtasun | Sparse probabilistic regression for activity-independent human pose inference[END_REF]. Our multiple-keyframe approach also exploits temporal information. As opposed to the mentioned approaches, we learn the reference models online in an unsupervised manner and do not require any preliminary step. Note anyway that our framework could also take advantage of already tracked sequences and generate the keyframe pool offline.

Overview

In this section we state our problem and give an overview of the proposed method. At every time frame t, a point cloud Y t = {y t i } i=1:Ny is reconstructed from silhouettes using EPVH [START_REF] Franco | Exact polyhedral visual hulls[END_REF]. y t i is a 6D vector that contains 3D spatial coordinates y and normal coordinates. The goal is to deform a reference model such that it fits the observations Y. Our model comprises a reference triangle surface and an intrinsic skeleton. We adopt the patch-based model proposed in [START_REF] Cagniart | Free-form mesh tracking: a patch-based approach[END_REF], where vertices are grouped into N p patches, as shown in Fig. 2(b). Our skeleton is a tree structure of N j nodes (3D joints) and the root is set at the pelvis, as shown in Fig. 2(c). The skeleton is rigged into the mesh using Pinocchio [START_REF] Baran | Automatic rigging and animation of 3d characters[END_REF], which gives the associations between vertices v and joints j. Deformations are parameterized with respect to: (i) the shape of the surface; (ii) the pose of the skeleton. The shape parameters Θ = {(R k , c k )} k=1:Np are the orientation and position pair for each patch k and encode the deformation of the reference mesh model. The pose parameters J = {x j } j=1:Nj are the 3D joint positions of the skeleton. Given the parametrization, the problem is formulated as the maximization of the joint probability distribution of the data and model:

max Θ,J P (Y, Θ, J ). (1)
This above distribution can be decomposed into P (Y|Θ) • P (J |Θ) • P (Θ), which represents respectively the likelihood of the shape given the observations, the probability of the pose given the shape, and the prior knowledge on shape deformations. Hence Eq. 1 can be rewritten as: where E r (Θ) = -ln P (Θ) is the rigidity energy in [START_REF] Cagniart | Probabilistic deformable surface tracking from multiple videos[END_REF] and E bone (Θ, J ) = -ln P (J |Θ) is the bone-binding energy in [START_REF] Huang | Robust human body shape and pose tracking[END_REF], both of which behave like regularization terms. The likelihood P (Y|Θ) is similar to [START_REF] Cagniart | Probabilistic deformable surface tracking from multiple videos[END_REF] and uses Gaussian Mixture Model (GMM) where every patch explains every observations y i according to:

min Θ,J [E r (Θ) + E bone (Θ, J ) -ln P (Y|Θ)], (2)
P (y i |Θ) = Np+1 k=1 Π k P (y i |z i = k, Θ). ( 3
)
z i is the latent variable for each y i : z i = k means that y i is generated by the mixture component associated with patch k. Π k = P (z i = k|y i , Θ) represents the probability that patch k explains observation y i . For each patch, when the closest vertex v k i with a compatible normal vector exists, the likelihood that y i is generated by the k-th component is modeled as a multivariate Gaussian with a mean located at the position x v k i of v k i and isotropic variance; otherwise the likelihood is a negligible number :

P (y i |z i = k, Θ) = N (y i |x v k i , σ 2 ) if v k i exists otherwise. ( 4
)
Solving Eq. 2 yields both the pose and the shape. In [START_REF] Cagniart | Probabilistic deformable surface tracking from multiple videos[END_REF], the tracking over a complete sequence is achieved by deforming the model on a frame-by-frame basis. That is, using (Θ t-1 , J t-1 ) as the initialization to solve Eq. 2 at frame t. To increase the robustness of this framework with respect to missing data and outliers better, we propose two methods that improve the deformation prior P (Θ) (Sec. 4) and the likelihood P (Y|Θ) (Sec. 5), respectively.

Multiple keyframe tracking framework

The rigidity energy E r (Θ) enforces neighboring patches to keep the original local configurations they have on the reference model, usually the shape at a given time t (e.g. t = 0). However, such local configurations do not always match with the current frame. This is particularly critical with missing data since patches without close observations tend to keep a possibly wrong prior reference configuration, as illustrated in Fig. 3.

Algorithm 1 Keyframe-based human motion tracking

1: F ← {0}, Ψ ← {(Θ 0 0 , J 0 )} 2: Overall shape parameters Θ t 0 ← {(I, c 0 k )} k=1:Np . 3: for t in timeFrames do 4:
Choose the reference model f ref based on Y t . 5:

Θ t-1 fref ← Θ t-1 0 * (Θ fref 0 ) -1 6:
With (Θ t-1 fref , J t-1 ) as initialization, solve Eq. 2 to obtain (Θ t fref , J t ).

7:

Θ t 0 ← Θ t fref * Θ fref 0 8:
if new keyframe detected then 9:

Update F and Ψ.

10:

end if 11: end for Therefore, to effectively handle missing data, we introduce a framework that exploits multiple reference models or keyframes. While already used in image tracking [START_REF] Vacchetti | Fusing online and offline information for stable 3d tracking in real-time[END_REF], keyframes have not yet, as far as we know, been applied to 3D human motion tracking problems. Multiple keyframes correspond to different instances of a shape that better represent the shape variability than a single pose at a given frame t. Our framework learns online keyframes and, at each frame, selects the best one as the reference model to be fitted to the observations. Let F = {f m } m=1:n f denote the keyframe pool where f m is the frame index and n f is the total keyframe number, and let Ψ = {(Θ fm 0 , J fm )} m=1:n f denote the corresponding parameter set. Our keyframebased tracking method is summarized in Alg. 1, where Θ f 0 corresponds to the accumulated rotation and translation from t = 0 to t = f , and (Θ f 0 ) -1 represents the inverse transformation. Allowing for different reference models enables different prior knowledge to be taken into account in the rigidity energy E r (Θ). Two crucial steps in Alg. 1 are:

• Line 4: how to select the best reference model from the keyframe pool?

• Line 8: when to add a new keyframe in the pool?

We tackle the former issue with shape dissimilarity and the latter with key pose detection. More details are elaborated in the following two subsections.

Keyframe detection

We first explain how new keyframes are added (i.e. Line 8 in Alg. 1). The essence of the multiple keyframe strategy lies in its ability to identify and record, during tracking, new local patch configurations. When the observed shape takes a pose at t that is very different from the reference pose, it is worth considering adding a new keyframe corresponding to frame t. Such analysis can be performed offline if knowledge on the shape poses is available prior to tracking. However, we consider here the more generic situation with little prior knowledge and where keyframes are detected online during tracking.

Pose descriptor. In order to incrementally identify distinctive poses, the previously-obtained skeleton poses J are explored. For each new pose J t , the pelvis of the skeleton is aligned to the global origin. The 3D coordinates of the remaining 14 joints are then concatenated to form a 42dimensional human body pose vector v. Aligning the pelvis to the origin cancels the global position and similar poses at different locations get similar representations. Here the skeleton is not rotated into a canonical direction and v still encodes orientation, which may not be a desirable attribute. We address this issue later. Similar poses yield vs that are close in R 42 whereas different poses correspond to vs that are distant. Using this descriptor, we cast the key pose detection as a clustering problem in R 42 . In this scenario, the number of clusters is supposed to be the output of the clustering technique and not a prior knowledge. Hence meanshift clustering naturally appears as a well adapted solution.

Mean-shift based key pose detector. Assume n f keyframes are already identified. Intuitively, there should be also n f clusters of poses. When a new pose vector v is obtained, mean-shift is performed and returns a number of clusters n c . If n c = n f , the number of clusters has not changed and we proceed to the next frame. If n c = n f + 1, it means that the shape pose has changed enough to justify a new cluster of poses. In general, poses in the same cluster also distribute closely in the time domain. Hence, starting from the current frame and going backward, the transition frame t where the new cluster starts is determined (e.g. t = 6 in Fig. 4(b)) and added as a new keyframe to F; its pose (Θ t 0 , J t) being added to Ψ. Fig. 4 illustrates this principle, where f last is the last element included in the keyframe set. The algorithm is sketched in Alg. 2, where Line 3 to 9 correspond to Line 8 to 10 in Alg. 1.

Note that, in general, t = t but t < t, which means that new keyframes are created online with some delay. For ex-Algorithm 2 Mean-shift-based key pose detector 1: Compute the bandwidth. Last keyframe f last ← 0 2: for each new incoming pose J t do 3:

De-pelvis the skeleton to obtain v t ∈ R 42 .

4:

De-pelvis all J f in Ψ and obtain V F = {v f }.

5:

Do mean-shift clustering on {v flast+1 . . . v t } ∪ V F .

6:

if the number of clusters = n f + 1 then Add transition frame t to F and (Θ t 0 , J t) to Ψ.

8:

f last ← t 9:
end if 10: end for ample, in Fig. 4(b), frame at t = 6 is detected as a keyframe with a 2 frame delay since in Fig. 4(a), v 6 is still in cluster no. 1. When a new keyframe at t is added, all the pose vectors before t, except the existing keyposes, are left out for further clustering (e.g. see Fig. 4(c)). This brings two advantages: first, the number of poses considered for clustering, i.e. (t-f last +n f ), is significantly reduced compared to the full set of poses; second, if a pose re-appears during tracking, the new collected v t is very likely to be clustered with existing keyframes, avoiding this way the occurrence of duplicated keyframes.

Mean-shift bandwidth. One concern with Alg. 2 is the bandwidth of mean-shift. A small bandwidth leads to many clusters while a large bandwidth gives few clusters. Since the intrinsic scale of the pose variation varies among sequences, an automatic way to determine this bandwidth is desirable. We achieve this using virtual pose vectors. Specifically, the de-pelvised skeleton model at t = 0 is rotated for 360 • with steps of 10 • , creating N s = 36 virtual pose vectors accordingly. Although they map to different points in R 42 , it is reasonable to cluster them together since they actually correspond to a single pose. We thus compute all Ns 2 pairwise distances and set the bandwidth as the half of their maximum. This way we ensure that they converge to the same mode with mean-shift. Recall that when a pose vector v is built from the estimated skeleton J , only the position is canceled but not the orientation. Using the above bandwidth, we expect pose vectors that differ only by a rotation to be clustered together hence canceling the rotation as well. More analysis on the influence of the bandwidth are presented in the experiment section.

Choosing the best keyframe

Given a new set of observations Y, the problem is now to determine the best keypose in F to be matched to these observations (Line 4 in Alg. 1). When the shape associated to such a keypose is close to the observed shape Y, it simplifies the estimation of the shape parameters and reduces the chances to fall in local minima. Therefore, we apply a shape similarity criterion to select the best keypose. This criterion uses shape histograms [START_REF] Ankerst | 3d shape histograms for similarity search and classification in spatial databases[END_REF] to describe 3D shapes and the L 2 distance between normalized histograms as the dissimilarity measure [START_REF] Huang | Shape similarity for 3d video sequences of people[END_REF]. The keyframe that presents the smallest dissimilarity with Y is chosen as the reference frame.

Patch-based outlier modeling

Besides missing data, sometimes point clouds contain false segmented foreground as the chair in Fig. 1 bottom row. In order to be robust to such outliers, care must be taken when designing the likelihood function P (Y|Θ). Note that the association of an observation y i to a patch, i.e. Eq. 4, applies only for z i = k ∈ [1, N p ] and that z i = N p + 1 is a special case introduced to model the outliers y i that are not explained by any patch. However, there is no physical outlier patch in the model to be associated to. In [START_REF] Cagniart | Probabilistic deformable surface tracking from multiple videos[END_REF], Cagniart et al. use a uniform distribution to model P (y i |z i = N p + 1, Θ) which basically assumes a certain proportion of the observations to be outliers and requires therefore some ad-hoc knowledge. Here we present a patch-based outlier modeling technique that takes into account spatial information and is based on the fact that the observations at frame t usually lie in the vicinity of the estimated surface at frame t -1.

Before modeling the general outlier event z i = N p + 1, we first consider the outlier event just for patch k, denoted as O k . The likelihood between y i and O k can be interpreted as how "bad" y i is explained by patch k. Since Eq. 4 expresses how well patch k explains y i under shape parameter Θ, we define the likelihood P (y i |O k , Θ) as:

P (y i |O k , Θ) ≡ 1 -P (y i |z i = k, Θ) = 1 -N (y i |x v k i , σ 2 ) if v k i exists 1 - otherwise. (5)
Eq. 5 also expresses, given shape parameter Θ and from the point of view of patch k, how likely y i is to be an outlier. Checking over all patches how poorly they explain y i , and assuming independence between O k , we can approximate the overall outlier likelihood as:

P (y i |z i = N p + 1, Θ) ≈ Np k=1 P (y i |O k , Θ). ( 6
)
From Eq. 6, we see that observations that are well explained by patches can not be outliers, where [START_REF] Cagniart | Probabilistic deformable surface tracking from multiple videos[END_REF] considers equal chances for each observation to be an outlier. 1. Sequences used for evaluation. We apply three different error measures, depending on the provided ground truth. A: silhouette overlap error. B: distances in R 3 between makers and associated vertices. C: distances in pixels with annotated joint positions. body parts because of the linear assumption of SVM. If we consider the SVM output for the outlier likelihood, as in Fig. 5(d), we cannot distinguish between chair and human observations. However, our patch-based outlier modeling is able to assign high values to chairs, as in Fig. 5(g), and to identify them by simple thresholding.

Optimization We follow the optimization framework in [START_REF] Cagniart | Probabilistic deformable surface tracking from multiple videos[END_REF][START_REF] Huang | Robust human body shape and pose tracking[END_REF] which alternates between estimating associations, i.e. Eq. 3, and solving for the model parameters, i.e. Eq. 2. These two steps correspond to the E-step and the M-step in the Expectation-Maximization framework [START_REF] Bishop | Pattern recognition and machine learning[END_REF]. The advantage of our outlier strategy is that it easily integrates into this method. In practice, outliers are not removed once and for all with hard thresholding, but Eq. 6 is substituted in Eq. 3. This means that outliers are estimated during the EM optimization, and that there is no need to add any other sophisticated learning-based method to improve outlier rejection.

Experiment results

The method was evaluated on 4 publicly available sequences as well as on 3 new sequences: WalkChair, Ham- [START_REF] Gall | Motion capture using joint skeleton tracking and surface estimation[END_REF] concern only their first stage results. Prev. are the results obtained when using the previous frame as the reference model. See Supplementary Material for more discussion.

merTable, and SideSit that contain static occlusions. These sequences were recorded with 9 cameras at 1000 × 1000 resolution. The occlusion objects are considered as foreground by the background subtraction and they remain in silhouettes and hence appear in the resulting point cloud observations, as shown in Fig. 8(c) and (d). We manually annotate the joint positions in 5 cameras to evaluate the poses of the skeleton. Due to the lack of realistic dynamic 3D surface ground truth, we use silhouette overlap error to evaluate the shape parameter estimation. If the occlusion objects are separated from the human body in silhouettes, we manually remove them and consider only the human parts as ground truth. These sequences are available at the 4D Repository 1 . To draw fair comparisons with other approaches, we do not refine the surfaces with silhouettes after tracking. These 7 sequences serve different purposes in the experiments and we summarize them in Table 1. Results are analyzed with respect to missing data and outliers, both qualitatively and quantitatively. In all the presented experiments, both the multiple keyframes strategy and the outlier rejection mechanism were used.

Robustness to missing data

Skirt and Dance demonstrate the effectiveness of the multiple keyframe strategy. The average silhouette overlap error for these sequences is shown in Table 2. In [START_REF] Gall | Motion capture using joint skeleton tracking and surface estimation[END_REF], Gall et al. refine the shape using silhouettes as a second stage of their method. Such refinement could fail if occlusion objects are close to the subject and appear in the silhouettes (e.g. Fig. 1 bottom row). We thus compare to their first stage results only. In WalkChair, missing data can be observed when the arms are too close to the torso. With the shape at t = 0 as the reference model result in Fig. 6(a obtained. The right arm stays in its configuration at t = 0 since not enough observations support the patches on the arm. Using another keyframe, as in Fig. 6(b), as the reference yields better result as shown in Fig. 6(c).

The influence of the mean-shift bandwidth. We report here on the influence of the mean-shift bandwidth on key pose detection and on the final tracking results. Tests were conducted on Skirt, HammerTable, and Fighting with varying bandwidths. Results on the keyframe numbers as well as on the corresponding errors are depicted in Fig. 7(a-c). The first two sequences were chosen for the repeating actions and the third one for its different numerical scale, see Fig. 7(d). In general, small bandwidths lead to more keyframes (green curves) and it appears that the errors (blue curves) also decline as the bandwidth decreases. However, small bandwidths have higher chances to identify tracking failures as keyframes and therefore accumulate errors. This explains why the error slightly rises when the bandwidth gets really small. We noticed that due to the different numerical scales, resulting from different recording setups, there is no fixed bandwidth that guarantees the best performance over the three sequences. This indicates that manually fixing the bandwidth is difficult. However, our strategy considers the numerical scale of the sequence and adjusts the bandwidth accordingly. This provides optimal or close to optimal performance (red dots). Note also that in Skirt, the subject raises up both arms and slowly rotates herself for a while, which leads to many similar poses only differing in orientations. As a result of our approach to cancel rotations, and the way we perform mean-shift with V F , duplicate key poses do not occur. Keyframe pools of all testing sequences are shown in the Supplementary Material.

Robustness to outliers

Static outliers. WalkChair, HammerTable and SideSit were used to demonstrate the robustness to outliers. Our approach was compared with [START_REF] Cagniart | Probabilistic deformable surface tracking from multiple videos[END_REF], [START_REF] Huang | Robust human body shape and pose tracking[END_REF] and the tracking framework proposed by Straka et al. ( [17]+ [START_REF] Straka | Simultaneous shape and pose adaption of articulated models using linear optimization[END_REF]). The silhouette overlap error as well as the discrepancies between the WlkChr. HmmrTbl. Cagniart [START_REF] Cagniart | Probabilistic deformable surface tracking from multiple videos[END_REF] 18482 fail Huang [START_REF] Huang | Robust human body shape and pose tracking[END_REF] 18063 fail Straka [START_REF] Straka | Skeletal graph based human pose estimation in real-time[END_REF][START_REF] Straka | Simultaneous shape and pose adaption of articulated models using linear optimization[END_REF] 12219 17285 ours 6803 3593 Table 3. Average silhouette overlap error of WalkChair (WlkChr.), and HammerTable (HmmrTbl.) from different approaches.

WlkChr.

HmmrTbl. SideSit Huang [START_REF] Huang | Robust human body shape and pose tracking[END_REF] 24.6 ± 10.7 fail 75 ± 40 Straka [START_REF] Straka | Skeletal graph based human pose estimation in real-time[END_REF][START_REF] Straka | Simultaneous shape and pose adaption of articulated models using linear optimization[END_REF] 20.6 ± 22.0 64.2 ± 53.9 84.4 ± 59.3 ours 15.9 ± 6.3 10.1 ± 3.0 19.3 ± 7.9 Table 4. Average joint 2D re-projection error in pixels of WalkChair (WlkChr.), HammerTable (HmmrTbl.), and SideSit. projected skeletons and the manually annotated joint positions in 5 cameras were measured. As reported in Table 3 and Table 4, our method attains consistently lowest errors in both pose and shape estimations. It is worth noting that in HammerTable, around 41% of the observations are not from the human subject, but we still get decent results (see Fig. 8(c) and (d)). We notice also that when the occlusions are closely touching the subject, it confuses the method in [START_REF] Straka | Skeletal graph based human pose estimation in real-time[END_REF] that deforms the skeleton model according to the observed skeletal graph in the point cloud, and thus the shape adaptation [START_REF] Straka | Simultaneous shape and pose adaption of articulated models using linear optimization[END_REF] cannot improve significantly the results. Please refer to the Supplementary Material and the accompanying video for more comparisons.

Dynamic outliers. Due to the lack of public datasets with dynamic outliers, we evaluate our approach with two multisubject sequences where we track only one subject and consider the others as outlier observations. For Basketball, the human subject is tracked against the ball observations. For Fighting, the subject with markers is tracked and the observations from the other subject are considered as outliers.

The results are shown in Fig. 8(a) and (b). First we con- firm that without the outlier rejection strategy presented in Sec. 5, the method in [START_REF] Cagniart | Probabilistic deformable surface tracking from multiple videos[END_REF] fails with this experimental setting, which is expected since it also fails to track against static outliers as in Table 3. Following the metric in [START_REF] Liu | Markerless motion capture of interacting characters using multiview image segmentation[END_REF], our approach attains 40.95 mm of average vertex position error with a standard deviation of 15.34 mm over 500 frames. Although the corresponding error in [START_REF] Liu | Markerless motion capture of interacting characters using multiview image segmentation[END_REF] is 29.61 mm and 25.50 mm, respectively, we would like to point out that their task is different from ours. Our objective is to robustly track in noisy environment whereas [START_REF] Liu | Markerless motion capture of interacting characters using multiview image segmentation[END_REF] simultaneously track two surfaces. Every observations is, in this case, associated to a patch or a vertex, and thus outliers are not to be considered. We observe anyway that, despite the presence of strong dynamic outliers (i.e. the ball and the second subject), our approach still provide reasonable results.

Figure 2 .

 2 Figure 2. (a-c): reference model. (d): observations. (a) is a surface with Nv vertices. In (b) the surface is decomposed into Np patches. (c) shows a rigged skeleton with Nj joints, and the associations between vertices and joints. (d) is a point cloud with Ny points. Typically, Nv ≈ 4K, Np ≈ 150, Nj = 15, and Ny ≈ 9K.

Figure 3 .

 3 Figure 3. Illustration of multi-keyframe advantages. (a) reference surface of Skirt [11] as the first keyframe. (b) second keyframe identified at t = 95. At t = 102, the left arm is missing in the observations, as in Fig. 1. Using (a) as the reference yields (c) while using (b) yields (d).

Figure 4 .

 4 Figure 4. Illustration of Algorithm 2. (a): time frame t = 6, F = {0}, and flast = 0. (b): t = 8, F becomes {0, 6}, and flast becomes 6. (c): t = 9, F = {0, 6}, and flast = 6. Only 5 instances (9 -6 + 2) are left for mean-shift.

Fig. 5

 5 illustrates this strategy. Fig. 5(a) and (b) shows the reference model represented in patches and body parts respectively. In [12], Huang et al. train a linear SVM on Fig. 5(b) obtained at t -1 to classify Y t into different rigid body parts, as shown in Fig. 5(c). Chair observations are classified as

Figure 5 .

 5 Figure 5. Comparison of outlier rejection in [12] (c-d) and ours (g-h). (a-b): reference surface colored in patches and body parts, respectively. (c): point cloud classified into body parts using the SVM trained on (b) from the previous frame, as suggested in [12]. (d): corresponding outlier likelihood from SVM. (e-f): two examples of Eq. 5. (g): point cloud colored based on Eq. 6. (h): points with outlier likelihood higher than 0.5 are colored in red.

Figure 6 .

 6 Figure 6. Results of frame t = 119 in WalkChair, with and without multiple keyframe strategy. Black dots are the observed point cloud where right arm gets merged into the torso. (a): estimated shape using surface at t = 0. (b): estimated shape at t = 32. (c): estimated shape using (b) as the reference surface.

Figure 7 .

 7 Fighting

Figure 8 .

 8 Figure 8. Results of (a) Fighting, (b) Basketball, (c) SideSit, and (d) HammerTable. Black dots represent the observed point clouds.

Table 2 .

 2 Average silhouette overlap error with different approaches. Image resolution: 1004 × 1004. Note that comparisons with

		Cagniart [6] Gall [11] Prev.	ours
	Skirt	7283	6900 7466 6715
	Dance	7881	7600	fail 6940

Acknowledgment. We would like to thank Matthias Straka from TU Graz for providing experimental results on their methods as well as valuable discussions.

Conclusion

We present an approach that captures human performances from multi-view video without markers. Considering realistic cases, we propose a multiple-keyframe-based tracking framework that uses mean-shift clustering to update a keyframe set online. A patch-based outlier modeling method is also presented to identify outliers more efficiently and effectively. Combining these two techniques into a surface deformation framework increases the robustness and enables the estimation of human shape and poses against missing data and outliers. The reliability of the proposed method is confirmed by the experiments on various public sequences as well as newly recorded sequences. Future di-rections include alleviating the requirement for background subtraction by considering photometric information.