Grounding Predicates through Actions

Symbols representing abstract states such as"dish in dishwasher"or"cup on table"allow robots to reason over long horizons by hiding details unnecessary for high-level planning. Current methods for learning to identify symbolic states in visual data require large amounts of labeled training data, but manually annotating such datasets is prohibitively expensive due to the combinatorial number of predicates in images. We propose a novel method for automatically labeling symbolic states in large-scale video activity datasets by exploiting known pre- and post-conditions of actions. This automatic labeling scheme only requires weak supervision in the form of an action label that describes which action is demonstrated in each video. We use our framework to train predicate classifiers to identify symbolic relationships between objects when prompted with object bounding boxes, and demonstrate that such predicate classifiers can match the performance of those trained with full supervision at a fraction of the labeling cost. We also apply our framework to an existing large-scale human activity dataset, and demonstrate the ability of these predicate classifiers trained on human data to enable closed-loop task planning in the real world.

I. INTRODUCTION
Enabling robots to perform long horizon tasks such as preparing meals or assembling furniture is a widely studied problem. Long horizon planning is rooted in early AI work that studied how to give robots the ability to reason through symbols [1]. Symbols allow robots to abstract away low-level details of the environment and perform logical reasoning at a higher level [2][3][4][5]. However, giving robots the ability to perceive symbols in real-world environments is still an unsolved problem. Without some form of sensory grounding, propositions such as "drawer is open" are simply a set of symbols that lack any actionable meaning for the robot. Thus, robots often execute symbolic plans without closed-loop visual feedback-if a robot fails to open a drawer, it has no way of knowing, because it does not know what "drawer is open" looks like. State-of-the-art methods for learning visual groundings of symbols require large amounts of annotated data [3,[6][7][8]. However, obtaining annotations of symbolic states is prohibitively expensive due to the sheer number of propositions in a single image. Furthermore, densely-labeled datasets cannot easily be transferred between domains, since different planning problems often require different symbols. The   More examples can be found on our project website: https://sites. google.com/stanford.edu/groundingpredicates. from table" changes the propositions "hand is empty" and "cup is on table" to "hand is holding cup" and "cup is not on table"). Labeling a dataset with actions is easier than with symbolic states: it requires first defining the preand post-conditions once for each action class, and then annotating each visual example with only its action. Then, partial labels of symbolic state come for free with the action pre-and post-conditions. For task planning applications, the action pre-and post-conditions will already be defined in the task planning domain.
With this novel partial labeling scheme, we train networks to infer symbolic states in images in 20BN Something Something v2 (20BN) [9], a large-scale human activity dataset. The result is a system that can identify symbolic states in real-world environments for robot manipulation (Fig. 1). These classifiers open up many opportunities for long horizon planning in the real world, such as closed-loop task planning or learning from demonstrations of sequential actions.
The main contributions of this paper are three-fold. 1) We provide a framework for automatically labeling symbolic arXiv:2109.14718v2 [cs.RO] 4 Mar 2022 states in real-world image frames through actions, using the logical formalism of Planning Domain Description Language (PDDL) [10]. 2) We evaluate this framework on two domains, 20BN and Gridworld, and analyze its advantages and disadvantages. 3) We demonstrate the ability of predicate classifiers that are only trained on human demonstrations in 20BN to be used for closed-loop task planning in a real-world robot environment.

II. RELATED WORK A. Visual Relationship Detection
The computer vision community has a recently growing body of work on symbolic state detection, under the names visual relationship detection and scene graph generation [3,[6][7][8][11][12][13][14][15]. A scene graph represents image scenes as graphs where nodes are objects or attributes like grass and green, and edges are relationships between nodes, e.g., grass − green [16]. The aim of visual relationship detection is to enable semantic scene understanding and to connect visual concepts with natural language. The three most common datasets for visual relationship detection are CLEVR (100k synthetic images with 3 objects and 16 predicates) [16], VRD (5k real world images with 100 objects and 70 predicates) [6], and Visual Genome (100k real world images with 34k objects and 110k predicates) [17]. A key challenge of visual relationship detection (i.e., symbolic state detection) is the immense number of possible states. Even for the relatively small VRD dataset, 70 binary predicates with 100 objects results in 70 * 100 * 99 = 693k possible propositions. VRD and Visual Genome both rely on manual annotation, and thus obtaining fully labeled symbolic states is infeasible. Our method for labeling symbolic states requires defining the pre-and post-conditions once for each action class, and then partial symbolic state labels are automatically generated for all examples in the dataset.
It is important to note that our method only works when before and after images of the actions being performed are provided; VRD and Visual Genome do not meet this requirement. Furthermore, learning from pre-and post-conditions of actions only works in domains where the symbolic state can be manipulated by actions. Propositions such as "sky is blue" may be difficult to learn through actions, for example. The Action Genome dataset [18] introduces spatio-temporal scene graphs with the purpose of using symbolic state changes across actions to improve action recognition in videos. However, Action Genome, like VRD and Visual Genome, relies on manual annotations, and furthermore does not include negative propositions. Action Genome is based on the Charades dataset for video activity recognition [19]. Charades is one of numerous activity recognition datasets, like ActivityNet [20], Epic-Kitchens [21], and Atomic Visual Actions (AVA) [22], where labeling symbolic states with our method could be possible. However, we choose to evaluate our method on the 20BN Something Something v2 dataset [9], a video activity recognition dataset that focuses on the manipulation of everyday objects and has been shown to be useful for learning robotic manipulation skills [23,24].
The focus of this work is providing an alternative framework for labeling symbolic states, not proposing a novel model architecture for visual relationship detection. To evaluate our labeling method, we therefore train a model whose architecture is based on a state-of-the-art model for visual relationship detection [25]. However, our approach is agnostic to the specific model architecture.

B. Visual Grounding for Planning
A long-term goal for this work is to help bridge the gap between perception and long horizon planning in robotics. A majority of task planning methods operate solely at the level of symbolic abstractions, assuming that in a real world application, there would be some way to perceive symbols in the environment. Some works attempt to learn the symbols themselves, so that perception and planning can work together in an end-to-end fashion [2,26]. In these works, the learned symbols are not easily interpretable, and guaranteeing correct behavior over a long horizon for complex domains is therefore difficult. Our labeling method benefits frameworks that rely on direct supervision of symbols to learn to plan over long horizons, such as [4,[27][28][29][30]. These systems contain submodules that require annotations of predicates in images. This dependency restricts these systems to simulated environments where symbolic states are easily obtainable or to real-world domains small enough that manually annotating symbolic states is feasible.
The work most directly related to ours is perhaps an imitation learning system that learns to ground predicates associated with the pre-and post-conditions of demonstrated actions [31]. With grounded predicates, their system can then generalize to new tasks using classical planning. Specifically, their Visuospatial Skill Learning module learns to classify the predicate far(a) for a cuboid block a using demonstrations from two actions: push(a) : ¬far(a) → far(a) and pull(a) : far(a) → ¬far(a). While they learn to ground predicates from actions, the complexity of their problem is much simpler: one predicate. Our method enables learning to ground any number of predicates simultaneously from actions with arbitrarily complex pre-and post-conditions.

C. Weak Supervision
Our framework might be considered a form of weak supervision, which deals with classification problems where humans define annotation functions (e.g. pre-and post-conditions of actions) that provide noisy labels for unlabeled datasets to alleviate the effort of full labeling. Although our dataset is not strictly unlabeled (we assume action labels are available), methods from this relatively new field could be applied to our framework in future work [32,33]. In Appx. A, we provide back-of-the-envelope calculations to demonstrate how much annotation time weak supervision can save. Our weak supervision framework makes it easier to train predicate classifiers for custom task planning domains, which often contain symbols not transferable from other domains.

III. PLANNING DOMAIN DESCRIPTION LANGUAGE
This paper uses the Planning Domain Description Language (PDDL) [10] to describe symbolic domains. A PDDL domain can be specified with a tuple (Φ, A), where Φ is the set of predicates and A is the set of actions. A PDDL problem is a tuple (O, s init , g), where O is the set of environment objects, s init is the initial state, and g is the goal-specified as a first-order logic formula-to be satisfied.

A. Predicates and Propositions
A predicate has a fixed number of parameters, each of which can be instantiated with an object in O to form a proposition. For example, the predicate in(a, b) has two parameters a and b, and instantiating the parameters with arguments box and hand results in the proposition in(box, hand). In this paper, we use "positive proposition" to refer to atomic formulas (e.g., in(box, hand)) and "negative proposition" to refer to negations of atomic formulas (e.g., ¬in(box, hand)). Let P be the set of all possible positive propositions and N be the number of these propositions (N = |P|).

B. States
Symbolic states are conjunctions (∧) of propositions, e.g., in(box, hand) ∧ above(box, table) ∧ ¬on(box , table). PDDL follows the closed-world assumption, which means that propositions that are not explicitly specified are assumed to be false. Thus, we can also represent states as the set of all true propositions, where propositions not in the set are false by default. Let s denote closed-world states, and let S be the set of all possible states. Note that |S| = 2 N .
Under the open-world assumption, propositions that are not explicitly specified in the state are unknown-neither true nor false. Here, symbolic states can be represented as a pair of signed states (s + , s − ), where s + is the set of positive propositions, s − is the set of negative propositions, and propositions that are excluded from s + and s − are unknown. Letŝ denote open-world states represented in this manner.
In practice, we represent closed-world states as boolean vectors s ∈ {0, 1} N and open-world states as boolean matricesŝ ∈ {0, 1} 2×N , where the first and second row correspond to positive and negative states, respectively.

C. Actions
Actions are defined by their pre-conditions-a firstorder logic formula that must be true before performing the action-and post-conditions-a formula that is applied to the symbolic state after the action. These formulas can include universal (forall), existential (exists), and conditional (when) quantifiers in PDDL. A symbolic planner aims to find a sequence of actions that starts at the initial state s init and ends up in a symbolic state that satisfies the goal g.

IV. PARTIAL STATE LABELS FROM ACTIONS
Our goal is to use a PDDL specification to automatically label symbolic states in an entire dataset of videos only annotated with actions. In this section, we formalize our method for obtaining partial state labels from pre-and post-conditions of actions and discuss how to use these partial labels to train symbolic state classifiers. An overview of the pipeline is provided in Fig. 2.

A. Predicate Classification Network
Let M be the maximum number of parameters of any predicate in P. Given I as an RGB image of any dimension, and M bounding boxes (b 1 , . . . , b M ) of M objects representing ordered predicate arguments, the network outputs a vectorỹ of P probabilities, one for each predicate in P: The benefit of using object bounding boxes is that we can query the predicate network for propositions with specific argument ordering, such as in(spoon, cup) or in(cup, spoon). The predicate network would then output a prediction based on its understanding of what it looks like for argument a to be in argument b. We assume that most robot manipulation applications will already rely on object tracking for planning and control, and therefore obtaining object bounding boxes for predicate detection would not incur any additional cost. However, using bounding boxes is an implementation choice that is orthogonal to our proposed method of labeling symbolic states with pre-and post-conditions. A classifier that predicts propositions without bounding boxes could also be used in place of the predicate classifier for environments where bounding boxes are difficult to obtain.
Details of the network architecture that we use in the experiments can be found in Appx. C.

B. Obtaining Partial State Labels from Actions
The logic formulas that describe pre-and post-conditions can be arbitrarily complex with deeply nested compound formulas. To make these nested formulas usable for computing a neural network loss for symbolic state predictions, we first convert them via first-order logic algebra into a flattened form called disjunctive normal form (DNF). DNFs are written as disjunctions (∨) of conjunctions (∧) of positive and negative propositions. Each conjunction in a DNF can be interpreted as a partially-specified state under the open-world assumption.
With this interpretation, a pre-condition DNF represents a set of possible partial states before performing the action, and a post-condition DNF represents a set of possible partial states after. These candidate states are partial because actions only specify conditions relevant to the action itself. For  pred(a, b)). The predicate classifier outputs C predictions each forỹpre,ỹpost ∈ [0, 1] P , where P is the number of predicates in the symbolic domain. These predicate vectors are transformed into symbolic state vectors ypre, ypost ∈ [0, 1] N , where N is the number of propositions. The ground truth labels used to evaluate these predictions come from partial state labelsŝpre,ŝpost ∈ {0, 1} 2×N , obtained from the PDDL definition of the action pick(a) on the right. DNF cross entropy is used to compare the predictions to the partial state labels. At test time, a single image is fed into the network along with bounding boxes for the predicate arguments, and the predicate classifier predicts the symbolic state of the image. example, the action pick(box) may require ¬in(box, hand) as a pre-condition, but does not care whether open(door) is true and is not even aware of door's existence. The candidate states in a pre-or post-condition DNF can encode varying amounts of information about the open world (i.e., specifying anywhere from 1 to N propositions), and some candidates may be supersets of others. A DNF is satisfied if and only if at least one of the candidate states in the DNF is true (Eq. 5). In other words, when an action is performed on a concrete symbolic state, at least one candidate state in the pre-condition DNF will be true before performing the action, and at least one in the post-condition DNF will be true after. Although it is impossible to know which candidate is true without more information, we do know the truth value of the propositions that appear in all the candidate states.
To label symbolic states using pre-and post-conditions, we therefore collapse each DNF into the largest single partial state that satisfies all its conjunctions. This is equivalent to the intersection of all the positive and negative states.
If collapsing a pre-condition DNF results in the empty set, this means that the pre-conditions are too general to provide labels for symbolic state. Post-conditions with empty collapsed DNFs means that there is a conditional effect (when) that may result in the action producing no changes. Both cases can be avoided by making the pre-conditions more descriptive. For example, if the action close(a) contains the conditional post-condition is-closable(a) =⇒ closed(a), the action can be redefined so that is-closable(a) is a pre-condition and the post-condition is simply closed(a).
Letŝ pre (a) : A → S 2 andŝ post (a) : A → S 2 be functions that return collapsed pre-and post-condition DNFs, respectively, for action a. These functions provide a way to label propositions whose values are guaranteed to be known before or after an action is performed. Any other proposition could be either true or false without violating a pre-or post-condition, so nothing can be said about their ground truth values. A formal proof for this statement can be found in Theorem 1.1.

C. DNF Cross Entropy Loss
As summarized in Fig. 2, a data point in the training set consists of an action a along with its pre-and postconditions, a pair of images (I pre , I post ) corresponding to the action's before and after state, and bounding boxes of the objects whose predicates one wishes to identify. Let y pre ∈ [0, 1] N and y post ∈ [0, 1] N be the predicted outputs of the network. The ground truth labels derived from the pre-and post-conditions of a can be represented as a pair of matricesŝ pre (a),ŝ post (a) ∈ {0, 1} 2×N .
To measure how much a network prediction output agrees with an action's pre-and post-conditions, we define a modified cross entropy loss to handle DNFs. Here, σ(y) is the sigmoid function.
The loss function used to train the network on a pre-and post-condition pair is: loss(y pre , y post ;ŝ pre ,ŝ post ) = CE DN F (y pre ,ŝ pre )

V. EXPERIMENTS
In the following experiments, we apply our framework to a large-scale real-world dataset, where obtaining complete ground truth labels is impractical. We then apply these trained classifiers to perform closed-loop task planning in a real robot environment. Finally, we evaluate the effectiveness of using partial labels of symbolic state from action preand post-conditions compared to complete ground truth labels in a toy Gridworld environment.

A. Learning Predicates from Large-Scale Datasets
The 20BN dataset [9] contains 220, 847 video examples of 174 manipulation actions. Although 20BN does not label bounding boxes for the action arguments, we use the bounding boxes from Something-Else [34]. We have defined the pre-and post-conditions for all but two of the actions-"putting # of sth onto sth" and "stacking # of sth"-since the # variable cannot be easily described as a PDDL object. We have also defined 35 predicates relevant to the 172 actions. The full PDDL description can be found on the project website.
The purpose of this experiment is to show that our framework makes it possible to train symbolic state classifiers on large-scale real world datasets without relying on expensive manual annotations. Back-of-the-envelope calculations in Appx. A indicate that using pre-and post-conditions to partially label the symbolic states in all 132, 853 videos of our subset of 20BN would take an expert roughly four 8-hour work days. By contrast, obtaining the equivalent partial state labels with manual annotation would take 690 8-hour work days.
1) Setup: We train the predicate classifiers using DNF cross entropy (DNF CE), as well as a weighted version of cross entropy (DNF WCE) using Class-Balanced Loss [35] to overcoming imbalances in the predicate distribution. There is no oracle because we do not have ground truth labels of the full symbolic state in this real-world dataset.
2) Results: Both DNF CE and DNF WCE are able to learn the predicates specified in the 20BN PDDL with 0.96 train F1 and 0.92 test F1. Although the two achieve the same test F1 scores (within 0.001), the predicate F1 scores in Table I reveal that overall, DNF WCE performs better on the less frequent predicates, increasing the average predicate F1 score from 0.60 to 0.65.
Given that each proposition in 20BN is positive 36% of the time on average, a random classifier that outputs positive propositions with 0.36 chance would be expected to get 0.36 F1. Thus, an F1 score of 0.65 is significantly better than random. Both DNF CE and DNF WCE achieve 0.93 test accuracy. Fig. 1 shows qualitative examples of DNF WCE.
3) Transfer to Real Robot Domains: To demonstrate the zero-shot transferability of predicate classifiers trained with our framework to robot manipulation domains, we directly apply the DNF WCE classifier trained on 20BN to a robot pick-and-place setting for closed-loop task planning (Fig. 3).
For our real robot experiment, we create a pick-and-place environment where the robot needs to clear fruit off the table and store them in appropriate containers. We use the YOLOv5 object detector [36] for predicate bounding boxes along with Mediapipe object tracking [37] for smoother tracking between frames. The predicate classifier then sends a symbolic state prediction to a task planner that searches for a sequence of actions that satisfies the PDDL goal. The whole pipeline (from perception to task planning) operates at 10Hz and therefore allows closed-loop task planning that can respond to failures and disturbances. A video of this demonstration is provided in the supplementary material.  In this closed-loop task planning example, a user first gives the task planner a high-level goal in the form of a logical formula such as in(banana, drawer). The predicate classifier then predicts the symbolic state of the environment from RGB image data and sends it to a task planner. The task planner performs symbolic tree search to produce a sequence of actions to accomplish the goal, using 20BN actions such as "open sth" or "put sth into sth". As the actions are executed by low-level controllers, the task planner continues to re-plan based on updated symbolic state predictions given by the predicate classifier.

B. Partial vs. Full State Labels in Gridworld
Full symbolic state labels are ideal for training predicate classifiers. However, obtaining such labels is only practical for small domains or simulated environments with direct access to the symbolic state. To evaluate the effectiveness of training on partial vs. complete state labels, we use the Gridworld environment (Fig. 4), where we have full control of the symbolic state.
1) Gridworld Environment: In this environment, an agent needs to obtain a trophy from inside a chest. There are 8 objects: agent, trophy, chest, chest key, door, door key, room a, and room b. There are 6 predicates: The actions are defined with complex pre-and postconditions, making use of universal (forall), existential (exists), and conditional (when) quantifiers.
To generate training data, we first sample a random symbolic state s 0 , where each proposition in s 0 has a 5% chance of being true. There are 172 propositions in Gridworld, so on average, a randomly sampled state will have 8 propositions. Then, given an action with its pre-and postconditionsŝ pre andŝ post , we generate two symbolic states: We then render s pre and s post to obtain images I pre , I post .
2) Setup: We compare an Oracle trained on full state labels to two models. The first, DNF CE, is trained with DNF cross entropy (Eq. 7). Because partial labels contain less information than full labels, we expect DNF CE to take longer to train than Oracle and/or perform slightly worse.
The third model is an ablation study to test whether the model benefits from seeing the visual change induced by an action in order to visually ground symbols. This ablation, Half DNF, trains DNF CE on only the pre-or post-conditions of any given action instance, but not both. With normal DNF training, a single data point comes with a pair of images: one before and one after the action. For the ablation, we feed the model with one image from the pair: either before or after. Across all the data, the model will see examples of both the pre-and post-conditions for each action, but never both for a single action instance. Because this ablation only sees half the data, we double the size of its dataset for a fair comparison. We expect Half DNF to perform worse than DNF CE.
Oracle and DNF CE are trained on 10, 000 examples of actions, and Half DNF is trained on 20, 000. All of them are evaluated against the full ground truth state (as opposed to the partial DNF state) on the same test set of 10, 000 examples. Using the full state allows us to get a clearer picture of the generalization abilities of these models.
3) Results: As shown in Fig. 4, all the models achieve a training F1 score of nearly 1, indicating that they finish learning by 20  DNF's worse performance indicates that seeing the visual changes induced by each action is more beneficial than simply receiving more data (twice the amount).
Increasing the dataset size from 10, 000 to 100, 000 results in perfect test scores for all the DNF models for all the predicates (results not shown). Because DNF training uses partial labels, it requires more data to match the effectiveness of training on complete labels. However, in practice, acquiring a large dataset of actions should be easier than acquiring a small dataset annotated with complete symbolic states (Appx. A). Furthermore, in robotics applications where symbolic state classifiers might be used to determine whether the pre-or post-conditions of an action are satisfied, being able to classify propositions outside of the contexts provided by DNF labels may not be necessary.

VI. CONCLUSION
In this work, we have presented a framework for extracting partial symbolic states from action pre-and post-conditions, which can be used to label large datasets with less effort. This new method drastically reduces the cost of labeling large-scale real-world datasets (690 vs. 4 work days for 20BN). Yet, predicate classifiers trained with our method are still able to nearly match the performance of models trained with full ground truth labels, as shown with our Gridworld experiment. Our closed-loop task planning example demonstrates that predicate classifiers trained on large-scale real-world datasets can be applied to real robot domains. In many cases, predicates learned from a general large-scale dataset may not be applicable to custom task planning domains, where accuracy is critical. However, our labeling framework would perhaps be most useful for these very applications where collecting new datasets is necessary. For these situations, our framework makes obtaining real-world symbolic state classifiers much more feasible.
This work opens up many opportunities for long horizon planning in the real world. Other interesting avenues of work include using natural language to learn the pre-and post-conditions of actions [38,39], or combining visual groundings of symbols with natural language [3,40].


extfieldsofstudy:
Computer Science
Mathematics

text:
Counterfactual Explanations in Sequential Decision Making Under Uncertainty

Methods to find counterfactual explanations have predominantly focused on one step decision making processes. In this work, we initiate the development of methods to find counterfactual explanations for decision making processes in which multiple, dependent actions are taken sequentially over time. We start by formally characterizing a sequence of actions and states using finite horizon Markov decision processes and the Gumbel-Max structural causal model. Building upon this characterization, we formally state the problem of finding counterfactual explanations for sequential decision making processes. In our problem formulation, the counterfactual explanation specifies an alternative sequence of actions differing in at most k actions from the observed sequence that could have led the observed process realization to a better outcome. Then, we introduce a polynomial time algorithm based on dynamic programming to build a counterfactual policy that is guaranteed to always provide the optimal counterfactual explanation on every possible realization of the counterfactual environment dynamics. We validate our algorithm using both synthetic and real data from cognitive behavioral therapy and show that the counterfactual explanations our algorithm finds can provide valuable insights to enhance sequential decision making under uncertainty.

Introduction
In recent years, there has been a rising interest on the potential of machine learning models to assist and enhance decision making in high-stakes applications such as justice, education and healthcare [1][2][3]. In this context, machine learning models and algorithms have been used not only to predict the outcome of a decision making process from a set of observable features but also to find what would have to change in (some of) the features for the outcome of a specific process realization to change. For example, in loan decisions, a bank may use machine learning both to estimate the probability that a customer repays a loan and to find by how much a customer may need to reduce her credit card debt to increase the probability of repayment over a given threshold. Here, our focus is on the latter, which has been often referred to as counterfactual explanations.
The rapidly growing literature on counterfactual explanations has predominantly focused on one-step decision making processes as the one described above [4][5][6]. In such settings, the probability that an outcome occurs is typically estimated using a supervised learning model and finding counterfactual explanations reduces to a search problem across the space of features and model predictions [7][8][9][10][11]. Moreover, it has been argued that, to obtain counterfactual explanations that are actionable and have the predicted effect on the outcome, one should favor causal models [12][13][14]. In this work, our goal is instead to find counterfactual explanations for decision making processes in which multiple, dependent actions are taken sequentially over time. In this setting, the (final) outcome of the process depends on the overall sequence of actions and the counterfactual explanation specifies an alternative sequence of actions differing in at most k actions from the observed sequence that could have led the process realization to a better outcome. For example, in medical treatment, assume a physician takes a sequence of actions to treat a patient but the patient's prognosis does not improve, then a counterfactual explanation would help the physician understand how a small number of actions taken differently could have improved the patient's prognosis. However, since there is (typically) uncertainty on the counterfactual dynamics of the environment, there may be a different counterfactual explanation that is optimal under each possible realization of the counterfactual dynamics. Moreover, since, in many realistic scenarios, a decision maker needs to decide among a small number of actions on the basis of a few observable covariates, which are often discretized into (percentile) ranges, in this work, we focus on decision making processes where the state and action spaces are discrete and low-dimensional.
The work most closely related to ours, which lies within the area of machine learning for healthcare, has achieved early success at using machine learning to enhance sequential (clinical) decisions [15][16][17][18]. However, rather than finding counterfactual explanations, this line of work has predominantly focused on using reinforcement learning to design better treatment policies. A notable exception is by Oberst and Sontag [18], which introduces an off-policy evaluation procedure for highlighting specific realizations of a sequential decision making process where a policy trained using reinforcement learning would have led to a substantially different outcome. While our work builds upon their modeling framework, we do not focus on off-policy evaluation but, instead, on the generation of (better) alternative action sequences as a means of explanation for the process's outcome. Therefore, we see their contributions as complementary to ours.
More broadly, our work is not the first to use counterfactual reasoning in the context of sequential decision making in the machine learning literature [19,20]. However, previous work has used counterfactual reasoning either to learn optimal policies using limited observational data [19] or to explain the action choices of a reinforcement learning agent [20]. In contrast, our work uses a causal model of the environment to compute alternative action sequences, close to the observed one, which, in retrospect, could have improved the outcome of the decision making process. Refer to Appendix A for a discussion of further related work.
Our contributions. We start by formally characterizing a sequence of (discrete) actions and (discrete) states using finite horizon Markov decision processes (MDPs). Here, we model the transition probabilities between a pair of states, given an action, using the Gumbel-Max structural causal model [18]. This model has been shown to have a desirable counterfactual stability property and, given a sequence of actions and states, it allows us to reliably estimate the counterfactual outcome under an alternative sequence of actions. Building upon this characterization, we formally state the problem of finding the counterfactual explanation for an observed realization of a sequential decision making process as a constrained search problem over the set of alternative sequences of actions differing in at most k actions from the observed sequence. Then, we present a polynomial time algorithm based on dynamic programming that finds a counterfactual policy that is guaranteed to always provide the optimal counterfactual explanation on every possible realization of the counterfactual transition probability induced by the observed process realization. Finally, we validate our algorithm using both synthetic and real data from cognitive behavioral therapy and show that the counterfactual explanations can provide valuable insights to enhance sequential decision making under uncertainty 1 .

Characterizing Sequential Decision Making using Causal Markov Decision Processes
Our starting point is the following stylized setting, which resembles a variety of real-world sequential decision making processes. At each time step t ∈ {0, . . . , T − 1}, the decision making process is characterized by a state s t ∈ S, where S is a space of n states, an action a t ∈ A, where A is a space of m actions, and a reward R(s t , a t ) ∈ R. Moreover, given a realization of a decision making process τ = {(s t , a t )} T −1 t=0 , the outcome of the decision making process o(τ ) = t R(s t , a t ) is given by the sum of the rewards.
Given the above setting, we characterize the relationship between actions, states and outcomes using finite horizon Markov decision processes (MDPs). More specifically, we consider an MDP M = (S, A, P, R, T ), where S is the state space, A is the set of actions, P denotes the transition probability P (S t+1 = s t+1 | S t = s t , A t = a t ), R denotes the immediate reward R(s t , a t ), and T is the time horizon. While this characterization Figure 1: Structural causal model C for a Markov decision process M. Green boxes represent endogenous random variables and pink boxes represent exogenous noise variables. The value of each endogenous variable is given by a function of the values of its ancestors in the structural causal model, as defined by Eq. 1. The value of each exogenous noise variable is sampled independently from a given distribution. An intervention do(A t = a) breaks the dependence of the variable A t from its ancestors (highlighted by dotted lines) and sets its value to a. After observing an event S t+1 = s t+1 , S t = s t , A t = a t , a counterfactual prediction can be thought of as the result of an intervention do(A t = a) in a modified SCM where U t takes values u t from a posterior distribution with support such that s t+1 = g S (s t , a t , u t ).
is helpful to make predictions about future states and design action policies [21], it is not sufficient to make counterfactual predictions, e.g., given a realization of a decision making process τ = {(s t , a t )} T −1 t=0 , we cannot know what would have happened if, instead of taking action a t at time t, we had taken action a = a t . To be able to overcome this limitation, we will now augment the above characterization using a particular class of structural causal model (SCM) [22,23] with desirable properties.
Let C be a structural causal model defined by the assignments where U t and V t are n-and m-dimensional independent noise variables, respectively, and g S and g A are two given functions, and let P C be the distribution entailed by C. Then, as argued by Buesing et al. [19], we can always find a distribution for the noise variables and a function g S so that the transition probability of the MDP of interest is given by the following interventional distribution over the SCM C: where do(A t = a t ) denotes a (hard) intervention in which the second assignment in Eq. 1 is replaced by the value a t . Under this view, given an observed realization of a decision making process τ = {(s t , a t )} T −1 t=0 , we can compute the posterior distribution P C | St=st,St+1=st+1,At=at (U t ) of each noise variable U t . Building on the conditional density function of this posterior distribution, which we denote as f C | St=st,St+1=st+1,At=at Ut (u), we can define a (non-stationary) counterfactual transition probability where, in (a), we drop the do(·) because U t and A t are independent in the modified SCM. Importantly, the above counterfactual transition probability allows us to make counterfactual predictions, e.g., given that, at time t the state was s t and, at time t + 1, the process transitioned to state s t+1 after taking action a t , what would have been the probability of transitioning to state s after taking action a = a t if the state had been s at time t. Refer to Figure 1 for a visual description of our model and the notion of counterfactual predictions. However, for state variables taking discrete values, the posterior distribution of the noise may be nonidentifiable without further assumptions, as argued by Oberst and Sontag [18]. This is because there may be several noise distributions and functions g S which give interventional distributions consistent with the MDP's transition probabilities but result in different counterfactual transition distributions. To avoid these non-identifiability issues, we follow Oberst and Sontag and restrict our attention to the class of Gumbel-Max SCMs, i.e., where U t,s ∼ Gumbel(0, 1) and P (· | S t , A t ) is the transition probability of the MDP. More specifically, this class of SCMs has been shown to satisfy a desirable counterfactual stability property, which goes intuitively as follows. Assume that, at time t, the process transitioned from state s t to state s t+1 after taking action a t . Then, in a counterfactual scenario, it is unlikely that, at time t, the process would transition from a state s to a state s other than s t+1 -the factual one-unless choosing an action a that decreases the relative chances of s t+1 compared to the other states. More formally, for a model satisfying counterfactual stability, given it holds that P τ,t (S t+1 = s | S t = s, A t = a) = 0. In practice, in addition to solving the non-identifiability issues, the use of Gumbel-Max SCMs allows for an efficient procedure to sample from the corresponding noise posterior distribution P C | St=st,St+1=st+1,At=at (U t ), described elsewhere [18,24], and given a set of d samples from the noise posterior distribution, we can compute an unbiased finite sample Monte-Carlo estimator for the counterfactual transition probability, as defined in Eq. 3, as follows:

Counterfactual Explanations in Markov Decision Processes
Inspired by previous work on counterfactual explanations in supervised learning [7,8], we focus on counterfactual outcomes that could have occurred if the alternative action sequence was "close" to the observed one. However, since in our setting, there is uncertainty on the counterfactual dynamics of the environment, we will look for a non-stationary counterfactual policy π that, under every possible realization of the counterfactual transition probability defined in Eq. 3, is guaranteed to provide the optimal alternative sequence of actions differing in at most k actions from the observed one. More specifically, let τ = {(s t , a t )} T −1 t=0 be an observed realization of a decision making process characterized by a Markov decision process (MDP) M = (S, A, P, R, T ) with a transition probability defined via a Gumbel-Max structural causal model (SCM), as described in Section 2. Then, to characterize the effect that any alternative action sequence would have had on the outcome of the above process realization, we start by building a non-stationary counterfactual MDP M τ = (S + , A, P + τ , R + , T ). Here, S + = S ×{0, . . . , T −1} is an enhanced state space such that each s + ∈ S + corresponds to a pair (s, l) indicating that the original decision making process would have been at state s ∈ S had already taken l actions differently from the observed sequence. Following this definition, R + denotes the reward function which we define as r + ((s, l), a) = R(s, a) for any (s, l) ∈ S + and a ∈ A, i.e., the counterfactual rewards remain independent of the number of modifications in the action sequence. Lastly, let P τ be the counterfactual transition probability, as defined by Eq. 3. Then, the transition probability P + τ for the enhanced state space is defined as: if (a = a t ∧ l = l) ∨ (a = a t ∧ l = l + 1) 0 otherwise, ALGORITHM 1: It samples a counterfactual explanation from the counterfactual policy π Input: counterfactual policy π, horizon T , counterfactual transition probability P τ , reward function R, initial state s 0 .
where note that the dynamics of the original states s are equivalent both under P + τ,t and P τ,t , however, under P + τ,t , we also keep track of the number of actions differing from the observed actions. Now, let π : S + × {0, . . . , T − 1} → A be a policy that deterministically decides about the counterfactual action a t that should have been taken if the process's enhanced state had been s + t = (s t , l t ), i.e., the counterfactual state at time t was s t after performing l t action changes. Then, given such a counterfactual policy π, we can compute the corresponding average counterfactual outcome as follows: where τ = {((s t , l t ), a t )} T −1 t=0 is a realization of the non-stationary counterfactual MDP M τ with a t = π((s t , l t ), t) and the expectation is taken over all the realizations induced by the transition probability P + τ and the policy π. Here, note that, if π((s t , 0), t) = a t for all t ∈ {0, . . . , T − 1}, thenō π (τ ) = o(τ ) matches the outcome of the observed realization.
Then, our goal is to find the optimal counterfactual policy π * τ that maximizes the counterfactual outcome subject to a constraint on the number of counterfactual actions that can differ from the observed ones, i.e., maximize πō π (τ ) subject to where a 1 , . . . , a T −1 is one realization of counterfactual actions and a 1 , . . . , a T −1 are the observed actions. The constraint guarantees that any counterfactual action sequence induced by the counterfactual transition probability P + τ and the counterfactual policy π can differ in at most k actions from the observed sequence. Finally, once we have found the optimal policy π * τ , we can sample a counterfactual realization of the process and the optimal counterfactual explanation using Algorithm 1.

Finding Optimal Counterfactual Explanations via Dynamic Programming
To solve the problem defined by Eq. 7, we break the problem into several smaller sub-problems. Here, the key idea is to compute the counterfactual policy values that lead to the optimal counterfactual outcome recursively by expanding the expectation in Eq. 6 for one time step. We start by computing the highest average cumulative reward h(s, r, c) that one could have achieved in the last r steps of the decision making process, starting from state S T −r = s, if at most c actions had been different to the observed ones in those last steps. For c > 0, we proceed recursively: and, for c = 0, we trivially have that: with s ∈ S, r ∈ {1, . . . , T }, c ∈ {1, . . . , k}, and h(s, 0, c) = 0 for all s and c. In Eq. 8, the first parameter of the outer maximization corresponds to the case where, at time T − r, the observed action a T −r is taken and the second parameter corresponds to the case where, instead of the observed action, the best alternative action is taken. By definition, we can easily conclude that h(s 0 , T, k) is the average counterfactual outcome of the optimal counterfactual policy π * τ , i.e., the objective value of the solution to the optimization problem defined by Eq. 7, and we can recover the optimal counterfactual policy π * τ by keeping track of the action chosen at each recursive step in Eq. 8 and 9. The overall procedure, summarized by Algorithm 2, uses dynamic programming-it first computes the values h(s, 1, c) for all s and c and then proceeds with the rest of computations in a bottom-up fashion-and has complexity O(n 2 mT k). Finally, we have the following proposition (proven by induction in Appendix B): The counterfactual policy π * τ returned by Algorithm 2 is the solution to the optimization problem defined by Eq. 7.

Experiments on Synthetic Data
In this section, we evaluate Algorithm 2 on realizations of a synthetic decision making process 2 . To this end, we first look into the average outcome improvement that could have been achieved if at most k actions had been different to the observed ones in every realization, as dictated by the optimal counterfactual policy. Then, we investigate to what extent the level of uncertainty of the decision making process influences the average counterfactual outcome achieved by the optimal counterfactual policy as well as the number of distinct counterfactual explanations it provides.
Experimental setup. We characterize the synthetic decision making process using an MDP with states S = {0, . . . , n − 1} and actions A = {0, . . . , m − 1}, where n = 20 and m = 10, and time horizon T = 20. For each state s and action a, we set the immediate reward equal to R(s, a) = s, i.e., the higher the state, the higher the reward. To set the values of the transition probabilities P (S t+1 | S t = s t , A t = a t ), we proceed as follows. First we pick one s ∈ S uniformly at random and we set a weight w s = 1 and then, for the remaining states s ∈ S \ s , we sample weights It is easy to check that, for each state-action pair (s t , a t ) at time t, s t+1 = s is most likely to be observed in the next timestep t + 1. The parameter α controls the level of uncertainty.
Then, we compute the optimal policy that maximizes the average outcome of the decision making process by solving Bellman's equation using dynamic programming [21] and use this policy to sample the (observed) realizations as follows. For each realization, we start from a random initial state s 0 ∈ S and, at each time step t, we pick the action indicated by the optimal policy with probability 0.95 and a different action uniformly at random with probability 0.05. This leads to action sequences that are slightly suboptimal in terms of average   : Influence that the level of uncertainty of a synthetic decision making process has on the average counterfactual outcomeō π * τ (τ ) achieved by the optimal counterfactual policy π * τ as well as on the number of distinct counterfactual explanations π * τ provides. In all panels, we set n = 20, m = 10 and d = 1,000 and, in each experiment, use 500 realizations from 10 different instances of the decision making process (50 realizations per instance), each with different w s . In panel (c), for each realization, we sample 100 counterfactual realizations and compute the average number of unique explanations across realizations. Shaded regions correspond to 95% confidence intervals. outcome 3 . Finally, to compute the counterfactual transition probability P τ,t for each observed realization τ , we follow the procedure described in Section 2 with d = 1,000 samples for each noise posterior distribution 4 .
Results. We first measure to what extent the counterfactual explanations provided by the optimal counterfactual policy π * τ would have improved the outcome of the decision making process. To this end, for each observed realization τ , we compute the relative difference between the average optimal counterfactual outcome and the observed outcome, i.e., (ō π * τ (τ ) − o(τ ))/o(τ ). Figure 2 summarizes the results for different values of k. We find that the relative difference between the average counterfactual outcome and the observed outcome is always positive, i.e., the sequence of actions specified by the counterfactual explanations would have led the process realization to a better outcome in expectation. However, this may not come as a surprise given that the counterfactual policy π * τ is optimal, as shown in Proposition 1. In this context, note that, in the worst case, the counterfactual policy π * τ would trivially repeat the observed action sequence, leading to an average counterfactual outcome equal to the observed outcome with probability 1. Moreover, we observe that, as the sequences of actions specified by the counterfactual explanations differ more from the observed actions (i.e., k increases), the improvement in terms of expected outcome increases.
Next, we investigate how the level of uncertainty α of the decision making process influences the average counterfactual outcome achieved by the optimal counterfactual policy π * τ as well as the number of distinct counterfactual explanations π * τ provides. Figure 3 summarizes the results, which reveal several interesting insights. As the level of uncertainty α increases, the average counterfactual outcome decreases, as shown in panel (a), however, the relative difference with respect to the observed outcome increases, as shown in panel (b). This suggest that, under high level of uncertainty, the counterfactual explanations may be more valuable to a decision maker who aims to improve her actions over time. However, in this context, we also find that, under high levels of uncertainty, the number of distinct counterfactual explanations increases rapidly with k. As a result, it may be preferable to use relatively low values of k to be able to effectively show the counterfactual explanations to a decision maker in practice.

Experiments on Real Data
In this section, we evaluate Algorithm 2 using real patient data from a series of cognitive behavioral therapy sessions. Similarly as in Section 5, we first measure the average outcome improvement that could have been achieved if at most k actions had been different to the observed ones in every therapy session, as dictated by the optimal counterfactual policy. Then, we look into individual therapy sessions and showcase how Algorithm 2, together with Algorithm 1, can be used to highlight specific patients and actions of interest for closer inspection 5 . Appendix D contains additional experiments benchmarking the optimal counterfactual policy against several baselines 6 .
Experimental setup. We use anonymized data from a clinical trial comparing the efficacy of hypnotherapy and cognitive behavioral therapy [25] for the treatment of patients with mild to moderate symptoms of major depression 7 . In our experiments, we use data from the 77 patients who received manualized cognitive behavioral therapy, which is one of the gold standards in depression treatment. Among these patients, we discard four of them because they attended less than 10 sessions. Each patient attended up to 20 weekly therapy sessions and, for each session, the dataset contains the theme of discussion, chosen by the therapist from a pre-defined set of themes (e.g., psychoeducation, behavioural activation, cognitive restructuring techniques). Additionally, a severity score is included, based on a standardized questionnaire [27], filled by the patient at each session, which assesses the severity of depressive symptoms. For more details about the severity score and the pre-defined set of discussion themes refer to Appendix C.
To derive the counterfactual transition probability for each patient, we start by creating an MDP with n = 5 states and m = 11 actions. Each state s ∈ S = {0, . . . , 4} corresponds to a severity score, where small numbers represent lower severity, and each action a ∈ A corresponds to a theme from the pre-defined list of themes that the therapists discussed during the sessions. Moreover, each realization of the MDP corresponds to the therapy sessions of a single patient ordered in chronological order and time horizon T ∈ {10, . . . , 20} is the number of therapy sessions per patient. Here, we denote the set of realizations for all patients as T .
In addition, to estimate the values of the transition probabilities, we proceed as follows. For every state-action pair (s i , a), we assume a n-dimensional Dirichlet(α i1 , . . . , α i1 ) prior on the probabilities p j | i,a = P (S t+1 = s j | S t = s i , A t = a), where α ij = 1 if j ∈ {i − 1, i, i + 1} and α ij = 0.01 otherwise. Then, if we observe c j transitions from state s i to each state s j after action a in the patients' therapy sessions T , we have that the posterior of the probabilities p j | i,a is a Dirichlet(α i1 + c 1 , . . . , α in + c n ). Finally, to estimate the value of the transition probabilities P (S t+1 | S t = s i , A t = a), we take the average of 100,000 samples from the posterior probability p j | i,a . This procedure sets the value of the transition probabilities proportionally to the number of times they appeared in the data, however, it ensures that all transition probability values are  Figure 4: Performance achieved by the optimal counterfactual policy π * τ in a series of real manualized cognitive behavioral therapy sessions T , where each realization τ ∈ T includes all the sessions of an individual patient sorted in chronological order. Panel (a) shows the distribution of the relative difference between the average counterfactual outcomeō π * τ (τ ) achieved by π * τ and the observed outcome o(τ ), i.e., (ō π * τ (τ ) − o(τ ))/o(τ ), for k = 3 (refer to Appendix E for additional results under other k values). Panels (b) and (c) show the average counterfactual outcomeō π * τ (τ ) achieved by π * τ and the average number of unique counterfactual explanations provided by each π * τ , averaged across patients, against the number of actions k differing from the observed ones. In panel (c), for each realization, the average number of unique counterfactual explanations provided by π * τ is estimated using 1,000 counterfactual realizations. In all panels, we set d = 1,000 and use data from 73 patients. Shaded regions correspond to 95% confidence intervals. non zero and transitions between adjacent severity levels are much more likely to happen. Moreover, we set the immediate reward for a pair of state and action (s, a) equal to R(s, a) = 5 − s ∈ {1, . . . , 5}, i.e., the lower the patient's severity level, the higher the reward. Here, if some state-action pair (s, a) is never observed in the data, we set its immediate reward to R(s, a) = −∞. This ensures that those state-action pairs never appear in a realization induced by the optimal counterfactual policy. Finally, to compute the counterfactual transition probability P τ,t for each realization τ ∈ T , we follow the procedure described in Section 2 with d = 1,000 samples for each noise posterior distribution.
Results. We first measure to what extent the counterfactual explanations provided by the optimal counterfactual policy π * τ would have improved each patient's severity of depressive symptoms over time. To this end, for each observed realization τ ∈ T corresponding to each patient, we compute the same quality metrics as in experiments on synthetic data in Section 5. Figure 4 summarizes the results. Panel (a) reveals that, for most patients, the improvement in terms of relative difference between the average optimal counterfactual outcomē o π * τ (τ ) and the observed outcome o(τ ) is rather modest. Moreover, panel (b) also shows that the absolute average optimal counterfactual outcomeō π * τ (τ ), averaged across patients, does not increase significantly even if one allows for more changes k in the sequence of observed actions. These findings suggest that, in retrospect, the choice of themes by most therapists in the sessions was almost optimal. That being said, for 20% of the patients, the average counterfactual outcome improves a ≥3 % over the observed outcome and, as we will later discuss, there exist individual counterfactual realizations in which the counterfactual outcome improves much more than 3%. In that context, it is also important to note that, as shown in panel (c), the growth in the number of unique counterfactual explanations with respect to k is weaker than the growth found in the experiments with synthetic data and, for k ≤ 4, the number of unique counterfactual explanations is smaller than 10. This latter finding suggests that, in practice, it may be possible to effectively show, or summarize, the optimal counterfactual explanations, a possibility that we investigate next.
We focus on a patient for whom the average counterfactual outcomeō π * τ (τ ) achieved by the optimal policy π * τ with k = 3 improves 9.5% over the observed outcome o(τ ). Then, using the policy π * τ , also with k = 3, and the counterfactual transition probability P τ , we sample multiple counterfactual explanations τ using Algorithm 1 and look at each counterfactual outcome o(τ ). Figure 5(a) summarizes the results, which show that, in most of these counterfactual realizations, the counterfactual outcome is greater than the observed outcome-if at most k actions had been different to the observed ones, as dictated by the optimal policy,  Figure 5: Insights on the counterfactual explanations provided by the optimal counterfactual policy π * τ for one real patient who received manualized cognitive behavioral therapy. Panel (a) shows the distribution of the counterfactual outcomes o(τ ) for the counterfactual realizations τ induced by π * τ and P τ . Panel (b) shows, for each time step, how frequently a counterfactual explanation changes the observed action as well as the observed severity level and the severity level in the counterfactual realization with the highest counterfactual outcome. Here, darker colors correspond to higher frequencies and higher severities. Panel (c) shows the action changes in the unique counterfactual explanations (green) provided by π * τ along with the mean of counterfactual outcomes (r) that each one achieves and how frequently (f) they appear across the counterfactual realizations. Here, the bottom row shows the observed actions that were changed by at least one of the counterfactual explanations. Refer to Appendix C for a definition of the actions (i.e., themes). In all panels, we set d = 1,000 and the results are estimated using 1,000 counterfactual realizations.
there is a high probability that the outcome would have improved. Next, we investigate to what extent there are specific time steps within the counterfactual realizations τ where π * τ is more likely to suggest an action change. Figure 5(b) shows that, for the patient under study, there are indeed time steps that are overrepresented in the optimal counterfactual explanations, namely t ∈ {10, 13, 16}. Moreover, the first of these time steps (t = 10) is when the patient had started worsening their depression after an earlier period in which they showed signs of recovery. Remarkably, we find that, in the counterfactual realization τ with the best counterfactual outcome, the worsening is mostly avoided. Finally, we look closer into the actual action changes suggested by the optimal counterfactual policy π * τ . Figure 5(c) summarizes the results, which reveal that π * τ recommends replacing some of the sessions on cognitive restructuring techniques (CRT) by behavioral activation (BHA) consistently across counterfactual realizations τ , particularly at the start of the worsening period. We discussed this recommendation with one of the researchers on clinical psychology who co-authored Fuhr et al. [25] and she told us that, from a clinical perspective, such recommendation is sensible since, whenever the severity of depressive symptoms is high, it is very challenging to apply CRT and instead it is quite common to use BHA. Appendix F contains additional insights about other patients in the dataset.

Conclusions
We have initiated the study of counterfactual explanations in decision making processes in which multiple, dependent actions are taken sequentially over time. Building on a characterization of sequential decision making processes using MDPs and the Gumbel-Max SCM, we have developed a polynomial time algorithm to find optimal counterfactual explanations. Using synthetic and real data from cognitive behavioral therapy, we have shown that the counterfactual explanations our algorithm finds can provide valuable insights to enhance sequential decision making under uncertainty.
Our work opens up many interesting avenues for future work, which may solve some of its limitations. For example, we have considered sequential decision making processes with discrete states and actions satisfying the Markov property. Although this setting may fit a plethora of real-world applications, it would be interesting to extend our work to continuous states and actions and/or semi-Markov processes. Moreover, we experimentally validated our method using a single real dataset. It would be valuable to evaluate counterfactual explanations generated by our algorithm using additional datasets from other (medical) applications. In this context, it would be worth to consider applications in which the true transition probabilities are due to a machine learning algorithm. Finally, the usefulness of the counterfactual explanations given by our algorithm crucially depends on the particular structural causal model we have focused on. To make such explanations more practical, it would be important to consider alternative structural causal models and carry out a user study in which the counterfactual explanations are shared with the human experts (e.g., therapists) who took the observed actions.

A Further related work
Our work also builds upon further related work on interpretable machine learning and counterfactual inference. In terms of interpretable machine learning, in addition to the work on counterfactual explanations for one-step decision making processes discussed in Section 1 [4][5][6], there is also a popular line of work focused on feature-based explanations [28][29][30]. Feature-based explanations highlight the importance each feature has on a particular prediction by a model, typically through local approximation. In this context, one particular type of feature-based explanations that have been relatively popular in explaining the action choices of reinforcement learning agents in simple gaming environments (e.g., Atari games) are those in the form of saliency maps [31][32][33]. In terms of counterfactual inference, the literature has a long history [34], however, it has primarily focused on estimating quantities related to the counterfactual distribution of interest such as, e.g., the conditional average treatment effect (CATE). More broadly, one may also draw connections between our work and automated planning [35,36] and off-policy evaluation in reinforcement learning [37][38][39].

B Proof of Proposition 1
Using induction, we will prove that the policy value π τ ((s, l), t ) set by Algorithm 2 is optimal for every s ∈ S, l ∈ {0, . . . , k}, t ∈ {0, . . . , T − 1} in the sense that following this policy maximizes the average cumulative reward h(s, r, c) that one could have achieved in the last r = T − t steps of the decision making process, starting from state S T −r = s, if at most c = k − l actions had been different to the observed ones in those last steps. Formally: subject to Recall that, a 1 , . . . , a T −1 are the observed actions and the counterfactual realizations a 1 , . . . , a T −1 are induced by the counterfactual transition probability P + τ and the policy π. We start by proving the induction basis. Assume that a realization has reached a state s + T −1 = (s, l) at time T − 1, one time step before the end of the process. If c = 0 (i.e., l = k), following Equation 9, the algorithm will choose the observed action π τ ((s, l), t ) = a T −1 and return an average cumulative reward h(s, 1, 0) = R(s, a T −1 ) + s ∈S P τ,T −1 (s | s, a T −1 ) h(s , 0, 0) = R(s, a T −1 ), where h(s , 0, 0) = 0 for all s ∈ S. Since no more action changes can be performed at this stage, this is the only feasible solution and therefore it is optimal.
If c > 0, since h(s , 0, c) = h(s , 0, c − 1) = 0 for all s ∈ S it is easy to verify that Equation 8 reduces to h(s, 1, c) = max a∈A R(s, a) and π τ ((s, l), t ) = argmax a∈A R(s, a) is obviously the optimal choice for the last time step. Now, we will prove that, for a counterfactual realization being at state s + t = (s, l) at a time step t < T − 1 (i.e., r = T − t , c = k − l), the maximum average cumulative reward h(s, r, c) given by Algorithm 2 is optimal, under the inductive hypothesis that the values of h(s , r , c ) already computed for r < r, c < c and all s ∈ S are optimal. Assume that the algorithm returns an average cumulative reward h(s, r, c) by choosing action π τ ((s, l), t ) = a while the optimal solution gives an average cumulative reward OP T s,r,c > h(s, r, c) by choosing an action a * = a. Here, by τ t = {((s t , l t ) , a t )} T −1 t=t we will denote realizations starting from time t with a t = π τ ((s t , l t ) , t) where π τ is the policy given by Algorithm 2 and we will use τ * t if the policy is optimal. Also, we will denote a possible next state at time t + 1, after performing action a, as (s , l ) where l = l + 1 if a = a t , l = l otherwise and, c = k − l . Similarly, after performing action a * , we will denote a possible next state as (s , l * ) where l * = l + 1 if a * = a t , l * = l otherwise and, c * = k − l * . Then, we get: where, in (a), we expand the expectation for one time step, in (b), we replace the average cumulative reward starting from time step t + 1 with h(s , r − 1, c ) and OP T s ,r−1,c * for the policy of Algorithm 2 and the optimal one respectively and, in (c), we replace OP T s ,r−1,c * with h(s , r − 1, c * ) due to the inductive hypothesis.
It is easy to see that, it can either be a * = a t with c * = c or a * ∈ A \ a t with c * = c − 1. If c = 0, following Equation 9, the algorithm will choose the observed action (i.e., a = a t ). This is the only feasible solution, since a * = a t would give c * = −1 and l * = k − c * = k + 1, which is not a valid state. Therefore, we get a = a * = a t , which is a a contradiction. If c > 0, because of the max operator in Equation 8, for the action a chosen by Algorithm 2, it necessarily holds that: which is clearly a contradiction.

C Additional Details about the Cognitive Behavioral Therapy Dataset
Each patient's severity of depression is measured using the standardized questionnaire PHQ-9 [27], which consists of 9 questions regarding the frequency of depressive symptoms (e.g., "Feeling tired or having little energy?") manifested over a period of two weeks. The patient has to answer each question by placing themselves on a scale ranging from 0 ("Not at all") to 3 ("Nearly every day"). The sum of those answers, ranging from 0 to 27, reflects the overall depression severity and it is usually discretized into five categories, corresponding to no depression (0 − 4), mild depression (5 − 9), moderate depression (10 − 14), moderately severe depression (15 − 19), severe depression (20 − 27) 8 . In our experiments, the states S = {0, . . . , 4} correspond to these five categories. Each session of cognitive behavioral therapy contains information about the topic of discussion between the patient and the therapist, among 24 pre-defined topics [25], with some of the topics having similar content. For example, there were 4 topics about "cognitive restructuring techniques" which, we observed that, some therapists merged and covered in 2 sessions. Here, we grouped the above topics into the following eleven broader themes: • STR -First session: Introduction, discussing expectations, getting to know each other, discussing the current symptoms / problems, current life situation.
• BIO -Biography: A look at biography, family and social frame of reference, school and professional development, emotional development, partnerships, important turning points or crises.
• PSE -Psychoeducation: Discuss symptoms of depression, recognize and understand connections between feelings, thoughts and behavior (depression triangle) based on a situation analysis from the current / last episode, causes of depression, develop a disease model, explain the treatment approach in relation to the model.
• BHA -Behavioural activation: Focus on behaviour, discuss the vicious circle (depression spiral), discuss list of pleasant activities, attention to life balance, if necessary improve the daily structure, recognizing and eliminating obstacles and problems.
• REV -Review: Review of the last sessions, collection of strategies learned so far, find suitable strategies for typical situations, draft a personal strategy plan, plan further steps.
• CRT -Cognitive restructuring techniques: Discuss influence of thoughts on feelings and actions, identify thought patterns, discuss influence of automatic thoughts / basic assumptions, check the validity of automatic thoughts.
• INR -Interactional competence: Self-assessment of your own self-confidence, discuss current interpersonal issues and derive goals, carry out role plays, transfer into everyday life.
• THP -Re-evaluation of thought patterns: Review, evaluate and rename basic assumptions, schemes and general plans.
• RLP -Relapse prevention: Explain the risk of relapse, discuss early warning symptoms, recognize risk situations, develop suitable strategies.
• END -Closing session: Finding a good end to the therapy, looking back on the last 5-6 months, parting ritual.
• EXT -Extra material: Sleep disorders, problem-solving skills, brooding module "When thinking doesn't help", discuss the influence of rumination on mood and impairments in everyday life, progressive muscle relaxation.
In our experiments, the actions A correspond to these broader themes. However, since the themes STR and END appeared only in the first (t = 0) and last (t = T − 1) time steps of each realization, we kept them fixed and we did not allow these themes to be used as action changes during the time steps t = {1, . . . , T − 2}. Average counterfactual outcome Random Greedy Noisy greedy Algorithm 2 Figure 6: Performance achieved by the optimal counterfactual policy π * τ given by Algorithm 2 and the baseline policies in the same series of cognitive behavioral therapy sessions T introduced in Section 6. The plot shows the average counterfactual outcome 1 T τ ∈Tō πτ (τ ) achieved by π * τ and the baseline policies, averaged over the set of observed realizations T , against the number of actions k differing from the observed ones. For each observed realization, the average counterfactual outcome is estimated using 1,000 counterfactual realizations. Here, we set d = 1,000 and use data from |T | = 73 patients. Shaded regions correspond to 95% confidence intervals.

D Performance Comparison with Baseline Policies
Experimental setup. In this section, we compare the average counterfactual outcome achieved by the optimal counterfactual policy, given by Algorithm 2, with that achieved by several baseline policies. To this end, we use the same experimental setup as in Section 6, however, instead of setting R(s, a) = −∞ for every unobserved pair (s, a), we set R(s, a) = 5 − s ∈ {1, . . . , 5}, similarly as for the observed pairs. This is because, otherwise, we observed that there were always realizations under the baselines policies for which the counterfactual outcome was −∞. In our experiments, we consider with the following baselines policies: • Random: At each time step t, the policy chooses the next action a * uniformly at random if l t < k and it chooses a * = a t otherwise.
• Greedy: At each time step t, being at state (s t , l t ), the policy chooses the next action a * greedily, i.e., if l t < k, then and, if l t = k, a * = a t .
• Noisy greedy: At each time step t, being at state (s t , l t ), it chooses the next action a * as follows. If l t < k, a * is given by Eq. 12 with probability 0.5 and a * = a t otherwise. If l t = k, a * = a t .
Results. Figure 6 shows the average counterfactual outcomes achieved by the optimal policy, as given by Algorithm 2, and the above baselines for different k values. The results show that, as expected, the optimal policy outperforms all the baselines across the entire range of k values and, moreover, the competitive advantage is greater for smaller k values. In addition, we also find that the performance of the random baseline policy drops significantly as k increases, since, as discussed in Section 6, the observed trajectories are close to optimal in retrospect and, differing from them causes the random policy to worsen the counterfactual outcome.  Figure 7: Empirical distribution of the relative difference between the average counterfactual outcomeō π * τ (τ ) achieved by the optimal counterfactual policy π * τ and the observed outcome o(τ ), i.e., P[(ō π * τ (τ ) − o(τ ))/o(τ )] for several values of k. Here, we consider the same series of cognitive behavioral therapy sessions T introduced in Section 6. In all panels, we set d = 1,000 and use data from 73 patients.

E Average counterfactual improvement for additional values of k
In this section, we use the same experimental setup as in Section 6 to measure to what extent the counterfactual explanations provided by the optimal counterfactual policy π * τ would have improved each patient's severity of depressive symptoms, under other values of k. To this end, for each observed realization τ ∈ T corresponding to each patient, we compute the relative difference between the average optimal counterfactual outcome and the observed outcome, i.e., (ō π * τ (τ ) − o(τ ))/o(τ ). Figure 7 summarizes the results, which show that, similarly as in Section 6, for most patients, the improvement in terms of relative difference between the average optimal counterfactual outcomeō π * τ (τ ) and the observed outcome o(τ ) is modest. That being said, we observe that, as the sequences of actions specified by the counterfactual explanations differ more from the observed actions (i.e., k increases), the improvement in terms of expected outcome presents a slight increase.

F Insights About Individual Patients
In this section, we provide insights about additional patients in the dataset. For each of these additional patients, we follow the same procedure in Section 6, i.e., we use Algorithm 1 with the policy π * τ , with k = 3, to sample multiple counterfactual explanations τ and look at the corresponding counterfactual outcomes o(τ ). Figure 8 summarizes the results, where each row corresponds to a different patient. The results reveal several interesting insights. For most of the patients, all of the counterfactual realizations lead to counterfactual outcomes greater or equal than the observed outcome (left column), however, the difference between the average counterfactual outcome and the observed outcome is relatively small. Notable exceptions are a few patients for whom there is a small probability that the counterfactual outcome is worse than the observed one (top row) as well as patients for whom the difference between the average counterfactual outcome and the observed outcome is high (bottom row). Additionally, we also find that the actual action changes suggested by the optimal counterfactual policies π * τ are typically concentrated in a few time steps across counterfactual realizations (right column), usually at the beginning or the end of the realizations. Action changes, severity vs. t Figure 8: Insights on the counterfactual explanations provided by the optimal counterfactual policy π * τ for five real patients who received manualized cognitive behavioral therapy. Each row corresponds to a different patient with an observed realization τ . The panels in the left column show the distribution of the counterfactual outcomes o(τ ) for the counterfactual realizations τ induced by π * τ and P τ . The panels in the right column show, for each time step, how frequently a counterfactual explanation changes the observed action as well as the observed severity level and the severity level in the counterfactual realization with the highest counterfactual outcome. Here, darker colors correspond to higher frequencies and higher severities. In all panels, we set d = 1,000, k = 3, and the results are estimated using 1,000 counterfactual realizations.


Mean-Based Breakpoint Selection on Circular Histogram

Circular histogram represents the statistical distribution of circular data; theH component histogram of HSI color model is a typical example of the circular histogram. When using H component to segment color image, a feasible way is to transform the circular histogram into a linear histogram, and then, the mature gray image thresholding methods are used on the linear histogram to select the threshold value. *us, the reasonable selection of the breakpoint on circular histogram to linearize the circular histogram is the key. In this paper, based on the angles mean on circular histogram and the line mean on linear histogram, a simple breakpoint selection criterion is proposed, and the suitable range of this method is analyzed. Compared with the existing breakpoint selection criteria based on Lorenz curve and cumulative distribution entropy, the proposed method has the advantages of simple expression and less calculation and does not depend on the direction of rotation.

Introduction
e data obtained from actual observations can be expressed in various measurement spaces, and the angles' space showing the angle change is one of the measurement spaces. Angles' space data processing belongs to the branch of the discipline of statistics: direction (circular) statistics [1,2]. Angle-based data are called direction data, and angles are commonly expressed as unit vectors. Different from the measurement based on the scale, the direction data have inherent periodic (cyclic) characteristics, which make the direction data have many unique and novel characteristics in modeling and statistical processing. e data that show the angle change of a single variable is called circular data [2], and one of its visual display methods is a point on the unit circle or a unit vector on a plane. A typical example of circular data is the H component in the HSI color model of color image [3]. e HSI color model is a mathematical image model proposed by the American colorist H. A. Munsell in 1915; it uses H (hue), S (saturation), and I (intensity) to describe color characteristics. e HSI color model is different from the commonly used RGB color models. e three components R(red), G(green), and B(blue) of RGB color model are linearly dependent, but the three components H, S, and I of the HSI color model are linearly independent. Since the HSI color model has a good capability of representing the colors of human perception, the color image segmentation in HSI color space has achieved good results [4][5][6][7][8][9][10].
As a typical example of circular data, Hue(H) represents the basic colors of the image (H ∈ [0, 2π)); it can be expressed as a circular histogram, due to its periodicity. At first, some scholars used the hue histogram to segment the color image without considering its periodicity [11,12]. Tseng et al. [13] proposed the thresholding method of circular histogram for color image segmentation for the first time in order not to lose the periodicity of the hue. Wu et al. [14] gave an iterative Otsu's algorithm based on the circular histogram, but this method is not the optimal method and cannot guarantee the convergence of the algorithm. Dimov et al. [15] gave the method of optimal thresholding and multithresholding of circular histogram through the symmetry constraint of threshold point pairs, but the calculation is very complicated. Utilizing the cyclic characteristic of a circular histogram, Lai and Posin [16] theoretically analyzed that when the circular histogram is expanded into a linear histogram and the Otsu method is adopted, only half of the points on the circle need to be searched to obtain the optimal threshold point pair, which successfully reduces the time complexity from O(N 2 ) to O(N). However, this method is not general and can only be applied to two-class threshold.
Lai and Posin's research [16] shows that it is a feasible way that first expand the circular histogram into a linear histogram and then use the mature linear histogram threshold methods (such as Otsu's [17], fuzzy entropy [18], and context sensitive [19] thresholding techniques) to obtain the threshold of a circular histogram. How to choose a suitable breakpoint to linearize the circular histogram becomes a key. For this reason, we have proposed two breakpoint selection criteria. One is the criterion based on the Lorenz curve [20]; we discussed the relation between the area difference and the expansion direction and gave the optimal breakpoint selection criterion in the anticlockwise or clockwise direction. e other is the criterion based on the cumulative distribution entropy [21], and we built a circular histogram expansion model based on the cumulative distribution entropy and discussed the optimal breakpoint selection criteria under different expansion directions.
ese two circular histogram expansion methods overcome the randomness of breakpoint selection. However, the computational complexity of the Lorenz curve and cumulative distribution entropy is relatively high and makes the selection of the optimal breakpoint spend much time.
Circular statistics, as a particular branch, generally deals with data composed of angles or directions. Due to the obvious periodicity of this data, it is necessary to distinguish between circular data and linear data. In circular statistics, angle mean is used to represent the average angle of a set of data on a circle; it is a circular statistical invariant on the circular histogram; it does not change with the rotation of the circular histogram. On the contrary, line mean represents the average of a linear set of data and is a linear statistical invariant on the linear histogram. Since the angle mean is a circular statistical invariant, the line mean is a linear statistical invariant; in view of this, this paper proposes a simple breakpoint selection criterion to minimize the distance between the angle mean on the circle and the point on the circle corresponding to the line mean of the expanded linear histogram. e proposed criterion can quickly and reasonably find the breakpoint that keeps the distribution unchanged after the circular histogram is expanded.
is paper proposes a fast method for breakpoint selection in circular histogram, which solves the problem of low efficiency in expanding circular histogram into linear histogram. It is organized as follows. Section 2 describes the angle mean of the circular histogram and the line mean of the linear histogram and gives the optimal breakpoint selection criterion. In Section 3, the suitable range of the proposed method is given by comparing it with the optimal breakpoint selection criteria of the Lorentz curve based and cumulative distribution entropy based. Section 4 summarizes the paper.

Criteria for Selection of Breakpoint in Circular Histogram
In this section, we use the H component histogram of the HSI color model to explain the circular histogram. Figure 1(a) shows the H component diagram in the HSI color model. e H component (H ∈ [0, 2π)) represents the periodic change of color in the anticlockwise direction. For example, red is 0, green is (2π/3), and blue is (4π/3). Taking into account the periodic changes of the H component, a circular histogram is used to represent the statistical distribution of the H component (Figure 1(b)). When we use H component to realize color image segmentation, a feasible approach is to transform the circular histogram into a linear histogram, and then, we use the threshold segmentation methods on the linear histogram to select the threshold. e distribution information carried by different linearized histograms produced by the same circular histogram at different cutting points may be different. Figure 2 shows the result of the circular histogram ( Figure 1(b)) expanded at two different points. Although they are derived from the same circular distribution, their linearized distributions are not similar. In order to keep the distribution of the linearized histogram as consistent as possible with the distribution of the circular histogram, a new breakpoint selection method is given below.

Angles Mean and Linearized Mean of Circular Histogram.
(1) e angles mean μ L on circular histogram is defined as where μ ′ � tan − 1 (S/C). e average direction (μ L ∈ 0, 2π) given in definition (2) is a statistic that describes the position state characteristics of the circular histogram. It does not depend on the starting point and the rotation direction, reflecting the center of the circular histogram [1,2]. e red line in Figure 3 represents the angles mean of the circular histogram. To show more clearly, Figure 3 uses the rose diagram to illustrate the circular histogram.
Suppose the circular histogram h(t) { } ( Figure 1 e corresponding point of the line mean μ L (t 0 ) on the circular histogram is formulated as

Breakpoint Selection
Criteria. e goal of linearizing the circular histogram is to be able to maintain the complete original distribution. To find the optimal breakpoint, considering the angles mean μ L is a circular invariant on a circular histogram and μ L (t 0 ) is a linear invariant of the linear histogram expanded at the breakpoint t 0 , it is hoped that the point μ L (t 0 ) on the circle corresponding to μ L (t 0 ) and μ L are as close as possible so that the linear histogram μ L (t 0 ) and μ L are points on the circle. Because of periodicity, the distance between them is different from the Euclidean distance, and more attention is paid to the difference in the direction of the two values. e cosine value of the angle between them can be used to measure the difference in the direction of them. e distance between μ L (t 0 ) and μ L can be measured by the cosine of the angles [1,2] and expressed as e value of D is only related to the angles of μ L (t 0 ) and μ L , D ∈ [0, 2]. When the two angles are the same, D � 0; when the directions of two angles are opposite, D � 2.
e mean-based selection criterion for the optimal breakpoint t * 0 is Obviously, when the circular histogram expands in the clockwise direction, the corresponding value of the line mean on the circular histogram is the same as the value obtained in equation (4). erefore, the method in this paper is unrelated to the expansion direction of the circular histogram.
It is important to emphasize that the idea of mean-based breakpoint selection criterion is different from the existing Lorenz curve-based and cumulative distribution entropybased breakpoint selection criteria [20,21]. e mean-based method uses the invariants of circular statistics and linear statistics. Lorenz curve-based and cumulative distribution entropy-based methods, using the cumulative distribution information of each linearized histogram, are related to the counterclockwise or clockwise direction of rotation. e algorithm of breakpoint selection on circular histogram is very simple and easy to implement. e algorithm of breakpoint selection with mean is illustrated in Algorithm 1.

Experiment Results and Analysis
e experiment is divided into two parts to evaluate the proposed method. e experiments are performed using Python3.8 on a PC with Intel Core 2.50 GHZ CPU and 8 GB RAM, under Windows 10 operating system. In circular models, the von Mises distribution (also known as the circular normal distribution) is the most important distribution. e status is equivalent to the normal distribution in the linear distribution. Many theories with applications in the circular statistics are often discussed for the von Mises distribution [1,2]. erefore, the first part shows the results of selecting breakpoint for different types of artificial von Mises distributions and discusses the influence of parameters b (the mean direction) and κ (the concentration parameter) of the bimodal von Mises distribution [1,2] on the proposed method. In the second part, the proposed meanbased breakpoint selection criterion is compared with the existing breakpoint selection criteria, including the breakpoint selection method based on the Lorenz curve [20] (Lorenz-based), cumulative distribution entropy [21] (CDFE-based), and artificial bee colony [22] (ABC-based) on the H component circular histogram corresponding to 8 images from the Berkeley dataset. For convenience, the quantitative level L of the H component in the experimental part is 360.
To more specifically illustrate the linearization effect of the mean-based method on the mixture distribution of the same κ, Figure 6 shows the relation between the percentage of the broken distribution (see the red box in Figure 5(f )) and the mean direction difference (|b 2 − b 1 |). Due to the symmetry of the circle, the mean direction difference only is selected from 0 to 180. e effect of the rest is equivalent to its symmetrical part. When the mean direction difference is closer to 180, the proportion of the broken distribution will suddenly increase, but the maximum will not exceed 0.5%. When the mean direction difference is less than 150, the linearization effect is similar and better.
Similarly, Figures 7-9 show the relation between the percentage of the broken distribution and the mean direction difference when κ 1 and κ 2 are the combinations of (5, 10), (10,15), and (10, 20), respectively. e maximum percentage of broken distribution is 20% in Figure 7, and it is 14% in Figure 8, which shows the lower the overall concentration parameter, the worse the overall linearization effect.
In Figure 9, the maximum percentage of broken distribution is 16%. Figures 8 and 9 show the small difference in the concentration parameter of the two distributions is conducive to the linearization of the circular histogram.
In Figures 7-9, the percentage of broken distribution in different concentration parameters is positively related to the mean direction difference, and it increases exponentially around 180. is exponential increase greatly reduces the linearization effect near 180.  Figure 5: Different mixture bimodal histograms and corresponding mean-based linearized histograms: (a-c) circular histograms of the mixture distributions of g(t; 0, 10) with g(t; 60, 10), g(t; 120, 10), and g(t; 180, 10), respectively, and (d-f ) linearized histograms of (a-c) with the mean-based breakpoint selection method, respectively.
Input: H-histogram Hist, Hue magnitude L Output: e optimal breakpoint t * t * ←0 en calculate the distance D(t * ) with Hist according to equation (5) for t � 1: L magnitude do Rotate historium to the right or left by t Calculate the distance D(t * ) with Hist according equation (5) if (D(t * ) < D(t * )) then (D(t * )←D(t * )) t * ←t end if end for return t * ALGORITHM 1: e breakpoint selection with mean.

Mathematical Problems in Engineering
In summary, from the results of Figures 6-9, it can be seen that the mean-based method is suitable for situations where the target is not far from the background in the circular histogram. e distance between the mean direction should generally not exceed 5π/6.

Real Circular Histograms.
To further illustrate the scope and effect of the mean-based breaking method, the H component circular histogram corresponding to the 8 color images in the Berkeley dataset is selected for breakpoint selection and compared with the Lorenz-based [20], CDFEbased [21], and ABC-based [22] breakpoint selection criteria. e linearization result of 8 images can be seen in Figures 10-17.
e variance and kurtosis have also been computed to fully compare the effects of 4 algorithms. Tables 1 and 2 depict the variance and kurtosis using the existing Lorenzbased, CDFE-based, ABC-based, and proposed meanbased histogram techniques. Equation (7) defines the calculation formula of variance. Variance represents the discrete trend of data distribution. When the data distribution is relatively scattered, the variance is large, and when the data distribution is relatively concentrated, the variance is small. Equation (8) defines the calculation formula of kurtosis. e lower limit of kurtosis will not be lower than 1, and the upper limit will not be higher than the number of data. e greater the kurtosis, the steeper the distribution: where σ 2 is variance, K is kurtosis, X i is the ith value of the variable X, X is the average of the variable X, L is the number of the variable X, and p i is the probability of the ith value of the variable X.
As we can see from Figures 10 and 11, when the distribution is unimodal (bimodal coincidence), the breaking results of Lorenz-based, CDFE-based, ABC-based, and mean-based methods are appropriate, which guarantees the integrity of the distribution.
It can be seen from the H component histograms in Figures 12 and 13 that the distance between the centers of the two distributions of the circular histogram is small. In terms of maintaining the integrity of circular distribution, the CDFE-based method and the mean-based method shows the better effect. e variance and kurtosis of the mean method is the best in Tables 1 and 2.
Comparing the results of Figures 14(c)-14(f ), the CDFEbased method shows the best effect and is the Lorenz-based method. e distance between the centers of the two distributions of the circular histogram (Figure 14       consistent with the conclusion obtained by the artificial circular histogram analysis. When the centers of the two distributions are far apart, the linearization effect of the mean-based method will deteriorate.
From Figure 15, we can see the linearization effect of the CDFE-based and mean-based methods is better than that of the Lorenz-based and ABC-based methods. e linearization result using the Lorenz-based method does not Mathematical Problems in Engineering maintain the integrity of the circular distribution; a small part becomes the right part of the linearized histogram. e linearization result using the ABC-based method completely destroys the distribution. For complex distributions in Figures 16 and 17 that most color types appear, the frequency is different. It can be seen from Tables 1-2 that the mean-based method has a slight advantage over the ABC-based method, the CDFE-based method, and the Lorenz-based method in variance and kurtosis.
On the whole, from the linearization results of the 8 circular histograms, the mean-based method is superior to CDFE-based, Lorenz-based, and ABC-based methods, when the difference between the center positions of the target and the background is not particularly large. e mean-based method is effective in suitable scenarios. Judging from the average of metrics of the 8 images shown in Tables 1 and 2, the mean-based method is the best among the four methods. Table 3 shows the time spent on the linearization of the above 8 images by the Lorenz-based, CDFE-based, ABCbased, and mean-based methods, respectively. e meanbased method has a great advantage in speed. Compared with the Lorenz-based method, on an average, it can save about 7 times of time. Compared with the CDFE-based and ABC-based methods, the improvement is even greater, on an average, shortening about 71 times the time.

Conclusions
For the linearization of circular histograms, we propose a new method to select breakpoint. e new method uses a simple mean operation to give the optimal breakpoint selection criterion. We discuss the applicable scenarios of this breakpoint selection criterion. Experiments show that the new method can guarantee the linearization effect of the circular histogram in suitable scenarios, reduce the computational complexity of the breakpoint selection, and provide a better way for the linearization of the circular histogram. In future, we will explore a new breakpoint selection method based on the mean-based method in this paper, which have a better linearization effect when the centers of the two distributions are far apart.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.