Sid-Ahmed-Ali Touati
email: sid-ahmed.touati@prism.uvsq.fr

Denis Barthou
email: denis.barthou@prism.uvsq.fr

On the Decidability of Phase Ordering Problem in Optimizing Compilation

Keywords: D [3]: 4 Algorithms, Performance, Theory Phase Ordering, Parameters Space Exploration, Iterative Compilation, Library Generation, Optimizing Compilation

We are interested in the computing frontier around an essential question about compiler construction: having a program P and a set M of non parametric compiler optimization modules (called also phases), is it possible to find a sequence s of these phases such that the performance (execution time for instance) of the final generated program P ′ is "optimal" ? We prove in this article that this problem is undecidable in two general schemes of optimizing compilation: iterative compilation and library optimization/generation. Fortunately, we give some simplified cases when this problem becomes decidable, and we provide some algorithms (not necessary efficient) that can answer our main question.

Another essential question that we are interested in is parameters space exploration in optimizing compilation (tuning optimizing compilation parameters). In this case, we assume a fixed sequence of optimization, but each optimization phase is allowed to have a parameter. We try to figure out how to compute the best parameter values for all program transformations when the compilation sequence is given. We also prove that this general problem is undecidable and we provide some simplified decidable instances.

INTRODUCTION

The notion of an "optimal" program is sometimes ambiguous in optimizing compilation. Using an absolute definition, an optimal program P * means that there is no other equivalent program P faster than P * , whatever be the input data. This is equivalent to state that the optimal program should run as fast as the longest dependence chain in its trace. This notion of optimality cannot exist in practice: Schwiegelshohn et al showed in [START_REF] Schwiegelshohn | On Optimal Parallelization of Arbitrary Loops[END_REF] that there are loops with conditional jumps for which no semantically equivalent time-optimal program exists on parallel machines, even with speculative execution 1 . More precisely, they showed why it is impossible to write a program that is the fastest for any input data. This is because the presence of conditional jumps makes the program execution paths dependent on the input data, so it is not guaranteed that a program shown faster for a considered input data set (i.e., for a given execution path) remains the fastest for all possible input data. Furthermore, Schwiegelshohn et al convinced us that "optimal" codes for loops with branches (with arbitrary input data) requires the ability to express and execute a program with unbounded speculative window. Since any real speculative feature is limited in practice2 , it is "impossible" to write an optimal code for some loops with branches on real machines.

In our work, we define the program optimality according to the input data. So, we say that a program P * is optimal if there is not another equivalent program P faster than P * considering the same input data. Of course, the optimal program P * related to the considered input data I * must still execute correctly for any other input data, but not necessarily in the fastest speed of execution. In other term, we do not try to build efficient specialized programs, i.e., we should not generate programs that execute only for a certain input data set. Otherwise a simple program that only prints the results would be sufficient for fixed input data.

With this notion of optimality, we can ask the general question: how to build a compiler that generates an optimal program given an input data set ? Such question is very difficult to answer, since we are not able till now to enumerate all the possible automatic program rewriting methods in compilation (some are present in the literature, others have to be set up in the future). So, we first address in this work another similar question: given a finite set M of compiler optimization modules, how to build an automatic method to combine them in a finite sequence that produces an optimal program? We mean by compiler optimization module a program transformation that rewrites the original code. Unless they are encapsulated inside code optimization modules, we exclude program analysis passes since they do not modify the code.

This paper provides a formalism for some general questions about phase ordering. Our formal writing allows us to give preliminary answers from the computer science perspective about decidability (what we can really do by automatic computation) and undecidability (what we can never do by automatic computation). We will show that our answers are tightly correlated to the nature of the models (functions) used to predict or evaluate the programs performances. Note that we are not interested in the efficiency aspects of compilation and code optimization: we know that most of the code optimization problems are inherently NP-complete. Consequently, the proposed algorithms in this paper are not necessarily efficient, and are written for the purpose of demonstrating the decidability of some problems. Proposing efficient algorithms for decidable problems is another research aspect outside the current scope.

This paper is organized as follows. Section 2 gives a short overview about some phase ordering studies in the literature, as well as some performance prediction modeling. Section 3 defines a formal model for the phase ordering problem that allows us to prove some negative decidability results. Next, in Section 4, we show some general optimizing compilation scheme in which the phase ordering problem becomes decidable. Section 5 explores the problem of tuning optimizing compilation parameters with a compilation sequence. Finally, we present our future work before concluding.

BACKGROUND

The problem of phase ordering in optimizing compilation is coupled to the problem of performance modeling, since the performance prediction/estimation may guide the search process. The two following subsections present a quick overview of related work.

Performance Modeling and Prediction

Program performance modeling and estimation on a certain machine is an old (and is still) an important research topic aiming to guide code optimization. The simplest performance prediction formula is the linear function that computes the execution time of a sequential program on a simple von-Neumann machine: it is simply a linear function of the number of executed instructions. With the introduction of memory hierarchy, parallelism at many level (instructions, threads, process), branch prediction and speculation, performance prediction becomes more complex than a simple linear formula. The exact shape or the nature of such function and the parameters that it involves are two unknown problems until now. However, there exist many articles that try to define approximated performance prediction functions:

• Statistical Linear Regression Models: the parameters involved in the linear regression are usually chosen by the authors. Many program executions or simulation through multiple data sets allow to build statistics that compute the coefficients of the model [START_REF] Alexander | Performance Prediction for Loop Restructuring Optimization[END_REF][START_REF] Eeckhout | Quantifying the Impact of Input Data Sets on Program Behavior and its Applications[END_REF].

• Static Algorithmic Models: usually, such models are algorithmic analysis methods that try to predict a program performance [START_REF] Callahan | Estimating Interlock and Improving Balance for Pipelined Architectures[END_REF][START_REF] Mangione-Smith | Approaching a Machine-Application Bound in Delivered Performance on Scientifique Code[END_REF][START_REF] Wang | Precise Compile-Time Performance Prediction for Superscalar-Based Computers[END_REF][START_REF] Theobald | On the Limits of Program Parallelism and its Smoothability[END_REF]. For instance, the algorithm counts the instructions of a certain type, or makes a guess of the local instruction schedule, or analyzes data dependencies to predict the longest execution path, etc.

• Comparison Models: instead of predicting a precise performance metric, some studies provide models that compare two code versions and try to predict the fastest one [START_REF] Kennedy | Static Performance Estimation in a Parallelizing Compiler[END_REF][START_REF] Triantafyllis | Compiler Optimization-Space Exploration[END_REF].

Of course, the best and the most accurate performance prediction is the Turing machine itself, since it executes the program and hence we can directly measure the performance. This is what is usually used in iterative compilation and library generation for instance.

The main problem with performance prediction models is their aptitude to reflect the real performance on the real machine. As well explained by Rai Jain [START_REF] Jain | The Art of Computer Systems Performance Analysis : Techniques for Experimental Design, Measurement, Simulation, and Modeling[END_REF], the common mistake in statistical modeling is to trust a model simply because it plots a similar curve compared to the real plot (a proof by eyes !). Indeed, this sort of experimental validation is not correct from the statistical science theory, and there exist formal statistical methods [START_REF] Jain | The Art of Computer Systems Performance Analysis : Techniques for Experimental Design, Measurement, Simulation, and Modeling[END_REF] that check if a model fits the reality. Until now, we have not found any study that validates a program performance prediction model using such formal statistical methods.

Some Attempts in Phase Ordering

Finding the best order in optimizing compilation is an old problem. The most common case is the dependence between register allocation and instruction scheduling in instruction level parallelism processors as shown in [START_REF] Freudenberger | Phase Ordering of Register Allocation and Instruction Scheduling[END_REF]. Many other cases of inter-phase dependencies exist, but it is hard to analyze all the possible interactions [START_REF] Whitfield | An Approach for Exploring Code-Improving Transformations[END_REF].

Click and Cooper in [START_REF] Click | Combining Analyses, Combining Optimizations[END_REF] present a formal method that combines two compiler modules to build a super-module that produces better (faster) programs than if we apply each module separately. However, they do not succeed to generalize their framework of module combination, since they prove it for only two special cases, which are constant propagation and dead code elimination.

In [START_REF] Almagor | Finding effective compilation sequences[END_REF], the authors use exhaustive enumeration of possible compilation sequences (restricted to a limited sequence size). They try to find if any "best" compilation sequence emerges. The experimental results show that, unfortunately, there is not a winning compilation sequence. We think that this is because such compilation sequence depends not only on the compiled program, but also on the input data and the underlying executing machine.

In [START_REF] Veldhuizen | Guaranteed Optimization: Proving Nullspace Properties of Compilers[END_REF], the authors target a similar objective as in [START_REF] Click | Combining Analyses, Combining Optimizations[END_REF]. They succeed to produce super-modules that guarantee performance optimization. However, they combine two analysis passes followed by a unique program rewriting phase. In our work, we try to find the best combination of code optimization modules, excluding program analysis passes (unless they belong to the code transformation modules).

In [START_REF] Zhao | A Model-Based Framework: An Approach for Profit-driven Optimization[END_REF], the authors evaluate by using a performance model the different optimization sequences to apply to a given program. The model determines the profit of optimization sequences according to register resource and cache behavior. Optimizations consider only scalars and the same optimizations are applied whatever be the values of the inputs. In our article, we assume on the contrary that the optimization sequence should depend on the value of the input (in order to be able to speak about the optimality of a program).

Finally, there is the whole field of iterative compilation. In this research activity, looking for a good compilation sequence requires to compile the program multiple times iteratively, and at each iteration, a new code optimization sequence is used [START_REF] Cooper | Adaptive Optimizing Compilers for the 21st Century[END_REF][START_REF] Triantafyllis | Compiler Optimization-Space Exploration[END_REF] until a "good" solution is reached. In such frameworks, any kind of code optimization can be sequenced, the program performance may be predicted or accurately computed via execution or simulation. There exist other attempts that try to combine a sequence of high level loop transformations [START_REF] Cohen | A Polyhedral Approach to Ease the Composition of Program Transformations[END_REF][START_REF] Wolf | Combining Loop Transformations Considering Caches and Scheduling[END_REF]. As mentioned, such methods are devoted to regular high performance codes and only use loop transformation in the polyhedral model.

In this paper, we give a general formalism for the phase ordering problem and its multiple variants that incorporate the work presented in this section.

TOWARDS A THEORETICAL MODEL FOR PHASE ORDERING PROBLEM

In this section, we give our theoretical framework about the phase ordering problem. Let M be a finite set of program transformations. We would like to construct an algorithm A that has three inputs: a program P, an input data I and a desired execution time T for the transformed program. For each input program and its input data set, the algorithm A must compute a finite sequence s = mn • mn-1 • • • • • m0, mi ∈ M * of optimization modules 3 . The same transformation can appear multiple times in the sequence, as it occurs already in real compilers (for constant propagation/dead code elimination for instance). If s is applied to P, it must generate an optimal transformed program P * according to the input data I. Each optimization module mi ∈ M has a unique input which is the program to be rewritten, and has an output P ′ = mi(P). So, the final generated program P * is (mn

• mn-1 • • • • • m0)(P).
We must have a clear concept and definition of a program transformation module. Nowadays, many optimization techniques are complex toolboxes with many parameters. For instance, loop unrolling and loop blocking require a parameter which is the degree of unrolling or blocking. Until Section 5, we do not consider such parameters in our formal problem. We handle them by considering, for each program transformation, a finite set of parameter values, which is the case in practice. Therefore loop unrolling with an unrolling degree of 4 and loop unrolling with a degree of 8 are considered as two different optimizations. Given such finite set of parameter values per program transformation, we can define a new compilation module for each pair of program transformation and parameter value. So, for the remainder of the text (until Section 5), a program transformation can be considered as a module without any parameter except the program to be optimized.

In order to check that the execution time has reached some value T , we assume that there is a performance evaluation function t that allows to precisely evaluate or predict the execution time (or other performance metrics) of a program P according to the input data I. Let t(P, I) be the predicted execution time. Thus, t can predict the execution time of any transformed program P ′ = m(P) when applying a program transformation c. If we apply a sequence of program transformations, t is assumed to be able to predict the execution time of the final transformed program, i.e., t(P ′ , I) = t((mn • mn-1 • • • • • m0)(P), I). t can be either the measure of performance on the real machine, obtained through execution of the program with its inputs, a simulator or a performance model. In this article, we do not make the distinction between the three cases and assume that t is an arbitrary computable function. Next, we give a formal description of the phase ordering problem in optimizing compilation.

PB. 1 (PHASE-ORDERING). Let t be an arbitrary performance evaluation function. Let M be a finite set of program transformations. ∀T ∈ N an execution time (in processor clock cycles), ∀P a program, ∀I input data, does there exist a sequence s ∈ M * such that t(s(P), I) < T ? In other words, if we define the set:

St,M(P, I, T ) = {s ∈ M * |t(s(P), I) < T } is the set St,M(P, I, T ) empty? 3 • denotes the symbol of function combination (concatenation).

Textually, the phase ordering problem tries to determine for each program and input whether there exists or not a compilation sequence s which results in an execution time lower than a bound T .

If there is an algorithm that decides the phase ordering problem, then there is an algorithm that computes one sequence s such that t(s(P), I) < T , provided that t always terminates. Indeed, enumerating the code optimization sequences in lexicographic order always finds an admissible solution to Problem 1. Deciding the phase ordering problem is therefore the key for finding the best optimization sequence.

Decidability Results

In our problem formulation, we assume the following characteristics:

1. t is a computable function. t(P, I) terminates when P terminates on the input I. This definition is compatible with the fact that t can be the measured execution time on a real machine;

2. each program transformation m ∈ M is computable, always terminates and preserves the program semantics;

3. program P always terminates;

4. the final transformed program P ′ = s(P) executes at least one instruction, i.e., the final execution time is strictly positive.

The phase ordering problem corresponds to what occurs in a compiler: whatever the program and input be given by the user (if the compiler resorts to profiling), the compiler has to find a sequence of optimizations reaching some (not very well defined) performance threshold. Answering the question of the phase ordering problem as defined in Problem 1 depends on the performance prediction model t. Since the function (or its class) t is not defined, Problem 1 cannot be answered as it is, and requires to have another formulation that slightly changes its nature. We consider in this paper a modified version, where the function t is not known by the optimizer. The adequation between this assumption and the real optimizing problem is discussed after the problem statement.

PB. 2 (MODIFIED PHASE-ORDERING).

Let M be a finite set of program transformations. For any performance evaluation function t, ∀T ∈ N an execution time (in processor clock cycles), ∀P a program, ∀I input data, does there exist a sequence s ∈ M * such that t(s(P), I) < T ? In other words, if we define the set:

SM(t, P, I, T ) = {s ∈ M * |t(s(P), I) < T },
is the set SM(t, P, I, T ) empty?

This problem corresponds to the case where t is not an approximates model but is the real executing machine (the most precise model). Let us present the intuition behind this statement: a compiler always has an architecture model of the target machine (resource constraints, instruction set, general architecture, latencies of caches,. . . ). This model is assumed to be correct (meaning that the real machine conforms according to the model) but does not take into account all mechanisms of the hardware. Thus in theory, an infinite number of different machines fit into the model, and we must assume the real machine is any of them. As the architecture model is incomplete and performance also depends usually on non-modeled features (conflict misses, data alignment, operation bypasses,...), the performance evaluation model of the compiler is incorrect. This suggests that the performance evaluation function of the real machine can be any performance evaluation function, even if there is a partial architectural description of this machine. Consequently, Problem 2 corresponds to the case of the phase ordering problem when t is the most precise performance model which is the real executing machine (or simulator): the real machine measures the performance of its own executing program (for instance, by using its internal clock or its hardware performance counters).

In the following lemma, we assume an additional hypothesis: there exists a program that can be optimized into an infinite number of different programs. This necessarily requires that there is an infinite number of different optimization sequences. But this is not sufficient. As sequences of optimizations in M are considered as words made of letters from the alphabet M, the set of sequences is always infinite, even with only one optimization in M. For instance, fusion and loop distribution can be used repetitively to build sequences as long as desired. However, this infinite set of sequences will only generate a finite number of different optimized codes (ranging from all fusioned loops, to all distributed loops). If the total number of possible generated programs is bounded, then it may be possible to fully generate them in a bounded compilation time: it is therefore easy to check the performance of every generated program and to keep the best one. In our hypothesis, we assume that the set of all possible generated programs (generated using the distinct compilation sequences belonging to M * ) is infinite. One simple optimization such as strip-mine, applied many times to a loop with parametric bounds, generates as many different programs. Likewise, unrolling a loop with parametric bounds can be performed an infinite number of times. Note that the decidability of Problem 2 when the cardinality of M * is infinite while the set of distinct generated programs is finite remains an open problem.

LEMMA 1. Modified Phase-Ordering is an undecidable problem if there exists a program that can be optimized into an infinite number of different programs.

PROOF. The intuition of the proof is the following: at least one program P0 can be optimized into an infinite number of different programs. Assume it is possible to enumerate these different programs, and number them with an integer. For this program P0, we therefore simplify the problem: instead of looking for an optimization sequence, we are looking now for an integer, that makes the value of t lower than some given bound (we replace the program parameter of t with the integer). As t is an input of the problem and is any computable function (terminating on any input), undecidability of this subproblem is a classical result in computation theory.

The detailed proof works by reduction of the following problem: Problem (Empty Set) Given L a recursive enumerable language, is L empty ? This problem is known undecidable by application of Rice theorem [START_REF] Hopcroft | Introduction to Automata Theory, Languages, and Computation[END_REF].

We first reformulate the phase-ordering problem as a problem on computable functions instead of a problem on optimization sequences and evaluation functions. Sequences s of optimizations in M are considered as words made of letters from the alphabet M. There exists a program P0 such that optimization sequences applied to P0 generate an infinite number of different programs. Thus, there exists an algorithm aM,P 0 that, given an integer i, enumerates optimization sequences in lexicographical order, finds the first i sequences si that generate different optimized programs and outputs si(P0). aM,P 0 (i) = si(P0). This function is a bijective mapping between integers and optimized versions of P0. Now, we can define, for any evaluation function t, a function t ′ : t ′ (m, I) = t(aM,P 0 (m), I) for all integers m, I. t ′ is a computable function that always terminate. This defines a bijective mapping between evaluation functions and terminating computable functions with two parameters. Only one of these parameters (m) is not an input of the problem. We introduce one more in order to be able to perform the reduction: as there exists a bijective function h mapping any pair of integers to an integer, we define a mapping between evaluation functions and terminating computable functions with three parameters: t ′′ (m, n, I) = t(aM,P 0 (h(m, n)), I). The Modified Phase-ordering problem is equivalent to deciding whether the set:

S ′ (t ′′ , I, T ) = {(m, n)|t ′′ (m, n, I) < T }
is empty or not, where t ′′ is any computable function that always terminate. Now, we build the reduction. Consider a recursive enumerable language L. There exists a computable function g in {0, 1} such that m ∈ L ⇔ g(m) = 1. We build the following function tg,T , based on g:

• Input: m, n and I

• Perform at most n steps of the computation of g(m).

• If the computation has finished and g(m) = 1 then return T -1

• Else return T + 1.
This function tg,T always terminates and is computable. Moreover, there exists n such that tg,T (m, n, I) < T iff g(m) = 1. Now, given some integer T and an input I:

• If S ′ (tg,T , I, T ) is not empty, then there exist m and n such that g(m) stops before n computation steps and g(m) = 1. Therefore L is not empty.

• If S ′ (tg,T , I, T ) is empty, then for all integer values m and n, tg,T (m, n, I) ≥ T . According to the definition of tg,T , it means that for all m and n, either g does not stop before n steps of computation, or g stops and g(m) = 0. That implies that for all m integer, m ∈ L: L = ∅.

This shows that there is a reduction from the problem Empty Set to the Modified Phase Ordering problem. As the problem Empty Set is undecidable, the modified phase ordering problem is also undecidable.

We provide here a variation on the modified phase ordering problem that corresponds to the library optimization issue: program and (possibly) inputs are known at compile-time, but the optimizer has to adapt its sequence of optimization to the underlying architecture/compiler. This is what happens in Spiral [START_REF] Pschel | SPIRAL: Code Generation for DSP Transforms[END_REF] and FFTW [START_REF] Frigo | A Fast Fourier Transform Compiler[END_REF]. If the input is also part of the unknowns, the problem has the same difficulty.

PB. 3 (PHASE ORDERING FOR LIBRARY OPTIMIZATION).

Let M be a finite set of program transformations, P the program of a library function, I some input and T an execution time. For any performance evaluation function t, does there exist a sequence s ∈ M * such that t(s(P), I) < T ? In other words, if we define the set: SP,I,M,T (t) = {s ∈ M * |t(s(P), I) < T } is the set SP,I,M,T (t) empty?

The decidability results of Problem 3 are stronger than those of Problem 2: here the compiler knows the program, its inputs, the optimizations to play with and the performance bound to reach. However, there is still no algorithm to find out the best optimization sequence, if the optimizations may generate a infinite number of different program versions.

LEMMA 2. Phase Ordering for library optimization is undecidable if optimizations can generate an infinite number of different programs for the library functions.

PROOF. The proof is the same as the previous one, as the proof does not depend neither on the input I nor on the bound T .

The next section gives other formulations of the Phase-Ordering problem that do not alter the decidability results proved in this section.

Another formulation of Phase Ordering Problem

Instead of having a function that predicts the execution time, we can consider a function g that predicts the performance gain or speedup. g would be a function with three inputs: the input program P, the input data I and a transformation module m ∈ M. The performance prediction function g(P, I, m) computes the performance gain if we transform the program P to m(P) and by considering the same input data I. For a sequence s = (mn • mn-1 • • • • m0) ∈ M * we define the gain g(P, I, s) = g(P, I, m0)×g(m0(P), I, m1

)ו • •×g((mn-1•• • ••m0)(P), I, mn).
Note that, since the gains (and speedups) are fractions, the whole gain of the final generated program is the product of the partial intermediate gains. The ordering problem in this case becomes the problem of computing a compilation sequence that results in a maximal speedup, formally written as follows. This problem formulation is equivalent to the initial one that tries to optimize the execution time instead of speedup.

PB. 4 (MODIFIED PHASE-ORDERING WITH PERF. GAIN).

Let M be a finite set of program transformations. For any performance gain function g, ∀k ∈ N a performance gain, ∀P a program, ∀I input data, does there exist a sequence s ∈ M * such that g(P, I, s) ≥ k? In other words, if we define the set: SM(g, P, I, k) = {s ∈ M * |g(P, I, s) ≥ k}, is the set SM(g, P, I, k) empty?

We can easily see that Problem 2 is equivalent to Problem 4. This is because g and t are dependent each other by the following usual equation of performance gain: g(P, I, m) = t(P, I)t(m(P), I) t(P, I)

EXAMPLES OF DECIDABLE SIMPLI-FIED CASES

In this section we give some decidable instances of the phase ordering problem. As a first case, we define another formulation of the problem that introduces a monotonic cost function. This formulation models the real existing compilation approaches. As a second case, we model generative compilation and show that phase ordering is decidable in this case.

Models with Compilation Costs

In Section 3, the phase ordering problem is defined using a performance evaluation function. In this section, we add another function c that models a cost. Such cost may be the compilation time, the number of distinct compilation passes inside a compilation sequence, the length of a compilation sequence, distinct explored compilation sequences, etc. The cost function has two inputs: the program P and a transformation pass m. Thus, c(P, m) gives the cost of transforming the program P to P ′ = m(P ). Such cost does not depend on input data I. The phase ordering problem including the cost function becomes the problem of computing the best compilation sequence with a bounded cost.

PB. 5 (PHASE-ORDERING WITH DISCRETE COST FUNCTION).

Let t be performance evaluation function that predicts the execution time of any program P given input data I. Let M be a finite set of optimization modules. Let c(P, m) be a function that computes the cost of transforming the program P to P ′ = m(P ), m ∈ M. Does there exist an algorithm A that solves the following problem ? ∀T ∈ N an execution time (in processor clock cycles), ∀K ∈ N a compilation cost, ∀P a program, ∀I input data, compute A(P, I, T ) = s such that s = (mn • mn-1 • • • • m0) ∈ M * and t(s(P), I) < T with c(P, m0) + c(m0(P), m1)

+ • • • + c((mn-1 • • • • • m0)(P), mn) ≤ K.
We see in this section that if the cost function c is a strictly increasing function, then we can provide a recursive algorithm that solves Problem 5. First, we define the monotonic characteristics of the function c. We say that c is strictly increasing iff ∀m, m ′ ∈ M, c(P, m) < c(s(P), m ′ ) That is, applying a program transformation sequence mn•mn-1 • • •• m0 ∈ M * to a program P has always a higher integer cost than applying mn-1 • • • • m0 ∈ M * . Such assumption is true for the case of function costs such as compilation time 4 , number of compilation passes, etc. Each practical compiler uses an implicit cost function.

Building an algorithm that computes the best compiler optimization sequence given a strictly increasing cost function is an easy problem because we can use an exhaustive search of all possible compilation sequences with bounded cost. Algorithm 1 provides a trivial recursive method: it first looks for all possible compilation sequences under the considered cost, then it iterates over all these compilation sequences to check whether we could generate a program with the bounded execution time. Such process terminates because we are sure that the cumulative integer costs of the intermediate program transformations will certainly reach the limit K.

As illustration, the work presented in [START_REF] Almagor | Finding effective compilation sequences[END_REF] belongs to this family of decidable problems. Indeed, the authors compute all possible compilation phase sequences, but by restricting themselves to a given number of phases in each sequence. Such number is modeled in our framework as a cost function defined as follows:

∀P a program , c(P, s) =  1 + c(P, (mn-1 • • • • • m0)) ∀(mn • • • • • m0) ∈ M * 1 ∀m ∈ M
Textually it means that we associate to each compilation sequence the cost which is simply equal to the number of phases inside the compilation sequence. The authors in [START_REF] Almagor | Finding effective compilation sequences[END_REF] limit the number of

One-Pass Generative Compilers

Generative compilation is a subclass of iterative compilation. In such simplified classes of compilers, the code of an intermediate program is optimized and generated in a one pass traversal of the abstract syntax tree. Each program part is treated and translated to a final code without any possible backtracking in the code optimization process. For instance, we can take the case of a program given as an abstract syntax tree. A set of compilation phases treats each program part, i.e. each sub-tree, and generates a native code for such part. Another code optimization module can no longer re-optimize the already generated program part, since any optimization module in generative compilation takes as input only This program optimization process as described by Algorithm 2 computes the best compilation phase greedily. Adding backtracking changes complexity but the process still terminates. More generally, generative compilers making the assumption that sequences of best optimized codes are best optimized sequences fit the onepass generative compiler description. For example, the SPIRAL project in [START_REF] Pschel | SPIRAL: Code Generation for DSP Transforms[END_REF] is a generative compiler. It performs a local optimization to each node. SPIRAL optimizes FFT formulae, from the formula level, by trying different decomposition of large FFTs. Instead of a program, SPIRAL starts from a formula, and the optimizations considered are decomposition rules. From a formula tree, SPIRAL recursively applies a set of program transformations at each node, starting from the leaves, generates C code, executes it and measures its performance. Using dynamic programming strategy 5 , composition of best performing formulae are considered as best performing compositions.

As can be seen, finding a compilation sequence in generative compilation that produces the fastest program is a decidable problem (Algorithm 2). Since the size of intermediate representation forms decreases at each local application of program transformation, we are sure that the process of program optimization terminates when all intermediate forms have been transformed to native In other terms, the number of possible distinct passes on a program becomes finite and bounded as shown in Algorithm 2: for each node of the abstract syntax tree, we apply locally a single code optimization (we iterate over all possible code optimization modules and we pick up the one that produces the best performance according to the chosen performance model). Furthermore, no code optimization sequence is searched locally (only a single pass is applied). Thus, if the total number of nodes in the abstract syntax tree is equal to ñ, then the total number of applied compilation sequences does not exceed |M| × ñ.

Of course, the decidability of one-pass generative compilers does not prevent them from having potentially high complexity: each local code optimization may be exponential (if it tackles NP-complete problem for instance). The decidability result only proves that, if we have a high computation power, we know that we can compute the optimal code after a bounded compilation time (possibly high).

This first part of the article investigates the decidability problem of phase ordering in optimizing compilation. Figure 1 synthesizes a whole view of the different classes of the investigated problems with their decidability results. The largest class of the phase ordering problem that we consider, denoted by C1, assumes a finite set of program transformations with possible optimization parameters (to explore). If the performance prediction function is arbitrary, typically if it requires program execution or simulation, then this problem is undecidable. The second class of the phase ordering problem, denoted by C2 ⊂ C1, has the same hypothesis as C1 except that the optimization parameters are fixed. The problem is undecidable too. However, we have identified two decidable classes of phase ordering problem which are C3 and C4 explained as follows. The class C3 ⊂ C2 considers one-pass generative compilation ; the program is taken as an abstract syntax tree (AST), and code optimization applies a unique local code optimization module on each node of the AST. The class C4 ⊂ C2 takes the same assumption as C2 plus an additional constraint which is the presence of a cost model: if the cost model is a discrete increasing function, and if the cost of the code optimization is bounded, then C4 is a class of decidable phase ordering problem.

The next section investigates another essential question in optimizing compilation, which is parameters space exploration.

COMPILER OPTIMIZATION PARAME-TERS SPACE EXPLORATION

Nowadays, many compiler optimization methods are parametrized. For instance, loop unrolling requires an unrolling degree; loop blocking requires a blocking degree as well, etc. Actually, the complexity of phase ordering problem does not allow to explore jointly the the best sequence of the compilation steps and the best combinations of modules parameters. Usually, the community tries to find the "best" parameter combination when the compilation sequence is fixed. This section is devoted to study the decidability of such problem.

Towards a Theoretical Model

First, we suppose that we have s ∈ M * a given sequence of optimizing modules belonging to a finite set M. We assume that s is composed of n compilation sequences.

We associate for each optimization module mi ∈ M a unique integer parameter ki ∈ N. The set of all parameters is grouped inside a vector -→ k ∈ N n , such that the i th component of -→ k is the parameter ki of the mi, the i th module inside the considered sequence s. If the sequence s contains multiple instances of the same optimization module m, the parameter of each instance may have a distinct value from those of the other instances.

For a given program P, applying a program transformation module m ∈ M requires a parameter value. Then, we write the transformed program as

P ′ = m(P, - → k ).
As in the previous sections devoted to the phase ordering problem, we assume here the existence of a performance evaluation function t that predicts (or evaluates) the execution time of a program P having I as input data. We denote t(P, I) the predicted execution time. The formal problem of computing the best parameter values of a given set of program transformations in order to achieve the best performance can be written as follows.

PB. 6 (BEST-PARAMETERS). Let t be a function that predicts the execution time of any program P given input data I. Let M be a finite set of program transformations and s a particular optimization sequence. Does there exist an algorithm At,s that solves the following problem ? ∀T ∈ N an execution time (in processor clock cycles), ∀P a program, ∀I input data, At,s(P,

I, T ) = - → k such that t(s(P, - → k ), I) < T .
This general problem cannot be addressed as it is, since the answer depends on the shape of the function t. In this paper, we assume that the performance prediction function is built by an algorithm a, taking s and P as parameters. Moreover, we assume the performance function t = a(P, s) built by a takes -→ k and I as parameters and is a polynomial function. Therefore, the performance of a program P with input I and optimization parameters -→ k is a(P, s)(I, -→ k ). We discuss about the choice of a polynomial model after the statement of the problem. We want to decide whether there are some parameters for the optimization modules that make the desired performance bound reachable: PB. 7 (MODIFIED BEST-PARAMETERS). Let M be a finite set of program transformations and s a particular optimization sequence of M * . Let a be an algorithm that builds a polynomial performance prediction function, according to a program and an optimization sequence. For all programs P, for all inputs I and performance bound T , we define the set of parameters as:

Ps,t(P, I, T ) = { - → k |a(P, s)( - → k , I) < T }.
Is Ps,t(P, I, T ) empty ?

As noted earlier, choosing an appropriate performance model is a central decision to define whether Problem 6 is decidable or not. For instance, Problem 7 considers polynomial functions, which are a family of usual performance models (arbitrary linear regression models for instance). Even a simple static model of complexity counting assignements evaluates usual algorithms with polynomials (n 3 for a straightforward implementation of square matrixmatrix multiply for instance). With such a simple model, any polynomial can be generated. It is assumed that a realistic performance evaluation function would be as least as difficult as a polynomial function. Unfortunately, the following lemma shows that if t is an arbitrary polynomial function, then Problem 7 is undecidable.

The following lemma states that Problem 7 is undecidable if there are at least 9 integer optimization parameters. In our context, this requires 9 optimizations in the optimizing sequence. Note that this number is constant when considering the best parameters, and is not a parameter itself. This number is fairly low compared to the number of optimizations found in state-of-the-art compilers (such as gcc or icc for instance). Now, if t is a polynomial and there are less than 9 parameters (the user has switched off most optimizations for instance): if there is only one parameter left, then the problem is decidable. For a number of parameters between 2 and 8, the problem is still open [START_REF] Matiyasevich | Elimination of quantifiers from arithmetical formulas defining recursively enumerable sets[END_REF] and Matiyasevich conjectured it as undecidable. LEMMA 3. The Modified Best-Parameters Problem is undecidable if the performance prediction function t = a(P, s) is an arbitrary polynomial and if there are at least 9 integer optimization parameters.

PROOF. The proof is based on a result published in 1982: given an arbitrary polynomial f with nine variables, Jones [START_REF] Jones | Universal Diophantine Equation[END_REF] proved that there is no recursive function which can determine whether f has a non-negative integer zero, in the sense that it finds an explicit zero or returns null otherwise.

Finding parameter values -→ k such that, for an arbitrary polynomial t and for some given constant value I, t( -→ k , I) < T is equivalent to finding the zeros of an arbitrary polynomial. Given a polynomial t, the polynomial T * t( -→ x , I) * t( -→ x , I) reaches a value lower than T for some -→ x only if -→ x is a zero of t. This shows that to find the values for the bound is as difficult as finding the zeros of a polynomial. If a generates arbitrary polynomials, according to the value of P, and there are at least 9 optimization parameters (I is not considered as a variable, as its value is constant), then Modified Best-Parameters Problem is undecidable.

Examples of Simplified Decidable Cases

Our formal problem Best-Parameters is the formal writing of library optimizations. Indeed, in such area of program optimizations, the applications are given with a training data set. Then, people try to find the best parameter values of optimizing modules (inside a compiler usually with a given compilation sequence) that holds in the best performance. In this section, we show that some simplified instances of Best-Parameters problem becomes easily decidable. A first example is the OCEAN project [START_REF] Aarts | OCEANS: Optimizing Compilers for Embedded Applications[END_REF], and a second one is the ATLAS framework [START_REF] Whaley | Automated Empirical Optimizations of Software and the ATLAS Project[END_REF].

The OCEAN project [START_REF] Aarts | OCEANS: Optimizing Compilers for Embedded Applications[END_REF] optimizes a given program for a given data set by exploring all combinations of parameter values. Potentially, such value space is infinite. However, OCEAN restricts the exploration to finite set of parameter intervals. Consequently, the number of parameter combinations becomes finite, allowing a trivial exhaustive search of the best parameter values: each optimized program resulting from a particular value of the optimization parameters is generated and evaluated. The one performing best is chosen. Of course, if we use such exhaustive search, the optimizing compilation time become very high. So, one can provide efficient heuristics for exploring the bounded space of the parameters [START_REF] Triantafyllis | Compiler Optimization-Space Exploration[END_REF]. Currently, this is outside the scope of our article.

ATLAS [START_REF] Whaley | Automated Empirical Optimizations of Software and the ATLAS Project[END_REF] is another simplified case of the Best-Parameter problem. In the case of ATLAS, the optimization sequence is known, the programs to optimize are known (BLAS variants), and it is assumed that the performance does not depend on the value of the input (independence w.r.t. the matrix and vector values). Moreover, there is a performance model for the cache hierarchy (basically, the size of the cache) that, combined to the dynamic performance evaluation, limits the number of program executions (i.e., performance evaluation) to do. For one level of cache and for matrix-matrix multiplication, there are three levels of blocking controlled by three parameters, bounded by the cache size and a small number of loop interchanges possible (for locality). Exhaustive enumeration inside admissible values enable to find the best parameter value.

Figure 2 synthesizes a whole view of the different classes of the investigated problems with their decidability results. The largest class of the best parameters exploration problem that we consider, denoted by C1, assumes a finite set of optimization parameters with unbounded values (infinite space); The compiler optimization sequence is assumed fixed. If the performance prediction function is arbitrary, then this problem is undecidable. The second class of the best parameters exploration problem, denoted by C2 ⊂ C1, has the same hypothesis as C1 except that the performance model is assumed as an arbitrary polynomial function. The problem is undecidable too. However, a trivial identified decidable class is the case of bounded (finite) parameters space. This is the case of the tools ATLAS (class C3) and OCEAN (class C4).

FUTURE WORK

The phase ordering problem studied in this article does not make any assumption about the kind or the family of the considered program transformations. Potentially, we can have an unbounded (but finite) number of optimizing modules inside a compiler, as long as they guarantee us the best performance. Consequently, the size of the compiler can be as large as we require. In a future work, we want to explore the phase ordering problem with an additional restriction which is the granted size to a compiler. For this purpose, and thanks to the results presented in [START_REF] Cohen | A Polyhedral Approach to Ease the Composition of Program Transformations[END_REF], we will restrict the family of program transformations to the polyhedral ones. Indeed, the authors in [START_REF] Cohen | A Polyhedral Approach to Ease the Composition of Program Transformations[END_REF] give a matrix coding of all polyhedral transformations and their possible combinations: the size of such matrix is finite and bounded, while its elements define all possible affine polyhedral program transformations. This matrix coding has the benefit for consuming a bounded space and allowing to ease the composition of program transformation.

Another future work to this article is to study the phase ordering problem with another kind of restriction. Instead of limiting the size of the optimizing compiler, we can put a limit on the size of the final transformed program. Does phase-ordering become decidable in this case ?

Finally, an important open problem remains the definition of a general family of performance prediction functions that makes the phase ordering problem decidable. In this paper, we proved that if such function requires the execution or the simulation of the considered program, then the phase ordering problem becomes undecidable. But what if the performance predictor does not require neither the execution nor the simulation of the program ? Of course, if the performance modeling is too trivial or too simple (see Section 4), it is highly probable that the phase ordering problem becomes decidable, but in this case the model would not fit the real program performance. So, we require to define a more general family of performance prediction functions, that are efficient enough to accurately model the real program performance while allowing to have a decidable phase ordering problem. As a first step, we will consider for instance linear regression models.

CONCLUSION

As far as we know, our article is the first formalisation of two known problems: the phase ordering in optimizing compilation and the compiler optimization parameters space exploration. Our article sets down the formal definition of the phase ordering problem in many compilation schemes such as static compilation, iterative compilation and library generation. Given an input data set for the considered program, the defined phase ordering problem is to find a sequence of code transformations (taken from a finite set of code optimizations) that increase the performance up to a fixed objective. Alternatively, we can consider too parametric code optimization modules, and then we can define the formal problem of best parameters space exploration. However in this case, the compilation sequence is fixed, and the searching process looks for the best code optimization parameters that increase the program performance up to a fixed objective.

We showed that the decidability of both these problems is tightly correlated to the function used to predict or to evaluate the program performance. If such function is an arbitrary polynomial function, or if it requires to execute a Turing machine (by simulation or by real execution on the considered underlying hardware), then both these problems are undecidable. This means that we can never have automatic solutions for them. We provided some simplified cases that make these problems decidable: for instance, we showed that if we include a compilation cost in the model (compilation time, number of generated programs, number of compilation sequences, etc.), then the phase ordering becomes obviously decidable. This is what all actual ad-hoc iterative compilation techniques really do. Also, we showed that if the parameters space is explicitly considered as bounded, then the best compiler parameter space exploration problem becomes trivially decidable too.

Our article proves then that the requirement to execute or to simulate a program is a major fundamental drawback for iterative compilation and for library generation in general. Indeed, they try to solve a problem that can never have an automatic solution. Consequently, it is impossible to bring a formal method that allows to accurately compare between the actual ad-hoc or practical methods of iterative compilation or for library generation [START_REF] Cooper | Adaptive Optimizing Compilers for the 21st Century[END_REF][START_REF] Almagor | Finding effective compilation sequences[END_REF][START_REF] Zhao | A Model-Based Framework: An Approach for Profit-driven Optimization[END_REF][START_REF] Triantafyllis | Compiler Optimization-Space Exploration[END_REF]. The experiments that can be made to highlight the efficiency of a method can never bring a guarantee that such iterative method would be efficient for other benchmarks. As a corollary, we can safely state that, since it is impossible to mathematically compare between iterative compilation methods (or between library generation tools) then we can consider that any proposed method is sufficiently "good" for only its set of experimented benchmarks and cannot be generalized as a concept or as a method.

Our article proves too that using iterative or dynamic methods for compilation is not fundamentally helpful for solving the general problem of code optimization. Such dynamic and iterative methods define distinct optimization problems that are unfortunately as undecidable as static code optimizations, even with fixed input data.

However, our article does not yet give information about the decidability of phase ordering or parameters space exploration if the performance prediction function does not require program execution. Simply because the answer depends on the nature of such function. If such function is too simple, then it is highly probable that the phase ordering becomes decidable but the experimental results would be weak (since the performance prediction model would be inaccurate). The problem of performance modeling then becomes the essential question. As far as we know, we did not find any model in the literature that has been formally validated by statistical fitting checks as explained in [START_REF] Jain | The Art of Computer Systems Performance Analysis : Techniques for Experimental Design, Measurement, Simulation, and Modeling[END_REF].

Figure 1 :

 1 Figure 1: Classes of Phase-Ordering Problems

Figure 2 :

 2 Figure 2: Classes of Best-Parameters Problems

  intermediate form. When a native code generation for a program part is carried out, there is no way to re-optimize such program portion, and the process continues for other sub-trees until finishing the whole tree. Note that the optimization process for each sub-tree is applied by a finite set of program transformations. In other words, generative compilers look for local "optimized" code instead of a global optimized program.

	Algorithm 2 Optimize Node(n)
	Require: an abstract syntax tree with root n
	Require: a finite set of program transformations M
	if n is not leaf then
	for all u child of n do
	Optimize N ode(u)
	end for
	/*Generate all possible codes and choose the best one*/
	best ← φ {best code optimization}
	time ← ∞ {best performance}
	for all m ∈ M do
	if t(n, m) ≤ time then
	best ← m
	time ← t(n, m)
	end if
	end for
	apply the best transformation to the node n without changing
	any child
	else {Generate all possible codes and choose the best one}
	best ← φ {best code optimization}
	time ← ∞ {best performance}
	for all m ∈ M do
	if t(n, m) ≤ time then
	best ← m
	time ← t(n, m)
	end if
	end for
	apply the best transformation to the node n
	end if
	program parts in

Indeed, the cited paper does not contain a formal detailed proof, but a persuasive reasoning.

If the speculation is static, the code size is finite. If speculation is made dynamic, the hardware speculative window is bounded.

The time on an executing machine is discrete since we have clock cycles.

The latest version of SPIRAL use more elaborate strategies, but still does no resort to exhaustive search/test.

phases (to 10 or 15 as example). Consequently, the number of possible combinations becomes bounded which makes the problem of phase ordering decidable. Algorithm 1 can be used to generate the best compilation sequence if we consider a cost function as a fixed number of phases.

The next section presents another simplified case in phase ordering, which is generative compilation.


Alp Guler

Nikolaos Kardaris

Siddhartha Chandra

Vassilis Pitsikalis

Christian Werner

Klaus Hauer

Costas Tzsafestas

Petros Maragos

Iasonas Kokkinos

Costas Tzafestas

Human Joint Angle Estimation and Gesture Recognition for Assistive Robotic Vision

établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The increase in elderly population is a fact worldwide [1]. In this context computer vision and machine learning research applied on human-robot-interaction from the perspective of assistive living has both scientific interest and social benefits. In this work we focus on two prominent directions and apply the respective methods in the context of a challenging assistive living human-robot interaction scenario. This involves a robotic rollator that interacts with the elderly subjects using visual sensors, assisting them in every-day activities. These directions involve the use of state of the art deep learning based approaches for human joint angle estimation for the future goal of subject stability estimation, as well as the application of action recognition methods to enable elderly subjects interact with the robot by means of manual gestures. Herein we focus on the visual processing pipelines of this interface, and show a variety of rich applications and experiments.

There has been a furore of activity on the pose estimation front in recent years. Pose estimation usually involves inferring the locations of landmarks or body parts and the quality of the prediction is measured by metrics that involve comparing the predicted and the ground truth locations in the image plane. In this work we address the problem of estimating human joint angles. The joint angle estimation task involves estimating the angles made by segments of the human body at the joint landmarks in world coordinates. More specifically, we are interested in estimating (a) the knee angles, that is, the angles between by the thigh and the shank segments of the left and right legs, and (b) the hip angles, or the angles made by the torso and thigh segments. These angles will later be used to determine if the user's pose is unstable. In case of instability, the robot can assist the user achieve a stable position by exerting a physical force of the optimal magnitude in the optimal direction.

Human action recognition is a very active research area with a multitude of applications, such as video retrieval, health monitoring as well as human-computer interaction for assistive robotics. Proposed approaches on action recognition span several directions, such as deep architectures [START_REF] Simonyan | Two-stream convolutional networks for action recognition in videos[END_REF] and local spatio-temporal features with the popular Bag-of-Visual-Words (BoVW) paradigm [START_REF] Niebles | A hierarchical model of shape and appearance for human action classification[END_REF][START_REF] Niebles | Unsupervised learning of human action categories using spatial-temporal words[END_REF][START_REF] Marszalek | Actions in context[END_REF], as well as one of the top performing approaches, the dense trajectories [START_REF] Wang | Action recognition by dense trajectories[END_REF][START_REF] Wang | Action recognition with improved trajectories[END_REF]. Application on human-robot interaction [START_REF] Fanello | Keep it simple and sparse: Real-time action recognition[END_REF] has received relatively little attention mainly due to the increased computational requirements of most action recognition systems, which prevents them from performing in real-time. We target the rocognition of elderly subjects' gestural commands by employing dense trajectories and exploring alternative encoding methods. We have implemented a real-time version of our gesture recognition system that uses an activity detector and is currently integrated in the robotic platform.

In this work we present our contributions on applying these approaches on data containing actual elderly subjects. We are guided by a vision whereby assistive humanrobot interaction is advanced by state-of-the-art results in mid and higher level vision applications. The assistive scenario involves a robotic rollator equiped with multiple sensors, as shown in Fig. 2, that is capable of estimating the pose and the joint angles, as well as recognizing manual gestures using the proposed approaches. We utilize a deep learning based 2D joint localization approach fused with 3D information acquired from RGB-D sensors to localize joints in 3D and estimate the angles between line segments(see Sec.3). Moreover, we apply our action-gesture recognition approach based on the recent dense trajectories features [START_REF] Wang | Action recognition by dense trajectories[END_REF], employing a variety of encoding schemes. We also feed mid-level vision information [START_REF] Chandra | Accurate human-limb segmentation in rgb-d images for intelligent mobility assistance robots[END_REF] to the higher level of actiongesture recognition. All cases are evaluated on rich scenarios of a new dataset and task with elderly subjects showing promising results. Finally, we examine practical issues concerning an online version of the system that is integrated in the robotic platform and has close-to-real-time performance.

The human-robot interaction assistive dataset

The experimental prototype used for data acquisition [START_REF] Fotinea | Data acquisition towards defining a multimodal interaction model for human-assistive robot communication[END_REF] consists of a robotic rollator equipped with sensors such as laser range sensors that scan the walking area for environment mapping, obstacle detection and lower limbs movement detection, force/torque handle sensors and visual sensors: an HD camera to record patient's upper body movements and two Kinect sensors. The first Kinect captures the torso, waist and hips and the second faces downwards towards the lower limbs. The recording process involved acquiring RGB-D videos using the open-source Robotics Operating System (ROS) software capturing human subjects in a set of predefined use-cases and scenarios.

The dataset we use for estimating joint angles (Fig. 2) consists of (a) colour and depth images captured by the ROS-Kinect system, and (b) human-joint landmark point Fig. 1: Sample gestures from Task-6a and illustration of the visual processing with dense trajectories for gesture recognition (Sec. 4). First row: "Come here". Second row: "I want to stand up". trajectories captured by a Qualisys motion capture (MoCap) system. This dataset has 9K images coming from 27 video recordings.

The data acquired for gesture recognition ("Task-6a") comprises recordings of the patient sitting in front of the rollator, which is placed at a distance of 2.5 meters. Task-6a includes 19 different gestural and verbal commands. Each command is performed by the 13 patients, 3 -6 times. Sample gestures are depicted in Fig. 1. The task is challenging, as mobility disabilities seriously impede the performance ability of a verbal and/or gestural command. Moreover, due to the cognitive disabilities of some users, in some cases we observe different pronunciations of a command even among multiple performances of the same user. In essence, we miss the consistency that is sometimes assumed in other datasets [START_REF] Escalera | Chalearn multi-modal gesture recognition 2013: grand challenge and workshop summary[END_REF]. In addition, background noise and other people moving in the scene make the recognition task even harder.

Human Joint Angle Estimation from Kinect RGB-D Images

Our first contribution is a deep learning framework for human pose estimation. More precisely, we are interested in estimating the hip and the knee angles. We define (a) the hip angle as the angle made at the hip joint by the shoulders and the knees, and (b) the knee angle as the angle made at the knee by the hips and the ankles. These angles give a fair indication of the human pose, and can be exploited to determine if the user's pose is unstable.

Pose Estimation Dataset

Our objective is to use the RGB-D images to estimate the joint-angles in the world coordinate system, and evaluate them against the angles computed from the pointtrajectories. We assume that any deviation between the actual joint angles and those computed from the point-trajectories, due to an offset between the motion capture markers and the human body, is neglegible. Our dataset suffers from several severe limitations. Firstly, the upper, lower kinect and MoCap systems are not aligned in time. Secondly, the upper and lower kinect sensors look at parts of the human body individually, Fig. 2: Example images from our dataset. Notice the scale changes, and noise in the depth images. Our dataset also suffers from occlusion and truncation of body parts. We also have multiple people in the videos.

and do not see the entire context. Finally, the depth images from the kinect are sparse and noisy (fig. 2). We alleviate the first two limitations by using preprocessing procedures described in the rest of this section. To cope with noisy depth estimates, we use a sparse coding based strategy to denoise the predictions. This is described in section 3.6.

Data Alignment As mentioned, our input data and our ground truth are recorded by two different systems, and are not aligned. To evaluate the accuracy of our predictions, we require that the images are aligned-in-time with the joint-angle trajectories. More specifically, we need to know which image-frame corresponds to which time-stamp in the MoCap trajectory. This alignment problem is non-trivial due to several factors. Firstly, the two recording systems do not start at the same time. Secondly, the two systems capture data at different frequencies. Besides, while the MoCap system has a constant frequency, the frequency of the ROS-kinect system varies with the system load. Due to this, the upper-body frames and lower-body frames coming from the two kinect sensors are not aligned. To fix the variable frame rate issue, we re-sample the kinect captured videos at 10 fps via nearest neighbour interpolation. To align the upper and lower body frames, we only use images from the time-stamps when both upper and lower frames are available, discarding the remaining frames. Finally, to align the MoCap trajectories with the kinect frames, we design an optimization framework. More specifically, we define variation-functions for the images and trajectories. The variation functions can be understood as how the data changes over time. Once we define these functions, we can align the data by ensuring that these variations are correlated in time, since each variation results from the same cause, namely the human subject moving.

We denote by f x (t) the variation function of the point trajectories. f x (t) is defined to be simply the L 2 norm of the position of center of the human body at each timestamp t in the point trajectories, where b(t) denotes the center of the body (the pelvis). We define, f y (t), the variation function of the images I as the L 2 distance, in the HOG feature space, of an image at timestamp t from the first image (t = 0) in the sequence.

f x (t) = ||b(t)|| 2 2 , (1)
f y (t) = ||I hog t -I hog 0 || 2 2 (2)
The alignment problem can now be solved by aligning the f x and f y signals. To achieve this, we minimize the cross-correlation between them by exhaustively sampling a delay parameter between them. The exhaustive search takes less than 2-3 minutes per video. Figure 3 shows the variation functions as the result of our alignment procedure for two different videos.

Image Stitching. In our experiments, we observed that images containing only the lower or the upper body are very hard samples for pose estimation methods because these lack a lot of informative context about the pose. To cope with this issue, we stitch the images from the lower and upper kinect to obtain a single image containing the full human body. This stitching is done by estimating a planar homography transformation between manually selected keypoints. Using this geometric transformation we overlay the upper and lower sensor images to get the mosaic image. However, the mosaic images look unrealistic due to perspective distortion. We fix this by performing perspective correction on the mosaic image by estimating the necessary affine transformation. Since the kinect sensors are stationary, this procedure is done only once, we use the same transformation matrices for all images in our dataset. Results of image stitching can be seen in fig. 2. In the following sections, we describe our approach.

Related Work

In the recent years, deep learning has emerged as the gold standard of machine learning on nearly all benchmarks, including the task of human pose estimation. A number of deep learning methods solve the pose estimation problem by first estimating the spatial locations of interesting landmarks, and then inferring the pose from these landmark locations. In a broad sense, these approaches can be classified into two categories. The first category consists of methods that directly recover spatial locations of interesting landmarks via regression. Inspired by the deep convolutional cascaded network of Sun et al. [START_REF] Sun | Deep convolutional network cascade for facial point detection[END_REF], the Deep Pose approach by Toshev et al. [START_REF] Toshev | Deeppose: Human pose estimation via deep neural networks[END_REF] treats pose estimation as a regression task, refining the predictions iteratively via a cascaded network architecture. The second class of methods estimate probability distributions of the presence of these landmarks, which can then be treated as unaries in an inference procedure that enforces spatial constraints to recover the likelihoods of these landmarks ( [START_REF] Tompson | Joint Training of a Convolutional Network and a Graphical Model for Human Pose Estimation[END_REF][START_REF] Wei | Convolutional pose machines[END_REF][START_REF] Belagiannis | Recurrent human pose estimation[END_REF][START_REF] Lifshitz | Human pose estimation using deep consensus voting[END_REF][START_REF] Haque | Viewpoint invariant 3d human pose estimation with recurrent error feedback[END_REF]). These approaches typically use fully convolutional neural networks to estimate a per pixel likelihood of joint locations. A number of these approaches [START_REF] Ramakrishna | Pose machines: Articulated pose estimation via inference machines[END_REF][START_REF] Carreira | Human pose estimation with iterative error feedback[END_REF] employ iterative error feedback, that is, a cascaded system to repeatedly refine the per pixel likelihoods, similar to recurrent networks. The DeeperCut approach [START_REF] Insafutdinov | Deepercut: A deeper, stronger, and faster multi-person pose estimation model[END_REF] exploits a very deep residual neural network and achieves promising results without using a sequential refinement scheme. More recently, Haque et al. [START_REF] Haque | Viewpoint invariant 3d human pose estimation with recurrent error feedback[END_REF] proposed a viewpoint invariant pose estimation framework for RGB-D images that utilizes Long Short Term Memory (LSTM) [START_REF] Hochreiter | Long short-term memory[END_REF] to incorporate the idea of error feedback.

Joint Angle Estimation from RGB-D Images

The joint angle estimation task can be addressed by either (a) directly regressing the angles from the images, or (b) first estimating the spatial locations of the human parts, and then estimating the joint angles using geometry. As described in section 3.2, we have two broad classes of methods that address the problem of spatial localization of landmarks. An advantage of methods that directly regress landmark locations from the image data is that they can learn pairwise interactions in a fully connected graph efficiently via fully connected layers, thereby exploiting the holistic context in the image.

However, these methods typically use an object detector at the testing stage [START_REF] Toshev | Deeppose: Human pose estimation via deep neural networks[END_REF] and then regress from the detected bounding boxes to the landmark locations. Object detection is necessary to prevent noise in the context around the object of interest from corrupting the features. While this is a reasonable strategy, this two step procedure can be undesirable for datasets containing low resolution images, or images captured in the wild where object detection can be inaccurate. Methods that estimate the probability distributions can work around this requirement by estimating the probability distributions over the entire image, and then inferring the landmark locations by enforcing spatial constraints. Our dataset has noisy, low resolution images, with significant occlusion, and dramatic scale changes. This motivates us to employ a fully convolutional network to estimate the spatial locations of the human parts in the image plane.

For our task of 3D localization of the human parts, it is natural to use the depth cues from the depth images. The depth field of view of the Microsoft Kinect sensor ranges from 25 centimeters to a few meters. This leads to a noisy and sparse reconstruction of the depth field in cases where the object of interest is close to the sensor. Our dataset has numerous instances where the human is too close to the kinect sensor. Consequently, the depth information is unreliable (or absent) in these cases. The RGB sensors do not suffer from this limitation. This motivates the usage of the clean RGB image for the estimation of joint locations, later to be combined with the depth information to reconstruct the 3d joint locations. We therefore, first estimate the part positions from the RGB images, then use the depth images to reconstruct the part positions in the world system, and finally estimate the joint angles via geometry. Our pipeline is described in Fig. 4. For estimating the part positions in 2D, we use the Convolutional Pose Machines [START_REF] Wei | Convolutional pose machines[END_REF] approach which achieves the state of the art results on the challenging MPII dataset [START_REF] Andriluka | 2d human pose estimation: New benchmark and state of the art analysis[END_REF]. This approach was trained using thirty thousand images from the MPII and Leeds Sports [START_REF] Andriluka | 2d human pose estimation: New benchmark and state of the art analysis[END_REF][START_REF] Johnson | Clustered pose and nonlinear appearance models for human pose estimation[END_REF] datasets.

Convolutional Pose Machines

As mentioned in section 3.3, our method is based on the convolutional pose machines [START_REF] Wei | Convolutional pose machines[END_REF] approach, that uses sequential prediction convolutional blocks that operate on image and intermediate belief maps and learn implicit image-dependent spatial models of the relationships between parts. We now briefly describe the prediction strategy.

The human part localization task is done in several stages. The first stage of the method predicts part beliefs from local image evidence via a fully convolutional network. The network has seven convolutional layers, trained to minimize the L 2 distance between the predicted beliefs and the ground truth beliefs. The ground truth beliefs are synthesized by putting Gaussian peaks at ground truth locations of parts. The second stage learns convolutional features from the image, and combines them with the beliefs from the first stage (via concatenation), and again learns a classifier to predict part beliefs from this combination of features and beliefs using the L 2 loss. The third stage predicts the part beliefs by combining features and the beliefs from the second stage and so on. While the network in [START_REF] Wei | Convolutional pose machines[END_REF] uses six stages, our network uses a cascade of four stages. The network is designed to learn interactions within increasingly larger neighbourhoods by repeatedly increasing the receptive fields of the convolutions in each subsequent stage of the cascade. Larger receptive fields are able to model long-range spatial interactions between the landmarks, and help impose geometrical constraints on the locations of the body parts, such as the knee, the hip, the shoulder, etc. The cascaded network predicts the per-pixel beliefs for fourteen landmarks, namely the neck, head, and the left and right ankle, knee, hip, wrist, elbow, shoulder and wrist.

Joint Angles from Part Locations

The method described in section 3.4 predicts human part locations in the image plane. However, our objective is to estimate joint angles in the world space. After the part beliefs are estimated in the mosaic image, we unstitch the mosaic image, and transfer the part beliefs to the upper and lower images. We then construct the 3D point cloud using the depth images and the intrinsic parameters of the kinect sensors, and transfer the part beliefs to the point clouds. This is followed by registering the point clouds of the upper and lower kinect sensors so that they are in the same space, and transferring the part beliefs to the common system. For each part, we choose the location to be the point with the maximum belief. If the maximum belief for a part is below a threshold , we conclude that this part is either occluded or truncated. This gives us 3D spatial locations of the body parts. Given the part locations in the world coordinate system, we compute the joint angles using vector algebra.

Denoising via Sparse Coding

As described in section 3.3, our dataset suffers from drastic scale changes causing the depth to be sparse and noisy. This causes the angle estimates to be noisy. To overcome this difficulty, we pose reconstruction of our predictions as a classic sparse coding based denoising approach. We estimate a dictionary D of codewords from the ground truth angles x by solving the following optimization

D * = argmin D,w ||x -Dw|| 2 2 + λ 1 ||w|| 1 + λ 2 ||w|| 2 2 . (3)
The reconstruction of noisy predictions x is done by estimating a sparse linear combination w * of the dictionary as

w * = argmin w ||x -D * w|| 2 2 + λ 1 ||w|| 1 + λ 2 ||w|| 2 2 y = D * w * (4)
This reconstruction suggests that the predictions can be represented sparsely over the dictionary, and helps tilt the distribution of the predicted angles towards the ground truth angles. We use this technique to denoise our hip angle predictions. We report results of denoising in section 3.7.

Experiments and Discussions

In this section we report empirical results for the estimation of hip, and knee angles. As described in section 3.1, our dataset has 27 videos, containing about 9K image frames.

We use the point trajectories of the shoulders, hips, knees and angles from the MoCap data to compute the ground truth angles. Our prediction pipeline is described in section 3.4. Our networks are trained on 30K images from the MPII and Leeds Action Dataset, as mentioned in section 3.4, and are the state of the art on these challenging datasets. The evaluation criterion is the Root Mean Square Error (RMSE) in degrees between the predicted angles and the ground truth angles. We also report the detection rate for each angle. An angle is said to be detected if the locations of all three parts constructing the angle are available. The detection rate for an angle is defined to be the ratio of frames in which the angle is detected.

We first study the effect of changing the threshold described in section 3.5. This threshold is applied to the maximum beliefs of parts, and introduces a measure of confidence in the part localization. If the maximum belief of a part is below we conclude that the part is absent. This absence occurs because of occlusion or truncation of parts of interest. Fig. 5 shows the variation of RMSE and detection rate for different values of . As the minimum accepted confidence increases, RMSE decreases substantially and naturally the detection rate drops due to increasing number of rejected angles. High confidence in part locations, which mainly depends on visibility of parts, leads to better estimates of joint angles. This increase in performance with large is observed despite the remaining erroneous joint locations, which would have been corrected if more context (and in many cases joints themselves) were visible. This provides evidence that a better placement of the cameras such that full body context is provided to the algorithm would lead to an increase in the performance of the algorithm. To further emphasize this, we exemplify estimated angles and groundtruth angles as a function of time for the same subject and same task under two different camera configurations in Fig. 6. In the configuration where the patient is closer to the sensors and the ankles and shoulders are occluded the algorithm is not able to estimate the angles, whereas in the setting where the subject is distanced from the sensor, the angle estimation is more accurate. For our quantitative results, we choose = 0.2, which gives a detection rate of 0.75 (figure 5). Our results are reported in table 1. We report the RMSE values for the estimation of the left and right knee and hip angles, averaged over the entire dataset. It can be seen that the errors in estimation of the hip angles are higher than those corresponding to the knee angles. This is due to two factors: (a) the shoulders are truncated when the human subject to too close to the camera (fig. 2 cluded by the elbows and hands (fig. 2, row 1). Table 1 also reports the estimation errors after the denoising process described in section 3.6. To estimate the dictionary (eq. 3), we use a leave-one-out strategy. The dictionary for denoising of a video is constructed using the dataset, leaving out this particular video. The parameters for the dictionary learning procedure are as follows: λ 1 = 10, λ 2 = 0.001, and the size of the codebook K = 200. We see an improvement in estimation of both the knee and hip angles after denoising, however the improvement in case of the hip angles is more prominent.

Dense Trajectories for Gesture Recognition

Gesture recognition allows the interaction of the elderly subjects with the robotic platform through a predefined set of gestural commands. Our gesture classification pipeline, depicted in Fig. 7, employs Dense Trajectories features along with the Bag-of-Visual-Words framework. We briefly describe the individual steps involved and present our experimental results on Task-6a. Dense trajectories [START_REF] Wang | Action recognition by dense trajectories[END_REF] consists in sampling feature points from each video frame on a regular grid and tracking them through time based on optical flow. Different descriptors are computed within space-time volumes along each trajectory. Descriptors are: the Trajectory descriptor, HOG [START_REF] Wang | Evaluation of local spatio-temporal features for action recognition[END_REF], HOF [START_REF] Laptev | Learning realistic human actions from movies[END_REF] and MBH [START_REF] Wang | Evaluation of local spatio-temporal features for action recognition[END_REF]. Features are encoded using BoVW, VLAD or Fisher vector to form a video representation. BoVW uses a codebook which is constructed by clustering a subset of ran-Fig. 7: Gesture classification pipeline.

domly selected training features into K = 4000 clusters. Each trajectory is assigned to its closest visual word and a histogram of visual word occurrences is computed, yielding a sparse K-dimensional representation. VLAD encodes first order statistics among features by computing differences between extracted features and the visual words. Fisher vector encodes first and second order statistics using a Gaussian Mixture Model (GMM) with K = 256 gaussians. To compensate for the high dimensionality of FV, we reduce the descriptors' dimensionality by a factor of 2 using Principal Component Analysis.

Videos encoded with VLAD or Fisher vector are classified using linear support vector machines (SVMs). Different descriptors are combined by concatenating their respective VLAD or FV vectors. In the BoVW case, videos are classified using SVMs with the χ 2 [START_REF] Wang | Evaluation of local spatio-temporal features for action recognition[END_REF] kernel. Different descriptors are combined in a multichannel approach, by computing distances between BoVW histograms as:

K (h i , h j ) = exp - c 1 A c D h c i , h c j , ( 5
)
where c is the c-th channel, i.e. h c i is the BoVW representation of the i-th video, computed for the c-th descriptor, and A c is the mean value of χ 2 distances D h c i , h c j between all pairs of training samples. Since we face multiclass classification problems, we follow the one-against-all approach and select the class with the highest score.

Experiments Gesture classification is carried out on a subset of the Task-6a dataset comprising 8 subjects and 8 gestures1 , without limiting the generality of results. Results are shown in Table 2. It is evident that the large variability of the gesture performance among patients has a great impact on performance. The combined descriptor performs consistently better, since it encodes complementary information extracted from the RGB channel. VLAD and Fisher vector further improve performance, since they encode rich information about the visual words' distribution. Fig. 8 depicts the mean confusion matrix over all patients computed for BoVW and the MBH descriptor. Naturally, gestures that are more easily confused by the elderly subjects are harder to classify correctly, e.g. "Help" and "PerformTask" both consist of a similar horizontal movement but in different height. To demonstrate the difficulty of the task, we use the same ges- 2).

Fig. 9: Sample frames from the gesture dataset acquired in [START_REF] Rodomagoulakis | A multimedia gesture dataset for human-robot communication: Acquisition, tools and recognition results[END_REF]. Table 3: Comparative results on Task-6a and the dataset acquired in [START_REF] Rodomagoulakis | A multimedia gesture dataset for human-robot communication: Acquisition, tools and recognition results[END_REF]. Mean classification accuracy (%) over all subjects is reported.

ture classification pipeline on the dataset acquired in [START_REF] Rodomagoulakis | A multimedia gesture dataset for human-robot communication: Acquisition, tools and recognition results[END_REF]. It includes 13 young subjects (Fig. 9) that perform the same 19 gestures as the ones in Task-6a under similar conditions. Training and testing is carried out with the same settings using BoW encoding.

Mean classification over all subjects is reported. Comparative results shown in Table 3 illustrate that variability in the execution of the gestures among healthy and cognitively intact subjects is effectively handled by our pipeline. This highlights the challenge that the development of a gestural communication interface for elderly people presents.

Filtering background trajectories

To further improve our gesture recognition results, we have worked towards the integration of the body part segmentation method introduced in [9] into our gesture recognition pipeline. Specifically, we use the subject's bounding box to reject noisy background trajectories (see Figure 10). We compute a separate BoVW histogram that corresponds to the area enclosed by the bounding box, which is used to augment the original BoVW vector. The mask shown in the Figure 10 is the output of the semantic segmentation algorithm introduced in [START_REF] Chandra | Accurate human-limb segmentation in rgb-d images for intelligent mobility assistance robots[END_REF] and applied on Task-6a data. Following the same experimental protocol as in the rest of the experiments described in this section, we obtained additional improvements on the related Task-6.a scenario, as shown in Table 4 below. Given the simplicity of the employed approach, results show remarkable improvement.

A holistic approach for gesture recognition, such as dense trajectories, can greatly benefit from exploiting mid-level information to remove background noise. Our current research plans include more efficient exploitation of the rich information contained in the output of the semantic segmentation pipeline. Table 4: Average classification accuracy (%) over all 8 patients using our baseline method (first column) and employing the foreground-background mask (second column). Results show a consistent improvement (third column) over multiple feature descriptors. Results are obtained using the BoW encoding.

On-line processing

Towards the realisation of a user interface that enables elderly subjects interact with the robotic platform, we have implemented an on-line version of our system that performs continuous gesture recognition. ROS is employed to provide the main software layer for interprocess communication.

Our overall system comprises two separate sub-systems: (a) the activity detector (AD node) that performs temporal localization of segments that contain visual activity and (b) the gesture classifier (GC node) that assigns each detected activity segment into a class. The AD node processes the RGB stream in a frame-by-frame basis and determines whether there is any visual activity in the scene, based on an activity "score", whose value is thresholded. The GC node is appropriately signaled at the beginning and the end of the activity segments. Small segments of activity are rejected to ensure that small spontaneous movements of the user are not processed.

The Gesture Classifier node processes video segments and assigns them to one of the pre-defined categories using the classification pipeline described previously. To reduce the total processing time we downsample RGB frames both in space and time by a factor of 2. The GC node caches input frames from the camera. When the appropriate signal is received, feature extraction begins immediately, starting from the indicated frame's timestamp. When the activity segment ends features are extracted from the remaining frames of the activity segment, the classification pipeline continues. The robot reacts to the recognized gestural command by either providing audio responses or moving in order to assist the elderly user. The system has been integrated in the robotic platform and used by actual elderly patients. It operates approximately at 20 fps on an Intel i7 CPU.

Conclusions

In this work we have examined multiple directions on an assistive human-robot interaction task, from the visual processing point of view. We have provided rich algorithmic and experimental evidence on how one may estimate joint angles of the elderly subjects with the given challenging setup and have signified the importance of the sensor localization that includes a full view of the subject. Additionally, we have shown that the adopted sparse coding based denoising approach increases performance. Another application concerns the automatic recognition of gesture commands, as well as a practical integrated system that online processes streaming video data, and accurately recognizes the commands with real-time performance. The systems are currently under validation studies conducted with elderly subjects in geriatric clinics, opening new perspectives in the field of robotic vision for assistive living. Ongoing and future plans involve the more deep integration of the multiple directions presented in this paper, by incorporating the information of the estimated pose and joint angles for gesture recognition, as well as the integration of gesture recognition in an overall network that jointly estimates the pose, and the higher level gesture concepts. Finally, further exploiting the computer vision interface applications with geriatric clinician experts can lead to a deeper understanding on how we should advance and jointly experiment with our approaches so that it would be for the best interest of the interdisciplinary research communities and above all for the benefit of the elderly population.

Fig. 3 :

 3 Fig.3: Results of our alignment procedure for 4 videos. The green curve shows the image variation function f x (t), the red curve shows the trajectory variation function f x (t). The human subjects in these videos performed the sit-stand exercise 3 times.

Fig. 4 :

 4 Fig.4: A schematic representation of our joint angle estimation pipeline. We have two stationary RGBD sensors (shown in the circle in the first column) on a moving rollator. We stitch the upper and lower RGB images, and perform part localization on the mosaic image using Convolutional Pose Machines. The mosaic image is then unstitched and the part beliefs are transferred to the registered point cloud. This gives us the spatial locations of the body parts in the world space. Given these part locations, we estimate the joint angles using vector algebra.

Fig. 5 :

 5 Fig. 5: Variation of Root Mean Square Error in degrees and detection rate with confidence threshold . As the confidence threshold increases (as the colours become hotter), the detection rate and the angle estimation errors decrease.

Fig. 6 :

 6 Fig. 6: Angle estimation on same subject and same task(sitting -standing) under different camera configurations. The left image corresponds to a time instance in the top row time series and the right image corresponds to an instance in the bottom row time series: Estimated (red) and the groundtruth angle (blue). The columns portray Left Hips, Right Hips, Left Knee and Right Knee angles respectively. Corresponding RMSE values are displayed on top of each time series.

1 Table 2 :

 12 5 56.0 44.1 51.5 58.3 56.3 45.8 81.3 57.0 69.9 66.3 combined 75.0 52.0 52.9 57.6 58.3 65.6 62.5 75.0 62.4 67.6 68.Classification accuracy (%) per patient on a subset of the Task-6a dataset that contains 8 gestures performed by 8 subjects. Results for different encoding methods are shown; "avg." stands for average accuracy over all patients.

Fig. 8 :

 8 Fig. 8: Mean confusion matrix over all subjects of the Task-6a dataset. The results are obtained with the MBH descriptor and BoW encoding (6th row of Table2).

Fig. 10 :

 10 Fig. 10: Combining semantic segmentation with feature extraction for gesture classification on the Task-6a. A separate BoVW histogram is computed for the whole frame and the area enclosed by the bounding box. The two vectors are combined into a single video representation.

Mobot System RGB Images Depth Images Stitched Image CNN Based Joint Location Estimation (CPM) Un-Stitching 3D Localization of Joints Point Cloud Generation Point Cloud Registration

Table 1 :

 1 Quantitative results of our Convolutional Pose Machines pipeline for the estimation of joint angles. The evaluation metric is the Root Mean Square Error in degrees. Denoising using a sparse coding strategy improves average performance by 5.4 degrees.

	Method	Left Knee Right Knee Left Hip Right Hip Average
	CPM	16.16	15.24	24.51	25.27	20.29
	CPM + Denoising 13.66	14.23	15.06	16.48	14.86

The 8 selected gestures are: "Help", "WantStandUp", "PerformTask", "WantSitDown", "ComeCloser", "ComeHere", "LetsGo", "Park".