Hugo Gilbert
email: hugo.gilbert@lip6.fr

Bruno Zanuttini
email: bruno.zanuttini@unicaen.fr

Paolo Viappiani
email: paolo.viappiani@lip6.fr

Paul Weng
email: paweng@cmu.edu

Esther Nicart
email: esther.nicart@unicaen.fr

Model-Free Reinforcement Learning with Skew-Symmetric Bilinear Utilities

In reinforcement learning, policies are typically evaluated according to the expectation of cumulated rewards. Researchers in decision theory have argued that more sophisticated decision criteria can better model the preferences of a decision maker. In particular, Skew-Symmetric Bilinear (SSB) utility functions generalize von Neumann and Morgenstern's expected utility (EU) theory to encompass rational decision behaviors that EU cannot accommodate. In this paper, we adopt an SSB utility function to compare policies in the reinforcement learning setting. We provide a model-free SSB reinforcement learning algorithm, SSB Q-learning, and prove its convergence towards a policy that is -optimal according to SSB. The proposed algorithm is an adaptation of fictitious play [Brown, 1951] combined with techniques from stochastic approximation [Borkar, 1997]. We also present some experimental results which evaluate our approach in a variety of settings.

INTRODUCTION

In problems of sequential decision-making under uncertainty (often represented as Markov Decision Problems-MDPs [Puterman, 1994]), an agent has to repeatedly choose according to her current state an action whose consequences are uncertain, in order to maximize a certain criterion in the long run. In most cases, the criterion chosen is the expectation of cumulated rewards, but it is not risk sensitive and fails to explain widely observed violations of axioms such as transitivity or von Neumann's independence axiom. One of the aims of decision theory is to provide criteria able to account for such behaviors.

Interestingly, the Skew-Symmetric Bilinear (SSB) utility theory [Fishburn, 1984] defines a family of decision criteria able to represent risk-averse and risk-seeking behaviors, intransitive choices and violations of the independence axiom. In particular it encompasses the expected utility (EU) model [START_REF] Von Neumann And Morgenstern | Theory of games and economic behavior[END_REF], which is the most popular risk sensitive criterion in decision theory. In SSB theory a binary functional ϕ over probability distributions is given, where the sign of ϕ(p, q) gives the preference between two distributions p and q. Furthermore, particular choices of ϕ allow decision criteria that only rely on ordinal pieces of information such as "this trajectory is preferred to this other trajectory" to be used. This property is of interest for MDPs as the optimal policy can be highly sensitive to the reward function, but specifying such a reward function is often difficult even for an expert user. Such preference-based approaches have received much attention lately [Akrour et al., 2012, Furnkranz et al., 2012, Busa-fekete et al., 2014[START_REF][END_REF], Wirth and Fürnkranz, 2013, Wirth and Neumann, 2015], and a special case of SSB utility function which optimizes the probability of yielding a preferred outcome has been investigated in various domains [START_REF] Busa-Fekete | [END_REF], Dudík et al., 2015, Rivest and Shen, 2010]. In reinforcement learning (RL), [START_REF] Busa-Fekete | [END_REF] provided a metaheuristic algorithm using evolutionary strategies to compute a "good" policy. Designing a learning scheme for an RL problem using an SSB utility function would therefore enable RL problems to be solved for a large class of criteria including "preference-based" criteria.

The possibility of intransitive preferences in SSB utility theory could be seen as a significant barrier to its use in automated decision making. However, the seminal work of Gilbert et al. [2015a] shows that an SSB-optimal strategy always exists as a mixture of policies, and furthermore, that in a finite horizon MDP where the model is known, such an SSB-optimal strategy can be computed using a double oracle approach. Unfortunately, this approach suffers two drawbacks: its time and space requirements might become prohibitive if the optimal mixture of policies is composed of too many policies, and it does not generalize to the case where the model is unknown, and thus cannot be used to derive a model-free RL algorithm.

We remedy this by providing a new algorithm which behaves optimally in an MDP with an SSB utility function, from which we derive a model-free RL algorithm in order to learn an SSB-optimal policy when the model is unknown. While the approach is still based on game-theoretic arguments, the algorithm differs in spirit and resembles an adaptation of a fictitious play algorithm [Brown, 1951] combined with Q-learning [START_REF] Watkins | [END_REF] using stochastic approximation techniques with two timescales [Borkar, 2008]. [START_REF] Kalathil | [END_REF] recently exploited a similar two-timescale technique in MDPs, but in a different context; using the average expected reward criterion with vector rewards, their goal is to learn a policy whose vectorial value approaches a fixed set.

The paper is organized as follows. In Section 2, we give the background elements and introduce our notations. In Section 3, we give a game theoretic view of the resolution of an RL problem with an SSB utility function and present a fictitious play approach to solve it. Lastly, in Section 4, we present some experimental results.

SSB MARKOV DECISION PROCESSES

We study in this paper episodic MDPs with finite state and action spaces. An MDP M is formally defined by a tuple (S, F, A, P, s 0 ) where: S is a finite collection of states; F = {f 1 , f 2 , . . . , f |F | } ⊂ S is a finite set of final states; A = {A s |s ∈ S} is a collection of finite sets of possible actions, one for each state; P is the transition function where P(s |s, a) is the probability that the state at time step t + 1 is s , given that the state at time step t was s and that the agent performed action a; s 0 ∈ S \ F is the initial state in which all episodes start. 1 Whereas preferences over state-action pairs are typically modeled with numerical rewards, in our framework we assume that the final states summarize the preference information. More precisely, we assume that the decision maker has a preference relation over possible final states, where f f means that ending in final state f is at least as good as ending in f . Note that we do not assume to be total or even transitive, thus accommodating a wide variety of preference behaviors, including those deviating from normative decision theory. A "standard" MDP (with rewards obtained at each time step) could still be represented in our setting with the notion of an augmented MDP [Gilbert et al., 2015a] at the cost of introducing additional states.

We assume (as standard in RL) that S, A, s 0 and F are known, that P is unknown and that at each step the agent knows exactly which state she occupies.

We call episode a succession of state-action pairs 1 A probability distribution P0 over initial states can easily be accommodated by using a dummy state s0 with one dummy action whose transitions to all the other states are governed by P0.

α i = 1, α i ≥ 0}
, where π is the mixed policy2 that randomly selects policy π i with probability α i at the beginning of each episode.

Example 1. As a running example, we consider a variant of the classical "Gardner dice" two-player game. Each player has three six-sided dice, written A, B, C and biased as shown in Table 1. Players simultaneously choose a die to throw, and whoever rolls the highest number wins.

It is easy to see that die A rolls higher than B most of the time, so die A should be preferred to B, but B mostly beats C, and C mostly beats A, hence the relation "more likely to win" is cyclic. The optimal policy for this problem is to play dice A, B and C with probabilities 3/13, 3/13 and 7/13, respectively [Gilbert et al., 2015a].

For illustrative purposes, we consider a variant where each player makes sequential decisions: she must first choose whether to throw die A (action a A ) or not (a BC ). If she chooses a BC , then she can choose to throw B (action a B ) or C (a C ). Clearly, this does not change the probabilities with which to throw each die in an optimal policy. This problem can be represented as an MDP where T max is 2, with an initial state s 0 where A s0 = {a A , a BC }; a state s BC (reached by choosing action a BC in s 0 ); and six final states {f 1 , . . . , f 6 } representing the numbers rolled. An example transition probability is P(f 3 |s BC , a B ) = 5/6, and an example episode is (s 0 , a BC , s BC , a C , f 5 ). An example (Markovian stationary deterministic) policy is π

B (s 0 ) = a BC , π B (s BC ) = a B .
Using similar notation for A and C, the optimal policy is the mixed policy π * =(π A |3/13,π B |3/13,π C |7/13), which dictates that the player draws one of π A ,π B ,π C at the start of an episode, following it for the whole episode. 3

Our aim is to find an optimal (defined in the next subsection) policy. Recall that the decision maker has a preference relation over possible final states; we want to compare policies by considering this preference relation. In other words, we want to lift the preference relation defined on final states to a preference relation defined on policies.

COMPARING POLICIES WITH AN SSB UTILITY FUNCTION

We assume throughout the paper that the agent's preferences between probability distributions are described by the SSB model as presented and axiomatized by Fishburn [1984]. In this model, an agent is endowed with a binary functional ϕ over ordered pairs (f, f ) ∈ F 2 of final states, indicating the intensity with which she prefers f to f , with f f ⇔ ϕ(f, f ) ≥ 0. Functional ϕ is assumed to be skew-symmetric, i.e., ϕ(f, f ) = -ϕ(f , f ) and is extended to the space of probability distributions over F by bilinearity (wrt the mixture operation on distributions).

The SSB criterion for comparing p and q is then written:

ϕ(p, q) = f,f ∈F p(f )q(f )ϕ(f, f ) (1)
where p(f ) (resp. q(f )) denotes the probability of reaching the final state f (resp. f ) in distribution p (resp. q). We have p q if ϕ(p, q) > 0 (strict preference), and p ∼ q if ϕ(p, q) = 0 (indifference).

Any policy π in an MDP induces a probability distribution p π over final states (reached after at most T max time steps); p π is referred to as the final state distribution of π. As comparing policies amounts to comparing their induced distributions, we write ϕ(π, π ) for ϕ(p π , p π ) to simplify notation and define the preference relation over policies:

π π ≡ ϕ(π, π ) ≥ 0 (2)
Example 2 (continued). The goal of beating the opponent's roll can be expressed as ϕ(f m , f n ) = 1 for m > n, -1 for m < n, and 0 for m = n. The deterministic policies π A , π B , π C amount to rolling the corresponding die, inducing the final state distributions p A , p B , p C of Table 1 3 Note the difference with the randomized (stationary) policy π which draws aA with probability 3/13 in s0, and, independently of this, draws aB (resp. aC ) with probability 3/10 (resp. 7/10) in sBC . The expectation of reaching each final state is nonetheless the same in π and π * (hence π is also optimal). over final states. We have (zero entries are irrelevant):

ϕ(p A , p B ) = p A (f 1 )p B (f 3 )ϕ(f 1 , f 3 ) + p A (f 1 )p B (f 6 )ϕ(f 1 , f 6 ) + p A (f 4 )p B (f 3 )ϕ(f 4 , f 3 ) + p A (f 4 )p B (f 6 )ϕ(f 4 , f 6 ) = 1 6 • 5 6 • (-1) + 1 6 • 1 6 • (-1) + 5 6 • 5 6 • 1 + 5 6 • 1 6 • (-1) = 14/36 > 0,
showing that die A should be preferred to die B.

If we were to define ϕ by ϕ(f m , f n ) = m -n, the strength of a victory (or defeat) would be taken into account. The policies in this case would be compared wrt the expectation of the roll (as explained below).

The SSB model is very general, as it can take into account choice intransitivity, which is widely observed in practice [Fishburn, 1991]. Moreover, the SSB model can represent different risk attitudes via an adequate choice of ϕ. For example, it represents risk-averse behavior (in the weak sense) if the certainty equivalent of a distribution p is less than or equal to its expected value, where the certainty equivalent of a distribution p is the element f such that ϕ(p, f ) = 0 (i.e., p ∼ f ). Nakamura [1989] shows how to design risk averse and risk seeking SSB utility functions.

The SSB model also encompasses many decision criteria, such as the expectation criterion ϕ(f, f ) = c(f ) -c(f ) (where c denotes a utility/cost function); the probability threshold criterion [START_REF] Yu | [END_REF]

ϕ(f, f ) = 1 c(f )≥τ - 1 c(f )≥τ , which states that p q if c(f )≥τ p(f ) > c(f )≥τ q(f ) for a threshold τ ∈ R; and the dominance relation ϕ(f, f ) = 1 (resp. 0, -1) if f f (resp. f ∼ f , f ≺ f ), which states that p q if f f p(f )q(f ) > f f p(f )q(f ) (as in Example 2).
In other words p is preferred to q if a final state generated according to p is more likely to be preferred to a final state generated according to q than the converse. This is called probabilistic dominance in the following.

Probabilistic Dominance (PD) is interesting as it only relies on ordinal pieces of information. Its axiomatic characterization was given by Blavatskyy [2006] and it has been explored lately in various domains such as RL [START_REF] Busa-Fekete | [END_REF], voting systems [START_REF] Rivest | [END_REF] and dueling bandits [START_REF] Dudík | [END_REF]. In the latter work, the authors adopt the name of von Neumann solution for the PD optimal solution. Indeed, as will be discussed in the next section, finding an SSB-optimal policy (and in particular, finding an optimal policy according to PD) is equivalent to finding a Nash equilibrium in a finite zero-sum two-player game. Thus the existence of a von Neumann solution is implied by von Neumann's minimax theorem.

In standard MDPs, the optimal policy can be highly sensitive to the reward function used, and yet designing a numerical reward function is often cognitively difficult, even for an expert user. This issue is tackled in preferencebased approaches [Akrour et al., 2012, Weng and Zanuttini., 2013, Gilbert et al., 2015b, Weng et al., 2013, Busafekete et al., 2014, Furnkranz et al., 2012[START_REF][END_REF], Wirth and Fürnkranz, 2013, Wirth and Neumann, 2015] by only considering ordinal pieces of information, such as feedbacks of the type "this trajectory is preferable to that one". We can distinguish two types of approach. The first [START_REF] Wirth | [END_REF], Wirth and Neumann, 2015, Weng et al., 2013] aims to recover a numerical reward function that explains most of the expressed preferences of the user and can be used with classic criteria. The second [START_REF] Busa-Fekete | [END_REF], Furnkranz et al., 2012] deals with purely ordinal criteria, e.g., [START_REF] Busa-Fekete | [END_REF] find a good policy (wrt to PD) using a metaheuristic algorithm based on evolutionary strategies in finite horizon continuous MDPs. Example 3. Consider a car racing against the clock. The driver aims to complete the race in the shortest time possible, but must find the right compromise between speed and the risk of being eliminated; driving too fast can cause the car to run off the track and be eliminated from the competition; the relation between the car's speed and going off track is stochastic (the faster, the more likely). This is easily modeled as an MDP where policies induce a distribution over the race's possible outcomes. The preference over final states is determined by two conditions: (i) race completion is always strictly preferred to an elimination, and (ii) a trajectory completing the race in time t 1 is strictly preferred to completion in time t 2 > t 1 .

In this model, adopting SSB with probabilistic dominance finds the "best" policy for a driver who wants to maximize her chance of winning a race (against other drivers facing the same MDP). In contrast, traditional approaches would solve this problem by setting a (large) negative reward r elim for an elimination and a small negative reward for each time step before completing the race, and then maximize expectation of rewards using, for example, Qlearning or any other classic algorithm. Different values of r elim would result in very different policies, and while it might be possible to find a good compromise value allowing for a competitive behavior, manually tuning the reward function would be difficult in more complex scenarios.

Consequently, designing and implementing an algorithm for solving MDPs in an RL setting with an SSB utility function also gives us a tool to compute optimal policies for a large class of criteria, including "preference based" ones. A first step towards the design of this algorithm is to give a game-theoretic view of the problem.

A GAME ON POLICIES

When an MDP is fixed, Equations 1 and 2 induce a zerosum two-player symmetric game where the set of strategies coincides with the set of possible policies. The players i ∈ {1, 2} simultaneously choose a strategy πi (pure or mixed). The resulting payoff is then given by ϕ(π 1 , π2 ). As emphasized by Gilbert et al. [2015a], an SSB optimal policy can be found by computing a Nash equilibrium of this game. Indeed, Nash equilibria (π * , π * ) are characterized by ∀π, ϕ(π * , π) ≥ 0.

Gilbert et al.

[2015a] also showed that in this game, a best response to a strategy π is given by a policy maximizing the expectation of cumulated rewards with reward function:

R p π (s) = 1 i Φp π if s = f i ∈ F 0 otherwise
where 1 i is the i th canonical vector, is the transpose operator and Φ is the SSB matrix (i.e., Φ i,j = ϕ(f i , f j )). Put another way, the reward obtained when arriving in the i th final state is given by the i th element of the vector Φp π (the reward is 0 for nonfinal states). Since this defines a standard MDP, there is always a stationary, Markovian and deterministic optimal policy (i.e., the best response to π).

Thus we can restrict ourselves to the finite game with Π s as the set of pure strategies. To summarize, a Nash equilibrium of this game will give an SSB-optimal policy, in the form of a mixed policy π * ∈ Πs .

Example 4 (continued). The game induced by our running example has pure strategies π A , π B , π C for both players, with payoff, e.g., ϕ(π A , π B ) = 14/36 (Example 2). Consider the mixed policy πAAB = (π A |2/3, π B |1/3); Table 1 gives its final state distribution p πAAB = (2/18, 0, 5/18, 10/18, 0, 1/18). The value of response π A to πAAB is given by Similarly, the values of π B , π C are -7/27 and 1/18 respectively, so a best response to πAAB is π A . The second component of vector Φ • p πAAB can be computed as follows. We know that rolling 2 is worse than 16/18 and better than 2/18 of the outcomes of p πAAB , and

p π A • Φ • p πAAB =( 1 6 , 0, 0, 5 6 , 0, 0) (- 16 18 , - 14
ϕ(f i , f j ) is always 1 for i > j. Therefore this component is -1 • 16/18 + 1 • 2/18 = -14/18.
In a finite zero-sum symmetric game, it is well-known that there exists a symmetric Nash Equilibrium (NE). We aim to compute this NE on policies characterized by payoff function ϕ (solving the game).

SOLVING THE GAME

Games in strategic forms can be solved by linear programming [START_REF] Chvátal | [END_REF], unfortunately, here, the game is too large to be solved directly (we remind the reader that the number of pure strategies of the game is equal to the number of deterministic stationary policies of the MDP, which is exponential in the number of states). 4 If the model was known, one could rely on the double oracle algorithm of Gilbert et al. [2015a]. Unfortunately, that approach suffers some drawbacks. Firstly, if the optimal mixture of policies π * is composed of too many policies, its time and space requirements might become prohibitive, as the double oracle algorithm would have to compute and store all polices that are in the support of π * . Secondly, this approach does not generalize to the case where the model is unknown. Thus it can not be used to derive a model-free RL algorithm. We therefore turn to a different approach, based on fictitious play and on a double timescale technique.

LEARNING SETTING

In RL, one typically expects convergence to playing an optimal policy at each step. Since in our case the optimal policy is mixed, this cannot be evaluated at each time step independently. Rather, after any number n of episodes, we consider the final state distribution, defined as

(f 1 |α 1 , . . . , f |F | |α |F | )
, where for all i, α i is the fraction of episodes so far in which the final state was f i .

We define the loss L(p) of a distribution over final states to be the value of its best response against it:

L(p) def = ϕ(p BR , p)
where p BR is the vector of final state frequencies of a policy which maximizes the expectation of cumulated rewards with respect to reward funtion R p (hence a best response to p). Since L(p π * ) = 0 characterizes SSB-optimal policies π * , it is natural to measure the quality of learning by the decrease in loss of the final state distribution so far p n , as n increases. Convergence to an SSB-optimal policy then amounts to L(p n )→0 with n → ∞.

Rephrasing, we expect that, considering the final states reached from the beginning, in retrospect their frequencies are approximately equivalent to those we would have obtained had we played a mixed optimal policy from the start (and more and more exactly as n increases). This corresponds to the standard "on-line" setting of RL, in which success is measured from the start.

FICTITIOUS PLAY

Fictitious Play is an algorithm that only needs a best response procedure to solve a game. The algorithm maintains for each player her mixed policy so far πn , defined as (π 1 |α 1 , . . . , π k |α k ), where for all i, α i is the fraction of episodes so far in which the stationary, deterministic policy π i has been played. At each time step, each player considers that πn perfectly represents the mixed strategy that is 4 Our running example does not illustrate this combinatorial explosion, but it clearly arises, e.g., in the "intransitive grid" and the "race against the clock" (see Section 4).

α i = n • α i /(n + 1) + 1/(n + 1) for π i = π n α i = n • α i /(n + 1) for π i = π n 6 π n+1 = BestResponseTo(π n+1 )
used by the adversary and plays a best response to it. The algorithm converges to a Nash equilibrium of the game (in the sense that L(p n ) converges to 0) when the game is a finite zero-sum game. Fictitious play is represented in Algorithm 1 for a symmetric two-player zero-sum game. As the game is symmetric, we only need to consider the mixed policy so far, πn , of a player playing against herself.

Example 5 (continued). Assume the agent chooses initial strategy π 0 = π B , then after one episode/game we have π1 = (π B |1). Now π A is a best response to π1 (Example 2), hence π 1 = π A . So the agent plays π A during the second episode, and we get π2 = (π A |1/2, π B |1/2). Now it is easy to see that π A is a best response to π2 , so the agent again plays π A and gets π3 = (π A |2/3, π B |1/3). From Example 4 we get that π A is a best response to π3 , and that the loss of π3 is L(p π3 ) = ϕ(π A , π3 ) = 7/54.

SSB Q-LEARNING

This subsection makes a first step towards the adaptation of fictitious play to solve the game induced by an MDP and an SSB utility function. Recall from Section 2.2 that the best responses to the mixed strategy so far, πn , are exactly the optimal policies in the (standard) MDP with reward function R p πn . In other words, best responses can be computed as a function of p πn only. Accordingly, instead of recording πn as such, which involves an exponential number k of pure policies π i in the worst case, as in Algorithm 1, it is enough to record the vector p n = p πn . Then we can rewrite Lines 4-6 of Algorithm 1 as:

p n+1 =n • p n /(n + 1) + p πn /(n + 1) π n+1 =BestResponseTo(p n+1 )
However, in practice, when policy π n is played, one observes only one drawing from p πn , and not the distribution itself. Hence our first adaptation of fictitious play is as given in Algorithm 2. Note that we do not know the model of the MDP, so we use BestResponseTo as an oracle. We find a better solution later in this section.

Example 6 (continued). Let π 0 = π B again. Hence the agent plays π B during one complete episode. If the die rolls 3, we get p 1 = (f 3 |1). Then the best response is computed as one to a strategy which always yields f 3 , and we get π 1 = π A . The agent therefore plays π A during one episode, Algorithm 2: Approximate Fictitious Play Data: MDP M, SSB function ϕ, arbitrary policy π 0 ∈ Π s 1 while True do 2 Play π n for one episode

3 f i = final state reached 4 p n+1 = n • p n /(n + 1) + 1 i /(n + 1) 5 π n+1 = BestResponseTo(p n+1 )
and if the die rolls 4, we get p 2 = (f 3 |1/2, f 4 |1/2), to which a best response is π 2 = π A . If the die then rolls 1, we get p 3 = (f 1 |1/3, f 3 |1/3, f 4 |1/3). In this case the best response is π A , just as in Example 4, but note that it has an estimated value p A • Φ • p 3 = 5/6 • 2/3 -1/6 • 2/3 = 4/9, instead of ϕ(π A , πAAB ) = 7/54 if p πAAB were directly observed.

The following theorem proves that observing only realizations of p πn does not prevent the convergence of (p n ) n∈N to the distribution of an SSB-optimal policy. Theorem 1. In Algorithm 2, L(p n ) tends to 0 as n → ∞.

Proof. Assume an optimal policy π n with respect to R pn is played during the (n+1)-th episode. At the end of the episode, a final state f i is reached and p n+1 is defined by:

p n+1 = n • p n n + 1 + 1 i n + 1 = p n + 1 n + 1 (1 i -p n )
We rewrite this equation in the following way :

p n+1 = p n + 1 n + 1 (p πn -p n + M n+1 ) (3)
where M n+1 = 1 i -p πn is a square integrable martingale difference sequence. Equation 3 is a standard single timescale process with continuous differential inclusion:

.

p(t) ∈ {p π -p(t) : π ∈ Π(p(t))} (4)
where Π(p) denotes the set of optimal policies with respect to reward R p . A similar differential inclusion can be obtained for the standard fictitious play. However, here, as p n is a distribution over final states (not over strategies) and as only a realization of p π is observed, we need to invoke a stochastic approximation argument.

Indeed, the best response correspondence is uppersemicontinuous, with closed and convex values. Hence the existence of at least one solution p(t) through each initial value p(0), which is Lipschitz continuous and defined for all positive times, is guaranteed [Hofbauer, 1995] 0))e -t , and hence L(p(t)) tends to 0 with t → ∞. Thus, the set of final state distributions of the optimal strategies is globally attracting for the best response dynamic, which implies the convergence of the discrete stochastic approximation (3) [Benaïm et al., 2006, Properties 1, 2].

If the model was known, one could use Algorithm 2. At each iteration of the while loop, the corresponding MDP would be solved using dynamic or linear programming to find the best response policy π n+1 . Unfortunately, as the model is unknown, this policy can not be computed directly and has to be learned. The "π n+1 = BestResponseT o(p n+1 )" line should therefore be replaced by, say, a Q-learning phase which converges asymptotically to the desired policy.

Recall that an agent using Q-learning (in a standard MDP) maintains an estimate of Q-values Q(s, a) using:

Q n+1 (s n , a n ) = Q n (s n , a n ) + α n (s n , a n )(r n+1 + max b {Q n (s n+1 , b)} -Q n (s n , a n ))
after taking action a n in state s n , ending up in state s n+1 , and observing the (numerical) reward r n+1 , and where α n (s n , a n ) is the value of the learning rate for (s n , a n ) at timestep n. Note that from now on, for simplicity, we use n to denote the number of time steps (while previously it denoted the number of episodes). Exploitation is realized by choosing a n = arg max a Q n (s n , a) at state s n . Qlearning is known to converge to an optimal policy in the limit, provided the reward function is stationary [START_REF] Watkins | [END_REF].

Hence the idea is to adapt Algorithm 2 to run Q-learning while keeping its target reward function fixed to R pn as when it started, and when it has converged (to an approximate best response to p n ), to update p n using what has been observed in this phase and start a new Q-learning phase. Unfortunately, the number of learning episodes required to learn an optimal or -optimal policy is unknown. Therefore, we want to avoid alternating the Q-learning phase and the reward update phase by using a technique from Borkar [1997] which interleaves both phases successfully using a two-timescale approach. In this approach, the following iterative equations are used concurrently:

p n+1 = p n + β n (1 i -p n ) when f i ∈ F is reached (5) Q n+1 (s n , a n ) = Q n (s n , a n ) + α n (s n , a n )∆ n+1 (6) ∆ n+1 = R pn (s n+1 ) + max b {Q n (s n+1 , b)} -Q n (s n , a n ) ∀s, a,              ∞ n=0 α n (s, a) = ∞ and ∞ n=0 β n = ∞ β n /α n (s, a) → 0 ∞ n=0 α n (s, a) 2 + β 2 n < ∞              (7)
The intuition is that Q n evolves more quickly than p n , giving time for Q n , and hence the best response, to adapt to changes in the target reward R pn .

Write η p (resp. η s,a ) for the number of times vector p n (resp. state action-pair (s, a)) has been updated (resp. experienced). In our case, β n is 1/(η p + 1) since Equation 5tracks the final state distribution so far, so we can set α n to 1/(η s,a + 1) 2/3 for instance. Note that this example works because p n is updated after at most T max steps, bounding the number of updates of α n (for any (s, a)) between two consecutive updates of β n . With these conditions, the two recursive equations on p n and Q n form a two-timescale stochastic approximation iteration, where p n is on a slower timescale than Q n . The following example gives the intuition behind the two time-scale technique; the convergence proof under Condition ( 7) is given in Theorem 2.

Example 7. Suppose that at some point the Q-values of a A and a BC in s 0 are 0.1 and 0.2, respectively, and that those of a B ,a C in s BC are 0.3,-0.7. Assume moreover p n =(1/5, 0, 1/5, 1/5, 1/5, 1/5), n = 80 (with 20 episodes of length 1 and 30 of length 2, hence η p 0.02) η s0,aBC = 30 (α n (s 0 , a BC ) 0.1) and η s BC ,a B = 20 (α n (s BC , a B ) 0.13).

Hence the agent will choose a BC in s 0 , reach s BC and update Q(s 0 , a BC ) to 0.2 + 0.1 • (0 + 0.3 -0.2) = 0.21. In s BC she will choose a B and, if she rolls 3, update

Q(s BC , a B ) to 0.3 + 0.13 • (1 3 • Φ • p n + 0 - 0.3) = 0.3 + 0.13 • (-2/5 -0.3) 0.
2, and p n to (10/51, 0, 11/51, 10/51, 10/51, 10/51).

Note that Q-values have evolved significantly more than p n (and hence than the target numerical reward R pn ).

HANDLING EXPLORATION

In order for Q-learning to converge, one needs to ensure that all state-action pairs are performed infinitely often. Usually, this exploration is guaranteed through some randomization, using an -greedy strategy, for instance. We present in this subsection an exploration strategy called Episodic--Greedy (EG for short) that guarantees that we converge to an -optimal SSB strategy using the algorithm described by Equations 5 and 6. During learning, an

= R pn (s n+1 ) 4 Q n+1 (s n , a n ) = Q n (s n , a n ) + α n (s n , a n )(r n+1 + max b {Q n (s n+1 , b)} -Q n (s n , a n )) 5 if s n+1 = f i ∈ F and exploration is off then 6 p n+1 = p n + 1 ηp+1 (1 i -p n )
episode will be generated using either the current best policy (defined by the Q-values), with probability (1 -), or the uniformly random policy, which we denote by π U , with probability . If an episode is generated using π U (i.e., the agent is exploring), then the update of Equation 5is not performed at the end of the episode. This guarantees that the convergence of p n is not biased by the exploration strategy.

The final proposed algorithm is presented in Algorithm 3.

We are now ready to prove the convergence of SSB Qlearning to an -optimal policy through two theorems.

Theorem 2. Under Conditions (7) with β n = 1/(η p + 1), in Algorithm 3, L(p n ) tends to 0 almost surely as n → ∞.

Proof. The idea is that p n can be viewed as quasi-static compared to Q n . Indeed, let π n be the greedy policy given by Q n . We can rewrite the equations as:

p n+1 = p n + α n ( n + M n ) (8) Q n+1 (s n , a n ) = Q n (s n , a n ) + α n (s n , a n )(T (Q n )(s n , a n ) -Q n (s n , a n ) + M n+1 )
where n = βn αn (p πnp n ) and M n = βn αn (1 i -p πn )

T (Q n )(s, a) = s P(s |s, a)(R pn (s )+max b {Q n (s , b)}) M n+1 =R pn (s n+1 ) + max b {Q n (s n+1 , b)} -T (Q n )(s n , a n )
Clearly n → 0 almost surely (||p πn || and ||p n || are bounded). Then (p n , Q n ) will converge to the internally chain transitive invariant set of the ODE [Borkar, 2008]:

.

p(t) = 0 . Q(t) = T (Q(t)) -Q(t)
Let Q * (p n ) denote the optimal Q-value function for reward function R pn . Therefore Q n -Q * (p n ) → 0 almost surely, which entails that (p n , π n ) converges to the set (p, π * (p)) with π * (p) a best response to p. We then rewrite (5) to:

p n+1 =p n + β n (p π * (pn) -p n + (p πn -p π * (pn) ) + (1 i -p πn ))
As p πnp π * (pn) → 0 almost surely, the asymptotic behavior is the same as in Theorem 1. Thus the loss of p n converges to 0 with n → ∞.

The result of Theorem 2 uses the fact that with the EG exploration strategy, exploration has no impact on p n . The drawback of this strategy is that p n does not truly represent the frequencies with which each final state has been obtained. If we let p real n represent the vector of true frequencies with which each final state has been obtained, then following theorem proves that p real n converges to the final state frequencies of an -optimal SSB-policy. Theorem 3. Under Conditions (7) with β n = 1/(η p + 1), when Algorithm 3 is run, p real n converges to the final state distribution of an -optimal SSB-policy almost surely.

Proof. Let denote the parameter of EG exploration and let p real = lim n→∞ p real n . Then asymptotically we have

p real = (1 -)p * + p π U
where p * is the final state distribution of an optimal SSBpolicy. Thus for any policy π:

ϕ(p real , p π ) = (1 -)ϕ(p * , p π ) ≥0 + ϕ(p π U , p π ) ≥-ϕmax ≥ -ϕ max with ϕ max = max f,f ϕ(f, f ).
Hence, with = /ϕ max , Algorithm 3 converges to an -optimal SSB strategy.

PROOF OF CONCEPT

Although SSB encompasses many different criteria, we focus here on the probabilistic dominance criterion as it is an important case for which no model-free algorithm that is provably correct has been proposed. We plot here the results of four experiments: "sequential Gardner dice", "who wants to be a millionaire", "intransitive grid" and "race against the clock" using the probabilistic dominance criterion (hence values are all in [-1, 1]). For all runs, 10, 000, 000 steps were performed in the MDP, was set to 0.1, and α n was set to 1/(η s,a + 1) 11/20 .

Gardner Dice. We first present the results on Gardner's dice problem as formalized in Example 1. Figure 1(a) shows the evolution of the frequencies (f A , f B , f C ) with which each die has been played for a representative run. The optimal frequency vector p * = (3/13, 3/13, 7/13) is shown as a green dot and the same vector biased by exploration p * = (1 -) * p * + * p π U by a red dot; the vector (f A , f B , f C ) tends towards p * , drawing triangles of decreasing surface around p * . Figure 1(b) presents the evolution of the Q-values of the three actions, a A , a B and a C . One can see that the best die alternates between the three dice and that max{Q A , Q B , Q C } tends towards 0. (The best response is always deterministic and so must be one of π A , π B , π C . However at convergence its value has to be 0.)

Who wants to be a millionaire. In this popular television game show, a contestant answers 15 multiple-choice questions (with four possible answers) of increasing difficulty, for increasingly large sums, roughly doubling the pot each question. At each time step, the contestant may decide to walk away with the money currently won. If she answers incorrectly, then all winnings are lost except what has been earned at a "guarantee point" (questions 5 and 10). The player is allowed 3 lifelines (50:50, removing two of the choices, ask the audience and call a friend for suggestions); each can only be used once. We used the first model of the Spanish 2003 version of the game presented by Perea and Puerto [2007]. The probability of answering correctly is a function of the question's number and increased by the lifelines used (if any).

Intransitive grid. In this domain we study an episodic grid MDP containing 9 states. The agent always starts an episode in the bottom-right corner of the grid. Three terminal states, f 1 , f 2 and f 3 can be attained at the three other corners of the grid. The agent can only go left and up. With a probability of 0.2, the agent makes a mistake and goes in the wrong direction. The preference relation between the final states is the following:

f 1 f 2 f 3 f 1 .
Race against the clock. Lastly, we discuss the domain discussed in Example 3. 5 The racing circuit is represented by 6 physical positions {p 1 , . . . , p 6 } plus one, p el , representing elimination. A state of the MDP is a triple (p, s, t) giving the current position p ∈ {p 1 , . . . , p 6 , p el }, the current speed of the car s ∈ {Slow, M edium, F ast} and the current time t ∈ N. In each state the agent can decide between 3 actions: accelerating, decelerating or keeping the same speed. At each time step t, the probability of running off the track is a function of s t , p t and a t , taking into account both the speed of the car and the difficulty of the current part of the circuit. Finally, the time spent between two positions decreases stochastically with the current speed.

Results. Figure 1(c)-(f) shows the evolution of L(p real n ) and L(p n ) (i.e., the values of the optimal policies regarding reward functions defined by Φp real n , Φp n ) for each domain. The results are averaged over 20 runs. As expected the value of L(p n ) tends towards 0 as the number of learning steps increase. The value of L(p real n ) decreases and is much lower than (= 0.1 in the experiments) on all figures.

Finally, in Table 2, we compare the "race against the clock" results obtained by our SSB Q-learning algorithm to the results obtained by a standard Q-learning algorithm launched with three different reward functions R 1 , R 2 and R 3 . For each reward function, a penalty of value -t is received by the agent each time the circuit is completed in time t. For reward function R i , an elimination results in a penalty of value r i elem ∈ {-10, -25, -40}. For each algorithm we give the frequency of elimination f elem and the average time T cc of circuit completion. The final columns show the probabilities with which the SSB Q-learning agent would beat (P ) and at least tie with (P ) each Q-learning agent.

As expected, f elem decreases and T cc increases with the penalty value of an elimination. For a Q-learning agent, this penalty would have to be tuned to give the best compromise. The SSB Q-learning agent does not face this problem and the last two columns show that this agent is more likely to produce a preferred episode.

Conclusion

Skew-Symmetric Bilinear (SSB) utility is a useful general decision model that encompasses many decision criteria (e.g., EU, threshold probability, probabilistic dominance, etc.). We designed a model-free reinforcement learning algorithm to compute an epsilon SSB-optimal policy and provided experimental results. Our current work can be extended in several natural ways. For instance, it would be interesting to tackle non-episodic problems. Another direction is to use more elaborate RL algorithms than Q-learning for computing best responses.

Our work is currently being applied in an industrial context with good results. An automated information extraction (IE) treatment chain is modelled as an MDP, and improved using a reward function balancing extraction quality and treatment time. SSB utility theory is used to formalise qualitative preferences expressed by human operators on the output of the treatments.

  π 1 |α 1 , . . . , π k |α k ) with 5

Algorithm 3 :

 3 SSB Q-learning Data: MDP M, SSB function ϕ 1 while True do 2 Choose a n using the EG exploration strategy 3 Play a n , observe s n+1 , and let r n+1

Figure 1 :

 1 Figure 1: For the Gardner dice domain, (a) convergence of the policy in the space of die frequencies, (b) evolution of Q-values, (c) evolution of loss; For (d) Who wants to be a millionnaire, (e) Intransitive grid and (f) Race against the clock, evolution of loss.

Table 1 :

 1 ,a 0 ,s 1 ,. . . ,s t-1 ,a t-1 ,s t ), starting in s 0 and ending in a final state s t ∈ F. When the current episode ends, a new episode starts in state s 0 . We further assume that the length of an episode is upper-bounded by a constant T max ∈ N.A policy π at horizon T indicates which action to perform in each nonfinal state for each time step t < T . A policy is Markovian if the action depends only on the current state and timestep (otherwise it may depend on all state-action pairs encoutered so far); deterministic if it prescribes exactly one action, or randomized if it prescribes a probability distribution over actions; stationary if the action prescribed does not depend on the timestep. We write Π s for the set of Markovian stationary deterministic policies.

		Probabilities for Each Gardner Die
		1	2	3	4	5	6
	p A 1/6	0	0	5/6	0	0
	p B	0	0	5/6	0	0	1/6
	p C	0	1/2	0	0	1/2	0
	(s 0

Importantly, given a set Π = {π 1 , π 2 , . . .} of policies, we define an enlarged set Π of policies, that denotes the set consisting of mixtures of policies, i.e., Π = {π = (π 1 |α 1 , π 2 |α 2 , . . .) : i

  As ζ(t) maximizes ϕ(., p(t)), the first term is null[Mas- Colell et al., 1995, pp. 964-965], and by linearity:

	d dt	L(p(t)) = ϕ(ζ(t), = ϕ(p(t),	. p(t))	. p(t))
			=0
		= ϕ(p(t),	. p(t)) + ϕ(p(t), p(t))
			=0
		= ϕ(p(t), ζ(t)) = -ϕ(ζ(t), p(t))
	as Φ is skew-symmetric (hence ϕ(p(t), p(t)) = ϕ( . p(t), . p(t)) = 0, and ϕ(p(t), ζ(t)) = -ϕ(ζ(t), p(t))).
	Thus, d dt L(p(t)) = -L(p(t)), L(p(t)) = L(p(

. Let p(t) be a solution of inclusion 4 and ζ(t) = p(t) + . p(t) = p π for a best response π ∈ Π(p(t)). By definition of L, we have L(p(t))= ϕ(ζ(t), p(t)), and by the envelope theorem: d dt L(p(t)) = ∂ϕ(ζ(t), p(t)) ∂ζ . ζ(t) + ∂ϕ(ζ(t), p(t)) ∂p . p(t) . p(t)) + ϕ( . p(t),

Table 2 :

 2 Comparisons of SSB Q-Learning with Three Q-Learning Algorithms (Results Averaged Over 20 Runs). with R 1 0.48 4.34 0.37 0.64 Q-learning with R 2 0.31 7.14 0.48 0.66 Q-learning with R 3 0.26 9.09 0.54 0.66

		f elem T cc	P	P
	SSB Q-learning	0.41 5.29	-	-
	Q-learning

Not to be confused with the notion of randomized policies.

The complete description is given as supplementary material at hugogilbert.pythonanywhere.com.

Acknowledgement

Work supported by the French National Research Agency through the Idex Sorbonne Universités, ELICIT project under grant ANR-11-IDEX-0004-02.


Tu-Hoa Pham

Abderrahmane Kheddar

Ammar Qammaz

Antonis A Argyros

Towards Force Sensing from Vision: Observing Hand-Object Interactions to Infer Manipulation Forces

We present a novel, non-intrusive approach for estimating contact forces during hand-object interactions relying solely on visual input provided by a single RGB-D camera. We consider a manipulated object with known geometrical and physical properties. First, we rely on model-based visual tracking to estimate the object's pose together with that of the hand manipulating it throughout the motion. Following this, we compute the object's first and second order kinematics using a new class of numerical differentiation operators. The estimated kinematics is then instantly fed into a second-order cone program that returns a minimal force distribution explaining the observed motion. However, humans typically apply more forces than mechanically required when manipulating objects. Thus, we complete our estimation method by learning these excessive forces and their distribution among the fingers in contact. We provide a full validity analysis of the proposed method by evaluating it based on ground truth data from additional sensors such as accelerometers, gyroscopes and pressure sensors. Experimental results show that force sensing from vision (FSV) is indeed feasible.

Introduction

Reliably capturing and reproducing human haptic interaction with surrounding objects by means of a cheap and simple set-up (e.g., a single RGB-D camera) would open considerable possibilities in computer vision, robotics, graphics, and rehabilitation. Computer vision research has resulted in several successful methods for capturing motion information. A challenging question is: to what extent can vision also capture haptic interaction? The latter is key for learning and understanding tasks, such as holding an object, pushing a chair or table, as well as enabling its reproduction from either virtual characters or physical (e.g., robotic) embodiments.

Contact forces are usually measured by means of hap-Figure 1: Using a single RGB-D camera, we track markerless hand-object manipulation tasks and estimate with high accuracy contact forces that are applied by human grasping throughout the motion.

tic technologies such as force transducers. The main drawback of such technologies is that they are obtrusive. Computer vision techniques would therefore be an ideal alternative to circumvent this issue. Yet, is it possible to estimate forces from visual observation? There is evidence that haptic perception can be induced through illusion and substitution dominated by vision, e.g. [START_REF] Lécuyer | Pseudo-haptic feedback: Can isometric input devices simulate force feedback[END_REF]. We aim at exploring computer vision to infer the forces exerted by humans on surrounding objects. In particular, we consider hand-object grasping and manipulation. The problem is extremely complex. Indeed, establishing that a hand-object contact has occurred is difficult because of occlusions and tracking inaccuracies. Nevertheless, the detection of events like an object being lifted or discontinuities in body motion may provide useful hints towards disambiguating discrete events. Additionally, even if contact positions can be determined efficiently, the estimation of the applied forces is still challenging because of the inherent multiplicity of solutions.

We demonstrate that, by solely using computer vision, it is possible to compute interaction forces occurring in handobject manipulation scenarios where object properties such as shape, contact friction, mass and inertia are known, along with human hand geometry. First, we monitor both the hand and the object motions by using model-based 3D tracking (other visual tracking techniques can also be used if they meet performance requirements). From the tracking data, we estimate hand-object contact points through proximity detection. Algebraic filtering computes the object's kinematics, i.e. velocity and acceleration. Contact force distributions explaining the kinematic observations are then resolved using conic optimization. When manipulating objects, humans typically apply more (internal) forces than what is required from the Newton-Euler dynamics. Thus, we improve our estimation method by using neural networks to learn the amount and distribution of these internal forces among the fingers in contact. The experimental results obtained on datasets annotated with sensors' groundtruth show the potential of the proposed method to infer hand-object contact forces that are both physically realistic and in agreement with the actual forces exerted by humans during grasping. To the best of our knowledge, this is the first time that this problem is addressed and solved based solely on markerless visual observations.

Related Work

Applications of force sensing in robotics, virtual reality and human physiological studies result in challenging requirements for transduction technologies [START_REF] Cutkosky | Springer Handbook of Robotics, chapter Force and Tactile Sensors[END_REF][START_REF] Dahiya | Tactile sensing -from humans to humanoids[END_REF]. Common drawbacks of mechatronic force sensing devices reside in their extensive need for calibration, time-varying accuracy (e.g., hysteresis) and cost. Additionally, they are required to be mounted onto the manipulated objects, thus impacting their shape or other physical properties, or onto the operator (e.g., force sensing gloves), thus obstructing the human haptic senses and limiting the natural range of motion.

An alternative to pressure sensors proposed in [START_REF] Mascaro | Photoplethysmograph fingernail sensors for measuring finger forces without haptic obstruction[END_REF] consists in instrumenting fingers with miniature LEDs and photodetectors to measure changes in fingernails coloration, that are then correlated to the touch force applied at fingertips. Later, this technology evolved to predict normal and shear forces, and even changes in posture, that appear to have different blood volume patterns [START_REF] Mascaro | Measurement of finger posture and three-axis fingertip touch force using fingernail sensors[END_REF]. In [START_REF] Sun | Predicting fingertip forces by imaging coloration changes in the fingernail and surrounding skin[END_REF], fingernail color and surrounding skin changes are monitored and processed using an external camera system to estimate contact fingertip forces [START_REF] Sun | Estimation of fingertip force direction with computer vision[END_REF][START_REF] Grieve | Force prediction by fingernail imaging using active appearance models[END_REF][START_REF] Urban | Computing grip force and torque from finger nail images using gaussian processes[END_REF]. Conversely, computer graphics models were developed to simulate fingertip appearance changes based on simulated forces [START_REF] Andrews | Data-driven fingertip appearance for interactive hand simulation[END_REF]. This result illustrates that computer vision can be used effectively to measure touch forces on fingertips. This approach is, however, limited to fingertip contacts and requires extensive calibration for each individual user. It is also limited by the necessity of having fingernails visible at all time and at high resolution, hence requiring appropriately mounted miniature cameras.

Used in conjunction with force sensing technologies, motion tracking can provide information on body dynam-ics to explain how humans interact with their environment. A setup combining marker-based motion capture and mechanical force sensors was used in [START_REF] Kry | Interaction capture and synthesis[END_REF] to estimate hand joint compliance and synthesize interaction animations. Although precise, marker-based motion capture is invasive and hardly applicable to real-world contexts.

Alternative tracking techniques were therefore developed to circumvent markers' intrusiveness. The markerless vision-based tracking of a hand was first treated in [START_REF] Rehg | Visual tracking of high dof articulated structures: an application to human hand tracking[END_REF] and has lately received renewed attention [START_REF] Oikonomidis | Efficient model-based 3d tracking of hand articulations using kinect[END_REF][START_REF] De La Gorce | Model-based 3d hand pose estimation from monocular video[END_REF][START_REF] Keskin | Hand pose estimation and hand shape classification using multi-layered randomized decision forests[END_REF][START_REF] Qian | Realtime and robust hand tracking from depth[END_REF][START_REF] Tang | Latent regression forest: Structured estimation of 3d articulated hand posture[END_REF][START_REF] Tompson | Real-time continuous pose recovery of human hands using convolutional networks[END_REF]. While these works show impressive results for tracking a hand in isolation, they are not applicable to our scenario, mostly due to severe mutual occlusions from strong handobject interaction, resulting in missing observations of both. For the 3D tracking of a hand in interaction with object(s), the existing methods can be classified into mostly bottom up, mostly top-down or hybrid, depending on the relative strength of the discriminative and the generative processes involved. The method in [START_REF] Romero | Nonparametric hand pose estimation with object context[END_REF] performed 3D tracking of a hand-object interaction by synthesizing its appearance based on a non-parametric discriminative model and by treating hand pose estimation as a classification problem.

In general, top-down methods do not treat occlusions as a distractor but rather as a source of information [START_REF] Oikonomidis | Full dof tracking of a hand interacting with an object by modeling occlusions and physical constraints[END_REF][START_REF] Ballan | Motion capture of hands in action using discriminative salient points[END_REF][START_REF] Oikonomidis | Tracking the articulated motion of two strongly interacting hands[END_REF][START_REF] Schmidt | Dart: Dense articulated real-time tracking[END_REF]. A generative process creates hypotheses of a joint hand-object model that are then evaluated against actual observations by solving a multi-parameter optimization problem. In the representative approach presented in [START_REF] Oikonomidis | Full dof tracking of a hand interacting with an object by modeling occlusions and physical constraints[END_REF] the availability of a joint 3D model facilitates the consideration of possible hand-object occlusions. Additionally, it makes possible the incorporation of strong priors in the optimization process, e.g., the fact that two different objects cannot share the same physical space. In this manner, even more complex problems such as tracking two strongly interacting hands [START_REF] Oikonomidis | Tracking the articulated motion of two strongly interacting hands[END_REF] or tracking two hands in interaction with an object [START_REF] Ballan | Motion capture of hands in action using discriminative salient points[END_REF] have been successfully addressed. To cope with hand-object occlusions, even stronger priors resulting from a physics-based simulation of the scene have also been incorporated [START_REF] Kyriazis | Physically plausible 3d scene tracking: The single actor hypothesis[END_REF].

The scalability issues of generative methods are dealt with in [START_REF] Kyriazis | Scalable 3d tracking of multiple interacting objects[END_REF], where it is shown that the 3D tracking of the interaction of two hands with several rigid objects can be performed accurately and efficiently. The approach taken in our work to support visual tracking of hand-object interactions is based on a variant of the method presented in [START_REF] Kyriazis | Scalable 3d tracking of multiple interacting objects[END_REF]. A representative hybrid method is proposed in [START_REF] Hamer | Tracking a hand manipulating an object[END_REF]. In that work, the 3D tracking of a hand manipulating an object is achieved by using bottom-up evidence to directly guide (rather than simply support) the formulation of hand pose hypotheses that is performed by generative means.

An inspiring use of motion tracking for force estimation is presented in [START_REF] Brubaker | Estimating Contact Dynamics[END_REF], where whole body contact forces and internal joint torques are estimated by solving an optimiza-tion problem linking contact dynamics and human kinematics. The goal of the work in [START_REF] Ye | Synthesis of detailed hand manipulations using contact sampling[END_REF] is the realistic synthesis of detailed hand manipulations. Towards this direction, motion capture data is used to formulate the synthesis problem taking into account contact points and contact forces. This method would benefit a lot from the unobtrusive visionbased estimation of actual contact forces we propose. Also aiming at realistic synthesis and motion reconstruction, the work in [START_REF] Zhao | Robust realtime physics-based motion control for human grasping[END_REF][START_REF] Wang | Video-based hand manipulation capture through composite motion control[END_REF] estimates the motion of a hand in interaction with an object together with contact points and exerted forces. Though this method successfully tracks challenging manipulation scenarios, the manipulation forces are only constructed to be compatible with visual observations, without aiming at matching the forces humans actually apply, as happens with the method we propose in this work.

Force Sensing from Vision (FSV)

Let S be a rigid body with mass m and inertia matrix J relative to its center of mass C. For any element e of its contacting environment (e.g. human hand, table), let µ e denote the corresponding Coulomb friction coefficient. We assume these quantities to be known as they can be obtained from the object of interest's CAD model or existing identification techniques [START_REF] Schedlinski | A survay of current inertia parameter identification methods[END_REF]. Interestingly, it has been shown that aspects of such information (e.g., mass distribution) can also be estimated by visual means [START_REF] Mkhitaryan | Visual estimation of object density distribution through observation of its impulse response[END_REF]. We then consider a scenario where S is grasped and manipulated by a hand, with possible contacts with the environment. Observing the scene with a sole RGB-D camera, which we suppose calibrated with the vertical direction known, our goal is to estimate the interaction forces between S and the user's hand, and between S and its environment when applicable. We decompose FSV into four subproblems as follows:

1. Track S and the hand and perform, for each step, vision-based proximity or collision detection to identify contacting fingers and corresponding contact points (Section 3.1). 2. Let X i be the estimated pose for S at instant i. Based on sequence (X i ) i∈[0,N ] , estimate for each frame the body's first and second-order kinematics, i.e. translational (resp. rotational) velocity v i (resp. ω i ) and acceleration a i (resp. α i ) (Section 3.2). 3. Compute a force distribution explaining the object's state computed at step 2 following Newton-Euler's laws of motion and Coulomb's friction model, using the contact points identified at step 1 (Section 3.4). 4. Learn and reproduce human internal force distributions among the fingers in contact (Section 3.6).

Each of these subproblems presents a number of challenges. First, the observation of manipulation tasks may be subject to hand-object occlusions. To overcome this issue, we address step 1 by means of model-based tracking as inspired by [START_REF] Kyriazis | Scalable 3d tracking of multiple interacting objects[END_REF]. Second, the limited camera acquisition frequency along with tracking errors can make the differentiation process of step 2 unstable. We tackle this issue by estimating derivatives using algebraic filtering derived from [START_REF] Mboup | Numerical differentiation with annihilators in noisy environment[END_REF]. Algebraic filtering was chosen for the sake of robustness, as it relies on no statistical assumption on the signal's noise. We then address step 3 by computing minimal force closure distributions as solutions of a secondorder cone program. Finally, step 4 stems from the fact that in contrast with [START_REF] Sun | Predicting fingertip forces by imaging coloration changes in the fingernail and surrounding skin[END_REF] where multiple photodetectors monitor each fingernail's blood flow individually, such microscopic features cannot be observed by a single RGB-D camera observing the whole scene. The object may indeed be grasped with more or less intensity without this being visible at a macroscopic scale. We tackle this statical indeterminacy with machine learning on usual human grasping practices.

Hand-object tracking

Our approach requires a good 3D pose estimate of the manipulated object together with that of the user's hand. To achieve this, we rely on a variant of the method proposed in [START_REF] Kyriazis | Scalable 3d tracking of multiple interacting objects[END_REF] that is tailored to our needs. In [START_REF] Kyriazis | Scalable 3d tracking of multiple interacting objects[END_REF], the model-based hand-object 3D tracking is formulated as an optimization problem, which seeks out the 3D object(s) pose and hand configuration that minimizes the discrepancy between hypotheses and actual observations. The optimization problem is solved based on PSO [START_REF] Eberhart | Swarm Intelligence (The Morgan Kaufmann Series in Evolutionary Computation)[END_REF].

Since this method estimates the generalized pose of a hand interacting with an object, it is straightforward to compute the 3D positions of the estimated fingertips in relation to the object's surface (i.e., contact points). Still, in our implementation of [START_REF] Kyriazis | Scalable 3d tracking of multiple interacting objects[END_REF], we have incorporated one important modification. The original 3D hand-object tracking framework provides solutions that are compatible with visual observations and are physically plausible in the sense that the hand and the object do not share the same physical space (i.e., the hand does not penetrate the modeled volume of the object). However, occluded fingers may have different poses that respect the above constraints, making the estimation of contact points an under-constrained problem. To overcome this issue, we assume that contact points do not change significantly when they cannot be observed. Time and space coherency is thus enforced by penalizing solutions in which hidden contact points are far from their last observed position.

Numerical differentiation for kinematics

In theory, velocity and acceleration can be estimated by numerical differentiation of poses obtained from tracking. However, this process is highly dependent on two factors: (a) the acquisition frequency of the RGB-D frames, and (b) the quality of the motion tracking. First, even a perfect tracking would result in poor velocity and acceleration estimates if performed over time steps far apart from each other, also depending on the way the hand moves. However, this is not a freely controllable parameter, as most commercial RGB-D cameras offer acquisition frame-rates capped between 30 and 60 fps. We present our results on a 30 fps SoftKinetic DepthSense 325 camera. Second, acceleration profiles occurring in manipulation tasks are naturally spiky (see for example Fig. 4). Therefore, numerical differentiation is challenging in that while the number of samples used for each derivative estimate must be sufficient to alleviate tracking errors, it must also be kept minimal to discern the sudden variations that acceleration profiles are subject to.

As an alternative to existing numerical differentiation methods, algebraic parameter estimation approaches [START_REF] Fliess | An algebraic framework for linear identification[END_REF] led to a new class of derivative estimators called algebraic numerical differentiators [START_REF] Mboup | Numerical differentiation with annihilators in noisy environment[END_REF]. The tracking errors resulting from the employed model-based tracking framework seem to follow a Gaussian distribution, yet they are not independent of one another, which rules out the white noise formalism. Subsequently, and in order to keep the kinematics estimation process unbiased by the use of a particular tracking method, we implement the so-called minimal (κ, µ) algebraic numerical differentiators, which do not assume prior knowledge of the signal errors' statistical properties.

From kinematics to dynamics

We suppose the rigid body S subject to n d non-contact forces (F d k ) k∈[1,n d ] applied at points (P d k ) k∈{1,...,n d } (e.g., gravitation, electromagnetism). We consider them fully known based on the object's properties. We seek to estimate n c contact forces (F c k ) k∈ [1,nc] applied at contact points with the hand or the environment (P c k ) k∈ [1,nc] that are obtained from tracking (Section 3.1). Using the object's kinematics as estimated in Section 3.2, its motion is governed by Newton-Euler equations. Therefore, the resulting net force F c and torque τ c due to the contact forces are such that:

F c = ma -F d τ c = J q • α + ω × (J q • ω) -τ d , (1)
with F d and τ d the resulting force and torque due to noncontact forces, and J q the inertia matrix with orientation q.

The contact forces are subject to friction, which we model using Coulomb's law. Let n k be the unit contact normal oriented inwards S at contact point P c k . Let then t x k and t y k be two unit vectors orthogonal to each other and to the normal n k , thus defining the tangent plane. Each contact force F c k is decomposed as follows:

F c k = f k n k + g k t x k + h k t y k , (2)
With µ k the friction coefficient at P c k , Coulomb's law reads:

g k t x k + h k t y k 2 ≤ µ k f k , (3)
which is a strict equality in the case of dynamic friction.

Nominal forces from cone programming

We address the estimation of the minimal contact forces responsible for the observed motion (i.e., nominal forces) as a second-order cone program (SOCP) [START_REF] Lobo | Applications of second-order cone programming[END_REF][START_REF] Boyd | Convex Optimization[END_REF][START_REF] Boyd | Fast computation of optimal contact forces[END_REF]:

min C(x) = 1 2 x T Px + r T x s. t. A j x + b j 2 ≤ c T x + d j , j = 1, . . . , m
and Fx = g.

(

) 4
As we track the object and the user's hand, we can determine, at each timeframe, newly established and broken contacts, and also those that remain still and those that slide. Therefore, we are explicitly considering static and kinetic (i.e. dynamic) friction in the constraints formulation. With n c,s and n c,k the respective numbers of friction forces and n c their sum, we choose the optimization vector as follows:

x = (f 1 , g 1 , h 1 , . . . , f nc,s , g nc,s , h nc,s , f nc,s+1 , f nc,s+2 , . . . , f nc,s+n c,k ) T (5)
The SOCP formulation in Eq. ( 4) then allows the direct handling of Coulomb static friction as inequality constraints and kinetic friction as equality constraints. Moreover, having each normal vector n k oriented inwards S, n c additional positivity constraints are added such that f k ≥ 0. Equality constraints enforcing that the resulting contact force distribution explains the observed kinematics stem from Newton-Euler's equations, as combining Eq. ( 1) with contact force expressions from Eq. ( 2) directly yields linear equations in x. We complete the SOCP with the objective function:

C SoS (x) = k∈F f 2 k + g 2 k + h 2 k = k∈F F c k 2 2 , ( 6
)
where F is the set of contacting fingers. As stated earlier in Section 3, there exists an infinity of possible force distributions for a given kinematics and set of contact points.

Using the contact forces' sum of squares as an indicator on the grasp intensity (i.e., its L 2 norm), the objective function C SoS allows to search for the optimal grasp in that sense. Numerical resolution is performed using CVXOPT [START_REF] Andersen | Cvxopt: A python package for convex optimization[END_REF].

Reproducing human grasping forces

Humans do not manipulate objects using nominal closures (i.e. minimal grasp forces). They tend to "over-grasp" and produce workless internal forces, i.e. firmer grasps than mechanically required. This human grasp property is described by considering finger forces as two sets [START_REF] Kerr | Analysis of multifingered hands[END_REF][START_REF] Yoshikawa | Manipulating and grasping forces in manipulation by multifingered robot hands[END_REF]: nominal forces responsible for the object's motion, and internal forces that secure the object through a firm grip but do not affect the object's kinematics as they cancel each other out [START_REF] Mason | Robot Hands and the Mechanics of Manipulation[END_REF][START_REF] Murray | A Mathematical Introduction to Robotic Manipulation[END_REF]. Studies showed that humans apply internal forces to prevent slip [START_REF] Johansson | Roles of glabrous skin receptors and sensorimotor memory in automatic control of precision grip when lifting rougher or more slippery objects[END_REF][START_REF] Flanagan | Hand movements[END_REF] to avoid muscle fatigue or damaging fragile objects [START_REF] Gorniak | Manipulation of a fragile object[END_REF][START_REF] Park | Optimality versus variability: effect of fatigue in multi-finger redundant tasks[END_REF].

(i) k )k∈F SOCP f (n) k , f (i)
We extend the formulation of the SOCP to this decomposition and learn from tactile sensors' measurements how humans apply internal forces when manipulating objects. Each finger force

F k is decomposed into a nominal com- ponent F (n)
k and an internal component

F (i) k : F k = F (n) k + F (i) k with F (n) k = f (n) k n k + g (n) k t x k + h (n) k t y k F (i) k = f (i) k n k + g (i) k t x k + h (i) k t y k . (7)
It is to mention that although both forces are decomposed along the same contact frame (n k , t x k , t y k ) as in Eq. ( 2), nothing constraints them to be similarly oriented. We subsequently gather the nominal and internal components into a new optimization vector x as in Eq. [START_REF] Boyd | Fast computation of optimal contact forces[END_REF].

By definition, nominal forces are responsible for the object's motion through the Newton-Euler equations while internal forces are neutral regarding its state of equilibrium:

       k∈F F (n) k = F c , k∈F --→ CP k × F (n) k = τ c k∈F F (i) k = 0, k∈F --→ CP k × F (i) k = 0. (8)
Equation [START_REF] Dahiya | Tactile sensing -from humans to humanoids[END_REF] provides a new set of constraints that we integrate into the SOCP of Section 3.4. Ensuring that the resulting distribution still obeys Coulomb's law of friction, we finally compute the distribution of nominal and internal forces that best match the tactile sensors' measurements ( f k ) k∈F , using a new objective function:

C dist (x) = k∈F F (n) k 2 2 + f (n) k +f (i) k -f k 2 . (9)
The reason why we do not directly identify internal forces as the differences between the measurements f k and the minimal forces resulting from the initial SOCP of Section 3.4 is that possible sensor measurement errors may lead them not to compensate each other. By integrating their computation into the SOCP, we ensure that the resulting internal forces f (i) k bridge the gap between f (n) k and measurements f k without perturbing the object's observed kinematics. We illustrate the decomposition process in Fig. 2(a).

Learning internal force distributions

Recent studies attempted to build mathematical models correlating grasp forces to kinematic data, yet limited to cyclic movement patterns and two-finger grasps [START_REF] Gao | Internal forces during object manipulation[END_REF][START_REF] Slota | Grip forces during object manipulation: experiment, mathematical model, and validation[END_REF], hence concealing the issue of statical indeterminacy. In contrast, our approach learns how humans apply internal forces by means of artificial neural networks (ANN). We first construct an experimental dataset by having human operators manipulate an instrumented box (see Section 4) over tasks such as pick-and-place, lift and release, rotations, and unguided compositions of these. Experiments were conducted over a pool of three female, including one left-handed, and three male operators using their preferred hand on different contact and object mass configurations. Executing 160 manipulation experiments of approximately 10 s duration, we perform motion tracking and record the tactile sensors' measurements to compute the best-matching decompositions f

(n) k , f (i) k
following the SOCP of Section 3.5. The next step is to learn the variations of internal forces f (i) k with motion and grasping features. We select the learning parameters as those that directly impact the force distributions through the Newton-Euler equations. Contact forces vary with the object's mass and acceleration, or more accurately including the contribution of gravity F c = m • (a -g). We can consider this dependence as twofold: on the magnitude of F c itself, and on the relative orientation of F c with the contact normals, as in [START_REF] Gao | Internal forces during object manipulation[END_REF]:

p 1 = F c 2 ( 10
)
p 2,k = n k • u F c , with u F c = F c F c 2 . ( 11
)
Similarly, we consider the case of rotational kinematics through the magnitude of the contact torque τ c of Eq. ( 1) and the individual torques each finger is able to generate:

p 3 = τ c 2 (12) p 4,k = --→ CP k × n k • u τ c , with u τ c = τ c τ c 2 . (13)
Finally, we learn internal forces as a function of kinematics and grasp parameters (p 1 , (p 2,k ) k∈F , p 3 , (p 4,k ) k∈F ) using two ANNs: a first network, ANN L 2 , estimates the amount of internal forces applied, quantified as the overall L 2 norm, while a second network, ANN decomp , jointly estimates the relative participation of each finger in the grasp's intensity. The outputs of ANN decomp are percentages constructed as the individual forces normalized with the overall L 2 norm. Note that, as that similar motions can stem from different force distributions, using a single ANN would mean linking similar inputs to highly varying individual forces. Yet, we observed that different grasp intensities still tend to be similarly shared among fingers, hence two ANNs to account for natural intensity variance but consistent decompositions. In order to avoid samples where FSV or tactile sensors are not reliable, we only use as learning data those where the resulting net forces are within a specified threshold from each other. In our experiments, setting this threshold to 1.5N yields a final dataset of 8200 samples, which we partition into training and validation datasets to construct and assess different ANN configurations by cross-validation. Performing numerical resolution with the neuralnet package for statistical analysis software R [START_REF] Fritsch | Package neuralnet, version 1.32[END_REF][START_REF] Core | R: A language and environment for statistical computing[END_REF], we choose ANN L 2 and ANN decomp with logistic neurons trained with resilient backpropagation and two hidden layers, with respectively 6 and 8 neurons in the first hidden layer, and 7 and 13 neurons in the second.

Experiments

In order to assess our approach, we perform manipulation experiments on a rectangular cuboid of dimensions 171mm × 111mm × 60mm. The simplified shape of this ground-truth object is chosen to meet sensing instrumentation constraints and offer several grasping possibilities. We instrument the box with two types of sensors. The first is an Xsens MTi-300 attitude and heading reference system (AHRS) motion sensor measuring reference rotational velocities and translational accelerations. Its purpose is to validate the numerical differentiation of tracking data by algebraic filtering, see Section 3.2. The second consists of five Honeywell FSG020WNPB piezoresistive one-axis force sensors that can be positioned at different predefined grasp spots on the box. We finally evaluate the contact forces estimated from the SOCP in Section 3.4 with the force sensors' measurements in terms of: (i) normal forces per finger, (ii) resulting net force, and (iii) sum of squares. We summarize the validation protocol in Fig. 3.

Kinematics from vision vs AHRS

We assess the validity of our approach by executing motions emphasizing each of the three coordinates of both translation accelerations and rotational velocities, and comparing the kinematics estimated from vision to measurements from the Xsens MTi-300 AHRS. Statistical analysis of the estimation errors shows that algebraic numerical differentiation is well suited for kinematics estimation (see Figure 3 Table 1). Though on translational acceleration the error is slightly higher than with Gaussian filtering, its variance is also considerably lower. Its performance on rotational kinematics is also the best of all three tested approaches. We summarize the results of the six-axis experiments in Fig. 4.

Nominal forces from vision-based kinematics

We now validate our vision-based force estimation framework using Honeywell FSG020WNPB sensors placed at pre-specified positions over the instrumented box. As a first validation step, contact points obtained from vision were compared to the expected contact points based on the sensors' locations and resulted in estimation errors of mean -1.55mm and standard deviation 6.13mm. Furthermore, we assessed the sensitivity of FSV to these uncertainties by comparing the force distributions obtained using either the contact points from vision or the tactile sensors' positions. We found that FSV is relatively robust to such estimation errors, resulting in force uncertainties of mean 0.216N and standard deviation 1.548N . Therefore, we rely solely on vision-based kinematics and contact points for the rest of this work. When performing experiments, we also observed that the force applied by the pinky finger was consistently below the sensitivity threshold of our force sensors, hence we present our results on four-finger experiments. Still, as the force distribution problem introduced in Section 3 becomes statically indeterminate from three fingers, using four fingers preserve the grasps' complexity and does not impact the generality of our results. We represent the force sensors' measurements along with FSV's outputs in Fig. 5. As mentioned in Section 3.4, the comparison of the normal components from vision and from tactile sensors shows that the latter's measurements are overall greater. This illustrates the fact that humans seize objects harder than the required force closure, in contrast with the L 2 -optimal grasp estimated from vision, which is visible in the sum of squares plot. Still, the resulting net forces are matching well, which demonstrates that FSV can successfully capture the object's motion characteristics and compute a reliable force distribution explaining the observed kinematics.

Reconstructing full contact force distributions

By recording new manipulation experiments, we extract as in Section 3.6 the kinematics and grasping parameters (p 1 , (p 2,k ) k∈F , p 3 , (p 4,k ) k∈F ) over time and use the trained ANNs to predict the internal forces the human operator most likely applied over the experiment, ( f (i) k ) k∈F . We finally construct the final contact force distributions using the variant of the SOCP described in Section 3.5 that features an objective function which aims at matching not the full contact forces but only their internal components: We illustrate the final estimation process in Fig. 2(b). By feeding the internal ANNs' predictions into the SOCP, we ensure that the final internal forces

C reconst (x) = k∈F F (n) k 2 2 + f (i) k -f (i) k 2 . (14
F (i) k k∈F
are not only consistent with natural grasping patterns but also physically correct and do not impact the object's observed kinematics through the resulting net force, as shown in Fig. 6.

Robustness analysis

We investigate the robustness of our approach to features that do not appear in the training dataset. To this end, we train another instance of the ANNs described in Section 3.6, not over the entire dataset but on a partial subset relative to a single operator, on a single grasp pose, and a single mass configuration. We then evaluate the resulting ANNs on datasets obtained with another user, another grasp, and/or a 10% mass increase. We report the relative errors with respect to the tactile sensors' measurements in Table 2 along with reference results from fully-trained ANNs.

First, it appears that ANNs trained over a single operator may be generalized to other users with no significant performance decrease, which suggests that humans tend to apply internal forces following similar patterns. Second, reasonable changes in mass do not seem to significantly impact the estimation accuracy either. This is allowed by the fact that in our problem formulation, mass is not a training variable by itself but is implicitly taken into account through the product F c = m • (a -g). Under this formalism, manipulating a heavy object with a given kinematics is analogous to manipulating a lighter object with a higher acceleration. Therefore, the ANNs may accomodate mass changes provided that they were trained over a sufficient variety of kinematics. In the end, ANNs seem most sensitive to grasp pose changes. This may be explained by the fact that placing fingers differently may substantially change their synergies. Still, the performance decrease remains reasonable while force distributions, by construction, still explain the observed motion. Eventually, the main sensitivity to grasp poses is comforted by the fact that also changing user and mass does not decrease the estimation accuracy further.

Summary and Discussion

Our work establishes that a single RGB-D camera can be used to capture interaction forces occurring in rigid object manipulation by a human hand without the need for visual markers or dedicated force sensing devices. Force sensing from vision is a novel and important contribution since it circumvents the intrusive instrumentation of object, environment and hands. Its exploitation can expand to the robotics field for daily on-line human activities monitoring, serving various purposes such as imitation learning.

Our method is validated with several experiments based on ground truth data and is able to estimate fairly accurately the force distributions applied during actual manipulation experiments. Although we confirmed that tracking noise is well mitigated by algebraic filtering, which produces truthful pose derivative estimates, guessing the hand-object contact points under strong occlusions remains a challenging, open problem in computer vision. We achieved this by using a state-of-the-art model-based tracking method under the somewhat practical assumption that occluded fingers remain at their last observed position until they are visible again. While this assumption is fairly valid in numerous interesting cases, it is not true when considering tasks such as dexterous manipulation with finger repositioning or sliding. Still, this limitation does not call into question the force estimation framework per se, and could be alleviated by extending the markerless tracking method to multi-camera inputs, which would remain non-intrusive and keep an edge over tactile sensors regarding usability and cost.

With respect to computational performance, SOCP and internal force predictions are performed in real-time, and only hand-object tracking is computationally expensive. Given the recent developments on GPGPU implementations of hand-object tracking [START_REF] Kyriazis | Scalable 3d tracking of multiple interacting objects[END_REF], our framework could be employed in real-time applications. This, combined with our reliance on a single camera, makes FSV suitable for daily observation and learning. Still, our approach is generic enough to accommodate any advance to the topic of 3D hand tracking and could be seamlessly extended to other methods, for instance when non real-time performance and a heavier setup are possible. Conversely, our framework could also be used as an implicit force model for physicsbased tracking and motion editing, as human-like forces could augment the pose search with biomechanical considerations such as muscle fatigue or energy expenditure.

Towards estimating contact forces from vision, we tackled the issue of static indeterminacy by applying machine learning techniques to internal forces. Rather than predicting new force distributions based on past observations, an alternative approach would be to formulate the evolution of the full contact forces following various objects and grasp taxonomies as an inverse optimal control problem. If invariants are found, they could be used to refine the cost function, which could result in more reliable contact forces than the nominal distributions computed by minimization of the grasp's L 2 -norm. Extending the ground truth force measurement setup with embedded three-axis or force-torque miniature sensors would also benefit both learning and optimal control approaches. Further work could also address the case of surface contact models in place of point contacts (as the fingertip is deforming), namely for dexterous manipulations, or make use of synergy properties of the hand for bimanual tasks. Finally, combining our approach with visual SLAM or automated camera calibration methods would allow it to be deployed in unknown, varying environments, e.g. on mobile robots.

Figure 2 :

 2 Figure 2: (a) Measurements from tactile sensors are used to estimate nominal and internal force decompositions from vision. (b) Full contact forces are reconstructed by feeding ANN predictions into a third SOCP variant.

Figure 4 :

 4 Figure 4: Comparison between vision-based kinematics and AHRS-embedded accelerometer and gyroscope.

Figure 5 :

 5 Figure 5: Contact forces from vision based on L 2 criterion are individually lower than tactile sensors' measurements but result in the same net force.

Figure 6 :

 6 Figure 6: Artificial neural networks used in conjunction with cone programming successfully predict force distributions that both explain the observed motion and follow natural human force distribution patterns.

  and control their magnitude

	Vision		Tactile Sensors	Vision
	Tracking		Measurements
	Kinematics,		Finger forces	Kinematics, contact points	ANNL2 , ANNdecomp
	contact points		( fk)k∈F
					Internal forces
					( f
	SOCP
	min x	Cdist(x)
	f k , f (n) k (i)	k∈F
		(a)

Table 1 :

 1 Kinematics estimation errors for central finite difference, Gaussian filtering, and algebraic filtering.

			: Validation
			Central Gaussian Algebraic
	Trans. acc. Avg. -0.029 -0.022	-0.024
	[m • s -2 ]	St.d.	1.686	1.627	0.904
	Rot. vel.	Avg.	0.084	0.070	0.052
	[rad • s -1 ] St.d.	1.559	1.294	1.241

Table 2 :

 2 Relative force estimation errors based on the exhaustivity of the training dataset. and × indicate features that respectively appear or not in the partial training dataset.

	)

Acknowledgements This work was partially supported by the FP7 EU RoboHow.Cog project.