Novel mathematical optimization Models for Group counterfactual Analysis

Programa de doctorado “Matemáticas” PHD DISSERTATION NOVEL MATHEMATICAL OPTIMIZATION MODELS FOR GROUP COUNTERFACTUAL ANALYSIS Author Miren Jasone Ramírez Ayerbe Supervisors Prof. Dr. Emilio Carrizosa Priego Prof. Dr. Dolores Romero Morales May 2024 Agradecimientos En primer lugar, y principalmente, tengo que darles las gracias por todo su apoyo y ayuda a mis dos directores, Emilio y Loli. Emilio, ya en la carrera, me animaste y me mostraste lo divertido, gratificante y emocionante que podía ser dedicarse a la investigación. Viste, antes que yo misma, que tenía las capacidades para desarrollar esta carrera. Te quiero dar las gracias por abrirme esa primera puerta, y por todas las horas de trabajo que han venido después. Por la rigurosidad y la originalidad de tus ideas. Por no olvidar nunca la parte humana. Loli, gracias por tu cercanía y por acogerme siempre, especialmente en Copenhague, que la siento ahora como una segunda casa. Tu pasión por este trabajo me inspira todos los días. Gracias por sacar tiempo de donde parece no haberlo, y por tus rigurosas correcciones. Gracias a los dos por cada conversación y cada debate. Por el apoyo en las etapas más duras. Si en un futuro, llego a ser una buena profesora y supervisora, será gracias a seguir vuestro ejemplo. Peter, working with you is always a pleasure. Thank you for your enthusiasm and for showing me new ways of looking at things and how to ask the right questions. I aspire to write in a manner as captivating and effective as you do. Tak. Andrea, thank you for welcoming me to New York and making it so easy to adapt to a new workplace. Thank you for giving me your time, even when it seemed nonexistent. I greatly enjoyed all of our talks throwing out ideas. Grazie. Gracias a todo el equipo de administración y secretaría del IMUS. Carmen, gracias por tu paciencia infinita y por hacer que coordinar cada estancia y cada viaje fuera lo más fácil posible. Gracias a todxs mis compañerxs del IMUS, a cada uno de ellxs. Estos cuatro años habrían sido mucho más duros si no hubiera podido compartir, no sólo las alegrías y los logros, sino, sobre todo, las frustraciones que son inevitables que ocurran en esta carrera de fondo. Gracias por el apoyo, las risas y el cariño. Cristina, no creo que hubiera podido terminar la tesis sin toda la ayuda que me has proporcionado. Siempre he dicho, que yo voy siguiendo el camino que ibas abriendo para mi. Gracias por la parte profesional, ayudándome con cada duda y papeleo, y por lo personal, por la alegría que contagias vayas a donde vayas. Nuria, me siento increíblemente afortunada de que empezáramos juntas en esto. Has sido mi apoyo más fundamental, mi compañera de viajes (y los que nos quedan). Gracias por estar II ahí, por entenderme y por las infinitas conversaciones. It won’t pass. Papá, mamá, sé que me diréis que no hace falta, pero gracias. III This research has been partially financed by research projects: • Optimización (FQM-329). Junta de Andalucía, Spain. • Data Science: from mathematical modelling to mathematical optimization (P18-FR2369). Junta de Andalucía, Spain. • Mathematical Optimization for Data-Driven Decision-Making (PID2019-110886RB-I00). MCIN/AEI/10.13039/501100011033, Spain. • Mathematical Optimization and Statistics for Explainable and Fair Machine Learning (PID2022-137818OB-I00). MCIN/AEI/10.13039/501100011033, Spain. • Mathematical Optimization in Data Science: Towards Interpretability in Classification, Regression and Dimensionality Reduction (US-1381178). Junta de Andalucía, Spain. • Research and Innovation Staff Exchange Network of European Data Scientists - NeEDS (H2020-EU RISE 822214). European Commission. • COST Action CA19130 - Fintech and Artificial Intelligence in Finance - Towards a transparent financial industry (FinAI). • Independent Research Fund Denmark (Grant 9038-00042 A) - “Benchmarking-based incentives and regulatory applications”. This support is gratefully acknowledged. V Resumen La inteligencia artificial y los algoritmos de aprendizaje automático son a menudo criticados por su falta de transparencia, siendo percibidos como “cajas negras” [Rudin, 2019]. Esta opacidad es especialmente indeseable en la toma de decisiones de alto riesgo, que llevan aparejadas decisiones críticas para los ciudadanos como la asignación de beneficios sociales, la concesión de préstamos, el diagnóstico médico o decisiones relacionadas con procedimientos legales [Azizi et al., 2018; Baesens et al., 2003; Zeng et al., 2022, 2017], con el peligro de producir resultados injustos para grupos en riesgo de sufrir discriminación [Besse et al., 2022; Miron et al., 2020, 2021]. La importancia de esta cuestión ha sido reconocida tanto por la industria como por las autoridades gubernamentales, incluida la Comisión Europea [European Commission, 2020; Goodman and Flaxman, 2017; Hupont et al., 2023]. En respuesta a esta necesidad, la Inteligencia Artificial Explicable (XAI) [Du et al., 2019; Goethals et al., 2022; Jung et al., 2020; Molnar et al., 2020; Rudin et al., 2022; Zhang et al., 2019] ha experimentado un notable auge, y ha crecido sustancialmente la producción de trabajos que desarrollan nuevos métodos para abordar estos retos. En la Inteligencia Artificial Explicable, los modelos de clasificación supervisada buscan un buen equilibrio entre precisión de predicción e interpretabilidad. Una vez que el clasificador ha sido entrenado, sería conveniente tener procedimientos para identificar cómo deberían modificarse las características de las instancias para ser clasificadas en la clase “buena”, por ejemplo, para ser clasificados como un buen pagador para un préstamo, o una persona sana para una prueba médica, o un profesional adecuado en un proceso de selección de personal, o un acusado clasificado como no reincidente. En esta tesis se estudian dichas soluciones modificadas, llamadas explicaciones contrafácticas, [Martens and Provost, 2014; Wachter et al., 2017]. Véase Artelt and Hammer [2019]; Guidotti [2022]; Karimi et al. [2022a]; Sokol and Flach [2019]; Stepin et al. [2021]; Verma et al. [2022] para revisiones bibliográficas recientes sobre Análisis Contrafáctico, y Browne and Swift [2020]; Freiesleben [2022]; Karimi et al. [2022a] para problemas relacionados que han aparecido en la literatura con nombres diferentes, como clasificación inversa o perturbaciones maliciosas. En esta tesis se estudian diferentes aspectos del problema de encontrar soluciones contrafácticas usando problemas de optimización. Más concretamente, se aborda el Análisis Contrafáctico Grupal, donde, en lugar de considerar una única instancia a la que se le asigna una VII sola explicación, consideramos un grupo de instancias para ser perturbadas y calculamos un conjunto de explicaciones para todas ellas. También estudiamos el problema cuando se tienen datos complejos y aplicamos el estudio más allá de la clasificación supervisada. A continuación, presentamos y motivamos los problemas abordados. El capítulo introductorio de esta tesis se basa en el trabajo de Carrizosa et al. [2024b]. Presentamos el problema de Análisis Contrafáctico Grupal, proporcionando un análisis riguroso sobre los elementos que definen este problema. También se identifican vínculos con problemas clásicos en la literatura de Análisis de Localización Continua y Discreta, como el de la mediana, el del centro, o el problema de la suma mínima de distancias al cuadrado, y destacamos los elementos novedosos, como las restricciones para garantizar que la probabilidad de ser clasificado en la clase positiva sea lo suficientemente alta. El Capítulo 2 se basa en tres trabajos diferentes, a saber, Carrizosa et al. [2024a], Carrizosa et al. [2024b], y Lodi and Ramírez-Ayerbe [2024]. En este capítulo, presentamos varios nuevos modelos de optimización que cubren todas las reglas de asignación posibles entre contrafácticos e instancias en el Análisis Contrafáctico Grupal. Primero, se presenta un modelo de asignación “one-for-one”, donde se incluyen diferentes restricciones que vinculan las explicaciones, como la imposición de dispersión (sparsity) global o restricciones de distancia entre las explicaciones contrafácticas. En segundo lugar, se formulan e ilustran modelos de optimización para reglas de múltiple asignación, a saber, “many-for-one”, “one-for-all” y “one-for-many”. Para este último, también se estudia una versión dispersa, donde se incluye un nivel adicional de interpretabilidad. Resolvemos el problema de minimizar el número de explicaciones contrafácticas necesarias para explicar un grupo de instancias. También se impone un número máximo de características que pueden cambiarse para todas las instancias que están vinculadas al mismo contrafáctico. Por lo tanto, no solo se detectan los contrafácticos representativos, sino también las características clave que necesitan cambiar en cada grupo. Se introduce un nuevo marco de generación de columnas diseñado para abordar eficazmente este problema, que, como se demuestra, tiene buenas propiedades de escalabilidad, rendimiento computacional y calidad de las soluciones. Los modelos se detallan para clasificadores lineales y bosques aleatorios, si bien el enfoque es válido para cualquier clasificador basado en puntuaciones (scores). El Capítulo 3 se basa en el trabajo de Carrizosa et al. [2023c]. En este capítulo, nos centramos en la naturaleza de los datos, específicamente desarrollamos un nuevo modelo de optimización matemática para datos funcionales, cuando se utiliza la regla de asignación “one-for- one”. El objetivo es identificar las muestras del conjunto de datos a partir del cual se construye la explicación contrafáctica, y cómo se deben combinar para que la instancia individual y su contrafáctico estén lo más cerca posible. Desarrollamos un modelo flexible que permite: incorporar varias medidas de distancia, incluida la popular distancia de deformación dinámica del tiempo (DTW); utilizar clasificadores basados en puntuaciones, específicamente bosques aleatorios; y lograr dos tipos de dispersión, en términos del número de atributos perturbados y en el número de muestras utilizadas para crear el contrafáctico. VIII Contents 1 Introduction 2 1.1 Ambientspace .................................. 6 1.2 Allocationrules.................................. 7 1.3 Constraints .................................... 8 1.4 Score-based classification models . . . . . . . . . . . . . . . . . . . . . . . . 10 1.5 Aggregating probabilities: Modelling P..................... 13 1.6 Modelling the cost criterion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2 Mathematical Models for Group Counterfactual Analysis with Different Allocation Rules 20 2.1 One-for-one allocation rule . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.1.1 Problem formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.1.2 Numerical illustration . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.2 Multiple allocations problems . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.2.1 A many-for-one allocation model . . . . . . . . . . . . . . . . . . . . 37 2.2.2 A one-for-many and one-for-all allocation model . . . . . . . . . . . . 38 2.2.3 A sparse one-for-many allocation rule model by column generation . . 39 2.3 Chaptersummary................................. 48 3 Counterfactual Analysis for Functional Data 70 3.1 Problemdefinition ................................ 73 3.2 AdditiveTreeModels............................... 75 3.3 Numericalillustration............................... 79 3.4 Chaptersummary................................. 82 4 Beyond classification: Counterfactual Analysis and Target Setting in Benchmarking 89 4.1 Introduction and motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 4.2 Background and literature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 4.3 TheSetting .................................... 97 4.4 Bilevel optimization for counterfactual analysis in DEA . . . . . . . . . . . . . 100 XV 4.5 Abankingapplication...............................106 4.5.1 Thedata..................................106 4.5.2 Counterfactual analysis of bank branches . . . . . . . . . . . . . . . . 107 4.6 Extensions.....................................112 4.6.1 Changing the returns to scale . . . . . . . . . . . . . . . . . . . . . . . 112 4.6.2 Changing the outputs . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 4.7 Chaptersummary.................................115 5 General Conclusions and Future Work 118 List of Figures 129 List of Tables 131 Appendix 133 References 138 XVI Chapter 1 Introduction 2 Finding counterfactual explanations amounts to solving a mathematical optimization model, whose ingredients will be detailed in this Chapter. While the most frequent version of Counterfactual Analysis found in the literature is the single-instance single-counterfactual case, this Chapter will deal with the more general case in which a group of counterfactuals is sought for a group of instances. We search for counterfactuals in a set X.We consider a binary classification problem, with labels in {−1,+1}(+1 being considered as the “good” class), and a classifier, identified by a function P:X −→ [0,1] such that P(x)is the probability of xbeing classified as positive, and thus if a deterministic classification is sought, the classifier labels as positive those xwith P(x)≥τfor some fixed threshold value τ, Kanamori et al. [2020]; Parmentier and Vidal [2021]. In the single-instance single-counterfactual case, we have at hand just one record x0∈ X where the classifier gives a low probability of being classified in the positive class, and we seek a so-called counterfactual instance xin the feasible set X(x0)⊂ X such that the cost C(x0,x) to perturb x0to yield xis low, while its probability P(x)of being classified as positive is high. In sum, in the single-instance single-counterfactual case, a counterfactual explanation for the instance x0can be found as an efficient solution of the following bi-objective optimization problem: min x∈X(x0)C(x0,x),−P(x).(CE) Needless to say, the same problem can be solved to change the prediction of a record from positive to negative. In this case, we have a record where the classifier gives a high probability of being classified in the positive class, and we seek a counterfactual explanation whose probability of being classified in the positive class is low, i.e., a high probability of being classified in the negative class. Different numerical solution approaches can be found in the literature to generate efficient solutions of Problem (CE). These include smooth optimization, e.g., Joshi et al. [2019]; Ramakrishnan et al. [2020], mixed integer optimization, e.g., Carrizosa et al. [2023c, 2024a]; Cui et al. [2015]; Fischetti and Jo [2018]; Kanamori et al. [2020, 2021]; Maragno et al. [2022]; Parmentier and Vidal [2021]; Russell [2019], multi-objective optimization, e.g., Dandl et al. [2020, 2023]; Del Ser et al. [2022]; Raimundo et al. [2022], robust optimization, e.g., Maragno et al. [2024], boolean satisfiability (SAT), e.g., Karimi et al. [2020], heuristic and metaheuristic approaches, e.g., Guidotti et al. [2019]; Poyiadzi et al. [2020]. In this thesis, we study the more general setting in which we are given a group of instances x0 s,s= 1,2...,S, hereafter denoted by x0= (x0 1,...,x0 S), to be perturbed to increase their probability of being classified in the positive class. Instead of finding one counterfactual instance for each instance x0 sindividually, we are interested in Group Counterfactual Analysis in which we seek for x0a group of Rcounterfactual instances x= (x1,...,xR)∈ X(x0)⊂ X:= XR. 4 From a stakeholder perspective, there are different reasons to perform a Group Counterfactual Analysis [Carrizosa et al., 2024a; Fernández et al., 2020, 2022; Rawal and Lakkaraju, 2020; Yousefzadeh and O’Leary, 2020; Yousefzadeh and O’Leary, 2022]. First, even if just one counterfactual is sought for each instance x0 s,s∈ {1,2, . . . , S}, linking constraints may exist, preventing solving the counterfactual models independently [Carrizosa et al., 2024a]. Such linking constraints appear, for example, when counterfactuals for records which are close should also be close, or when the statistical distribution of the counterfactuals should resemble the one of the original instances [Slack et al., 2021]. Second, several counterfactuals may be sought for the same instance, sufficiently far from each other [Wachter et al., 2017], and thus the procedure would yield for each instance a collection of counterfactuals that are hopefully diverse. Third, stakeholders may be in search of just a few counterfactual instances, to be seen as benchmarks for the group {x0 s}S s=1, and hence several instances will share the same counterfactual [Lodi and Ramírez-Ayerbe, 2024]. To end, stakeholders may want to detect a small subset of features such that perturbing these features can increase the probability P(x0 r)of being classified in the positive class, for all r, Eckstein et al. [2021]; Ley et al. [2023]; Piccialli et al. [2023]; Sharma et al. [2020]; Warren et al. [2023]. The Group Counterfactual Analysis yields the following optimization problem: min x∈X(x0)C(x0,x),−P(x),(GroupCE) where C(x0,x)measures the cost incurred when x0is perturbed to yield xand P(x)is defined as a component-wise nondecreasing function of the probabilities P(x)of the counterfactuals. To generate efficient solutions of a bi-objective optimization problem such as (GroupCE), two main methods are usually found in the literature of Multi-Objective Optimization, and they differ in the way the second objective Pis controlled. If imposed as a hard constraint, i.e., P must be above a threshold value ν∈[0,1], we obtain (GroupCEhard), min x∈X(x0)C(x0,x)(GroupCEhard) s.t. P(x)≥ν. (1.1) Alternatively, we may optimize a convex combination of both objectives, as done in Wachter et al. [2017] for the single-instance single-counterfactual case. However, it is well-known that, as opposed to (GroupCEhard), optimizing weighted sums may not generate all the efficient solutions of Problem (GroupCE), unless both the objective function and the feasible region are convex, assumptions which, as shown below, are unlikely to be fulfilled in this problem. See Carrizosa and Fliege [2002]; Carrizosa and Romero-Morales [2001]; Ogryczak [2001]; Ruiz et al. [2009] and references therein for other ways to obtain solutions to a bi-objective optimization model. 5 If we leave aside for a moment the Pcriterion, Problem (GroupCE) as well as its scalarized version (GroupCEhard) strongly resemble facility location problems, Drezner and Hamacher [2004]; Laporte et al. [2020]. Indeed, the Sinstances x0 1,...,x0 Sin x0may be seen as the set of users, and Rnew facilities, x1,...,xR, are sought in a region X(x0)to minimize the transportation costs from the users (the records) to the facilities (counterfactuals). The criterion P, concerning the probabilities of the counterfactuals being classified as positive, can be seen as related to the fixed cost of opening facilities. As for facility location problems, different types of optimization problems are obtained, depending on the choices of the ingredients defining Problem (GroupCE). We will describe the different ingredients required to define the mathematical optimization problems arising in Counterfactual Analysis, namely, the ambient space, i.e., the space where counterfactuals are taken from; the allocation rule, i.e., how counterfactual explanations are assigned to instances; the constraints imposed to build the counterfactuals; the classifier used and how the probabilities given by such classifier to the different counterfactuals to be in the positive class are aggregated; and the cost criterion measuring how difficult is for the instances to be perturbed to yield their counterfactuals. 1.1 Ambient space The first modeling decision affects the set X, hereafter the ambient space, containing each of the Rcounterfactuals. As in the literature in Location Analysis, in which this same question induces a taxonomy between Discrete Location, Mirchandani and Francis [1990], and Continuous Location, Plastria [1995], here we will have endogenous and exogenous counterfactuals. The set Xcan be a finite collection of observed datapoints, yielding the so-called endogenous counterfactuals, as in Ramon et al. [2020]; Wexler et al. [2019]. Using endogenous explanations has the advantage that the counterfactual explanations obtained actually exist in reality [Keane and Smyth, 2020]. Alternatively, and this is the most popular approach in the literature, Guidotti [2022], the counterfactuals can be synthetically built, yielding in this case the so-called exogenous explanations. When searching for exogenous explanations, it is common to assume that Xis a finite-dimensional vector space of dimension J, which is the case when the instances are characterized by numerical features. Categorical features are typically represented as a binary vector (in the so-called one-hot encoding, a categorical feature with v categories is represented by a binary vector in {0,1}v−1), and thus can also be seen as points in a vector space. For more complex types of data, a common practice consists of mapping the data into a vector space via a vector embedding [Van Looveren and Klaise, 2021], such as word embeddings for text data [Tolkachev et al., 2022] or node embeddings for graph data [Prado-Romero et al., 2022]. While such embeddings are very useful since they allow us to use the powerful machinery of vector analysis, one needs to revert to the original space (of words, graphs, etc.) the counterfactual solution obtained in the embedded vector space. Finally, and 6 to prevent the infinite-dimensional nature of the problems arising in the presence of functional features, the counterfactuals can be built as linear or convex combinations of observed datapoints, yielding thus a finite-dimensional X, making optimization much easier [Ates et al., 2021; Carrizosa et al., 2023c; Delaney et al., 2021; Karlsson et al., 2020]. Chapter 3 provides a detailed framework for conducting Counterfactual Analysis specifically tailored for functional data. In terms of the type of mathematical optimization problems that are obtained, finding endogenous counterfactual explanations amounts to solving combinatorial problems, with a very similar structure to discrete facility location problems, such as the classic p-median problem, Avella et al. [2007]; García et al. [2011]; Grötschel and Wakabayashi [1989]; Marín and Pelegrín [2019]; Mladenovi´ c et al. [2007], whereas searching for exogenous counterfactual explanations will lead to either continuous or mixed integer optimization problems, such as the classic minimum-sum-of-squared-distances problem, Aloise et al. [2012]; Liberti and Manca [2022]; Piccialli et al. [2022], the Weber problem, Chandrasekaran and Tamir [1990]; Plastria [1992]; Weiszfeld and Plastria [2009], or its extensions to the multi-facility case, namely, the multi-source Weber problem, Brimberg et al. [2000], and the minisum multi-facility Weber problem, Lefebvre et al. [1991]. This means that some of the existing numerical optimization tools developed in Discrete as well as Continuous Location Analysis can be easily tailored to tackle Group Counterfactual Analysis problems successfully, and, conversely, new facility location models are obtained motivated by this emerging field of application. 1.2 Allocation rules Most of the literature addresses the problem of finding one counterfactual explanation xfor one instance x0, as in Problem (CE). Even if it is assumed that there is a tuple of instances x0= (x0 1,...,x0 S), and a counterfactual for each individual is sought, Problem (CE) is solved separately for the different individuals. A user perspective is then followed. However, when a stakeholder perspective is followed, and a tuple xof counterfactuals is built for the tuple x0of individuals, the question of how to allocate instances to counterfactuals arises. (a) One-for-one (b) Many-for-one (c) One-for-all (d) One-for-many Figure 1.1: Illustration of the different allocation rules between instances (squares in red) and their counterfactual explanations (circles in blue) in Group Counterfactual Analysis. Consider a bipartite graph as the ones in Figure 1.1, where the origin nodes (squares in red) represent the instances, and the destination nodes (circles in blue) represent the labels for the 7 Rcounterfactuals. Observe that the value xr∈ X corresponding to each destination node is a decision variable. Moreover, in some cases, the edges of this bipartite graph are given, and in some cases, they will also be decision variables. Let us now introduce some notation for the edges of the bipartite graph linking instances to the labels of their counterfactuals. For each s∈ {1,2, . . . , S}, let Rsbe the set of indices r∈ {1,2, . . . , R}such that counterfactuals xr are associated with instance x0 s.Conversely, for each r∈ {1,2, . . . , R},let Srbe the set of indices s∈ {1,2, . . . , S}such that instances x0 sare associated with counterfactual xr.Observe that, by construction, r∈ Rsiff s∈ Sr.Depending on whether the edges in this bipartite graph are fixed or are decision variables, and depending on the geometry of such connections, different models are obtained. We briefly discuss four of them in what follows, which are easy to relate to the Location Analysis literature and are versatile enough to encompass the existing literature on Counterfactual Analysis. In the first allocation rule under consideration, for each instance, x0 sexactly one counterfactual xsis sought and vice versa, as in Figure 1.1a, i.e., Rsand Srare both a singleton and known in advance, yielding what we call the one-for-one allocation model, in which just the locations of the R=Scounterfactuals are sought. In the second allocation rule, the many-for-one allocation model, as in Figure 1.1b, all edges are known in advance too, satisfying that Sris a singleton for each r∈ {1,2, . . . , R}and each Rs=∅: each instance x0 shas associated its tuple (xr)r∈Rsof counterfactuals, Mothilal et al. [2020], fixing the concern raised in Wachter et al. [2017], where it is stated that one single counterfactual instance could be too restrictive and not take into account the user’s personal circumstances. In the third allocation rule, the one-for-all allocation model, as in Figure 1.1c, R= 1 and S1={1,2, . . . , S}: all instances x0 s,s= 1,2, . . . , S, share the same counterfactual x1, to be seen as a benchmark for the group. Therefore, in this allocation rule, all edges are known in advance too. Finally, in the fourth allocation rule, the one-for-many allocation model, as in Figure 1.1d, each Rsis a singleton and each Sr=∅. Here, the edges are decision variables too, and thus, one seeks both the location xrof the counterfactuals, r∈ {1,2, . . . , R}, and also a partition of {1,2, . . . , S}such that each subset {x0 s:s∈ Sr}shares the same counterfactual. Different models for all allocation rules are detailed in Chapter 2. 1.3 Constraints We discuss in what follows the type of constraints the counterfactual explanations x∈ X := XRfor the group x0are expected to satisfy. Most are valid for both models with endogenous and exogenous counterfactuals, whereas some require a certain structure (e.g., a vector space) on the ambient space X.Some can be found in the literature on single-instance singlecounterfactual models, and some appear naturally when a stakeholder’s perspective is followed. 8 [2019]; Plastria and Carrizosa [2001, 2012]. A gauge in RJis a function σ:RJ−→ R+ which is definite positive (σ(u)≥0∀u, σ(u)=0iff u=0), it is positively homogeneous (σ(τu) = τσ(u)∀τ≥0) and it is subadditive (σ(u+u′)≤σ(u) + σ(u′)).The unit ball of σ, i.e., the lower level set {u:σ(u)≤1}is a convex compact set, containing the origin in its interior, and indicates the moves which cost just 1unit. When σis absolutely homogeneous (σ(τu) = |τ|σ(u)∀τ),i.e., when the unit ball of σis symmetric with respect to the origin, σ is a norm. A natural choice for Dissimilarity(x0,x)is Dissimilarity(x0,x) = S X s=1 X r∈Rs ωsπ(σs(xr−x0 s)),(1.9) where ωs>0,πis a convex increasing function in R+,and σsis a gauge in RJ. Under these conditions, (1.9) is convex in x, and thus the objective function is convex if Rsis fixed and not a decision variable, as in the one-for-one, the many-for-one or the one-for-all cases. The convexity of πimplies that movements are penalized more than linearly. Using different weights ωsfor different instances x0 sallows stakeholders to ask for smaller perturbations for some individuals. This may be relevant in the framework of fairness, in which we may have records split into two groups, namely, those belonging to a sensitive group and the remaining ones, and a higher weight is given to individuals x0 sin the sensitive group. With respect to the choice of the gauges σs,norms, i.e., symmetric gauges, such as ℓp norms, are systematically used in the single-instance single-counterfactual models, Kanamori et al. [2020]; Russell [2019]; Wachter et al. [2017]. However, in practice, increasing in δunits one feature may not be as costly as decreasing δunits the very same feature, Karimi et al. [2021], and thus there is a need to depart from the state-of-the-art and address models with asymmetric gauges. Plausible asymmetric gauges for this problem with good structural properties are, nevertheless, easy to build. A well-known family of asymmetric gauges are the so-called skewed norms, Plastria [1992], namely, gauges σof the form σ(u) = σ0(u) + ηu,(1.10) where σ0is a norm, and ηis a fixed vector with σ0 0(η)<1,where σ0 0is the dual norm of σ0. Observe that σis asymmetric unless η= 0. See Plastria [1992] for a model through a skewed norm of the work expended when moving along an inclined plane, and Drezner and Drezner [2021]; Plastria [1992] for the effort of flying under steady wind conditions. Let us discuss two particular cases of skewed norms, namely the skewed ℓ1and ℓ2, which yield much more realistic models than those obtained with symmetric norms, and are as tractable as their symmetric counterparts. Taking in (1.10) as σ0the ℓ2norm in RJ,and ηsuch that σ0 0(η) = σ0(η)<1, the unit ball of σis an ellipsoid whose center is the origin only for η= 0 15 [Plastria, 1992]. Another plausible alternative in (1.10) is to take σ0as the ℓ1norm in RJ,and ηwith σ0 0(η) = maxj|ηj|<1, yielding σ(u) = J X j=1 |uj|+ J X j=1 ηjuj= J X j=1 ((1 + ηj) max(uj,0) + (1 −ηj) max(−uj,0)) ,(1.11) and thus, the cost of increasing δunits the feature jis δ(1 + ηj), while the cost of decreasing δunits the same feature is δ(1 −ηj).Observe that (1.11) is used as loss function in quantile regression [Yu et al., 2003]. Hence, taking in (1.9) πaffine and all σsin the form of (1.11), Dissimilarity(x0,·)is piecewise linear, and its optimization can be done, after adding auxiliary variables, by optimizing a linear objective. This way we can address with the same methods as in the weighted ℓ1 norm, Russell [2019]; Wachter et al. [2017], a much more realistic case in which perturbations increasing vs decreasing a feature are not equally costly. Another class of (asymmetric) gauges with good structural properties but seemingly unexplored in the literature of counterfactual analysis is obtained as an extension of (1.11), and called hereafter quantile gauges. Let us recall that a norm σ0in RJis said to be absolute if, for all u∈RJ,one has σ0(u) = σ0(|u|),where |u|is the vector (|u1|,...,|uJ|),[Bauer et al., 1961]. Observe that ℓpnorms (and their convex combinations) are absolute norms. Given an absolute norm σ0in RJand a vector η∈RJwith maxj|ηj|<1,define σηas ση(u) = σ0(|u|+Dηu),(1.12) where Dηis the diagonal matrix with ηin its diagonal. As before, observe that σηis asymmetric unless η= 0,and also it is a gauge. Note also that, if we take as σ0the ℓ1norm, we obtain the skewed norm in (1.11). Using in (1.9) each σsas a quantile gauge as (1.12), Dissimilarity(x0,·)is a sum of (increasing convex functions of) gauges, its minimization being a slight variant of the problem of minimizing a sum of ℓpnorms, for which Xue and Ye [1997, 2000] give polynomial time procedures based on the construction of logarithmically homogeneous self-concordant barrier functions. Moreover, if we take as σ0the ℓ2norm and π(t) = t2, one obtains Dissimilarity(x0,x) = (1.13) = S X s=1 X r∈Rs ωs J X j=1 (1 + ηj) max(xrj −x0 sj,0) + (1 −ηj) max(x0 sj −xrj,0)2. If the allocations defining Rsare fixed, as in the one-for-one, the many-for-one or the onefor-all models, the function in (1.13) is convex piecewise quadratic, and can be optimized as the squared Euclidean distance ℓ2 2has been in the single-instance single-counterfactual case in 16 the literature, but now addressing the asymmetry which penalizes differently the increase and decrease of the different features. Moreover, in the one-for-many rule, in which allocations (and thus Rs) are decision variables, one obtains (1.13) is a slight generalization of the minimumsum-of-squares-clustering problem, Aloise et al. [2012]; Liberti and Manca [2022]; Piccialli et al. [2022]. As a final remark on our discussion on quantile gauges, we should mention that they can be easily extended to accommodate causality relations between features [Karimi et al., 2022b; Mahajan et al., 2019; Pearl, 2009]. In its simpler form, assume that, if feature iis perturbed from xito xi+δ, then feature jis perturbed automatically (for free) from xjto xj+hijδ. In other words, paying for a perturbation δwould move the coordinates from x0to x0+δ(I+H), where Idenotes the identity matrix. Assuming I+His regular, this implies paying for δ:= (x−x0)(I+H)−1would move the point from x0to x,and thus in the cost function one should replace quantile gauges σsby σHs,defined as σHs(u) = σsu(I+Hs)−1.(1.14) The discussion above applies for the common case in which the ambient space Xis a J- dimensional vector space. When Xhas a different nature, other distance measures are to be used, see Kaufman and Rousseeuw [2009]. In the presence of mixed data, including numerical as well as categorical data, Gower distance, Gower [1971], a weighted sum of a distance measure on the numerical features and Hamming-type measures for categorical features is the most popular choice, see also Brughmans et al. [2023]; Wilson and Martinez [1997]. When Xconsists of time series and functional data, Esling and Agon [2012]; Xing et al. [2010], the dissimilarity between functions can be measured by some integrated Euclidean distance (easily extended to integrated quantile gauges), and, in case functions are inspected at a different speed, by using the Dynamic Time Warping distance [Carrizosa et al., 2023c]. If features represent frequencies of a discrete event, dissimilarities such as the chi-squared distance, Carrizosa et al. [2023a], can be used. When dealing with text data [Ramon et al., 2020; Tolkachev et al., 2022], an encoder may be used to map the instances into RJ, where counterfactuals are sought. The models of norms or gauges discussed above for data in RJmay be no longer meaningful to measure dissimilarities between mapped instances and counterfactuals, and instead, proximity measures such as the cosine similarity can be used. Once the counterfactual has been calculated in the embedding space, a decoder is needed in order to obtain the final counterfactual instance in the original space. Finally, for image data, see Vermeire et al. [2022] for an extensive survey. In addition to the dissimilarity, the cost function Ccan also capture how complex the perturbation to move from x0to xis. In this case, the cost function may have the form: C(x0,x) = Dissimilarity(x0,x) + λcComplexity(x0,x),(1.15) with λc>0. 17 When the ambient space Xis a J-dimensional vector space, complexity is usually measured through the number of features one needs to perturb, and it is, therefore, the complement of sparsity. Sparsity ensures more interpretability and it has been argued that people prefer explanations where fewer features are changed [Miller, 2019]. Sparsity can be modeled at the instance level or group level. For a given instance x0, this can be done directly through the ℓ0norm, ∥x0−x∥0, counting the number of features changed, or, in order to retain convexity, one can use the ℓ1norm instead. If sparsity at the group level is sought, one may minimize the total number of features changed for at least one instance in x0. In other words, we may minimize γ0defined as: γ0(x0,x) =    max s|x0 sj −xsj|J j=1   0 ,(1.16) where xsj denotes the value of feature jin xs. Notice that maxs|x0 sj −xsj|J j=1 is a vector of Jcomponents, where each component is the maximum change in feature jacross all the instances. Then, we count the number of features changed globally with the ℓ0norm. Note that more sophisticated forms of sparsity may be required in the presence of complex data, as in Carrizosa et al. [2017, 2021a, 2022b] for categorical features. When Xis defined by linear or convex combinations of observed datapoints, sparsity can be measured as the ℓ0norm of the vector of these coefficients, as in Carrizosa et al. [2023c]. In the extreme case, the counterfactuals built in this way are observed datapoints themselves. To have a smooth formulation of these measures of complexity, tractable with mixed integer optimization solvers, binary decision variables are introduced indicating whether a given feature can be perturbed, and constraints are expressed either via the usual big-M method or other strategies to address indicator constraints such as those in Belotti et al. [2016]; Bomze and Peng [2023]; Han et al. [2023]; Wei et al. [2022]. 18 Chapter 2 Mathematical Models for Group Counterfactual Analysis with Different Allocation Rules 20 In this Chapter we present several optimization models for Group Counterfactual Analysis, covering all possible allocation rules between counterfactuals and instances, namely the onefor-one, the many-for-one, the one-for-all and the one-for-many allocation rules. The remainder of the Chapter is organized as follows. In Section 2.1, we formalize a onefor-one allocation model, introducing the optimization problem that obtains counterfactual explanations for score-based classifiers. We detail this formulation for a specific cost function, accounting for the perturbations as well as the number of features changed, and for two families of classifiers, namely, a logistic regression model and an additive tree model. In Section 2.2, we introduce the optimization models for multiple allocation problems, specifically for the allocation rules many-for-one, one-for-all and one-for-many. In the case of the one-for-many allocation rule, two optimization models are developed: a first case where no sparsity is introduced and the number of counterfactuals sought is previously known, and a second case, where sparsity is included and the number of counterfactual explanations is minimized. For the latter, a column generation framework is developed. Several numerical illustrations for real-world datasets can be found in both sections. We end the Chapter with Section 2.3, where a summary is provided. 2.1 One-for-one allocation rule In this Section, we formalize the problem of finding counterfactual explanations when the one-for-one allocation rule is considered. We first discuss the main ingredients of the model, namely, the admissible perturbations, the way perturbations costs are measured, and the baseline classifier. We then address an easy but relevant case, namely, the separable case, which reduces the problem to solving a series of single-instance single-counterfactual problems. Then, for the non-separable case, a mixed integer programming formulation is given. Let us consider a binary classification problem, with classes denoted by +1,−1, in a set X ⊂ RJ. Having the data in RJdoes not mean that all features considered are numerical ones. Indeed, if we have categorical features, we can transform them into numerical ones using the standard one-hot encoding. For all the models described in this section we will use the probability criterion (1.7) as described in Section 1.5 from Chapter 1, to ensure that each counterfactual explanation has a high enough probability of being classified in the positive class. Thus, for score-based classifiers, see Section 1.4, xis assigned to the positive class if f(x)≥φ−1(ν), for a given threshold value ν. There may be some instances that do not behave like the rest of the group, and by perturbing them we would incur very high costs. These instances may be outliers to the group, and it is important to identify them to be treated separately. Whereas outliers’ identification is a well studied topic in data analysis, e.g. Boukerche et al. [2020], as far as the authors are aware, this issue has not been addressed so far in the literature on counterfactual explanations. As we 22 are dealing with groups of instances, we will allow for some of them to not be assigned an explanation, identifying thus outliers in the group of interest. Let us formalize these ideas and express the search of counterfactuals with the one-for- one allocation rule as a mathematical optimization problem. Let x0be the set of Sinstances, which form the group of interest, and let S∗⊂ {1, . . . , S}be the sought subset of cardinality S∗corresponding to the indices of the instances to be perturbed so that the classifier would classify them all in the positive class. With respect to the cost function used, although our approach is valid for more general settings as described in Section 1.6, we will consider in what follows a parametrized family of cost functions C: C(x0,x) = S X s=1 ∥x0 s−xs∥2 2+λind S X s=1 ∥x0 s−xs∥0+λglobγ0(x0,x).(2.1) The first term in (2.1) is an overall measure of the perturbations (measured through the Euclidean distance) to move from x0to x.Thus, we will consider in this case a symmetric gauge, the ℓ2norm, ωs= 1 for all instances and πas the quadratic function in (1.9). Regarding the complexity term, the second term in the cost function measures, weighted by λind ≥0, the overall number of perturbations in features needed to move from x0to x., and γ0is given by (1.16). Hence, the third term measures, weighted by λglob ≥0,the number of features perturbed for at least one instance to move from x0to x. 2.1.1 Problem formulation Given the score-based classifier induced by the score function f, the feasible set X0and the cost function Cas in (2.1), the mathematical optimization formulation of finding the minimal cost perturbations of the features to classify S∗instances in the positive class reads as follows: min x,S∗C(x,x0)(ColCE) s.t. f(xs)≥φ−1(ν)∀s∈ S∗(2.2) xs=x0 s∀s /∈ S∗(2.3) |S∗|=S∗(2.4) x∈ X0.(2.5) Let us first address an easy yet very relevant case. When λglob = 0,i.e., when the perturbation cost is measured by the sum of the squared Euclidean distances and the ℓ0norm between each x0 sand its corresponding individual counterfactual xs,then Cis separable. Moreover, when there are no linking constraints in X0, i.e., when X0=QS s=1 X0 s, then Problem (ColCE) is separable as well. This way, the problem can be split into Sindividual counterfactual problems, linked just by the constraint (2.4). This yields the following: 23 Proposition 2.1.1. Suppose λglob = 0,and X0=QS s=1 X0 s. For s= 1, . . . , S, let x∗ sbe an optimal solution to the single-instance single-counterfactual problem min xs ∥x0 s−xs∥2 2+λind∥x0 s−xs∥0(CEsep) s.t. f(xs)≥φ−1(ν)(2.6) xs∈ X0 s. Then an optimal solution for Problem (ColCE) is obtained by sorting the costs ∥x0 s−x∗ s∥2 2+ λind∥x0 s−x∗ s∥0and selecting the S∗instances with the lowest value of such cost. In what follows, we address the case in which the separability assumptions above do not hold. Let us express then (ColCE) in a more tractable form, by first rewriting the objective function as a quadratic convex function, and then rewriting the constraints (2.2)-(2.4), so that, when X0is polyhedral (eventually with some coordinates in integer numbers), a mixed integer convex quadratic program is obtained. In order to linearize the ℓ0norm and γ0in (1.16), binary variables are introduced. For every feature jand every instance s, we define the binary variable ξsj with ξsj = 1 if feature jof instance sis perturbed, i.e., xsj =x0 sj. The ℓ0norm is linearized by adding to the model the following constraints: −Msjξsj ≤x0 sj −xsj ≤Msjξsj j= 1, . . . , J, s = 1, . . . , S (2.7) ξsj ∈ {0,1},(2.8) where each Msj is a big constant, which exists thanks to the bounded nature of X0.In particular, any Msj ≥ |x0 sj|+ max z∈X0 s |zj|(2.9) is valid. A similar linearization can be done for the global sparsity γ0. In that case, a binary variable ξ∗ jis introduced. For each feature j,ξ∗ jtakes the value 1 if xsj takes a value different from x0 sj for some instance s= 1, . . . , S. The following constraints are added to our formulation: ξ∗ j≥ξsj s= 1, . . . , S, j = 1,...J (2.10) ξ∗ j∈ {0,1}j= 1, . . . , J. (2.11) Moreover, in order to indicate whether an instance is perturbed or not, new binary variables yare introduced. For s= 1, . . . , S,ys= 1 if instance sis allowed to be perturbed, and ys= 0 otherwise. Then constraints (2.2) and (2.4) are expressed respectively as ys(f(xs)−φ−1(ν)) ≥0s= 1, . . . , S. (2.12) 24 (a) Logistic regression (b) Random forest Figure 2.1: Counterfactual explanations for instances x0 sin Table 2.2 from the Boston housing dataset for two classifiers, the logistic regression and the random forest. The explanations have been calculated solving model (CEsep) with Cs(x0 s,xs) = ∥x0 s−xs∥2 2+ λind∥x0 s−xs∥0for λind = 0.02. The feature perturbations are displayed. classified in the negative class by the logistic regression model, i.e. that comply with wx0 s+b < 0,∀s= 1, . . . , S, (2.41) the key features that globally would need to change in order to change its class. We consider S∗=Sand ν= 0.5. To obtain the Pareto front, we will consider the global sparsity term as a hard constraint, instead of in the objective function. Thus, instead of solving (ColCELR), we will consider the following optimization problem: min x,y,ξ,ξ∗,u S X s=1 ∥x0 s−xs∥2 2 s.t. J X j=1 ξ∗ j≤Fmax (2.5),(2.7),(2.8) (2.10),(2.11),(2.13) −(2.20). (ColCELRhard) To describe the Pareto front, we solve this problem for all the possible values of the global sparsity that can be achieved, i.e., when Fmax takes values 2 to 13. Note that for Fmax = 1, the optimization problem is infeasible. The features that are changed are displayed in Figure 2.3. In this case, we can see a smooth transition in the increase of the perturbed features. Thirdly, we look at S∗, which can detect outliers from a group. Considering the same group as before, i.e., the instances that meet condition (2.41), we calculate the counterfactual 31 (a) Logistic regression (b) Random forest Figure 2.2: Counterfactual explanations for instances x0 sin Table 2.2 from the Boston housing dataset for two classifiers, the logistic regression and the random forest. The explanations have been calculated solving models (ColCELR) and (ColCEATM) with cost function (2.1) with λind = 0,λglob = 0.2. The feature perturbations are displayed. explanations for 95% of them, leaving 5%. We do this for cost function (2.1) with λind = 0 and two values of λglob = 0.1,10. The counterfactuals can be seen in Figure 2.4. In this case, for both values of λglob = 0.1,10 considered the same outliers are outlined. The index of the outliers, i.e., the instances that are not perturbed and therefore do not have a counterfactual associated, are 141, 373, 374, 384, 385, 386, 387, 388, 398, 412, 414, 438, 489, 490. Until now we have considered that the separability is destroyed because of the cost function, but it can also happen if we impose, for example, the Lipschitz continuity constraint through a dissimilarity d. Recall from Section 1.3 in Chapter 1 that this constraint had the following form: d(xi,xj)≤τd(x0 k,x0 l),∀i∈ Rk, j ∈ Rl,∀k, l ∈ {1,2, . . . , S},(2.42) for some threshold value τ. In the following, we provide some illustrations for the logistic regression model and a random forest with T= 100 trees, maximum depth 3 and wt=1 T, for all t= 1,...,100. All the parameters and hyper-parameters used in the following illustrations are detailed in https://github.com/jasoneramirez/GroupCE for its reproducibility. We now consider the case where the cost function is separable, with λind = 0.01 and λglob = 0, but the Lipschitz continuity constraint (2.42) with τ= 10 is imposed, linking the counterfactuals. This is calculated for instances x0 1and x0 2from Table 2.3, and for the random forest. The difference between imposing the Lipschitz continuity constraint or not is shown in Figure 2.5, where both the perturbations and feature values are displayed. The distance between the counterfactual explanations without imposing the Lipschitz continuity constraint is 0.51, whereas when the constraint is 32 Figure 2.3: Features that need to be perturbed (in red) for all the instances in the Boston housing dataset that comply with (2.41) when solving Problem (ColCELRhard) for all the values of Fmax to describe the Pareto front. The classifier considered is a logistic regression model. Fmax = 1 yields an infeasible problem. imposed the distance between explanations reduces to 0.46. sCRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX PTRATIO B LSTAT 1 0.001554 0.125 0.271628 0.0 0.286008 0.500287 0.959835 0.438387 0.173913 0.236641 0.276596 1.000000 0.480684 2 0.002303 0.125 0.271628 0.0 0.286008 0.396628 1.000000 0.450354 0.173913 0.236641 0.276596 0.974104 0.778146 3 0.001840 0.125 0.271628 0.0 0.286008 0.468097 0.854789 0.496731 0.173913 0.236641 0.276596 0.974305 0.424117 4 0.002457 0.125 0.271628 0.0 0.286008 0.539567 0.941298 0.474416 0.173913 0.236641 0.276596 0.988956 0.516556 5 0.001249 0.125 0.271628 0.0 0.286008 0.469055 0.823893 0.463503 0.173913 0.236641 0.276596 1.000000 0.318433 6 0.000983 0.125 0.271628 0.0 0.286008 0.446062 0.371782 0.392956 0.173913 0.236641 0.276596 0.983862 0.385762 7 0.007007 0.000 0.281525 0.0 0.314815 0.457559 0.606591 0.325355 0.130435 0.229008 0.893617 1.000000 0.180188 8 0.007099 0.000 0.281525 0.0 0.314815 0.485725 0.840371 0.303022 0.130435 0.229008 0.893617 0.957436 0.235375 9 0.006981 0.000 0.281525 0.0 0.314815 0.435524 0.552008 0.306359 0.130435 0.229008 0.893617 0.996772 0.185982 10 0.008743 0.000 0.281525 0.0 0.314815 0.465415 0.811535 0.284471 0.130435 0.229008 0.893617 0.974406 0.357064 Table 2.3: Feature values of x0= (x0 1,...,x0 10)from the Boston housing dataset used for the one-for-one models in Figures 2.5, 2.6 and 2.19. To end the illustrations for this dataset, we illustrate how the lower bound νimposed on the probabilities P(xs)affects the choice of the counterfactuals. In all the counterfactuals calculated before, the value of νimposed has been 0.5. We show in Figure 2.6a the Pareto frontier when the value of νchanges for the case where λind = 0,λglob = 0.1, the classifier is the logistic regression, and no other linking constraints have been imposed. The counterfactuals have been calculated for the 10 instances in Table 2.3. Figure 2.6b shows the features that need to be changed for the different values of ν. We can see how the larger the value of νis, the more features we perturb. Some of the experiments described for the Boston Housing dataset, have been carried 33 out in the other 3 additional datasets in Table 2.1. In these datasets, there are some nominal categorical features that have been transformed into numerical ones using the standard one-hot encoding, with one dummy variable for each category. Therefore, in the heatmaps each category is represented. However, one must take into account that a change from one category to another, although visualized as two dummy variables changes in the heatmaps (one is activated and changes from 1 to 0, whereas another is deactivated and goes from 0 to 1), in reality only one categorical feature is being perturbed. We will now describe the experiments for the COMPAS,Taiwanese Credit and Students performance datasets. COMPAS The COMPAS dataset was collected by ProPublica [Angwin et al., 2019] as part of their analysis on recidivism decisions in the United States. Following Karimi et al. [2020]; Mothilal et al. [2020]; Parmentier and Vidal [2021], 5 features are considered. The description of the dataset is detailed in Table 5.2 in the Appendix. First, counterfactual explanations for 10 individuals that were classified originally in the negative class by both classifiers, logistic regression and random forest, are calculated. The original values of the 10 instances can be seen in Table 2.4. Race AgeGroup Sex PriorsCount ChargeDegree African-American 2 0 7 1 African-American 2 0 16 1 African-American 2 0 18 1 African-American 1 0 3 1 African-American 1 0 1 1 African-American 1 0 3 1 Caucasian 3 0 28 1 African-American 3 0 38 1 Caucasian 3 1 28 1 African-American 2 0 19 1 Table 2.4: Feature values of x0= (x0 1,...,x0 10)from the COMPAS dataset used for the models in Figures 2.7-2.8. We solve models (CEsep), (ColCELR) and (ColCEATM), considering different combinations of the values of λind and λglob. For this dataset, the counterfactuals obtained also make sense. As one should expect, to not be labeled as recidivist, the number of prior criminal records has to decrease as well as the charge degree. Also, the COMPAS dataset has been used to exemplify unfairness in classification. This can be detected, in the separable case, for the random forest. Indeed, for one instance, if we were allowed to change their race from African-American to Caucasian, decreasing their charge degree would result in a positive classification. 34 Comparing the results in the separable case (Figure 2.7) with the non-separable case (Figure 2.8) we see that 4 features need to be changed globally (i.e., Race, AgeGroup, PriorsCount and ChargeDegree) in the separable case in the random forest case, whereas only 2 features need to be changed (AgeGroup and PriorsCount) in the non-separable case. For the logistic regression case, the perturbed features are the same in the separable and non-separable cases. Looking at the results, it can be seen that the features AgeGroup and PriorsCount are the features that have the greatest impact on changing the instances from the negative to the positive class for this group of instances and for both models. It is worth noticing that previous work [Dressel and Farid, 2018] about this dataset reached the conclusion that despite having the original dataset a collection of 137 features, the same accuracy could be achieved with a simple linear classifier with only two features. These two features where the age and the priors counts. Next, we display for all the instances in the COMPAS dataset that were classified in the negative class by the logistic regression model, i.e. that comply with (2.41), the key features that globally would need to be modified in order to change its class. We solve Problem (ColCELRhard) for Fmax = 1,2,3,4,5. For Fmax = 1, the optimization problem is infeasible. For the rest of the values of Fmax, the features that are changed are displayed in Figure 2.9. We reach the same conclusion as before: the only critical features to change the predicted class are AgeGroup and PriorsCount. Finally, for this dataset, we look at S∗, which can detect outliers from a group. The counterfactuals can be seen in Figure 2.10. For both values of λ= 0.1,10, the two features that are perturbed are the same. We obtain two optimal sets of outliers with the same total euclidean distance between the original samples and the counterfactuals. Taiwanese Credit The Taiwanese Credit dataset [Yeh and Lien, 2009] was used to predict the probability of default payments in Taiwan. The data was preprocessed as in Karimi et al. [2020]; Parmentier and Vidal [2021]. The description of the dataset is detailed in Table 5.3 in the Appendix. First, counterfactual explanations for 10 individuals that were classified originally in the negative class by both classifiers, are calculated. The original values of these 10 instances can be seen in Table 2.5. We solve models (CEsep), (ColCELR) and (ColCEATM), considering different combinations of the values of λind and λglob. For this dataset, to not be labeled as default, the number of total months overdue and the most recent payment amount has to decrease, whereas the number of months with zero balance has to increase. Comparing the changes in the separable case (Figure 2.11) with the non-separable case (Figure 2.12), we can see a big difference on the number of features needed to be perturbed. We see that one must change 4 features globally in the separable case, both in the logistic 35 regression and random forest case, whereas it is reduced to 2 and 1 feature in the non-separable case, respectively. Looking at the results, it can be seen that the critical features are different for each classifier for these 10 instances. For the logistic regression, the features that have the greatest impact are MostRecentPayment and TotalMonthsOverdue, whereas for the random forest is the feature about the number of months with zero balance over the last 6 months. Next, we display for all the instances in the Taiwanese Credit dataset that comply with (2.41), the key features that globally would need to to be modified in order to change its class. We solve Problem (ColCELRhard) for the different values of Fmax from 1 to 14. The features changed are displayed in Figure 2.13. We can see that the maximum number of features needed to be perturbed to change the prediction are 5. Finally, for this dataset, we look at S∗, which can detect outliers from a group. The counterfactuals can be seen in Figure 2.14. From the 144 outliers, 113 are in common for both cost functions. Students performance The Students performance dataset includes student grades, demographic, social and school related features and aims at predicting whereas a student will pass a course or not. The description of the dataset is detailed in Table 5.4 in the Appendix. First, counterfactual explanations for 10 individuals that were classified originally in the negative class by both classifiers, are calculated. The original values of these 10 instances can be seen in Table 2.6. We solve models (CEsep), (ColCELR) and (ColCEATM), considering different combinations of the values of λind and λglob. For this dataset, to be classified in the positive class, i.e., passing the course, the features that have to be increased are mostly the mother’s and father’s education, the study time, the free time, the quality of family relations and the health status; whereas the features that have to be mostly decreased are the age, the number of failures, the number of absences, travel time and how much they go out with friends. Comparing the changes in the separable case (Figure 2.15) with the non-separable case (Figure 2.16), one must perturb 7 and 15 features globally in the separable case in the logistic regression and random forest, respectively. This is reduced to 5 and 11 features in the nonseparable case. Focusing in the non-separable case ,we clearly see in the figures that the features that have the greatest impact are the age of the students, the mother’s education, the study time, the number of failures and amount of free time that the students have. In the case of the random forest there seems to be more features that impact on the decision made by the classifiers. 36 Besides the one already mentioned, the father’s education, the health status, the travel time and the number of failures are also perturbed for most records. Next, we display for all the instances in the Students performance dataset that comply with (2.41), the key features that globally would need to be perturbed in order to change its class. We solve Problem (ColCELRhard) for different values of the parameter Fmax,ranging from 1 to 30. For values 1 to 5, the problem was infeasible. The features changed are displayed in Figure 2.13. We can see that the maximum number of features needed to be perturbed to change the prediction are 14. Finally, for this dataset, we look at S∗, which can detect outliers from a group. The counterfactuals can be seen in Figure 2.18. For λglob = 0.1the indices of the outliers are 28, 57, 62, 69, 79, 85, 91, 121, 151, 174, 198, 225, 254, whereas for λglob = 10 the indices of the outliers are 27, 28, 31, 46, 50, 122, 124, 140, 146, 167, 174, 202, 216. There are only two in common. 2.2 Multiple allocations problems In this section, we focus now on the following allocation rules: the many-for-one, the one-for- all and the one-for-many allocation rules. 2.2.1 A many-for-one allocation model For the many-for-one allocation rule, we focus on the case where there are no linking constraints between counterfactuals associated with different instances. Thus, it boils down to the problem of calculating the counterfactuals for each instance, separately. When considering more than one explanation to an instance, some type of diversity is sought between the different explanations. We will illustrate one possible choice of diversity, see Section 1.3 for other types of diversity. Considering a specific feature of interest, we will impose that each of the counterfactuals has different values on this feature. In this way, we ensure that we have different options at hand to increase the probability, without relying on this feature to be modified to a specific value, that may be hard to obtain. The values allowed in each case for the considered feature are imposed in the set X(x0 s), that changes for each explanation. Thus solving the problem of finding Rcounterfactuals explanations for an instance, is equivalent to solving Roptimization problems as follows: min xr∈Xr(x0 s)∥x0 s−xr∥2 2+λind∥x0 s−xr∥0(2.43) s.t. f(xr)≥φ−1(ν),(2.44) where Xr(x0 s)corresponds to the feasible set for the r-th explanation associated to x0 s. 37 For the Boston housing dataset, we choose LSTAT as the feature of interest. We will calculate for each instance R= 3 counterfactuals, imposing LSTAT to be in the first quartile, between the first and third quartile, and above the third quartile, respectively. Specifically, we add to Xr(x0 s)one of the following constraints: (xr)LSTAT ≤Q1or Q1<(xr)LSTAT ≤Q3or (xr)LSTAT > Q3.(2.45) We illustrate this for the logistic regression model, the first two instances in Table 2.3 and R= 3. In Figure 2.19a the three counterfactuals for instance x0 1are displayed. One can see that imposing in the third explanation a high value for LSTAT, i.e., not allowing this feature to be decreased too much, results then in the necessity of moving other features to reach the desired probability of being classified in the positive class. The same happens for the counterfactuals for instance x0 2, displayed in Figure 2.19b. 2.2.2 A one-for-many and one-for-all allocation model We will now focus on the one-for-many and one-for-all allocation rules models. Whereas in the latter the only decision variable is the location of the counterfactual, in the former, both the location of counterfactuals and the assignment of counterfactuals to instances are to be decided. In this case, the constraints imposed on the explanations cannot depend on specific instances x0, thus X(x0) = X. Of course, there may be constraints imposed on all counterfactuals generally, and also, cannot link constraints could be added as well. As before, one way to define Pis as (1.7), ensuring in this way that each counterfactual explanation has a high enough probability. As cost function we take the ℓ2 2. To formulate the one-for-many allocation rule, binary variables ysr are introduced, where ysr = 1 if instance x0 sis assigned to counterfactual xrand 0otherwise. For a linear classifier such as the logistic regression model, we can formulate Problem (GroupCEhard) as follows: min x∈X,y R X r=1 S X s=1 ysr∥x0 s−xr∥2 2(2.46) s.t. wxr+b≥φ−1(ν)∀r= 1,2, . . . , R (2.47) R X r=1 ysr = 1 ∀s= 1,2, . . . , S (2.48) ysr ∈ {0,1} ∀s= 1,2, . . . , S ∀r= 1,2, . . . , R. (2.49) Constraint (2.47) guarantees that each counterfactual explanation has at least a probability of being classified in the positive class of ν, while constraints (2.48) and (2.49) ensure that each instance is assigned to exactly one counterfactual explanation. To solve this problem an alternating algorithm can be used, similar to Lloyd’s algorithm 38 [Lloyd, 1982], where two phases arise: the allocation of the instances to the counterfactual instance that minimizes the cost function, and the location of the explanations, where the counterfactual explanation is calculated in the case where the clusters are already known. Since we are dealing with exogenous counterfactuals, the chosen cost function is the squared distance and we are considering a logistic regression model, this problem has a very similar structure to the classical minimum-sum-of-squared-distances problem, with the addition of linear constraints (2.47). The problem becomes more complex when, for instance, other classifiers are chosen, such as an additive tree model. Setting R= 1 in (2.46)-(2.49), one obtains the model for the one-for-all allocation rule. In such case, the variables ysr are not needed, as all the instances are assigned to the same counterfactual explanation, and there is no allocation phase. To visualize the output of model (2.46)-(2.49), we consider x0 s,s= 1,...,295,to be all the instances in the Boston housing dataset that were given by the logistic regression model a probability of belonging to the positive class below 0.5, i.e., φ(f(x0 s)) <0.5,∀s= 1,2, . . . , S. (2.50) We calculate for them R= 3 counterfactual explanations in the one-for-many case, and one single explanation in the one-for-all case. We impose different values of ν,ν= 0.5,0.7,0.9. Figures 2.20-2.22 display the feature values of the clusters and the explanations. For the one-for-many allocation rule, for ν= 0.5, we note that, the first cluster is characterized by high values of the feature RM and AGE and lower values of features DIS, RAD and TAX, whereas cluster 2 is characterized by high values of RAD, TAX and PTRATIO. Cluster 3 can be characterized by lower values of most features except B. For the one-for-all allocation rule, the explanation that defines the benchmark for the positive class resembles the counterfactual associated with the first cluster in the case R= 3, but with higher values of the features RAD and TAX. Similar conclusions can be derived for other choices of ν, as shown in Figures 2.21 and 2.22. 2.2.3 A sparse one-for-many allocation rule model by column generation In the previous section for the one-for-many model, the cost function considered was just the Euclidean norm and the number of counterfactuals needed to explain the group of instances was a predefined parameter. In this subsection, we will: • Directly consider sparsity by limiting the number of features that can be changed in any instance and linking the features that are changed in any specific group, i.e., we implicitly detect subgroups of the population that share the same perturbed features in their counterfactual explanation. Doing so, we ease the use of the counterfactuals in the contexts in which interpretability is a requirement. 39 • Minimize the number of counterfactuals that can explain all the negative instances. Doing so, we significantly increase interpretability, because having fewer counterfactual explanations significantly contributes characterizing the set used to cover all the original instances. By limiting the number of features allowed to change in each subgroup, we are considering in our biobjective problem the cost function as a constraint. This allows us then to minimize the number of counterfactuals needed to explain the group. Because there is a limitation on how many features are allowed to be perturbed, we have to ensure that each of the instances assigned to their counterfactual moves no more than the limit imposed, thus, the minimum number of explanations needed to explain all the instances is not previously known. This is achieved by developing a sophisticated mixed-integer programming formulation based on a novel column generation framework. We compare this formulation with a simple adaptation of the mixed-integer programming formulation by Carrizosa et al. [2024a] and shown in Section 2.1. We adapt it to suit cases where the cost function solely integrates the global sparsity term, while also incorporating the one-for-many allocation rule, and we computationally show the significant advantages in terms of scalability, computational performance and quality of the solutions of the column generation framework. The framework incorporates as a black-box any classifier that can be formulated (after training) as a mixed-integer programming problem, for example, neural networks or logistic regression. Problem definition Recall that we have a probabilistic binary classification on a set X ⊂ {0,1}|J |. We will assume that all features xjare binary. In the case where there are numerical or categorical features, discretization followed by one-hot encoding will be considered. This assumption is made to be able to add sparsity in the one-for-many allocation rule model. We want to obtain the minimum number of counterfactual explanations that cover a set of Sinstances, while bounding the features that are changed in each subgroup. We want global explanations for the whole group and we want to detect the key features to flip the classifier’s decision for each of the groups. So, for the set of instances assigned to one counterfactual, the perturbed features will be the same, thus resulting in a much more interpretable rule. We will bound the number of features changed in each cluster by a given parameter Jmax, while our objective will be to minimize the number of total explanations used for the whole group. The above problem can be formulated through mixed-integer linear programming techniques. We will use the notation in Table 2.7. 40 the two solutions have the same value. Evaluation and analysis In Tables 2.10-2.15, the evaluation of the column generation framework is detailed. Specifically, for different values of the number of instances, |S|, and the maximum number of features allowed to change, Jmax, we report the smallest number of explanations obtained by the CG framework, y∗, the gap between the CG and the MILP solutions, and the respective runtimes. We can see in Table 2.10 for |S| = 10 that for small values of |S|, when dealing with a dataset with few features like COMPAS, the MILP is faster. However, for higher values of |S|, for the rest of the datasets and especially with NN as classifier, using the column generation framework is not only considerably faster, but also, when using the MILP we are not able to reach the same quality of the solution within the time limit. Notice that in those cases, the gap is always4negative, i.e., not only solving the MILP does not certify optimality but the number of explanations is higher than the value obtained with the CG methodology. For example, in Table 2.10 for |S| = 100 and Jmax = 2, a gap of -7 indicates that the solution of the MILP reached the time limit of 3600 CPU seconds with a set of 16 explanations, while the CG certifies that 9 are enough. Overall, the results in Tables 2.10-2.15 clearly indicates that (i) the CG is an effective and scalable approach to compute small, structured counterfactual sets, and (ii) a naive approach based on the extension of a MILP framework from the literature is significantly outperformed by the CG one. The effect of the score-based model used within the optimization framework (both CG and MILP) to classify the counterfactual explanations (see Section 2.2.3) is less clear from the experiments. • On the one side, the differences of accuracy between NNs and LR are so minor that it is hard to define a trend: ideally, a more accurate model could require less counterfactuals and this intuition seems to be confirmed by a couple of results in the tables. Namely, for the German credit dataset for |S| = 20 and Jmax = 10, we obtain 8 explanations for the LR whereas only 5 are needed by the NN that is slightly more accurate (see Table 2.9). Conversely, for the Students performance dataset for |S| = 20 and Jmax = 10, the slightly more accurate LR model needs 10 explanations versus the 11 of the NN. • On the other side, from an optimization standpoint, modeling NNs requires more constraints and variables, i.e., larger optimization models and Tables 2.10-2.15 clearly show that the CG approach does not suffer from scalability issues as much as the MILP does. 4The only exception is in Table 2.10 for |S| = 100 and Jmax = 3, where the best solution obtained by the MILP is one unit smaller than the one of the CG algorithm. 47 Thus, one can assume that for datasets where more advanced Machine Learning techniques (like NNs) are required to obtain satisfactory classification, CG would be the only option. Drawing explainability It is informative to picture the obtained counterfactuals, so as to make apparent that the explainability of the approach. For the COMPAS dataset, considering the NN and Jmax = 3, we show the counterfactual instances for different sizes of S, namely, in Figure 2.23, we show the three explanations for the |S| = 10. In the heatmap, each row is an explanation, and each column a feature. In green, we show the variables with a value of 1, and in white the ones that are 0. We visualize the same information for |S| = 20 in Figure 2.24. Finally, it is also worthwhile to see the most perturbed features. We can visualize for each of the datasets in the NN case, the number of total times one of the features is perturbed for all values of Jmax and |S| in Figures 2.25-2.27. 2.3 Chapter summary In this Chapter, models for each possible allocation rule have been outlined and illustrated. First, for the one-for-one allocation rule, an optimization model is proposed where minimal cost perturbations on the instances are sought. Perturbations costs are modeled taking into account both the magnitude of the change and also the number of predictor variables that need to be changed, either individually or globally. The model is expressed as a mathematical optimization problem, detailed for the case of a linear classifier (as e.g. in logistic regression or linear SVMs) and for ATM (as e.g. in random forests or XGBoost). Second, a model for the many-for-one allocation rule is outlined. We have detailed it for a specific type of diversity, where it is imposed that each of the counterfactuals has a different value for a specific feature. Third, two different models have been presented for the one-for-many allocation rule. Considering as a parameter the number of explanations needed and not taking into account sparsity, we have formulated the problem for linear models and ATM. On the other hand, we have addressed explainability altogether (i) by minimizing the number of explanations and at the same time (ii) by considering sparsity in the number of features allowed to change, with a novel column generation approach that has favorable properties in terms of scalability, computational performance and quality of the solutions. 48 (a) λglob = 0.1(b) λglob = 10 Figure 2.4: Counterfactual explanations for all the instances in the Boston housing dataset that comply with (2.41). The classifier considered is a logistic regression model. The explanations have been calculated solving model (ColCELR) with λind = 0,λglob = 0.1(left), 10 (right), and S∗=⌈0.95S⌉. The feature perturbations are displayed. 49 (a) Perturbations without (2.42) (b) Perturbations with (2.42) (c) Feature values without (2.42) (d) Feature values with (2.42) Figure 2.5: One-for-one counterfactual explanations for instances x0 5and x0 6in Table 2.3 from the Boston housing dataset where the classifier is the random forest. The explanations have been calculated solving model (ColCEATM) with λind = 0.01,λglob = 0 and ν= 0.5. Features perturbations are displayed in the two pictures above, without the Lipschitz continuity constraint (2.42) for τ= 10 and with this constraint, respectively, whereas, in the two pictures below, the corresponding features values are displayed. (a) Pareto frontier (b) Features perturbed Figure 2.6: On the left, the objective function of Problem (ColCELR) vs threshold value ν, when λind = 0,λglob = 0.1, for the logistic regression model. On the right, the features used in the counterfactual explanations. The explanations have been calculated for the 10 instances from the Boston housing dataset in Table 2.3. 50 (a) Logistic regression (b) Random forest Figure 2.7: Counterfactual explanations for instances x0 sin Table 2.4 from the COMPAS dataset for two classifiers, the logistic regression and the random forest. The explanations have been calculated solving model (CEsep) with Cs(x0 s,xs) = ∥x0 s−xs∥2 2+λind∥x0 s−xs∥0for λind = 0.02. The feature perturbations are displayed. (a) Logistic regression (b) Random forest Figure 2.8: Counterfactual explanations for instances x0 sin Table 2.4 from the COMPAS dataset for two classifiers, the logistic regression and the random forest. The explanations have been calculated solving models (ColCELR) and (ColCEATM) with cost function (2.1) with λind = 0,λglob = 0.2. The feature perturbations are displayed. 51 Figure 2.9: Features that need to be perturbed (in red) for all the instances in the COMPAS dataset that comply with (2.41) when solving Problem (ColCELRhard) for all the values of Fmax to describe the Pareto front. The classifier considered is a logistic regression model. Fmax = 1 yields an infeasible problem. 52 (a) λglob = 0.1(b) λglob = 10 Figure 2.10: Visualization of 100 of the counterfactual explanations for all the instances in the COMPAS dataset that comply with (2.41). The classifier considered is a logistic regression model. The explanations have been calculated solving model (ColCELR) with λind = 0, λglob = 0.1(left), 10 (right), and S∗=⌈0.95S⌉. The feature perturbations are displayed. 53 isMale isMarried AgeGroup EducationLevel MaxBillAmountOverLast6Months MaxPaymentAmountOverLast6Months MonthsWithZeroBalanceOverLast6Months MonthsWithLowSpendingOverLast6Months MonthsWithHighSpendingOverLast6Months MostRecentBillAmount MostRecentPaymentAmount TotalOverdueCounts TotalMonthsOverdue HasHistoryOfOverduePayments 1 1 3 4 0.1704 0.0107 0 0 6 0.2659 0.0357 1 22 1 1 0 2 3 0.0632 0.0119 0 0 5 0.0754 0.0396 1 32 1 1 0 2 2 0.1431 0.0169 0 0 6 0.2217 0.0201 1 12 1 0 1 2 4 0.2153 0.0095 0 0 6 0.3341 0.0318 1 12 1 1 0 2 4 0.1724 0.0072 0 0 6 0.2832 0.0240 2 8 1 1 1 3 4 0.0811 0.0087 0 0 5 0.1175 0.0292 2 8 1 1 0 2 4 0.1079 0.0089 0 0 5 0.1623 0.0298 1 12 1 1 1 3 2 0.1559 0.0086 0 0 6 0.2689 0.0220 1 12 1 1 0 2 4 0.1071 0.0086 0 0 6 0.1599 0.0285 1 8 1 1 0 2 3 0.1785 0.0078 0 0 6 0.2958 0.0259 1 12 1 Table 2.5: Feature values of x0= (x0 1,...,x0 10)from the Taiwanese credit dataset used for the models in Figures 2.11–2.12. 54 (a) Logistic regression (b) Random forest Figure 2.11: Counterfactual explanations for instances x0 sin Table 2.5 from the Taiwanese credit dataset for two classifiers, the logistic regression and the random forest. The explanations have been calculated solving model (CEsep) with Cs(x0 s,xs) = ∥x0 s−xs∥2 2+λind∥x0 s− xs∥0for λind = 0.02. The feature perturbations are displayed. (a) Logistic regression (b) Random forest Figure 2.12: Counterfactual explanations for instances x0 sin Table 2.5 from the Taiwanese credit dataset for two classifiers, the logistic regression and the random forest. The explanations have been calculated solving models (ColCELR) and (ColCEATM) with cost function (2.1) with λind = 0,λglob = 0.2. The feature perturbations are displayed. 55 Figure 2.13: Features that need to be perturbed (in red) for all the instances in the Taiwanese Credit dataset that comply with (2.41) when solving Problem (ColCELRhard) for all the values of Fmax to describe the Pareto front. The classifier considered is a logistic regression model. 56 (a) Cluster 1 and its counterfactual (b) Cluster 2 and its counterfactual (c) Cluster 3 and its counterfactual (d) All instances and its counterfactual Figure 2.21: One-for-many and one-for-all analysis for x0from the Boston housing dataset as in (2.50) for the logistic regression and R= 3 and R= 1 respectively. We consider Problem (2.46)-(2.49) with ν= 0.7. The plot represents the feature values of the counterfactual explanations (below) and their associated clusters (above). (a) Cluster 1 and its counterfactual (b) Cluster 2 and its counterfactual (c) Cluster 3 and its counterfactual (d) All instances and its counterfactual Figure 2.22: One-for-many and one-for-all analysis for x0from the Boston housing dataset as in (2.50) for the logistic regression and R= 3 and R= 1 respectively. We consider Problem (2.46)-(2.49) with ν= 0.9. The plot represents the feature values of the counterfactual explanations (below) and their associated clusters (above). 63 Table 2.12: Evaluation of the sparse one-for-many allocation rule model with the column generation framework for the Students performance dataset for the LR model for different values of |S| and Jmax. |S| Jmax y∗Gap CG time MILP time 10 10 7 0 29.726 587.94 20 10 10 -1 689.34 3600.00 50 10 33 -15 3600.00 3600.00 100 10 77 -23 3600.00 3600.00 10 15 4 0 34.20 3071.42 20 15 5 0 648.29 3600.00 50 15 12 -6 3600.00 3600.00 100 15 43 -57 3600.00 3600.00 10 20 2 0 4.48 5.85 20 20 2 0 87.41 40.16 50 20 3 0 3600.00 3600.00 100 20 4 0 3600.00 3600.00 Table 2.13: Evaluation of the sparse one-for-many allocation rule model with the column generation framework for the Students performance dataset for the NN for different values of |S| and Jmax. |S| Jmax y∗Gap CG time MILP time 10 10 7 0 31.337 3600.00 20 10 11 0 916.76 3600.00 50 10 33 -15 3290.40 3600.00 100 10 79 -21 3600.00 3600.00 10 15 4 0 35.900 3600.00 20 15 5 -1 930.42 3600.00 50 15 12 -36 3600.00 3600.00 100 15 45 -54 3600.00 3600.00 10 20 2 0 20.22 48.21 20 20 2 0 223.40 2887.83 50 20 3 -46 3600.00 3600.00 100 20 4 -92 3600.00 3600.00 Table 2.14: Evaluation of the sparse one-for-many allocation rule model with the column generation framework for the German credit dataset for the LR model for different values of |S| and Jmax. |S| Jmax y∗Gap CG time MILP time 10 5 9 0 5.31 548.38 20 5 19 0 21.62 3600.00 50 5 38 -12 1676.81 3600.00 100 5 90 -10 3600.00 3600.00 10 10 4 0 13.32 109.08 20 10 8 -1 432.72 2874.25 50 10 15 -35 3600.00 3600.00 100 10 58 -42 3600.00 3600.00 10 15 2 0 18.55 14.84 20 15 4 0 482.71 3600.00 50 15 6 -1 3600.00 3600.00 100 15 11 -80 3600.00 3600.00 Table 2.15: Evaluation of the sparse one-for-many allocation rule model with the column generation framework for the German credit dataset for the NN for different values of |S| and Jmax. |S| Jmax y∗Gap CG time MILP time 10 5 9 0 6.40 3600.00 20 5 19 0 67.52 3600.00 50 5 38 -12 2139.81 3600.00 100 5 90 -10 3600.00 3600.00 10 10 4 0 23.49 3600.00 20 10 5 -1 499.60 3600.00 50 10 15 -21 3600.00 3600.00 100 10 64 -36 3600.00 3600.00 10 15 2 0 10.22 79.35 20 15 4 0 459.90 3600.00 50 15 6 -2 3600.00 3600.00 100 15 11 -87 3600.00 3600.00 64 Figure 2.23: Counterfactual explanations for the COMPAS dataset for a NN model when solving a sparse one-for-many allocation model with |S| = 10 and Jmax = 3.Each row represents a counterfactual instance and each column a feature. The value of the counterfactual features are illustrated, in green if its value is 1, in white if it is 0. Figure 2.24: Counterfactual explanations for the COMPAS dataset for a NN model when solving a sparse one-for-many allocation model with |S| = 20 and Jmax = 3.Each row represents a counterfactual instance and each column a feature. The value of the counterfactual features are illustrated, in green if its value is 1, in white if it is 0. 65 Figure 2.25: Number of times each of the features is perturbed for the COMPAS dataset for a NN, when solving the sparse one-for-many allocation model for |S| ∈ {10,20,50,100}and Jmax = 2,3. 66 Figure 2.26: Number of times each of the features is perturbed for the Students performance dataset for a NN, when solving the sparse one-for-many allocation model for |S| ∈ {10,20,50,100}and Jmax = 10,15,20. 67 Figure 2.27: Number of times each of the features is perturbed for the German credit dataset for a NN, when solving the sparse one-for-many allocation model for |S| ∈ {10,20,50,100}and Jmax = 5,10,15. 68 Chapter 3 Counterfactual Analysis for Functional Data 70 In this Chapter, we specifically address the problem of interpretability when data are functions. This type of data arises in important domains such as econometrics, energy, or marketing [Jank and Shmueli, 2006; Sood et al., 2009; Sunar and Swaminathan, 2021]. There is rather extensive literature on the use of machine learning to analyse functional data, e.g., adapting Support Vector Machine models to functional data [Blanquero et al., 2019; Chaovalitwongse et al., 2008], using regression trees to detect critical intervals [Blanquero et al., 2023] or novel forms of interpretability when dealing with functional data [Carrizosa et al., 2022a; Martín- Barragán et al., 2014]. See also Aneiros et al. [2022]; Ramsay [2006] for an overview of methods for functional data analysis. While the literature on machine learning to analyse functional data is extensive, this is not the case for Counterfactual Analysis. Most of the work on counterfactual explanations focuses on tabular, image data or text data [Karimi et al., 2022a; Ramon et al., 2020; Tolkachev et al., 2022], and much less on functional data. In principle, one could apply the methods developed for tabular data also to functional data, just by considering that each feature is the measurement of the function at a time instant. However, doing so, fundamental information such as the autocorrelation structure along consecutive time instants would be lost. For this reason, some works on counterfactual explanations exploiting the functional nature of data have been suggested, e.g., Ates et al. [2021]; Delaney et al. [2021], but, as far as the authors know, none of them uses the structure and properties of the machine learning model. Moreover, when working with functional data, other types of distance measures may appear, such as the Dynamic Time Warping distance [Xing et al., 2010], which, with our methodology, we can consider. Regarding the allocation rules used, we will consider only the one-for-one allocation rule where no other linking constraints exists, thus, we will study the separable case in the previous Chapter in Section 2.1.1, but now considering that data are functional. An instance x∈ X ⊂ FJis defined as a vector of Jfunctional features. The counterfactual explanation x∈ X0⊂ X of the instance x0is a hypothetical instance generated by combining existing instances in the dataset, heareafter prototypes, so that the cost C(x,x0)of perturbing the features in x0to yield xis minimal. With this, we achieve certain interpretability goals. First, we can deal with multivariate functional data, i.e., our data are functions taking values in some RJ.Second, we are able to identify the instances from the dataset that generate the counterfactual for each instance, controlling how sparse the counterfactual explanation is, in terms of both the number of prototypes used to create the counterfactual xand the number of functional features changed to move from x0to x. Third, we can model the cost function Cby means of different distance measures, including popular distances in functional analysis such as the Dynamic Time Warping distance. We will show that, under mild assumptions, obtaining counterfactual explanations reduces to solving a Mixed Integer Convex Quadratic Model with linear constraints, which can be solved with standard optimization packages. The remainder of the Chapter is organized as follows. In Section 3.1, we model the pro- 72 the warping paths W, yielding min x,z,α,ξ,uλ0 J X j=1 ξj+λ2 Q X r=1 J X j=1 xj(q1r)−x0 j(q2r)2 s.t. (3.12) −(3.23),(3.24) −(3.25) (q11, q21)→(q12, q22)→. . . →(q1Q, q2Q)∈ W (CEFDTW) Notice how for a fixed warping path in W, constraints (3.12)-(3.25) are all linear, while we have both binary and continuous variables. Hence, if X0is again defined by linear constraints, since the objective function is quadratic, Problem (CEFDTW) is a Mixed Integer Convex Quadratic Model with linear constraints, that can be solved using standard optimization packages. For this reason we propose an alternating heuristic to solve Problem (CEFDTW), detailed in Algorithm 1. Algorithm 1: Algorithm to calculate counterfactual explanations with the DTW- based cost function (3.6). 1Initialisation: Let π(0) be the warping path (t1, t1)→(t2, t2). . . →(tn, tn); 2r←0; 3Solve Problem (CEFDTW) with π(r)as warping path, and obtain the counterfactual instance xand DTW distance δ(r); 4Find the optimal warping path π∗and its corresponding δ∗by solving: δ∗= min PQ r=1 PJ j=1 xj(q1r)−x0 j(q2r)2 s.t. (q11, q21)→(q12, q22)→. . . →(q1Q, q2Q)∈ W if δ∗=δ(r)then 5stop; 6else 7update π(r+ 1) = π∗; 8r←r+ 1; 9go to Step 3; 10 end 11 Output: counterfactual instance x. 3.3 Numerical illustration We will illustrate our methodology in two real-world datasets with functional data, one univariate and another multivariate, from the UCR archive [Dau et al., 2019]. For a given instance, we are able to identify the individuals of the dataset from which the corresponding counterfactual is made up and what their contribution is. Furthermore, we show the two different sparsities that we can obtain with our model, namely, the number of prototypes used for the counterfac- 79 tual and the number of functional features that change. The use of different distances, i.e., the Euclidean and the DTW distances, is also displayed. All the mathematical optimization problems have been implemented using Pyomo optimization modeling language [Hart et al., 2017, 2011] in Python 3.8. As solver, we have used Gurobi 9.0 [Gurobi Optimization, 2021]. A value of ϵ= 1e−6has been imposed in (3.12) and (3.13). The values of the big-Min (3.12) and (3.13) are node dependent, and they have been tightened following (2.25) and (2.26). For all the computational experiments, the classification model considered has been a Random Forest with T= 200 trees, maximum depth of 4 and wt=1 T, t = 1, . . . , T. and we have imposed ν= 0.5. Our experiments have been conducted on a PC, with an Intel R CoreTM i7-1065G7 CPU @ 1.30GHz 1.50 GHz processor and 16 gigabytes RAM. The operating system is 64 bits. The source code to reproduce all results is available at https://github.com/jasoneramirez/FunctionalCE. The first dataset, ItalyPowerDemand [Keogh et al., 2006], has one functional feature. There are 1096 instances and each instance is a time series of length 24, representing the power demand in Italy in six months. The binary classification task is to distinguish days from October to March (response value −1) from April to September (response value +1). The second dataset, NATOPS [Ghouaiel et al., 2017], has 24 functional time series of length 51 representing the X, Y, and Z coordinates of the left and right hand, wrist, thumb and elbows as captured by a Kinect 2 sensor. There are 260 instances and we chose two classes of the 6 that there are in the dataset. The binary classification task is thus to distinguish the gesture “All Clear" (response value −1) from “Not Clear" (response value +1). ItalyPowerDemand We present the counterfactual for an instance x0of the dataset ItalyPowerDemand in Figure 3.2. In each case, we represent the original curve, the prototypes, and the final counterfactual. The first cost model analysed is the squared Euclidean model (3.5) with λ0= 0 (since we have only one feature, λ0>0is meaningless). Different values of Bmax have been used. The smaller the value of Bmax is, the more sparse the counterfactual is in terms of prototypes, while the larger the value of Bmax is, the higher the freedom to use prototypes and therefore the smaller the distances obtained. In Figure 3.3a we plot the relation between the number of prototypes and the distance. It is illustrated how using more than one prototype may be beneficial, but using more than 4 prototypes gives us less sparsity without smaller distances. To show the flexibility of our model, the same experiments have been carried out but changing the cost based on DTW distances (3.6), again with λ0= 0. The counterfactual solutions have been calculated with the heuristic procedure described in Algorithm 1. The results are depicted in Figure 3.4. As before, one can see how the objective function decreases as the number of prototypes Bmax increases. However, in this case, it is sufficient to use 2 prototypes, as 3 or more will not improve much the objective function, see Figure 3.3b. 80 (a) Bmax = 1 (b) Bmax = 2 (c) Bmax = 3 (d) Bmax = 4 Figure 3.2: Counterfactual explanations for x0of the ItalyPowerDemand dataset. Different values of Bmax, i.e., the number of prototypes used for the convex combination, have been imposed. The cost function is (3.5) with λ0= 0, λ2= 1, and ν= 0.5. NATOPS We present now the counterfactual for an instance x0of the multivariate dataset NATOPS. The cost function used has been of the form is the squared Euclidean model (3.5) with λ0= 1, λ2= 0.005.As we are interested in detecting the critical functional features to flip the classifier’s decision, we give more weight in the cost function to the ℓ0as an illustration. In Figure 3.5 the counterfactual instance xfor x0for Bmax = 1 is shown. As the cost function Ccontains as its first term the ℓ0norm, we obtain a sparse solution in the sense of the features we need to change to move from x0to x. Indeed, to change its class, only three functional features have to be modified. In Figure 3.6 the changed features are presented. As in the univariate case, we can impose different values of Bmax. In Figure 3.7 we show the counterfactual explanation for Bmax = 2 and for the same cost function. Note how giving the flexibility to use more than one prototype, results in only having to change two features, 81 (a) Euclidean distance (b) DTW distance Figure 3.3: Distance obtained vs number of prototypes used for the convex combination for the counterfactual explanation xfor x0of the ItalyPowerDemand dataset. see Figure 3.8. 3.4 Chapter summary In this Chapter, we have proposed a novel approach to build counterfactual explanations when dealing with multivariate functional data in classification problems by means of mathematical optimization. With our method, we ensure feasible and sparse explanations, controlling not only the number of prototypes of the dataset used to create the counterfactuals, but also the number of features that need to be changed. Our model is also flexible enough to be used with different distance measures, e.g., the Euclidean distance or the DTW distance. Moreover, our methodology is applicable to score-based classifiers, including additive tree models, such as random forest or XGBoost models, as well as linear models, such as logistic regression and linear support vector machines. We have illustrated our methodology on various real-world datasets. 82 (a) Bmax = 1 (b) Bmax = 2 (c) Bmax = 3 (d) Bmax = 4 Figure 3.4: Counterfactual explanations for x0of the ItalyPowerDemand dataset. Different values of Bmax, i.e., the number of prototypes used for the convex combination, have been imposed. The cost function is (3.6) with λ0= 0, λ2= 1, and ν= 0.5. 83 Figure 3.5: Counterfactual explanations for x0of the NATOPS dataset for a Random Forest. Bmax = 1 prototype has been imposed to create the counterfactual. The cost function is (3.5) with λ0= 1, λ2= 0.005, and ν= 0.5. 84 (a) Feature 7 (b) Feature 11 (c) Feature 13 Figure 3.6: Changed features in the counterfactual explanation for x0of the NATOPS dataset for a Random Forest with Bmax = 1. The cost function is (3.5) with λ0= 1, λ2= 0.005, and ν= 0.5. 85 Figure 3.7: Counterfactual explanations for x0of the NATOPS datase for a Random Forest. Bmax = 2 prototypes has been imposed to create the counterfactual. The cost function is (3.5) with λ0= 1, λ2= 0.005, and ν= 0.5. 86 (a) Feature 11 (b) Feature 22 Figure 3.8: Changed features in the counterfactual explanation for x0of the NATOPS dataset for a Random Forest with Bmax = 2. The cost function is (3.5) with λ0= 1, λ2= 0.005, and ν= 0.5. 87 and is expressed in percentages of the existing inputs and outputs. The resulting best practice alternative is illustrated in the “Benchmark” column. We see that the DMU in this example expresses an interest in reducing Supervision and Credit personnel but simultaneously seeks to increase the number of personal loan accounts. Figure 4.2: Directional search in Interactive Benchmarking software. Real-world dataset of bank branches in Section 4.5. Applications of interactive benchmarking have typically been in settings where the DMU in a trial-and-error like process seeks alternative production plans. Such processes can certainly be useful in attempts to learn about and gain trust in the modeling, to guide decision making and to find the performance enhancing changes that a DMU may find relatively easy to implement. From the point of view of Multiple Criteria Decision Making (MCDM) we can think of such processes as based on progressive articulation of preferences and alternatives, cf. e.g. the taxonomy of MCDM methods suggested in Rostami et al. [2017]. It is clear from this small example, however, that the use of an interactive process guided solely by the DMU may not always be the best approach. If there are more than a few inputs and outputs, the process can become difficult to steer towards some underlying optimal compromise between the many possible changes in inputs and outputs. In such cases, the so-called prior articulation of preferences may be more useful. If the DMU can express its preferences for different changes, e.g., as a cost of change function C((x,y),(x∗,y∗)) giving the cost of moving from the present production plan (x,y)to any new production plan (x∗,y∗), then a systematic search for the optimal change is possible. Our approach is based on this idea. We consider a class of cost functions and show how to find optimal changes in inputs and outputs using bilevel optimization. In this sense, it corresponds to endogenizing the directional choice so as to make the necessary changes in inputs and outputs as small as possible. Of course, by varying the parameters of the cost function, one can also generate a reasonably representative set of alternative production plans that the DMU can then choose from. This would correspond to the idea of a prior articulation of alternatives approach in the MCDM taxonomy. The idea of introducing a cost of change function to guide the search for alternative production plans is closely related to the idea of targets in DEA. At a general level, a target is 95 here understood as an alternative production plan that a DMU should move to.1There has been a series of interesting DEA papers on the determination of targets using the principle of least action, see for example Aparicio et al. [2014] and the references in here. In Aparicio et al. [2014], the authors explicitly introduce the principle of least action referring to the idea in physics that nature always finds the most efficient course of action. The general argument underlying these approaches is that an inefficient firm should achieve technical efficiency with a minimum amount of effort. Different solutions have been proposed using different distance measures or what we call the cost of change. In many papers, this corresponds to minimizing the distance to the efficient frontier in contrast to the traditional efficiency measurement problem, where we are looking for the largest possible savings or the largest possible expansions of the services provided. A good example is Aparicio and Pastor [2013]. In this sense, our idea of finding close counterfactuals fits nicely into the DEA literature. The choice of targets has also been discussed in connection with the slack problem in radial DEA measures. A Farrell projection may not lead to a Pareto efficient point and in a second stage, it is therefore common to discuss close alternatives that are fully Pareto efficient. Again, different solutions – using for example constraints or the determination of all full facets – have been proposed, cf. Aparicio et al. [2017]. Identifying all facets and measuring distance to these like in Silva Portela et al. [2003] is theoretically attractive but computationally cumbersome in most applications. Our approach can be seen as a generalization of the literature on targets. In particular, if E∗= 1 and the considered cost function coincides with a norm or with a typical technical efficiency measure, then the previous DEA target approaches are particular cases of the general approach introduced in this Chapter. We formulate the target setting problem for general cost-of-change functions using a bilevel program2, and we reformulate the constraints to get tractable mathematical optimization problems. Using combinations of ℓ0,ℓ1and ℓ2norms, as we do in the illustrations, the resulting problems are Mixed Integer Convex Quadratic Problems with linear constraints. It is worthwhile to note also, that we do not necessarily require the target to be Pareto efficient, allowing for the possibility that a DMU may not seek to become fully efficient but, for example, just 90% Farrell efficient which also implies that targets may be on non-full facets. In the following, we combine the ideas of DEA, least action targets, and counterfactual explanations. We formulate and solve bilevel optimization models to determine “close” alternative production plans or counterfactual explanations in DEA models that lead to desired 1In Aparicio et al. [2007], the authors distinguish between setting targets and benchmarking in the sense that targets are the coordinates of a projection point, which is not necessarily an observed DMU, whereas benchmarks are real observed DMUs. This distinction can certainly be relevant is several contexts, but it is not how we use targets here. We use benchmarking as the general term for relative performance comparison, and target to designate an alternative production plan that a DMU should choose to improve performance in the easiest possible way. 2The idea of using a bilevel linear programming approach has also appeared in the DEA literature. It should be noted in particular that Aparicio et al. [2017] proposed to resort to a bilevel linear programming model when strictly efficient targets are to be identified using the Russel output measure to capture cost-of-change. 96 relative performance levels and also take into account the strategic preferences of the entity. 4.3 The Setting We consider K+ 1 DMUs (indexed by k), using Iinputs, xk= (xk 1, . . . , xk I)⊤∈RI +, to produce Ooutputs, yk= (yk 1, . . . , yk O)⊤∈RO +. Hereafter, we will write (xk,yk)to refer to production plan of DMU k,k= 0,1, . . . , K. Let Tbe the technology set, with T={(x,y)∈RI +×RO +|xcan produce y}. We will initially estimate Tby the classical DEA model. It determines the empirical reference technology T∗as the smallest subset of RI +×RO +that contains the actual K+ 1 observations, and satisfies the classical DEA regularities of convexity, free-disposability in inputs and outputs, and Constant Returns to Scale (CRS). It is easy to see that the estimated technology can be described as: T∗(CRS) = {(x,y)∈RI +×RO +| ∃λ∈RK+1 +:x≥ K X k=0 λkxk,y≤ K X k=0 λkyk}. To measure the efficiency of a firm, we will initially use the so-called Farrell input-oriented efficiency. It measures the efficiency of a DMU, say DMU 0, as the largest proportional reduction E0of all its inputs x0that allows the production of its present outputs y0in the technology T∗. Hence, it is equal to the optimal solution value of the following LP formulation min E,λ0,...,λKE(DEA) s.t. Ex0≥ K X k=0 λkxk y0≤ K X k=0 λkyk 0≤E≤1 λ∈RK+1 +. This DEA model has K+ 2 decision variables, Ilinear input constraints and Olinear output constraints. Hereafter, we will refer to the optimal objective value of (DEA), say E0, as the efficiency of DMU 0. In the following, and assuming that firm 0 with production plan (x0,y0)is not fully efficient, E0<1, we will show how to calculate a counterfactual explanation with a desired efficiency level E∗> E0, i.e., the minimum cost changes needed in the inputs of the firm, 97 x0, in order to obtain an efficiency E∗. Given a cost function C(x0,ˆ x)that measures the cost of moving from the present inputs x0to the new counterfactual inputs ˆ x, and a set X(x0) defining the feasible space for ˆ x, the counterfactual explanation for x0is found by solving the following optimization problem: min ˆ xC(x0,ˆ x) s.t. ˆ x∈ X(x0) (ˆ x,y0)has at least an efficiency of E∗. With respect to C(x0,ˆ x), different norms can be used to measure the difficulty of changing the inputs. A DMU may, for example, be interested to minimize the sum of the squared deviations between the present and the counterfactual inputs. We model this using the squared Euclidean norm ℓ2 2. Likewise, there may be an interest in minimizing the absolute value of the deviations, which we can proxy using the ℓ1norm, or the number of inputs changed, which we can capture with the ℓ0norm. When it comes to X(x0), this would include the non-negativity of ˆ x, as well as domain knowledge specific constraints. With this approach, we detect the most important inputs in terms of the impact they have on the DMU’s efficiency, and with enough flexibility to consider different costs of changing depending on the DMU’s characteristics. In the next section, we will show that finding counterfactual explanations involves solving a bilevel optimization problem of minimizing the changes in inputs and solving the above DEA problem at the same time. In Section 4.6, we will also discuss how the counterfactual analysis approach can be extended to other technologies and to other efficiency measures like the output-oriented Farrell efficiency and other DEA technologies. Before turning to the details of the bilevel optimization problem, it is useful to illustrate the idea of counterfactual explanations in benchmarking using a small numerical example. Suppose we have four firms with the inputs, outputs, and Farrell input efficiencies as in Table 4.1. The efficiency has been calculated solving the classical DEA model with CRS, namely (DEA). In this example, firms 1 and 2 are fully efficient, whereas firms 3 and 4 are not. Firm x1x2y E 1 0.50 1 1 1 2 1.50 0.50 1 1 3 1.75 1.25 1 0.59 4 2.50 1.25 1 0.50 Table 4.1: Inputs, outputs and corresponding Farrell input-efficiency of 4 different firms when solving the classical DEA model with CRS. First, we want to know the changes needed in x3for firm 3 to have a new efficiency E∗of at least 80%. Since we only have two inputs, we can illustrate this graphically as in Figure 4.3a. The results are shown in Table 4.2 for different cost functions. It can be seen that we, in all 98 cases, get exactly 80% efficiency with the new inputs. We see from column ℓ2 2that the Farrell solution, i.e., solving (DIR) with dx=xand dy=0, is further away from the original inputs than the counterfactual solution based on the Euclidean norm. To the extent that difficulties of change is captured by the ℓ2 2norm, we can conclude that the Farrell solution is not ideal. Moreover, in the Farrell solution one must by definition change both inputs, see column ℓ0. Using a cost function combining the ℓ0norm and the squared Euclidean norm, denoted by ℓ0+ℓ2, one penalizes the number of inputs changed. With this we detect the one input that should change in order to obtain a higher efficiency, namely the second input. In contexts like negotiations with various input suppliers, it is often more practical to focus negotiations on just one or a select few inputs, instead of dealing with all inputs at the same time. Cost function ˆx1ˆx2y E ℓ2 2ℓ0 Farrell 1.29 0.92 1 0.8 0.32 2 ℓ0+ℓ21.75 0.69 1 0.8 0.31 1 ℓ21.53 0.80 1 0.8 0.25 2 Table 4.2: Counterfactual explanations for firm 3 in Table 4.1 imposing E∗= 0.8and different cost functions. Let us now focus on firm 4 and again find a counterfactual instance with at least 80% efficiency. The results are shown in Table 4.3 and Figure 4.3b. Notice how in the Farrell case one obtains again the farthest solution and also the least sparse from the three of them. As for the counterfactual explanations with our methodology, the inputs nearest to the original DMU that give us the desired efficiency are in a non full-facet of the efficiency frontier. By using the Farrell to measure the desired efficiency level, we only need to change one input, namely, the second input, and can have “slack" in the first input. We here deviate from Aparicio et al. [2017], in which the authors look for targets on the strongly efficient frontier, i.e., without slack. Cost function ˆx1ˆx2yEℓ2 2ℓ0 Farrell 1.56 0.78 1 0.8 1.10 2 ℓ0+ℓ22.50 0.63 1 0.8 0.39 1 ℓ22.50 0.63 1 0.8 0.39 1 Table 4.3: Counterfactual explanations for firm 4 in Table 4.1 imposing E∗= 0.8and different cost functions. In Figure 4.3, the space where we search for the counterfactual explanation is shaded. Although in these illustrations the frontier is explicitly given, in general, the frontier points are convex combinations of the observed DMUs, and obtaining the isoquants is not easy. Therefore, when finding counterfactual explanations a bilevel optimization is needed to search for “close” inputs in the complement of a convex set. 99 (a) Explanations for firm 3 (b) Explanations for firm 4 Figure 4.3: Counterfactual explanations for firms 3 and 4 in Tables 4.2 and 4.3 respectively imposing E∗= 0.8and different cost functions. 4.4 Bilevel optimization for counterfactual analysis in DEA Suppose DMU 0 is not fully efficient, i.e., the optimal objective value of Problem (DEA) is E0<1. In this section, we formulate the counterfactual explanation problem in DEA, i.e., the problem that calculates the minimum cost changes in the inputs x0that make DMU 0 have a higher efficiency. Let ˆ xbe the new inputs of DMU 0that would make it at least E∗efficient, with E∗> E0. With this, we have defined the counterfactual instance as the one obtained changing the inputs, but in the same sense, we could define it by changing the outputs. This alternative output-based problem will be studied in Section 4.6. Since the values of the inputs are to be changed, the efficiency of the new production plan (ˆ x,y0)has to be calculated using Problem (DEA). The counterfactual explanation problem in DEA reads as follows: min ˆxC(x0,ˆx)(4.1) s.t. ˆx∈RI +(4.2) E≥E∗(4.3) E∈arg min ¯ E,λ0,...,λK {¯ E:¯ Eˆ x≥ K X k=0 λkxk,y0≤ K X k=0 λkyk,¯ E≥0,λ∈RK+1 +}, (4.4) where in the upper level problem in (4.1) we minimize the cost of changing the inputs for firm 0,x0, to ˆ x, ensuring nonnegativity of the inputs, as in constraint (4.2), and that the efficiency is at least E∗, as in constraint (4.3). The lower level problem in (4.4) ensures that the efficiency of (ˆ x,y0)is correctly calculated. Therefore, as opposed to Counterfactual Analysis in 100 interpretable machine learning and what we have studied in previous Chapters, here we are confronted with a bilevel optimization problem. Notice also that to calculate the efficiency in the lower level problem in (4.4), the technology is already fixed, and the new DMU (ˆ x,y0) does not take part in its calculation. In what follows, we reformulate the bilevel optimization problem (4.1)-(4.4) as a singlelevel model, by exploiting the optimality conditions for the lower-level problem. This can be done for convex lower-level problems that satisfy Slater’s conditions, e.g., if our lower-level problem was linear. In our case, however, not all the constraints are linear, since in (4.4) we have the product of decision variables ¯ Eˆ x. To be able to handle this, we define new decision variables, namely, F=1 Eand βk=λk E, for k= 0, . . . , K. Thus, (4.1)-(4.4) is equivalent to: min ˆx,F C(x0,ˆx) s.t. ˆx∈RI + F≤F∗ F∈arg max ¯ F,β {¯ F:ˆ x≥ K X k=0 βkxk,¯ Fy0≤ K X k=0 βkyk,¯ F≥0,β∈RK+1 +}.(4.5) This equivalent bilevel opimization problem can now be reformulated as a single-level model. The new lower-level problem in (4.5) can be seen as the ˆ x-parametrized problem: max F,βF(4.6) s.t. ˆ x≥ K X k=0 βkxk(4.7) Fy0≤ K X k=0 βkyk(4.8) F≥0(4.9) β≥0.(4.10) The Karush-Kuhn-Tucker (KKT) conditions, which include primal and dual feasibility, stationarity and complementarity conditions, are necessary and sufficient to characterize an optimal solution. Thus, we can replace problem (4.5) by its KKT conditions. Primal feasibility is given by (4.7)-(4.10). Dual feasibility is given by: γI,γO, δ, µ≥0,(4.11) for the dual variables associated with constraints (4.7)-(4.10), where γI∈RI +,γO∈RO +, 101 δ∈R+,µ∈RK+1 +. The stationarity conditions are as follows: γ⊤ Oy0−δ= 1 (4.12) γ⊤ Ixk−γ⊤ Oyk−µk= 0 k= 0, . . . , K. (4.13) Lastly, we need the complementarity conditions for all constraints (4.7)-(4.10). For constraint (4.7), we have: γi I= 0 or ˆxi− K X k=0 βkxk i= 0 i= 1, . . . , I. (4.14) In order to model this disjunction, we will introduce binary variables ui∈ {0,1},i= 1, . . . , I, and the following constraints using the big-M method: γi I≤MIui,ˆxi− K X k=0 βkxk i≤MI(1 −ui), i = 1, . . . , I, (4.15) where MIis a sufficiently large constant. The same can be done for the complementarity condition for constraint (4.8), introducing binary variables vo∈ {0,1},o= 1, . . . , O, big-M constant MO, and constraints: γo O≤MOvo,−Fy0 o+ K X k=0 βkyk o≤MO(1 −vo), o = 1, . . . , O. (4.16) The complementarity condition for constraint (4.10) would be the disjunction βk= 0 or µk= 0. Using the stationarity condition (4.13) and again the big-M method with binary variables wk∈ {0,1},k= 0, . . . , K, and big-M constant Mf, one obtains the constraints: βk≤Mfwk,γ⊤ Ixk−γ⊤ Oyk≤Mf(1 −wk), k = 0, . . . , K. (4.17) Finally, for constraint (4.9) the complementarity condition yields F= 0 or δ= 0. Remember that F= 1/E,0≤E≤1, thus Fcannot be zero by definition and we must impose δ= 0. Using stationarity condition (4.12), this yields: γ⊤ Oy0= 1.(4.18) We now reflect on the meaning of these constraints. Notice that constraints (4.15) and (4.16) model the slacks of the inputs and outputs respectively, while constraint (4.17) models the firms that define the frontier, i.e., the firms with which DMU 0 is to be compared. If binary variable ui= 1, then there is no slack in input i, i.e., ˆxi=PK k=1 βkxk i, whereas if ui= 0 that means there is. The same happens with binary variable vo, namely, it indicates whether there is a slack in output o. On the other hand, when wk= 1, then the equality of the dual constraint will hold γ⊤ Ixk=γ⊤ Oyk, i.e., firm kis fully efficient and it is used to define the efficiency of 102 the counterfactual instance. If wk= 0 then βk= 0, and firm kis not being used to define the efficiency of the counterfactual instance. Let us go back to the example in the previous section with four firms with 2 inputs and 1 output and several choices of cost function Cof changing the inputs. When C=ℓ2 2, we can see that firm 3 is compared against firms 1 and 2, while firm 4 is compared against firm 2 only. Notice that µkis only present in (4.13), thus it is free. In addition, we know that δ= 0. Therefore, we can transform the stationarity conditions (4.12) and (4.13) to γ⊤ Oy0= 1 (4.19) γ⊤ Ixk−γ⊤ Oyk≥0k= 1, . . . , K (4.20) γI,γO≥0,(4.21) that are exactly the constraints in the dual DEA model for the Farrell output efficiency. The new reformulation of the counterfactual explanation problem in DEA is as follows: min ˆx,F,β,γI,γO,u,v,wC(x0,ˆx) s.t. F≤F∗ ˆx∈RI + u,v,w∈ {0,1} (4.7) −(4.10) primal (4.19) −(4.21) dual (4.15) −(4.16) slacks (4.17) frontier. So far, we have not been very specific about the objective function C(x0,ˆ x). Different functional forms can be introduced, and this may require the introduction of further variables to implement these. In Section 4.3, we approximated the firm’s cost-of-change using combinations of the ℓ0 norm, the ℓ1norm, and the squared ℓ2norm. As a starting point, we therefore propose the following objective function: C(x0,ˆ x) = ν0∥x0−ˆ x∥0+ν1∥x0−ˆ x∥1+ν2∥x0−ˆ x∥2 2,(4.22) where ν0, ν1, ν2≥0. Taking into account that there may be specific product input prices and output prices or that inputs may have varying degrees of difficulty to be changed, one can consider giving different weights to the deviations in each of the inputs. In order to have a smooth expression of objective function (4.22), additional decision vari- 103 ables and constraints have to be added to the counterfactual explanation problem in DEA. To linearize the ℓ0norm, binary decision variables ξiare introduced. For input i,ξi= 1 models x0 i= ˆxi,i= 1, . . . , I. Using the big-M method the following constraints are added to our formulation: −Mzeroξi≤x0 i−ˆxi≤Mzeroξi, i = 1, . . . , I (4.23) ξi∈ {0,1}, i = 1, . . . , I, (4.24) where Mzero is a sufficiently large constant. For the ℓ1norm we introduce continuous decision variables ηi≥0,i= 1, . . . , I, to measure the absolute values of the deviations, ηi=x0 i−ˆxi, which is naturally implemented by the following constraints: ηi≥x0 i−ˆxi, i = 1, . . . , I (4.25) −ηi≤x0 i−ˆxi, i = 1, . . . , I (4.26) ηi≥0, i = 1, . . . , I. (4.27) Thus, the counterfactual explanation problem in DEA with cost function Cin (4.22), hereafter (CEDEA), reads as follows: min ˆx,F,β,γI,γO,u,v,w,η,ξν0 I X i=1 ξi+ν1 I X i=1 ηi+ν2 I X i=1 η2 i(CEDEA) s.t. F≤F∗ ˆx∈RI + u,v,w∈ {0,1} (4.7) −(4.10),(4.15) −(4.17),(4.19) −(4.21), (4.23) −(4.24),(4.25) −(4.27). Notice that in Problem (CEDEA) we assumed X(x0) = RI +as the feasible space for ˆx. Other relevant constraints for the counterfactual inputs could easily be added, e.g., bounds or relative bounds on the inputs, or inputs that cannot be changed in the short run, say capital expenses, or that represent environmental conditions beyond the control of the DMU. In the case where only the ℓ0and ℓ1norms are considered, i.e., ν2= 0, the objective function as well as the constraints are linear, while we have both binary and continuous decision variables. Therefore, Problem (CEDEA) can be solved using an Mixed Integer Linear Programming (MILP) solver. Otherwise, when ν2= 0, Problem (CEDEA) is a Mixed Integer Convex Quadratic model with linear constraints, which can be solved with standard optimization packages. When all three norms are used, Problem (CEDEA) has 3I+K+2+Ocontinuous variables and 2I+O+K+ 1 binary decision variables. It has 7I+ 3O+ 3K+ 5 constraints, 104 (a) C=ℓ0+ (ℓ2)(b) C=ℓ0+ℓ2(c) C=ℓ2 Figure 4.5: The inputs that change (in red) when solving Problem (CEDEA) with desired efficiency E∗= 1 for all the inefficient firms in the Canadian bank dataset. Figure 4.6: Counterfactual Explanations for DMU 238 from the Canadian bank dataset when solving Problem (CEDEA) with desired efficiency E∗= 0.8. 111 Cost function Teller Typing Acc. and Ledgers Supervision Credit Mean ℓ0(x−ˆx) ℓ0+ (ℓ2)0.01 0.30 0.18 0.08 0.44 1.00 ℓ0+ℓ20.03 0.32 0.19 0.09 0.47 1.11 ℓ20.13 0.43 0.51 0.18 0.63 1.88 Table 4.8: Average results on how much the inputs (personnels) change when solving Problem (CEDEA) with desired efficiency E∗= 0.8for all the inefficient firms in the Canadian bank dataset. Cost function Teller Typing Acc. and Ledgers Supervision Credit Mean ℓ2(r) ℓ0+ (ℓ2)0.28 0.56 0.28 0.41 0.27 0,3708 ℓ0+ℓ20.12 0.52 0.25 0.39 0.26 0.3707 ℓ20.05 0.41 0.11 0.25 0.21 0.3702 Table 4.9: Average results on how much the inputs (personnels) change when solving Problem (CEDEA) with desired efficiency E∗= 0.8for all the inefficient firms in the Canadian bank dataset. can see now an increasing number of rows with no markings compared to Figure 4.5, belonging to the firms that had already an efficiency of 0.8. 4.6 Extensions Here, we extend the analysis in Section 4.4 by investigating alternative returns to scale and by investigating changes in the outputs rather than the inputs. In both extensions, we will consider a combination of the ℓ0,ℓ1and ℓ2norms as in objective function (4.22). 4.6.1 Changing the returns to scale In Section 4.4, we have considered the DEA model with constant return to scale (CRS), where the only requirement on the values of λis that they are positive, i.e., λ∈RK+1 +, but we could consider other technologies. In that case, to be able to transform our bilevel optimization problem to a single-level one, we should take into account that for each new constraint derived from the conditions on λ, a new dual variable has to be introduced. We will consider the varying return to scale (VRS) model as it is one of the most preferred one by firms [Bogetoft, 2012], but extensions to other models are analogous. Consider the input case. With the same transformation as before, we have: 112 (a) C=ℓ0+ (ℓ2)(b) C=ℓ0+ℓ2(c) C=ℓ2 Figure 4.7: The inputs that change (in red) when solving Problem (CEDEA) with desired efficiency E∗= 0.8for all the inefficient firms in the Canadian bank dataset. min ˆx,F C(x0,ˆx) s.t. ˆx∈RI + F≤F∗ F∈arg min ¯ F,λ0,...,λK {¯ F:ˆx≥ K X k=0 βkxk,ˆ Fy0≤ K X k=0 βkyk,¯ F≥0,β∈RK+1 +, K X k=0 βk=F}. Notice that the only difference is that we have a new constraint associated with the technology, namely, PK k=0 βk=F. Let κ≥0be the new dual variable associated with this constraint. Then, the following changes are made in constraints (4.19) and (4.20): γ⊤ Oy0+κ= 1 (4.28) γ⊤ Ixk−γ⊤ Oyk−κ≥0k= 0, . . . , K. (4.29) The single-level formulation for the counterfactual problem for VRS DEA is as follows: min ˆx,F,β,γI,γO,u,v,w,κ,η,ξν0 I X i=1 ξi+ν1 I X i=1 ηi+ν2 I X i=1 (x0 i−ˆxi)2(CEVDEA) 113 s.t. F≤F∗ K X k=0 βk=F ˆx∈RI + κ≥0 u,v,w∈ {0,1} (4.7) −(4.10),(4.15) −(4.17) (4.21),(4.23) −(4.29). 4.6.2 Changing the outputs We have calculated the counterfactual instance of a firm as the minimum cost changes in the inputs in order to have a better efficiency. In the same vein, we could consider instead changes in the outputs, leaving the same inputs. Again, suppose firm 0 is not fully efficient, E0<1. Now, we are interested in calculating the minimum changes in the outputs y0that make it to have a higher efficiency E∗> E0. Let ˆ ybe the new outputs of firm 0 that make it to be at least E∗efficient. We have, then, the following bilevel optimization problem: min ˆy,E C(y0,ˆy) s.t. ˆy∈RO + E≥E∗ E∈arg min ¯ E,λ0,...,λK {¯ E:¯ Ex0≥ K X k=0 λkxk,ˆ y≤ K X k=0 λkyk,¯ E≥0,λ∈RK+1 +}. Following similar steps as in previous section, the single-level formulation for the counterfactual problem in DEA in the output case is as follows: min ˆy,E,λ,γI,γO,u,v,w,η,ξν0 O X o=1 ξo+ν1 O X o=1 ηo+ν2 O X o=1 (y0 o−ˆyo)2(CEODEA) s.t. ˆy∈RO + E≥E∗ Ex0≥ K X k=0 λkxk ˆy≤ K X k=0 λkyk 114 γ⊤ Ix0= 1 γ⊤ Oyk−γ⊤ Ixk≤0k= 0, . . . , K γi I≤MIuii= 1, . . . , I Ex0 i− K X k=0 λkxk i≤MI(1 −ui)i= 1, . . . , I γo O≤MOvoo= 1, . . . , O −ˆyo+ K X k=0 λkyk o≤MO(1 −vo)o= 1, . . . , O λk≤Mfwkk= 0, . . . , K γ⊤ Oyk−γ⊤ Ixk≤Mf(1 −wk)k= 0, . . . , K −Mzeroξo≤y0 o−ˆyoo= 1, . . . , O y0 o−ˆyo≤Mzeroξoo= 1, . . . , O ηo≥y0 o−ˆyoo= 1, . . . , O −ηo≤y0 o−ˆyoo= 1, . . . , O E, λ,γI,γO,η≥0 u,v,w,ξ∈ {0,1}. As in the input model, depending on the cost function, we either obtain an MILP model or a Mixed Integer Convex Quadratic model with linear constraints. This model could be formulated analogously for the VRS case. 4.7 Chapter summary In this Chapter, we have proposed a collection of optimization models to setting targets and finding counterfactual explanations in DEA models, i.e., the least costly changes in the inputs or outputs of a firm that leads to a pre-specified (higher) efficiency level. With our methodology, we are able to include different ways to measure the proximity between a firm and its counterfactual, namely, using the ℓ0,ℓ1, and ℓ2norms or a combination of them. Calculating counterfactual explanations involves finding “close" alternatives in the complement of a convex set. We have reformulated this bilevel optimization problem as either an MILP or a Mixed Integer Convex Quadratic Problem with linear constraints. In our numerical section, we can see that for our banking application, we are able to solve this model to optimality. DEA models can capture very complex relationships between multiple inputs and outputs. This allows more substantial evaluations and also offers a framework that can support many operational, tactical and strategic planning efforts. However, there is also a risk that such a model is seen as a pure black-box which in turn can lead to mistrust and some degree of model 115 or algorithm aversion. By looking at counterfactuals, a firm can get a better understanding of the production space and is more likely to trust in the modeling. Counterfactuals in DEA can also help a firm choose which changes to implement. It is not always enough to simply think of a strategy and which factors can easily be changed, say a direction in input space. It is also important how the technology looks like and therefore how large such changes need to be to get a desired improvement in efficiency. In this way, the analysis of close counterfactuals can help endogenize the choice of both desirable and effective directions to move in. By varying the parameters of the cost function, the firm can even get a menu of counterfactuals, from which it can choose, having thus more flexibility and leading the evaluated firm to gain more trust in the underlying model. Note also that by calculating the counterfactual explanations for all firms involved, as we did in our banking application, one can determine which combinations of inputs and outputs that most commonly shall be changed to improve efficiency. This is interesting from an overall system point of view. Society at large - or for example a regulator tasked primarily with incentivizing natural monopolies to improve efficiency - may not solely be interested that everyone becomes efficient. It may as well be important how the efficiency is improved, e.g. by reducing the use of imported or domestic resources or by laying off some particular types of labor and not other types. 116 Chapter 5 General Conclusions and Future Work 118 In this thesis, following a stakeholder’s perspective we have extended Counterfactual Analysis in several ways. In Chapter 1, we have provided a critical discussion on the different ingredients defining the problem of generating Group Counterfactual Explanations. In general, finding efficient solutions to this problem calls for solving Mixed Integer Nonlinear Programming formulations. We have related those to classic problems in the Continuous and Discrete Location Analysis literature, such as the p-median, the p-center, or the minimumsum-of-squared-distances problems. The link established between Counterfactual Analysis and Location Analysis may allow to transfer to the world of Counterfactual Analysis the well developed algorithmic machinery of Location Analysis. Conversely, unexplored versions of Location Analysis problems appear motivated by their counterparts in Counterfactual Analysis. In Chapter 2, we have formulated the problem of Group Counterfactual Analysis with all possible allocation rules as mathematical optimization models. We develop flexible and novel methodologies for score-based models, specifically detailed for liner models, additive trees models and neural networks. By calculating Group Counterfactual Analysis, we are able to detect: the critical features to perturb in a group setting to change the classifier’s prediction, outliers that do not behave like the rest of the group, and representative counterfactual explanations that cluster the original instances, when assigned to the same explanation. In Chapter 3, we extend the analysis of counterfactuals when data are functions. For both univariate and multivariate functional data, we are able to generate counterfactuals by combining existing instances in the dataset. With our methodology, we detect both the features that are to be perturbed, and the instances from the dataset that generate the explanations. In Chapter 4, we extend counterfactual explanations beyond classification. Specifically, we focus on a benchmarking framework, where data envelopment analysis is used to determine how efficient different firms are when comparing them to each other. We extend the concept of targets in benchmarking using the ideas of Counterfactual Analysis. By defining a bilevel optimization model, transformed into an equivalent single-level one employing the optimality conditions of the DEA model, we are capable of deploying the changes needed in order to improve the efficiency of firms, with the flexibility of including different ways to measure proximity. There are several promising interesting lines of future research for the Operations Research community. We now outline some of them below. The practical success of Counterfactual Analysis in Supervised Classification as developed in this thesis will be conditioned to three premises: the counterfactuals xbuilt for the tuple of instances x0should be feasible, the perturbations suggested on x0to yield xshould be doable, and counterfactuals, used as benchmarks for records labeled by the classifier used in the negative class should be indeed in the positive class -a condition which may differ from being classified as positive by the classifier -. Chapter 1 has given an overview of different judicious choices for the model ingredients to fulfill the three above-mentioned premises. However, if 120 2.14 Visualization of 100 of the counterfactual explanations for all the instances in the Taiwanese Credit dataset that comply with (2.41). The classifier considered is a logistic regression model. The explanations have been calculated solving model (ColCELR) with λind = 0,λglob = 0.1(left), 10 (right), and S∗=⌈0.95S⌉. The feature perturbations are displayed. . . . . . . . . . . . . 57 2.15 Counterfactual explanations for instances x0 sin Table 2.6 from the Students performance dataset for two classifiers, the logistic regression and the random forest. The explanations have been calculated solving model (CEsep) with Cs(x0 s,xs) = ∥x0 s−xs∥2 2+λind∥x0 s−xs∥0for λind = 0.02. The feature perturbations are displayed. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 2.16 Counterfactual explanations for instances x0 sin Table 2.6 from the Students performance dataset for two classifiers, the logistic regression and the random forest. The explanations have been calculated solving models (ColCELR) and (ColCEATM) with cost function (2.1) with λind = 0,λglob = 0.2. The feature perturbations are displayed. . . . . . . . . . . . . . . . . . . . . . . . 59 2.17 Features that need to be perturbed (in red) for all the instances in the Students performance dataset that comply with (2.41) when solving Problem (ColCELRhard) for all the values of Fmax to describe the Pareto front. The classifier considered is a logistic regression model. Fmax = 1,2,3,4,5yields an infeasible problem. 60 2.18 Counterfactual explanations for all the instances in the Students performance dataset that comply with (2.41). The classifier considered is a logistic regression model. The explanations have been calculated solving model (ColCELR) with λind = 0,λglob = 0.1(left), 10 (right), and S∗=⌈0.95S⌉. The feature perturbations are displayed. . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 2.19 Many-for-one counterfactual explanations with R= 3 for instances x0 1and x0 2in Table 2.3 from the Boston housing dataset where the classifier is the logistic regression model. The explanations have been calculated solving model (2.43)-(2.44) with λind = 0.01 and ν= 0.5. Each row is an explanation, calculated by adding one constraint from (2.45) in Xr(x0 s). Features perturbations are displayed. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 2.20 One-for-many and one-for-all analysis for x0from the Boston housing dataset as in (2.50) for the logistic regression and R= 3 and R= 1 respectively. We consider Problem (2.46)-(2.49) with ν= 0.5. The plot represents the feature values of the counterfactual explanations (below) and their associatedclusters(above). ............................... 62 127 2.21 One-for-many and one-for-all analysis for x0from the Boston housing dataset as in (2.50) for the logistic regression and R= 3 and R= 1 respectively. We consider Problem (2.46)-(2.49) with ν= 0.7. The plot represents the feature values of the counterfactual explanations (below) and their associatedclusters(above). ............................... 63 2.22 One-for-many and one-for-all analysis for x0from the Boston housing dataset as in (2.50) for the logistic regression and R= 3 and R= 1 respectively. We consider Problem (2.46)-(2.49) with ν= 0.9. The plot represents the feature values of the counterfactual explanations (below) and their associatedclusters(above). ............................... 63 2.23 Counterfactual explanations for the COMPAS dataset for a NN model when solving a sparse one-for-many allocation model with |S| = 10 and Jmax = 3. Each row represents a counterfactual instance and each column a feature. The value of the counterfactual features are illustrated, in green if its value is 1, in whiteifitis0. .................................. 65 2.24 Counterfactual explanations for the COMPAS dataset for a NN model when solving a sparse one-for-many allocation model with |S| = 20 and Jmax = 3. Each row represents a counterfactual instance and each column a feature. The value of the counterfactual features are illustrated, in green if its value is 1, in whiteifitis0.................................... 65 2.25 Number of times each of the features is perturbed for the COMPAS dataset for a NN, when solving the sparse one-for-many allocation model for |S| ∈ {10,20,50,100}and Jmax = 2,3. ........................ 66 2.26 Number of times each of the features is perturbed for the Students performance dataset for a NN, when solving the sparse one-for-many allocation model for |S| ∈ {10,20,50,100}and Jmax = 10,15,20................... 67 2.27 Number of times each of the features is perturbed for the German credit dataset for a NN, when solving the sparse one-for-many allocation model for |S| ∈ {10,20,50,100}and Jmax = 5,10,15. .................. 68 3.1 Comparison between different warping paths between the functions x(in blue) and x′(inorange). ................................ 75 3.2 Counterfactual explanations for x0of the ItalyPowerDemand dataset. Different values of Bmax, i.e., the number of prototypes used for the convex combination, have been imposed. The cost function is (3.5) with λ0= 0, λ2= 1, and ν= 0.5..................................... 81 3.3 Distance obtained vs number of prototypes used for the convex combination for the counterfactual explanation xfor x0of the ItalyPowerDemand dataset. . 82 128 3.4 Counterfactual explanations for x0of the ItalyPowerDemand dataset. Different values of Bmax, i.e., the number of prototypes used for the convex combination, have been imposed. The cost function is (3.6) with λ0= 0, λ2= 1, and ν= 0.5..................................... 83 3.5 Counterfactual explanations for x0of the NATOPS dataset for a Random Forest. Bmax = 1 prototype has been imposed to create the counterfactual. The cost function is (3.5) with λ0= 1, λ2= 0.005, and ν= 0.5. .......... 84 3.6 Changed features in the counterfactual explanation for x0of the NATOPS dataset for a Random Forest with Bmax = 1. The cost function is (3.5) with λ0= 1, λ2= 0.005, and ν= 0.5. ........................ 85 3.7 Counterfactual explanations for x0of the NATOPS datase for a Random Forest. Bmax = 2 prototypes has been imposed to create the counterfactual. The cost function is (3.5) with λ0= 1, λ2= 0.005, and ν= 0.5.............. 86 3.8 Changed features in the counterfactual explanation for x0of the NATOPS dataset for a Random Forest with Bmax = 2. The cost function is (3.5) with λ0= 1, λ2= 0.005, and ν= 0.5. ........................ 87 4.1 Directional search for an alternative production plan to (x,y)for a DMU along (dx,dy)using (DIR) in a benchmarking context. . . . . . . . . . . . . . . . . 94 4.2 Directional search in Interactive Benchmarking software. Real-world dataset of bank branches in Section 4.5. . . . . . . . . . . . . . . . . . . . . . . . . . 95 4.3 Counterfactual explanations for firms 3 and 4 in Tables 4.2 and 4.3 respectively imposing E∗= 0.8and different cost functions. . . . . . . . . . . . . . . . . . 100 4.4 Counterfactual Explanations for firm 238 from the Canadian bank dataset when solving Problem (CEDEA) and desired efficiency E∗= 1.........108 4.5 The inputs that change (in red) when solving Problem (CEDEA) with desired efficiency E∗= 1 for all the inefficient firms in the Canadian bank dataset. 111 4.6 Counterfactual Explanations for DMU 238 from the Canadian bank dataset when solving Problem (CEDEA) with desired efficiency E∗= 0.8........111 4.7 The inputs that change (in red) when solving Problem (CEDEA) with desired efficiency E∗= 0.8for all the inefficient firms in the Canadian bank dataset.113 129 List of Tables 2.1 Datasets used for the numerical experiments in Section 2.1.2 and 2.2, for the one-for-one, many-for-one, one-for-all and one-for-many allocation rule models. 29 2.2 Feature values of x0= (x0 1,...,x0 10)from the Boston housing dataset used for the models in Figures 2.1–2.2. . . . . . . . . . . . . . . . . . . . . . . 30 2.3 Feature values of x0= (x0 1,...,x0 10)from the Boston housing dataset used for the one-for-one models in Figures 2.5, 2.6 and 2.19. . . . . . . . . . . 33 2.4 Feature values of x0= (x0 1,...,x0 10)from the COMPAS dataset used for the modelsinFigures2.7-2.8. ............................ 34 2.7 Notation used to define the sparse one-for-many allocation rule model. . . . . . 41 2.9 Structure and Accuracy of the classifiers used for the numerical illustrations for the sparse one-for-many allocation rule model. . . . . . . . . . . . . . . . . 46 2.5 Feature values of x0= (x0 1,...,x0 10)from the Taiwanese credit dataset used for the models in Figures 2.11–2.12. . . . . . . . . . . . . . . . . . . . . 54 2.6 Feature values of x0= (x0 1,...,x0 10)from the Students performance dataset used for the models in Figures 2.15–2.16. . . . . . . . . . . . . . . . . 58 2.8 Summary of the datasets used for the experiments for the sparse one-for-many allocationrulemodel................................ 60 2.10 Evaluation of the sparse one-for-many allocation rule model with the column generation framework for the COMPAS dataset for the LR model for different values of |S| and Jmax. .............................. 60 2.11 Evaluation of the sparse one-for-many allocation rule model with the column generation framework for the COMPAS dataset for the NN for different values of |S| and Jmax................................... 60 2.12 Evaluation of the sparse one-for-many allocation rule model with the column generation framework for the Students performance dataset for the LR model for different values of |S| and Jmax..................... 64 2.13 Evaluation of the sparse one-for-many allocation rule model with the column generation framework for the Students performance dataset for the NN for different values of |S| and Jmax. ..................... 64 130 2.14 Evaluation of the sparse one-for-many allocation rule model with the column generation framework for the German credit dataset for the LR model for different values of |S| and Jmax. ......................... 64 2.15 Evaluation of the sparse one-for-many allocation rule model with the column generation framework for the German credit dataset for the NN for different values of |S| and Jmax............................ 64 3.1 Parameters and decision variables used to model Problem (3.1) for additive tree models, when data are functions. . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.1 Inputs, outputs and corresponding Farrell input-efficiency of 4 different firms when solving the classical DEA model with CRS. . . . . . . . . . . . . . . . . 98 4.2 Counterfactual explanations for firm 3 in Table 4.1 imposing E∗= 0.8and differentcostfunctions............................... 99 4.3 Counterfactual explanations for firm 4 in Table 4.1 imposing E∗= 0.8and differentcostfunctions............................... 99 4.4 Descriptive statistics of the Canadian bank dataset in Schaffnit et al. [1997]. 107 4.5 Value of the parameters ν0, ν1and ν2in (4.22) for the different cost functions used in the calculation of counterfactuals for the Canadian bank dataset. . . 107 4.6 Average results on how often the inputs (personnels) change when solving Problem (CEDEA) with desired efficiency E∗= 1 for all the inefficient firms in the Canadian bank dataset...........................110 4.7 Average results on how much the inputs (personnels) change when solving Problem (CEDEA) with desired efficiency E∗= 1 for all the inefficient firms in the Canadian bank dataset. ........................110 4.8 Average results on how much the inputs (personnels) change when solving Problem (CEDEA) with desired efficiency E∗= 0.8for all the inefficient firms in the Canadian bank dataset. .....................112 4.9 Average results on how much the inputs (personnels) change when solving Problem (CEDEA) with desired efficiency E∗= 0.8for all the inefficient firms in the Canadian bank dataset. .....................112 5.1 Description of the features of the Boston housing dataset. . . . . . . . . . 135 5.2 Description of the features of the COMPAS dataset. ...............135 5.3 Description of the features of the Taiwanese Credit dataset. . . . . . . . 135 5.4 Description of the features of the Students performance dataset. . . . . 136 5.5 Description of the features of the German Credit dataset. . . . . . . . . . . 136 131 Appendix 133 Description of datasets Variable Definition Type CRIM per capita crime rate by town numerical ZN proportion of residential land zoned for lots over 25,000 sq.ft numerical INDUS proportion of non-retail business acres per town numerical CHAS Charles River dummy variable (1 if tract bounds river; 0 otherwise) binary NOX nitric oxides concentration (parts per 10 million) numerical RM average number of rooms per dwelling numerical AGE proportion of owner-occupied units built prior to 1940 numerical DIS weighted distances to five Boston employment centres numerical RAD index of accessibility to radial highways numerical TAX full-value property-tax rate per $10,000 numerical PTRATIO pupil-teacher ratio by town numerical B1000(Bk −0.63)2where Bk is the proportion of blacks by town numerical LSTAT % lower status of the population numerical MEDV +1 if Median value of owner-occupied homes in $1000’s over the 50th percentile, −1otherwise target Table 5.1: Description of the features of the Boston housing dataset. Variable Definition Type Race race of defendants: African-American, Asian, Caucasian, Hispanic, Native American, Other categorical (nominal) AgeGroup age category of defendants: 1 (<25), 2 (25-45), 3 (>45) categorical (ordinal) Sex sex of defendants: 0 (male) or 1 (female) binary PriorsCount prior criminal records of defendants numerical ChargeDegree charge degree of defendants: 0 (Misdemeanor), 1 (Felony) binary TwoYearRecid binary variable indicate whether defendant is rearrested at within two years: 1 (no), -1 (yes) target Table 5.2: Description of the features of the COMPAS dataset. Variable Definition Type isMale sex of individual: 0 (female), 1(male) binary isMarried marital status: 0 (not married), 1 (married) binary AgeGroup age category: 1 (<25), 2 (25-39), 3 (40-59), 4 (≥60) categorical (ordinal) EducationLevel education level: 1 (other), 2 (highschool), 3 (university), 4 (graduate) categorical (ordinal) MaxBillAmountOverLast6Months maximum bill amount over the last 6 months numerical MaxPaymentAmountOverLast6Months maximum payment amount over the last 6 months numerical MonthsWithZeroBalanceOverLast6Months months with zero balance over the last 6 months numerical MonthsWithLowSpendingOverLast6Months months with low spending (bill<0.2) over last 6 months numerical MonthsWithHighSpendingOverLast6Months months with high spending (bill>0.8) over last 6 months numerical MostRecentBillAmount most recent bill amount numerical MostRecentPaymentAmount most recent payment amount numerical TotalOverdueCounts total overdue counts numerical TotalMonthsOverdue total months overdue numerical HasHistoryOfOverduePayments binary variable indicating whethere the individual has a history of overdue payment: yes(1), no (0) binary NoDefaultNextMonth binary variable indicating whether the individual will default next month (0) or not (1) target Table 5.3: Description of the features of the Taiwanese Credit dataset. 135 Variable Definition Type school student’s school: 0 (Gabriel Pereira), 1 (Mousinho da Silveira) binary sex student’s sex: 0 (male), 1 (female) binary age student’s age numerical address student’s home address type: 0 (rural), 1(urbal) binary famsize family size: 0 (less or equal to 3), 1 (greater than 3) binary Pstatus parent’s cohabitation status: 0 (living together), 1 (apart) binary Medu mother’s education: 0 (none), 1 (primary education), 2 (5th to 9th grade), 3 (secondary education), 4 (higher education) categorical (ordinal) Fedu father’s education: 0 (none), 1 (primary education), 2 (5th to 9th grade), 3 (secondary education), 4 (higher education) categorical (ordinal) Mjob mother’s job: ’teacher’, ’health’ care related, civil ’services’ (e.g. administrative or police), ’at_home’ or ’other’ categorical (nominal) Fjob father’s job: ’teacher’, ’health’ care related, civil ’services’ (e.g. administrative or police), ’at_home’ or ’other’ categorical (nominal) reason reason to choose this school: close to ’home’, school ’reputation’, ’course’ preference or ’other’ categorical (nominal) guardian student’s guardian: ’mother’, ’father’ or ’other’ categorical (nominal) traveltime home to school travel time: 1 (<15 min) , 2 (15 to 30 min), 3 (30 min to 1 hour), 4 (>1 hour) categorical (ordinal) studytime weekly study time: 1 (<2 hours), 2 (2 to 5 hours), 3 (5 to 10 hours), or 4 (>10 hours) categorical (ordinal) failures number of past class failures numerical schoolsup extra educational support: 0 (no), 1 (yes) binary famsup family educational support: 0 (no), 1 (yes) binary paid extra paid classes within the course subject: 0 (no), 1 (yes) binary activities extra-curricular activities: 0 (no), 1 (yes) binary nursery attended nursery school: 0 (no), 1 (yes) binary higher wants to take higher education : 0 (no), 1 (yes) binary internet Internet access at home: 0 (no), 1 (yes) binary romantic with a romantic relationship: 0 (no), 1 (yes) binary famrel quality of family relationships: from 1 (very bad) to 5(excellent) categorical (ordinal) freetime free time after school: from 1 (very low) to 5 (very high) categorical (ordinal) goout going out with friends: from 1 (very low) to 5 (very high) categorical (ordinal) Dalc workday alcohol consumption: from 1 (very low) to 5 (very high) categorical (ordinal) Walc weekend alcohol consumption: from 1 (very low) to 5 (very high) categorical (ordinal) health current health status: : from 1 (very low) to 5 (very high) categorical (ordinal) absences number of school absences numerical Class 1 if the student passes the math course, -1 otherwise target Table 5.4: Description of the features of the Students performance dataset. Variable Definition Type status_account status of existing account: A11 (<0 DM), A12 (0-200 DM), A13 (≥200 DM), A14 (no checking account) categorical credit_history credit history: A30 (no credits taken), A31 (all credits paid back duly), A32 (existing credits paid back), A33 (delay in paying off), A34 (critical account) categorical purpose purpose categorical savings savings: A61 (<100DM), A62 (100-500 DM), A63 (500-1000 DM), A64 (≥1000 DM), A65 (unknown) categorical employment present employment since A71 (unemployed), A72 (<1year), A73 (1-4 years), A74 (4-7 years), A75 (≥7years) categorical status_sex personal status and sex: A91 (male and divorced), A92 (female and divorced or married), A93 (male and single), A94 (male and married), A95 (female and single) categorical guarantors guarantors: A101 (none), A102 (co-applicant), A103 (guarantor) categorical property property: A121(real estate), A122 (life insurance), A123 (car or other), A124 (no property) categorical installment_plans other installment plans: A141 (bank), A142 (stores), A143 (none) categorical housing housing: A151 (rent), A152 (own), A153 (for free) categorical job job: A171 (unemployed), A172 (unskilled), A173 (skilled), A174 (highly qualified) categorical duration duration in months numerical amount credit amount numerical installment_rate installment rate in percentage of disposable income numerical residence_since present residence since numerical age age in years numerical number_credits number of existing credits at this bank numerical maintance_obligors number of people being liable to provide maintenance for numerical tel telephone: A191 (none), A192 (registered) binary foreign foreign worker: A201 (yes), A202 (no) binary Table 5.5: Description of the features of the German Credit dataset. 136