Thomas Ricatte
email: thomas@ricatte.fr

Rémi Gilleron
email: remi.gilleron@univ-lille3.fr

Marc Tommasi
email: marc.tommasi@univ-lille3.fr

Learning with Hypernode Graphs Hypernode Graphs for Learning from Binary Relations between Groups in Networks *

Keywords: Hypergraphs, Graph Laplacians, Graph Kernels, Spectral Learning, Semisupervised Learning, Skill Rating Algorithms

The aim of this paper is to propose methods for learning from interactions between groups in networks. We propose a proper extension of graphs, called hypernode graphs as a formal tool able to model group interactions. A hypernode graph is a collection of weighted relations between two disjoint groups of nodes. Weights quantify the individual participation of nodes to a given relation. We define Laplacians and kernels for hypernode graphs and prove that they strictly generalize over graph kernels and hypergraph kernels. We then proceed to prove that hypernode graphs correspond to signed graphs such that the matrix D -W is positive semi-definite. As a consequence, homophilic relations between groups may lead to non homophilic relations between individuals. We also define the notion of connected hypernode graphs and a resistance distance for connected hypernode graphs. Then, we propose spectral learning algorithms on hypernode graphs allowing to infer node ratings or node labelings. As a proof of concept, we model multiple players games with hypernode graphs and we define skill rating algorithms competitive with specialized algorithms.

Introduction

Networks are commonly modeled by graphs where nodes correspond to individual objects and edges correspond to binary relationships between individuals. This is for instance the case for social networks with friendships between users, or computer networks with connections between machines. But in many situations, only interactions between groups are observed without any information on individuals besides group membership. As an example, in team sports and online games only outcomes of games between teams of multiple players are filed. Also, group interactions are observed when dealing with clusters in computer networks or communities in social networks. However in such contexts, the simultaneous evaluation of functions on groups and on individuals play an important role for applications. This is exactly the case for functions that measure CPU or cluster loads, individual players or team skills, users or communities ratings. The two evaluation levels of functions on items and functions on groups are however related by assuming an additive model. For instance, in multiple player games, the team skill is assumed to be the weighted sum of player skills as in [START_REF] Herbrich | TrueSkill TM : A Bayesian Skill Rating System[END_REF]. A task that can be defined under these assumptions is to evaluate the function values of all groups and individuals when only a few values are observed. In the case of multiple player games, skills must be evaluated from a limited number of game outcomes, thus allowing to predict game outcomes for new games or to do matchmaking with well-balanced games.

In this paper, we propose a new model called hypernode graphs as a formal tool able to model group interactions. Hypernode graphs allow to consider node valuation problems as well as group valuation problems by assuming an additive model. We define a hypernode graph to be a set of hyperedges. A hyperedge is a pair of disjoint hypernodes, where a hypernode is a set of nodes. Every node of a hypernode is given a non negative real-valued weight and weights in a hyperedge satisfy an equilibrium condition. Roughly speaking, a hypernode models a group, a hyperedge models a relation between two groups, and individual weights correspond to the contribution of each individual to the relation between the two groups. An example of hypernode graph is shown in Figure 1. Hypernode graphs generalize over graphs because there is a one to one correspondence between undirected graphs and hypernode graphs in which all hypernodes are singleton sets.

Then, a learning task is to infer node function values over a hypernode graph. We consider semi-supervised setting where function values are given for only a few nodes or a few hypernodes. For semi-supervised learning on hypernode graphs, we assume that connected hypernodes tend to have similar values. For real-valued node functions over discrete structures, the gradient is a discrete derivative that allows to measure such a similarity. Therefore, we define a gradient for hypernode graphs. If we denote by G the gradient of a hypernode graph h, we define the Laplacian ∆ of h to be ∆ = G T G, and we define the kernel ∆ † to be the Moore-Penrose pseudo-inverse of the Laplacian ∆. Then, the quantity f T ∆f measures the smoothness of a real-valued node function f on the hypernode graph h, i.e. f T ∆f is small when connected hypernodes have close values. As the Laplacian is positive semi-definite, (graph) spectral learning algorithms (see [START_REF] Von | A tutorial on spectral clustering[END_REF][START_REF] Zhou | Learning from labeled and unlabeled data on a directed graph[END_REF][START_REF] Zhu | Semi-supervised learning using gaussian fields and harmonic functions[END_REF] can be applied for learning in hypernode graphs.

As a proof-of-concept, we consider the skill rating problem for multiple players games, i.e. how to infer individual skills of players from game outcomes in multiple players games. We consider competitive games between two teams where each team is composed of an arbitrary number of players. We model a set of games by a hypernode graph where hypernodes are teams of players. A player skill rating function is thus a function on nodes of the hypernode graph. Then, in order to model game outcomes, we add outcome nodes. An outcome node is added to every losing team in a hyperedge for a game. Such an outcome node will be given a positive value in order to model the relation between team skills given by the game outcome. Thus, the inductive bias used to learn skill ratings is to choose a smooth function on the hypernode graph such that outcome nodes values are fixed. Therefore, the skill rating problem is equivalent to a semi-supervised learning problem over a hypernode graph. The SSL problem can be solved by minimizing Ω(s) = s T ∆s where ∆ is the Laplacian of the hypernode graph and where function values for game outcome nodes are given. The inferred skill rating function allows to compute team ratings and to predict game outcomes for new games. We consider real datasets of multiple players games (double tennis and online games) and apply our method in a batch setting. Experimental results show that we obtain very competitive results compared to specialized algorithms such as Elo duelling (see [START_REF] Emrick | The Rating of Chess Players, Past and Present[END_REF] and TrueSkill (see [START_REF] Herbrich | TrueSkill TM : A Bayesian Skill Rating System[END_REF].

One question raised by the introduction of our new framework is whether or not problems on hypernode graphs can be solved using graphs. Before answering this question, let us consider the case of hypergraphs, a popular generalization of graphs where a hyperedge is a set of nodes. Hypergraphs have been studied from a machine learning perspective and have been used in many applications (see [START_REF] Klamt | Hypergraphs and Cellular Networks[END_REF][START_REF] Zhang | The automatic construction of a view-independent relational model for 3-D object recognition[END_REF][START_REF] Cai | End-to-end coreference resolution via hypergraph partitioning[END_REF]. It has been proved in [START_REF] Agarwal | Higher Order Learning with Graphs[END_REF] that the hypergraph Laplacians introduced so far do not provide a strict gain of expressiveness with respect to graph Laplacians since we can reduce any learning problem based on a hypergraph Laplacian to a learning problem based on a graph Laplacian using an adequate graph construction. This is not true for hypernode graphs. Indeed, we show that no graph construction allows to define a set of smooth functions over the graph which coincides with the set of smooth functions over a given hypernode graph. Thus, we state that hypernode graphs Laplacians strictly generalize over graph Laplacians.

It can be noted that the class of hypernode graph Laplacians strictly contains the class of graph Laplacians and also contains the class of graph kernels (a graph kernel is the Moore-Penrose pseudo-inverse of a graph Laplacian). We can also observe that a convex linear combination of graph kernels (used in multiple graph kernel learning as in [START_REF] Argyriou | Combining Graph Laplacians for Semi-Supervised Learning[END_REF]) is a hypernode graph kernel (i.e., the Moore-Penrose pseudo-inverse of a hypernode graph Laplacian). We also show that, for every hypernode graph, we can construct a signed graph with adjacency matrix W such that D -W is the Laplacian of the hypernode graph. The fact that W potentially contains negative entries should be related to the strict gain of expressiveness provided by our new framework. Let us recall that, for arbitrary signed graphs, the matrix D -W can be indefinite. Thus, many works have studied the notion of Laplacian for signed graphs. [START_REF] Herbster | Exploiting cluster-structure to predict the labeling of a graph[END_REF] has considered the class of symmetric diagonally dominant matrices with non negative diagonal entries. While others (see [START_REF] Yao | Bounds for the least Laplacian eigenvalue of a signed graph[END_REF][START_REF] Andrew B Goldberg | Dissimilarity in graph-based semi-supervised classification[END_REF][START_REF] Kunegis | Spectral analysis of signed graphs for clustering, prediction and visualization[END_REF] have used a modified definition of the Laplacian to obtain a positive semi-definite Laplacian. Last, [START_REF] Koren | ACE: a fast multiscale eigenvectors computation for drawing huge graphs[END_REF] have considered the class of signed graphs such that D -W is positive semi-definite for the problem of drawing graphs. A signed graph in this class can therefore be represented by a hypernode graph with the same Laplacian.

We finally study the definition of a distance between nodes in hypernode graphs. Indeed, it has been shown that, for connected graphs, a resistance distance can be defined from the graph kernel and that it can be expressed in term of the commute time distance in the graph (see [START_REF] Klein | Resistance distance[END_REF][START_REF] Ashok | The electrical resistance of a graph captures its commute and cover times[END_REF]. We define a pseudo-distance between nodes of a hypernode graph using the hypernode graph kernel. We show that it is a distance for a well chosen notion of connected hypernode graphs. But, it should be noted that we leave open the question of finding an algorithmic definition of connected components in hypernode graphs. Also, we introduce a diffusion operator for hypernode graphs and we show that the distance can be expressed from the differences of potentials. As a special case, we get the result on the commute time distance for graphs. But, when the hypernode graph is not a graph, the question of finding an interpretation of the distance between two nodes with random walks is left open because negative terms are involved.

Hypernode Graphs

The following definition is our contribution to the modeling of binary relationships between sets of entities.

Definition 1 A hypernode graph h = (V, H) is a set of nodes V and a set of hyperedges H. Each hyperedge h ∈ H is an unordered pair {s h , t h } of two non empty and disjoint hypernodes (a hypernode is a subset of V ). Each hyperedge h ∈ H has a weight function w h mapping every node i in s h ∪ t h to a positive real number w h (i) (for i / ∈ s h ∪ t h , we define w h (i) = 0). Each weight function w h of h = {s h , t h } must satisfy the Equilibrium Condition defined by

i∈t h w h (i) = i∈s h w h (i) .
In the following, we denote by n the number of nodes, we denote by p the number of hyperedges, and we assume an ordering of the set of nodes and of the set of hyperedges. We say that a node i belongs to a hyperedge h, that we denote by i ∈ h, if w h (i) = 0. We define the degree of a node i by deg(i) = h w h (i). The degree of a node is positive when it belongs to at least one hyperedge. We define the diagonal degree matrix by D = diag(deg(1), . . . , deg(n)) and the volume of the hypernode graph by Vol(h) = i∈N deg(i) = Tr(D).

An example of hypernode graph is shown in Figure 1. The blue hyperedge h blue links the sets {1, 2, 3} and {4, 5}. The weights of h blue satisfy the Equilibrium condition:

1/2 + 1/2 + √ 2 = √ 2 + √ 2.
The green hyperedge h green links the sets {4, 5} and {1, 6} and the Equilibrium condition is obviously satisfied. The red hyperedge h red links the two singleton sets {3} and {6} and the weights of h red are equal, thus the hyperedge h red can be viewed as an edge with edge weight 1/3.

As noted in the previous example, when a hyperedge h is an unordered pair {{i}, {j}} of two nodes i, j, the Equilibrium Condition states that the weights w h (i) and w h (j) are equal. Therefore, every hypernode graph such that all hyperedges are unordered pairs of singleton nodes can be viewed as a graph with adjacency matrix W defined by W i,j = W j,i = w h (i) = w h (j) for every hyperedge {{i}, {j}}, and 0 otherwise. Conversely, a graph can be viewed as hypernode graph where every hypernode is a singleton set.

We are mainly interested in evaluating real-valued functions on nodes and on hypernodes. We assume a weighted linear model such that any real-valued node function f can be extended to a hypernode u of a hyperedge h by f (u) = i∈u f (i) w h (i). Moreover, we will consider an homophilic assumption stating that connected hypernodes tend to have similar values. In order to implement this assumption, we introduce a gradient for hypernode graphs by Definition 2 Let h = (V, H) be a hypernode graph and f be a real-valued node function, the (hypernode graph) unnormalized gradient of h is a linear application, denoted by grad, that maps every real-valued node function f into a real-valued hyperedge function grad(f ) defined, for every

h = {s h , t h } in H, by grad(f )(h) = f (t h ) -f (s h ) = i∈t h f (i) w h (i) - i∈s h f (i) w h (i) , ( 1
)
where an arbitrary orientation of the hyperedges has been chosen. We denote by G ∈ R p×n the matrix of grad.

Because of the Equilibrium Condition, the gradient of every constant node function is the zero-valued hyperedge function. The square n × n real valued matrix ∆ = G T G is defined to be the unnormalized Laplacian of the hypernode graph h. When the hypernode graph is a graph (all hypernodes are singleton sets), then the hypernode graph Laplacian is equal to the unnormalized graph Laplacian. The Laplacian ∆ does not depend on the arbitrary orientation of the hyperedges used for defining the gradient. By definition, the Laplacian ∆ is positive semi-definite. We define the smoothness of a real-valued node function f over a hypernode graph h to be Ω(f ) = f T ∆f . Last, we define the hypernode graph kernel of a hypernode graph h to be the Moore-Penrose pseudo-inverse ∆ † of the hypernode graph Laplacian ∆.

Because a hypernode graph Laplacian is positive semi-definite, we can leverage the spectral learning algorithms defined in Von Luxburg (2007), [START_REF] Zhou | Learning from labeled and unlabeled data on a directed graph[END_REF], and [START_REF] Zhu | Semi-supervised learning using gaussian fields and harmonic functions[END_REF] from graphs to hypernode graphs. In the next section we show how to use spectral learning algorithms over hypernode graphs for the skill rating problem in multiple players games.

Skill Rating for Multiple player Games

We consider competitive games between two teams where each team is composed of an arbitrary number of players. A first objective is to compute the skill ratings of individual players from a batch of games with their outcomes. A second objective is to predict a game outcome for a new game. In this section, we show how to model a batch of games by a hypernode graph and we apply learning algorithms over the hypernode graph to infer player skills that allow to predict game outcomes of new games.

Multiplayer Games

Let us consider a set of players P = {1, . . . , n} and a set of games Γ together with their respective game outcome. Each game is between two teams of an arbitrary number of players. Let us also consider that a player i contributes to a game γ with a positive real value c(i). Let us assume that each player has a skill s(i), we suppose an additive model as in [START_REF] Herbrich | TrueSkill TM : A Bayesian Skill Rating System[END_REF] stating that the skill of a team is the weighted sum of the skills of the players in the team. More formally, given two teams of players A = {a 1 , a 2 , . . . , a } and

B = {b 1 , b 2 , . . . , b k } playing game γ, then A is predicted to be the winner if i=1 c(a i )s(a i ) > k i=1 c(b i )s(b i ) .
(2)

Equivalently, one can rewrite this inequality by introducing a positive real number o on the right hand side such that

i=1 c(a i )s(a i ) = o + k i=1 c(b i )s(b i ) , ( 3
)
where the real number o quantifies the game outcome. In the case of a draw, there is an equality between the team skills, i.e. the game outcome o is set to 0. Given a set of games with their respective outcome, the goal is to infer a skill rating function s ∈ R n that respects equations 3 as much as possible. We define the cost of a game γ with outcome o for a skill function s by

C γ (s) = i=1 c(a i )s(a i ) -(o + k i=1 c(b i )s(b i )) 2 .
Consequently, given a set of games Γ and the corresponding game outcomes, the goal is to find a skill rating function s * that minimizes the sum of the different costs, i.e. search for

s * = arg min s γ∈Γ C γ (s) . (4)

Modeling Games with Hypernode Graphs

We introduce the general construction by considering an example. Let us consider a game γ between two teams A = {1, 2, 3} and B = {4, 5}. Let us also assume that all the players contribute to the game with weight 1. Such a game can be modeled by a hyperedge between sets of nodes {1, 2, 3} and {4, 5} with weights equal to 1. Now, let us suppose that B wins the game, then the skill rating function s must satisfy Equation (3), that is s(1) + s( 2 then the skills must satisfy s(1) + s(2) + s(3) + 2 is equal to s(4) + s( 5). And, if we constraint a skill function on the hyperedge to satisfy s(O) = 2 and s(Z) = 0, then the homophilic assumption on the hyperedge states that s(1) + s(2) + s(3) + 2 must be close from s(4) + s( 5) which expresses the expected relations between player skills given the game outcome.

In the general case, let us consider a set of players P = {1, . . . , n} and a set of games Γ. Each game γ is between two teams A and B of multiple players, the winning team is known as well as the game outcome o. Let us also consider that a player i contributes to a game γ with a non negative weight c(i). We can define, for every game γ a hyperedge h as follows 1. The players of A define one of the two hypernodes of h, the weight of a player node i is defined to be c(i) 2 ;

2. do the same construction for the second team B;

3. add a outcome node O to the set of player nodes corresponding to the losing team, its weight is set to 1;

4. add a lazy node Z to the set of player nodes corresponding to the winning team, its weight is chosen in order to ensure the Equilibrium condition for the hyperedge h.

We define the hypernode graph h = (V, H) as the set of all hyperedges h for all games γ in Γ. Now, a skill rating function is the restriction to player nodes of a real-valued node function over the hypernode graph. In order to model the game outcomes, we fix the skill rating function values over the additional nodes. A skill function s over h modeling Γ must satisfy: for every lazy node Z, the function value s(Z) must be 0, and, for every outcome node O, the function value s(Z) must be the outcome value o. Let ∆ be the unnormalized Laplacian of h, and let s be a real-valued node function on h. Then, it is easy to show that the skill rating problem (4) is equivalent to find a real valued node function s over h solution of minimize s s T ∆s subject to for every outcome node O, s(O) = o; for every lazy node Z , s(Z) = 0

(5)

Regularizing the hypernode graph

When the number of games is small, the number of hyperedges is small and some node players belong to only one hyperedge. Thus, player skills can be defined independently while satisfying the constraints and it will be irrelevant to compare them. In order to solve this issue, we introduce in Equation ( 5) a regularization term based on the standard deviation σ(s p ) of players skills where s p = (s(1), . . . , s(n)). This allows us to reduce the spread of the player skills and leads us to find a real valued node function s over h solution of where µ is a regularization parameter. That is, we control the spread of s p avoiding extreme values for players involved in a small number of games.

We show in Appendix A that the regularized optimization problem 6 can be rewritten as find a real valued node function s over h µ solution of where ∆ µ is the Laplacian of the hypernode graph h µ obtained from the hypernode graph h by: add a regularizer node R; for every player node, add an hyperedge between the player node and the regularizer node R and set the (node) weights to µ/n. The added hyperedges can be viewed as edges with weight µ/n.

Note that, assuming unitary contributions and a Gaussian distribution for skill ratings, we can prove that the value of µ/n should have the same order of magnitude than the average number of games played by a player. Indeed, this will allow us to obtain similar expected values for the two terms s T ∆s and µσ(s p ) 2 . We follow this guideline in the experiments presented below.

Inferring Skill Ratings and Predicting Game Outcomes

The problem (7) can be viewed as a semi-supervised learning problem on the hypernode graph h µ because the question is to predict player node values (skill ratings) from known values for lazy nodes and outcome nodes. The Laplacian ∆ µ is positive semi-definite because it is a hypernode graph Laplacian. Therefore, we can use the semi-supervised learning algorithm presented in [START_REF] Zhu | Semi-supervised learning using gaussian fields and harmonic functions[END_REF]. This algorithm was originally designed for graphs and solves exactly the problem (7) by putting hard constraints on the outcome nodes and on the lazy nodes. It should be noted that the algorithm is parameter free. We denote this learning method by H-ZGL.

Another approach to solve the semi-supervised learning problem is to use a regression algorithm. For this, we consider the hypernode graph kernel ∆ † µ (defined as the Moore-Penrose pseudo-inverse of the Laplacian ∆ µ ), we train a regression support vector machine using outcome nodes values and lazy nodes values, and we predict node player values. The main parameter is the soft margin parameter. We denote this method by H-SVR.

Using one of these two methods, we can infer players skills, then deduce teams skills, and finally predict game outcomes for new games. For this, we suppose that we are given a training set of games Γ l with known outcomes together with a test set of games Γ u for which game outcomes are hidden. The goal is to predict game outcomes for the test set Γ u . We use the following Algorithm 1 Predicting game outcomes Input: Training set of games Γ l , set of testing games Γ u 1: Build, as described in Sections 3.2 and 3.3, the regularized hypernode graph h µ from Γ l 2: Compute an optimal skill rating s * using H-ZGL or H-SVR.

3: Compute the mean skill s among players in Γ l 4: for each game in Γ u do 5:

Assign skills given by s * to players involved in Γ l , and s otherwise 6:

Evaluate the inequality (2) and predict the winner 7: end for

Experiments

In this section, we report experimental results for the inference of player skills and the prediction of game outcomes for different datasets in the batch setting described above. Note that other works have considered the online setting as in [START_REF] Herbrich | TrueSkill TM : A Bayesian Skill Rating System[END_REF] or Elo (1978).

Tennis Doubles

We consider a dataset of tennis doubles collected between January 2009 and September 2011 from ATP tournaments (World Tour, Challengers and Futures). Tennis doubles are played by two teams of two players. Each game has a winner (no draw is allowed). A game is played in two or three winning sets. The final score corresponds to the number of sets won by each team during the game. The dataset Γ consists in 10028 games with 1834 players.

Along the experimental process, given a proportion ρ, we draw randomly a training set Γ l of size ρ% of the number of games in Γ. The remaining games define the test set Γ u . We present in Figure 3 several statistics related to this process on the Tennis dataset. First, it can be noticed that many players have played a small number of games. Second, when the number of games in the training set is small, half of the players are involved in games in the test set. Therefore, the skill rating problem and the game outcome prediction problem become far more difficult to solve when few games are used for learning.

The experimental setting is defined as follows. Let ρ be a proportion varying from 10% to 90% by 10%. For every value of ρ, we repeat ten times: draw a random training set Γ l of size ρ% of the number of games in Γ; let the remaining games define the test set Γ u ; infer player skills and predict game outcomes following Algorithm 1. For the definition of the hypergraph, we set all node player weights to 1 because we do not have information on player contributions. The node outcome values o are set to be the difference between the number of sets won by the two teams, thus the possible outcome node values are 1, 2 or 3. Therefore, there are 3 outcome nodes and one lazy node shared by all hyperedges. The resulting hypernode graph has at most 1839 nodes: at most 1834 player nodes, 1 lazy node, 3 outcome nodes, and 1 regularizer node. In order to complete the definition of the We report in Figure 4 prediction results using: algorithm 1 and H-ZGL, algorithm 1 and H-SVR, Elo Duelling, and Trueskill 1 . It can be noted that Elo Duelling performs poorly but Elo was designed for one against one games. It should be noted that we have done similar experiments for tennis singles. The results are not reported here but they show that Elo Duelling and True skill obtain similar results but are outperformed by H-ZGL and H-SVR. For the dataset of tennis doubles, it must be noted that H-ZGL and H-SVR outperform Trueskill for a small number of training games. Also, H-ZGL outperforms Trueskill whatever is the number of training games. On this set of experiments H-ZGL outperforms H-SVR but it should be noted that H-ZGL is parameter free and that we use H-SVR with the default soft margin parameter value.

Xbox Title Halo2

The Halo2 dataset was generated by Bungie Studio during the beta testing of the XBox title Halo2. It has been notably used in [START_REF] Herbrich | TrueSkill TM : A Bayesian Skill Rating System[END_REF] to evaluate the performance of the Trueskill algorithm. We consider the Small Teams dataset with 4992 players and 27536 games opposing up to 12 players in two teams which can have a different size. Each game can result in a draw or a win of one of the two teams. The proportion of draws is 22.8%. As reported in [START_REF] Herbrich | TrueSkill TM : A Bayesian Skill Rating System[END_REF], the prediction of draws is challenging and it should be noted that Trueskill and our algorithm fail to outperform a random guess for the prediction of draw.

We consider the same experimental process. For the construction of the hypernode graph in Algorithm 1, we fix all players contributions in games to 1 and we again set the parameter 1. TrueSkill and Elo implementations are from [START_REF] Hamilton | PythonSkills: Implementation of the TrueSkill, Glicko and Elo Ranking Algorithms[END_REF]. Results were double-checked using Lee (2013b) and Lee (2013a) value µ/n to 16. In the optimization problem 7, the game outcomes o are set to 1 when the game has a winner and 0 otherwise because game scores vary depending on the type of game. Therefore in the hypernode graph there are one lazy node and two outcome nodes (for 0 and 1) shared by all hyperedges. We again compare the skill rating algorithms H-ZGL, H-SVR, Elo Duelling and Trueskill. The number of prediction errors over game outcomes is computed assuming that a draw can be regarded as half a win, half a loss as in [START_REF] Lasek | The predictive power of ranking systems in association football[END_REF]. We present the experimental results in Figure 5. For a proportion of 10% of games in the training set, H-ZGL, H-SVR and Trueskill give similar results while with larger training sets, our hypernode graph learning algorithms outperform Trueskill. Contrary to the previous experiment, H-SVR (with default parameter value) outperforms H-ZGL.

Hypernode Graphs, Graphs and Signed Graphs

In this section, we study the class of hypernode graph Laplacians. In particular, we show that any hypernode graph Laplacian can be seen as the Laplacian matrix of a signed graph (i.e., a graph that potentially contains negative edge weights).

The Class of Hypernode Graphs

Create a new hyperedge h

= {s h , t h } with s h = t h = ∅ 5:
for each node i ∈ V do

6:

Define w h (i) = r(i) 2

7:

if r(i) < 0 then Add i to s h else if r(i) > 0 Add i to t h end if 8:

end for

9:

Add h to H 10: end for 11: return The hypernode graph h = (V, H) From a square root decomposition G T G of M , for each line of G, we define a hyperedge h = {s h , t h } where nodes in s h (respectively in t h ) have positive (respectively negative) values in the line, and node weights are chosen to be squares of values in the line. Then, the Equilibrium condition 1 is satisfied because 1 ∈ Null(M ). And it is easy to check that G is a gradient matrix of h and, consequently, that M = G T G is the Laplacian of h.

As a consequence of Proposition 3, the class of hypernode graph Laplacians

• is closed under convex linear combination,

• is closed under pseudo-inverse,

• and strictly contains the class of graph Laplacians and the class of graph kernels (a graph kernel is the Moore-Penrose pseudo-inverse of a graph Laplacian).

Linear combinations of graph kernels have been used for semi-supervised learning (see for instance [START_REF] Argyriou | Combining Graph Laplacians for Semi-Supervised Learning[END_REF]. Another consequence of Proposition 3 is that a convex linear combination of graph kernels is a hypernode graph kernel (also a hypernode graph Laplacian). To conclude this section, we propose a conjecture Conjecture 4 The class of hypernode graph Laplacians is the smallest class closed by convex linear combinations which contains all graph kernels.

The difficult part is to prove that every hypernode graph Laplacian is a convex linear combination of graph kernels. This is because a graph kernel, which is defined to be the pseudo-inverse of a graph Laplacian, has no simple analytic form.

Hypernode Graph Laplacians Strictly Generalize Graph Laplacians

As said above, the class of hypernode graph Laplacians strictly contains the class of graph Laplacians. But this does not allow us to claim that hypernode graph Laplacians provide a gain of expressiveness over graph Laplacians. Indeed, it could be the case that through a well chosen graph construction, the set of smooth functions defined on hypernode graphs could be made to coincide exactly with the set of smooth functions defined on graphs. The class of hypergraph Laplacians has been studied in that perspective of expressiveness. Indeed, it has been shown in [START_REF] Agarwal | Higher Order Learning with Graphs[END_REF] that all hypergraph Laplacians defined so faramong them the Laplacians ∆ B from [START_REF] Bolla | Spectra, euclidean representations and clusterings of hypergraphs[END_REF], ∆ R from Rodríguez ( 2003) and ∆ ZHS from [START_REF] Zhou | Learning with hypergraphs: Clustering, classification, and embedding[END_REF] -can be defined as (restrictions of) graph Laplacians using a graph construction such as the clique expansion (where each hyperedge is replaced by a clique graph with uniform weights) or the star expansion (where, for every hyperedge, a new node is added and is linked with all the nodes in the hyperedge).

While one can think of similar constructions for the case of hypernode graphs, we prove that there does not exist a (finite) graph expansion of a hypernode graph which defines the same set of smooth functions over the original set of nodes. The proof is based on the very simple hypernode graph h shown in Figure 6. Note that a function is smooth over h if and only if it satisfies

(f (1) + f (2) -f (3) -f (4)) 2 = 0 (C)
We show by contradiction in Appendix B that there does not exist a finite graph g whose node set contains {1, 2, 3, 4} and that satisfies the following conditions: 1. Every smooth function over g satisfies (C) 2. Every function satisfying (C) can be extended into a smooth function over g. The first condition ensures that smooth functions on g are related to smooth functions over the original hypernode graph h. The second condition ensures that smooth functions on the original hypernode graph can be defined from smooth functions over extended finite graphs. The combination of these two conditions ensure that the problem of finding smooth functions on h can be rephrased as a graph problem (finding the smooth functions on g and compute their restrictions to {1, 2, 3, 4}). The fact that no graph g can satisfy both conditions at the same time for a simple hypernode graph h shows that the smoothness conditions defined by the hypernode graphs can not be defined with undirected extended graphs.

G = 1 1 -1 -1 0 0 -1 1 1 2 3 4 0.5 0.5 2 0.5 0.5 2 G = - √ 2 2 - √ 2 2 √ 2 0 - √ 2 2 - √ 2 2 0 √ 2
Figure 7: two L-equivalent hypernode graphs with their respective gradient which are square root of the Laplacian ∆ shown in Figure 8.

L-equivalent Hypernode Graphs and Signed Graphs

Let us consider a symmetric positive semi-definite real-valued matrix M such that 1 ∈ Null(M ), M is the Laplacian of a hypernode graph by Proposition 3. But, as the square root decomposition is not unique, there are several hypernode graphs with the same Laplacian that we called L-equivalent. Examples of L-equivalent hypernode graphs are given in Figure 7. In order to study the L-equivalent relation and to obtain results showing whether a hypernode graph is L-equivalent to a graph, we first show that hypernode graphs are related to signed graphs.

It is known that the Laplacian matrix ∆ of a graph can be written D -W where W is the adjacency matrix of the graph and D is the corresponding degree matrix. Let us consider a hypernode graph h with Laplacian ∆, we define the pairwise weight matrix W of h by W i,j = -∆ i,j if i = j, and 0 otherwise. Note that the pairwise weight matrix coincides with the classic adjacency matrix when the hypernode graph is a graph. Let us define the degree of a node i to be deg(i) = j∈V W i,j and let us denote by D the diagonal matrix of all deg(i). Then, because of the property 1 ∈ Null(M ) in Proposition 3, it is immediate that, for every i, ∆ i,i = j∈V W i,j . As a consequence, we have Proposition 5 Let h = (V, H) be a hypernode graph, let W be the pairwise weight matrix of h, and let D be the diagonal degree matrix of h. Then, the Laplacian of h is ∆ = D -W .

An example is shown in Figure 8. As a consequence of the above proposition, we can leverage the pairwise weight matrix to characterize L-equivalent hypernode graphs by Proposition 6 Two hypernode graphs are L-equivalent if and only if they have the same pairwise weight matrix. Two L-equivalent hypernode graphs have the same degree matrix.

∆ =      1 1 -1 -1 1 1 -1 -1 -1 -1 2 0 -1 -1 0 2      ; W =      0 -1 1 1 -1 0 1 1 1 1 0 0 1 1 0 0      ; D =      1 0 0 0 0 1 0 0 0 0 2 0 0 0 0 2     
.

Figure 8: Laplacian matrix, pairwise weight matrix, and degree matrix for the two Lequivalent hypernode graphs shown in Figure 7. The pairwise weight matrix W contains in general negative weights and can be thus interpreted as the adjacency matrix of a signed graph. Therefore, we define the reduced signed graph of a hypernode graph h to be the signed graph g with adjacency matrix W . An example is shown in Figure 9. Then, we can show that • the reduced signed graph of a graph g (viewed as a hypernode graph) is the graph g,

• a hypernode graph h is L-equivalent to a graph if and only if its reduced signed graph is a graph,

• a signed graph with adjacency matrix W is the reduced signed graph of a hypernode graph if and only if the matrix D -W is positive semi-definite.

The definition of the pairwise weight matrix may seem ad hoc to mimic the definition of the Laplacian in the graph case. But, we can show that the pairwise weight matrix can be defined directly from the hypernode graph using the following formula

∀i = j, W i,j = h∈H P (h, i, j) w h (i) w h (j) , ( 8
)
where P (h, i, j) = 1 if i and j belongs to two different ends of h, P (h, i, j) = -1 if i and j belongs to the same end of h, and 0 otherwise. And the formula can be interpreted as follows. The pairwise weight W i,j is computed as a sum over all the hyperedges, then it can be understood as an aggregation of all the information over hyperedges involving the nodes i and j. For every term in the sum, the quantity w h (i) can be viewed as the cost of entering the hyperedge h at node i, the quantity w h (j) as the cost of exiting the hyperedge h at node j, and P (h, i, j) as a sign depending whether i and j are in the same end or in different ends of the hyperedge.

Resistance Distance and Random Walks in Hypernode Graphs

In this section, we study whether a distance can be defined between nodes of a hypernode graph and how such a distance can be interpreted in the hypernode graph. Throughout the section, we consider a hypernode graph h = (V, H) with Laplacian ∆.

Defining a distance in Hypernode Graphs

Let ∆ † be the Moore-Penrose pseudo-inverse of ∆, let us define d by, for every i, j in V ,

d(i, j) = ∆ † i,i + ∆ † j,j -2∆ † i,j . ( 9
)
Because ∆ † is symmetric positive semi-definite, we have Proposition 7 d defines a pseudo-metric on h, i.e., it is positive, symmetric and satisfies the triangle inequality (for all i, j, k the inequality d(i, j) ≤ k) + d(k, j) holds).

For the pseudo-metric d to be a distance, d should satisfy also the coincidence axiom: d(i, j) = 0 ⇒ i = j. This is not true in general as d(1, 2) = 0 for the hypernode graph in Figure 6. The intuition is that the nodes 1 and 2 can not be distinguished because the smoothness condition is on the sum f (1) + f (2). Nevertheless we can show that Proposition 8 When Null(∆) = Span(1), d defines a metric (or distance) on h.

Let us recall that, in the graph case, Proposition 7 holds and that Proposition 8 holds when the graph is connected. By analogy with the case of graphs, the condition Null(∆) = Span(1) can be viewed as an algebraic definition of a connected hypergraph. So far, we have not found an alternative algorithmic definition of connected components in hypernode graphs.

Finally, let us note that, in the graph case, for connected graphs, both d and d 2 are metrics while, for hypernode graphs, d 2 does not satisfy the triangle inequality even if Null(∆) = Span(1).

Resistance Distance and Diffusion in Hypernode Graphs

Let us suppose that h satisfies Null(∆) = Span(1). Our goal is to study whether d 2 can be written in terms of a diffusion function in the hypernode graph. First, let us consider the Poisson equation ∆f = In that models the diffusion of an input charge In through a system associated with the Laplacian operator ∆. Let us consider a node j called sink node, we consider the input function In j defined by In j (j) = deg(j) -Vol(h) and In j (i) = deg(i) if i = j. We can prove that Lemma 9 The solutions of ∆f = In j are the functions f = µ1 + ∆ † In where µ ∈ R.

The compete proof is presented in Appendix C. It is based on properties of pseudoinverse matrices and on the hypothesis Null(∆) = Span(1). Then, for every pair of nodes (k, ) in V , we define V j (k, ) as the difference of potential between k and , i.e., we define V j (k, ) = f (k) -f ( ) where f is a solution of the Poisson equation ∆f = In j . Lemma 9 allows to show that the definition of V j (k, ) does not depend on the choice of the solution of the Poisson equation and that an equivalent definition is

V j (k, ) = (e k -e ) T ∆ † In j ,
where e i is the unit vector with 1 in component i. We can now relate the distance d 2 and the diffusion potential V j by Proposition 10 For every i, j in V , we have Vol(h)d 2 (i, j) = V j (i, j) + V i (j, i), where V j (i, i) = 0, and for i = j,

V j (i, j) = h|i∈h w h (i) deg(i)   1 + k∈h,k =i P (h, i, k) w h (k) w h (i) V j (k, j)   . ( 10
)
The proof is given in Appendix D.

Resistance Distance and Random Walks in Hypernode Graphs

Proposition 10 provides an expression of d 2 in terms of diffusion potentials. We first show that the proposition generalizes the result of [START_REF] Klein | Resistance distance[END_REF][START_REF] Ashok | The electrical resistance of a graph captures its commute and cover times[END_REF] relating the resistance distance and random walks in graphs. Indeed, let us consider a hypernode graph h where all hypernodes are singleton sets. The hypernode graph can be viewed as a graph and every hyperedge h that contains a node i is an edge {{i}, {k}} with k ∈ N . Thus, we have w h (i) = w h (k) = W i,k and Equation ( 10) can be rewritten as

V j (i, j) = k∈N W i,k deg(i) (1 + V j (k, j)) .
Thus, V j (i, j) can be interpreted as the hitting-time distance from i to j (average number of steps needed by a random walker to travel from i to j). Consequently, Proposition 10 states that, when h is a graph, the distance d 2 (i, j) between two nodes i and j is equal to the commute-time distance between i and j divided by the overall volume, which is the expression found in [START_REF] Klein | Resistance distance[END_REF][START_REF] Ashok | The electrical resistance of a graph captures its commute and cover times[END_REF]. Now, let us consider the general case where the hypernode graph is not a graph. Let us define p(h|i) = w h (i) deg(i) and p(k|h, i) = P (h, i, k) w h (k) w h (i) . Then, we can rewrite Equation (10) as

V j (i, j) = h|i∈h p(h|i)   1 + k∈h,k =i p(k|h, i)V j (k, j)   .
The term p(h|i) can be interpreted as a jumping probability from node i to the hyperedge h because p(h|i) is non-negative and h p(h|i) = 1. Let us now consider the term p(k|h, i), it can be shown that k p(k|h, i) = 1. But the term p(k|h, i) is negative when i and k belong to the same end of h. Consequently, the quantity p(k|h, i) can not be interpreted as a jumping probability from node i to node k with the hyperedge h. Therefore, there is no easy interpretation of d 2 in terms of random walks in the hypernode graph because of possible negative values in the expression of d 2 (i, j).

Conclusion

We have introduced hypernode graphs for learning from binary relations between groups in networks. We have defined a spectral theory allowing to model homophilic relations between groups. Properties of Laplacians and kernels have allowed us to define learning algorithms and an application to skill rating for multiple players games has been presented. This paper opens many research problems. From a machine learning perspective, we hope that the model will open the way to solving new learning problems in networks. From a fundamental perspective, many questions must be investigated. For instance, we have seen that connectivity properties and random walks in hypernode graphs need more research. Also, the definition of cuts in hypernode graphs and the generalization of the Max-Flow Min-Cut theorem are open questions.

Let us now consider the component C of g e that contains the node 1. The indicator function 1 C is smooth on g e but we have (1 C (1) + 1 C (2) -1 C (3) -1 C (4)) 2 = (1 + 0 -0 -0) 2 = 0 .

Consequently, 1 C does not satisfy (C), which violates the second condition and concludes the proof.

V j (i, j) + V i (j, i) = Vol(h)(∆ † i,i + ∆ † j,j -∆ † i,j -∆ † j,i ) = Vol(h)(∆ † i,i + ∆ † j,j -2∆ † i,j ) = Vol(h)d 2 (i, j) , which concludes the first part of the proof.

It remains to show that V j satisfies Equation (10). By definition of V j , we get that, for every node i, V j (i, i) = 0. Let us now consider i = j and let us consider f in S j , i.e., a solution of ∆f = In j . As i = j, we have d(i) = e T i In j . Since f ∈ S j , we have In j = ∆f = G T Gf so we can rewrite the previous equality as d(i) = e T i G T Gf = (Ge i ) T (Gf ). Now, because of (1), we have G1 = 0, thus Gf = G(f -f (j)1). This leads to

d(i) = (Ge i ) T (G(f -f (j)1) = h (Ge i ) (h) • (G(f -f (j)1) (h) (13)
For any oriented hyperedge h = (s h , t h ), we define h (i) to be equal to 1 if i ∈ t h and -1 if i ∈ s h (0 if the node i. does not belong to the hyperedge). We can observe that we always have P (h, i, j) =h (i) h (j) . We have for all i, P (h, i, i) = -1 and h w h (i) = d(i) so the previous equality can be rewritten under the form

d(i) = V j (i, j)d(i) + h|i∈h w h (i)    k∈h,k =i (-P (h, k, i)) w h (k)V j (k, j)    .
Hence, we get the linear system V j (i, j) = 1 + h|i∈h w h (i)

d(i)    k∈h,k =i P (h, k, i) w h (k)V j (k, j)    = h|i∈h w h (i) d(i)    1 + k∈h,k =i P (h, k, i) w h (k) w h (i) V j (k, j)    ,
which concludes the proof.

h

  Figure 1: A hypernode graph with three hyperedges. A hyperedge is represented by a rectangle. Nodes in one of the two hypernodes are connected to the same side of the rectangle.

Figure 2 :

 2 Figure 2: Hyperedge h for a game A = {1, 2, 3} against B = {4, 5} where B wins. The outcome node O is added to A (the loser). The lazy node is added to B (the winner). All weights are set to 1 except the weight for Z chosen in order to ensure the Equilibrium condition.

  minimize s s T ∆s + µσ(s p ) 2 subject to for every outcome node O, s(O) = o; for every lazy node Z , s(Z)

  minimize s s T ∆ µ s subject to for every outcome node O, s(O) = o; for every lazy node Z , s(Z) = 0 (7)

Figure 3 :

 3 Figure 3: [left] Distribution of the number of players against the number of played games;[right] Average percentage of players in Γ u which are involved in some game in Γ l

  . Parameters of Elo and TrueSkill are the default parameters of Hamilton (2012) (K = 32 for Elo, µ0 = 25, β = 12.5, σ = 8.33 and τ = 0.25 for TrueSkill).

Figure 5 :

 5 Figure 4: Predictive error for Double Tennis dataset

Laplacians Proposition 3

 3 The class of hypernode graph Laplacians is the class of symmetric positive semi-definite real-valued matrices M such that 1 ∈ Null(M ), where Null(M ) is the null space of M . Proof It is easy to verify from the Laplacian definition that a Laplacian M is a symmetric positive semi-definite real-valued matrix such that 1 ∈ Null(M ). Conversely, let us consider the following algorithm Algorithm 2 Construction of a hypernode graph Input: M ∈ R n×n symmetric positive semi-definite such that 1 ∈ Null(M ) 1: Set V to {1, . . . , n}; set H to ∅ 2: Compute a square root decomposition M = G T G 3: for each row r of G do 4:

Figure 6 :

 6 Figure 6: Hypernode graph h and a candidate star expansion.

Figure 9 :

 9 Figure 9: The reduced signed graph with adjacency matrix W of Figure 8 of the two L-equivalent hypernode graphs shown in Figure 7.

  , we can develop the expression of d(i) from Equation (13) and writed(i) = h w h (i) h (i) •   k∈h w h (k) h (k)(f (k) -f (j)) h, k, i)) w h (k)(f (k) -f (j)) (h, k, i)) w h (k)V j(k, j)

Appendix A. Regularized Hypernode Graph for Skill Rating

Formally, we assume a numbering of V such that V = {1, . . . , N } where N is the total number of nodes, the first n nodes are the player nodes followed by the t lazy nodes, then followed by the outcome nodes, that is, V = {1, . . . , n}∪{n+1, . . . , n+t}∪{n+t+1, . . . , N }. Let ∆ be the unnormalized Laplacian of h, and let s be a real-valued node function on h, s can be seen as a real vector in R N where the first n entries represent the skills of the n players. Then, Problem 6 can be rewritten as minimize s∈R N s T ∆s + µσ(s p ) 2 subject to ∀n + 1 ≤ j ≤ n + t, s(j) = 0 (for lazy nodes)

where µ is a regularization parameter and s p denotes the vector of player skills (s(1), . . . , s(n)). In order to apply graph-based semi-supervised learning algorithms using hypernode graph Laplacians, we now show that the regularized optimization problem can be rewritten as an optimization problem for some hypernode graph Laplacian. For this, we will show that it suffices to add a regularizer node in the hypernode graph h. First, let us recall that if s is the mean of the player skills vector s p = (s(0), . . . , s(n)), then, for all q ∈ R, we have

Thus, in the problem 6, we can instead minimize s T ∆s + µ n n i=1 (s(i) -q) 2 over s and q. We now show that this can be written as the minimization of r T ∆ µ r for some vector r and well chosen hypernode graph Laplacian ∆ µ . For this, let us consider the p × N gradient matrix G of the hypernode graph h associated with the set of games Γ, and let us define the

, where B is the n × (N + 1) matrix defined by, for every 1 ≤ i ≤ n, B i,i = -1, B i,N +1 = 1, and 0 otherwise.

The matrix G µ is the gradient of the hypernode graph h µ obtained from the hypernode graph h by: add a regularizer node N + 1; add, for every player node, a new hyperedge between the player node and the regularizer node N + 1 with node weights µ/n. Note that such a hyperedge can be viewed as an edge with edge weight µ/n.

Let us denote by r the vector (s(0), . . . , s(N ), q), then since ∆ = G T G, we can write r T ∆ µ r subject to ∀n + 1 ≤ j ≤ n + t, r(j) = 0 (for lazy nodes)

Appendix B. Expressiveness of the hypernode graph Laplacian

Let us consider the simple hypernode graph h first introduced in Figure 6 and recalled in Figure 10 below.

Proposition 11 There do not exist a finite graph g whose node set contains {1, 2, 3, 4} and that satisfies the conditions: 1. All the smooth functions on g satisfy (C).

Any function that satisfies (C) can be extended to a smooth function on g.

Proof Let us define S = {1, 3} and denote by ∆ the Laplacian matrix of h. The indicator vector 1 S is in Null(∆) and, thus, define a smooth function on h. Let us consider a graph g e = (N e , E e ) whose nodeset N e contains the nodeset N = {1, 2, 3, 4} and that satisfies the two conditions presented in Proposition 11. We denote by ∆ e the Laplacian matrix of g e . Because of the first condition, the function 1 S can be extended to a smooth function f e on g e .

Since g e is a graph, we know that Null(∆ e ) is spanned by the indicator vectors of the connected components of g e (see for instance Von Luxburg (2007)). Since f e (1) = 1 S (1) = f e (2) = 1 S (2), 1 and 2 must be in different components in g e . Note that the same holds with f e (1) = f e (4), f e (3) = f e (2) and f e (3) = f e (4). By following a similar reasoning with S = {1, 4}, we eventually deduce that the four original nodes 1, 2, 3 and 4 must be in distinct components in g e .

Appendix C. Proof of Lemma 9

Proof

Let us consider a hypergraph h with Laplacian ∆ such that Null(∆) = Span(1). Therefore, R n is the direct sum of the space Span(1) and of the space Null(∆) ⊥ . Let us consider f ∈ R n , it can be written f = µ1 + g where g ∈ Null(∆) ⊥ . Thus, we have ∆f = ∆(µ1 + g) = ∆g.

Let us suppose that f satisfies ∆f = In j , we deduce that ∆g = In j . As seen above, we know that the operator ∆ † ∆ is the orthogonal projector operator on Null(∆) ⊥ . Since g ∈ Null(∆) ⊥ , we have g = ∆ † ∆g = ∆ † In j . Hence we can write f under the form µ1+∆ † In j .

Conversely, let us consider f = µ1+∆ † In j with µ ∈ R. We have ∆f = ∆ µ1 + ∆ † In j = µ∆1 + ∆∆ † In j . From the definition of In j , we deduce that In j ∈ Null(∆) ⊥ . Since ∆ is symmetric, ∆∆ † is also the orthogonal projector on Null(∆) ⊥ and thus ∆∆ † In j = In j . Since ∆1 = 0, we get ∆f = In j which concludes the proof.

Appendix D. Proof of Proposition 10

Proof First, we show that Ω(i, j) =

V j (i,j)+V i (j,i) Vol(h)

. For that, let us develop the expression of V j (i, j) obtained above V j (i, j) =(e i -e j ) T ∆ † In j

where

We can observe that, for every node i and j in N , we have R i (i, j) + R j (j, i) = 0. Thus, we can write


Daniel Meneveaux

Benjamin Bringier

Emmanuelle Tauzia

Mickaël Ribardière

Lionel Simonot

Isotropic and Anisotropic Interfaced Lambertian Microfacets

Specular microfacet distributions have been successfully employed by many authors for representing materials glossiness, and they are generally combined with a Lambertian term that accounts for the colored aspect. Such a representation makes use of the Fresnel reflectance factor at the interface, but the transmission factor is often ignored. In addition, the generalization to microfacet distributions with a more general reflectance is known to be complex, since it requires to solve an angular integral that has no analytical solution. This paper proposes a complete framework for physically handling both reflection and transmission with microfacet distributions. First, we show how transmission a↵ects reflectance of an interfaced Lambertian model, and provide an analytical description of an individual microfacet reflectance. Second, we describe a method for handling distributions of such microfacets in any physically based Monte-Carlo rendering systems, as well as an approximate representation (which includes light multiple reflexions based on V-cavities representation) that can be used for GPU applications or for measured data fitting. Our approach generalizes several previous models, including flat Lambertian materials as well as specular and Lambertian microfacets. The result section illustrates the wide range of materials that can be possibly taken into account with this representation.

Introduction

Realistic and physically plausible representations of reflectance have been studied for a long time in optics and computer graphics. Yet, many questions still remain open because the mathematical models dedicated to Bidirectional Reflectance Distribution Functions (BRDFs) must fulfill various requirements, such as physical plausibility, ability to naturally represent existing/measured materials, visual realism, computation time e ciency for producing computer generated images, material control for users and artists, etc.

The wide variety of empirical models, introduced for instance by Phong [START_REF] Phong | Illumination for computer generated pictures[END_REF], Blinn-Phong [START_REF] Blinn | Models of light reflection for computer synthesized pictures[END_REF], or the revised version proposed by Lewis [START_REF] Lewis | Making shaders more physically plausible[END_REF], o↵ers a competitive and intuitive way to define surface optical properties, but the parameters used do not correspond to actual physical characteristics and they have limited expressiveness. Other approaches favor function bases whose coe cients can be obtained from measurements using projection or parameter fitting, e.g. with wavelets [START_REF] Claustres | BRDF measurement modelling using wavelets for e cient path tracing[END_REF], spherical harmonics [START_REF] Westin | Predicting reflectance functions from complex surfaces[END_REF], or Phong lobes [START_REF] Lafortune | Non-linear approximation of reflectance functions[END_REF]. Unfortunately with such representations, both the nature and the large number of coe cients make it di cult for users to tune material appearances manually. Finally, physically based BRDF models [START_REF] Cook | A reflectance model for computer graphics[END_REF][START_REF] Ward | Measuring and modeling anisotropic reflection[END_REF][START_REF] Oren | Generalization of Lambert's reflectance model[END_REF][START_REF] Ashikhmin | A microfacet-based BRDF generator[END_REF][START_REF] Walter | Microfacet models for refraction through rough surfaces[END_REF][START_REF] Jakob | A comprehensive framework for rendering layered materials[END_REF][START_REF] Kurt | An anisotropic brdf model for fitting and monte carlo rendering[END_REF][START_REF] Wu | Microfacet-based interference simulation for multilayer films[END_REF] are often preferred because they are defined by parameters that correspond to physical aspects of materials, such as surface roughness, chromaticity, or refractive indices. Many of them rely on a statistical distribution of microfacets associated with their reflectance property.

BRDFs based on microfacet distributions have been introduced in computer graphics for the representation of glossy reflections on rough surfaces [START_REF] Cook | A reflectance model for computer graphics[END_REF][START_REF] Torrance | Theory for o↵-specular reflection from roughened surfaces[END_REF]. The formulation corresponds to a product of Fresnel reflectance, distribution and attenuation factors. Materials can be metallic (complex refractive index), highly reflecting light [START_REF] Callet | Physically based rendering of metallic paints and coated pigments[END_REF][START_REF]From fundamental optical properties of materials to visual appearance: metals, alloys, gilts, glasses, paints, photonics and biophotonics[END_REF], or dielectric (real refractive index) which only partly reflect light. In this latter case, transmission should be accounted for. For instance, Walter et al. [START_REF] Walter | Microfacet models for refraction through rough surfaces[END_REF] propose to handle refraction of rough transparent dielectrics. With paints, plastics, ceramics, and so on, light goes through the dielectric interface and scatters within the material body before re-exiting the surface. In most models, this e↵ect is approximated using a Lambertian term that does not depend on the specular aspect of the surface.

However, as pointed out by several authors, the matte component of rough materials is not constant and often exhibit backscattering [START_REF] Oren | Generalization of Lambert's reflectance model[END_REF][START_REF] Van Ginneken | Di↵use and specular reflectance from rough surfaces[END_REF][START_REF] Kienle | 250 years lambert surface: does it really exist?[END_REF]. The Lambertian term should at least account for BRDF specular albedo to deduce the remaining body reflectance behavior [START_REF] Kelemen | A microfacet based coupled specular-matte BRDF model with importance sampling[END_REF][START_REF] Shirley | A practitioners' assessment of light reflection models[END_REF]. Though physically plausible, this latter representation does not explicitly correlate reflectance with the substrate roughness, and backscattering is not handled.

Simonot [START_REF] Simonot | A photometric model of di↵use surfaces described as a distribution of interfaced Lambertian facets[END_REF] introduces interfaced Lambertian microfacets for physically handling the transmission through a Fresnel interface and multiple reflections with a Lambertian substrate. This theoretical model has been originally defined with a Gaussian distribution of normals, without taking masking and shadowing into account, and numerical integration is required because microfacets are not only specular.

In this paper, we extend Simonot's work and make it possible to use interfaced Lambertian microfacets in rendering systems. We render the e↵ects of this model, and confront the obtained results with other microfacet-based BRDFs. Figure 1 illustrates some results produced with the approaches described in this paper. While Jakob et al. [START_REF] Jakob | A comprehensive framework for rendering layered materials[END_REF] define materials explicitly with rough translucent layers, interfaced Lambertian microfacets are closer to Oren and Nayar's approach [START_REF] Oren | Generalization of Lambert's reflectance model[END_REF], and handle backscattering. In our case, contrary to Jakob et al., explicit lighting simulation within layers is not required since interreflections between the substrate and the interface are analytically described. Our representation can be directly integrated in any Monte-Carlo based lighting simulation renderer, and we also propose an approximate model which includes light multiple reflections within V-cavities.

More specifically, the contributions of this paper consist of:

• a consistent framework dedicated to microfacet BRDFs, that integrates a homogeneous matte body with a rough interface; it extends previous models for opaque materials, including flat specular or Lambertian, as well as a range from specular to Lambertian rough surfaces, with various distributions;

• a practical solution for straightforwardly integrating such models within Monte-Carlo-based rendering systems, including importance sampling and anisotropy;

• a method based on path tracing for handling light multiple reflection contributions within Vcavities;

• an approximate model that can be employed for direct lighting, GPU, and measurement fitting, based on the Oren-Nayar representation, and which handles light inter-reflections between microfacets.

Our results show that this model accounts for many reflection phenomena, including backscattering, anisotropy, and BRDF darkening at incoming and outgoing grazing angles, that can be observed on some measured materials. These e↵ects are naturally included in the model since it physically handles balance between specular and body reflections.

The remainder of the paper is organized as follows. Section 2 presents the theory of microfacet BRDF models and introduces notations. Section 3 describes interfaced Lambertian microfacet BRDFs. Section 4 explains how scattering microfacets can be employed in various configurations in rendering systems. Section 5 presents our results for various cases, including comparisons and fittings from measured BRDFs. Conclusion and future work are provided in Section 6.

Microfacet-Based BRDF

This paper focuses on microfacet-based models for opaque materials, as described by Torrance and colleagues [START_REF] Cook | A reflectance model for computer graphics[END_REF][START_REF] Torrance | Theory for o↵-specular reflection from roughened surfaces[END_REF], Oren and Nayar [START_REF] Oren | Generalization of Lambert's reflectance model[END_REF], Ashikhmin et al. [START_REF] Ashikhmin | A microfacet-based BRDF generator[END_REF], Bagher et al. [START_REF] Bagher | Accurate fitting of measured reflectances using a Shifted Gamma micro-facet distribution[END_REF], Dupuy et al. [START_REF] Dupuy | Linear e cient antialiased displacement and reflectance mapping[END_REF], or Heitz [START_REF] Heitz | Understanding the masking-shadowing function in microfacet-based brdfs[END_REF]. They have been successful at representing various families of materials with only a few parameters clearly related to surface observations and intuitive to manipulate. So far, they have been described using (i) a microfacet distribution, (ii) an elementary reflectance (purely specular or Lambertian), and (iii) geometric factors accounting for light shadowing and masking.

Definitions and Notations

The notations used in this paper are given in Table 1, they correspond to those proposed by Nicodemus et al. [START_REF] Nicodemus | Geometrical consideration and nomenclature for reflectance[END_REF]. Any direction u is defined with a polar angle ✓ u and an azimutal angle ' u expressed according to the macrosurface normal n.

Let us consider a surface sample of normal n lit by a collimated light source, from direction i, as illustrated in Figure 1. The radiance reflected toward an observer o is given by the BRDF f (i, o, n). It is defined as the ratio between the elementary radiance dL(i, o, n) reflected by the surface in direction o and the incident irradiance dE(i, n) coming from direction i:

f (i, o, n) = dL(i, o, n) dE(i, n) . (1)
A BRDF should be physically plausible, i.e. should satisfy Helmoltz reciprocity, i.e.

f (i, o, n)=f (o, i, n) and energy conservation, i.e. R ⌦ + f (i, o, n)|on|d! o  1, 8i 2 ⌦ + , with |on|= o•n if o•n > 0
and 0 otherwise. The microfacet representation of a BRDF f (i, o, n) corresponds to a statistical description of microfacet distribution. Given the BRDF f µ (i, o, m) of an individual microfacet associated with normal m, its contribution is weighted by the distribution D(m) and a geometric attenuation factor (or GAF) G(i, o, m). The distribution function D(m) defines the surface roughness, indicating the proportion of microfacets oriented according to a given orientation m; the attenuation factor G(i, o, m) determines the portion of a microfacet of normal m visible from both the light source and the observer, with a major influence at grazing angles. Many authors have studied the use of various distributions and geometric attenuation factors [7, 9-11, 15, 20, 23, 24, 27], which have to be carefully chosen together [START_REF] Ashikhmin | A microfacet-based BRDF generator[END_REF][START_REF] Heitz | Understanding the masking-shadowing function in microfacet-based brdfs[END_REF]:

f (i, o, n) = Z ⌦ + |im| |in| f µ (i, o, m) |om| |on| D(m)G(i, o, m)d! m . (2)
All microfacets are supposed to be oriented in the upper hemisphere, with D(m) = 0 if m• n0; the projected areas of all microfacets have to be equal to the macroscopic surface, i.e.

R ⌦ + D(m)|mn|d! m = 1.
In the case of purely specular microfacets, Equation 2 simplifies to [START_REF] Cook | A reflectance model for computer graphics[END_REF][START_REF] Walter | Microfacet models for refraction through rough surfaces[END_REF][START_REF] Torrance | Theory for o↵-specular reflection from roughened surfaces[END_REF]:

f (i, o, n) = F (i, h)D(h)G(i, o, h) 4|ih||oh| , (3)
where h = i+o ||i+o|| is the half angle vector between i and o, and F (i, h) corresponds to Fresnel reflectance, depending on the medium refractive index n i (see Table 1). This equation defines the glossy aspect of the surface.

Discussion

On the one hand, rough metallic surfaces are almost opaque (Fresnel reflectance values are close to 1), and n 1 is defined by a complex value, with reflectance variations according to wavelengths [START_REF] Callet | Physically based rendering of metallic paints and coated pigments[END_REF].

f (i, o, n) Macroscopic BRDF f µ (i, o, m) Microfacet BRDF D(m) Microfacet normal distribution x , y
Roughness parameters given in local coordinate system for dealing with anisotropy

G(i, o, m) Geometric attenuation factor n 0 External refractive index n 1
Internal refractive index

n i = n 1 n 0 Ratio of refractive indices F (i, m)
Fresnel term, interface reflectance

T (i, m) Interface transmittance T (i, m) = 1 -F (i, m)
Such materials are handled with specular microfacets directly since the resulting chromaticity and glossiness are completely handled by the Fresnel term and the chosen distribution. On the other hand, dielectrics correspond to real values of n 1 , with negligible wavelength dependency (the wavelength dependency is only related to body reflection). For common values of n i (i.e. 1  n i < 3), the Fresnel reflectance is low except at grazing angles. With energy conservation, the interface transmittance is significant although it is generally ignored [START_REF] Cook | A reflectance model for computer graphics[END_REF] or roughly approximated [START_REF] Kelemen | A microfacet based coupled specular-matte BRDF model with importance sampling[END_REF]. The model presented in this paper handles this latter type of material and deals with transmission of light through the interface represented by a real refractive index and multiple light reflections between a Lambertian substrate and the interface.

The strong assumption of considering only pure Lambertian or pure specular microfacets restricts the gamut of materials that can be reproduced. For instance, the objects represented by the Oren-Nayar model [START_REF] Oren | Generalization of Lambert's reflectance model[END_REF] cannot exhibit any glossy appearance, while the Cook-Torrance model [START_REF] Cook | A reflectance model for computer graphics[END_REF] only employs a constant Lambertian term, that does not account for backscattering with rough materials [START_REF] Oren | Generalization of Lambert's reflectance model[END_REF]. These models should not be combined for taking into account simultaneously glossy and matte aspects for three main reasons: first, Oren and Nayar use a distribution which is very rarely used in other models; second, interfaced Lambertian microfacets exhibit chromaticity variations according to the interface refractive index discontinuity, that cannot be handled using a simple combination of models; third, energy conservation is not handled with such a combination since body reflection should also depend on n

i . The choice of D(m), G(i, o, m), and f µ (i, o, m) may provide a wide range of physically plausible materials, as long as the combinations are carefully chosen [START_REF] Ashikhmin | A microfacet-based BRDF generator[END_REF][START_REF] Walter | Microfacet models for refraction through rough surfaces[END_REF][START_REF] Heitz | Understanding the masking-shadowing function in microfacet-based brdfs[END_REF][START_REF] Bourlier | One and two-dimensional shadowing functions for any height and slope stationary uncorrelated surface in the monostatic and bistatic configurations[END_REF].

Interfaced Lambertian BRDF

The model presented in this paper defines rough surfaces (Figure 2) made up with a Lambertian substrate of parameter K d (which depends on the wavelength ), covered with a flat interface corresponding to a refractive index discontinuity n i (a real index for dielectrics) between these two media. It accounts for both surface and body reflections.

Single Microfacet (or flat surface)

Interfaced Lambertian microfacets are associated with a BRDF f µ (i, o, m) defined by a pure specular interface reflection f µ s (i, o, m) combined with the material volume body di↵use reflection f µ b (i, o, m):

f µ (i, o, m) = f µ s (i, o, m)+ f µ b (i, o, m). (4)
Since microfacet interfaces are considered as perfectly flat, the surface reflection component

f µ s (i, o, m)
is [START_REF] Walter | Microfacet models for refraction through rough surfaces[END_REF]:

f µ s (i, o, m) = F (i, h) ! m (h, m) 4|ih| 2 , ( 5
)
where h is the bissector direction of i and o, ! m is the dirac distribution associated with the elementary solid angle d! m around m, so that for any solid angle ⌦:

Z ⌦ ! m (s, m)d! m = ( 1 if s 2 ⌦ 0 otherwise (6)
and F (i, h) is the reflectance of the interface (for unpolarized irradiance), corresponding to the Fresnel's law between an outer medium of index n 0 and a medium of index n 1 :

F (i, h) = 1 2 (g -c) 2 (g+ c) 2 ⇢ 1+ [c(g+ c) -1] 2 [c(g -c)+ 1] 2 , (7)
where

c = i•h and g 2 = n 2 i + c 2 -1.

Lambertian substrate

Fresnel interface

Incoming light

First specular reflection Scattering after multiple reflections Substrate-interface multiple reflections Body scattering (Figure 3) should account for the first interface transmission T (i, m) = 1 -F (i, m) of light, followed by a Lambertian reflection due to the substrate (of reflectance K d ), inner multiple interactions between interface and substrate 1/(1 -K d r i ), and final transmissions toward the outgoing direction, of transmittance T (o, m):

f µ b (i, o, m) = 1 ⇡n 2 i T (i, m)T (o, m) K d (1 -K d r i ) , (8)
where r i is the internal reflectance on the flat interface lit by a Lambertian source coming from the medium of index n 1 :

r i = Z ⇡ 2 ✓u = 0 F i (u, m) sin(2✓ u )d✓ u , (9)
with F i the microfacet inner Fresnel reflectance. We use the analytical expression

n 2 i (1 -r i ) = 1 -r e
with r e provided by Molenaar et al. [START_REF] Molenaar | Determination of kubelka-munk scattering and absorption coe cients by di↵use illumination[END_REF]:

r e = 1 2 - 2n 3 i (n 2 i + 2n i -1) (n 2 i + 1)(n 4 i -1) + (n i -1)(3n i + 1) 6(n i + 1) 2 + ( 10
)
8n 4 i (n 4 i + 1) (n 2 i + 1)(n 4 i -1) 2 ln(n i )+ n 2 i (n 2 i -1) 2 (n 2 i + 1) 3 ln ✓ n i -1 n i + 1 ◆
This model considers an interface between two media, i.e. with a refractive index discontinuity, thus not subject to interferences.

With interfaced Lambertian microfacets in Equation 8, 1/n 2 i expresses an extension of the light beam from a medium of index n 1 to a less refringent medium of index n 0 < n 1 . The multiple reflections between the interface and the substrate lead to chromatic variations depending on parameter n i . Figure 4a actually varies according to n i , providing darker and more saturated colors. Figure 4c illustrates colored disks with these variations, according with n i . With achromatic objects, some variations are also noticeable, since increasing values of n i also tend to provide higher transmission, and consequently darker body reflection. Figure 5 shows the curves and a rendered 3D object with materials corresponding to a flat interfaced Lambertian surface, with varying values of n i . Light inter-reflections between the interface and the substrate a↵ect the object appearance depending on K d /(1 -K d r i ). Reflected radiance deacreases when n i increases, and when observation angles become grazing, contrary to Lambertian or rough Lambertian materials.

Rough Interfaced Lambertian BRDF

Using f µ from Equation 4 in Equation 2, the global BRDF remains an expression composed of a glossy and a matte components (resp. f s and f b ):

f (i, o, n) = f s (i, o, n)+ f b (i, o, n). (11)
The specular component that provides glossiness, corresponds to Equation 3 (coming from Equations 2 and 5), while from Equations 2 and 8, the body reflection becomes:

f b (i, o, n) = K d ⇡n 2 i (1 -r i K d ) ⇥ (12) Z ⌦ + T (i, m)T (o, m)D(m)G(i, o, m) |im| |in| |om| |on| d! m .
With this representation, several noticeable characteristics can be underlined:

• when K d = 0, this formulation corresponds to purely specular microfacets; • when n i = 1, then T = 1 and r i = 0, it corresponds to purely Lambertian microfacets; • when = 0, the BRDF is a flat interfaced Lambertian surface;

• when n i = 1 and = 0, the model is equivalent to a flat Lambertian material. In this paper, we have employed various distributions and geometric attenuation factors from the literature.

Gaussian distribution of normals is essentially used in the Oren-Nayar model [START_REF] Oren | Generalization of Lambert's reflectance model[END_REF]. This distribution is employed in this paper only for comparison purposes, with the following isotropic formulation:

D(m) = C e (-✓ 2 m /2 2 ) / |mn|, (13)
where is the roughness parameter (expressed in radians) and

1/C = 2⇡ R ⇡/2 ✓m = 0 e (-✓ 2 m /2 2
) sin ✓ m d✓ m . Using such a distribution remains di cult in practice because of the estimation of C. The correspondance with Beckmann or GGX distributions requires 0 = p 2 . We employ this distribution as Oren and Nayar [START_REF] Oren | Generalization of Lambert's reflectance model[END_REF], in association with the Torrance-Sparrow geometric attenuation factor [START_REF] Torrance | Theory for o↵-specular reflection from roughened surfaces[END_REF], which is mathematically well-posed [START_REF] Heitz | Understanding the masking-shadowing function in microfacet-based brdfs[END_REF]. This factor assumes V-cavities much longer than large:

G(i, o, m) = max ✓ 0, min ✓ 1, 2|in||mn| |im| , 2|on||mn| |om| ◆ ◆ . ( 14
)
However, this assumption is somehow physically unrealistic and prone to discontinuities at grazing angles.

Beckmann distribution is often referenced in the literature [START_REF] Cook | A reflectance model for computer graphics[END_REF][START_REF] Ashikhmin | A microfacet-based BRDF generator[END_REF][START_REF] Kelemen | A microfacet based coupled specular-matte BRDF model with importance sampling[END_REF]. It corresponds to a Gaussian distribution of microfacet slopes. Ward [START_REF] Ward | Measuring and modeling anisotropic reflection[END_REF] proposes a derived anisotropic distribution, and Kurt et al. [START_REF] Kurt | An anisotropic brdf model for fitting and monte carlo rendering[END_REF] provides the following formulation, that includes the correct normalization:

D(m) = e ✓ tan 2 ✓m ✓ cos 2 'm 2 x + sin 2 'm 2 y ◆◆ ⇡ x y |mn| 4 (15)
where x and y correspond to the roughness parameters (standard deviations of slopes) associated with a local coordinate system. The well-known Cook-Torrance model corresponds to the isotropic form of this distribution (for x = y ), associated with the Torrance-Sparrow GAF. Several authors have discussed attenuation factors [START_REF] Ashikhmin | A microfacet-based BRDF generator[END_REF][START_REF] Heitz | Understanding the masking-shadowing function in microfacet-based brdfs[END_REF][START_REF] Smith | Geometrical shadowing of a random rough surface[END_REF] since the Torrance-Sparrow model makes the assumption of Vcavities that cannot represent actual physical surfaces. One of the most used representations [START_REF] Ashikhmin | A microfacet-based BRDF generator[END_REF][START_REF] Walter | Microfacet models for refraction through rough surfaces[END_REF][START_REF] Heitz | Understanding the masking-shadowing function in microfacet-based brdfs[END_REF] has been proposed by Smith [START_REF] Smith | Geometrical shadowing of a random rough surface[END_REF], making the assumption of a separable form, so that masking and shadowing are independent:

G(i, o, m) = G 1 (i, m)G 1 (o, m). (16)
Given that the GAF should be independent from the normal m of a microfacet, i.e. G 1 (u, m) does not correlate with any other microfacet, even nearby, it leads to the following expression:

G 1 (u, m) = ( G 1 (u) when u•m > 0 0 w h e nu•m  0, (17)
with

G 1 (u) = u•n R ⌦ + u D(m)|um|d! m . (18)
Ashikhmin et al. [START_REF] Ashikhmin | A microfacet-based BRDF generator[END_REF] suggest a prior numerical integration for all directions u. In the case of Beckmann distribution, this integral has an analytical expression [START_REF] Heitz | Understanding the masking-shadowing function in microfacet-based brdfs[END_REF] with G 1 (u) = 1/(1+ ⇤(u)), where:

⇤(u) = erf(a) -1 2 + 1 2a p ⇡ e -a 2 , (19)
and

a = 1/(↵ u tan ✓ u ), ↵ u = q 2 x cos 2 ' u + 2 y sin 2 ' u .
In the remainder of this paper, we will make the distinction between Gaussian and Beckmann distributions as described in the two previous paragraphs, although Beckmann also corresponds to a Gaussian.

GGX distribution, described in [START_REF] Walter | Microfacet models for refraction through rough surfaces[END_REF][START_REF] Trowbridge | Average irregularity representation of a rough surface for ray reflection[END_REF], has been extended to anisotropic distribution [START_REF] Heitz | Understanding the masking-shadowing function in microfacet-based brdfs[END_REF]:

D(m) = 1 ⇡↵ x ↵ y |mn| 4 ⇣ 1+ tan 2 ✓ m ⇣ cos 2 'm 2 y + sin 2 'm 2 y ⌘⌘ 2 . (20)
The associated GAF is also obtained using G 1 (u) = 1/(1+ ⇤(u)), with:

⇤(u) = 1 + p 1+ 1/a 2 2 , ( 21
)
where

a = 1/(↵ u tan ✓ u ), ↵ u = q 2 x cos 2 ' u + 2 y sin 2 ' u .

Light Inter-reflections

Light multiple reflections between microfacet interface and body (Figure 3) are analytically handled by our model (Equations 9 and 8). In addition, light reflections between microfacets should also be accounted for. Equation 2 stands for single-bounce reflections; to the best of our knowledge, multiplebounce inter-reflections with microfacets have only been handled by Oren and Nayar [START_REF] Oren | Generalization of Lambert's reflectance model[END_REF], who consider the influence of two bounces of reflections with pure Lambertian microfacets, under the assumption of V-cavities (Figure 6).

This problem is di cult since it requires to solve the infinite dimensional integral rendering equation. We have implemented a numerical Monte-Carlo integration for estimating light inter-reflections with Vcavities and interfaced Lambertian microfacets given in Equation 4. Given an incident light direction i, the BRDF is the sum of all the contributions

f 1 = f 1 + f 2 + • • •
, where f k corresponds to the observed radiance (or reflectance with unitary irradiance) after k reflections: where dL 1 is the total reflected radiance in direction o for a collimated light beam coming from direction i. This formulation takes into account all light interactions with all possible microfacets complying with distribution D. R 1 (i, o, m) = R S L(x)dx is the total reflected radiance for all light paths x (from the path domain S and L(x) is the path contribution function) with at least one bounce inside a V-cavity represented by its normal m. L(x) is estimated using path tracing.

f 1 (i, o, n) = dL 1 (i, o, n) |in| = Z ⌦+ R 1 (i, o, m) |in| |om| |on| D(m)d! m , (22)
We propose to evaluate Equation 22 using Monte-Carlo quadrature. First, a microfacet orientation is chosen according to importance sampling, depending on the probability density function (pdf) D(m)|mn|; second, a point x i is uniformly chosen on the microfacet; third, a light path is built similarly to the photon mapping rendering method [START_REF] Jensen | Realistic Image Synthesis Using Photon Mapping[END_REF], from a randomly chosen direction i. This path hits the point x i from direction i, reflects several times within the V-cavity and finally gets out of the microsurface. The path contribution is collected on a virtual hemispherical sensor, subdivided such that all cells correspond to the same solid angle (Figure 6). With this simulation method, shadowing and masking terms are automatically handled. Moreover, similarly to Oren and Nayar [START_REF] Oren | Generalization of Lambert's reflectance model[END_REF], V-cavities are considered as much longer than large.

Inside each sensor cell, density estimation is used to capture f 1 (i, o, n) where o is the center direction of the cell solid angle:

f 1 (i, o, n) ⇡ P N j = 1 1 |in| L(x j ) |ox j m j | |ox j n| D(m j ) N tot ! o , (23)
where N tot is the total number of samples used to estimate Equation 22, N is the number of paths contributing to the cell, o xj is the outgoing direction of path xj , m j is the sampled microfacet normal and d! o the cell's solid angle. Figure 7 illustrates the e↵ects of taking into account inter-reflections (two and multiple bounces). Simulations are done with one billion light paths and a sensor containing 3240 cells. When the microfacet BRDF is close to Lambertian, inter-reflections play a role in the total reflection model in case of a rough surface. On the other hand, pure specular microfacets produce almost no additional reflected light. Two light bounces are significant for Lambertian microfacets, while further reflections produce only very low additional reflected radiance in most configurations (as already explained in [START_REF] Oren | Generalization of Lambert's reflectance model[END_REF]). Some more results can be found in the joined additional material file. This protocol makes use of V-cavities, with the advantage that it does not require to explicitly construct a surface composed of fixed microfacets, and enables the computation without loss of generality. Using other types of representation complying with Smith's GAF would require some more experiments, which are beyond the scope of this paper.

BRDF Integration and Sampling

Interfaced Lambertian microfacet BRDFs are composed of two terms f s and f b (Section 3.2) that can be processed independently with the rendering equation, due to integration linearity:

L o (x, o, n) = Z ⌦ + L i (x, i, n)f (i, o, n)|in|d! i , (24)
where x is the considered surface element location, L o (x, o, n) corresponds to the outgoing radiance,

L i (x, o, n) is the incident radiance, coming from direction ! i , and f (i, o, n) = f s (i, o, n)+ f b (i, o, n).
The glossy term corresponds to the well-known formulation provided by Equation 3, while the body part is directly derived from Equation 2, using the microfacet body contribution f µ b (Equation 12):

f b (i, o, n) = Z ⌦ + |im| |in| f µ b (i, o, m) |om| |on| D(m)G(i, o, m)d! m . (25)
Unfortunately, Equation 25 requires numerical integration. The next section shows how it can be handled straightforwardly in any Monte-Carlo rendering system, producing correct results without making use of any approximation. We also provide in this paper an approximate isotropic model, that can be used without Monte-Carlo integration, for instance for GPU applications or data fitting.

Monte-Carlo Integration

Most renderers based on ray-tracing approaches make use of Monte-Carlo integration for estimating the rendering equation. In those cases, BRDF importance sampling is used for reducing variance, as pointed out by many authors [START_REF] Walter | Microfacet models for refraction through rough surfaces[END_REF][START_REF] Kelemen | A microfacet based coupled specular-matte BRDF model with importance sampling[END_REF][START_REF] Lafortune | Mathematical models and monte carlo algorithms for physically based rendering[END_REF][START_REF] Veach | Robust monte carlo methods for light transport simulation[END_REF].

Sampling BRDFs requires integrating Equation 25(corresponding to the body contribution), using a stochastic process that only requires choosing microfacet orientations. Although microfacet reflectance and the resulting BRDF are smooth, uniform sampling is not e cient due to the distribution function D that might introduce high frequencies (thin peaks and long tails), resulting in noise. Importance sampling should thus also be used, based on D(m)|mn|, as for path tracing importance sampling.

The specular component (Equation 3) can be processed as in any rendering system, without changing the importance sampling strategy, as it corresponds to the usual Cook-Torrance formulation (with potentially di↵erent distributions and/or attenuation factors).

Choosing between surface and body sampling is the first step. Ideally, f s and f b should be integrated over hemispherical incidence directions to determine each term weight. However, this process would be ine cient in practice. Instead, we propose to apply the process as for a Lambertian surface with a flat interface (Equations 5 and 8). The resulting analytical formulation depends on K d and n i ; it is physically exact with a plane interface, and this assumption is e cient with rough surfaces since the balance between the two terms remains approximately identical.

Let us define R s = r e the total specular reflectance provided by Molenaar, in Equation 11, and the total body reflectance R b = 1

n 2 i (1 -r e ) 2 K d 1 -K d r i ;
We use the proportion between R s and R b for choosing the contribution between specular and body directions. Given a uniform random value ⇠

0 2 [0, 1[, if ⇠ 0 < R s /(R s + R b )
, then sample a specular direction, otherwise, sample a body direction.

Importance sampling the specular/glossy term has been discussed by many authors; Walter et al. [START_REF] Walter | Microfacet models for refraction through rough surfaces[END_REF] provide an overview on the subject. Based on the microfacet distribution (see below for an anisotropic representation), the pdf is given by D(m)|mn|. The weight associated with the outgoing direction is:

w(i) = f s (i, o, m)D(m)G(i, o, m) 4|in||on|p o (i) , (26)
with:

p o (i) = pdf(m) @! h @! o = pdf(m) 1 4|om| . ( 27
)
Body term sampling is more di cult since all microfacet contributions should be accounted for.

We propose to transfer BRDF integration in the rendering system directly. The Monte-Carlo process can be applied to the whole rendering equation, including the analytical BRDF:

L o (o, n) = Z ⌦ + L i (i, n) Z ⌦ + fb (i, o, n)|in|d! m d! i , (28)
where

fb (i, o, n) = |im| |in| f µ b (i, o, m) |om| |on| D(m)G(i, o, m
), as stated in Equation 25. This equation exhibits a four dimensional integral that can also be solved according to Monte-Carlo importance sampling. Therefore, the process consists in selecting the incoming direction i according to a first importance sampling process, based on a pdf equal to |in|/⇡, as well as a microfacet orientation m with on a second importance sampling process based on a pdf equal to D(m)|mn|. Anisotropy is straightforwardly handled by the chosen distribution of D(m).

Using the anisotropic distributions for Beckmann and GGX (Equations 15 and 20), importance sampling of a single microfacet is given by the following sampling. In the following, (⇠ 1 , ⇠ 2 ) are two uniform random variables in [0, 1[ 2 ; the value of ' m is identically sampled for all distributions:

' m = arctan ✓ y x tan(2⇡⇠ 2 ) ◆ . (29)
The value for ✓ m depends on the chosen distribution, reversing the cumulative density function (cdf) associated with pdf(m) = D(m)|mn|.

With the Beckmann distribution:

✓ m = arctan 0 B @ v u u t log(⇠ 1 ) cos 2 'm 2 x + sin 2 'm 2 y 1 C A . (30)
With the GGX distribution:

✓ m = arctan 0 B @ v u u t ⇠ 1 (1 ⇠ 1 ) ⇣ cos 2 'm 2 x + sin 2 'm 2 y ⌘ 1 C A . (31)
Note that importance sampling a Gaussian distribution is not straightforward since the cdf cannot be analytically represented and thus inverted.

Results and Discussion

Our model has been added as a plugin to the Mitsuba renderer [START_REF] Jakob | Mitsuba renderer[END_REF], and the approximate version has been derived as a glsl shader. Table 2 provides the configurations of distributions and attenuation factors we have used in this paper. In the following tables and figures, IL refers to our interfaced Lambertian microfacets. Influence of interface on body reflection is illustrated in Figure 8, according to the choice of functions D(m) and G(m). The discontinuity corresponding to the Torrance-Sparrow GAF V-cavities essentially influences the BRDF at grazing incidence and/or observation angles. In addition, when is low and/or when n i is close to 1, the di↵erence between Torrance-Sparrow and Smith GAFs tends to decrease. Note that the higher tail for GGX distribution produces increasing values at grazing observation angles (Figure 8.c).

Conclusion

This paper presents an analysis and implementation methods for interfaced Lambertian microfacet BRDFs. The presented model relies on very few parameters (refractive index, substrate reflectance and surface roughness) and generalizes the Lambert model, as well as specular microfacets and Lambertian microfacets. It is physically plausible and can be straightforwardly integrated in any Monte-Carlobased rendering system using importance sampling with any distribution and geometric attenuation factor. The proposed methodology can be extended to any type of microfacet-based reflectances.

We have also proposed an approximate version (corresponding to two distributions) of this model that can be employed with interactive rendering systems, or fitting BRDF data, for avoiding numerical integration that remains time consuming. It also integrates light inter-reflections between microfacets.

Our model and the methods we propose in this paper do not correspond to a subset of Jakob et al.'s approach [START_REF] Jakob | A comprehensive framework for rendering layered materials[END_REF] due to the interface management which is di↵erent from multi-layers; for instance, pure Lambertian microfacet BRDFs are inherently included in our model, and no light simulation between layers is required.

A limitation of our method concerns the mean to estimate light inter-reflections between microfacets.

The corresponding contributions are handled by our approximate model, with V-cavities approximations, but including this computation within path tracing remains a di cult issue, a deeper analysis would be needed.

An interesting question concerns the correlation between the interface and the substrate roughnesses, since a single-roughness value constrains both specular peaks and body reflection, which is not always observed in measured materials. We also wish to investigate such a representation.

Figure 1 :

 1 Figure 1: Geometry of reflection with our notation.

Figure 2 :

 2 Figure 2: Surface built up with interfaced Lambertian microfacets. The substrate scatters light while the interface provides brightness.

Figure 3 :

 3 Figure 3: Interfaced microfacet body reflection: light is both reflected and transmitted by the interface. Our model also accounts for light multiple scattering.

  provides three spectral reflectances of the Macbeth color checker corresponding to CIE xy values given in Figure 4b for normalized illuminant D65. This diagram shows that chromaticity Spectral reflectance K d ( ) for three colored samples. Chromaticity variations for the three color checker samples given above and the color gamut PAL/SECAM, according to ni for values equal to ni ={1.0, 1.1, 1.2, 1.3, 1.4, 1.5}, with ✓i = 65 . The arrows illustrate the variation direction when ni is increased. ni = 1.0 ni = 1.1 ni = 1.2 ni = 1.3 ni = 1.4 ni = 1.5 (c) Visual aspect of chomaticity variations according to n i for the purple sample, with PAL/SECAM color gamut and ✓ i = 65 .

Figure 4 :

 4 Figure 4: K d ( ) and chromaticity variations according to n i , for three sample colors from the Macbeth color checker.

( a )

 a Body reflection according to refractive index (the specular reflection peak is not illustrated here since it only corresponds to a dirac). ni = 1 ni = 1.2 ni = 1.33 ni = 1.5 (b) Rendered images of a flat interfaced Lambertian BRDF, with various values of ni.

Figure 5 :

 5 Figure 5: (a) BRDF of a flat interfaced Lambertian surface with K d = 0.6 for each component r,v,b, and n i ={1.0, 1.2, 1.33, 1.5}; (b) rendered images of the same BRDFs applied on a 3D object.

Figure 6 :

 6 Figure6: Light inter-reflections within V-cavities: (a) Light path starts using a given direction i and interacts with the microfacets before exiting in direction o; (b) its contribution is collected on a hemispherical sensor subdivided in equal solid angle cells; f 1 (i, o, n) is computed in each cell using density estimation based on all paths contributing to the cell.

6 Figure 7 :

 67 Figure 7: V-cavities reflectance with inter-reflections, Beckmann distribution: f 1 corresponds to a single bounce; f 2 corresponds to the BRDF with two light bounces; f 1 includes all bounces, using the Russian roulette strategy.

6 Figure 8 :

 68 Figure 8: Influence of distributions and GAFs for several values of n i and , illustrated for ✓ i = 60 o (log scale).

Figure 9 :

 9 Figure 9: Comparison between Oren-Nayar coupled with specular microfacets (Gaussian distribution and Torrance-Sparrow GAF) and interfaced Lambertian (IL, same distribution and GAF), with K d = 0.6, n i = 1.5, and various values of (log scale).

Figure 10 :

 10 Figure 10: Isotropic IL elephants (with both the same parameters), rendered using path tracing: (top) n i = 1.5, for glossy materials and various roughnesses, clearly visible on the specular highlight; (bottom) n i = 1.1, with matte objects and backreflections (progressively darkening the right part of the bigger elephant).

Figure 11 :

 11 Figure 11: Appearance variations according to n i and , with Beckmann distribution and Smith GAF. The bottom line corresponds to pure Lambertian microfacets, and the bottom-left image is an almost flat Lambertian surface.

Figure 12 :

 12 Figure 12: Anisotropic interfaced Lambertian microfacets, with (top) n i = 1.5 and (bottom) n i = 1.0.

Table 1 :

 1 Notations used in this paper.

	n	Macroscopic normal unit vector
	m	Microfacet normal unit vector
	i	Incident light unit vector direction
	o	Observer unit vector direction
	h	Half angle unit vector, bisector between i and o
	✓	u	Polar angle associated with direction u
	'

u Azimutal angle associated with direction u d! u Elementary solid angle around direction u ⌦ + Upper hemisphere of direction n

Table 2 :

 2 Configurations of distributions and attenuation factors used in this paper.

	Notation Distribution (D)	Attenuation factor
			(G)
	Gauss, TS Gaussian normals	Torrance-Sparrow
	Beckm, TS Beckmann (slopes)	Torrance-Sparrow
	Beckm, SB Beckmann (slopes)	Smith-Bourlier
	GGX	GGX	GGX GAF