See pages 1 of title_page.pdf

# Accelerated RRT* and its evaluation on Autonomous Parking

Jiri Vlasak12 , Michal Sojka2 and Zdeněk Hanzálek2

(1Faculty of Electrical Engineering, Czech Technical University in Prague
2Czech Institute of Informatics, Robotics and Cybernetics, Czech Technical
University in Prague
{jiri.vlasak.2, michal.sojka<EMAIL_ADDRESS>)

###### Abstract

Finding a collision-free path for autonomous parking is usually performed by
computing geometric equations, but the geometric approach may become unusable
under challenging situations where space is highly constrained. We propose an
algorithm based on Rapidly-Exploring Random Trees Star (RRT*), which works
even in highly constrained environments and improvements to RRT*-based
algorithm that accelerate computational time and decrease the final path cost.
Our improved RRT* algorithm found a path for parallel parking maneuver in
$95\text{\,}\%$ of cases in less than
$0.15\text{\,}\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}\mathrm{s}$.

Keywords: Autonomous parking, Rapidly-Exploring Random Trees, Reeds and Shepp
steering, Dijkstra optimization, Nearest neighbor heuristics.

##  1 INTRODUCTION

Modern cars are commonly equipped with parking assistants that can perform
parallel or perpendicular parking maneuvers. Parking is a relatively easy task
as the movement is slow and the car dynamics might be neglected. Usually,
geometric equations are used for planning these maneuvers. A geometric
approach has limitations when applied in unexpected environments or when more
than a simple parking maneuver has to be planned. In this paper we address the
cases, when more advanced planners need to be used, and one of the problems
experienced by those complex planners is their computational complexity.

For this paper, we define the parking problem as finding a collision-free path
from an initial car position (i.e., $x$, $y$, and $heading$) to the goal
position under the presence of an arbitrary number of known static obstacles.
The path may consist of an arbitrary number of path segments alternating
forward and backward drives of the car. We are interested in a close to
optimal parking maneuver path in the sense of path length respecting the
kinematic constraints of the car.

In this paper, we propose an RRT*-based algorithm to solve the autonomous
parking problem, which we define more formally in Section 2. Contrary to well-
known A* algorithm, RRT* algorithm does not need space discretization. Also,
it handles nonholonomic constraints by design. The RRT* algorithm searches the
state space by creating a tree structure that represents possible paths.
RRT*-based algorithms were successfully applied to a wide range of planning
problems from the robot, vehicle, and aerial domains. However, they were also
used in not such apparent problems as tunnel detection in proteins from the
field of molecular biology.

Our algorithm uses Reeds and Shepp curves for particular path segments when
building the tree and Euclidean distance as a metric for the nearest neighbor
search. We complemented the RRT* algorithm with an optimization procedure
based on the Dijkstra algorithm used to reduce the number of the path segments
and to lower the cost of the path connecting initial and goal pose.

The main contributions of this paper are:

  * •

Minimization of the path cost with an optimization procedure based on the
shortest path by Dijkstra algorithm.

  * •

Speed up of the RRT* path search with the nearest neighbor heuristics.

In our experiments (see Section 5), we compare multiple cost functions of the
nearest neighbor search and show that the fastest approach to find the path is
to use the Euclidean distance as the cost function in the nearest neighbor
search (see Figure 4). We also evaluate the effectiveness of our optimization
procedure based on the Dijkstra algorithm and show (see Figure 5) that it
significantly improves the cost of the path even when compared to other
algorithms such as RRT*-Smart. In Section 6 we summarize our results. The
source code of our algorithm is available
111http://rtime.felk.cvut.cz/gitweb/hubacji1/iamcar.git.

###  1.1 Related Works

A common approach to solve a parking problem is to split the task to the
environment detection, the path planning, and the path execution. In this
paper, we consider the path planning part.

Typical parking problems can be classified into two classes: parallel parking
and perpendicular parking. Some publications consider only parallel parking
[Gupta et al., 2010], [Cheng et al., 2013], [Vorobieva et al., 2013], or only
perpendicular parking [Petrov et al., 2015]. In this paper, we propose a
universal method which considers obstacles of arbitrary shape.

Many published approaches use Reeds and Shepp curves [Reeds and Shepp, 1990]
for path planning [Lee et al., 2006] without considering obstacles. In
[Fraichard and Scheuer, 2004], the authors present Continuous-Curvature Paths
that extend the Reeds and Shepp line segments and circular arcs with clothoid
arcs. Resulting paths have continuous curvature, so a car that follows a path
does not have to stop to change orientation of the wheels. Continuous-
Curvature Paths have been used in [Muller et al., 2007], [Vorobieva et al.,
2013], [Cheng et al., 2013], and [Yi et al., 2017]. In [Kim et al., 2010], the
authors use two basic motions to create a set of motions. Finally, [Hsu et
al., 2008], [Gupta et al., 2010], and [Liang et al., 2012] describe parking
using paths generated with two circles geometry.

However, in real-life situations, a typical parking scenario may be disturbed
by sloppy parked neighbor car, temporary parked bike, non-standard parking
slot shape, or other unspecified constraints. Therefore, when a parking slot
is detected, evaluation of a situation may fail, and an approach based on
geometric equations may become unusable in such a case.

In this paper, we propose RRT*-based algorithm which can handle complex
parking situations. Rapidly-Exploring Random Trees (RRT) [LaValle, 1998] is a
randomized algorithm that can handle nonholonomic constraints. Although RRT is
probabilistically complete (with probability $1\text{\,}$, the algorithm
converges to the solution, as time tends to infinity), it is not
asymptotically optimal [Karaman and Frazzoli, 2011]. Therefore, Karaman and
Frazzoli proposed the RRT* algorithm, which converges to an optimal solution
as time tends to infinity. In [Islam et al., 2012], the authors improved the
RRT* algorithm by using path optimization and intelligent sampling and named
the resulting algorithm RRT*-Smart. After the initial path is found,
RRT*-Smart converges to the optimum faster than RRT*. In our approach, we stop
the RRT* algorithm when a path is found, and then we optimize the path by
Dijkstra algorithm.

##  2 THE PARKING PROBLEM

In this section, we define the problem and terminology used throughout this
paper.

Figure 1: Example scenario with the init pose (1), the goal pose (2), two
circle obstacles (3), the obstacle compound of line segment obstacles (4), and
the final path (5).

A pose is a triplet $p=(x,y,\theta)$, where $x,y$ are cartesian coordinates
and $\theta$ is a heading.

A search space is a set of poses $S=\\{(x,y,\theta)\mid
x\in[XMIN,XMAX],y\in[YMIN,YMAX],\theta\in[0,2\pi)\\}$, where $XMIN$, $XMAX$,
$YMIN$, and $YMAX$ are borders of search space.

A scenario is a quintuple $s=(S,p_{init},p_{goal},O_{C},O_{S})$, where $S$ is
a search space, $p_{init},p_{goal}$ are init and goal poses, and $O_{C}$ resp.
$O_{S}$ are sets of circle obstacles resp. line segment obstacles. We can see
an example scenario with the final path connecting initial and goal pose in
Figure 1. Example scenario also demonstrates segment obstacles (borders),
circle obstacles, and the complex obstacle of arbitrary shape (compound of
line segment obstacles).

Circle obstacle is a triplet $o_{c}=(x,y,r)$, where $x,y$ are cartesian
coordinates of the center, and $r$ is the radius. Line segment obstacle is a
quadruple $o_{s}=(x_{1},y_{1},x_{2},y_{2})$, where $x_{1},y_{2}$ are
coordinates of the line segment start and $x_{2},y_{2}$ are coordinates of the
line segment end.

A car is a quadruple $c=(l,w,R,b)$, where $l$ is a length of the car, $w$ is a
width of the car, $R=\frac{1}{\kappa}$ is car minimum turning radius, $\kappa$
is curvature, and $b$ is car wheelbase (the distance between front and rear
axles). In Figure 1, red crosses represent $x,y$ coordinates of $init$ and
$goal$ poses. The U-Shape frame represents length $l$, width $w$, and pose
heading. Finally, example obstacles are hatched.

A path from pose $a$ to pose $b$ is a sequence of poses $P_{a,b}=\\{p_{i}\mid
i\in\\{0,1,...,n-1\\},p_{0}=a,p_{n-1}=b\\}$, such that $P$ satisfies kinematic
constrains given by car $c$.

The $collide(p,O)$ function returns $True$ when a car $c$ positioned at pose
$p$ is inside arbitrary obstacle $o\in O$, or the frame of car $c$ collides
with this obstacle. Otherwise, the function returns $False$.

The $collide(P,O)$ function returns $True$ when for any pose $p\in P$ the
$collide(p,O)$ returns $True$. Otherwise, the function returns $False$.

The $cost(P)$ is a path cost defined in Equation 1, where $RSDist(a,b)$ is
Reeds & Shepp distance from pose $a$ to pose $b$.

$cost(P)=\sum_{i=0}^{i=n-2}RSDist(p_{i},p_{i+1})$ (1)

We define final path $P_{F}$ in Equation 2, where $P_{all}=\\{P_{a,b}\mid
a=p_{init},b=p_{goal},\neg collide(P,O)\\}$. An example of the final path is
in Figure 1.

$P_{F}=\underset{P\in P_{all}}{{\arg\min}_{P}}\ cost(P)$ (2)

##  3 RRT*

Rapidly-Exploring Random Tree Star (RRT*) is an asymptotically optimal
randomized algorithm to solve path planning problems, such as the parking
problem defined in Section 2.

RRT* uses a tree data structure that represents poses and paths, it handles
nonholonomic constraints and can hold general restrictions on $p_{init}$ and
$p_{goal}$ poses, or obstacles. Therefore, the RRT* should be able to solve
the unpredictable, real-life scenarios. We can see basic RRT* pseudocode
(lines 4 to 19) as part of complete RRT*-based Algorithm 1.

The fundamental element of RRT* is a $node$. The $node$ is a pose extended
with $parent$ (the pointer to the predecessor $node$), $children$ (the array
of successor $nodes$), and cumulative cost $ccost=cost(P_{p_{init},node})$. As
$node$ is extension to pose, we may update our definition of path
$P_{a,b}=\\{p_{i}\mid i\in\\{0,1,...,n-1\\},p_{0}=a,p_{n-1}=b\\}$, such that
$a$, $b$, and $p_{i}$ are nodes, where $p_{i}$ is parent of $p_{i+1}$. We use
a $path$ as the sequence of poses or $nodes$ interchangeably.

In RRT* algorithm, all $nodes$ are stored in tree data structure
$\mathcal{T}=(root,V,E)$, where $root$ node corresponds to $p_{init}$ pose,
$V$ is set of nodes, $E$ is set of edges, and $\forall n_{1},n_{2}\in
V:\\{n_{1},n_{2}\\}\in E\Leftrightarrow n_{1}\text{ is the parent of }n_{2}$.

Figure 2: Example scenario with the init pose (1), the goal pose (2), two
circle obstacles (3), the obstacle compound of line segment obstacles (4), the
final path (5), the final path before optimization (6), and line segments and
circle segments (gray).

###  3.1 Basic Procedures

In this section, we describe the basic procedures of RRT* used to build
$\mathcal{T}$ data structure.

RandomSample procedure returns a node with a pose from search space $S$, where
$x$, $y$, and $\theta$ are randomly generated.

Cost($nn,rs$) function is a metric used in RRT*.

NearestNeighbor($rs$) procedure searches for a $node$ with the lowest
Cost($node,rs$) in $\mathcal{T}$.

Steer($nn,rs$) procedure returns a path $P_{nn,rs}$. We can see the results of
Steer procedure in Figure 2 (gray).

NearNodes($ns,dist$) procedure returns a set of nodes $nns$ from
$\mathcal{T}$, such that $\forall n\in nns:\textsc{Cost}(n,ns)<dist$.

Connect($ns,nns$) procedure searches in near nodes ($nns$) for the best
candidate node to expand $\mathcal{T}$ towards the $ns$. The best candidate
node is the node in $\mathcal{T}$ that minimizes the cumulative cost of $ns$
when it becomes the parent of the $ns$. The path from the best candidate node
to the $ns$ must be free of collisions. If the best candidate node is found,
the $na$ is added as a child of the best candidate node, and Connect($ns,nns$)
returns $True$. Otherwise, the procedure returns $False$.

Rewire($ns,nns$) procedure checks if for any $node$ in $nns$ there is a path
with lower cumulative cost via $ns$. And swaps parents if so. This procedure
along with Connect($ns,nns$) ensures the asymptotical optimality of RRT*.

GoalFound returns $True$ if $p_{goal}\in\mathcal{T}$ and $False$ otherwise.

Collides($na,ns$) returns $True$ if the path $P_{na,ns}$ collides with any
obstacle of scenario, and $False$ otherwise.

Algorithm 1 Accelerated RRT*

1:Input:

  * •

initial pose

  * •

goal pose

  * •

array of obstacles

2:Output:

  * •

$True$ if goal pose reached, $False$ otherwise

  * •

array of paths connecting initial and goal pose

3:procedure RRT*

4: while $\textsc{Elapsed}<TMAX$ do

5: $rs\leftarrow\textsc{RandomSample}$

6: $nn\leftarrow\textsc{NearestNeighbor}(rs)$

7: $pn\leftarrow nn$

8: $newNodes\leftarrow\emptyset$

9: for $ns\leftarrow\textsc{Steer}(nn,rs)$ do

10: $nns\leftarrow pn\cup\textsc{NearNodes}(ns,dist)$

11: if $\textsc{Connect}(ns,nns)$ then

12: $\textsc{Rewire}(ns,nns)$

13: $newNodes\leftarrow newNodes\cup ns$

14: if GoalFound then

15: break while

16: end if

17: $pn\leftarrow ns$

18: end if

19: end for

20: for $na\leftarrow newNodes$ do

21: $pn\leftarrow na$

22: for $ns\leftarrow\textsc{Steer}(na,goal)$ do

23: if $\textsc{Collide}(pn,ns)$ then

24: break

25: end if

26: $pn.children\leftarrow pn.children\cup ns$

27: if GoalFound then

28: break while

29: end if

30: $pn\leftarrow ns$

31: end for

32: end for

33: end while

34: if GoalFound then

35: OptPath

36: end if

37: return GoalFound

38:end procedure

###  3.2 Implementation

Our RandomSample procedure samples randomly from the whole space $S$
(including obstacles). We use OMPL [Sucan et al., 2012] implementation of
Reeds and Shepp [Reeds and Shepp, 1990] optimal paths for Steer and Cost
functions. NearNodes, Connect and Rewire procedures work the same as in
[Karaman and Frazzoli, 2011].

For two nodes we implemented auxiliary IsNear($n_{1},n_{2}$) function that
returns $True$ if $n_{1}$ is within the predefined Euclidean distance from
$n_{2}$ ($GFDIST$) and the difference between headings of $n_{1}$ and $n_{2}$
is less than the specified angle ($GFANGLE$). We use this function to specify
if the goal was found, the Steer procedure reached $rs$, or if two nodes are
the same. For computational experiments in Section 5, we used $GFDIST=0.05$
and $GFANGLE=\frac{\pi}{32}$.

In each iteration of RRT*-based algorithm, there is an expansion of
$\mathcal{T}$ towards the $p_{goal}$ (see lines 20 to 32 in Algorithm 1) as
used in [Kuwata et al., 2008]. We added path optimization procedure to
RRT*-based algorithm (see line 35 in Algorithm 1) that is run when the goal is
found as explained in Section 4.2. We can see an example of optimized final
path (5) and final path before optimization (6) in Figure 2.

##  4 RRT* IMPROVEMENTS

In this section, we introduce our improvement to nearest neighbor search and
details about path optimization procedure.

###  4.1 Nearest Neighbor

Algorithm 2 Nearest neighbor search

1:$IYSIZE$$\triangleright$ nn structure size

2:$IYSTEP$$\triangleright$ increment distance

3:$nodes[IYSIZE]$$\triangleright$ array of lists of nodes

4:

5:Input:

  * •

$node$ to be added to data structure

6:Output:

  * •

data structure of nodes

7:procedure AddIY($node$)

8: $iy\leftarrow\lfloor\frac{node.y}{IYSTEP}\rfloor$

9: $nodes[iy]\leftarrow nodes[iy]\cup node$

10:end procedure

11:

12:Input:

  * •

$node$ to be searched

13:Output:

  * •

the nearest neighbor of $node$

14:procedure NearestNeighbor($node$)

15: $iy\leftarrow\lfloor\frac{node.y}{IYSTEP}\rfloor$

16: $nn\leftarrow NULL$$\triangleright$ nearest neighbor

17: $c_{min}\leftarrow\infty$$\triangleright$ minimum cost

18: $as\leftarrow 0$$\triangleright$ array step

19: while $c_{min}>as\cdot IYSTEP$ do

20: $i\leftarrow\max(iy-as,0)$

21: $j\leftarrow\min(iy+as,IYSIZE-1)$

22: for $n\in nodes[i]\cup nodes[j]$ do

23: if $\textsc{EDist}(n,node)<c_{min}$ then

24: $c_{min}\leftarrow\textsc{EDist}(n,node)$

25: $nn\leftarrow n$

26: end if

27: end for

28: $as\leftarrow as+1$

29: end while

30:end procedure

Because the nearest neighbor procedure returns a node with the lowest cost,
such a node is a good candidate for tree expansion. The pseudocode of the
nearest neighbor search is outlined in Algorithm 2. To improve the performance
of finding the nearest neighbor, we use a $nodes$ data structure (the array of
linked lists of nodes) defined in line 3. The $nodes$ data structure allows us
to split search space $S$ along the $y$-axis ($y-axis$ suits better for
parallel parking scenario we experimented with in Section 5), so we can
compare nodes within multiples of $IYSTEP$ (increment distance based on
$nodes$ data structure) constant first.

Lines 7 to 10 describes how a $node$ is added to $nodes$. First, we compute
the index of $nodes$ array ($iy$) where the $node$ should be stored. Then, the
$node$ is added to the list of nodes at that $iy$ index.

When looking for the nearest neighbor of the $node$ in the indexing structure
(lines 14 to 30), we compute $iy$ index again. Then, we search the list of
nodes stored in the array $nodes$ on index $iy$ ($nodes[iy]$). Finally, we
repeatedly widen the interval of indexes to be investigated while the minimum
cost is higher than half of the interval width times $IYSTEP$ and search the
lists of nodes stored in the array on indexes corresponding to the widened
interval.

We use Euclidean distance as the cost function in the nearest neighbor search
in contrast to Reeds and Shepp path length as the cost function for building
RRT*. This approach speeds up the process but does not influence the final
path cost as discussed in Section 5.

###  4.2 Path Optimization

Algorithm 3 Path optimization

1:Input:

  * •

path connecting initial and goal pose

2:Output:

  * •

lower cost path connecting initial and goal pose

3:procedure OptPath

4: $tips\leftarrow$ cusp nodes$\triangleright$ array

5: $pq\leftarrow\emptyset$$\triangleright$ priority queue

6: $pq\leftarrow pq\cup tips[0]$

7: while $|pq|\neq 0$ do

8: $n_{i}\leftarrow\textsc{Pop}(pq)$

9: if $n_{i}=tips[\textsc{Size}(tips)-1]$ then

10: break

11: end if

12: for all $j>i$ do

13: $n_{j}\leftarrow tips[j]$

14: $P_{n_{i},n_{j}}\leftarrow\textsc{Steer}(n_{i},n_{j})$

15: $c\leftarrow n_{i}.ccost+\textsc{Cost}(n_{i},n_{j})$

16: if $\textsc{Collide}(n_{i},n_{j})$ then

17: continue

18: end if

19: if $c<n_{j}.ccost$ then

20: $n_{j}.ccost\leftarrow c$

21: $n_{j}.parent\leftarrow i$

22: if $n_{j}.visited=False$ then

23: $n_{j}.visited\leftarrow True$

24: $pq\leftarrow pq\cup n_{j}$

25: end if

26: end if

27: end for

28: end while

29: $opath\leftarrow\emptyset$$\triangleright$ new optimized path

30: $i\leftarrow\textsc{Size}(tips)-1$

31: while $i>0$ do

32: $opath\leftarrow opath\cup tips[i]$

33: $i\leftarrow tips[i].parent$

34: end while

35: $opath\leftarrow opath\cup tips[0]$

36: if better cost of $opath$ then

37: return True

38: end if

39: return False

40:end procedure

The path optimization procedure is run when the goal is found. Even that RRT*
is asymptotically optimal, it converges to the optimal solution very slowly.
When the $p_{goal}$ is reached for the first time a final path $P_{F}$ is
probably far from optimum in the sense of cost (we use Reeds and Shepp path
length as cost). The purpose of path optimization procedure is to decrease the
final path cost.

A final path $P_{F}$ consists of topologically ordered nodes (see definition
of a path in Section 2). We select tip nodes from the final path that are also
topologically ordered. In our case, tip nodes are cusp nodes (nodes where the
direction of movement changes) along with $p_{init}$ and $p_{goal}$.

In Algorithm 3 we initialize tip nodes and priority queue in lines 4 to 6. In
lines 7 to 28, we use Dijkstra algorithm to find the shortest path from the
first tip node ($p_{init}$) to the last one ($p_{goal}$). From the priority
queue, we pop the node $n_{i}$ (where $i$ is the index of $n_{i}$ node in
$tips$ array) with the lowest cumulative cost. Then, we call
Steer($n_{i},n_{j}$) procedure from $n_{i}$ to all $n_{j}$ for $j>i$ (see
lines 12 to 27) that returns path $P_{n_{i},n_{j}}$. If $P_{n_{i},n_{j}}$ is
collision free and cumulative cost of $n_{j}$ is smaller when reached via
$P_{n_{i},n_{j}}$ then parent and cumulative cost of $n_{j}$ are updated, and
$n_{j}$ is pushed to the priority queue if not visited already. The process
repeats until the priority queue is empty or $n_{i}=p_{goal}$.

Optimized path found by Dijkstra is retrieved in lines 29 to 35. If the
cumulative cost of $p_{goal}$ is better, OptPath procedure returns $True$ and
$False$ otherwise.

##  5 COMPUTATIONAL EXPERIMENTS AND EVALUATION

We present the results of computational experiments for parallel parking
scenario with no obstacle in Section 5.1, and the results of computational
experiments for parallel parking scenarios with circle obstacle in Section
5.2.

We are interested in the nearest neighbor search and path optimization
procedures. Specifically, we are interested in how does the cost function,
used in the nearest neighbor search, influences algorithm computation time. We
experimented with the following implementations of the nearest neighbor
search:

  * •

Nearest neighbor search with the cost based on Reeds and Shepp path length.

  * •

Nearest neighbor search with the cost based on Reeds and Shepp path length but
with the heading of nodes temporarily set to the same value.

  * •

Nearest neighbor search with the cost based on Euclidean distance.

Also, we would like to know if the path optimization procedure influences the
cost of the final path. We tested the following path optimization
possibilities:

  * •

No path optimization.

  * •

Path optimization from [Islam et al., 2012].

  * •

Path optimization described in Algorithm 3.

The car we use for experiments is $1.625\text{\,}\mathrm{m}$ wide and
$3.760\text{\,}\mathrm{m}$ long. The minimum turning radius of the car is
$10.820\text{\,}\mathrm{m}$ and wheelbase is $2.450\text{\,}\mathrm{m}$.

We run computational experiments on a single core of Intel(R) Core(TM)
i7-5600U CPU @ $2.60\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ with MemTotal:
$16\,322\,516\text{\,}\mathrm{k}\mathrm{B}$.

###  5.1 Scenario with No Obstacle

We tested RRT*-based algorithm on parallel parking scenario shown in Figure 3.
The parking lot is $2.2\text{\,}\mathrm{m}$ wide and $6.5\text{\,}\mathrm{m}$
long [CSN 73 6056, 2011]. The width of the street is
$2.75\text{\,}\mathrm{m}$.

We let the Algorithm 1 to run for up to 10 seconds. When the RRT*-based
algorithm finds the goal, the OptPath procedure optimizes the final path. We
repeated the experiment $10\,000\text{\,}$ times for this scenario.

Figure 3: Parallel parking scenario with no obstacle. On the left, there is a
final path before optimization (orange) and optimized final path (blue). On
the right, there is a complete tree of all paths (gray).

####  5.1.1 Nearest Neighbor Search

We compare the computation times when the algorithm found the final path for
different cost functions used in the nearest neighbor search implementations.
We can see the results in the histogram with the logarithmic scale in Figure
4.

Figure 4: Histogram of time to find the path. Vertical dashed lines represent
$95\text{\,}\%$ percentile (red is 0.13, blue is 1.16, orange is 5.98).

For the nearest neighbor search implementation with the Reeds and Shepp cost
function (the same cost function used for building $\mathcal{T}$, orange in
Figure 4), the algorithm did not find the goal in all runs. On the other hand,
for the nearest neighbor search implementation where we used the Euclidean
distance as the cost function (red in Figure 4), the goal was found in
$100\text{\,}\%$ of runs. For comparison purposes, we run the experiment for
the nearest neighbor search implementation with Reeds and Shepp cost function,
where the heading of the nodes was temporarily set to the same value (blue in
Figure 4).

####  5.1.2 Path Optimization

We also compared the final path costs for different path optimization
procedures. We can see the results in the histogram with the logarithmic scale
in Figure 5.

Figure 5: Path cost histogram. Vertical dashed lines represent $95\text{\,}\%$
percentile (red is 33.29, blue is 37.56, orange is 40.31).

We can see the improvement over no path optimization (orange) when algorithm
from [Islam et al., 2012] is used (blue). And we can see that the path
optimization from Algorithm 3 (red) has the best results.

###  5.2 Scenario with Circle Obstacle

Further, we tested RRT*-based algorithm on parallel parking scenarios shown in
Figure 6. The parking lot is $2.2\text{\,}\mathrm{m}$ wide and
$6.5\text{\,}\mathrm{m}$ long [CSN 73 6056, 2011]. The width of the street is
$2.75\text{\,}\mathrm{m}$. There is a random circle obstacle with diameter of
$0.5\text{\,}\mathrm{m}$ laying on the street near the parking lot.

We let the Algorithm 1 to run for up to 10 seconds. When the RRT*-based
algorithm finds the goal, the OptPath procedure optimizes the final path. We
repeated the experiment $10\,000\text{\,}$ times for this scenario.

Figure 6: Parallel parking scenarios with circle obstacle. Scenarios 1 and 2
differ in the position of circle obstacle. There is the final path before
optimization (orange) and the optimized final path (blue).

The results are similar to the results in Section 5.1. The cost based on the
Euclidean distance speeds up the algorithm computation time ($95\text{\,}\%$
percentile), dependent on the obstacle position, to $0.12\text{\,}\mathrm{s}$
for Scenario 1 and to $0.66\text{\,}\mathrm{s}$ for Scenario 2. The path
optimization procedure decreases the cost of the final path ($95\text{\,}\%$
percentile) by $4\text{\,}\%$ concerning No optimization case for both
scenarios.

##  6 CONCLUSION

We proposed the RRT*-based algorithm for planning parking paths and
experimented with the nearest neighbor search and path optimization
procedures.

RRT*-based algorithm without improvements uses the cost function based on the
Reeds and Shepp path length in the nearest neighbor search (as well as for
building the $\mathcal{T}$ data structure), and no optimization procedure. Our
improvements include the cost function based on Euclidean distance in the
nearest neighbor search and optimization procedure based on the Dijkstra
algorithm.

We have shown that when we use the cost function based on the Reeds and Shepp
path length for building the $\mathcal{T}$ data structure and the cost
function based on the Euclidean distance in the nearest neighbor search, there
is a significant acceleration in algorithm computation time.

Additionally, we have shown that the path optimization procedure based on the
Dijkstra algorithm for the shortest path search can optimize the final path to
$63\text{\,}\%$ of the original cost in $95\text{\,}\%$ of cases, for the
parallel parking scenario without obstacles, which is a better result than the
optimization procedure used in [Islam et al., 2012]. However, for parallel
parking scenario with circle obstacle, the optimized cost is only
$96\text{\,}\%$ of the original cost in $95\text{\,}\%$ of cases.

Finally, from the experiments we can see that for parallel parking scenario
with no obstacle, RRT*-based algorithm with improvements tends to
significantly faster computation time as well as to lower final path cost.
However, for parallel parking scenario with circle obstacle, RRT*-based
algorithm with improvements tends to significantly faster computation time but
about $10\text{\,}\%$ to $20\text{\,}\%$ worse final path cost then RRT*-based
algorithm without improvements.

In our future work, we are going to experiment with the improvements presented
in this paper. Particularly, the recognition and selection of tip nodes seem
to be interesting. Also, the bidirectional RRT* algorithms, such as [Jordan
and Perez, 2013] and [Klemm et al., 2015], could lead to significant
improvements in the matter of computational time.

## ACKNOWLEDGEMENTS

This work was supported by the Technology Agency of the Czech Republic under
the Centre for Applied Cybernetics TE01020197.

## References

  * [Cheng et al., 2013] Cheng, K., Zhang, Y., and Chen, H. (2013).  Planning and control for a fully-automatic parallel parking assist system in narrow parking spaces.  In Proc. IEEE Intelligent Vehicles Symp. (IV), pages 1440–1445.
  * [CSN 73 6056, 2011] CSN 73 6056 (2011).  Parking areas for road vehicles.  Technical report, Praha.
  * [Fraichard and Scheuer, 2004] Fraichard, T. and Scheuer, A. (2004).  From reeds and shepp’s to continuous-curvature paths.  20(6):1025–1035.
  * [Gupta et al., 2010] Gupta, A., Divekar, R., and Agrawal, M. (2010).  Autonomous parallel parking system for ackerman steering four wheelers.  In Proc. IEEE Int. Conf. Computational Intelligence and Computing Research, pages 1–6.
  * [Hsu et al., 2008] Hsu, T., Liu, J., Yu, P., Lee, W., and Hsu, J. (2008).  Development of an automatic parking system for vehicle.  In Proc. IEEE Vehicle Power and Propulsion Conf, pages 1–6.
  * [Islam et al., 2012] Islam, F., Nasir, J., Malik, U., Ayaz, Y., and Hasan, O. (2012).  RRT*-Smart: Rapid convergence implementation of RRT* towards optimal solution.  In Proc. IEEE Int. Conf. Mechatronics and Automation, pages 1651–1656.
  * [Jordan and Perez, 2013] Jordan, M. and Perez, A. (2013).  Optimal bidirectional rapidly-exploring random trees.
  * [Karaman and Frazzoli, 2011] Karaman, S. and Frazzoli, E. (2011).  Sampling-based algorithms for optimal motion planning.  The international journal of robotics research, 30(7):846–894.
  * [Kim et al., 2010] Kim, D., Chung, W., and Park, S. (2010).  Practical motion planning for car-parking control in narrow environment.  IET Control Theory Applications, 4(1):129–139.
  * [Klemm et al., 2015] Klemm, S., Oberländer, J., Hermann, A., Roennau, A., Schamm, T., Zollner, J. M., and Dillmann, R. (2015).  Rrt $\ast$-connect: Faster, asymptotically optimal motion planning.  In Proc. IEEE Int. Conf. Robotics and Biomimetics (ROBIO), pages 1670–1677.
  * [Kuwata et al., 2008] Kuwata, Y., Fiore, G. A., Teo, J., Frazzoli, E., and How, J. P. (2008).  Motion planning for urban driving using rrt.  In Proc. IEEE/RSJ Int. Conf. Intelligent Robots and Systems, pages 1681–1686.
  * [LaValle, 1998] LaValle, S. M. (1998).  Rapidly-exploring random trees: A new tool for path planning.
  * [Lee et al., 2006] Lee, K., Kim, D., Chung, W., Chang, H. W., and Yoon, P. (2006).  Car parking control using a trajectory tracking controller.  In Proc. SICE-ICASE Int. Joint Conf, pages 2058–2063.
  * [Liang et al., 2012] Liang, Z., Zheng, G., and Li, J. (2012).  Automatic parking path optimization based on bezier curve fitting.  In Proc. IEEE Int. Conf. Automation and Logistics, pages 583–587.
  * [Muller et al., 2007] Muller, B., Deutscher, J., and Grodde, S. (2007).  Continuous curvature trajectory design and feedforward control for parking a car.  15(3):541–553.
  * [Petrov et al., 2015] Petrov, P., Nashashibi, F., and Marouf, M. (2015).  Path planning and steering control for an automatic perpendicular parking assist system.  In 7th Workshop on Planning, Perception and Navigation for Intelligent Vehicles, PPNIV, volume 15, pages 143–148.
  * [Reeds and Shepp, 1990] Reeds, J. and Shepp, L. (1990).  Optimal paths for a car that goes both forwards and backwards.  Pacific journal of mathematics, 145(2):367–393.
  * [Sucan et al., 2012] Sucan, I. A., Moll, M., and Kavraki, L. E. (2012).  The open motion planning library.  IEEE Robotics Automation Magazine, 19(4):72–82.
  * [Vorobieva et al., 2013] Vorobieva, H., Minoiu-Enache, N., Glaser, S., and Mammar, S. (2013).  Geometric continuous-curvature path planning for automatic parallel parking.  In Proc. SENSING AND CONTROL (ICNSC) 2013 10th IEEE INTERNATIONAL CONFERENCE ON NETWORKING, pages 418–423.
  * [Yi et al., 2017] Yi, Y., Lu, Z., Xin, Q., Jinzhou, L., Yijin, L., and Jianhang, W. (2017).  Smooth path planning for autonomous parking system.  In Proc. IEEE Intelligent Vehicles Symp. (IV), pages 167–173.


# Counterfactuals for the Future

Lucius E.J. Bynum, 1 Joshua R. Loftus, 2 Julia Stoyanovich 1

###### Abstract

Counterfactuals are often described as ‘retrospective,’ focusing on
hypothetical alternatives to a realized past. This description relates to an
often implicit assumption about the structure and stability of exogenous
variables in the system being modeled — an assumption that is reasonable in
many settings where counterfactuals are used. In this work, we consider cases
where we might reasonably make a different assumption about exogenous
variables, namely, that the exogenous noise terms of each unit do exhibit some
unit-specific structure and/or stability. This leads us to a different use of
counterfactuals — a ‘forward-looking’ rather than ‘retrospective’
counterfactual. We introduce “counterfactual treatment choice,” a type of
treatment choice problem that motivates using forward-looking counterfactuals.
We then explore how mismatches between interventional versus forward-looking
counterfactual approaches to treatment choice, consistent with different
assumptions about exogenous noise, can lead to counterintuitive results.

##  1 Introduction

Counterfactuals are often described as retrospective, focusing on hypothetical
alternatives to a realized past. Here, we explore how this view of
counterfactuals as retrospective connects to implicit assumptions about how
exogenous (noise) variables behave in the system being modeled. We frame our
discussion in terms of graphical approaches to causality (Pearl 2009; Spirtes
et al. 2000; Peters, Janzing, and Schölkopf 2017). We characterize assumptions
about such exogenous variables as being of two kinds: (1) assumptions about
how a unit looks exogenously _compared to itself over time_ , and (2)
assumptions about how a unit looks exogenously _compared to other units_. We
refer to the first kind as assumptions about ‘stability’, and to the second as
assumptions about ‘structure.’ In characterizing these assumptions, we explore
how they lead us to a ‘forward-looking’ rather than a ‘retrospective’ use of
counterfactuals.

The questions we explore in this work are motivated by a distinction between
two types of approaches to maximizing social welfare (or minimizing
inequality) in a target population by learning a treatment rule. We now
contrast these two approaches with an example.

Suppose you want to allocate resources within a particular population to
improve an outcome. For a concrete example, consider allocating tutoring to
students at a school in order to improve test performance. Further, suppose
that unobserved external factors about the students (e.g., family income,
encouragement from parents) explain at least some of the variation in the
outcome across units. How should we allocate resources to these students?

We frame this as a computational problem involving a few steps. The input is a
dataset describing units, in the style of an observational study at one time
step, along with assumptions about causal structure (which could also be
learned from the data). The output is a treatment assignment across units, for
example, which students get tutoring. The optimization problem operates over
an exponential space of possible treatment assignments and involves estimating
the impact of a given assignment on the outcome distribution. In this work, we
are agnostic to the particular optimization objective. However, we focus on
objectives that would require the entire outcome distribution, which is of
particular interest if we care about mitigating inequality rather than
maximizing utility. Our work here concerns the following question: _Which
approach — interventional or counterfactual — should we use to estimate the
outcome distribution?_ (See Section 2 for precise definitions of
interventional and counterfactual distributions.)

_Approach 1._ We could use our model of the data to identify which subgroups
of students to treat based on their covariate values only. After learning an
optimal treatment rule, we would apply it to our observed students. We could
describe this as a treatment choice problem focused on _interventional
distributions_ — a problem studied extensively in econometrics. We may take an
approach to tackling this problem such as Empirical Welfare Maximization (EWM)
that, in short, maximizes the social welfare in a target population using a
sampled population to learn a treatment assignment rule (Kitagawa and Tetenov
2018). Notice that here, in contrast to much work in econometrics, our sample
_is_ our entire population of interest.

_Approach 2._ We could use our model of the data to identify which students to
treat based on their covariate values _in addition to_ any student-specific
idiosyncrasies not captured by the covariates we have measured. This gives
rise to a problem focused on _counterfactual distributions_ , which are
particular to the observed students and their modeled exogenous variables,
i.e., noise.

Importantly, _Approach 1_ and _Approach 2_ will in general produce a different
estimate of the outcome distribution and, crucially, lead us to tutor a
different set of students. In other words, these two approaches will result in
two different policies. In this work, we characterize how accurate or
inaccurate each approach would be depending on our assumptions about exogenous
structure and stability (i.e., what do we assume about family income and
encouragement from parents year-to-year?). We show how mismatched assumptions
between interventional and counterfactual approaches to treatment choice can
lead to counterintuitive results.

#### Contributions.

  * •

We make a case for using counterfactuals in an explicitly forward-looking
rather than retrospective manner in particular policy choice settings, and
connect the use of counterfactual (versus interventional) distributions to
assumptions about exogenous noise (Section 3).

  * •

We introduce a distinct type of decision problem focused on inequality —
_counterfactual treatment choice_ — that motivates the use of forward-looking
counterfactuals (Section 4).

  * •

We study empirically the exogenous conditions under which interventional
approaches versus forward-looking counterfactual approaches to treatment
choice problems are more appropriate (Section 5).

##  2 Background & Related Work

###  2.1 Structural Causal Models

Let capital letter $X$ denote a random variable, where lowercase letter $X=x$
denotes the value it obtains. Let boldface capital letter
$\mathbf{X}=\\{X_{1},\ldots,X_{n}\\}$ denote a set of random variables, with
value $\mathbf{X}=\mathbf{x}$. Capital $P_{X}$ denotes the cumulative
distribution function of variable $X$ and lowercase $p_{X}$ denotes the
density (or mass) function. Let $P_{Y\mid X=x}$ denote the conditional
distribution of $Y$ given $X=x$ and $P_{Y\mid X}$ denote the collection of
$P_{Y\mid X=x}$ for all $x$, i.e., the conditional of $Y$ given $X$. Next we
introduce the structural causal model framework and associated
notation.111Notation in this section is a combination of notation from Peters,
Janzing, and Schölkopf (2017) with language and notation from Buesing et al.
(2019) and Oberst and Sontag (2019).

######  Definition 1 (Structural Causal Model (SCM)).

An SCM $\mathfrak{C}$ is a four-tuple
$\left<\mathbf{U},\mathbf{X},\mathbf{F},P_{\mathbf{U}}\right>$ where
$\mathbf{U}$ is a set of independent random noise variables
$\mathbf{U}=\\{U_{1},\ldots,U_{n}\\}$ with distribution $P_{\mathbf{U}}$,
$\mathbf{X}$ is a set of random variables
$\mathbf{X}=\\{X_{1},\ldots,X_{n}\\}$, and $\mathbf{F}$ is a set of functions
$\mathbf{F}=\\{f_{1},\ldots,f_{n}\\}$. For all $i$,
$X_{i}=f_{i}(\mathbf{PA}_{X_{i}},U_{i})$, where
$\mathbf{PA}_{X_{i}}\subseteq\mathbf{X}\setminus X_{i}$ is the set of parents
of $X_{i}$ in the causal DAG $\mathcal{G}$. The prior distribution
$P_{\mathbf{U}}$ and functions $\mathbf{F}$ determine the distribution
$P^{\mathfrak{C}}$.

An SCM entails an _observational distribution_ $P_{\mathbf{X}}$ as well as
distributions after _interventions_.

######  Definition 2 (Interventions).

An intervention

$I=do\left(X_{i}\coloneqq\tilde{f}(\tilde{\mathbf{PA}}_{X_{i}},\tilde{U}_{i})\right)$

corresponds to replacing the structural mechanism
$f_{i}(\mathbf{PA}_{X_{i}},U_{i})$ with
$\tilde{f}_{i}(\tilde{\mathbf{PA}}_{X_{i}},\tilde{U}_{i})$. This definition
encompasses atomic interventions, denoted $\text{do}(X_{i}=x)$. We denote an
SCM after intervention $I$ as $\mathfrak{C}^{\text{do}(I)}$, and the resulting
interventional distribution as $P^{\mathfrak{C};\text{do}(I)}$.

An SCM also entails _counterfactual distributions_.

######  Definition 3 (Counterfactuals).

A counterfactual distribution is an interventional distribution specific to an
observed assignment $\mathbf{X}=\mathbf{x}$ (over any set of observed
variables). The counterfactual distribution, denoted
$P^{\mathfrak{C}\mid\mathbf{X}=\mathbf{x};\text{do}(I)}$, is the distribution
entailed by the SCM $\mathfrak{C}^{\text{do}(I)}$ using the posterior
distribution $P_{\mathbf{U}\mid\mathbf{X}=\mathbf{x}}$ instead of the prior
$P_{\mathbf{U}}$, i.e., the posterior distribution
$P_{\mathbf{U}\mid\mathbf{X}=\mathbf{x}}$ is passed through the modified
structural mechanisms in $\mathfrak{C}^{\text{do}(I)}$ to obtain the
counterfactual distribution specific to assignment $\mathbf{X}=\mathbf{x}$.
Note in the continuous case, we would condition instead on $\mathbf{X}\in A$
with $P(\mathbf{X}\in A)>0$ rather than $\mathbf{X}=\mathbf{x}$.

###  2.2 Welfare Maximization in Econometrics

Problems that assign treatment based on some function of people’s outcomes
(e.g., a measure of welfare) are often referred to as statistical treatment
rules, policy choice, or policy learning problems, and relate to a large body
of work in econometrics and causal inference (Athey and Imbens 2017; Manski
2003). Typically in this literature, policies are learned to maximize a
utilitarian (additive) welfare criterion rather than a criterion focused
explicitly on inequality. However, a growing collection of works look at
treatment choice with a focus on inequality, such as Kitagawa and Tetenov
(2019) and Kasy (2016), and inequality research itself has a long history.
Though these and other related works on welfare maximization have a similar
focus on inequality, in this work we show how econometric treatment choice
problems are often what we define as interventional treatment choice problems,
with a focus on estimating interventional rather than forward-looking
counterfactual distributions.

###  2.3 Algorithmic Fairness and Inequality

The distributional impacts of decisions, tensions between fairness and
inequality, and other questions of social welfare are widely discussed in
algorithmic fairness and algorithmic decision-making literature (Barocas,
Hardt, and Narayanan 2019; Bynum, Loftus, and Stoyanovich 2021; Green and Chen
2019; Heidari, Ferrari, Gummadi, and Krause 2019a, b; Hu and Chen 2020;
Kannan, Roth, and Ziani 2019; Kasy and Abebe 2021; Khan, Manis, and
Stoyanovich 2021; Kusner, Russell, Loftus, and Silva 2019; Liu, Dean, Rolf,
Simchowitz, and Hardt 2018; Madras, Pitassi, and Zemel 2018; Mullainathan
2018; Nabi, Malinsky, and Shpitser 2019). Specific examples of “counterfactual
treatment choice problems,” which we define generally here, have come up in
different instances in this literature, e.g., the “impact remediation problem”
introduced by Bynum, Loftus, and Stoyanovich (2021) and “discriminatory impact
problem“ introduced by Kusner et al. (2019). To the best of our knowledge,
these or other types of problems have not been discussed with connections to
forward-looking counterfactuals nor assumptions about exogenous variables.

##  3 Forward-looking Counterfactuals

We can demonstrate forward-looking counterfactuals through an example.
Consider a set of $n$ schools for which we observe data
$\\{Z_{0}^{(i)},Y_{0}^{(i)}\\}_{i=1}^{n}$ for covariates $Z,Y$ at time $t=0$,
described by SCM $\mathfrak{C}$ where

$\displaystyle Z^{(i)}_{0}$ $\displaystyle=U^{(i)}_{Z},$ $\displaystyle
Y^{(i)}_{0}$ $\displaystyle=Z^{(i)}_{0}+U^{(i)}_{Y},$

with exogenous noise variables $U^{(i)}_{Z},U^{(i)}_{Y}\sim\mathcal{N}(0,1)$.
At time $t=1$, we can intervene on $Z$ to determine the value of $Z_{1}$ with
the goal of changing future outcomes $Y_{1}$. Suppose that our intervention is
to add constant $\delta$ to the value of $Z$, i.e.,
$Z^{(i)}_{1}=Z^{(i)}_{0}+\delta\cdot w(i)$, with $w(i)\in\\{0,1\\}$ indicating
the treatment decision for individual $i$.

Consider using an interventional distribution to characterize the distribution
of $Y_{1}$ after performing an intervention $I$ that determines which set of
individuals to treat $\\{i:w(i)=1\\}$. The following modified SCM
$\mathfrak{C}^{\text{do}(I)}$,

$\displaystyle Z^{(i)}_{1}$ $\displaystyle=Z^{(i)}_{0}+\delta\cdot w(i),$
$\displaystyle Y^{(i)}_{1}$ $\displaystyle=Z^{(i)}_{1}+U^{\prime(i)}_{Y},$

with a sample from the prior $U^{\prime(i)}_{Y}\sim\mathcal{N}(0,1)$, induces
our interventional distribution of interest
$P_{Y_{1}}^{\mathfrak{C};\text{do}(I)}$. This distribution describes the
distribution of $Y_{1}$ for a group of individuals with the same post-
intervention covariates $Z_{1}$ as the observed sample.

Now consider using a counterfactual distribution to characterize the
distribution of $Y_{1}$ after performing intervention $I$. Instead of using
prior $U^{(i)}_{Y}\sim\mathcal{N}(0,1)$ in $\mathfrak{C}^{\text{do}(I)}$, we
have

$\displaystyle Z^{(i)}_{1}$ $\displaystyle=Z^{(i)}_{0}+\delta\cdot w(i),$
$\displaystyle Y^{(i)}_{1}$ $\displaystyle=Z^{(i)}_{1}+\tilde{U}^{(i)}_{Y},$

with posterior $\tilde{U}^{(i)}_{Y}=U^{(i)}_{Y}$, where, rather than sample a
new value of $U_{Y}$ for each individual, we set $\tilde{U}_{Y}=U_{Y}$ during
abduction. The entailed distribution
$P_{Y_{1}}^{\mathfrak{C}\mid\mathcal{X};\text{do}(I)}$ describes the
distribution of $Y_{1}$ for this specific set of individuals using their
observed exogenous noise terms from time $t=0$.

Which approach should we use to estimate $P_{Y_{1}}$? This depends on what we
assume about the structure and stability of our exogenous noise term $U_{Y}$.
Notice that from our modeling perspective, ‘noise’ in the structural equation
for $Y$ can also be described as unobserved factors outside of our model that
explain some of the variation in $Y$. Assumptions about ‘noise,’ then, can
also be framed as assumptions about any variables we have not measured that
explain some of the variation in our outcome of interest.

In this work, we discuss how units ‘look exogenously’ with a focus only on an
outcome of interest $Y$ and ignore exogenous variation in other variables for
now. We also focus for now on just two time steps, with observation at time
$t=0$ and intervention at time $t=1$.

If we assume that a unit’s exogenous variables are constant over time, e.g.,
that a unit will look exactly the same next year as it does this year, then
for that unit a counterfactual estimate of $Y$ would be appropriate to
estimate $Y_{1}$, while an interventional estimate (if there is any
variability in exogenous variables across units) would not. By contrast, if we
assume that a unit’s exogenous variables change significantly over time, a
counterfactual estimate of $Y$ might be an appropriate retrospective estimate
of $Y_{0}$ after intervention, but not an appropriate estimate of $Y_{1}$. We
refer to the use of a counterfactual to predict a _future_ rather than
hypothetical past outcome as a _forward-looking counterfactual_.

Whether or not assumptions about structure and stability matter depends also
on _who_ our population of interest is at time $t=1$. If at time $t=1$ we are
concerned with a new sample of individuals governed by the same SCM, it does
not matter how exogenous noise behaves for the original sample — we would not
use counterfactual distributions anyway.

For what type of problem, then, would we want to use forward-looking
counterfactuals, and when would these assumptions about exogenous noise
matter? In Section 4, we define a type of problem where a forward-looking
counterfactual approach would be a natural fit, and in Section 5, we explore
assumptions about exogenous stability and structure in such settings
empirically.

##  4 Counterfactual Treatment Choice

In this section, we define a type of inequality-focused treatment choice
problem that motivates the use of forward-looking counterfactuals.

Consider a decision maker $\mathbbm{D}$ who assigns interventions, and a
population at time $t=0$ characterized by the joint distribution of
$(Z^{(i)}_{0},X^{(i)}_{0},Y^{(i)}_{0}(0),Y^{(i)}_{0}(1))$ with covariates
$X^{(i)}_{0}\in\mathcal{X}\subset\mathbb{R}^{d_{x}}$, treatments/interventions
$Z^{(i)}_{0}\in\Omega_{Z}$, and potential outcomes
$Y^{(i)}_{0}(0),Y^{(i)}_{0}(1)\in\mathbb{R}$. The data are a size $n$ sample,
$\\{Z^{(i)}_{0},X^{(i)}_{0},Y^{(i)}_{0}\\}_{i=1}^{n}$.222Here we have used
notation building on the setup of (Kitagawa and Tetenov 2019). Assume the
variables $\\{Z,X,Y\\}$ follow a data generating process described by SCM
$\mathfrak{C}$ with corresponding causal graph (DAG) $\mathcal{G}$. Denote the
set of endogenous variables $\mathbf{V}=\\{Z,X,Y\\}$, the set of exogenous
variables $\mathbf{U}=\\{U_{Z},U_{X},U_{Y}\\}$, and the observed data
$\mathcal{X}=\\{\mathbf{v}^{(i)}\\}_{i=1}^{n}$. At time step $t=1$, decision
maker $\mathbbm{D}$ performs an intervention, changing the values
$\\{Z^{(i)}_{1}\\}_{i=1}^{n}$ to affect the distribution of outcomes $Y_{1}$
at time step $t=1$.

We define a _counterfactual treatment choice problem_ as a treatment choice
problem characterized by the following criteria:

  1. 1.

The sample is the entire population of interest.

  2. 2.

There are unobserved variables that are specific to each unit and that explain
some of the variation in outcomes. These unobserved variables also exhibit
some stability over time.

In other words, individuals $i=1,\ldots,n$ are the entire population of
interest, and exogenous variables $U_{Y}$ are not identical across units and
exhibit some stability over time (characterized more explicitly in Section 5).
Thus, in a counterfactual treatment choice problem, decision maker
$\mathbbm{D}$ focuses on the following counterfactual distribution:

$p_{Y_{1}}^{\mathfrak{C}\mid\mathcal{X};\text{do}(I)}\equiv\frac{1}{n}\sum_{i=1}^{n}p_{Y_{1}}^{\mathfrak{C}\mid\mathbf{V}=\mathbf{v}^{(i)};\text{do}(I)},$
(1)

with interventions of the form $I=\text{do}(Z_{1}^{(1:n)}=\vv{z})$, denoting a
vector of treatments across variables $Z^{(1)}_{1},\ldots,Z^{(n)}_{1}$.
$\mathbbm{D}$ wants to maximize a measure of welfare
$W(P_{Y_{1}}^{\mathfrak{C}\mid\mathcal{X};\text{do}(I)})$ as a function of
interventions $\vv{z}\in\bigtimes_{i=1}^{n}\Omega_{Z}$.

By contrast, in an _interventional treatment choice problem_ , we are instead
focused on the following distribution:

$p_{Y_{1}}^{\mathfrak{C};\text{do}(I)},$ (2)

where decision maker $\mathbbm{D}$ wants to maximize a measure of welfare
$W(P_{Y_{1}}^{\mathfrak{C};\text{do}(I)})$ in a target population different
from the current sample but whose data generating process follows the same SCM
$\mathfrak{C}$, thus sharing the same interventional distributions. Because
the concern is an interventional rather than counterfactual distribution in a
target population different from the current sample, the treatment assignment
in the interventional case is a decision problem based on covariate values $X$
rather than the specific units $\\{1,\ldots,n\\}$ in the sample. As in
Kitagawa and Tetenov (2018), this means that the set of treatment rules can be
indexed by their _decision sets_ $G\subset\mathcal{X}$ of covariate values,
which determine the group of individuals $\\{X\in G\\}$ to whom treatment 1 is
assigned.

In other words, in a counterfactual treatment choice problem, we do not assign
treatments for an individual based on their covariate values $X$, because we
are not learning a function that we can use to assign treatment to new
individuals, so we do not need to predict outcomes for unseen individuals. In
a decision problem over the sample, two individuals with the same covariate
values $X$ can have different values for their outcomes $Y$ and thus might
receive different treatments. The decision sets then are the _units_ in the
sample rather than subsets of covariate values.

###  4.1 Sociological Implications for Fairness

A focus on counterfactuals rather than interventions also has significant
sociological implications for reasoning about fairness as well as social
categories like race and gender. In general, operationalizing social
categories in the context of causal modeling is a nuanced technical and
sociological problem; see Bynum, Loftus, and Stoyanovich (2021) and references
therein for a detailed discussion. Were we to use racial categories in an
interventional treatment choice problem, we would be deciding to treat future
individuals with unseen outcomes based on their racial categories. Kitagawa
and Tetenov (2019), for example, suggest legality issues in assigning
treatments in an econometric setting based on social categories such as race.
By contrast, in a counterfactual treatment choice problem, racial categories
might _reveal_ disparate outcomes, but it is the realized, factual, already-
observed disparate outcomes that drive our treatment assignment rather than
assumptions or predictions about future disparities based on racial
categories. Such a problem would aim to tackle “pre-existing bias” as
discussed by Friedman and Nissenbaum (1996). In short, a counterfactual
treatment choice problem focuses on whichever grouping of people needs the
most help as measured by their outcomes, regardless of what the grouping is.

###  4.2 Example: EWM as Interventional Treatment Choice

To illustrate this distinction concretely, we draw a direct contrast with
Kitagawa and Tetenov (2019) as a representative interventional treatment
choice problem, referred to as EWM. We will refer to the counterfactual
approach as CF. Consider the following SCM $\mathfrak{C}$.

$\displaystyle\begin{split}X&=U_{X},\\\ Z&=U_{Z},\\\ Y&=X+Z+U_{Y},\end{split}$
(3)

where $U_{X},U_{Z}\sim\text{Bern}(0.5)$ and $U_{Y}\sim\text{U}(\\{0,1,2\\})$
with four observations — Units 1, 2, 3, and 4 — whose observed $X,Z,Y$ values
are shown in Table 3, along with the posterior values of the exogenous
variables and predicted potential outcomes $Y(0),Y(1)$ for possible treatments
$Z=0$ and $Z=1$. Assume for simplicity in this comparison that $X,\mathbf{U}$
are constant over time, so we do not need to consider time steps (we will
relax this assumption in Section 5).

Unit | X | Z | Y | $U_{X}$ | $U_{Z}$ | $U_{Y}$ | $Y(0)$ | $Y(1)$
---|---|---|---|---|---|---|---|---
1 | 0 | 0 | 1 | 0 | 0 | 1 | 1 | 2
2 | 0 | 0 | 2 | 0 | 0 | 2 | 2 | 3
3 | 1 | 0 | 1 | 1 | 0 | 0 | 1 | 2
4 | 1 | 0 | 2 | 1 | 0 | 1 | 2 | 3
Table 1: Observed data ($X,Z,Y$), inferred exogenous variables
($U_{X},U_{Z},U_{Y}$), and predicted potential outcomes ($Y(0),Y(1)$)
corresponding to the SCM in Equation 5.

From the observed data, we infer a posterior
$P_{\mathbf{U}\mid\mathbf{X}=\mathbf{x}}^{\mathfrak{C}}$ and (through
‘abduction, action, prediction’) compute counterfactual outcomes for each
individual for all treatment values $z\in\\{0,1\\}$. In this example, the
exogenous variables can be uniquely reconstructed from the observed data, so
the corresponding counterfactual distributions are point masses.333Point mass
posteriors are inspired by Example 6.18 from Peters, Janzing, and Schölkopf
(2017). For example,
$p_{Y}^{\mathfrak{C}\mid\mathbf{V}=\mathbf{v}_{2};\text{do}(Z=0)}(y)$ is a
point mass at $y=2$. The two possible point masses for each unit correspond to
potential outcome values $Y(0)$ or $Y(1)$, shown as two columns in Table 3.

Assumption | Model |  Interpretation
---|---|---
(A1) Exogenous factors are constant over time. | $\sigma_{U}=0$ |  Among the relevant variables we have not measured, each unit looks exactly the same next year as it does this year.
(A2) Exogenous factors vary over time. | $\sigma_{U}>0$ |  Among the relevant variables we have not measured, each unit looks somewhat the same next year as it does this year. Similarities grow weaker with larger $\sigma_{U}$ values.
(A3) Exogenous factors exhibit unstructured variation. | $\sigma_{\mu}=0$ |  Among the relevant variables we have not measured, each unit looks the same as any other unit, apart from random variability with time.
(A4) Exogenous factors exhibit structured (unit-specific) variation. | $\sigma_{\mu}>0$ |  Among the relevant variables we have not measured, there are units that look unlike other units, in addition to random variability with time. Units look less like each other with larger $\sigma_{\mu}$.
Table 2: A table of different possible assumptions about the structure and
stability of exogenous factors that influence the outcome of interest. Each
assumption is stated, along with its corresponding interpretation and
parameterization in Equation 4.

Imagine that we can afford at most two interventions. Given only one binary
covariate $X\in\\{0,1\\}$, three possible EWM decision sets we might consider
are $G_{\emptyset}\equiv\emptyset$, $G_{0}\equiv\\{0\\}$, and
$G_{1}\equiv\\{1\\}$, where we treat no units, units with $x=0$, or units with
$x=1$ respectively. In other words, consider the set of feasible policies
$\mathcal{G}\equiv\\{G_{\emptyset},G_{0},G_{1}\\}$.

Equality-minded EWM aims to maximize a rank-dependent social welfare function
in the overall population using the sample analog social welfare after
treatment. Consider the standard Gini social welfare function

$W(P_{Y})=\int_{0}^{\infty}(1-P_{Y}(y))^{2}dy.$

Assume we can perfectly recover the post-intervention outcome distribution, so
we can in turn perfectly estimate the population social welfare function for
each feasible policy. From SCM $\mathfrak{C}$ above, we can obtain the post-
treatment cumulative distribution functions
$P_{G_{\emptyset}}(y),P_{G_{0}}(y),P_{G_{1}}(y)$ for each policy
$G\in\mathcal{G}$ and compute the corresponding post-treatment population
social welfare values: $W(P_{G_{\emptyset}})=\frac{35}{36}$;
$W(P_{G_{0}})=\frac{56}{36}$; and $W(P_{G_{1}})=\frac{46}{36}$.444See
Technical Appendix A for a more detailed exposition of this example with
additional calculations shown. Thus, our social-welfare-maximizing treatment
policy from EWM is

$G^{*}_{\textsf{EWM}}\equiv\underset{G\in\mathcal{G}}{\text{sup
}}W(P_{G})=G_{0}.$

Now consider this problem from a perspective focused instead on the sample at
hand. We are now interested in the potential distributions of outcome $Y$
after different sets of interventions across _these four individuals_ rather
than the wider population. We can think of the distribution (probability mass
function) of outcome $Y$ after intervention on each of these units as an
average of the four corresponding point masses, given by Equation 1. Our
optimal treatment rule according to EWM, treatment rule $G_{0}$, would treat
Unit 1 and Unit 2, giving us cumulative distribution function
$P_{Y}^{\mathfrak{C}\mid\mathcal{X};\text{do}(Z^{(1:4)}=[1,1,0,0])}$ and post-
intervention social welfare in the sample
$W(P_{Y}^{\mathfrak{C}\mid\mathcal{X};\text{do}(Z^{(1:4)}=[1,1,0,0])})=\frac{26}{16}$.

Now consider instead treating Unit 1 and Unit 3. Observe that this treatment
rule does not determine treatment assignment for individuals by their
covariate values $X$ only (under an EWM treatment rule, we could not treat
both Units 1 and 3 and satisfy our budget constraint). Instead, our treatment
rule is based on selecting specific units driven by their realized outcome
values. Observe also that this treatment rule has the same ‘cost’ of treating
only two units. Under this intervention, we obtain cumulative distribution
function $P_{Y}^{\mathfrak{C}\mid\mathcal{X};\text{do}(Z^{(1:4)}=[1,0,1,0])}$
and post-intervention social welfare in the sample
$W(P_{Y}^{\mathfrak{C}\mid\mathcal{X};\text{do}(Z^{(1:4)}=[1,0,1,0])})=\frac{32}{16}$.

This is one of the possible treatment assignments we would consider from a
counterfactual (CF) perspective. Even without enumerating other possible
assignments, we can conclude that $G^{*}_{\textsf{CF}}\neq
G^{*}_{\textsf{EWM}}$, and, of more significance,
$W(G^{*}_{\textsf{CF}})>W(G^{*}_{\textsf{EWM}})$, so EWM and CF lead to
different optimal policies (as well as different candidate policies), and the
EWM-optimal policy is not necessarily the best choice for the welfare of the
given sample.

In summary, if we focus on the whole population, as in EWM, then we should
treat all units with $X=0$. But if we instead focus on our particular sample
of four units, as in CF, then we should treat Units 1 and 3, whose covariates
$X$ are different from each other and identical to those of other untreated
units.

##  5 Assumptions About Exogenous Variables

(a)

(b)

Figure 1: (a) Interventional, forward-looking counterfactual, and true
densities for $Y_{1}$ after intervention on $Z_{1}$. (b) Scatter-plots of
($Y^{(i)}_{0}$, $Y^{(i)}_{1}$) for each individual $i$ with the line
$Y_{0}=Y_{1}$ overlaid in gray.

In this work, we do not propose a catch-all characterization of exogenous
stability, and we do not propose a structure that applies to all models.
Instead, we focus on the following illustrative parameterization in order to
demonstrate insights that are applicable more generally.

Consider an outcome $Y$ and intervention $Z$ at two time steps (years), $t=0$
and $t=1$. We observe data across $n$ individuals
$\\{Z^{(i)}_{0},Y^{(i)}_{0}\\}_{i=1}^{n}$ in order to learn the relationship
between $Z$ and $Y$, and we intervene at time step $t=1$, determining the
value $Z_{1}$ in order to affect the distribution of $Y_{1}$. Intervention on
individual $i$ with value $Z^{(i)}_{0}$ increases their $Z$ value by constant
amount $\delta$, i.e., $Z^{(i)}_{1}=Z^{(i)}_{0}+\delta\cdot w(i)$, with
$w(i)\in\\{0,1\\}$ indicating the treatment decision for individual $i$.

#### Ground Truth.

The true relationship between $Z$ and $Y$ at each time step is as follows.

$\displaystyle Z^{(i)}_{0}$ $\displaystyle=U^{(i)}_{Z}$ (4) $\displaystyle
Y^{(i)}_{0}$ $\displaystyle=Z^{(i)}_{0}+U^{(i)}_{0}$ $\displaystyle
Z^{(i)}_{1}$ $\displaystyle=Z^{(i)}_{0}+\delta\cdot w(i)$ $\displaystyle
Y^{(i)}_{1}$ $\displaystyle=Z^{(i)}_{1}+U^{(i)}_{1}$

where $U^{(i)}_{Z}\sim\mathcal{N}(\mu_{Z},\sigma^{2}_{Z})$ and
$U^{(i)}_{0},U^{(i)}_{1}\overset{iid}{\sim}\mathcal{N}(\mu^{(i)}_{U},\sigma^{2}_{U})$,
with $\mu^{(i)}_{U}\sim\mathcal{N}(0,\sigma^{2}_{\mu})$. With this setup, we
can explore assumptions about exogenous stability and structure through
parameters $\sigma_{U}$ and $\sigma_{\mu}$. Table 2 shows 4 possible
assumptions and their parameterization in this example using $\sigma_{U}$ and
$\sigma_{\mu}$. Assumptions (A1) and (A2) are about stability (how a unit
looks exogenously compared to itself over time), while Assumptions (A4) and
(A5) are instead about structure (how a unit looks exogenously compared to
other units).

#### Model Perspective.

From the perspective of the modeler, who sees only
$\\{Z^{(i)}_{0},Y^{(i)}_{0}\\}_{i=1}^{n}$, a correct specification of the data
generating process at time $t=0$, given the available information, would be

$\displaystyle Z^{(i)}_{0}$ $\displaystyle=U^{(i)}_{Z}$ $\displaystyle
Y^{(i)}_{0}$ $\displaystyle=Z^{(i)}_{0}+U^{(i)}_{0}$

where $U^{(i)}_{Z}\sim\mathcal{N}(\mu_{Z},\sigma^{2}_{Z})$ and
$U^{(i)}_{0}\sim\mathcal{N}(0,\sigma^{2}_{\mu}+\sigma^{2}_{U})$. With this
model, the modeler would make the following interventional and forward-looking
counterfactual estimates of $P_{Y_{1}}$ after intervention $I$ that determines
treatment decisions $w(i)$.

#### Interventional Estimate of $P_{Y_{1}}$.

For an interventional estimate of $P_{Y_{1}}$, we would draw $U_{1}$ as a new
value from $\mathcal{N}(0,\sigma^{2}_{\mu}+\sigma^{2}_{U})$, i.e.,

$\displaystyle Z^{(i)}_{1}$ $\displaystyle=Z^{(i)}_{0}+\delta\cdot w(i)$
$\displaystyle U^{\prime(i)}_{1}$
$\displaystyle\sim\mathcal{N}(0,\sigma^{2}_{\mu}+\sigma^{2}_{U})$
$\displaystyle Y^{(i)}_{1}$ $\displaystyle=Z^{(i)}_{1}+U^{\prime(i)}_{1}$

with entailed distribution $P_{Y_{1}}^{\mathfrak{C};\text{do}(I)}$.

#### Forward-looking Counterfactual Estimate of $P_{Y_{1}}$.

A forward-looking counterfactual estimate of $P_{Y_{1}}$, rather than drawing
a new value of $U$, would instead set $U_{1}=U_{0}$ during abduction:

$\displaystyle Z^{(i)}_{1}$ $\displaystyle=Z^{(i)}_{0}+\delta\cdot w(i)$
$\displaystyle\tilde{U}^{(i)}_{1}$ $\displaystyle=U^{(i)}_{0}$ $\displaystyle
Y^{(i)}_{1}$ $\displaystyle=Z^{(i)}_{1}+\tilde{U}^{(i)}_{1}$

with entailed distribution
$P_{Y_{1}}^{\mathfrak{C}\mid\mathcal{X};\text{do}(I)}$.

###  5.1 Empirical Results

Figure 1(a) compares the true post-intervention density $p_{Y_{1}}$ to the
interventional estimate $p_{Y_{1}}^{\mathfrak{C};\text{do}(I)}$ and the
counterfactual estimate $p_{Y_{1}}^{\mathfrak{C}\mid\mathcal{X};\text{do}(I)}$
for $n=1000$, $\mu_{Z}=0$, $\sigma_{Z}=1$, $\delta=1$, and varying values of
$\sigma_{U},\sigma_{\mu}\in\\{0,0.5,5\\}$. We consider a setting where higher
values of outcome $Y$ are more desirable. The treatment rule is to treat all
individuals with outcomes $Y_{0}<0$.

If, among unobserved factors influencing $Y$, each unit looks exactly the same
next year as it does this year ($\sigma_{U}=0$), and each unit looks exactly
like every other unit ($\sigma_{\mu}=0$), then the interventional,
counterfactual, and true distribution for $Y_{1}$ all line up. The
corresponding scatter plot in Figure 1(b) shows that, for those with
$Y_{0}<0$, $Y_{1}$ is exactly $Y_{0}+\delta$, and for those untreated with
$Y_{0}>0$, $Y_{1}$ is exactly $Y_{0}$.

If each unit, among the exogenous factors influencing $Y$, looks exactly like
itself next year ($\sigma_{U}=0$), and units’ exogenous variables are not
exactly the same ($\sigma_{\mu}>0$), then the counterfactual distribution
lines up with the distribution of $Y_{1}$, but the interventional distribution
does not. In these cases in Figure 1(b), though the range of $Y_{0}$ values is
wider, $Y^{(i)}_{1}$ is again exactly $Y^{(i)}_{0}+\delta\cdot w(i)$.

Figure 2: Estimated Kullback–Leibler divergence values between true densities
for $Y_{1}$ and both interventional and forward-looking counterfactual
densities with individual treatment effect $\delta=1$, shown for
$\sigma_{U}\in[0,5]$ in increments of 0.1 across a few settings of
$\sigma_{\mu}\in\\{0,0.5,1,5\\}$. Figure 3: Estimated Kullback–Leibler
divergence values between true densities for $Y_{1}$ and both interventional
and forward-looking counterfactual densities with individual treatment effect
$\delta=5$, shown for $\sigma_{U}\in[0,12]$ in increments of 0.5 across a few
settings of $\sigma_{\mu}\in\\{0,0.5,1,5\\}$.

If, instead, among unobserved factors influencing $Y$, each unit will vary
next year compared to this year ($\sigma_{U}>0$), and each unit looks exactly
like every other unit apart from this random variation with time
($\sigma_{\mu}$ = 0), then the true distribution for $Y_{1}$ can lie between
the counterfactual and interventional distributions. For some values of
$\sigma_{U}>0$, the true distribution for $Y_{1}$ is visibly closer to the
interventional distribution (e.g., $\sigma_{U}=0.5$). For large enough values
of $\sigma_{U}$, the counterfactual distribution is again closer to the true
distribution of $Y_{1}$ (e.g., $\sigma_{U}=5$). In the corresponding row in
Figure 1(b), as $\sigma_{U}$ increases, for any given value of $Y_{0}$, the
value of $Y_{1}$ starts to vary enough to drown out the effect of the
intervention. As variation in exogenous factors across units increases
(increasing $\sigma_{\mu}$), forward-looking counterfactual estimate
$p_{Y_{1}}^{\mathfrak{C}\mid\mathcal{X};\text{do}(I)}$ matches $p_{Y_{1}}$
more closely than does $p_{Y_{1}}^{\mathfrak{C};\text{do}(I)}$.

Figure 2 gives a higher granularity view of changes in
$p_{Y_{1}}^{\mathfrak{C}\mid\mathcal{X};\text{do}(I)}$ and
$p_{Y_{1}}^{\mathfrak{C};\text{do}(I)}$ relative to $p_{Y_{1}}$ as
$\sigma_{U}$ increases for a few different values of $\sigma_{\mu}$.
Kullback–Leibler divergence values
$D_{\text{KL}}(p_{Y_{1}}||p_{Y_{1}}^{\mathfrak{C}\mid\mathcal{X};\text{do}(I)})$
and $D_{\text{KL}}(p_{Y_{1}}||p_{Y_{1}}^{\mathfrak{C};\text{do}(I)})$ are
estimated using the FNN package in R (Beygelzimer et al. 2022) to quantify how
close the interventional and counterfactual distributions are to the ground
truth. There are regions, when $\sigma_{U}$ is smaller and $\sigma_{\mu}$ is
sufficiently small, in which, by this measure,
$p_{Y_{1}}^{\mathfrak{C};\text{do}(I)}$ is a better estimate of $p_{Y_{1}}$
than $p_{Y_{1}}^{\mathfrak{C}\mid\mathcal{X};\text{do}(I)}$. We can focus on
the behavior in these regions by amplifying the individual treatment effect
$\delta$. Figure 3 shows the same information as Figure 2, now with $\delta=5$
instead of $\delta=1$ and for a correspondingly wider range of $\sigma_{U}$
values.

In summary, the two settings (1) where exogenous factors are sufficiently
stable over time and (2) where units’ exogenous factors are sufficiently
dissimilar (regardless of stability over time) are both appropriate for
forward-looking counterfactuals.

##  6 Discussion and Future Work

A standard story about counterfactuals is that we envision an unreal past in
order to understand the present. In this paper, we have explored how and when
we could alternatively use counterfactuals to envision an unreal past in order
to draw conclusions about our future. We introduced assumptions about the
stability and structure of exogenous variables that lead us to this forward-
looking rather than retrospective use of counterfactuals.

We also introduced counterfactual treatment choice problems as a setting that
motivates the use of forward-looking counterfactuals and studied empirically
the behavior of forward-looking counterfactual distributions as well as their
interventional counterparts under different exogenous conditions, showing that
each approach leads to different estimates.

What happens when decision maker $\mathbbm{D}$ makes a decision based on
assumptions consistent with the interventional case but exogenous variables
are actually stable, or vice-versa? A key implication is that these
distributional differences can be crucial for downstream policy decisions. To
illustrate this point, let us return to the case in Figure 1, where
$\sigma_{U}=5$, $\sigma_{\mu}=5$, and $\delta=1$. The variance $\mathbb{V}$
for each distribution is shown in the following table:

$\mathbb{V}[P_{Y_{0}}]$ | $\mathbb{V}[P_{Y_{1}}]$ | $\mathbb{V}[P_{Y_{1}}^{\mathfrak{C}\mid\mathcal{X};\text{do}(I)}]$ | $\mathbb{V}[P_{Y_{1}}^{\mathfrak{C};\text{do}(I)}]$
---|---|---|---
12.2 | 9.36 | 9.61 | 51.1

If, for example, our aim is to allocate treatments in order to decrease
variance in $Y$, the post-intervention variance we get from an interventional
estimate (51.1) is an order of magnitude larger than the actual value (9.36).
With this approach, we would expect an intervention for individuals with
$Y_{0}<0$ to increase variance overall, when in fact the opposite is true.
Conversely, in an unstable setting, exogenous noise could cause more
variability in actual outcomes than a forward-looking counterfactual approach
would predict, leading $\mathbbm{D}$ to disappointment if they expect their
chosen interventions to reduce $\mathbb{V}$.

Our work here also suggests several avenues for future exploration, including
extending results to more than two time steps with deeper connections to
latent variable and time series cross-sectional models, where exogenous
components of variables could be time-varying in different _estimable_ ways.
For example, mixed-effects models are closely related in settings where
estimation is possible. Adding additional generality to our characterization
of exogenous stability and structure in Table 2 could enable tests for
stability in different settings, as well as support reasoning about exogeneity
for multiple variables in arbitrary DAGs. Counterfactual treatment choice
problems also have their own set of technical challenges to investigate. For
example, if we can only reliably estimate intervals for counterfactuals (Lei
and Candès 2021), that would induce a policy choice problem with intervals
rather than point estimates.

Lastly, our exploration here of how mismatched assumptions about widely-used
inferential objects (like interventional distributions) can lead to
counterintuitive policy results suggests the creation of educational material
in these areas (e.g., Bynum et al. (2022)) as well as cross-collaboration with
decision and policy makers as important additional areas of future work.

## References

  * Athey and Imbens (2017) Athey, S.; and Imbens, G. W. 2017.  The State of Applied Econometrics: Causality and Policy Evaluation.  _Journal of Economic Perspectives_ , 31(2): 3–32.
  * Barocas, Hardt, and Narayanan (2019) Barocas, S.; Hardt, M.; and Narayanan, A. 2019.  _Fairness and Machine Learning_.  http://www.fairmlbook.org.
  * Beygelzimer et al. (2022) Beygelzimer, A.; Kakadet, S.; Langford, J.; Arya, S.; Mount, D.; and Li, S. 2022\.  _FNN: Fast Nearest Neighbor Search Algorithms and Applications_.  R package version 1.1.3.1.
  * Buesing et al. (2019) Buesing, L.; Weber, T.; Zwols, Y.; Racanière, S.; Guez, A.; Lespiau, J.-B.; and Heess, N. M. O. 2019.  Woulda, Coulda, Shoulda: Counterfactually-Guided Policy Search.  _ArXiv_ , abs/1811.06272.
  * Bynum et al. (2022) Bynum, L. E.; Khan, F. A.; Konopatska, O.; Loftus, J. R.; and Stoyanovich, J. 2022\.  An Interactive Introduction to Causal Inference.  _VISxAI: Workshop on Visualization for AI Explainability_.  https://lbynum.github.io/interactive-causal-inference.
  * Bynum, Loftus, and Stoyanovich (2021) Bynum, L. E.; Loftus, J. R.; and Stoyanovich, J. 2021.  Disaggregated Interventions to Reduce Inequality.  _Equity and Access in Algorithms, Mechanisms, and Optimization_.
  * Friedman and Nissenbaum (1996) Friedman, B.; and Nissenbaum, H. 1996.  Bias in Computer Systems.  _ACM Transactions on Information Systems_ , 14: 330–347.
  * Green and Chen (2019) Green, B.; and Chen, Y. 2019.  Disparate Interactions: An Algorithm-in-the-Loop Analysis of Fairness in Risk Assessments.  In _Proceedings of the Conference on Fairness, Accountability, and Transparency_ , FAT* ’19, 90–99. New York, NY, USA: Association for Computing Machinery.  ISBN 9781450361255.
  * Heidari et al. (2019a) Heidari, H.; Ferrari, C.; Gummadi, K. P.; and Krause, A. 2019a.  Fairness Behind a Veil of Ignorance: A Welfare Analysis for Automated Decision Making.  arXiv:1806.04959.
  * Heidari et al. (2019b) Heidari, H.; Loi, M.; Gummadi, K. P.; and Krause, A. 2019b.  A Moral Framework for Understanding Fair ML through Economic Models of Equality of Opportunity.  In _Proceedings of the Conference on Fairness, Accountability, and Transparency_ , FAT* ’19, 181–190. New York, NY, USA: Association for Computing Machinery.  ISBN 9781450361255.
  * Hu and Chen (2020) Hu, L.; and Chen, Y. 2020.  Fair Classification and Social Welfare.  In _Proceedings of the 2020 Conference on Fairness, Accountability, and Transparency_ , FAT* ’20, 535–545. New York, NY, USA: Association for Computing Machinery.  ISBN 9781450369367.
  * Kannan, Roth, and Ziani (2019) Kannan, S.; Roth, A.; and Ziani, J. 2019.  Downstream Effects of Affirmative Action.  In _Proceedings of the Conference on Fairness, Accountability, and Transparency_ , FAT* ’19, 240–248. New York, NY, USA: Association for Computing Machinery.  ISBN 9781450361255.
  * Kasy (2016) Kasy, M. 2016.  Partial Identification, Distributional Preferences, and the Welfare Ranking of Policies.  _Review of Economics and Statistics_ , 98: 111–131.
  * Kasy and Abebe (2021) Kasy, M.; and Abebe, R. 2021.  Fairness, Equality, and Power in Algorithmic Decision-Making.  _Proceedings of the 2021 ACM Conference on Fairness, Accountability, and Transparency_.
  * Khan, Manis, and Stoyanovich (2021) Khan, F. A.; Manis, E.; and Stoyanovich, J. 2021.  Fairness as Equality of Opportunity: Normative Guidance from Political Philosophy.  _CoRR_ , abs/2106.08259.
  * Kitagawa and Tetenov (2018) Kitagawa, T.; and Tetenov, A. 2018.  Who should be Treated? Empirical Welfare Maximization Methods for Treatment Choice.  _Econometrica_ , 86 2: 591–616.
  * Kitagawa and Tetenov (2019) Kitagawa, T.; and Tetenov, A. 2019.  Equality-Minded Treatment Choice.  _Journal of Business & Economic Statistics_, 39: 561 – 574.
  * Kusner et al. (2019) Kusner, M.; Russell, C.; Loftus, J.; and Silva, R. 2019.  Making Decisions that Reduce Discriminatory Impacts.  In Chaudhuri, K.; and Salakhutdinov, R., eds., _Proceedings of the 36th International Conference on Machine Learning_ , volume 97 of _Proceedings of Machine Learning Research_ , 3591–3600. San Diego, California: PMLR.
  * Lei and Candès (2021) Lei, L.; and Candès, E. J. 2021.  Conformal Inference of Counterfactuals and Individual Treatment Effects.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_.
  * Liu et al. (2018) Liu, L. T.; Dean, S.; Rolf, E.; Simchowitz, M.; and Hardt, M. 2018.  Delayed Impact of Fair Machine Learning.  In Dy, J.; and Krause, A., eds., _Proceedings of the 35th International Conference on Machine Learning_ , volume 80 of _Proceedings of Machine Learning Research_ , 3150–3158. PMLR.
  * Madras, Pitassi, and Zemel (2018) Madras, D.; Pitassi, T.; and Zemel, R. 2018.  Predict Responsibly: Improving Fairness and Accuracy by Learning to Defer.  In Bengio, S.; Wallach, H.; Larochelle, H.; Grauman, K.; Cesa-Bianchi, N.; and Garnett, R., eds., _Advances in Neural Information Processing Systems_ , volume 31. Curran Associates, Inc.
  * Manski (2003) Manski, C. F. 2003.  Statistical Treatment Rules for Heterogeneous Populations.  _Econometrica_ , 72: 1221–1246.
  * Mullainathan (2018) Mullainathan, S. 2018.  Algorithmic Fairness and the Social Welfare Function.  In _Proceedings of the 2018 ACM Conference on Economics and Computation_ , EC ’18, 1. New York, NY, USA: Association for Computing Machinery.  ISBN 9781450358293.
  * Nabi, Malinsky, and Shpitser (2019) Nabi, R.; Malinsky, D.; and Shpitser, I. 2019.  Learning Optimal Fair Policies.  _Proceedings of machine learning research_ , 97: 4674–4682.
  * Oberst and Sontag (2019) Oberst, M.; and Sontag, D. A. 2019.  Counterfactual Off-Policy Evaluation with Gumbel-Max Structural Causal Models.  _ArXiv_ , abs/1905.05824.
  * Pearl (2009) Pearl, J. 2009.  _Causality_.  Cambridge University Press.
  * Peters, Janzing, and Schölkopf (2017) Peters, J.; Janzing, D.; and Schölkopf, B. 2017.  _Elements of Causal Inference: Foundations and Learning Algorithms_.  MIT Press.
  * Spirtes et al. (2000) Spirtes, P.; Glymour, C. N.; Scheines, R.; and Heckerman, D. 2000.  _Causation, Prediction, and Search_.  MIT Press.

##  7 Acknowledgments

This work is supported in part by NSF Awards No. 1922658, 1934464, and
1916505. We are grateful to Rajesh Ranganath for clarifying discussions.

##  Appendix A EWM as Interventional Treatment Choice

Consider the following SCM $\mathfrak{C}$.

$\displaystyle\begin{split}X&=U_{X},\\\ Z&=U_{Z},\\\ Y&=X+Z+U_{Y},\end{split}$
(5)

where $U_{X},U_{Z}\sim\text{Bern}(0.5)$ and $U_{Y}\sim\text{U}(\\{0,1,2\\})$
with four observations — Units 1, 2, 3, and 4 — whose observed $X,Z,Y$ values
are shown in Table 3, along with the posterior values of the exogenous
variables and predicted potential outcomes $Y(0),Y(1)$ for possible treatments
$Z=0$ and $Z=1$. Assume for simplicity in this comparison that $X,\mathbf{U}$
are constant over time, so we do not need to consider time steps.

Unit | X | Z | Y | $U_{X}$ | $U_{Z}$ | $U_{Y}$ | $Y(0)$ | $Y(1)$
---|---|---|---|---|---|---|---|---
1 | 0 | 0 | 1 | 0 | 0 | 1 | 1 | 2
2 | 0 | 0 | 2 | 0 | 0 | 2 | 2 | 3
3 | 1 | 0 | 1 | 1 | 0 | 0 | 1 | 2
4 | 1 | 0 | 2 | 1 | 0 | 1 | 2 | 3
Table 3: Observed data ($X,Z,Y$), inferred exogenous variables
($U_{X},U_{Z},U_{Y}$), and predicted potential outcomes ($Y(0),Y(1)$)
corresponding to the SCM in Equation 5.

Denote the set of endogenous variables $\mathbf{V}=\\{X,Z,Y\\}$, with observed
data
$\mathcal{X}=\\{\mathbf{v}_{1},\mathbf{v}_{2},\mathbf{v}_{3},\mathbf{v}_{4}\\}$,
and the set of exogenous variables $\mathbf{U}=\\{U_{X},U_{Z},U_{Y}\\}$. From
the observed data, we infer a posterior $p_{\mathbf{U}}^{\mathfrak{C}}$ and
(through ‘abduction, action, prediction’) compute counterfactual outcomes for
each individual for all treatment values $z\in\\{0,1\\}$. In this example, the
exogenous variables can be uniquely reconstructed from the observed data, so
the corresponding counterfactual distributions are point masses.555Point mass
posteriors are inspired by Example 6.18 from Peters, Janzing, and Schölkopf
(2017). For example,

$p_{Y}^{\mathfrak{C}\mid\mathbf{V}=\mathbf{v}_{2};\text{do}(Z=0)}(y)=\begin{cases}1&y=2\\\
0&\text{otherwise}\end{cases}.$

The two possible point masses for each unit correspond to potential outcome
values $Y(0)$ or $Y(1)$, shown as two columns in Table 3.

Imagine that we can afford at most two interventions. Given only one binary
covariate $X\in\\{0,1\\}$, three possible EWM decision sets we might consider
are $G_{\emptyset}\equiv\emptyset$, $G_{0}\equiv\\{0\\}$, and
$G_{1}\equiv\\{1\\}$, where we can treat no units, units with $x=0$, or units
with $x=1$ respectively. In other words, consider the set of feasible policies
$\mathcal{G}\equiv\\{G_{\emptyset},G_{0},G_{1}\\}$.

Equality-minded EWM aims to maximize a rank-dependent social welfare function
in the overall population using the sample analog social welfare after
treatment. Consider the standard Gini social welfare function

$W(P_{Y})=\int_{0}^{\infty}(1-P_{Y}(y))^{2}dy.$

Assume we can perfectly recover the post-intervention outcome distribution, so
we can in turn perfectly estimate the population social welfare function for
each feasible policy. From SCM $\mathfrak{C}$ above, we have

$\displaystyle Y(1)\mid(X=1)$ $\displaystyle\sim\text{U}(\\{2,3,4\\})$
$\displaystyle Y(1)\mid(X=0)$ $\displaystyle\sim\text{U}(\\{1,2,3\\})$
$\displaystyle Y(0)\mid(X=1)$ $\displaystyle\sim\text{U}(\\{1,2,3\\})$
$\displaystyle Y(0)\mid(X=0)$ $\displaystyle\sim\text{U}(\\{0,1,2\\})$

from which we obtain the following post-treatment cumulative distribution
functions for each policy $G\in\mathcal{G}$ via

$\displaystyle P_{G}(y)$
$\displaystyle\equiv\int_{\mathcal{X}}\left[\vphantom{\int}P_{Y(0)\mid
X=x}(y)1\\{x\notin G\\}\right.$ $\displaystyle\left.+P_{Y(1)\mid
X=x}(y)1\\{x\in G\\}\vphantom{\int}\right]dP_{X}(x).$

Computing these distribution functions from the corresponding conditional
distributions for each feasible policy, we have

$P_{G_{\emptyset}}(y)=\begin{cases}0&y<0\\\ \frac{1}{6}&0\leq y<1\\\
\frac{3}{6}&1\leq y<2\\\ \frac{5}{6}&2\leq y<3\\\ 1&3\leq y\end{cases},\
P_{G_{0}}(y)=\begin{cases}0&y<1\\\ \frac{1}{3}&1\leq y<2\\\ \frac{2}{3}&2\leq
y<3\\\ 1&3\leq y\end{cases},$ $P_{G_{1}}(y)=\begin{cases}0&y<0\\\
\frac{1}{6}&0\leq y<1\\\ \frac{2}{6}&1\leq y<2\\\ \frac{4}{6}&2\leq y<3\\\
\frac{5}{6}&3\leq y<4\\\ 1&4\leq y\end{cases},$

and the corresponding post-treatment population social welfare values:
$W(P_{G_{\emptyset}})=\frac{35}{36}$; $W(P_{G_{0}})=\frac{56}{36}$; and
$W(P_{G_{1}})=\frac{46}{36}$. Thus our social-welfare-maximizing treatment
policy from EWM is

$G^{*}_{\textsf{EWM}}\equiv\underset{G\in\mathcal{G}}{\text{sup
}}W(P_{G})=G_{0}.$

Now consider this problem from a perspective focused instead on the sample at
hand. We are now interested in the potential distributions of outcome $Y$
after different sets of interventions across _these four individuals_ rather
than the wider population. We can think of the distribution (probability mass
function) of outcome $Y$ after intervention on each of these units as an
average of the four corresponding point masses. Our optimal treatment rule
according to EWM, treatment rule $G_{0}$, would treat Unit 1 and Unit 2,
giving us

$\displaystyle p_{Y}$
${}^{\mathfrak{C}\mid\mathcal{X};\text{do}(Z^{(1:4)}=[1,1,0,0])}$
$\displaystyle=\frac{1}{4}\sum_{\mathbf{v}_{i}\in\mathcal{X}}p_{Y}^{\mathfrak{C}\mid\mathbf{V}=\mathbf{v}_{i};\text{do}(Z^{(1:4)}=[1,1,0,0])}$
$\displaystyle=\begin{cases}\frac{1}{4}&y=1\\\ \frac{1}{2}&y=2\\\
\frac{1}{4}&y=3\\\ 0&\text{otherwise}\end{cases}$

with cumulative distribution function

$P_{Y}^{\mathfrak{C}\mid\mathcal{X};\text{do}(Z^{(1:4)}=[1,1,0,0])}=\begin{cases}0&y<1\\\
\frac{1}{4}&1\leq y<2\\\ \frac{3}{4}&2\leq y<3\\\ 1&3\leq y\end{cases}$

and post-intervention social welfare in the sample

$W(P_{Y}^{\mathfrak{C}\mid\mathcal{X};\text{do}(Z^{(1:4)}=[1,1,0,0])})=\frac{26}{16}.$

Now consider treating Unit 1 and Unit 3. Notice that this treatment rule does
not determine treatment assignment for individuals by their covariate values
$X$ only (under an EWM treatment rule, we could not treat both Units 1 and 3
and satisfy our budget constraint). Instead, our treatment rule is based on
selecting specific units driven by their realized outcome values. Notice also
that this treatment rule has the same ‘cost’ of treating only two units. Under
this intervention, we have

$\displaystyle p_{Y}$
${}^{\mathfrak{C}\mid\mathcal{X};\text{do}(Z^{(1:4)}=[1,0,1,0])}$
$\displaystyle=\frac{1}{4}\sum_{\mathbf{v}_{i}\in\mathcal{X}}p_{Y}^{\mathfrak{C}\mid\mathbf{V}=\mathbf{v}_{i};\text{do}(Z^{(1:4)}=[1,0,1,0])}$
$\displaystyle=\begin{cases}1&y=2\\\ 0&\text{otherwise}\end{cases}$

with cumulative distribution function

$P_{Y}^{\mathfrak{C}\mid\mathcal{X};\text{do}(Z^{(1:4)}=[1,0,1,0])}=\begin{cases}0&y<2\\\
1&2\leq y\end{cases}$

and post-intervention social welfare in the sample

$W(P_{Y}^{\mathfrak{C}\mid\mathcal{X};\text{do}(Z^{(1:4)}=[1,0,1,0])})=\frac{32}{16}.$

This is one of the possible treatment assignments we would consider from a CF
perspective. Even without enumerating other possible assignments, we can
conclude that $G^{*}_{\textsf{CF}}\neq G^{*}_{\textsf{EWM}}$ and of more
significance, $W(G^{*}_{\textsf{CF}})>W(G^{*}_{\textsf{EWM}})$, so EWM and CF
lead to different optimal policies (as well as different candidate policies),
and the EWM-optimal policy is not necessarily the best choice for the welfare
of the given sample.

In summary, if we focus on the whole population, as in EWM, then we should
treat all units with $X=0$. But if we instead focus on our particular sample
of four units, as in CF, then we should treat Units 1 and 3, whose covariates
$X$ are different from each other and identical to those of other untreated
units.