Finding Counterfactually Optimal Action Sequences in Continuous State Spaces
Authors: Stratis Tsirtsis
Manuel Gomez Rodriguez
Venue: NeurIPS 2023 poster (2023)
Keywords: Counterfactual reasoning
Markov decision process
Structural causal model
A* search
TL;DR: We propose a method to compute an action sequence that, in retrospect, would have led a particular observed episode of a continuous-state MDP to a better outcome.

Abstract: Whenever a clinician reflects on the efficacy of a sequence of treatment decisions for a patient, they may try to identify critical time steps where, had they made different decisions, the patient's health would have improved. While recent methods at the intersection of causal inference and reinforcement learning promise to aid human experts, as the clinician above, to *retrospectively* analyze sequential decision making processes, they have focused on environments with finitely many discrete states. However, in many practical applications, the state of the environment is inherently continuous in nature. In this paper, we aim to fill this gap. We start by formally characterizing a sequence of discrete actions and continuous states using finite horizon Markov decision processes and a broad class of bijective structural causal models. Building upon this characterization, we formalize the problem of finding counterfactually optimal action sequences and show that, in general, we cannot expect to solve it in polynomial time. Then, we develop a search method based on the A* algorithm that, under a natural form of Lipschitz continuity of the environment’s dynamics, is guaranteed to return the optimal solution to the problem. Experiments on real clinical data show that our method is very efficient in practice, and it has the potential to offer interesting insights for sequential decision making tasks.

Finding Counterfactually Optimal Action Sequences in Continuous State Spaces
Stratis Tsirtsis
Manuel Gomez-Rodriguez
Max Planck Institute for Software Systems Kaiserslautem, Germany stsirtsisdmpi-sws.
Max Planck Institute for Software Systems Kaiserslautern, Germany manuel gr grammi-rmssessee
Abstract
Whenever clinician reflects on the efficacy of sequence of treatment decisions for patient, they may try to identify critical time steps where, had they made different decisions, the patient's 's health would have improved. While recent methods at the intersection of causal inference and reinforcement learning promise to aid human experts, as the clinician above, to retrospectively analyze sequential decision making processes, they have focused on environments with finitely many discrete states. However, in many practical applications, the state of the environment is inherently continuous in nature. In this paper, we aim to fill this gap. We start by formally characterizing sequence of discrete actions and continuous states using finite horizon Markov decision processes and a broad class of bijective structural causal models. Building upon this characterization, we formalize the problem of finding counterfactually optimal action sequences and show that, in general, we cannot expect solve it in polynomial time. Then, we develop search method based on the A algorithm that, under a natural form of Lipschitz continuity of the environment's dynamics, is guaranteed to return the optimal solution the problem. Experiments on real clinical data show that our method is very efficient in practice, and it has the potential offer interesting insights for sequential decision making tasks.
Introduction
Had the chess player moved the king one round later, would they have avoided losing the game? Had the physician administered antibiotics one day earlier, would the patient have recovered? The process of mentally simulating alternative worlds where events of the past play out differently than they did in reality is known counterfactual reasoning []] Thoughts of this type of human decisions and they are tightly connected to the way we attribute causality and responsibility to events and others' actions [2]. The last decade has seen rapid development of reinforcement learning agents, presenting (close to) human-level performance making tasks, such as gaming B autonomous driving (5) and clinical decision support 66 7 In conjunction with the substantial progress made in the field of causal inference [8]) (9) this has led to a growing interest in machine learning methods that employ elements of counterfactual reasoning to improve or to retrospectively analyze decisions in sequential settings [1.155
common by-product
a variety of sequential decision
In the context of reinforcement learning, sequential decision making is typically modeled using Markov Decision Processes (MDPs) [6]] Here, we consider MDPs with a finite horizon where each episode (i.e., each sequence of decisions) consists of finite number of time steps. As example, consider a clinician treating a patient in an intensive care unit (ICU). At each time step, the clinician observes the current state of the environment (e.g., the patient's 'ss vital signs) and they choose among set of potential actions (e.g., standardized dosages of a drug). Consequently. the chosen action causes the environment transition (stochastically) into state, and the clinician earns reward (e.g..
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

satisfaction inversely proportional to the patient's 's severity). The process repeats until the horizon is met and the goal of the clinician is to maximize the total reward.
In this work, our goal is to aid the retrospective analysis of individual episodes as the example above. For each episode, aim to find an action sequence that differs slightly from the one taken in reality but, under the circumstances of that particular episode, would have led to a higher counterfactual reward. In our example above, assume that the patient's is condition does not improve after a certain period of time. A counterfactually optimal action sequence could highlight to the clinician a small set of time steps in the treatment process where, had they administered different drug dosages, the patient's s severity would have been lower. In turn, manual inspection of those time steps could provide insights to the clinician about potential ways to improve their treatment policy.
To infer how a particular episode would have evolved under a different action sequence than the one taken in reality, one needs to represent the stochastic state transitions of the environment using a structural causal model (SCM) 8 17]. This has been a key aspect of a line of work at the intersection of counterfactual reasoning and reinforcement learning, which has focused on methods to either design better policies using offline data [10 [2] or to retrospectively analyze individual episodes m 113. Therein, the work most closely related to ours is by Tsirtsis et al. [3, which
introduces a method to compute counterfactually optimal action sequences in MDPs with discrete states and actions using a Gumbel-Mas
SCM to model the environment dynamics [T]] However, in many practical applications, such as in critical care, the state of the environment is inherently continuous in nature [18]. In our work, we aim to fill this gap by designing method to compute counterfactually optimal action sequences in MDPs with continuous states and discrete actions. Refer to Appendix [A] for a discussion of further related work and to Pearl g for an overview of the broad field of causality.
Our contributions. We start by formally characterizing sequential decision making processes with continuous states and discrete actions using finite horizon MDPs and a general class of bijective SCMs 19. Notably, this class of SCMs includes multiple models introduced in the causal discovery literature 20 25). Building on this characterization, we make the following contributions:
We formalize the problem of finding counterfactually optimal action sequence for a particular episode in environments with continuous states under the constraint that it differs from the observed action sequence in at most A actions.
(ii) We show that the above problem is NP-hard using novel reduction from the classic parti- tion problem (26). This is in contrast with the computational complexity of the problem in environments with discrete states, which allows for polynomial time algorithms [13]. (iii) We develop search method based on the A" algorithm that, under a natural form of Lipschitz continuity of the environment's 's dynamics, is guaranteed to return the optimal solution to the problem upon termination.
Finally, we evaluate the performance and the qualitative insights of our method by performing a series of experiments using real patient data from critical care
2 A causal model of sequential decision making processes
At each time step & E [T' - 1] == {0, , 1}. where I' is a time horizon, the decision making process is characterized by D-dimensional vector state s; E S RD, an action a E A, where A is a finite set of N actions, and reward R(st, 04) E IR associated with each pair of states and actions. Moreover, given an episode of the decision making process,
{(*t, (t), [-1 the process's outcome o(T) = Σ, R(st, a4) is given by the sum of the rewards. In the remainder, we will denote the elements of a vector s as St,1,...,
Further, we characterize the dynamics of the decision making process using the framework of structural causal models (SCMs).
general, an SCM is consisted of four parts: (i) a set of endogenous variables (ii) a set of exogenous (noise) variables (iii) a set of structural equations assigning values to the endogenous variables, and (iv) a set of prior distributions characterizing the exogenous variables [8). In our setting, the endogenous variables of the SCM C are the random variables representing the states So.... ST and the actions Ag....,AT-I The action A, at time
I 'Our code is accessible at httt:://ghhhhh..oo/////os www/Networks-I www.eeetorrss--Lannnnooottteeeee Lamaiaaggaoannnfftttloottttnnnnn Taable] Li in Appendix B summarizes the notation used throughout the paper.

---

step t is chosen based on the observed state S and is given by structural (policy) equation
4 : 94 (S,, Z.),
(1)
where Z, Z is vector-valued noise variable, to allow some level of stochasticity in the choice of the action, and its prior distribution pe (Z₂) is characterized by a density function fz. Similarly, the state S1+1 in the next time step given by a structural (transition) equation
S g As, U.),
(2)
where U, E U is vector-valued noise variable with its prior distribution pc (Ut) having density function fü. and we refer to the function Is as the transition mechanism. Note that, in Eq. 2 the noise variables {U=) are mutually independent and, keeping the sequence of actions fixed, they are the only source of stochasticity in the dynamics of the environment. In other words, a sampled sequence of noise values {u=)r=o and fixed sequence of actions {a)) result into single (deterministic) sequence of states {st \T= This implicitly assumes that the state transitions are stationary and there are no unobserved confounders. Figure 4 Appendix B depicts the causal graph G corresponding to the SCM defined above.
The above representation of sequential decision making using an SCM C is a more general refor- mulation of Markov decision process, where a (stochastic) policy :(a 8) is entailed by Eq. and the transition distribution (i.e., the conditional distribution of S [.,, A,) is entailed by Eq. Specifically, the conditional density function of S4+1 St,, Ar is given by
p² (S = = = St, A₁ = az) = =p -piidoAA (S =s =s S. =8)
1 gs at, u)] fū (u)du,
(3)
Jusu
where do( (2 = ae) denotes a (hard) intervention on the variable A,, whose value is set to a Here. the first equality holds because S1+1 and At are d-separated by S; in the sub-graph obtained from G after removing all outgoing edges of A and the second equality follows from Eq. 2
Moreover, as argued elsewhere [1 13, by using an SCM to represent sequential decision making, instead of a classic MDP, we can answer counterfactual questions. More specifically, assume that, at time step t. observed the state S, St. took action A, a, and the next state was S = 8t+1. Retrospectively, we would like to know the probability that the state S would have been s' if, at time step we had been in a state and we had taken an action a, (generally) different from St, Og. Using the SCM C, we can characterize this by a counterfactual transition density function
oe 11 11111mm A mo do( (0 () (S = 8'|S = =)
1 = gs(s, a, u)] 1t. 51 (u)du, (4) Meu
where [ 1:.:1= is the posterior distribution of the noise variable U, with support such that 8c+1 = gs(Be; 0r; u).
In what follows, we will assume that the transition mechanism gs is continuous with respect to its last argument and the SCM C satisfies the following form of Lipschitz Lipschitzccontinnity:
Definition 1. An SCM C is Lipschitz- Lipschit iff the transition mechanism gs and the reward R are Lipschitz-comlimmons with respect to their first argument, i.e., for each a E A, a E U, there exists a Lipschitz constant Ka,u E RR_ such that, for any 8, s' &S.,Igst &, a, a,u.) - gs(s', a, u)ll S Ka ||s - s'1,, and, for each a E A, there exists a Lipschitz constant C. E <R44, such that, for any 8, 8' E S, (R(s,a) R(s', a)) E C. || s'l). In both cases, ||-|] denotes the Euclidean distance.
s'1l.
Note that, although they are not phrased in causal terms, similar Lipschitz continuity assumptions for the environment dynamics are common in prior work analyzing the theoretical guarantees of reinforcement learning algorithms (27 Moreover, for practical applications (e.g., in healthcare), this is a relatively mild assumption to make. Consider two patients whose vitals s and s' are similar
In general, the do operator also allows for soft interventions (i.e., setting probability distribution for As).
"Thi follows directly from the rules of do-calculus. For further details, refer to Chapter 3 of Pearl (8).
3

---

at certain point in time, they receive the same treatment and every unobserved factor u that may affect their health is also the same. Intuitively, Definition iimplies that their vitals will also evolve similarly in the immediate future, Le., the values (s, a, u) and gs a, u) will not differ dramatically. In this context, it is worth mentioning that, when the transition mechanism gs is modeled by a neural network, it is possible to control its Lipschitz constant during training, and penalizing high values can be seen as a regularization method [36 37].
Further, we will focus on bijective SCMs [9) a fairly broad class of SCMs, which subsumes multiple models studied in the causal discovery literature, such as additive noise models [20], post-nonlinear causal models (2I), location-scale noise models [22 and more complex models with neural network
components (2)).
Definition 2. An SCM C is bijective iff the transition mechanism gs is bijective with respect to its last argument, i.e., there is well-defiedd inverse function 9s" SS x A x S + U such that, for every ( Ct, Ut), it holds that u. = 9s (st, at 81+1)-
combination of St+1, St, at, Un with 8:41
Importantly, bijective SCMs allow for a more concise characterization of the counterfactual transition density given in Eq. 4 More specifically, after observing an event S1+1 S1+1: S. 8t, A, = Ut. the value u, of the noise variable U, can only be such that ut = 0s (8t, 0z, 8t+1), ie.,, the posterior distribution of U2 is a point mass and its density is given by
Jt.
el S: (u) = 1[u 9s (8t, at 81+1)],
(5)
where [[:] denotes the indicator function. Then, for given episode of the decision making process, we have that the (non-stationary) counterfactual transition density is given by
(S+11 = s'1 Ss 8, A₁ a) == p [S++1 |S++1=8441 88+1,,5 S.. ma do( Ar aa) (S =s'1 8)
[s = as(
a, u] ] 1 1/ue = 95 ( at, 11t1)]du
well
=1 [ = gs (s, 05 (st, ac, 81+1))]
Since this density is also a point mass, the resulting counterfactual dynamics are purely deterministic. That means, under a bijective SCM, the answer to the question "What would have been the state at time t +++ 1, had we been at state s and taken action a at time t, given that, in reality, we were at 8t, we took at and the environment transitioned to S1+1?" ?" is just given by s' gs (8,0.93 -1 (st, at, St-1)).
On the counterfactual identifiability of bijective SCMs. Very recently, Nasr-Esfahany and Kici- man (88 have shown that bijective SCMs are in general not counterfactually identifiable when the exogenous variable U, is multi-dimensional. In other words, even with access to an infinite amount of triplets (st, de, $t+1) sampled from the true SCM C, it is always possible to find SCM M * C with transition mechanism hi and distributions PM (U,) that entails the same transition distributions as C (i.e., it fits the observational data perfectly), but leads to different counterfactual predictions. Although our subsequent algorithmic results do not require the SCM C to be counterfactually identifi- able, the subclass of bijective SCMs we will use in our experiments in Section 5 is counterfactually identifiable. The defining attribute of this subclass, which we refer to as element- wise bijective SCMs, is that the transition mechanism Is can be decoupled into D independent mechanisms Is.; such that S.+ = gs,((S; (S,, Ar, Uz,)) for E {l,..., D}. This implies S 1 Ut,1] Uz.,, St, Ar for j 4 i, however, Ur. Ú.. do not need to be independent. Informally, we have the following identifiability result (refer to Appendix C for a formal version of the theorem along with its proof, which follows a similar reasoning to proofs found in related work [2 19):
Theorem (Informal). Ler C and M be two dlement-wise bijective SCMs such that their entailed transition distributions for S4n1 given any value of St,, A, are always identical. Then, all their counterfactual predictions based on an observed transition (8t, at, 8(+1) will also be identical.
On the assumption of no unobserved confounding. The assumption that there are no hidden confounders is a frequent assumption made by work at the intersection of counterfactual reasoning and reinforcement learning [3]] and, more broadly. in the causal inference literature (39 43 That said, there is growing interest in developing off- off-policy methods for partially observable MPDs (POMDPs) that are robust to certain types of confounding [44146), and in learning dynamic treatment regimes in sequential settings with mon-Markovian structure [47 48 Moreover, there is a line of work focusing on the identification of counterfactual quantities in non-sequential confounded environments [4). In that context, we consider the computation of (approximately) optimal counterfactual action sequences under confounding as very interesting direction for future work.

---

3 Problem statement
Let r an observed episode of decision making process whose dynamics are characterized by a Lipschitz-continouus bijective SCM. To characterize the counterfactual outcome that any alternative action sequence would have achieved under the circumstances of the particular episode, we build upon the formulation of Section 2 and we define a non-stationary counterfactual MDP = (S+, A,, F Rtt, (tt,7 T) with deterministic transitions. Here, S S x [T - 1] is an enhanced state space such that each st E St is a pair (s, I) indicating that the counterfactual episode would have been at state 8 E S with / action changes already performed. Accordingly, Rt is reward function which takes the form R+ (ts, I), a) = R(s, a) for all (s, E S a E A. i.e., it does not change depending the number of action changes already performed. Finally, the time-dependent transition function F S+ x A - S+ is defined as
-
P ((s,). = {(s (8, a.a,, (8t, Gt 81+1) ,1+1) if (a
[(s (8,0t, 0c.95 95 (8,, at 8(41)) .1) otherwise.
(7)
Intuitively, here the transition function according to the point mass of the counterfactual transition density given in Eq. 6 and we use the second coordinate to keep track of the changes that have been performed in comparison to the observed action sequence up the time step
Now, given the initial state 8o of the episode T and any counterfactual action sequence {a)) compute the corresponding counterfactual episode r' {($) 4), df Its sequence of states is given recursively by
4) = Fo ((+0, 0), d() and (8/41) = P, ((8) (4, a/) for E
(8)
and its counterfactual outcome is given by o+ (T') (T') E, R+ (($), (4, d() E, R (s),
Then, similarly as in Tsirtsis et al. [3]. our ultimate goal is to find the counterfactual action sequence [a) that, starting from the observed initial state su. maximizes the counterfactual outcome sub- ject to constraint the number of counterfactual actions that can differ from the observed ones, i.e.
T
maximize 4g s...
o o (T')
subject to s/ =0 and Σ 1aa + sk
1=0
where 0o,-... aT are the observed actions. Unfortunately, using a reduction from the classic partition problem (26), the following theorem shows that cannot hope to find the optimal action sequence in polynomial time.
Theorem 4. The problem defined by Eq. is NP-Hard.
The proof of the theorem relies on reduction from the partition problem (26), which is known to be NP-complete, to problem, defined in Eq. 9 At a high-level, we map any instance of the partition problem to an instance of our problem, taking special care to construct a reward function and an observed action sequence, such that the optimal counterfactual outcome ot (r*) takes specific value and only if there exists a valid partition for the original instance. The hardness result of Theorem motivates our subsequent focus the design of method that always finds the optimal solution to our problem at the expense of a potentially higher runtime for some problem instances.
4 Finding the optimal counterfactual action sequence via A search
To deal with the increased computational complexity of the problem, we develop an optimal search method based the classic A" algorithm [52). which we have found to be very efficient in practice.
Our starting point is the observation that, the problem of Eq. presents an optimal substructure, i.e. its optimal solution can be constructed by combining optimal solutions smaller sub-problems. For an observed episode т, let (= t) be the maximum counterfactual reward that could have been achieved in counterfactual episode where, at time the process is (counterfactual) state s, and there are so far / actions that have been different in comparison with the observed action sequence. Formally,
T-1 Σ R(s de)
(s, = max aly
subject to s and Σ laa + d' I.
'The supporting proof' of all Theorems, Lemmas and Propositions can be found in Appendix
5

---

% O
8,
R 1e. +1
Otserved action sequencc
R 1
(a) Scarch graph
(b) Hcuristic function computation
Figure 1: Main components of our search method based on the A algorithm. Panel (a) shows the search graph for a problem instance with AA| = 2. Here, each box represents a node the graph, and each edge represents a counterfactual transition. Next each edge, we include the action a E A causing the transition and the associated reward. Panel (b) shows the heuristic function computation, where the two axes represent a (continuous) state space S IR= and the two levels on the z-axis correspond to differences in the (integer) values (1, () and (la, + 1). Here, the blue squares correspond to the finite states in the anchor set S and (8a, = F7., ((s,1),a).
(s. 1, t) of
Then, it is easy to see that the quantity (=, can be given by the recursive function
(,,,, max {R(s. a) + V ta, + 1)}, for all s E S, << l and t < -,,, (10) aEA
where (Sa, F ((s, 1), ,a). In the base case of = k (i.e., all allowed action changes are already
performed), we have (s, k t) = R(s, a.) + V2 (Sa, la,, L + 1) for all E S and < T 1, and V= k. -1) R(s, aT-1) for 1 = T 1. Lastly, when t = - and < k, we have V= - = m R(s,a) for all e S.
Given the optimal substructure of the problem, one may be tempted to employ a typical dynamic programming approach to compute the values (=,,, f) in bottom-up fashion. However, the complexity of the problem lies in the fact that, the states 8 are real-valued vectors whose exact values depend the entire action sequence that led to them. Hence, to enumerate all the possible values that s might take, one has to enumerate all possible action sequences in the search space, which is problem with brute force search. In what follows,
equivalent to solving method to find optimal solutions using the A' algorithm, with the caveat that its runtime varies depending on the problem instance, and it can be equal to that of brute force search in the worst case.
present our proposed
Casting the problem as graph search. We represent the solution space of our problem as a graph, where each node corresponds to tuple (s, l,t) with s E S, E [k] and E [T]. Every node v = (s, I, t) with < k and < T 1 has |A| outgoing edges, each one associated with an action E A, carrying a reward R(s, a). and leading to a node va = (8a, la, + 1) such that (s., La) F+ ((s. ()s,,),=a .a). In the case of / = k: the node has exactly one edge corresponding to the observed action 0t at time f. Lastly, when
1, the outgoing edge(s) lead(s) to a common node uT (sp, k. I') which we call the goal node, and it has zero outgoing edges itself. Note that, the exact value of 88 is irrelevant, and we only include it for notational completeness.
Let 8n be the initial state of the observed episode. Then, it is easy to notice that, starting from the root node 0o = (so, O, 0), the first elements of each node 2 on path u.,..., Ui...., UT form a sequence of counterfactual states, and the edges that connect those nodes are such that the corresponding counterfactual action sequence differs from the observed one at most i actions. That said, the
counterfactual outcome (+( = T- H(s(, af) is expressed as the sum of the rewards associated with each edge in the path, and the problem defined by Eq. 9 is equivalent to finding the path of maximum total reward that starts from Up and ends in Ur. Figure Ta illustrates the search graph for simple instance of our problem. Unfortunately, since the states are vectors of real values, even enumerating all the graph's 's nodes requires time exponential in the number of actions |.A], which makes classic algorithms that search over the entire graph non-practical.
To address this challenge, we resort to the A algorithm, which performs a more efficient search over the graph by preferentially exploring only parts of it where we have prior information that they are

---

Algorithm 1: It computes upper bounds (s, for the values (s,
Input: States S, actions A, observed action sequence {0e)==0 1-T-1 horizoo T, transition function Fi, reward function R, constraint k anchor set Sr. Initialize: (s, T 1) + max.EA a)a
8 8 S 1= O, .,..,, -.
V, (s, k 1) H(s,at), s S.
for = T - .,....,0 do
for 1 = k ....,, 0 do
avoilable_ _cctiinssstssst + az if 1 == k else A
for E S do bounds 0
for E available octions do a) Va +- mi, (V. (st, la, + 1) bounds 4- bounds U (R(s,a) a) + V.} end
82, la + F7,,
/- Get the min bound for (a (t, 1)
T441 ||st sall3
based on S */
(s,, + max(bounds); -
/+ Get the max bound over the actions */
end end end
return (8,1,)) for all E St, E [k]] E [T -1]
more likely to lead to paths of higher total reward. Concretely, the algorithm proceeds iteratively and maintains a queue of nodes to visit, initialized to contain only the root node vo. Then, at each step, it selects one node from the queue, and it retrieves all its children nodes in the graph which are subsequently added to the queue. It terminates when the node being visited the goal node vT. Refer to Algorithm 2 in Appendix E for pseudocode implementation of the A algorithm.
The key element of the A algorithm is the criterion based on which it selects which node from the queue to visit next. Let e = (s;, 4, t) be a candidate node in the queue and Tu, be the total reward of the path that the algorithm has followed so far to reach from og to U/. Then, the A" algorithm visits next the node i; that maximizes the sum Top + (s,,, 4,,4, where V, is a heuristic function that aims to estimate the maximum reward that can be achieved via any path starting from vs (s., L; I) and ending in the goal node uT, i.e., gives an estimate for the quantity V- (8; li, l). Intuitively, the heuristic function can be thought of as "eye into the future" of the graph search, that guides the algorithm towards nodes that are more likely to lead to the optimal solution and the algorithm'ss 's performance depends on the quality of the approximation of (s; 4, f) by V. (s, 4. t). Next, we will look for heuristic function that satisfies consistency to guarantee that the A algorithm as described above returns the optimal solution upon termination [52].
Computing a consistent heuristic function. We first propose algorithm that computes the function's values (s, 1, t) for a finite set of points such that E [k], E [T [ 1]. S₁ is a pre-defined finite set of states-an an anchor set-whose construction we discuss later. Then, based on the upper bounds of the corresponding values V= 1, and we expand the definition of the heuristic function V- over all 8 E S by expressing it in terms of those upper bounds. Finally, we prove that the function resulting from the aforementioned procedure is consistent.
St C S, where
of the SCM C, we show that these computed values of V, are valid
To compute the upper bounds V/z, we exploit the observation that the values V=(,,,, satisfy form of Lipschitz-comtinuiyy as stated the following Lemma.
Lemma 5. Let u, = 9g (st, dt; 8t+1), 1), Ku, maxae A Kayuy C max C. and the sequence Lo....,Ly LT-1 E R2, be such that LT-1 = C and L = C + LL+1Ku, for t E [ 2]. Then, it V, (s', 1,t) S L; [|s - s'11, for all [T- 1], [k] and 8, s' S.
holds that (.
Based on this observation, our algorithm proceeds in a bottom-up fashion and computes valid upper bounds of the values V, (s, () for all / [k]] [T - 1] and s in the anchor set St. To get the
A heuristic function V, is consistent iff, for nodes v (,) % = (sa, la, 1) connected with an edge associated with action o, it satisfies (, 2 R(s, a) + V= la, +11 53
7

---

intuition, assume that, for given t, the values (= 1) are already computed for all & E St. I E [k] and they are indeed valid upper bounds of the corresponding V=(s, + 1). Then, let (a) = F7,, ((s,1), a) for some 8 E St and E [K]. Since Sa itself may not belong the finite anchor set S the algorithm uses the values (st, la, 1-41111 +11) of all anchors 8 E S in combination with their distance Sa. and it sets the value of V, (s, in way that it is also guaranteed to be (maximally tight) upper bound of V (8, 1, 1). Figure [b] iluutuusssssss the above operation. Algorithm summarizes the overall procedure, which is guaranteed to return upper bounds, as shown by the following proposition:
Proposition 6. For all S 1 [ [T -11, it holds that (s 2 V+(s,l,t), where (+ are the values of the heuristic function computed by Algorithm 7
Next, we use the values V (s, 1, computed by Algorithm Ito expand the definition of V, over the entire domain as follows. For some E S, a E A, let (8a, (a) = Fr. ([ a), then, we have that
0
t= t t -1
max R(s,a) aEA'
(,
(11)
(=, L141 ||s sull [|lll
max "EA'
R(s,a) + min of ES
otherwise,
where A'=
= {a=} for and A for k. Finally, the following theorem shows that the resulting heuristic function V, is consistent:
Theorem 7. For any nodes v = (8,1,), U2 = (Say lay + 1) with < T connected with an edge associated with action a, holds that V=(,, 2 R(s,a) + V- (Sa, la, + 1). Moreover, for any node I, 1) and edge connecting it to the goal node uT (sp. k, T), holds that (s, 1) 2 R(s,a) a) + (, T).
Kick-starting the heuristic function computation with Monte Carlo anchor sets. For any Si, whenever compute V, using Eq. 11 the resulting value is set based the value V=(8,,la, 1) of some anchor 8t. increased by a penalty term L++1 || - sall- Intuitively, this allows us to think of the heuristic function V, as an upper bound of the function V, whose looseness depends on the magnitude of the penalty terms encountered during the execution of Algorithm [ and each subsequent evaluation of Eq. 11 To speed up the A" algorithm, note that, ideally, one would want all penalty terms to be zero, i.e., an anchor set that includes all the states s of the nodes
= (s, that are going to appear in the search graph. However, as discussed in the beginning of Sec. an enumeration of those states requires a runtime exponential in the number of actions.
To address this issue, we introduce a Monte Carlo simulation technique that adds to the anchor set the observed states (so...., 8T-1) and all unique states {s/,..... s resulting by M randomly sampled counterfactual action sequences dbg- ......... .4. Specifically, for each action sequence, we first sample number A'' of actions to be changed and what those actions are going to be, both uniformly at random from (1...., k} and AA respectively. Then, we sample from {0, .... -......TT1. 1} the & time steps where the changes take place, with each time step having a probability L./ E.. Lu to be selected. This biases the sampling towards earlier time steps, where the penalty terms are larger due to the higher Lipschitz constants. As we will see in the next section, this approach works well in practice, and it allows to control the runtime of the A' algorithm by appropriately adjusting the number of samples M. We experiment with additional anchor set selection strategies in Appendix
5 Experiments using clinical sepsis management data
Experimental setup. To evaluate our method, we use real patient data from MIMIC-II
54 freely
accessible critical care dataset commonly used in reinforcement learning for healthcare 6 55 57].
We follow the preprocessing steps of Komorowski et al. [6] to identify a cohort of 20,926 patients treated for sepsis [58). Each patient record contains vital signs and administered treatment information in time steps of 4-hour intervals. As an additional preprocessing step, we discard patient records whose associated time horizon T is shorter than 10, resulting in a final dataset of 15,992 patients with horizons between 10 and 20.
8

---

2.5
50
2.2
30
2.1
EBB
anecagg
20
2.0
10
200
1.0
0.8
0.9
1.0
1.1
1.2
1.3
500 1000 1500 2000 2500 9000 Morre Carlo samplos, M
Tipechito cometand,
adtion dhangge, A0
(a) Efficiency vs. La
(b) Efficiency vs.
(c) Efficiency vs.
Figure 2: Computational efficiency of method under different configurations, measured by the effective branching factor (pink-left axis) and the runtime of the A algorithm (green-right axis). In Panel (a), we set M 2000 and 3. In Panel (b), we set T., 1.0 and k = 3. In Panel (c), set La = 1.0 and M = 2000. In all panels, we set Lg 0.1 and error bars indicate 95% confidence intervals over 200 executions of the A algorithm for 200 patients with horizon T 12.
To form state space 8 = RD we use D = 13 features. Four of these features are demographic contextual and thus we always set their counterfactual values to the observed ones. The remaining
δ features are time- varying and include the SOFA score [59] a standardized score of organ failure aate-along with eight vital signs that are required for its calculation. Since SOFA scores positively correlate with patient mortality [60], we assume that each s gives a reward R(s) equal to the negation of its SOFA value. Here, it is easy to see that this reward function is just a projection of s, fherefore, it is Lipschitz continuous with constant C. = for all a E A. Following related work [6] [53 [7]. consider action space A that consists of 25 actions, which correspond to 5 x 5 levels of administered vasopressors and intravenous fluids. Refer to Appendix G for additional details on the features and actions.
To model the transition dynamics of the time-narying features, we consider an SCM C whose transition mechanism takes a location-scale form gs (s, A,, U,) = h(S,, Ar) + o(Sr At) @ U.. where : S model is glement-wise bijective and hence it is counterfactually identifiable, as shown in Section 2 Moreover, we use neural networks to model the location and scale functions and and enforce their Lipschitz constants to be L,, and Lgs respectively. This results in a Lipschitz continuous SCM C with Ka,u = LA + Lg max; [u... Further, we assumee that the noise variable U, follows a multivariate Gaussian distribution with zero mean and allow its covariance matrix be (trainable) parameter.
R.D,. and @ denotes the element- wise multiplication (2 24 Notably. this
We jointly train the weights of the networks and 6 and the covariance matrix of the noise prior on the observed patient transitions using stochastic gradient descent with the negative lodelikd of each transition as a loss. our experiments, if not specified otherwise, we use an SCM with Lipschitz constants L2, = 1.0. Lg = 0.1 that achieves log-likeling lig-liielihood only 6% lower to that of the best model trained without any Lipschitz constraint. Refer to Appendix for additional details on the network architectures, the training procedure and the way we enforce Lipschitz continuity
Results. We start by evaluating the computational efficiency of our method against (i) the Lipschitz constant of the location network Ln. (ii) the number of Monte Carlo samples M used to generate the anchor set Sr, and (iii) the number of actions k that can differ from the observed ones. We measure efficiency using running time and the effective branching factor (EBF) [2]. The EBF is defined as real number 2 1 such that the number of nodes expanded by A equal to 1 + b b2 +...+ 6 where T is the horizon, and values close to 1 indicate that the heuristic function is the most efficient in guiding the search. Sigureejsmmmrizze the results, which show that our method maintains overall fairly low running time that decreases with the number of Monte Carlo samples M used for the generation of the anchor set and increases with the Lipschitz constant La, and the number of action changes k. That may not come as surprise since, as La increases, the heuristic function becomes more loose, and as & increases, the size of the search space increases exponentially. To put things in perspective, for a problem instance with LP. 1.0, k 3 and horizon T 12. the A" search led by our hcuristic function is effectively equivalent to an exhaustive search over a full tree with
All experiments were performed using internal cluster of machines equipped with 16 Intel(R) Xeon(R) 3. 3.20GH CPU cores, 512GBs of memory and 2 NVIDIA A40 48GB GPUs.

---

Acenge ef (injeoooenn
Nombhe of quliints
2030
100
15
10-
8%%
1000
SOFX
2076
055
3
066
1095
20%
3005
40%
10
20
30
40
3 of action changge,
Counteefactual Iregrovmeen
Time (bours)
(a) Average c/f improvement vs. (b) Distribution of ef improvement
(c) Observed vs. c/f episode
Figure 3: Retrospective analysis of patients' episodes. Panel (a) shows the average counterfactual improvement as a function of for a set of 200 patients with horizon T = 12, where error bars indicate 95% confidence intervals. Panel (b) shows the distribution of counterfactual improvement across all patients for A = 3, where the dashed vertical line indicates the median. Panel (c) shows the observed (solid) and counterfactual (dashed) SOFA score across time for a patient who presents 19.9% counterfactual improvement when A = 3. Upward (downward) arrows indicate action changes that suggest a higher (lower) dosage of vasopressors (V) and fluids (F). In all panels, we set M 2000.
2.1 7,355 leaves while the corresponding search space of our problem consists of more than million action sequences more than 3 million paths to reach from the root node to the goal node.
Next, investigate what extent the counterfactual action sequences generated by method would have led the patients in our dataset to better outcomes. For each patient, we measure their counterfactual improvement-thh relative decrease in cumulative SOFA score between the counter- factual and the observed episode. Figures [3a (3b] summarize the results, which show that: (i) the average counterfactual improvement shows a diminishing increase as k increases; (ii) the median counterfactual improvement is only 5%, indicating that, the treatment choices made by the clinicians for most of the patients were close to optimal, even with the benefit of hindsight; and (iii) there are 176 patients for whom our method suggests that a different sequence of actions would have led to an outcome that least 15% better. That said, we view patients at the tail of the distribution as "interesting cases" that should be deferred to domain experts for closer inspection, and we present one such example in Fig. Be| In this example, our method suggests that, had the patient received an early higher dosage of intravenous fluids while some of the later administered fluids where replaced by vasopressors, their SOFA score would have been lower across time. Although we present this case as purely anecdotal, the counterfactual episode is plausible, since there are indications of decreased mortality when intravenous fluids are administered at the early stages of a septic shock [6]]
6 Conclusions
In this paper, we have introduced the problem of finding counterfactually optimal action sequences in sequential decision making processes with continuous state dynamics. We showed that the problem is NP-hard and, to tackle it, we introduced a search method based on the A" algorithm that is guaranteed to find the optimal solution, with the caveat that its runtime instance. Lastly, using real clinical data, we have found that our method is very efficient in practice, and it has the potential to offer interesting insights to domain experts by highlighting episodes and time-steps of interest for further inspection.
vary depending on the problem
Our work opens up many interesting avenues for future work. For example, it would be interesting to develop algorithms with approximation guarantees that run in polynomial time, at the expense of not achieving strict counterfactual optimality. Moreover, since the practicality of methods like ours relies on the assumption that the SCM describing the environment is accurate, it would be interesting to develop methods to learn SCMs that align with human domain knowledge. Finally, it would be interesting to validate our method using real datasets from other applications and carry out user studies in which the counterfactual action sequences found by our method are systematically evaluated by the human experts (e.g., clinicians) who took the observed actions.
Acknowledgements. Tsirtsis and Gomez- Rodreg-Rdr acknowledge support from the European Re- search Council (ERC) under the Furopean Union's 's Horizon 2020 research and innovation programme (grant agreement No. 945719).
10

---

References
[1] Ruth MJ Byrne Counterfactual thought. Annual review of psychology, 67:135-157, 2016.
[2] David A Lagnado, Tobias Gerstenberg, and Ro'i Zultan. Causal responsibility and counterfactu- als. Cognitive science, (76)):036-1073, 2013.
[3] David Silver, Aja Huang, Chris Maddison, Arthur Guez, Laurent Sifre, George Van Den Driess- che, Julian Schrittwieser, loannis Antonoglou, Veda Panneershelvam, Marc Lanctot, et al. Mas- tering the game of go with deep neural networks and tree search. nature, 529(7587):484 2016.
[4] Lukasz Kaiser, Mohammad Babaeizadeh, Piotr Milos, Blažej Osiński, Roy H Campbell, Konrad Czechowski, Dumitru Erhan, Chelsea Finn, Piotr Kozakowski, Sergey Levine, Afroz Mohiuddin, Ryan Sepassi, George Tucker, and Henryk Michalewski. Model based reinforcement learning for atari. In International Conference on Learning Representations, 2020.
[5] B Ravi Kiran, Ibrahim Sobh, Victor Talpaert, Patrick Mannion, Ahmad A Al Sallab, Senthil Yogamani, and Patrick Pérez. Deep reinforcement learning for autonomous driving: A survey. IEEE Transactions on Intelligent Transportation Systems, 23(6):4909-499) 4926, 2021.
[6] Matthieu Komorowski, Leo Celi, Omar Badawi, Anthony C Gordon, and A Aldo Faisal. The artificial intelligence clinician learns optimal treatment strategies for sepsis in intensive care. Nature medicine, 24(11):1716-125) 2018.
[7] Chao Yu, Jiming Liu, Shamim Nemati, and Guosheng Yin. Reinforcement learning in healthcare: A survey, ACM Computing Surveys (CSUR), 55(1):1-36, 2021.
[8] Judea Pearl. Causality. Cambridge university press, 2009.
[9] Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Elements of causal inference: founda- tions and learning algorithms. The MIT Press, 2017.
[10] Lars Buesing, Theophane Weber, Yori Zwols, Nicolas Heess, Sebastien Racaniere, Arthur Guez, and Jean Baanttaptiste Lespiau. Woulda, coulda, shoulda: Counterfactually-gune policy search. In International Conference Learning Representations, 2019.
[11] Michael Oberst and David Sontag. Counterfactual off-policy evaluation with gumbel-max structural causal models. In International Conference on Machine Learning, pages 4881-488 PMIR, 2019.
[12] Chaochao Lu, Biwei Huang, Ke Wang. José Miguel Hernández- Kun Zhang, and Bernhard Schölkopf. Sample-efficient reinforcement learning via counterfactual-be based data augmentation. arXiv preprint arXiv:2012.09092, 2020.
[13] Stratis Tsirtsis. Abir De, and Manuel Gomez-Rodrigee Counterfactual explanations in sequential decision making under uncertainty. Advances in Neural Information Processing Systems, 34:30127-30133 2021.
[14] Jonathan Richens, Rory Beard, and Daniel H Thompson. Counterfactual harm. In Advances in Neural Information Processing Systems, volume 35, pages 36350-36365, 2022.
[15] Kimia Noorbakhsh and Manuel Gomez-Rodriguez. Counterfactual temporal point processes. Advances Neural Information Processing Systems, 0124810-24822, 2022.
[16] Richard Sutton and Andrew Barto. Reinforcement learning: An introduction. MIT press, 2018.
[17] Elias Bareinboim, Juan D Correa, Duligur Theling, and Thomas leard. On pearl's hierarchy and the foundations of causal inference. In Probabilistic and causal inference: the works of judea pearl, pages 507-556. 2022.
[18] Malcolm Elliott and Alysia Coventry. Critical care: the eight vital signs of patient monitoring. British Journal of Nursing, 21(10):621-62, 2012.

---

[19] Arash Nasr-Eaae Fasr-Esfahany, Mohammad Alizadeh, and Devavrat Shah. Counterfactual identifiability of bijective causal models. In Proceedings of the 40th International Conference on Machine Learning, ICMI'23. JorR.org, 2023.
[20] Patrik Hoyer, Dominik Janzing, Joris M Mooij. Jonas Peters, and Bernhard Schölkopf. Nonlinear causal discovery with additive noise models. Advances in neural information processing systems, 21, 2008.
[21] Kun Zhang and Aapo Hyvärinen. On the identifiability of the post-nonlinear causal model. In Proceedings of the Twenty-Fifth Conference Uncertainty in Artificial Intelligence, UAI '09, page 647-655, Arlington, Virginia, USA, 2009. AUAI Press.
[22] Alexander Immer, Christoph Schultheiss, Julia E Vogt, Bernhard Schölkopf, Peter Bühlmann, and Alexander Marx. On the identifiability and estimation of causal location-scale noise models. In International Conference on Machine Learning, pages 14316-14332. PMLR, 2023.
[23] Nick Pawlowski, Daniel Coelho de Castro, and Ben Glocker. Deep structural causal models for tractable counterfactual inference. Advances in Neural Information Processing Systems, 33:857-869 -69, 2020.
[24] Ilyes Khemakhem, Ricardo Monti, Robert Leech, and Aapo Hyvarinen. Causal autoregressive flows. In International conference on artificial intelligence and stafistics, pages 3520-3528. PMLR, 2021.
[25] Pedro Sanchez and Sotirios A. Tsaftaris. Diffusion causal models for counterfactual estimation. In First Conference on Causal Learning and Reasoning, 2022.
[26] Richard M Karp. Reducibility among combinatorial problems. In Complexity of computer compulations, pages 85-103. Springer, 1972.
[27] Dimitri Bertsekas. Convergence of discretization procedures in dynamic programming. IEEE Transactions on Automatic Control, 20(3):415-419 1975.
[28] Karl Hinderer. Lipschitz continuity of value functions in markovian decision processes. Mathe- matical Methods of Operations Research, 62:3-22, 2005.
[29] Emmanuel Rachelson and Michail Lagoudakis. On the locality of action domination in sequential decision making. 2010.
[30] Jason Pazis and Ronald Parr. Pac optimal exploration in continuous space markov decision processes, In Proceedings of the AAAI Conference on Artificial Intelligence, volume 27, pages 774-781, 2013.
[31] Matteo Pirotta, Marcello Restelli, and Luca Bascetta. Policy gradient in lipschitz markov decision processes. Machine Learning, 100:255-28.: 2015.
[32] Felix Berkenkamp, Matteo Turchetta, Angela Schoellig, and Andreas Krause. Safe model-based reinforcement learning with stability guarantees. Advances in neural information processing systems, 30, 2017.
[33] Kavosh Asadi, Dipendra Misra, and Michael Littman. Lipschitz continuity in model-based reinforcement learning. In International Conference on Machine Learning, pages 264-273. PMLR, 2018.
[34] Ahmed Touati, Adrien Ali Taiga, and Mare G Bellemare. Zooming for efficient model-free reinforcement learning in metric spaces. arXiv preprint arXiv:2003.04006 2020.
[35] Omer Gottesman, Kavosh Asadi, Cameron Allen, Sam Lobel, George Konidaris, and Michael Littman. Coarse-grained smoothness for rl in metric spaces. arXiv preprint arXiv:21 10. /2276, 2021.
[36] Cem Anil, James Lucas, and Roger Grosse. Sorting out lipschitz function approximation. In International Conference Machine Learning, pages 291-301. PMLR, 2019.
12

---

[37] Henry Gouk, Eibe Frank, Bernhard Pfahringer, and Michael J Cree. Regularisation of neural networks by enforcing lipschitz continuity. Machine Learning. 110:393-416, 2021.
[38] Arash Yasr-Esfahany and Emre Kiciman. Counterfactual (non- jidenifiiiiif of learned struc- tural causal models. arXiv preprint arXiv:2301.0900 1,9031, 2023.
[39] Scott Sussex, Caroline Uhler, and Andreas Krause. Near-optimal multi multi-perturbation experimen- tal design for causal structure learning. Advances in Neurul Information Processing Systems, 34:777-788, 2021.
[40] loana Bica, Daniel Jarrett, and Mihaela van der Schaar. Invariant causal imitation learning for generalizable policies. Advances in Neural Information Processing Systems, 1-3952239944, 2021.
[41] Chris Cundy, Aditya Grover, and Stefano Ermon. Bed nets: Scalable variational approaches for bayesian causal discovery. Advances in Neural Information Processing Systems, 34:7095-7110, 2021.
[42] Virginia Aglictti, Neil Dhir, Javier González, and Theodoros Damoulas. Dynamic causal bayesian optimization. Advances in Neural Information Processing Systems, 34:10549-10560, 2021.
[43] Xinshi Chen, Haoran Sun, Caleb Ellington, Eric Xing. and Le Song. Multi-task Ieariing of order- consistent causal graphs. Advances Neural Information Processing Systems, 34:11083-11095, 2021.
[44] Guy Tennenholtz, Uri Shalit, and Shie Mannor. Off-policy evaluation in partially observable environments. In Proceedings of the AAAI Conference Artificial Intelligence, volume 34, pages 10276-10283, 2020.
[45] Hongseok Namkoong, Ramtin Keramati, Steve Yadlowsky, and Emma Brunskill. Off-policy policy evaluation for sequential decisions under unobserved confounding. Advances Neural Information Processing Systems, 33:188/ 18889188811, 2020.
[46] Lingxiao Wang, Zhuoran Yang, and Zhaoran Wang. Provably efficient causal reinforcement learning with confounded observational data. Advances in Neural Information Processing Systems, 34:21 164-21175, 2021.
[47] Junzhe Zhang and Elias Bareinboim. Near-optimal reinforcement learning in dynamic treatment regimes. Advances in Neural Information Processing Systems, 32, 2019.
[48] Junzhe Zhang and Elias Bareinboim. Designing optimal dynamic treatment regimes: A causal reinforcement learning approach. In International Conference on Machine Learning, pages 11012-11022. PMLR, 2020.
[49] Ilya Shpitser and Eli Sherman. Identification of personalized effects associated with causal pathways. Uncertainty in artificial intelligence: proceedings of the... conference. Conference on Uncertainty in Artificial Intelligence, volume 2018. NIH Public Access, 2018.
[50] Juan Correa, Sanghack Lee, and Elias Bareinboim. Nested counterfactual identification from arbitrary surrogate experiments. Advances in Neural Information Processing Systems, 34:6856- 6867, 2021.
[51] Junzhe Zhang, Jin Tian, and Elias Bareinboim. Partial counterfactual identification from observational and experimental data. International Conference on Muchine Learning, pages PMLR, 2022.
[52] Stuart Russel, Peter Norvig, et al. Artificial intelligence: a modern approach, volume 256. Pearson Education Limited London, 2013.
[53] Judea Pearl. Heuristics: intelligent search strategies for computer problem solving. Addison- Wesley Longman Publishing Co., Inc., 1984.
13

---

[54] Alistair EW Johnson. Tom Pollard, Lu Shen, Li-wei H Lehman, Mengling Feng, Mohammad Ghassemi, Benjamin Moody. Peter Szolovits, Leo Anthony Celi, and Roger G Mark. Mimic-iii, a freely accessible critical care database. Scientific data, 3(1):1-9, 2016.
[55] Aniruddh Raghu, Matthieu Komorowski, Imran Ahmed, Leo Celi, Peter Szolovits, and Marzyeh Ghassemi. Deep reinforcement learning for sepsis treatment. arXiv preprint arXiv:1711.0900 2017.
[56] Siqi Liu, Kay Choong See, Kee Yuan Ngiam, Leo Anthony Celi, Xingzhi Sun, and Mengling Feng. Reinforcement learning for clinical decision support in critical care: comprehensive review, Journal of medical Internet research, 22(7)=e18477, 2020.
[57] Taylor W Killian, Haoran Zhang, Jayakumar Subramanian, Mchdi Fatemi, and Marzych Ghassemi. An empirical study of representation learning for reinforcement learning in healthcare. arXiv preprint arXiv:2011.1122 111323,,, 2020.
[58] Mervyn Singer, Clifford S Deutschman, Christopher Warren Seymour, Manu Shankar-Hae Harki, Djillali Annane, Michael Bauer, Rinaldo Bellomo, Gordon R Bernard, Jean-Daniel Chiche, Craig M Coopersmith, et al. The third international consensus definitions for sepsis and septic shock (sepsis-3). Jama, 31 315(())801-810, 2016.
[59] Simon Lambden, Pierre Francois Laterre, Mitchell M Levy, and Bruno Francois. The sofa score developmen,,, utility and challenges of accurate assessment in clinical trials. Critical Care, 23(1):1-9, 2019..
2019.
[60] Flavio Lopes Ferreira, Daliana Peres Bota, Annette Bross, Christian Mélot, and Jean-Louis Vincent. Serial evaluation of the sofa score to predict outcome in critically ill patients. Jama, 286(14):1554-1558 2001.
[61] Jason Waechter, Anand Kumar, Stephen Lapinsky, John Marshall, Peter Dodek, Yaseen Arabi, Joseph E Parrillo, R Phillip Dellinger, Allan Garland, Cooperative Antimicrobial Therapy of Septic Shock Database Research Group, et Interaction between fluids and vasoactive agents on mortality in septic shock: multicenter, observational study, Critical care medicine, 42(10):2158-21 2014.
[62] Lihong Li and Michael L Littman. Lazy approximation for solving continuous finite-horizon mdps. In AAAI, volume 5, pages 1 1175-1180, 2005.
[63] Richard Dearden, Nicholas Meuleau, Richard Washington, and Zhengzhu Feng. Dynamic programming for structured continuous markov decision problems. In Twentieth Conference on Uncertainty Artificial Intelligene (UAI-04). 2004.
[64] Prashan Madumal, Tim Miller, Liz Sonenberg, and Frank Vetere. Explainable reinforcement learning through a causal lens. In Proceedings of the AAAI conference on artificial intelligence, volume 34, pages 2493-2500, 2020.
[65] Ioana Bica, Daniel Jarrett, Alihan Hüyük, and Mihaela van der Schaar. Learning" " what-if" explanations for sequential decision- arXiv preprint arXiv:2007.13531, 2020.
[66] Amir- moss-Hossein Karimi, Julius Von Kügelgen, Bernhard Schölkopf, and Isabel Valera. Algorith- mic recourse under imperfect causal knowledge: probabilistic approach. Advances in neural information processing systems, 33:265-277, 2020.
[67] Stratis Tsirtsis and Manuel Hameez-Rodriguez
Decisions, counterfactual explanations and strategic behavior. Advances in Neural Information Processing Systems, 33:16749-16760, 2020.
14

---

A Further related work
Counterfactual reasoning and reinforcement learning. As mentioned in Section there closely related line of work [1013] that focuses on the development of machine learning methods that employ elements of counterfactual reasoning to improve to retrospectively analyze decisions in sequential settings. Buesing et al. [10] use SCMs to express the transition dynamics in Partially Observable MDPs (POMDPs), and they propose a method to efficiently compute policy based on counterfactual realizations of logged episodes. Lu et al. [2] adopt a similar modeling framework, and they propose a counterfactual data augmentation approach to speed up standard Q-learning. Oberst and Sontag [TT] introduce the Gumbel-Max SCM to express the dynamics of an arbitrary discrete POMPD, and they develop a method for counterfactual off-policy evaluation to identify episodes where a given alternative policy would have achieved higher reward. However, none of these works aims to find an action sequence, close to the observed sequence particular episode, that is counterfactually optimal.
Planning in continuous-state MDPs. Our work has additional connections to pieces of work that aim to approximate the optimal value function in an MDP with continuous states and a finite Therein, the work most closely related to ours the one by Bertsekas 227] shows that, under Lipschitz continuity assumption on the transition dynamics, a value function computed via value iteration in an MDP with discretized states, converges to the optimal value function of the original (continous-
horizon [27 62 63
MDP the discretization becomes finer. Although some of the proof mechanics of this work are similar to ours, the contributions are orthogonal, as we do not employ any form of discretization, and we leverage the Lipschitz continuity assumption to compute optimal action sequences in continuous states using the A algorithm.
Counterfactual reasoning and explainability. Our work has ties to pieces of work that use forms of counterfactual reasoning as a tool towards learning explainable machine learning models. For example, Madumal et al. [64] propose to express the action selection process of reinforcement learning agent as a causal graph, and they use generate explanations for the agent's chosen actions. Bica et al. [65] introduce an inverse reinforcement learning approach to learn an interpretable reward function from expert demonstrated behavior that is expressed terms of preferences over potential (counterfactual) outcomes. Moreover, our work is broadly related to the work on counterfactual explanations (in classification) that aims to find minimally different set of features that would have led to different outcome, in settings where single decisions are taken 66 67]
15

---

B Supporting table and figure
Table It summarizes the most important notation used in the main body of the paper. Symbol/Definition
Description
Time horizon
S, S, A E A
State and action (random variables) at time
81, ai
State and action (observed values) at time
::S. A I
Reward function
T
{(s,. 0z) 1=)
Observed episode
EE R(s(,0))
Outcome of episode
U. EU
Transition noise (random variable) at time
111
Transition noise (inferred value) at time /
9s 8 x Axu +8
Transition mechanism (def. in Eq. 2
Ka,u
Lipschitz constant of gs under A, = a and U, = u
C.
Lipschitz constant of R under A,
S = SX T-1
Enhanced state space
18,1) E S+
Enhanced state representing c/f state 8 after I action changes
F7 : St X A -ASS* S+
Time-dependent transition function (def. in Eq.
(8 d.), 1
Counterfactual episode
Ei a()
Counterfactual outcome of episode r' Max. c/f reward achievable for from enhanced state (s. time
(s,
Sa 1
Fr,,((s,))a)
Enhanced state resulting from (s. () with action
v (s,1,t)
Node in the search graph of the A algorithm
Up 0,)
Root node
UT
8g, T))
Goal node
Reward accumulated along the path that A followed from vo to Lipschitz constant V,, (8,l,t)
Lt
(1,,,
Heuristic function approximating (,, Anchor set
S
S
S
A,
U
Z
U
Figure 4: Directed acyclic graph G' of the SCM C, modeling sequential decision making process (figure adapted from (I3). Green boxes represent endogenous random variables and red boxes represent exogenous noise variables. All exogenous random variables are root nodes of the graph, and each one independently sampled from respective prior distribution. Each endogenous variable an effect of its ancestors in the graph G, and takes its value based on Eqs. do(A; = a) breaks the dependence of the variable A, from its ancestors (highlighted by dotted lines) and sets its value to After observing prediction corresponds to an intervention do(A, = a) in a modified SCM where U, takes values u from posterior distribution with support such that $1+1 ( a, u.).
[ and An intervention
event S 8t+1, S,
A, ae, a counterfactual
16

---

C
Counterfactual identifiability of element-wise bijective SCMs
In this section, we show that element-wise bijective SCMs, a subclass of bijective SCMs which we formally define next, are counterfactually identifiable.
Definition 8. An SCM C is element-wise bijective iff it is bijective and there exist functions gs, R x A x R- - I. with i E {l...., 8141 gs 8t, aq, u.), it holds that $1+1,5 95, a2 Ut,s) for i DI)....D..
......,,, D} such that, for every combination of 8(+1 8t, Qe Ue with
Under our assumption that the transition mechanism gs is continuous with respect to its third argument, it is easy to see that, for any element-wise bijective SCM, the functions gs. are always strictly monotonic functions of the respective ut Based on this observation, we have the following theorem of counterfactual identifiability:
Theorem 9. Let C and M be two element-wit wise bijective SCMs with transition mechanisms gs and h.s respectively, and for any observed transition (St, Gt, 8 ), let ut = 9s (8t, 0t, 8141) and ün hi (8r, ae; 1.+1). Moreover, given any 8 E SS a A. let s' gs(8, Uz) and s" = hi((,, u.). If pe (S S: = 8, A: = a) = (S S: 8, A, a) for all E S,Sa E A, it must hold that
Proof. We prove the theorem by induction, starting by establishing the base case s = sy. Without loss of generality, assume that both gs.1 and hg,1 are strictly increasing with respect to their third argument. Since the two SCMs entail the same transition distributions, we have that
pe (St+1. S S841111 St = 8t, At A =ba=aa a.) (a)= = PM (St+1. S 8841,11 St 82, Ar = at) (
pe (gs,1 (st; ac, U211 S 98,1 0e, (4,1)) pM (hs,1 (st, ar, U2,1) S h (8t, at, (n,1)) (xx)
pe (U1,1 S 1111) PM (U1,1 S i1,1),
where (+) holds because both SCMs are element-wise bijective, and (**) holds because IS,1 and h5,1 are increasing with respect to their third argument. Similarly, we have that
pC (S141,1 S 8 S, = 8, A, = a) pC (gs,1 (8, a, U,,) S gs,1
a, (s,a,
pe (U1,1 S $1,1)
pm (U1,1 S ut,1)
pm (hs,1 (8, U1,1 - hs (s,a, (.,1)) PM (St+1,1 S 5 S = A, a)
pc (S++1.1 S 8 S 8, a)
where in (+), (+) we have used the monotonicity of gs and hg, and (**) follows from the previous result. The last equality implies that s and s correspond to the same quantile of the distribution for Sp-1.1 S.
A, a. Therefore, it is easy to see that s s since the opposite would be in contradiction to Is.1 being bijective. Note that, we can reach that conclusion irrespective of the direction of monotonicity of g5,1 and hs,1, since any change in the direction of the inequalities happening at step (**) is reverted at steps (*) and (t).
Now, starting from the inductive hypothesis that si S/ for all [I....,,. with n D, we show the inductive step, i.e., 5y14 = S"+1- Again, without loss of generality, assume that both 98.n+1 and hg.441 are strictly increasing with respect their last argument. Note that, the two SCMs entail the same transition distributions, i.e. the same joint distributions for S St, Ar. Following from the law of total probability, they also entail the same conditional distributions for Si41n+1| S1+1.5n: St. Ae, where we use the notation I<n to refer to vector that contains the first n elements of D-dimensional vector æ. Therefore, we have that
PC (S1+1,441 S 8c+1,141 S411.5 = 81+1,5n; S. = 8t, A: = az) =
pd (S S S141.411 S5 = 81-1,5n, S. = St, A. = =a,) (
17

---

pe (gs,t41 (st, az Ut,+1) S g8.++1 (st, Qr. u.,+1) Uisa = U1,5n)
= PP (hs, (st, on Utn+1) S hgg (st, Ge, dr,n+1) Ut,5 un (*))
pe (ULn+1 & ut,++1 U1,5n = u. = Pd (Utn+1 S ugn+1 Ui,5n ún,5n),
where for the first equality have used the inductive hypothesis, (*) holds because both SCMs are element-wise bijective, and (+*) holds because 98,++1 and hs, are increasing with respect to their third argument. Similarly, get that
pe (5441 S s S4115m S141,5n: =
pe (gs,-1 (8, a, Ulm+1) S g8.++1 (8, a, ve,+ 195.5n (8, a, Ut,5n) = ga.sn (8, a, u.
pe (Unn+1 uen+1 Utinn uc,gn)
PM (Ut,t+1 S Waat1 U1,5n ün,5a)
1 Pm (hs,tt1 (8, Ut,1) & h8,+t1 (8,0, (.,1 hss (8, Ut,gn) ha,ga (8, (s,
P. (St+1,41 8 S1t1,gn = S111,5 St 8, A a)
p (St41,41 S 4/ St+1,5n 8411.n, S. 8, A, = a):
where (+), (t) we have used the invertibility and monotonicity of gs and hs, and (**) follows from the previous result. With the same argument in the base case, the last equality implies that Su+1 = " That concludes the proof.
18

---

D Proofs
D.1 Proof of Theorem
We prove the hardness of our problem as defined in Eq. by performing a reduction from the partition problem 26), which is known to be NP-Complete. In the partition problem, given multiset of B positive integers V = [ui...., FB} and the goal is to decide whether there is a partition of V into two subsets V1, 1/2 with V12 n V/2 = 0 and V1 U V/2 V, such that their sums are equal, i.e.,
E Evv, 1/-
E0,ev, Vi-
Consider instance of our problem where S = U R² A contains 2 actions ddift' Anall and the horizon is = B + 1. Let be an element-wise bijective SCM with arbitrary prior distributions pe(), (U,) such that their support is on I and transition mechanism Is such that
St, S + U. and
S 0
gs(St, Odiff, Ui)
( Onull: U2
+
Ur.
(12)
Moreover, assume that the reward function is given by
(V)
sum (V)
R(Si, 4air) R(Si, (manl) = max 0, Sr,1
2 S1,2
2
(13)
sum(V)
sum (V)
max 0, 2
-S,,1
S1,
2
where sVm is the sum of all elements continuous suggested by the following Lemma (refer to Appendix I.I.J for a proof).
Σ
1=1 Us Note that, the SCM C defined above is Lipschitz-
Lemma 10. The SCM defined by Equations 12 13 Lipschitz-
according to Definition ?
Now, assume that the counterfactual action sequence can differ in an arbitrary number of actions from the action sequence in the observed episode i.e., k = T and, let the observed action sequence
be such that a, = amall for t E {0.... T - 1}. Lastly, let the initial observed state be so = [0, vil.
the observed states {8) be such that 8; = Σi. U/ ut+1 for t E {I...., T - 2} and the
last observed state be -T-1 [sum(V), O]. Then, is easy to see that the noise variables U, have posterior distributions with point mass on the respective values
e++1 Pr+2
for t E (O,...,T--3 3} and Ur = PT
Up
Note that, for all E (1..... T 2}. we have S $t,1 < sVm(V) and 81, 2 1, hence the immediate reward according to Eq. [3 is equal to O. Consequently, the outcome of the observed episode is o (7) R( ) max(0, *um(()) max (0, sam(V)
svm
2
Next, we will characterize the counterfactual outcome o(T') of counterfactual episode r' with sequence of states {s()7-) resulting from alternative sequence of actions {a TT)) Let NY denote the set of time steps until time f, where the actions taken in a counterfactual episode r' are aair and auull respectively. Formally, D/- (t' E {o. ar daa)) N² {t' E {o...... ar = amall). Then, as an intermediate result, we get the following Lemma (refer to Appendix D.1.2 for proof).
Lemma 11. holds that s',1 = =EVEN; for all E (l....T-- 13.
Following from that, we get that 0 < s/,1 S sum(V) for all E f1...., 1}. Moreover, we can observe that the transition mechanism given in Eq. 12 is such that ( Ar, U222) = U, for all E (O....,
, T 2), independently of Sr and A, Therefore, holds that $4,2 4 for (1. ......... T 2). and $6,2 S0.2 = 01 1. a direct consequence, is easy to see that R(s a() for all e {O...., T 2}. and the counterfactual outcome is given by o + (T') = R(
(14)
In addition that, we have that ur 2,2 0, hence
T Σier U+t1
(15)
0
19

---

Now, we will show that, if we can find the action sequence {až |2=0 that gives the optimal coun- terfactual outcome o* (T*) for the aforementioned instance polynomial time, then we can make a decision about the corresponding instance of the partition problem, also in polynomial time. To this end, let {s:) rt=
T-1 be the sequence of states in the optimal counterfactual realization and, let
Di-2- {t {o...... T 2} az = Gaif): N {t' (O...., 2] ag. = Qgull))
From Eq. [4| we get that the optimal counterfactual outcome is o (7*) R(s 1:: af-1), and it is easy to see that the reward function given in Eq. 13 is always less equal than zero. If ((())) 0, has to hold that
max 0, 07.57-11
sum(V)
sum (V)
)
ST
2
2
sum(V)
sumn(V) (V)
max
0,
2
-8-1111
ST
1.2
2
Σ MENT
v++1
sum V)
0
and
())
Σ
1i+1
< =
2
2
AEN
1/41 swm(V)
Σ
1EN
where (+) follows from Eq. 15 As a consequence, the subsets V1 = {us i 1 E N and 1/2 = (1; 2) partition V and their sums are equal.
On the other hand, if o (r*) 0, as we will show, there is no partition of into two sets with equal sums. For the sake of contradiction. assume there are two sets V1, 1/2 that partition V. with (V)) = sum((2)
sum(V)/2, and let N7-2 {t E {0. (o. 2} Ut+1 E V1 and DI-2 = {t E {0, (O....,7 - 2} Ut+1 V2). Then, consider the counterfactual episode 7' with an action sequence (a])==) such that its elements take values Onull and Oar based on the sets 24 22 respectively, with -1 taking an arbitrary value. It is easy that
o+('') R(s4-1,
R { EiENg
U++1
4T-1
0
swm(V)
sum(V)
Σ
max
0,
1t+1
max
0,
2e+1
2
2
15
CEN;
sum(V)
sum (V)
surn(V)
sum(()
max
0,
max
0.
2
2
2
= = 00>00(m) (r*),
which is a contradiction. This step concludes the reduction and, therefore, the problem given Eq. cannot be solved in polynomial time, unless P NP.
D.1.1 Proof of Lemma 10
easy to see that, for all u E 24 and for all s, E S, the function gS Onall, 2) satisfies (s Gmall, u) - gs(s' anull, u)ll < ||s ---s s'Il, and therefore Ka satisfies Definition For the case of A, = daim, we have that
(s Gan, u) gs(s', Oamr uill
81 82
1 5/
(51 - s/1 + (5/2 82)
0
0
[(s 8) (5/ -22) S [1 $1 + 182 s/gb
and therefore (s adift, u) gs((', adiff, u)ll² < (s) - 212 (52 2/)2 2|s1 1/1152 s/g. We also have that
2 |s s'll = v2,/ (s1 - s1)² (s2 2²22 + 2 s s' 2( 81)² + 2(32 -2-2
20

---

By combining these, we get
2 || - s111 lgs (8, Chiir, uA) - (s(s Gairt, uill 2 (8) = (82 x/2 2 8/11 2/21 +
2 0| 8 8'11" - 1las(8, adift, u) gs( adift, 2) 2 (Is1 - 51 - + - s'1)²² 2 0.
Hence, we can easily see that K Kagruu v satisfies Definition
Next, we need to show that, for all E A there exists C. E R4 such that, for all s, E S. it holds [R(s, a) R('s,, a)l S Ca ||s - 8'11. Note that, to show that a function of the form max(0, f())) with f I2a + R is Lipschitz continuous, it suffices to show that f(s) is Lipschitz continuous, since the function max(0, x) with x E R has a Lipschitz constant equal to
We start by showing that the function f()) (V) is positive constant. For an arbitrary pair s, s' E S. we have that
- a is Lipschitz continuous, where
Is()) - f(s') = 181 a1=220 - 2/21 & |s1 s11 + 0:/8 - 2/1 +
If(s) - f(s'))2 (s1 - -1) + (s2 - =2) + 2a|s1 sil|s2 - $221
We also have that
VI VI+a a ||s s' = VI + a / (s1 51 (32 - 5/2²²
(1 + 8 s 2 (1 a) - -11² + (1 a) 222)²
By combining these, we get
(1+a) ||s s'12 1 f(s) f(s 2 a(s1 s1 + a(s2 2222²
22 $1 sills2 2'2 +
a1+ ||s - - (88) /(x') 2 a ([s1 - 811 |82 $ 1) 2 0.
Hence, we arrive to s(x)) < vI+ 1 8/1. and the function f is Lipschitz continuous. is easy to see that the function ((s) = a - 81 822 a is also Lipschitz continuous with the proof being almost identical. As a direct consequence, the reward function given in Equation 13 satisfies
Definition 1 with Comn = Cass 2 1 swm(V)
This concludes the proof of the lemma.
2
D.1.2 Proof of Lemma 11
We will prove the lemma by induction. For the base case of t 1, we distinguish between the cases
a = dairr and a Gmall- In the first case, we have $ = u0, + So, $0,2 = + 0 - 0 and N₁ 0 and, therefore, the statement holds. In the second case, we have $1.1 uoo1 + 80,1
"1 + 0 = 01 .²¹ = {0) and EvENY U1 U1. Therefore, the statement also holds.
For the inductive step (t 1), assume that 5/ا1, EvEN; ur+1 and will show that
K, = ΣveNy, w0+1- Again, we distinguish between the cases d/ = car and d/ = auill
However, note that, in both cases, 6/-1.2 = U1- + = wr. Therefore, in the case a = ddiff. we get
1 - U1-1, +1-11115 166111+6111111111
Σ
U++1 U2
Σ MEN:
Up41
Σ
vP+1s
V'NN
MENL
where the last equation holds because d/- Qai and, therefore, N No 2. In the case of d' = Gmall we get
1/11 U1-1,1 + 1/11,11
Σ
U++1 Σ 1841
MEN
PENT
where the last equation holds because 0/-1 Onull and, therefore, No = Ni_ U ft (t-1)) 1}.
D.2 Proof of Lemma 5
We will prove the proposition by induction, starting from the base case, where
First,
Let = T 1 and = k. It is easy to see that, if the process is at state E S in the last
21

---

time step with no action changes left, the best reward that can be achieved is R(s. aT) as already discussed after Eq. 10 Therefore, holds that |V,(s, - 1) V- k. - [[(s, aT- R(s', -ar-i)
Car- s'1 I C s'Il, where the last step holds because C = A Ca. Now, consider the case of = 7 - with taking an arbitrary value in {O...., k- 1}. Let 8, s' be two states in S and a" be the action that gives the maximum immediate reward at state 8, that is, a" = arg [[(((())))))))) Then, we get
(V,, (8 1)1 max{ max{ a)}l
aEA aEA
aEA
| RRss' a Ca. ||s s'1l s s'1l
where (x) follows from the fact that R(s', a²) I maX R(s',a)}. Therefore, for any 1 E (0. ........h kk] and 8, 8 E S, it holds that (V, -1 V ,,, 1)] LT- | 8'I, where
LT-1= CC
Now, we will proceed to the induction step. Let L < T - 1, 1 < k and, as an inductive hypothesis, assume that Lt+1 E R as defined in Lemma is such that, for all / {O.....k) .k and 8, s' E S, L ||s - s'll- Additionally, let (s):la), (s'2 la) denote the enhanced states that follow from (s,1), (5', after taking an action a, i.e. (8a: (a) = Fte ((s, 1) a) and (s'a (a) = Pt, ((6', (),))) a). Lastly, let a" be the action that maximizes the future total reward starting from state 8, Le., a" = argmax [R(s, have that
it holds that (v, + 1) (,', 1, 1)
(R(s, a) + V. (8a, la I + 1)}. Then, we
(o,10,))) - Vi(6', (1,)) |mmax R(ss, + Vz (sa: la; 1 + 1)] -max max( a) + V, la, 1)31 aEd
aEA
|1(5, a") + +V (sa-
R(s's a²) V, (5' la-
|R(s a²) -R(s', a²) + IV~ haa+, 1) V, (5 la-, 1)
(xx)
% C. ||s - 8'1l + Le+1 |1s0-
8/-11
Saa ||s - -s'l + L1+1 Ka me ||8 8'll
(1x*)
I C||s- C-|s-6 s'll + L Ku, ||8-
= (C + Ta+1 Kw,) ||s - s'll = Toe || 8'1l-
In the ahove, (+) holds due to R(s',a") V- (a'a-
max R(s', a)
Vr (8 lay 1 + 1)}} (**) follows from the inductive hypothesis, and (* (***) holds because C m C. and Ku, maxas Ka,u, is easy see that, similar arguments hold for the simple case of = k, therefore, we omit the details. This concludes the inductive step and the proof of Lemma 5
D.3 Proof of Proposition 6
We will prove the proposition by induction, starting from the base case, where (. = 1- the algorithm initializes V. (s, 1,T 1) to max R(s, a) for all St, ( E (0, k 1} and V. (s. ks T 1) to R(s, -1) It easy to see that those values are optimal, as already discussed 1) follows trivially.
after Eq.
10 Therefore, the base case (s, (, - 1) 2 V, (s, Now, we will proceed to the induction step. Let f < T T'11 and, as an inductive hypothesis, assume that V. (s 2 V= 1 1) for all s E Sy. E {O..... k). Our goal is to show that V. (s. I, t) > V= L, for all s E S, E {O...., k). First, let I < k. For a given point E Sy. Algorithm I finds the next state sa that would have occurred by taking each action a, i.e., (sa, la) = Fr,, ((s, (),a))))) a).
and it computes the associated value V. = mingfes, (V, (8t, t + 1) L || sall). Then,
22

---

simply sets (s, equal to max (R(ss VaJ. We have that
Va min
(V. (8, ++, L141 ||8 8all} #,E8
min (V (V l, s1Ett
Le+1 |8 - 8all}
(xx)
2 min {V, (Sa + 1)}
st5t
V, (sa t.,
where (+) follows from the inductive hypothesis, and (**) is consequence of Lemma Then, we get
V=(8,1,t) = max (R(s, + V.J 2 max (ixx + (8 + 1)} = (8,,,) aEA
aEA
Additionally, when = k we have k,) a(s,a:) + ming, es, (V, k,t t + 1) + L++1 [|8 8a,ll |} and V k,, a(s,a;) + (8,,, k, + 1). Therefore, the proof for (s k,t 2 (s k, t) is almost identical.
D.4 Proof of Theorem
We start from the case where t = 1. v= (8, 1) and, consider an edge associated with action a" connecting the goal node vr (sp, k. T) that carries reward R(s, a²). Then, have
max R(s,a) 2 R(s, a R(s,a") k.
aEA'
and the base case holds.
For the more general case, where < T 1, first establish the following intermediate result, whose proof is given in Appendix D.4.1
Lemma 12. For every 8, 8' S, (O...., E {0, T 1}, holds that (V((,,,,) V-(s', 1,1)1 S L, ||s s'Il, where L, is as defined in Lemma
That said, consider edge associated with action a' connecting node (8,4t) node Ua- = (8ar, la+, 1). Then, we have
V=(,, = max "EA''
R(s,a) + min e,E,,
(5) la, 1) + L1t1 ||s| - sull
2 R(s,a")+ + min
Vy (8 la, + 1) + Tst1 || - 80-11 180-81
9,58,
2 R(s, a'ssaa)+ + min e,58,
lar, tor +1)
R(s, a²) (=(su+ +++
That concludes the proof and, therefore, the heuristic function V. is consistent.
D.4.1 Proof of Lemma 12
Without loss of generality, we will assume that / < k, since the proof for the case of | = k is similar and more straightforward. We start from the case where f T 1 and, for states s, s' have
- 1)]
R(s, a) max R(s'.a)
aEA
aEA
V=
V, 1) < C ||s -s LT |s -s(]
where the last inequality follows from Lemma
23

---

Now, consider the case t < T 1, and let (sa, (a) denote the enhanced state that follows from (s, 1)
after taking an action at time t, i.e., (sa, (a) = F7., ((s, () a). Then, we have
(s, V. (,',',,,t))
la, + + L ||s| sall}}
= max aed
R(s + min stESt
(Vr la, + Let1 || ||+- - s'ull 1 <6ll}} (16)
R(s',a) + min max aEA
stE3,
Let a" be the action A that maximizes the first part of the above subtraction, i.e..
[((s,)) min syEs
(Vr ta, + 1) Le+1 ||st sa
arginxx aEA
Then, Eq. 16 implies that
-V (s), R(s,a') + min 4,E8,
(oo, 1) + Let1 [|s| sa- 8a-1]
-R(s' a²) min (1 1) + T441 ||st 86-11
t158
[(ss.a)) - R(s' a")l
{ (s la 1) Ltt1 ||s sa. s..11 [|}
min #4ES,
(Vr(s + 1) L (4illss - s¹a. s/-
min 4,E8,
(17)
Now, let & be the 8 E St that minimizes the second part of the above subtraction, i.e.
8 = argmin { (8t 1) + T4+1 || - 8'a-11
a,ES,
As consequence and in combination with Eq. we get
(s,1,,))) V=(6,,,,0) < IR(8 a²) R(ss'a")
+ (i la- + L1+1 || 8a-
V.(8, 1a- 1) + T441 | - -
1 || - 11-0-'0-11 = R[(s R(s'sa' aa") ||+ L141 ||| Su ||
(0)
Ca ||s s'1 + L1+1 ||sa-
8/--11
(+x)
C||s s'll
L441ku, ||s s'll L,
s'
where in (+) use the triangle inequality and the fact that the SCM C is Lipschitz-comtinuuus, and (**) follows from Lemma
24

---

Algorithm 2: Graph search via A*
Input: States S, actions A, observed action sequence {at Jr=0 1==7-1 horizon T, transition function F7, reward function R, constraint k, initial state so, heuristic function Vy. Initialize: NODE
{"tuple": (so, O, 0), "rwd" 0, "par" Null, "act" /Nul]. STACK action_sedd setiionseqqence t-
QUEUE Q -- {root} SET explored +
while True do
argmax.ee {v' .rudd (v tuple)): Q Q (o: /* Next node to visit */ v.tuple = (s, 8 T) then
while v.par * Null do
action action_sequenee sequence.pushi(e.cct
/* Retrieve final action sequence */
== v.par
end
return acciion_sequene end
explored + explored U {v}
/* Set node as explored */
it = then availabe_actin + { else
cvailabe_actions + end
for
a E available acailableeactions
actions do
(s, v.tuple
(8a, - Fit,, ((s,1),a)
/* Identify w's children nodes */
+ {"tuple (8a, la. 1), "rwd" v.rwd + R(s, a), "par" "act" a)
c, $ Q and va & explored then
Q QU(va)
/+ Add them to the queue if unexplored */
end
end
end
E A algorithm
Algorithm 2 summarizes the step-by-step process followed by the A" algorithm. Therein, represent each node by an object with attributes: (i) the "tuple" (s, I, t) of the node, (ii) the total reward "rwd" of the path that has led the search from the root node up the node v, (iii) the parent node "par" from which the search arrived to and (iv) the action "act" associated with the edge connecting the node v with its parent. In addition the queue of nodes to visit, the algorithm maintains set of explored nodes, and it adds new node to the queue only if it has not been previously explored. The algorithm terminates when the goal node is chosen to be visited, i.e. the "tuple" attribute the format (*, *, 7), where * denotes arbitrary values. Once the goal node ur has been visited, the algorithm reconstructs the action sequence that led from the root node to the goal node and returns it output.
has
25

---

3.10
1300
2.0
5.05
1.9
19
300
10000
15000
20000
20000
20000
10000
10000
20000
25000
30000
12000
18000
20000
20000
20000
Size of onchor
anchor
Size
(8))
(a) MC-Lipschitz
(b) MC-Uniform
(c) Facility-ILocation
Figure 5: Computational efficiency our method, measured by means of the Effective Branching Factor (EBF), under three different anchor set selection strategies against the size of the anchor set Sp. In all panels, we set LA 1.0, La 0.1 and
Error bars indicate 95% confidence intervals
over 200 executions of the A algorithm for 200 patients with horizon T 12.
F Experimental evaluation of anchor set selection strategies
In this section, we benchmark the anchor set selection strategy presented in Section against two alternative competitive strategies using the sepsis management dataset and the same experimental More specifically, we consider the following anchor set selection strategies:
setup as in Section 5
(i) MC-Lipschitz: This is the strategy described in depth in Section 4 based Monte Carlo simulations of counterfactual episodes under randomly sampled counterfactual action se- quences. Notably, the time steps where each counterfactual action sequence differs from the observed one are sampled proportionally to the respective Lipschitz constant L, of the SCM's transition mechanism. To ensure fair comparison with other strategies, instead of controlling the number of sampled action sequences M, we fix the desired size of the anchor set St, and we repeatedly sample counterfactual action sequences until the specified size is met.
(ii) MC-Uniform: This strategy is a variant of the previous strategy where we sample the time steps where each counterfactual action sequence differs from the observed uniformly random, rather than biasing the sampling towards time steps with higher Lipschitz constants Lr.
(iii) FacilityLLccatio: Under this strategy, the anchor set is the solution minimax facility location problem defined using the observed available data. Let S, be the union of all state vectors observed in all episodes T in a given dataset. Then, choose an anchor set S C S, of fixed size |S; b, such that the maximum distance of any point in S, to its closest point S is minimized. Here, the rationale that counterfactual states resulting from counterfactual action sequences for one observed episode are likely to be close to the observed states of some other episode in the data. Formally.
St argmin max min {|ls s'11} S'CS.:(8"
(18)
(ESS ES
Although the above problem is known to be NP-Complete, farthest-point clustering algorithm, which is known to have an approximation factor equal and runs in polynomial time. The algorithm starts by adding one point from S, to S at random. Then, proceeds iteratively and, at each iteration, it adds to S the point from S, that is the furthest from all points already in S,, i.e., S, S₁ U 8, where max (mina,es, || s+1]]- The algorithm terminates after b iterations.
find solution using the
Results. We compare the computational efficiency of our method under each of the above anchor set selection strategies for various values of the size of the anchor set [Sy]- Figure 5 Summarrzes the results. We observe that the Facility-Location selection strategy performs rather poorly compared to the other two strategies, achieving an effective branching factor (EBF) higher than 3. contrast, the MC-Lipschitz and MC-Uniform strategies achieve an EBF close to 2, which decreases rapidly as the size of the anchor set increases. Among these two strategies, the MC-Lipschitz strategy, which we use in our experiments in Section 5 achieves the lowest EBF.
26

---

Table 2: Levels of vasopressors and intravenous fluids corresponding the 25 actions in A
Vasopressors (mcg/kg/min)
Intravenous fluids (mL/4 hours)
0.00
0
0.04
30
0. 113
80
0.225
279
0.788
850
G A Additional details on the experimental setup
G.1 Features and actions in the sepsis management dataset
R,,, where D 13 is the number of features. We distinguish between three types of features: (i) demographic features, whose values remain constant across time, (ii) contextual features, for which we maintain their observed (and potentially varying) values throughout all counterfactual episodes and, (ii) time-varying features, whose counterfactual values are given by the SCM The list of features is as follows:
As mentioned in Section our state space
Gender (demographic)
Re-admission (demographic)
Age (demographic)
Mechanical ventilation (contextual)
FiO₂ (time-varying)
PaO₂ (time-varying)
Platelet count (time-varying)
Bilirubin (time-varying)
Glasgow Coma Scale (time-varying)
Mean arterial blood pressure (time-varying)
Creatinine (time-varying)
Urine output (time-varying)
SOFA score (time-varying)
To define of actions A follow related work 6
and consider 25 actions corresponding to 5 x levels of administered vasopressors and intravenous fluids. Specifically, for find all non-zero values appearing in the data, and divide them
55 57
both vasopressors and fluids, into 4 intervals based on the quartiles of the observed values. Then, we set the 5 levels to be the median values of the intervals and 0. Table G.I shows the resulting values of vasopressors and fluids.
G.2 Additional details on the network architecture & training
We represent the location and scale functions and ø of the SCM C using neural networks with 1 hidden layer, 200 hidden units and tanh activation functions. The mapping from state s and an action a the hidden vector takes the form z tanh(W, + W.a). where a is 2-D vector representation of the respective action. The mapping from the hidden vector z to the network's 's output done via fully connected layer with weights W.. enforce a network have a Lipschitz constant I. with respect to the state input, apply spectral normalization to the weight matrices W, and W., so that their spectral norms are [Wall, | = 1. Additionally, we add intermediate layers between the input and the hidden layer and between the hidden layer and the output layer, each one multiplying its respective input by constant VT. Since it is known that the lanh activation function has a Lipschitz constant of 1, it is easy to see that, by function composition, the resulting network is guaranteed to be Lipschitz continuous with respect to its state input with constant L. Note
27

---

0.01
0.02
0.05
-0%
0.1
0.2
0.5
1.0
-0%
1.1
1,2
1.5
0.01
0.02
0.05
2
Ln
SCM C. measured by means of the percentage decrease in log likelihood of the data comparison with SCM trained without Lipschitz-continiig constraints. The and axes correspond to different enforced values for the Lipschitz constants Ln- Lø of the location and scale networks h and respectively. Darker values indicate that the learned SCM achieves significantly lower log likelihood than the unconstrained SCM.
Figure 6: Goodness of fit of the
that, since the matrix Wa is not normalized, the network's Lipschitz constant with respect to the action input can be arbitrary.
To train the SCM C, for each sample, we compute the negative log-likelihood of the observed transition under the SCM's current parameter values (i.e., network weight matrices & covariance matrix of the multivariate Gaussian prior), and we use that as a loss. Subsequently, we optimize those parameters using the Adam optimizer with learning rate of 0.001, batch size of 256, and we train the model for 100 epochs.
We train the model under multiple values of the Lipschitz constants Ln, Ls of the location and scale networks, and we evaluate the log-likelihone likeliieeoodd of the data under each model using 5-fold cross- validation. Specifically, for each configuration of L, and Lap we randomly split the dataset into training and a validation set (with a size ratio 4-to-1), we train the corresponding SCM using the training set, and we evaluate the log-likelihood of the validation set based on the trained SCM. This results in the log-likelihood always being measured on different set of data points than the one used for training. For each configuration of Li, and La. repeat the aforementioned procedure 5 times and we report the average log- log-likelihood achieved on the validation set. addition, we train an unconstrained model without spectral normalization, which can have arbitrary Lipschitz constant.
Figure 6 shows the decrease in logeliiildoodd of the respective constrained model as percentage of the ligliiklidood achieved by the unconstrained model, under various values of the Lipschitz constants Ln, La- We observe that, the model's performance is predominantly affected by the Lipschitz constant of the location network T,,, and its effect is more pronounced when L,, takes values smaller than 1. Additionally, we can see that the scale network's Lipschitz constant L has milder effect on performance, especially when L, is greater or equal than 1. Since we are interested in constraining the overall Lipschitz constant of the SCM C, in our experiments in Section 5 we set La = 1 and La 0.1, which achieves log-likelihood only 6% lower to that of the best model trained without any Lipschitz constraint.
28


Reviews:
[{"review_id": "A7jZx1GB5x", "reviewer": "Reviewer_R23o", "summary": "The paper studies the problem of finding a counterfactual action sequence to maximize the outcome of a trajectory in an MDP characterized by a structural causal model. The focus is on continuous state spaces under a set of Lipschitz constraints. The authors show that this problem is NP hard and propose an algorithm based on A* search, which they test empirical on a dataset on clinical decision making. ", "questions": "- What could be baselines to compare the empirical results to?\n   - What impact does the chosen heuristic have on the performance of the search algorithm?\n   - How much worse would be a method that discretizes the environment?\n\n- How can we be confident that the present dataset is a good benchmark to evaluate the approach?", "limitations": "Yes", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "- The paper studies a relevant well-motivated problem that will be of interest to the NeurIPS community\n- The authors do a good job in clarifying and formalizing the problem\n- Showing that the proposed problem is NP-hard is a valuable insight to guide future theoretical and empirical work\n- The proposed algorithm is a natural choice and it's formulation is sound\n- The empirical evaluation is performed on a realistic dataset, and the results are promising", "weaknesses": "- The discussion of related work is too limited. I think there is a lot of related work about related problems in SCM (outside RL) that should be mentioned. I'm not very familiar with this literature myself, which made it difficult to judge the potential impact of the contribution the present paper makes.\n\n- Theorem 3 (NP-hardness) is a key contribution of the paper. However, it is only stated in the paper and almost not discussed. It would be valuable to give a high level picture of the reduction used to prove the result, as well as a discussion of the implications of this result and how it results in the choice of algorithm later in the paper.\n\n- The empirical results are difficult to interpret due to a lack of baselines. While the efficiency results in Figure 1 and performance results in Figure 2 look as expected qualitatively, there is no point of comparison.\n   - For efficiency it would be useful to compare to some naive baselines, like a search based method without a specific heuristic.\n   - The counterfactual performance improvements in Figure 2 are pretty small which makes me wonder if the dataset is a good enough benchmark to evaluate the present method. Maybe the trajectories in the dataset are too close to optimal?\n   - It could be interesting to compare to an approach the discretizes the environment and uses prior methods for discrete state spaces.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their careful and insightful comments, which will help improve our paper. Please, find a point-by-point response below.\n\n**[Related work]** We agree with the reviewer that there is a rich literature on structural causal models (SCMs) that is not connected to reinforcement learning (RL). Due to space limitations, we chose to focus the discussion only on the most closely related work [10-13] in the main body of our paper and include a further discussion of related work in Appendix A. In the revised version of our paper, we will bring some of the related work from Appendix A to the main body, and we will add a footnote explicitly pointing the reader to Pearl [8] for a more complete overview of prior work on causal inference based on SCMs.\n\n**[NP-Hardness]** The main implication of the NP-Hardness result is that there is no hope in designing a polynomial time algorithm that solves the problem. Following the reviewer’s advice, after Theorem 3, we will include a paragraph where we will give a high level picture of the reduction used in the proof of the Theorem and discuss the importance and implications of the Theorem in the design of our algorithm.\n\n**[Lack of baselines]** We thank the reviewer for their concrete experimental suggestions. To the best of our knowledge, the only natural baseline to compare against would be the method introduced by Tsirtsis et al., NeurIPS 2021 [13] which solves a closely related problem in discrete state spaces, as we discuss in lines 52-54. Unfortunately, their method has a quadratic complexity with respect to the number of discrete states and thus does not scale to continuous multidimensional vector states as those used in our experiments. For example, discretizing the 9 continuous features we consider in our experiments into ten discrete levels each---a rather coarse-grained discretization---would lead to 1 billion discrete states.\n\nIn terms of naive baselines, we believe it is non-trivial to come up with methods other than an exhaustive search that simultaneously work without a heuristic function and guarantee finding the optimal solution. We do not directly compare our method with an exhaustive search because it is too computationally expensive to run in a reasonable amount of time. That said, our experimental results indicate that our method does perform better than an exhaustive search since, otherwise, it would be exploring the entire search space independently of the exact form of our computed heuristic function. For example, varying the number of Monte Carlo samples $M$ would have no effect on efficiency, i.e., Figure 1(b) would be a flat line.\n\n**[Quality of the sepsis management dataset]** We decided to perform experiments using the MIMIC-III sepsis management dataset since it is a commonly used dataset in the literature on reinforcement learning for healthcare (see lines 297-298) and it is relevant to our problem motivation. Note that, as an implication of Theorem 6, our method is guaranteed to *always* find the optimal solution, and we do not have reasons to expect that it would perform poorly using a different dataset. However, we would like to clarify that the reviewer’s comment that *``maybe the trajectories in the dataset are too close to optimal''* is part of our experimental findings (see lines 347-349) and we do not view this as an indicator that the dataset is not a good benchmark for the evaluation of our method."}, {"review_id": "0XcNRjxxtQ", "reviewer": "Reviewer_MQT2", "summary": "This paper studies a question that is very natural for a sequential decision maker to ask itself: \"how could I best improve the return I obtained in a trajectory by only changing a fixed (k) number of actions from the sequence of actions I just executed while keeping the rest of the actions fixed\". Authors study this question in a context where the state space is continuous but the action space is discrete. The paper generalizes the result of previous work which studied a similar question, but in the context of discrete states. The setting is very nicely motivated and the problem is quite relevant and important.\n\nThis formulation seems to be novel, general, and useful. However, the paper shows that answering this question is in fact an NP-hard, one that we can only solve in the worst case using exhaustive search. Therefore, the paper goes on to find an admissible heuristic to be used in conjunction with the A* algorithm to reduce the complexity of the search space.\n\nTo propose a consistent/admissible heuristic, one needs to find reasonable upper bounds for the value of unseen states which are found using counterfactual reasonings. The paper proposes to get this done using the notion of Lipschitz continuity where intuitively, the difference between the value of the two states is upper bounded by a constant times the norm of the difference of the two states. This may be a stringent requirement generally, but it seems to make sense in the domain that is of interest in the paper, namely in medical trials. Experiments show that using this heuristic does result in significantly reducing the average runtime of the search.\n\nOverall, this is a fair contribution in the intersection of causal reasoning and planning/reinforcement learning.", "questions": "I would like to see a comparison with the discrete case, in particular the \"Woulda, coulda, shoulda: Counterfactually-guided policy search\" paper of Buesing et al, 2018. I do understand that this paper is dealing with the continuous state space. But one natural baseline to compare against could for example be the paper mentioned above where one can simply discretize the continuous state space and then apply their algorithm to compare against.\n\nHow important is the kick-starting part with the monte carlo anchoring? In particular, it would have been nice to see ablations where the monte-carlo anchoring is performed more or less to understand its effect on performance.", "limitations": "Authors maintain that they will release code for their experiments if accepted, but they seem to have access to confidential patient data. I hope and trust that authors will they take sufficient precautions to ensure the anonymity of patients whose information is used in these experiments.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The competency that the paper attempts to develop, namely being able to answer counterfactual questions in the context of planning and reinforcement learning, is quite interesting. Formulating this problem and showing that the original problem is NP-hard also sheds some light on the complexity of asking such questions. Some of the assumptions made, such as using a bijective causal model does limit the scope of the result, but I still find the formulation very interesting.\n\n", "weaknesses": "We know that Lipchitz continuity can be a pretty lose upper bound for the function value. Especially in the context of this paper, when data is sparse, the upper bounds derived in Lemma 4 can be pretty bad. Moreover, recently there have been more efforts to develop Lipschitz-like tools that are more conducive to RL algorithms by ensuring that the new tool are Coarse-Grained. See for example \"Coarse-Grained Smoothness for Reinforcement Learning in Metric Spaces\", by Gottesman et al, 2023 or \"Zooming for Efficient Model-Free Reinforcement Learning in Metric Spaces\" by Touati et al 2020. This paper though relies on more classical notions of smoothness (vanilla Lipchitz continuity) which is not the most effective in the context of RL.\n\n\nAnother point is that under the assumptions made in this paper (namely Lipschitz reward and transition dynamics), the value function itself becomes Lipschitz, so when learning the value function one can restrict the set of feasible value functions to be Lipschitz also. Is there any reason why this is not leveraged? In the same vein, we clearly know that the value function is always between (1-\\gamma) R_max and -(1-\\gamma) R_max. Can you elaborate if the bounds that are obtained by the Lipschitz interpolation would be able to provide much tighter bounds or that the bounds become vacuous? Can you also elaborate on the amount of data needed in your medical domains before these bound become non vacuous?\n\nIn practice, we also need to make sure that we do not underestimate the Lipschitz constant, because otherwise the heuristic used in A* will no longer be admissible. The estimate also needs to be tight enough so as to make sure that it is effective when used as a heuristic. But I am not sure how this trade-off will be maintained in problems where the Lipschitz constant is unknown.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their careful and insightful comments, which will help improve our paper. Please, find a point-by-point response below.\n\n**[Lipschitz-like tools for RL algorithms]** We would like to thank the reviewer for bringing these papers to our attention, which we will cite in the revised version of our paper. However, we would like to highlight that, in our paper, we do not learn a model-free RL policy nor study the relationship between smoothness, Lipschitz continuity and the amount of data needed to learn a model-free RL policy. Rather, given an observed sequence of actions and states, we search for an alternative sequence of actions that would have retrospectively maximized the (counterfactual) outcome under the counterfactual transition dynamics given by a (trained) structural causal model (SCM)---this is stated verbally in lines 40-42 and formalized in lines 180-183. In other words, our paper studies an *algorithmic* problem and not a *learning* problem. Therefore, it is hard to make a direct comparison with the papers brought up by the reviewer. That said, we would also like to clarify that the bounds obtained in Lemma 4 do *not* depend on the sparsity of the available data and they are an implication of the Lipschitz continuity of the SCM (see Definition 1).\n\n**[Lipschitz value function]** In our paper, to guarantee that the computed heuristic function $\\hat{V}_ \\tau$ is provably consistent, we do leverage the assumption that the SCM is Lipschitz continuous (see Lemma 4, Proposition 5, Theorem 6 and the related proofs in Appendix C). Also, we would like to clarify that our approach does not involve *learning* a value function. Note that the function $V_ \\tau$ does not denote the value function, as typically defined in the RL literature, but denotes the counterfactual reward that could have been achieved in a counterfactual episode where, at time $t$, the process is at a (counterfactual) state $s$, and there are so far $l$ actions that have been different in comparison with the observed action sequence in the observed episode $\\tau$ (see lines 192-195).\n\nTo our best understanding, the upper bound $(1- \\gamma ) R_{max}$ mentioned by the reviewer is not an upper bound relevant to our function $V_\\tau$ but rather a bound for the value function in an RL problem with an infinite horizon. Note that, in our work, we consider a finite horizon (see lines 31, 60 & 75), and the outcome of the decision making process is the sum of the rewards $o(\\tau)=\\sum_t R(\\mathbf{s}_t, a_t)$, which implies that $\\gamma=1$. Therefore, together with the explanation given in the previous paragraph, we believe that the bound given by the reviewer is not applicable to our problem.\n\nFinally, we would also like to clarify that, to compute a better heuristic function $\\hat{V}_\\tau$, we do not require more data to get better bounds, but rather more anchor points, which we get via Monte Carlo simulations (see lines 277-295).\n\n**[Underestimation of Lipschitz constant]** In our experiments, we do not estimate the Lipschitz constant but rather we train an SCM that is Lipschitz-continuous by design and whose Lipschitz constant we can control. For a detailed description of this process and the overall model architecture, please refer to Appendix F2 in the supplementary.\n\n**[Comparison with \"Woulda, coulda, shoulda…\"]** We would like to clarify that “Woulda, coulda, shoulda: Counterfactually-guided policy search” by Buesing et al. [10] solves a different problem and thus it is incomparable with our method. In this context, a natural baseline to compare against would be the method introduced by Tsirtsis et al., NeurIPS 2021 [13] which solves a closely related problem in discrete state spaces, as we discuss in lines 52-54. Unfortunately, their method has a quadratic complexity with respect to the number of discrete states and thus does not scale to continuous multidimensional vector states as those used in our experiments. For example, discretizing the 9 continuous features we consider in our experiments into ten discrete levels each---a rather coarse-grained discretization---would lead to 1 billion discrete states.\n\n**[Monte Carlo anchoring]** The kick-starting approach we describe in lines 277-295 is an important part of our method because it generates the anchor set $\\mathcal{S}_ \\dagger$ required to compute the heuristic function $\\hat{V}_ \\tau$. Note that, in Figure 1(b), we vary the amount of Monte Carlo anchoring to investigate its effect on performance and, in Appendix E in the supplementary, we evaluate alternative selection strategies for the anchor set.\n\n**[Patient data]** We will release the code we used in our experiments but not the *anonymized* patient dataset (i.e., the MIMIC-III dataset). In this context, we would like to clarify that we are not the data owners, and we have gained access to the MIMIC-III dataset by submitting an application to the official owners (Physionet). As part of the process, we have passed a short online class in ethics. Please, refer to the MIMIC-III website for more details."}, {"review_id": "6QmE76r7hn", "reviewer": "Reviewer_DTwm", "summary": "The paper tackles the problem of finding optimal action sequences in domains with continuous state spaces with counterfactual reasoning.  They propose an A* based search approach and show the efficacy of the approach on a clinical sepsis management problem. ", "questions": "1. Graphical implication of an intervention on a variable $X$ (ie.  $do(X=x)$ ) is that all the incoming edges to $X$ are removed. In other words, as we are intervening on X, the influence from the $PA_X$  (parent nodes of X) on the node $X$ becomes zero. Refer to paragraph after Def 3.2.1 in Chapter 3 of Pearl. In Equation 3, it seems that all the arrows pointing out of $A_t$ are removed to adjust for $do(A_t = a_t)$ (as mentioned in line 106). So, the d-separation assumed does not follow from the do-calculus. \n2. Similar d-separation assumptions are made in Eq. 4 and 6.  \n\t1. Equation 1. suggests that two variables have a parent-child relationship: $S_t \\rightarrow A_t$.  \n\t2. Equation 2 suggests that their variable form a collider: $S_t \\rightarrow S_{t+1} \\leftarrow A_{t}$.\n\t3. On performing do-operation on $A_t$  would remove the parent-child relationship. But still maintain the collider  $S_t \\rightarrow S_{t+1} \\leftarrow A_{t}$. So, $S_{t+1}$ and $A_t$ are not d-separated by $S_t$. \n3. What are the implications of modifying the d-separation assumption on the proposed algorithms?\n4. How was the generated counterfactual action sequence evaluated? The SOFA score requires vital signs for calculations, how were the vital signs obtained?\n5. Figure 1 needs a legend. It is not clear what the pink & green lines represent.\n6. Did any physician analyze the counterfactual action sequences? Was any physian involved at any stage of the process?\n7. I recommend authors provide a reference table for notations in the appendix, it was really difficult to keep track of all the notations used in the paper.", "limitations": " ", "rating": 6, "confidence": 2, "soundness": 1, "presentation": 1, "contribution": 2, "strengths": "* The paper focuses on a very significant research problem\n* The experimental evaluation on a  real-world clinical dataset demonstrates the strong utility of the approach", "weaknesses": "* I have some concerns regarding the causal d-separation assumption (see Q1 & Q2).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their careful and insightful comments, which will help improve our paper. Please, find a point-by-point response below.\n\n**[d-separation]** Although the reviewer is correct that, in general, a do() operation on a variable (in our case, $A_t$) removes all *incoming* edges to that variable (in our case, $\\mathbf{S}_ t \\rightarrow A_t$), we respectfully disagree that Eq. 3 does not follow from the do-calculus. On the contrary, the first equality of Eq. 3 is a direct application of the 2nd rule of the do-calculus (action/observation exchange), as given in Theorem 3.4.1 of Pearl [8]. Specifically, the rule (applied to our case) states that one can write the (observational) probability $p^\\cal{C}(\\mathbf{S}_ {t+1} = \\mathbf{s} | \\mathbf{S}_ t = \\mathbf{s}_ t, A_t = a_t)$ as an interventional probability $p^{\\cal{C} ; do(A_t=a_t)} (\\mathbf{S}_ {t+1} = \\mathbf{s} | \\mathbf{S}_ t = \\mathbf{s}_ t)$ and vice versa if $\\mathbf{S}_ {t+1}$ is conditionally independent of $A_t$ given $\\mathbf{S}_ t$ in the graph resulting after deleting all *outgoing* edges of $A_t$. Note that this is equivalent to what lines 104-106 of our submission state—in the graph $A_t \\leftarrow \\mathbf{S}_ t \\rightarrow \\mathbf{S}_ {t+1}$, where the original outgoing edge $A_t \\rightarrow \\mathbf{S}_ {t+1}$ is deleted, $\\mathbf{S}_ t$ acts as a confounder and conditioning on it d-separates $A_t$ and $\\mathbf{S}_ {t+1}$.\n\nIn the revised version of our paper, we will include a figure illustrating the causal graph, and we will rephrase lines 104-106 to clarify this point. In light of the aforementioned explanation, we do not see any technical issues with Eqs. 3, 5 and 6, and we would be grateful if the reviewer could reconsider their score.\n\n**[Evaluation based on SOFA score]** As discussed in lines 344-346, the generated counterfactual action sequences were evaluated based on the counterfactual improvement they would have provided according to the SCM $\\mathcal{C}$, i.e., the relative decrease in cumulative SOFA score between the counterfactual and the observed episode. The SOFA score and the eight vitals required to compute it form the dimensions of the state vector, as stated in lines 305-307, and their counterfactual values are given by the (trained) SCM.\n\n**[Figure 1]** To avoid occlusions, we specify what the pink and green lines represent in the caption of Figure 1. The pink line represents the effective branching factor (EBF) and the green line represents the  average runtime of the A* search (in seconds). Here, we also used matching colors in the left and right y-axis of each panel to further indicate that the pink line corresponds to EBF and the green line corresponds to A* average running time.\n\n**[Evaluation by physicians]** The counterfactual action sequences were not analyzed by physicians for the purposes of this submission. In our work, we focus on formalizing and tackling algorithmically the problem of finding counterfactually optimal action sequences for episodes of sequential decision making processes with continuous states. However, as discussed in lines of 370-373, performing a user study with a systematic evaluation of counterfactual action sequences by human experts is a very interesting direction for future work.\n\n**[Notation]** Following the reviewer’s advice, in the revised version of our paper, we will include a reference table for notations in the Appendix."}, {"review_id": "qIEaxT5IAf", "reviewer": "Reviewer_qYbP", "summary": "This paper tackles the problem of finding counterfactual action sequences in sequential decision making problems. The main difference to previous methods is that this paper regards a continuous state space instead of a discrete one, which renders previous solution methods infeasible. The authors formalize the problem as an SCM and show that the solution in general is NP hard. As a solution, they present a new method based on the A-star algorithm, which yields for many problems an efficient solution. They evaluate the method on clinical data.", "questions": "How important are the assumptions that the SCM is bijective and Lipschitz-continuous? You need this assumption for the bounds in your algorithm. Can you still approach the problem if one condition is violated? How much do you think do the conditions restrict the field of applications of your method?\n\nAny thoughts about the conditions of the problem under which the method works efficiently?", "limitations": "Limitations and directions for future work are adequately addressed.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The paper is well-written and easy to follow.\n\nThe problem of efficiently finding counterfactual action sequences for the continuous state space seems new to me and relevant to the community.\n\nEven if the SCM formulation follows closely the work of Tsirtsis et al. (2021), a likewise dynamic programming approach would require enumerating all possible action sequences. From my point of view, the main contribution lies therefore in the A-star-based search method (and the anchor set selection), which requires further technical details such as bijectivity and Lipschitz continuity of the SCM for the heuristics. Therefore, the technical contribution seems good to me.\n\nThe method is evaluated with regard to the influence of different parameters.", "weaknesses": "The worse case complexity of the method is the same as for brute-force search. The efficiency of the method depends on the number of samples of the Monte Carlo samples in the anchor set. It is not clear to me, which conditions the problem must have in order for the method to work effiently and when it falls back to worse case complexity.\n\nMinor (typos):\nline 190 comma before \"such that\"", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their careful and insightful comments, which will help improve our paper. Please, find a point-by-point response below.\n\n**[Efficiency of our method]** The efficiency of our method depends on the tightness of the bounds $\\hat{V}_ {\\tau}$, which depends (i) on the number of the Monte Carlo samples in the anchor set, and (ii) on the Lipschitz constant of the SCM. Therefore, in general, our method will work more efficiently the lower the Lipschitz constant of the SCM, as shown in Figure 1(a). However, one cannot just decrease the Lipschitz constant arbitrarily to increase efficiency because this would degrade the goodness of fit of the SCM with respect to observational data. In Figure 4 in Appendix F, we investigate the goodness of fit of the SCM under different values of the Lipschitz constants.\n\n**[Bijectivity and Lipschitz-continuity]** The assumptions that the SCM is bijective and Lipschitz continuous are both necessary to conclude that the solution returned by our algorithm is optimal and that our algorithm is efficient. Since the bijectivity assumption is satisfied by many classes of SCMs studied in the literature, as discussed in lines 132-135, and the Lipschitz continuity is a natural assumption, as discussed in lines 122-131, we do not think these two assumptions significantly restrict the field of applications of our method."}, {"review_id": "QZd9pq2Cbx", "reviewer": "Reviewer_rmB3", "summary": "The authors present a method for finding c/f (optimal) action sequences in sequential decision making problems with uncertainty, with the novelty that they consider continuous state dynamics. They apply their method to the interesting setting of sepsis treatment at the end of the paper and therein demonstrate impressive results for their method (compared to the actual, observed, action sequences). ", "questions": "## Abstract\n\n- I re-read the abstract multiple times and I unfortunately struggle with the first half. Consider changing the sentence construction or just simplify the first half. It is very unclear what you are trying to say and I fear that the positioning of the paper suffers because of it. One suggestion could be to anchor the paper with one or two real examples and then take the justification from there.\n- The first paragraph of the introduction does a much better job of framing your paper, I would suggest adapting that approach for the first half of the abstract.\n\n## Introduction \n\n- Spell out ICU before you use it; it is common-enough term in the West but may be unfamiliar to a non-native speaker (or indeed anyone who does not frequently use English).\n- Excellent second paragraph.\n- Line 49: an SCM does not have the concept of state or transitions, merely a distribution P() over the endogenous and exogenous variables in the model (as defined originally by Pearl). Consequently how you have chosen to introduce it here does not quite follows. Consider re-phrasing. Moreover, SCMs are, as originally defined, state-less representations of a causal environment.\n- It would be helpful to the reader and yourselves if you could show us (the readers) a quadrant of work which has considered MDPs with: (continuous states, discrete actions), (continuous states, continuous action) … and so on — so that we understand what gap you are filling with this work (certainly you spell it out, but I would argue that a quadrant would be more useful and visually more powerful). \n- Line 60: please explain what a bijective SCM is. It is not enough to simply cite the original paper it if it is going to be key to your method.\n\n##\tA Causal Model of Sequential Decision Making Processes \n\n- Line 76: presumably $a_t \\in \\mathcal{A} = \\mathbb{Z}_+$?\n- How come you are not considered the discounted cumulative reward?\n- I think you should differentiate your setting more than what you are currently doing. At present it sounds as if eq (2) is part of the standard SCM definition. Perhaps introduce a definition paragraph for your SCM? That way the reader is clear that your are formalising a new concept and there is now ambiguity between yours and Pearl’s definition.\n- I would like to see a much longer discussion on why you make the causal sufficiency assumption (Line 97) in this work (no unobserved confounders). It is a common assumption as you say, but merely saying that other works make the same assumption, is not good enough justification to make it a tenable assumption. It is a huge assumption and highly unrealistic in most realistic scenarions. It makes life computationally far easier but trades-off usability down the line. Consequently, please discuss why you make this assumption. \n- Line 104: the do() express does not have to constitute a hard (atomic) intervention, it can also be a soft intervention (and others). Please consider or at least mention these settings too and why you settle for hard interventions.\n- Definition 1: is this definition somehow different from the standard one? Perhaps you ought to also give it a source since this is common fare but it may differ from the one readers are used to (and would also elucidate on the matter if it the standard definition or not).\n- The paragraph at the bottom of page 3 is great, very interesting application of Lipschitz continuity on a real problem.\n- Define the indicator function in equation 5.\n- Lines 149-161: I wonder if it may be not better for you to formally have a small lemma here, possibly even a theorem, explicating the results for which you have a proof. It currently reads as if it was merely a by-product of your method, whereas I imagine it is rather more important than that?\n\n## Problem statement\n\n- I would recommend revising the structure of your paper given that the problem statement appears on page four at conference with a page limit of nine. The paper will read far better if you frame the problem early on, allowing the reader to understand your angle of attack from the very start.\n- Line 163: you already introduced the episode on line 78.\n- Where is practically relevant to know that eq 9 is NP hard? A trivial question no doubt, but I would be interested to hear what the authors have in mind with regards to this theorem (is it really a theorem?)\n\n## Finding the optimal counterfactual action sequence via A* search\n\n- Line 210: I think a diagram here would be very helpful, explaining algorithm 1; a lot of page six could gainfully be explained with a figure rather than paragraphs of (rather dense) text.\n\n## Experiments using clinical sepsis management data\n\n- Looking at figure 2(a) I would have thought that as k increased we would have seen a much larger increase in the y-axis but this does not seem to be the case, the increase is fairly modest. Can you explain? \n- According to panel 2(b), 50% of patients see less than or equal to ~5% c/f improvement? Is that correctly read?\n- I enjoyed this section and the example is good but it is long and arduous, and we have no evidence to suggest that your method works for more (even synthetic) settings. It is an old trope of reviewing that we (reviewers) always want more experiments but in this case I think it is a fair request. Simply because I am struggling to understand how your ideas generalise across different settings and applications (even synthetic) nor do we have a comparison against other methods (are there any to which you can gainfully compare against?) \n", "limitations": "Please see the Questions section for the full review.", "rating": 6, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 3, "strengths": "Please see the Questions section for the full review.", "weaknesses": "Please see the Questions section for the full review.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their careful and insightful comments, which will help improve our paper. Please, find a point-by-point response below.\n\n**[Miscellaneous comments on presentation]** We would like to thank the reviewer for all their concrete suggestions regarding presentation/organization of the paper. We will address all of them in the revised version of the paper.\n\n**[Line 49]** We would like to clarify that we do not refer to states or transitions as concepts inherent to the original SCM definition. As lines 48-49 mention, we *use* an SCM (as originally defined by Pearl) to *represent* the stochastic state transitions of a sequential decision making environment. For more details, please refer to our response to **[SCM definition]**.\n\n**[Bijective SCM]** Please note that we formally define what a bijective SCM is in Definition 2 (lines 136-138).\n\n**[Line 76]** $\\mathcal{A}$ is a discrete *finite* set of actions, not the set of positive integers.\n\n**[Discounted cumulative reward]** Since we focus on a finite horizon setting, we think it is more meaningful to consider the undiscounted cumulative reward.\n\n**[SCM definition]** We would like to highlight that, in lines 82-84, where we first introduce the notion of an SCM, we do use Pearl’s (general) definition, i.e., that an SCM is consisted of four parts: (i) a set of endogenous variables (ii) a set of exogenous (noise) variables (iii) a set of structural equations assigning values to the endogenous variables, and (iv) a set of prior distributions characterizing the exogenous variables. Subsequently, in lines 84-92, we formally introduce a particular type of SCM $\\mathcal{C}$ that we use to represent sequential decision making processes—an instance of Pearl’s general definition. In our SCM, each state $S_t$ and each action $A_t$ in a finite sequence of actions and states $\\\\{S_t, A_t\\\\}_ {t=0}^{T-1}$ corresponds to a different endogenous variable (see lines 84-86), the variables $\\\\{Z_t, U_t\\\\}_ {t=0}^{T-1}$ are exogenous variables characterized by prior distributions $\\\\{P^\\mathcal{C} (Z_t), P^\\mathcal{C} (U_t)\\\\}_ {t=0}^{T-1}$ (see lines 88-89 & 91-92), and the value of each variable $A_t$ ($S_t$) is given by a structural equation $g_A$ ($g_S$), as stated in Equation 1 (2). In the revised version of our paper, we will rephrase lines 82-84 to clarify that they refer to the general definition of an SCM given by Pearl.\n\n**[No unobserved confounders]** We agree with the reviewer that the assumption of no unobserved confounding may be more or less realistic, depending on the exact domain and application at hand. In our work, our focus is on developing a method to solve the problem of finding counterfactually optimal action sequences in continuous state spaces. Since this problem has not been studied before, for simplicity, we have decided to tackle the problem under the assumption that there are no unobserved confounders, similarly as others in the literature have done when addressing a problem for the first time. However, we do not think that assumption nullifies our contribution, and we hope that future research will build upon our work to develop methods that work under unobserved confounders. Following the reviewer’s advice, in the revised version of our paper, we will include a separate paragraph at the end of Section 2, where we will discuss the implications of this assumption and highlight related work that focuses on environments where unobserved confounding exists (currently discussed in Appendix A in the supplementary).\n\n**[Hard interventions]** We settle on hard interventions because they fit better the problem formulation we consider. In the revised version of the paper, we will clarify that it would be interesting to consider variations of our problem using soft interventions.\n\n**[Definition 1]** To the best of our knowledge, the specific definition of Lipschitz-continuity for the type of SCMs we consider has not been used before.\n\n**[NP-hardness]** It is practically relevant to know that Eq. 9 is NP-hard because it tells us that there is no hope in designing a polynomial time algorithm that solves the problem we focus on. Since the problem has not been proven to be NP-hard before and the reduction we use to prove NP-hardness is novel, we think it is reasonable to include such a result as a Theorem. In the revised version of our paper, after Theorem 3, we will include a paragraph where we will discuss the importance and implications of the Theorem in the design of our algorithm.\n\n**[Figure 2(a)]** Figure 2(a) shows that, as $k$ increases, the marginal gains that could have been achieved in terms of the total counterfactual reward are diminishing. That finding follows our intuition since it indicates that, in retrospect, a small number of actions in each episode had the most significant effect on the episode’s outcome.\n\n**[Figure 2(b)]** The reviewer reads panel 2(b) correctly, 50% of patients see less or equal to ~5% c/f improvement. As discussed in lines 348-349, this indicates that the treatment choices made by the clinicians for most of the patients were close to optimal, even with the benefit of hindsight.\n\n**[Other experimental settings, other methods]** As an implication of Theorem 6, our method is guaranteed to *always* find the optimal solution to the problem we study, and we do not have reasons to expect that it would perform poorly in a different dataset or application. That said, we would like to share that we did consider performing synthetic experiments prior to the submission of the paper, and we chose not to go forward because we believe there are no additional insights one could gain about our method in a synthetic environment, other than the ones already presented in the paper. Due to the rebuttal's space constraints, we kindly ask the reviewer to refer to the response **[Lack of baselines]** to Reviewer **R23o** for more details re. other methods (or lack thereof) to compare our method against."}]

Author rebuttal:

Meta-review: The reviewers unanimously agreed on acceptance.
Decision: Accept (poster)