Computing Shapley values in the plane

We consider the problem of computing Shapley values for points in the plane, where each point is interpreted as a player, and the value of a coalition is defined by the area of usual geometric objects, such as the convex hull or the minimum axis-parallel bounding box. For sets of $n$ points in the plane, we show how to compute in roughly $O(n^{3/2})$ time the Shapley values for the area of the minimum axis-parallel bounding box and the area of the union of the rectangles spanned by the origin and the input points. When the points form an increasing or decreasing chain, the running time can be improved to near-linear. In all these cases, we use linearity of the Shapley values and algebraic methods. We also show that Shapley values for the area and the perimeter of the convex hull or the minimum enclosing disk can be computed in $O(n^2)$ and $O(n^3)$ time, respectively. In this case the computation is closely related to the model of stochastic point sets considered in computational geometry, but here we have to consider random insertion orders of the points instead of a probabilistic existence of points.


Introduction
One can associate several meaningful values to a set P of points in the plane, like for example the area of the convex hull or the area of the axis-parallel bounding box. How can we split this value among the points of P ? Shapley values are a standard tool in cooperative games to "fairly" split common cost between different players. Our objective in this paper is to present algorithms to compute the Shapley values for points in the plane when the cost of each subset is defined by geometric means.
Formally, a coalitional game is a pair (N, v), where N is the set of players and v : 2 N → R is the characteristic function, which must satisfy v(∅) = 0. Depending on the problem at hand, the characteristic function can be seen as a cost or a payoff associated to each subset of players.
In our setting, the players will be points in the plane, and we will use P ⊂ R 2 for the set of players. Such scenario arises naturally in the context of game theory through modeling: each point represents an agent, and each coordinate of the point represents an attribute of the agent. We will consider characteristic functions defined through the area, such as the area of the convex hull. As mentioned before, such area could be interpreted as a cost or a payoff associated to each subset of the points P , depending on the problem.
Shapley values are probably the most popular solution concept for coalitional games, which in turn are a very common model for cooperative games with transferable utility. The objective  Figure 1: Different costs associated to a point set that are considered in this paper. The cross represents the origin. In all cases we focus on the area.
Note that in all the problems we consider monotone characteristic functions: whenever Overview of our contribution. We show that the Shapley values for the AreaConvexHull and AreaEnclosingDisk games can be computed in O(n 2 ) and O(n 3 ) time, respectively. These problems resemble the models recently considered in stochastic computational geometry; see for example [1,2,12,13,14,22,30]. However, there are some key differences in the models. In the most basic model in stochastic computational geometry, sometimes called unipoint model, we have a point set and, for each point, a known probability of being actually present, independently for each point. Then we want to analyze a certain functional, like the expected area of the minimum enclosing disk, the expected area of the smallest bounding box, or the probability that some point is contained in the convex hull.
In our scenario, we have to consider random insertion orders of the points and analyze the expected increase in the value of the characteristic function after the insertion of a fixed point p. Thus, we have to consider subsets of points constructed according to a different random process. In particular, whether other points precede p or not in the random order are not independent events. In our analysis, we condition on properties of the shape before adding the new point p. In the case of the minimum enclosing disk we use that each minimum enclosing disk is determined by at most 3 points, while in the case of the convex hull we condition separately on each single edge of the convex hull before the insertion of p. A straightforward application of this principle gives polynomial-time algorithms. To improve the running time we carry out this idea but split it into two computations: the expected value after the insertion of p, if the insertion of p changed the shape; the expected value before the insertion of p, if the insertion of p changed the shape. Finally, we use arrangements of lines and planes to speed up the computation by an additional linear factor.
For the AreaAnchoredRectangles, AreaBoundingBox, and AreaAnchoredBound-ingBox games we show that Shapley values can be computed in O(n 3/2 log n) time. In the special case where the points form a chain (increasing or decreasing y-coordinate for increasing x-coordinate), the Shapley values of those games can be computed in O(n log 2 n) time. We refer to these games as axis-parallel games.
It is relative easy to compute the Shapley values for these axis-parallel games in quadratic time using arrangements of rectangles and the linearity of Shapley values. We will discuss this as an intermediary step towards our solution. However, it is not obvious how to get subquadratic time. Besides using the linearity of Shapley values, a key ingredient in our algorithms is using convolutions to evaluate at multiple points some special rational functions that keep track of the ratio of permutations with a certain property. The use of computer algebra in computational geometry is very recent, and there are very few results [3,4,18] using such techniques in geometric problems. (In contrast, there is a quickly growing body of works in discrete geometry using real algebra.) At the basement of our algorithms, we need to count the number of permutations with certain properties. For this we use some simple combinatorial counting. In summary, our results combine fundamental concepts from several different areas and, motivated by classical concepts of game theory, we introduce new problems related to stochastic computational geometry and provide efficient algorithms for them.
Related work. The chapter by Deng and Fang [6] gives a summary of computational aspects of coalitional games. The book by Nisan et al. [20] provides a general overview of the interactions between game theory and algorithms.
Puerto, Tamir and Perea [23] study the the Minimum Radius Location Game in general metric spaces. When specialized to the Euclidean plane, this is equivalent to the PerimeterEn-closingDisk. The paper also considers the L 1 -metric, which is in some way related to the axis-parallel problems we consider here. However, the focus of their research is on understanding the core of the game, and do not discuss the computation of Shapley values. In particular, they show that the Minimum Radius Location Game in the Euclidean plane has nonempty core. Puerto, Tamir and Perea [24] also discuss the Minimum Diameter Location Game, which can be defined for arbitrary metric spaces, but then focus their discussion on graphs.
Faigle et al. [10] consider the TSP coalitional game in general metric spaces and specialize some results to the Euclidean plane. They also analyze an approximate allocation of the costs.
The computation of Shapley values have been considered for several games on graphs. The aforementioned Airport problem can be considered a shortest spanning-path game in a (graphtheoretic) path. Megiddo [17] extended this to trees, while Deng and Papadimitriou [7] discuss a game on arbitrary graphs defined by induced subgraphs. There is a very large body of follow up works for graphs, but we could not trace other works considering the computation of Shapley values for games defined through planar objects, despite being very natural.
Assumptions. We will assume general position in the following sense: no two points have the same x or y coordinate, no three points are collinear, and no four points are cocircular. In particular, the points are all different. The actual assumptions depend on the game under consideration. It is simple to consider the general case, but it makes the notation more tedious.
We assume a unit-cost real-RAM model of computation. In a model of computation that accounts for bit complexity, time bounds may increase by polynomial factors (even if the input numbers are integers, the outputs may be rationals with large numerators and denominators).
Organization. We start with preliminaries in Section 2, where we set the notation, define Shapley values, explain some basic consequences of the Airport game, discuss our needs of algebraic computations, and count permutations with some properties. The section is long because of our use of tools from different areas.
Then we analyze different games. The AreaConvexHull and AreaEnclosingDisk games are considered in Section 3 and Section 4, respectively. In Section 5 we discuss the AreaAnchoredRectangles game, while in Section 6 we discuss the AreaBoundingBox and AreaAnchoredBoundingBox games. To understand Section 6 one has to go through Section 5 before. The remaining sections (and games) have no dependencies between them and can be read in any order. We conclude in Section 7 with some discussion.

Preliminaries
Here we provide notation and background used through the paper.
Geometry. In most cases we will consider points in the Euclidean plane R 2 . The origin is denoted by o. For a point p ∈ R 2 , let R p be the axis-parallel rectangle with one corner at the origin o and the opposite corner at p. As already mentioned in the introduction, for each Q ⊂ R 2 we use bb(Q) for the (minimum) axis-parallel bounding box of Q, med(Q) for a smallest (actually, the smallest) disk that contains Q, and CH(Q) for the convex hull of Q.
We try to use the word rectangle when one corner is defined by an input point, while the word box is used for more general cases. All rectangles and boxes in this paper are axis-parallel, and we drop the adjective axis-parallel when referring to them. An anchored box is a box with one corner at the origin.
For a point p ∈ R 2 we denote by x(p) and y(p) its x-and y-coordinate, respectively. We use x(Q) for {x(q) | q ∈ Q}, and similarly y(Q). A set P of points is a decreasing chain, if x(p) < x(q) implies y(p) > y(q) for all p, q ∈ P . A set P of points is an increasing chain if x(p) < x(q) implies y(p) < y(q) for all p, q ∈ P .

Shapley values.
We provide a condensed overview suitable for understanding our work. We refer to some textbooks in Game Theory ( [11,Chapter IV.3], [19,Section 9.4], [21,Section 14.4]) for a comprehensive treatment that includes a discussion of the properties of Shapley values.
Consider a coalitional game (P, v). For us, a player is usually denoted by p ∈ P . We denote by n the number of players and by [n] the set of integers {1, . . . , n}.
A permutation of P is a bijective map π : P → [n]. Let Π(P ) be the set of permutations of P . Each permutation π ∈ Π(P ) defines an ordering in P , where π(p) is the position of p in that order. We will heavily use this interpretation of permutations as defining an order in P . For each element p ∈ P and each permutation π ∈ Π(P ), let P (π, p) be the elements of P before p in the order defined by π. Thus P (π, p) = {q ∈ P | π(q) ≤ π(p)}.
We can visualize P (π, p) as adding the elements of P one by one, following the order defined by π, until we insert p. The increment in v(·) when adding player p is The Shapley value of player p ∈ P in the game (P, v) is It is straightforward to note that where π is picked uniformly at random from Π(P ). This means that Shapley values describe the expected increase in v(·) when we insert the player p in a random permutation.
Since several permutations π define the same subset P (π, p), the Shapley value of p is often rewritten as We will use a similar principle in our algorithms: find groups of permutations that contribute the same to the sums. For non-axis-parallel games we will also use other equivalent expressions of the Shapley value that follow from rewriting the sum differently. It is easy to see that Shapley values are linear in the characteristic functions. Indeed, since for any two characteristic functions v 1 , v 2 : 2 P → R and for each λ 1 , λ 2 ∈ R we have In fact, (a weakened version of) linearity is one of the properties considered when defining Shapley values axiomatically.
It is easy to obtain the Shapley values when the characteristic function v is constant over all nonempty subsets of P . In this case φ(p, v) = v(P )/n.
We say that two games (P, v) and (P , v ), where P and P are point sets in Euclidean space, are isometrically equivalent if there is some isometry that transforms one into the other. That is, there is some isometry ρ : R 2 → R 2 such that P = ρ(P ) and v = v • ρ. Finding Shapley values for (P, v) or (P , v ) is equivalent because φ(ρ(p), v ) = φ(p, v) for each p ∈ P .

Airport and related games.
We will consider the following games, some of them 1-dimensional games.
Airport game: The set of players P is on the positive side of the real line and the characteristic function is v air (Q) = max(Q) for each nonempty Q ⊂ P .
LengthInterval game: The set of players P is on the real line and the characteristic function is v li (Q) = max(Q) − min(Q) for each nonempty Q ⊂ P .
AreaBand game: The set of players P is on the plane and the characteristic function is v ab (Q) = (max(y(P )) − min(y(P )) · v li (x(Q) for each nonempty Q ⊂ P .
PerimeterBoundingBox game: The set of players P is on the plane and the characteristic function is v pbb (Q) = per(bb(Q)) for each nonempty Q ⊂ P .
PerimeterAnchoredBoundingBox game: The set of players P is on the plane and the characteristic function is v pabb (Q) = per(bb(Q ∪ {o})) for each nonempty Q ⊂ P . Proof. The Airport game was solved by Littlechild and Owen [15], as follows. If P is given by 0 < x 1 < · · · < x n , then This shows the result for the Airport game because it requires sorting.
The LengthInterval game can be reduced to the Airport problem using simple inclusionexclusion on the line. Consider the values α = min(P ) and β = max(P ) and define the characteristic functions They define games that are isometrically equivalent to Airport games or a constant game, and we have v li (Q) = v 1 (Q) + v 2 (Q) + v 3 (Q) for all nonempty Q ⊂ P . The result follows from the linearity of the Shapley values. 6 For the AreaBand game we just notice that v ab (Q) = v li (Q) · (max(y(P )) − min(y(P ))) and use again the linearity of Shapley values. (The value max(y(P )) − min(y(P )) is constant.) For the PerimeterBoundingBox game we note that v pbb (Q) = 2 · (max(x(P )) − min(x(P ))) + 2 · (max(y(P )) − min(y(P ))), which means that we need to compute the Shapley values for two LengthInterval games. For the PerimeterAnchoredBoundingBox we use that v pabb (Q) = 2 · (max(x(P ), 0) − min(x(P ), 0)) + 2 · (max(y(P ), 0) − min(y(P ), 0)), which is a linear combination of a few Airport games.
The perimeter for the anchored bounding box is the same as for the union of anchored rectangles. So this solves the perimeter games for all axis-parallel objects considered in the paper. Remark. We can now point out the convenience of using the real-RAM model of computation. In the Airport game with x i = i for each i ∈ [n], we have φ(x n , v air ) = i∈[n] 1 i , the n-th harmonic number. Thus, even for simple games it would become a challenge to carry precise bounds on the number of bits.

Algebraic computations.
For the axis-parallel problems we will be using the fast Fourier transform to compute convolutions. Assume we are given values a 0 , . . . , a n and b 0 , . . . , b n , and define a i and b i equal zero for all other indices. Using the fast Fourier transform we can compute the convolutions c k = i+j=k a i b j for all integer k in O(n log n) operations. (The value is nonzero for at most 2n + 1 indices.) Our use of this is encoded in the following lemma, which provides multipoint evaluation for a special type of rational functions. Note that using the more generic approach of Aronov, Katz and Moroz [4,18] for multipoint rational functions gives a slightly worse running time.
Lemma 2. Let b 0 , . . . , b n , ∆ be real numbers and consider the rational function Given an integer > −∆, possibly negative, and a positive integer m, we can evaluate R(x) at all the integer values x = , + 1, . . . , + m in O((n + m) log(n + m)) time.
Proof. Set a i = 1/(∆ + + m − i) for all i ∈ {0, . . . , m}. Note that the assumption ∆ + > 0 implies that the denominator is always positive. All the other values of a * and b * are set to 0. Define the convolutions c k = i a i b k−j for all k ∈ {0, . . . , m} and compute them using the fast Fourier transform in O((n + m) log(n + m)) time. Then, for each integer x we have Thus, from computing the convolution of a * and b * , we get the values R( ), . . . , R( + m)

Permutations.
We will have to count permutations with certain properties. The following lemmas encode this. Proof. Permute the α elements of A arbitrarily and put them after x. Similarly, permute the β elements of B arbitrarily and put them before x. Finally, insert the remaining |N \(A∪B∪{x})| = n − α − β − 1 elements arbitrarily, one by one. Each permutation that we want to count is constructed exactly once with the procedure, and thus there are Lemma 4. Let N be a set with n elements. Let A, B and C be pairwise disjoint subsets of N with α, β and γ elements, respectively.
• Consider a fixed element a ∈ A. There are n! · permutations of N such that: (i) a is the first element of A and (ii) a is before all elements of B or before all elements of C. (Thus, the whole B or the whole C is behind a.) Proof. We start showing the first item. Fix an element a ∈ A. For a set X, disjoint from A, let Π X be the set of permutations that have a before all elements of (A \ {a}) ∪ X. We can count these permutations as follows. First, permute the elements of (A \ {a}) ∪ X and fix their order. Then insert a at the front. Finally, insert the elements of N \ (A ∪ X) anywhere, one by one. All permutations of Π X are produced with this procedure exactly once. Therefore, there are This finishes the proof of the first item. Now we prove the second item. Fix an element b ∈ B. We construct the desired permutations as follows. Select one element c ∈ C and place c before b. Now insert all the elements of A ∪ B \ {b} after b in arbitrary order. Then, insert all the elements of C \ {c} after c, one by one. Finally, insert the remaining elements in any place, one by one. Each of permutations under consideration is created exactly once by this procedure. Therefore, there are permutations, which is exactly n! · γ (α+β)(α+β+γ) . Then we use that

Convex hull
In this section we consider the area and the perimeter of the convex hull of the points. We will focus on the area, thus the characteristic function v ch , and at the end we explain the small changes needed to handle the perimeter. Consider a fixed set P of points in the plane. For simplicity we assume that no three points are collinear. : Left: triangulating the difference between CH(P (π, p)) and CH(P (π, p) \ {p}). Right: in order for pqq to be in T (π, p), q and q must appear before p, which in turn must appear before all other points in the halfplane H(q, q ).
Proof. For each q, q ∈ P (q = q ), let H(q, q ) be the open halfplane containing all points to the left of the directed line from q to q . Define level(q, q ) to be the number of points in P ∩ H(q, q ). We can decompose the difference between CH (P (π, p)) and CH (P (π, p) \ {p}) into a set T (π, p) of triangles (see Figure 2 (left)), where q and q appear before p in π, and no points before p in π lie in H(q, q )}.
In other words, pqq ∈ T (π, p) if and only if p ∈ H(q, q ), and among the level(q, q ) + 2 points in H(q, q ) ∪ {q, q }, the two earliest points are q and q , and the third earliest point is p. See Figure 2 (right). (Note that if CH (P (π, p)) = CH (P (π, p) \ {p}), then T (π, p) is empty.) For fixed p, q, q ∈ P with p ∈ H(q, q ), Lemma 3 with x = p, B = {q, q } and A = H(q, q ) ∩ P \ {p} tells that the probability that pqq ∈ T (π, p) with respect to a random permutation π is exactly .
It follows that the Shapley value of p is From the formula, we can immediately compute φ(p, v ch ) for any given p ∈ P in O(n 2 ) time, if all ρ(q, q ) values have been precomputed. Each value ρ(q, q ) can be computed from level(q, q ) using O(1) arithmetic operations. Thus, precomputing ρ(q, q ) requires precomputing level(q, q ) for all O(n 2 ) pairs q, q . In the dual, this corresponds to computing the levels of all O(n 2 ) vertices in an arrangement of n lines. The arrangement of n lines can be constructed in O(n 2 ) time [5,Chapter 8], and the levels of all vertices can be subsequently generated by traversing the arrangement in O(1) time per vertex.
Naively applying Lemma 5 to all points p ∈ P gives O(n 3 ) total time. We can speed up the algorithm by a factor of n: Proof. Let p = (x, y) ∈ P . Observe that for fixed q, q ∈ P (q = q ), if p ∈ H(q, q ), then area( pqq ) is a linear function in x and y and can thus be written as a(q, q )x+b(q, q )y +c(q, q ). Let A(q, q ) = a(q, q ) · ρ(q, q ), B(q, q ) = b(q, q ) · ρ(q, q ), and C(q, q ) = c(q, q ) · ρ(q, q ). (As noted earlier, we can precompute all the ρ(q, q ) values in O(n 2 ) time from the dual arrangement of lines.) By (1), We describe how to compute A(p), B(p), and C(p) for all p ∈ P in O(n 2 ) total time. Afterwards, we can compute φ(p, v ch ) for all p ∈ P in O(n) additional time.
The problem can be reduced to 3 instances of the following: Given a set P of n points in the plane, and given O(n 2 ) lines each through 2 points of P and each assigned a weight, compute for all p ∈ P the sum of the weights of all lines below p (or similarly all lines above p).
In the dual, the problem becomes: Given a set L of n lines in the plane, and given O(n 2 ) vertices in the arrangement, each assigned a weight, compute for all lines ∈ L the sum S( ) of the weights of all vertices below .
To solve this problem, we could use known data structures for halfplane range searching, but a direct solution is simpler. First construct the arrangement in O(n 2 ) time. Given , ∈ L, define S( , ) to be the sum of the weights of all vertices on the line that are below . For a fixed ∈ L, we can precompute S( , ) for all ∈ L \ { } in O(n) time, since these values correspond to prefix or suffix sums over the sequence of weights of the O(n) vertices on the line . The total time for all lines ∈ L is O(n 2 ). Afterwards, for each ∈ L, we can compute S( ) in O(n) time by summing S( , ) over all ∈ L \ { } and dividing by 2 (since each vertex is counted twice). The total time for all ∈ L is O(n 2 ). Proof. We adapt the algorithm for area to handle the perimeter. Let E(π, p) be the set of directed edges of CH (P (π, p)) that are not in CH (P (π, p) \ {p}), where edges are oriented in clockwise order. (If CH (P (π, p) = CH (P (π, p) \ {p}), then E(π, p) is empty; otherwise, it has exactly two edges because of general position.) Then (q, p) ∈ E(π, p) if and only if q appears before p in π and no points before p in π lie in H(q, p); in other words, among the level(q, p) + 2 points in H(q, p) ∪ {p, q}, the earliest point is q and the second earliest point is p. We can use Lemma 3 to bound the probability of this event. Thus, for fixed p, q ∈ P , the probability that (q, p) ∈ E(π, p) (with respect to a random permutation π) is exactly ρ (q, p) = 1 (level(q, p) + 2)(level(q, p) + 1) . Similarly, the probability that (p, q) ∈ E(π, p) is exactly ρ (p, q). It follows that the expected total length of the edges in E(π, p) is On the other hand, a modification of the proof of Lemma 5 shows that the expected total length of the edges in CH (P (π, p) \ {p}) that are not in CH (P (π, p)) is . We can compute φ − (p) for all p ∈ P as in the proof in Theorem 6, in O(n 2 ) total time (in fact, the algorithm is a little simpler, since linear functions are not required).

Smallest enclosing disk
In this section we consider the area and the perimeter of the smallest enclosing disk of the points. Recall that v ed is the characteristic function. For simplicity, we assume general position in the following way: no four points are cocircular and circles through three input points do not have a diameter defined by two input points. We use P for the set of points.
First we explain how to compute the Shapley values for the area of the minimum enclosing disk.
Proof. For each subset of points Q ⊂ P , let X(Q) be the subset of points of Q that lie on the boundary of med(Q). We now recall some well known properties; see for example [27] or [5,Section 4.7]. The disk med(Q) is unique and med(X(Q)) = med(Q). If a point p is outside med(Q), then p is on the boundary of med(Q ∪ {p}), that is, p ∈ X(Q ∪ {p}). Because of our assumption on general position, X(Q) has at most three points. When X(Q) has two points, then they define a diameter of med(Q). We have |X(Q)| ≤ 1 only when |Q| ≤ 1. See Figure 3.  Figure 4: Left: in order for the shown disk to be med(P (π, p) \ {p}), the points q, q , q must appear before p, which in turn must appear before all other points outside the disk. Right: in order for the shown disk to be med(P (π, p)), the points q, q must appear before p, which in turn must appear before all points outside the disk.
earliest point is p. See Figure 4 (left). Lemma 3 implies that, for a fixed B and a fixed p, the probability that X(P (π, p) \ {p}) = B (with respect to a random permutation π) is exactly .
Let I(π, p) be the indicator variable that is 1 if p ∈ med(P (π, p) \ {p}), and 0 otherwise. It follows that the expected value of area(med(P (π, p) \ {p})) · I(π, p) is On the other hand, for a basis B and a point p ∈ B, we have X(P (π, p)) = B if and only if all points of B \ {p} appear before p in π, and no points before p in π lie in out(B). In other words, among the level(B) + |B| points in out(B) ∪ B, the |B| − 1 earliest points are the points of B \ {p} and the next earliest point is p. See Figure 4 (right). We can bound the probability of this event using Lemma 3. Thus, for a fixed B and p ∈ B, the probability that X(P (π, p) \ {p}) = B (with respect to a random permutation π) is exactly ρ (B) = (|B| − 1)! (level(B) + |B|) · · · (level(B) + 1) .
It follows that the expected value of area(med(P (π, p)) · I(π, p) is By linearity of expectation, the Shapley value of p is φ(v ed , p) = φ + (p) − φ − (p). From the formulas (2)  Proof. We already know how to compute φ + (p) in O(n 2 ) time for each p ∈ P , and thus in O(n 3 ) total time. It suffices to focus on computing φ − (p). In the formula (2), terms for bases of size 2 can be handled in O(n 2 ) time for each p; so it suffices to focus on bases of size 3. The problem can be formulated as follows: Given a set P of n points in the plane, and given O(n 3 ) disks each with 3 points of P on the boundary and each assigned a weight, compute for all p ∈ P the sum of the weights of all disks not containing p.
By the standard lifting transformation, the problem reduces to the following: Given a set H of n planes in 3 dimensions, and given O(n 3 ) vertices in the arrangement, each assigned a weight, compute for all h ∈ H the sum S(h) of the weights of all vertices below h.
To solve this problem, we could use known data structures for halfspace range searching, but an approach as in the proof of Theorem 6 is simpler.

Union of anchored rectangles
In this section we consider the AreaAnchoredRectangles game defined by the characteristic function v ar . It is easy to see that one can focus on the special case where all the points are in a quadrant; see Figure 5. Our discussion will focus on the case where the points of P are on the positive quadrant of the plane.
Consider a fixed set P of points in the positive quadrant. In the notation we will drop the dependency on P . For simplicity, we assume general position: no two points have the same x-or y-coordinate. We first introduce some notation that will be used in this section and in Section 6. q NE(q) NW(q) Figure 6: Notation for axis-parallel problems. Left: the quadrants to define NW(q), NE(q) and SE(q). Right: cells of A.

Notation for axis-parallel problems
For each point q of the plane, we use the "cardinal directions" to define subsets of points in quadrants with apex at q: We use lowercase to denote their cardinality: nw(q) = | NW(q)|, ne(q) = | NE(q)| and se(q) = | SE(q)|. See Figure 6, left. Let x 1 < . . . < x n denote the x-coordinates of the points of P , and let y 1 < . . . < y n be their y-coordinates. We also set x 0 = 0 and y 0 = 0. For each i, j ∈ [n] we use w i = x i − x i−1 (for width) and h j = y j − y j−1 (for height).
Let L be the set of horizontal and vertical lines that contain some point of P . We add to L both axis. Thus, L has 2n + 2 lines. The lines in L break the plane into 2-dimensional cells (rectangles), usually called the arrangement and denoted by A = A(L). More precisely, a (2-dimensional) cell c of A is a maximal connected component in the plane after the removal of the points on lines in L. Formally, the (2-dimensional) cells are open sets whose closure is a rectangle, possibly unbounded in some direction. We are only interested in the bounded cells, and with a slight abuse of notation, we use A for the set of bounded cells. We denote by c i,j the cell between the vertical lines x = x i−1 and x = x i and the horizontal lines y = y j−1 and y = y j . Note that c i,j is the interior of a rectangle with width w i and height h j . See Figure 6, right.
Since NE(q) is constant over each 2-dimensional cell c of A, we can define NE(c), for each cell c ∈ A. The same holds for NW(c) and SE(c) and their cardinalities, ne(c), nw(c) and se(q). See Figure 7.
A block is a set of cells B = B(i 0 , i 1 , j 0 , j 1 ) = {c i,j | i 0 ≤ i ≤ i 1 , j 0 ≤ j ≤ j 1 } for some indices i 0 , i 1 , j 0 , j 1 , with 1 ≤ i 0 ≤ i 1 ≤ n and 1 ≤ j 0 ≤ j 1 ≤ n. The number of columns and rows in B is i 1 − i 0 + 1 + j 1 − j 0 + 1 = O(i 1 − i 0 + j 1 − j 0 ). A block B is empty if no point of P is on the boundary of at least three cells of B. Equivalently, B is empty if no point of P is in the interior of the union of the closure of the cells in B. See Figure 8 for an example.
We will be using maximal rows and columns within a block B to compute some partial information. Thus, for each block B and each index i, we define the vertical slab V (i, B) = {c i,j | 1 ≤ j ≤ n, c i,j ∈ B}. Similarly, for each block B and each index j, we define the horizontal slab H(j, B) = {c i,j | 1 ≤ i ≤ n, c i,j ∈ B}. Such slabs are meaningful only for indices within the range that defines the slab. We call them the slabs within B.

Interpreting Shapley values geometrically
First, we reduce the problem of computing Shapley values to a purely geometric problem. See Figure 7 for the relevant counters considered in the following result. For each subset C of cells of A, we define (We will only consider sets C of cells with ne(c) > 0 for all c ∈ C.) Note that we want to compute σ(·) for the sets of cells contained in the rectangles R p for all p ∈ P . Using standard tools in computational geometry we can compute the values φ(p, v ar ) for all p ∈ P in near-quadratic time. Here is an overview of the approach. An explicit computation of A takes quadratic time in the worst case. Using standard data structures for orthogonal range searching, see [28] or [5,Chapter 5], we can then compute ne(c) for each cell c ∈ A. Finally, replacing each cell c by a point q c ∈ c with weight w c = area(c)/ ne(c), we can reduce the problem of computing φ(p, v ar ) to the problem of computing qc∈Rp w c , which is again an orthogonal range query. An alternative is to use dynamic programming across the cells of A to compute first ne(c) and then partial sums of the weights w c .
Our objective in the following sections is to improve this result using the correlation between adjacent cells. It could seem at first glance that segment trees [5, Section 10.3] may be useful. We did not see how to work out the details of this; the problem is that the weights are inversely proportional to ne(c).

Handling empty blocks
In the following we assume that we have preprocessed P in O(n log n) time such that ne(q) can be computed in O(log n) time for each point q given at query time [28]. This is a standard range counting for orthogonal ranges and the preprocessing has to be done only once.
When a block B is empty, then we can use multipoint evaluation to obtain the partial sums σ(·) for each vertical and horizontal slab of the block.
Lemma 11. Let B be an empty block with k columns and rows. We can compute in O(k log n) time the values σ(C) for all slabs C within B.
We look into the first vertical slab V (i 0 , B) and make groups of cells depending on their value ne(·). More precisely, for each we define J( ) = {j | j 0 ≤ j ≤ j 1 , ne(c i 0 ,j ) = }. Let 0 and 1 be the minimum and the maximum such that J( ) = 0, respectively.
We set up the following rational function with variable x: Setting t = − 0 , b t = j∈J( 0 +t) h j and ∆ = 0 , we have Thus, this is a rational function of the shape considered in Lemma 2 with 1 − 0 ≤ j 1 − j 0 + 1 ≤ k terms. The coefficients can be computed in O(k log n) time because we only need the values h j and ne(c i 0 ,j ) for each j. These latter values ne(·) are obtained from range counting queries. Note that A similar statement holds for all the other vertical slabs within B. We make the statement precise in the following. Consider two consecutive vertical slabs V (i, B) and V (i + 1, B) within the block B. Because the block B is empty, the difference ne(c i+1,j )−ne(c i,j ) is independent of j. See Figure 9. It follows that, for each index i with i 0 ≤ i ≤ i 1 , there is an integer δ i such that ne(c i,j ) = ne(c i 0 ,j ) + δ i for all j with j 0 ≤ j ≤ j 1 . Moreover, for each i with i 0 ≤ i ≤ i 1 and each with 0 ≤ ≤ 1 , the value of ne(c i,j ) is constant over all j ∈ J( ). Therefore, for each j ∈ J( ) we have ne(c i,j ) = + δ i .
Each value δ i can be obtained using that δ i = ne(c i,j 0 ) − ne(c i 0 ,j 0 ). This means that the values δ i 0 , . . . , δ i 1 can be obtained in O(k log n) time. Now we note that, for each index i with i 0 ≤ i ≤ i 1 , we have i, B)).
We use Lemma 2 to evaluate the After this, we get each value σ(V (i, B)) = w i · R(δ i ) in constant time.

Chains
In this section we consider the case where the points are a chain. As discussed before, it is enough to consider that P is in the positive quadrant. After sorting, we can assume without loss of generality that the points of P are indexed so that 0 < x(p 1 ) < · · · < x(p n ). We start with the easier case: increasing chains. The problem is actually an Airport game in disguise.
Thus, z i is the area of the region R p i \ R p i−1 divided by ne(c), for some c ⊂ R p i \ R p i−1 . See Figure 10. Because of Lemma 10, we have for each i ∈ [n] Therefore φ(p 1 , v ar ) = z 1 and, for each i ∈ [n] \ {1}, we have φ(p i , v ar ) = z i + φ(p i−1 , v ar ). The result follows.
(Actually this game is just the 1-dimensional Airport game if each point p i is represented on the real line with the point with x-coordinate area(R p i ).) It remains the more interesting case, when the chain is decreasing. See Figure 11 for an example. It is straightforward to see that, if i + j > n + 1, then ne(c i,j ) = 0, and if i + j ≤ n + 1, we have ne(c i,j ) = n + 2 − i − j. So in this case we do not really need data structures to obtain ne(c i,j ) efficiently. (The proof of Lemma 11 can be slightly simplified for this case because J( ) contains a single element due to the special structure of the values ne(c). See Figure 12.) We use a divide-and-conquer paradigm considering certain empty blocks defined by two indices and r, where < r. Since the indexing of rows is not the most convenient in this case, it is better to introduce the notation B ,r for the block B( + 1, m, n − r + 2, n − m + 1), where  X(a, B ,r ). In this example, 3n/16 < a < n/4.
For each point p a of P , we can compute the set I(a) and then the sum ( ,r)∈I(a) σ(X(a, B ,r )) = φ(p a , v ar ) in O(log n) time. Over all points of P we spend O(n log n) time in this last computation.
When the point set is a chain over different quadrants, we can reduce it to a few problems over the positive quadrant and obtain the following.

General point sets
We consider now the general case; thus the points do not form a chain. See Figure 7 for an example. Like before, we restrict the discussion to the case where P is in the positive quadrant. In this case, ne(c i,j ) is not a simple expression of i, j anymore, As mentioned in Section 5.3, we preprocess P in O(n log n) time such that ne(q) can be computed in O(log n) time [28].
In this scenario we consider horizontal bands. A horizontal band B is the block between two horizontal lines. Thus, B = {c i,j | j 0 ≤ j ≤ j 1 } for some indices 1 ≤ j 0 ≤ j 1 ≤ n. See Figure 16. We keep using the notation introduced for blocks. Thus, for each i ∈ [n], let V (i, B) be the vertical slab with the cells c i,j ∈ B. Let P B be the points of P that are the top-right corner of some cell of B. We use k B = |P B |. Because of our assumption on general position, k B = j 1 −j 0 +1 and thus k B is precisely the number of horizontal slabs in the band B. Furthermore, for each point p ∈ P B we define the rectangle R(p, B) as the cells of B to the left and bottom of p.
Lemma 15. For a band B with k B rows we can compute in O(((k B ) 2 +n) log n) time σ (V (i, B)), for all i ∈ [n], and σ (R(p, B)), for all p ∈ P B . Proof. Assume that B is defined by the row indices j 0 ≤ j 1 . Let q 1 , . . . , q k B be the points of P B sorted by increasing x coordinate. Take also q 0 as a point on the y-axis and q k B +1 as a point on the right boundary. For each t denote by m t the number of vertical slabs between the vertical lines through q t−1 and q t . See Figure 16. We divide the band B into blocks B 1 , B 2 , . . . using the vertical lines though the points of P B . We get k B + 1 blocks (or k B if the rightmost point of P belongs to P B ), and each of them is empty by construction. Since the block B t has m t vertical slabs and k B horizontal slabs we can compute σ(V (·, B t )) and σ(H(·, B t )) for all the slabs within B t in O((k B + m t ) log n) time using Lemma 11. Using that the O(k B ) blocks B 1 , B 2 , . . . are pairwise disjoint, which means that t m t = n, we conclude that in O(((k B ) 2 + n) log n) time we can compute all the values σ(V (·, B t )) and σ(H(·, B t )) for all slabs within all blocks B t .
Since for each vertical slab V (i, B) of B there is one block B t(i) that covers it, we then have σ(V (i, B)) = σ(V (i, B t(i) )) for such index t(i). This means that we have already computed the values σ(V (i, B)) for all i. Now we explain the computation of σ(R(p, B)) for all p ∈ P B . Note that within each block B t we have computed σ(H(j, B t )) for all j. With this information we can compute the values σ(R(p, B)). Namely, for each block B t and each j with j 0 ≤ j ≤ j 1 , we define Using that we compute all the values σ(H ≤ (j, B t )) in O((k B ) 2 ) time. Since R(p, B) = H ≤ (j(t), B t(p )) for some indices j(p) and t(p) that we can easily obtain, the lemma follows. Proof. As discussed earlier, it is enough to consider the case when P is in the positive quadrant because the other quadrants can be transformed to this one. The cells inside a rectangle R p now correspond to the disjoint union where B p is the band that contains p (perhaps on its top boundary) and the index i p is such that the cell c ip,j has p on its top right corner. See Figure 17.

Area of the bounding box
In this section we are interested in the AreaBoundingBox game defined by the characteristic function v bb . The structure of this section is similar to the structure of Section 5. In particular, we keep assuming that all points have different coordinates and we keep using the notation introduced in Section 5.1. First we note that it is enough to consider the AreaAnchoredBoundingBox problem and assume that the points are in one quadrant. Recall that in this problem the origin o has to be included in the bounding box. Thus, it uses the characteristic function v abb (Q) = area(bb(Q ∪ {o})). Proof. We use inclusion exclusion, as indicated in Figure 18. There are characteristic functions v 1 , . . . , v 9 such that, for each Q ⊆ P , we have v bb (Q) = 9 i=1 v i (Q). Moreover, each characteristic function v i is either a constant value game (first summand in Figure 18), isometrically equivalent to an AreaBand game (last four terms in Figure 18), or isometrically equivalent to an AreaAnchoredBoundingBox problem with all points in one quadrant (the remaining four summands in Figure 18). Finally, note that if the points form a chain, they also form a chain in each of the cases possibly exchanging the increasing/decreasing character. The following lemma is summarized in Figure 19. permutations that fulfill this criteria. This means that, for each p ∈ NE(c) we have

Interpreting Shapley values geometrically
Consider now a point p ∈ NW(c). For a permutation π ∈ Π(P ), we have ∆(v c , π, p) = area(c) if and only if p is the first point of NW(c) in the permutation π, all the points of NE(c) are after p in the permutation π, and at least one point of SE(c) is before p in the permutation π. permutations that fulfill this criteria. This means that, for each p ∈ NW(c) we have φ(p, v c ) = area(c) · ψ NW (c).
The case for a point p ∈ SE(c) is symmetric to the case p ∈ NW(c), where the roles of se(c) and nw(c) are exchanged. Therefore we conclude As a sanity check it good to check that ne(c) · ψ NE (c) + nw(c) · ψ NW (c) + se(c) · ψ SE (c) = 1. This is the case.
Noting that we have, For each subset of cells C of A, we define Like in the case of anchored boxes, our objective here is to compute these values for several vertical and horizontal slabs.

Handling empty blocks
In the following we assume that we have preprocessed P in O(n log n) time such that ne(q), nw(q) and se(q) can be computed in O(log n) time for each point q given at query time [28]. We use again multipoint evaluation to compute the values σ * (·) for each vertical and horizontal slab of an empty block and for each * ∈ {NE, NW, SE}. However, the treatment has to be a bit more careful now because we have to deal with different rational functions.
Lemma 19. Let B be an empty block with k columns and rows. We can compute in O(k log n) time the values σ * (C) for all slabs C within B and each * ∈ {NE, NW, SE}.
Proof. The proof is very similar to the proof of Lemma 11, but some adaptation has to be made. Assume that B is the block B(i 0 , i 1 , j 0 , j 1 ). We explain how to compute the values σ NE (V (i, B)) for all i 0 ≤ i ≤ i 1 . The technique to compute σ NE (H(j 0 , B)), . . . , σ NE (H(j 1 , B)) and for NW and SE instead of NE is similar.
We set up the following fractional constant and rational functions with variable x: Each of them is the sum of at most k fractions. As discussed in the proof of Lemma 11, R 2 (x) and R 3 (x) can be rewritten so that we can use Lemma 2 to evaluate them at consecutive points. The coefficients can be computed in O(k log n) time.
Consider two consecutive vertical slabs V (i, B) and V (i + 1, B) within the block B. Figure 20 shows the changes in the counters. Because the block B is empty we have the following properties: • The sum ne(c i,j ) + nw(c i,j ) is independent of i.  • The difference ne(c i+1,j ) + se(c i+1,j ) − ne(c i,j ) + se(c i,j ) is always −1. Therefore we have for all relevant i and j. In particular, for each j ∈ J( ) we have ne(c i,j )+se(c i,j ) = −(i−i 0 ).
• The difference ne(c i+1,j ) + nw(c i+1,j ) + se(c i+1,j ) − ne(c i,j ) + ne(c i,j ) + se(c i,j ) depends on whether the point with x(p) = x i is above or below B. Thus, this difference is independent of the row j. This means that there exist, for each index i, a value δ i such that ne(c i,j ) + nw(c i,j ) + se(c i,j ) = ne(c i 0 ,j ) + nw(c i 0 ,j ) + se(c i 0 ,j ) + δ i for all relevant j. In particular, for each j ∈ J ( ) we have ne(c i,j )+nw(c i,j )+se(c i,j ) = +δ i . Each single value δ i can be computed as δ i = se(c i,j ) − se(c i 0 ,j ) in O(log n) time using range counting queries.
Note that for each relevant i we have (V (i, B)).
Thus, we need to evaluate R 2 (x) at the values x = 0, −1, . . . , i 0 − i 1 and we have to evaluate R 3 (x) at the values δ i for i = i 0 , . . . , i 1 . According to Lemma 2, this can be done in O(k log n) time. After this, we get each value σ NE (V (i, B)) in constant time.

Algorithms for bounding box
The same techniques that were used in Sections 5.4 and 5.5 work in this case. We compute partial sums σ * (·) for several vertical bands and horizontal bands. The only difference is that we have to cover the whole bounding box bb(P ∪ {o}), instead of only the portion that was dominated by some point. In the case of an increasing chain in the positive quadrant, the anchored bounding box and the union of anchored rectangles is actually the same object. In the case of a decreasing chain, we have to work above and below the chain to cover the whole bounding box. See Figure 21. This is twice as much work, so it does not affect the asymptotic running time. For arbitrary point sets, we again use the bands with roughly √ n horizontal slabs, and break each slab into empty boxes, as it was done in Lemma 15. For the partial sums that we consider (like σ(H ≤ (·)) and σ(V ≤ (·))) we also construct the symmetric versions σ * (H ≥ (·)) and σ * (V ≥ (·)). With this information we can recover the sums that we need in each of the three quadrants with an apex in p ∈ P ; recall Figure 19. In the case of a decreasing chain, we get an extra case that is handled by noting that η a = c∈A, pa∈SE(c) area(c) · ψ SE (c) is η a = η a−1 + n j=a+1 area(c a,j ) · ψ SE (c a,j ) + a−1 i=1 area(c i,n−a+1 ) · ψ SE (c i,n−a+1 ).
The last two terms are a vertical and a horizontal band around the cell c a,n−a+1 . The sums for ψ NW are handled similarly. We omit the details and summarize. Because of Lemma 17 we also obtain the following.

Conclusions
We have discussed the efficient computation of Shapley values, a classical topic in game theory, for coalitional games defined for points in the plane and characteristic functions given by the area of geometric objects. For axis-parallel problems we used computer algebra to get faster algorithms. For non-axis-parallel problems we provided efficient algorithms based on decomposing the sum into parts and grouping permutations that contribute equally to a part.
In game theory, quite often one considers coalitional games where some point, say the origin, has to be included in the solutions that are considered. For example, we could use the minimum enclosing disk that contains also the origin. We do this for the anchored versions of the games. This setting is meaningful in several scenarios, for example, when we split the costs to connect to something. Our results also hold in this setting through an easy adaptation.
The problems we consider here are a new type of stochastic problems in computational geometry. The relation to other problems in stochastic computational geometry is unclear. For example, computing the expected length of the minimum spanning tree (MST) in the plane for a stochastic point set is #P-hard [13]. Does this imply the same for the coalitional game based on the length of the Euclidean MST? Is there a two-way relation between Shapley values and the stochastic version for geometric problems? In particular, is there a FPRAS for computing the Shapley values of the game defined using the Euclidean MST? For the stochastic model this is shown by Kamousi et al. [13]. Note that the length of the Euclidean spanning tree is not monotone: adding points may reduce its length. Thus, some Shapley values could be potentially 0, which, at least intuitively, makes it harder to get approximations.
Finally, it would be worth to understand whether there is some relation between Shapley values in geometric settings and the concept of depth in point sets, in particular for the Tukey depth. For stochastic points, this relation has been explored and exploited by Agarwal et al. [1].
Our algorithms generalize to higher dimensions, increasing the degree of the polynomial describing the running time. Are there substantial improvements possible for higher dimensions? In dimension d, computing the volume of the bounding box can be done in O(dn) time, while the obvious algorithms to compute the Shapley values for the corresponding game require n Θ(d) time. Is this linear dependency on d in the degree of the polynomial actually needed, under some standard assumption in complexity theory?