text:
\section{Introduction and related work}\label{sec:introduction}

Invariance, or equivariance to certain transformations of data, has proven essential in numerous applications of machine learning (ML). Invariance leads to better generalization capabilities of ML models~\citep{arjovsky2019invariant, chen2020simple, bloem2020probabilistic}.
For instance, in image recognition tasks, predictions ought to remain unchanged under scaling, translation, or rotation of the input image. Data augmentation is one of the earliest heuristics developed to promote this kind of invariance that has become indispensable in training successful models, like deep neural networks (DNNs)~\citep{shorten2019survey, NEURIPS2020_44feb009}. State-of-the-art architectures can handle invariance and equivariance requirements by design. Well-known examples include convolutional neural networks (CNNs) for translation invariance of image data~\citep{NIPS2012_c399862d}, group equivariant CNNs for other group transformations~\citep{cohen2016group}, recurrent neural networks (RNNs) and transformers for sequential data~\citep{vaswani2017attention}, DeepSet~\citep{NIPS2017_f22e4747} for sets, and graph neural networks (GNNs) for different types of geometric structures~\citep{battaglia2018relational}.

Many real-world applications in modern ML, however, call for an arguably stronger notion of invariance, called \emph{counterfactual invariance}. This is the case for image classification, algorithmic fairness~\citep{NIPS2016_9d268236, mitchell2021algorithmic}, robustness \citep{buhlmann2020invariance}, and out-of-distribution  generalization \citep{lu2021invariant}. These applications require predictors to exhibit invariance with respect to hypothetical manipulations of the data generating process (DGP)~\citep{peters2016causal,heinze2018invariant,rojas2018invariant,arjovsky2019invariant,buhlmann2020invariance}. In image classification, for instance, we want to build a model that ``would have made the same prediction, if the object position had been different with everything else being equal''. This model ought to be invariant to the DGP of some attributes, such as an object position and light condition. Similarly, in algorithmic fairness \citet{kilbertus2017avoiding,DBLP:conf/nips/KusnerLRS17} introduce notions of interventional and counterfactual fairness. These notions are based on the causal relationships between the observed variables of a model, and they require invariance to the DGP of protected attributes.

Inspired by problems in natural language processing (NLP), \citet{DBLP:conf/nips/VeitchDYE21} study a notion of counterfactual invariance. In their work, \citet{DBLP:conf/nips/VeitchDYE21} provide a method to enforce their notion of counterfactual invariance, based on the maximum mean discrepancy (MMD). However, this method promotes only a \emph{necessary}, but not sufficient observable imprint or proxy of counterfactual invariance during training. Furthermore, their work can only handle discrete random variables, and it applies to specific causal graphs. Motivated by this line of research, we study the following question:
\begin{center}
\emph{How do we learn predictors that are invariant\\under counterfactual changes of certain covariates?}
\end{center}

We propose a new definition of counterfactual invariance that generalizes the notion studied by \citet{DBLP:conf/nips/VeitchDYE21}. Furthermore, we propose a new method that directly enforces counterfactual invariance during the learning process. Unlike existing approaches, the proposed method is \emph{sufficient} for counterfactual invariance, and it works well for both categorical and continuous covariates and outcomes, as well as in multivariate settings

\paragraph{Technical overview.} Our notion of counterfactual invariance inevitably relies on strong assumptions about the DGP, as counterfactuals by their very nature cannot be observed. This is a major challenge, since counterfactual invariance cannot be checked, unless strong prior knowledge of the DGP is available. We circumvent this problem by deriving a connection between counterfactual invariance and conditional independence, under certain observability assumptions in the causal graph. The resulting conditional independence statements can then be enforced via a regularization term based on the flexible Hilbert-Schmidt Conditional Independence Criterion (\textsc{HSCIC}{}) \citep{Park20:CME}.

\section{Preliminaries}\label{sec:preliminaries}

\paragraph{Counterfactual invariance.}
We introduce structural causal models as in \citet{pearlj}.

\begin{definition}[Structural causal model (SCM)]
\label{def:SCM}
A structural causal model is a tuple $(\mathbf{U}, \mathbf{V}, F, \mathbb{P}_{\mathbf{U}})$ such that
\begin{enumerate}
\item $\mathbf{U}$ is a set of background variables that are exogenous to the model;
\item $\mathbf{V}$ is a set of observable (endogenous) variables; \item $F=\{f_{V}\}_{V\in\mathbf{V}}$ is a set of functions from (the domains of) $\mathsf{pa}(V) \cup U_V$ to (the domain of) $V$, where $U_V \subset \mathbf{U}$ and $\mathsf{pa}(V) \subseteq \mathbf{V} \setminus \{V\}$ such that $V = f_V(\mathsf{pa}(V), U_V)$; (iv) $\mathbb{P}_{\mathbf{U}}$ is a probability distribution over the domain of $\mathbf{U}$.
Further, the subsets $\mathsf{pa}(V) \subseteq \mathbf{V} \setminus \{V\}$ are chosen such that the graph $\mathcal{G}$ over $\mathbf{V}$ where the edge $V' \to V$ is in $\mathcal{G}$ if and only if $V' \in \mathsf{pa}(V)$ is a directed acyclic graph (DAG).
\end{enumerate}
\end{definition}
We always denote with $\mathbf{Y} \subset \mathbf{V}$ the outcome (or prediction target), and with $\hat{\mathbf{Y}}$ a predictor for that target. The predictor $\hat{\mathbf{Y}}$ is not strictly part of the SCM, because we get to tune it. However, since it takes inputs from $\mathbf{V}$, we often treat it as an observed variable in the SCM.
Each SCM implies a unique \emph{observational distribution} over $\mathbf{V}$ \citep{pearlj}.
However, SCMs allow us to also define interventional distributions. Given a variable $A \in \mathbf{V}$, an intervention $A \gets a$ amounts to replacing $f_A$ in $F$ with the constant function $A = a$. This yields a new SCM, which induces the \emph{interventional distribution}.
Similarly, we can intervene on multiple variables $\mathbf{V} \supseteq \mathbf{A} \gets \mathbf{a}$.
We then write $\mathbf{Y}^*_{\mathbf{a}}$ for the outcome in the intervened SCM, also called potential outcome. Note that the interventional distribution $\mathbb{P}_{\mathbf{Y}^*_{\mathbf{a}}}(\mathbf{y})$ differs in general from the conditional distribution $\mathbb{P}_{\mathbf{Y}\mid \mathbf{A}}(\mathbf{y} \mid \mathbf{a})$, due to possible unobserved confounding effects.
We can also condition on a set of variables $\mathbf{W} \subseteq \mathbf{V}$ in the (observational distribution of the) original SCM before performing an intervention, which we denote by $\mathbb{P}_{\mathbf{Y}^*_{\mathbf{a}} \mid \mathbf{W}} (\mathbf{y} \mid \mathbf{w})$.
This is a counterfactual distribution: ``Given that we have observed $\mathbf{W} = \mathbf{w}$, what would $\mathbf{Y}$ have been had we set $\mathbf{A} \gets \mathbf{a}$, instead of the values $\mathbf{A}$ had actually taken?''.
Note that the sets $\mathbf{A}$ and $\mathbf{W}$ need not be disjoint. We can now define counterfactual invariance.
\begin{definition}[Counterfactual invariance]\label{def:counterfactual_invariance}
Let $\mathbf{A}$, $\mathbf{W}$ be (not necessarily disjoint) sets of nodes in a given SCM. A predictor $\hat{\mathbf{Y}}$ is counterfactually invariant in $\mathbf{A}$ with respect to $\mathbf{W}$ if $\mathbb{P}_{\hat{\mathbf{Y}}^*_{\mathbf{a}} \mid \mathbf{W}}(\mathbf{y}\mid \mathbf{w}) = \mathbb{P}_{\hat{\mathbf{Y}}^*_{\mathbf{a}'} \mid \mathbf{W}}(\mathbf{y}\mid \mathbf{w})$ almost surely, for all $\mathbf{a}, \mathbf{a}'$ in the domain of $\mathbf{A}$, and all $\mathbf{w}$ in the domain of $\mathbf{W}$.\footnote{With a mild abuse of notation, if $\mathbf{W} = \emptyset$ then the requirement of conditional counterfactual invariance becomes $\mathbb{P}_{\hat{\mathbf{Y}}^*_{\mathbf{a}}}(\mathbf{y}) = \mathbb{P}_{\hat{\mathbf{Y}}^*_{\mathbf{a}'}}(\mathbf{y})$ almost surely, for all $\mathbf{a}, \mathbf{a}'$ in the domain of $\mathbf{A}$.}
\end{definition}
A predictor that is counterfactually invariant can thus be viewed as robust to changes of $\mathbf{A}$, in the sense that the (conditional) post-interventional distribution of $\hat{\mathbf{Y}}$ does not change for different values of the intervention.
In other words, the classification would not change had we gone back and changed the position of the object.
This definition is stronger than previous definitions of counterfactual invariance (see Definition 1.1 by \citet{DBLP:conf/nips/VeitchDYE21}). In fact, the invariance notion by \citet{DBLP:conf/nips/VeitchDYE21} only requires $\mathbb{P}_{\hat{\mathbf{Y}}^*_{\mathbf{a}}}(\mathbf{y}) = \mathbb{P}_{\hat{\mathbf{Y}}^*_{\mathbf{a}'}}(\mathbf{y})$ almost surely, for all $\mathbf{a}, \mathbf{a}'$ in the domain of $\mathbf{A}$.

\paragraph{Kernel mean embeddings and conditional measures.}
Our regularization term heavily relies on kernel mean embeddings (KMEs).
We now highlight the main concepts pertaining KMEs and refer the reader to \citet{smola2007hilbert,learning_kernels_book,prob_stats_book,Muandet17:KME} for more details.
Fix a measurable space $\mathscr{Y}$ with respect to a $\sigma$-algebra $\mathcal{F}_{\mathscr{Y}}$, and consider a probability measure $\mathbb{P}$ on the space $(\mathscr{Y}, \mathcal{F}_{\mathscr{Y}})$. Let $\mathcal{H}$ be a reproducing kernel Hilbert space (RKHS) with a bounded kernel $k_\mathbf{Y} \colon \mathscr{Y} \times \mathscr{Y} \rightarrow \mathbb{R}$, i.e., $k_\mathbf{Y}$ is such that $\sup_{\mathbf{y} \in \mathscr{Y}}k(\mathbf{y},\mathbf{y}) < +\infty$. The kernel mean embedding $\mu_{\mathbb{P}}$ of $\mathbb{P}$ is defined as the expected value of the function $k(\ \cdot \ ,\mathbf{y})$ with respect to $\mathbf{y}$, i.e., $\mu_{\mathbb{P}} \coloneqq \expect{}{k(\ \cdot \  ,\mathbf{y})}$.

The definition of KMEs can be extended to conditional distributions \citep{DBLP:journals/jmlr/FukumizuSG13,DBLP:conf/icml/GrunewalderLGBPP12,DBLP:conf/icml/SongHSF09,DBLP:journals/spm/SongFG13}.
Consider two random variables $\mathbf{Y}$, $\mathbf{Z}$, and denote with $(\Omega_\mathbf{Y}, \mathcal{F}_{\mathbf{Y}})$ and $(\Omega_\mathbf{Z}, \mathcal{F}_{\mathbf{Z}})$ the respective measurable spaces.
These random variables induce a probability measure $\mathbb{P}_{\mathbf{Y}, \mathbf{Z}}$ in the product space $\Omega_\mathbf{Y} \times \Omega_\mathbf{Z}$.
Let $\mathcal{H}_\mathbf{Y}$ be a RKHS with a bounded kernel $k_\mathbf{Y} (\cdot, \cdot )$ on $\Omega_\mathbf{Y}$.
We define the KME of a conditional distribution $\mathbb{P}_{\mathbf{Y}\mid \mathbf{Z}}(\cdot \mid \mathbf{z})$ via $\mu_{\mathbf{Y}\mid \mathbf{Z} = \mathbf{z}} \coloneqq \expect{}{k_\mathbf{Y}(\ \cdot \ ,\mathbf{y})\mid \mathbf{Z} = \mathbf{z}}$.
Here, the expected value is taken over $\mathbf{y}$.
KMEs of conditional measures can be estimated from samples. To illustrate this, consider i.i.d.\ samples $(\mathbf{y}_1, \mathbf{z}_1), \dots, (\mathbf{y}_n, \mathbf{z}_n)$.
Denote with $\hat{K}_\mathbf{Y}$ the kernel matrix with entries $[\hat{K}_\mathbf{Y}]_{i,j} \coloneqq k_\mathbf{Y}(\mathbf{y}_i, \mathbf{y}_j)$.
Furthermore,  let $k_\mathbf{Z}$ be a bounded kernel on $\Omega_\mathbf{Z}$.
Then, $\mu_{\mathbf{Y}\mid \mathbf{Z} = \mathbf{z}}$ can be estimated as
\begin{equation}\label{eq:cond_measure_estimand}
 \hat{\mu}_{\mathbf{Y}\mid \mathbf{Z} = \mathbf{z}} \coloneqq  \sum_{i = 1}^n \hat{w}_{\mathbf{Y}\mid \mathbf{Z}}^{(i)}(\mathbf{z})k_\mathbf{Y}(\cdot, \mathbf{y}_i)\; \text{ with }\; \hat{w}_{\mathbf{Y}\mid \mathbf{Z}}(\cdot) \coloneqq ( \hat{K}_\mathbf{Z} - n \lambda I )^{-1} \left [k_\mathbf{Z}(\cdot,\mathbf{z}_1), \cdots , k_\mathbf{Z}(\cdot,\mathbf{z}_n)\right ]^T
\end{equation}
where, $I$ is the identity matrix and $\lambda$ is a regularization parameter.
Here, $\hat{w}_{\mathbf{X}\mid \mathbf{A}}^{(i)}(\cdot)$, the $i$-th entry of $\hat{w}_{\mathbf{X}\mid \mathbf{A}}(\cdot)$, are the coefficients of kernel ridge regression \citep{DBLP:conf/icml/GrunewalderLGBPP12}.

\section{Counterfactually invariant predictors}\label{sec:ci}
\paragraph{Counterfactual invariance from conditional independence.}
In this section, we establish a simple graphical criterion to express the desired counterfactual invariance as a conditional independence in the observational distribution, rendering it estimable from i.i.d.\ data.
We first repeat the notion of \emph{blocked paths} \citep{pearlj}.
\begin{restatable}{definition}{blockedpath}
\label{def:blocked_path}
Consider a path $\pi$ of causal graph $\mathcal{G}$. A set of nodes $\mathbf{Z}$ blocks $\pi$, if $\pi$ contains a triple of consecutive nodes connected in one of the following ways: $N_i \rightarrow Z \rightarrow N_j $, $N_i \gets Z \rightarrow N_j $, or $N_i \rightarrow M \gets N_j $, with $N_i, N_j \notin \mathbf{Z}$, $Z \in \mathbf{Z}$, and neither $M$ nor any descendent of $M$ is in $\mathbf{Z}$.
\end{restatable}
Using \ref{def:blocked_path}, we give the following criterion for counterfactual invariance.
\begin{restatable}{theorem}{condcounterfactual}
\label{lemma:conditiona_counterfactual}
Let $\mathcal{G}$ be a causal diagram, and consider two (not necessarily disjoint) sets of nodes $\mathbf{A}$, $\mathbf{W}$.
Let $\mathbf{Z}$ be a set of nodes that blocks all non-directed paths from $\mathbf{A}\cup \mathbf{W}$ to $\mathbf{Y}$.
Then, for any SCM compatible with $\mathcal{G}$, any predictor $\hat{\mathbf{Y}}$ that satisfies $\hat{\mathbf{Y}}\perp \!\!\! \perp \mathbf{A}\mid \mathbf{Z}$ is counterfactually invariant in $\mathbf{A}$ with respect to $\mathbf{W}$.
\end{restatable}

\begin{figure}[t]
  \centering
\begin{tikzpicture}[node distance=7mm and 10mm, thick, main/.style = {draw, circle, minimum size=0.8cm}, >={triangle 45}]
    \node[main] (1) {$\scriptstyle \mathbf{A}$};
    \node[main] (2) [right =of 1] {$\scriptstyle \mathbf{X}$};
    \node[main] (3) [right =of 2] {$\scriptstyle \mathbf{Y}$};
    \draw[->] (1) -- (2);
    \draw[->] (2) -- (3);
    \draw[->] (1) to [out=320, in=220, looseness=1] (3);
    \draw[<->, dashed] (1) to [out=45, in=135, looseness=1] (2);
    \node (label) [above =of 1] {\textbf{(a)}};
  \end{tikzpicture}%
  \hfill
  \begin{tikzpicture}[node distance=7mm and 10mm, thick, main/.style = {draw, circle, minimum size=0.8cm}, >={triangle 45}]
    \node[main] (1) {$\scriptstyle \mathbf{A}$};
    \node[main] (2) [right =of 1] {$\scriptstyle \mathbf{X}$};
    \node[main] (3) [right =of 2] {$\scriptstyle \mathbf{Y}$};
    \node[main] (4) [above =of 2] {$\scriptstyle \mathbf{C}$};
    \draw[->] (1) -- (2);
    \draw[->] (2) -- (3);
    \draw[->] (4) -- (2);
    \draw[->] (4) -- (3);
    \draw[->] (1) to [out=320, in=220, looseness=1] (3);
    \node (label) [left =of 4] {\textbf{(b)}};
  \end{tikzpicture}%
  \hfill
    \begin{tikzpicture}[node distance=7mm and 10mm, thick, main/.style = {draw, circle, minimum size=0.8cm}, >={triangle 45}]
    \node[main] (1) {$\textbf{A}$};
    \node[main, dashed] (3) [right =of 1] {$\scriptstyle \textbf{U}$};
    \node[main] (5) [above =of 3] {$\scriptstyle \textbf{C}$};
    \node[main] (4) [right =of 3] {$\scriptstyle \textbf{Y}$};
    \node[main] (6) [above =of 4] {$\scriptstyle \textbf{X}$};
    \draw[->, dashed] (1) -- (3);
    \draw[->] (5) -- (1);
    \draw[->, dashed] (5) -- (3);
    \draw[->] (5) -- (4);
    \draw[->, dashed] (3) -- (4);
    \draw[->] (6) -- (4);
    \draw[->] (1) to [out=320, in=220, looseness=1] (4);
    \node (label) [left =of 5] {\textbf{(c)}};
  \end{tikzpicture}%
  \caption{%
    \textbf{(a) An example for \ref{lemma:conditiona_counterfactual}}, where $\mathbb{P}_{\mathbf{Y}^*_\mathbf{a}}(\mathbf{y})$ is unidentifiable, due to the confounding effect denoted by the dashed double arrow. However, by \ref{lemma:conditiona_counterfactual} any predictor $\hat{\mathbf{Y}}$ such that $\hat{\mathbf{Y}}\perp \!\!\! \perp \mathbf{A}\mid \mathbf{Z}$ with $\mathbf{Z} = \{\mathbf{X}\}$ is counterfactually invariant in $\mathbf{A}$. Hence, the criterion in \ref{lemma:conditiona_counterfactual} is weaker than identifiability.
    \textbf{(b) Assumed causal structure for the Adult dataset}, where $\mathbf{A}$ consists of the protected attributes \emph{gender} and \emph{age}. A predictor is learned for the income $\mathbf{Y}$ that is invariant in $\mathbf{A}$, which corresponds to counterfactual fairness (see \ref{def:counterfactual_fairness}).
    \textbf{(c) Causal structure for the constructed dSprites ground truth}, where $\mathbf{A} = \{\text{\small{Pos.X}}\}$, $\mathbf{U} = \{\text{\small{Scale}}\}$, $\mathbf{C} = \{\text{\small{Shape}}, \text{\small{Pos.Y}}\}$,
    $\mathbf{X} = \{\text{\small{Color}}, \text{\small{Orientation}}\}$, and $\mathbf{Y} = \{\text{\small{Outcome}}\}$.
   $\mathbf{U}$ is unobserved.
    We want a predictor for $\mathbf{Y}$, which is counterfactually invariant in $\mathbf{A}$ with respect to all other observed attributes (see \ref{sec:counterfactual_invariance}).
    }\label{fig:graphs}
\end{figure}

The proof is deferred to Appendix \ref{proof:conditiona_counterfactual}, where we derive the graphical criterion by studying a graphical model similar to the counterfactual graph \citep{DBLP:journals/jmlr/ShpitserP08}.
Our key observation is that the set $\mathbf{Z}$ as in \ref{lemma:conditiona_counterfactual} acts as a valid adjustment set for certain random variables in this joint pre- and post-interventional graphical model.
Using this observation, we derive \ref{theorem:identifiability_app} with similar techniques to covariate adjustment. However, our proof does \emph{not} achieve identification of the counterfactual distributions by simply applying the do-calculus \citep{pearlj}.
Crucially, \emph{the existence of a measurable set $\mathbf{Z}$ as in Theorem \ref{lemma:conditiona_counterfactual} does not imply the identifiability of counterfactual distributions $\mathbb{P}_{\mathbf{Y}^*_{\mathbf{a}}} (\mathbf{y})$} (see Figure \ref{fig:graphs}(a) for a counterexample). In particular,  Theorem \ref{lemma:conditiona_counterfactual} does not rule out hidden confounding in the model. Hence, our method is applicable when (certain parameters of) $\mathbb{P}_{\mathbf{Y}^*_{\mathbf{a}}} (\mathbf{y})$ or $\mathbb{P}_{\mathbf{Y}^*_{\mathbf{a}} \mid \mathbf{W}} (\mathbf{y})$ cannot be learned from observational data.

\paragraph{Counterfactually invariant predictors.}
We now propose our novel regularization term denoted by $\hscic{\hat{\mathbf{Y}},\mathbf{A}}{\mathbf{Z}}$ that promotes counterfactual invariance by encouraging the conditional independence $\hat{\mathbf{Y}} \perp \!\!\! \perp \mathbf{A} \mid \mathbf{Z}$ (see Theorem \ref{lemma:conditiona_counterfactual}). Given a loss function $\mathcal{L}(\hat{\mathbf{Y}})$ for the predictive performance of $\hat{\mathbf{Y}}$, we propose to minimize the following loss
\begin{equation}\label{eq:total_loss}
    \mathcal{L}_{\textsc{total}}(\hat{\mathbf{Y}}) = \mathcal{L}(\hat{\mathbf{Y}}) + \gamma \cdot\textsc{HSCIC}{}(\hat{\mathbf{Y}},\mathbf{A}\mid \mathbf{Z})\:,
\end{equation}
with a regularization parameter $\gamma \ge 0$.
We will now develop and justify the regularization term (see  Corollary~\ref{cor:hscic_independence}) for which $\textsc{HSCIC}(\hat{\mathbf{Y}},\mathbf{A}\mid \mathbf{Z}) = 0$ if and only if $\hat{\mathbf{Y}}\perp \!\!\! \perp \mathbf{A}\mid \mathbf{Z}$.

\paragraph{Hilbert-Schmidt Conditional Independence Criterion (\textsc{HSCIC}{}).}
Consider two random variables $\mathbf{Y}$ and $\mathbf{A}$, and denote with $(\Omega_\mathbf{Y}, \mathcal{F}_{\mathbf{Y}})$ and $(\Omega_\mathbf{A}, \mathcal{F}_{\mathbf{A}})$ the respective measurable spaces. Suppose that we are given two RKHSs $\mathcal{H}_{\mathbf{Y}}$, $\mathcal{H}_{\mathbf{A}}$ over the support of $\mathbf{Y}$ and $\mathbf{A}$ respectively. The tensor product space $\mathcal{H}_{\mathbf{Y}} \otimes \mathcal{H}_{\mathbf{A}}$ is defined as the space of functions of the form $(f \otimes g)(\mathbf{y}, \mathbf{a}) \coloneqq f(\mathbf{y}) g(\mathbf{a})$, for all $f \in \mathcal{H}_{\mathbf{Y}}$ and $g \in \mathcal{H}_{\mathbf{A}}$. The tensor product space yields a natural RKHS structure, with kernel $k$ defined by $k(\mathbf{y} \otimes \mathbf{a}, \mathbf{y}' \otimes \mathbf{a}') \coloneqq k_{\mathbf{Y}}(\mathbf{y}, \mathbf{y}') k_{\mathbf{A}}(\mathbf{a}, \mathbf{a}')$. We refer the reader, i.e., to \citet{DBLP:journals/jmlr/SzaboS17} for more details on tensor product spaces.
With this notation we define:
\begin{definition}[\textsc{HSCIC}{}]
\label{def:definition_hscic}
For (sets of) random variables $\mathbf{Y}$, $\mathbf{A}$, $\mathbf{Z}$ the \textsc{HSCIC}{} between $\mathbf{Y}$ and $\mathbf{A}$ given $\mathbf{Z}$ is defined as the real-valued random variable $\hscic{\mathbf{Y},\mathbf{A}}{\mathbf{Z}} = H_{\mathbf{Y},\mathbf{A}\mid \mathbf{Z}} \circ \mathbf{Z}$. Here, $H_{\mathbf{Y},\mathbf{A}\mid \mathbf{Z}}$ is a real-valued deterministic function, defined as
\[
H_{\mathbf{Y},\mathbf{A}\mid \mathbf{Z}}(\mathbf{z}) \coloneqq \left \| \mu_{\mathbf{Y},\mathbf{A}\mid \mathbf{Z} = \mathbf{z}} - \mu_{\mathbf{Y}\mid \mathbf{Z} = \mathbf{z}}\otimes \mu_{\mathbf{A}\mid \mathbf{Z} = \mathbf{z}} \right \|,
\]
with $\left \| \cdot \right \|$ the metric induced by the inner product of the tensor product space $\mathcal{H}_{\mathbf{X}}\otimes \mathcal{H}_{\mathbf{A}}$.
\end{definition}
The \textsc{HSCIC}{} has the following important property.
\begin{restatable}{theorem}{hscicind}
\label{thm:hscic_independence}
With the notation introduced above, suppose that the kernel $k$ of the tensor product space $\mathcal{H}_{\mathbf{X}}\otimes \mathcal{H}_{\mathbf{A}}$ is characteristic.\footnote{The tensor product kernel $k$ is characteristic if the mapping $\mathbb{P}_{\mathbf{Y}, \mathbf{A}} \mapsto \expect{}{k(\ \cdot \ , \mathbf{y}\otimes \mathbf{a})}$ is injective. Here, the expected value is taken over the pair $(\mathbf{y}, \mathbf{a})$.} Then, $\hscic{\mathbf{Y},\mathbf{A}}{\mathbf{Z}} = 0$ almost surely if and only if $\mathbf{Y}\perp \!\!\! \perp \mathbf{A} \mid \mathbf{Z}$.
\end{restatable}
We provide a proof in Appendix \ref{appendix:hscic_independence}.
The \textsc{HSCIC}{} as a criterion for conditional independence is not a novel concept \citep{Park20:CME}.
However, our definition differs from \citet{Park20:CME}, who define the \textsc{HSCIC}{} via KMEs based on the Bochner conditional expected value.
In contrast, we define the \textsc{HSCIC}{} using the embeddings of conditional measures $\mu_{\mathbf{Y}\mid \mathbf{Z} = \mathbf{z}}$. This allows us to bypass some technical assumptions required by \citet{Park20:CME}, while retaining generality over cross-covariance based methods \citep{DBLP:conf/nips/FukumizuGSS07}.
We refer readers to Appendix \ref{appendix:ind_bochner} and \ref{appendix:ind_cross_covariance} for a comparison with previous approaches.
Combining  Theorem \ref{lemma:conditiona_counterfactual} with Theorem \ref{thm:hscic_independence}, we can now use the \textsc{HSCIC}{} to promote counterfactual invariance.
\begin{corollary}
\label{cor:hscic_independence}
Consider an SCM with causal diagram $\mathcal{G}$ and fix two (not necessarily disjoint) sets of nodes $\mathbf{A}$, $\mathbf{W}$. Let $\mathbf{Z}$ be a set of nodes in $\mathcal{G}$ that blocks all non-directed paths from $\mathbf{A}\cup \mathbf{W}$ to $\mathbf{Y}$. Then, any predictor $\hat{\mathbf{Y}}$ that satisfies $\hscic{\hat{\mathbf{Y}},\mathbf{A}}{\mathbf{Z}} = 0$ almost surely is counterfactually invariant in $\mathbf{A}$ with respect to $\mathbf{W}$.
\end{corollary}

\paragraph{Estimating the \textsc{HSCIC}{} from samples.}
We show how to estimate $\hscic{\mathbf{Y},\mathbf{A}}{\mathbf{Z}}$ from finite samples using eq.~\eqref{eq:cond_measure_estimand}.
Given $n$ samples $\{(\mathbf{y}_i, \mathbf{a}_i, \mathbf{z}_i) \}_{i = 1}^n$, denote by $\hat{K}_{\mathbf{Y}}$ and $\hat{K}_{\mathbf{A}}$ the corresponding kernel matrices for $\mathbf{Y}$ and $\mathbf{A}$. We then estimate the $H_{\mathbf{Y},\mathbf{A}\mid \mathbf{X}}$ as \citep{Park20:CME}
\begin{align}
    \hat{H}^2_{\mathbf{Y},\mathbf{A}\mid \mathbf{Z}}(\cdot) & = \hat{w}^T_{\mathbf{Y},\mathbf{A}\mid \mathbf{Z}}(\cdot) \left ( \hat{K}_{\mathbf{Y}} \odot \hat{K}_{\mathbf{A}}\right ) \hat{w}_{\mathbf{Y},\mathbf{A}\mid \mathbf{Z}}(\cdot) - 2 \left ( \hat{w}^T_{\mathbf{Y} \mid \mathbf{Z}}(\cdot) \hat{K}_{\mathbf{Y}} \hat{w}_{\mathbf{Y},\mathbf{A} \mid \mathbf{X}}(\cdot)  \right )\label{eq:estimate_hscic} \\
    & \cdot \left ( \hat{w}^T_{\mathbf{A} \mid \mathbf{X}}(\cdot) \hat{K}_{\mathbf{A}} \hat{w}_{\mathbf{Y},\mathbf{A} \mid \mathbf{Z}}(\cdot)  \right ) + \left ( \hat{w}^T_{\mathbf{Y} \mid \mathbf{Z}}(\cdot) \hat{K}_{\mathbf{Y}} \hat{w}_{\mathbf{Y} \mid \mathbf{Z}}(\cdot)  \right ) \left ( \hat{w}^T_{\mathbf{A} \mid \mathbf{Z}}(\cdot) \hat{K}_{\mathbf{A}} \hat{w}_{\mathbf{A} \mid \mathbf{Z}}(\cdot)  \right )\;, \nonumber
\end{align}
where $\odot$ is element-wise multiplication.
As mentioned, the vectors $\hat{w}_{\mathbf{Y}\mid \mathbf{Z}}(\cdot)$, $\hat{w}_{\mathbf{A}\mid \mathbf{Z}}(\cdot)$, and $\hat{w}_{\mathbf{Y},\mathbf{A}\mid \mathbf{Z}}(\cdot)$ are found via kernel ridge regression.
Theorem 1 by \citet{DBLP:journals/focm/CaponnettoV07} shows that if the regularization parameters for the kernel ridge coefficients decay
to $0$ sufficiently quickly as $n \to \infty$, then $\hat{H}^2_{\mathbf{Y},\mathbf{A}\mid \mathbf{Z}}(\cdot)$ convergences to $H^2_{\mathbf{Y},\mathbf{A}\mid \mathbf{Z}}(\cdot)$ at a rate of $n^{-c}$, where $c$ is a parameter depending on the complexity of the embeddings and the effective dimensions of $\mathcal{H}_{\mathbf{Y}}$ and $\mathcal{H}_{\mathbf{A}}$.
In practice, computing the \textsc{HSCIC}{} approximation $\hat{H}^2_{\mathbf{Y},\mathbf{A}\mid \mathbf{Z}}(\cdot)$ may be time-consuming.
To speed it up, we can use random Fourier features to approximate the matrices $\hat{K}_{\mathbf{Y}}$ and $\hat{K}_{\mathbf{A}}$ as above \citep{DBLP:conf/nips/RahimiR07,DBLP:conf/icml/AvronKMMVZ17}.
We refer readers to Appendix \ref{appendix:random_fourier_features} for a discussion of the use of random Fourier features.

\paragraph{Measuring counterfactual invariance.}
Besides predictive performance, e.g., mean squared error (\textbf{MSE}) for regression, our key metric of interest is the level of counterfactual invariance achieved by the predictor $\hat{\mathbf{Y}}$.
We quantify the overall counterfactual variance as a single scalar
\begin{equation}\label{eq:vcf}
  \textbf{VCF}(\hat{\mathbf{Y}}) = \mathbb{E}_{\mathbf{w} \sim\mathbb{P}_{\mathbf{W}}}\left[ \var_{\mathbf{a}' \sim \mathbb{P}_{\mathbf{A}}}[\mathbb{E}_{\hat{\mathbf{Y}}^*_{\mathbf{a}'} \mid \mathbf{W}}(\hat{\mathbf{y}} \mid \mathbf{w})]\right]\:.
\end{equation}
For deterministic predictors (i.e., point estimators), which we use in all our experiments, the variance term is zero if and only if counterfactual invariance holds at $\mathbf{w}$ for (almost) all $\mathbf{a}'$ in the support of $\mathbb{P}_{\mathbf{A}}$ (i.e., the prediction remains constant).
Since the variance is non-negative, $\mathbf{VCF}$ (we often drop the argument) is zero if and only if counterfactual invariance holds almost surely.
However, the outer expectation of $\mathbf{VCF}$ may still be sensitive to outliers in its arguments, which is why it can be instructive to also analyze the variances for different conditioning values $\mathbf{Z}=\mathbf{z}$.
To estimate \textbf{VCF} in practice, we pick $d$ datapoints $(\mathbf{w}_i)_{i=1}^d$ from the observed data, and for each compute the counterfactual outcomes $\hat{\mathbf{Y}}^*_{\mathbf{a}'}\mid \mathbf{w}_i$ for $k$ different values of $\mathbf{a}'$.
The inner expectation is simply the predictor output.
We use empirical variances with $k$ examples for each of the $d$ chosen datapoints, and the empirical mean
of the $d$ variances for the outer expectation.
Finally, we note that by \ref{cor:hscic_independence}, \textsc{HSCIC}{} itself also serves as a useful measure of how well counterfactual invariance has been enforced.

\begin{wrapfigure}{R}{0.45\textwidth}
\vspace{-3mm}
\centering
\includegraphics[width=0.35\textwidth]{dsprites_pic.png}
\vspace{1mm}
\caption{\textbf{Examples of images from the dSprites dataset} \citep{dsprites17}. This dataset consists of 2D shapes generated from 6 ground truth independent latent factors.}
\label{fig:dspritesexamples}
\vspace{-5mm}
\end{wrapfigure}

\section{Applications of counterfactual invariance}\label{sec:counterfactual_invariance}
To clarify the benefits and challenges of counterfactual invariance, we now provide brief explanations of potential application areas, which we will subsequently study empirically with real-world data.

\paragraph{Image classification.} Counterfactual invariance serves as a strong notion of robustness in high-dimensional settings, such as image classification: ``Would the truck have been classified correctly had it been winter in this exact situation instead of summer?'' For concrete demonstration, we use the dSprites dataset \citep{dsprites17}, which consists of relatively simple, yet high-dimensional, square black and white images of different shapes (squares, ellipses, etc.), sizes, and orientations (rotation) in different xy-positions (see Figure \ref{fig:dspritesexamples}).
Since this dataset is fully synthetic and labelled (we know all factors for each image), we consider a causal model as depicted in Figure \ref{fig:graphs}(c).
The full structural equations are provided in Appendix \ref{app:sec:imageresults}.
The main goal in this application is to evaluate our method for high-dimensional data, including mixed categorical (e.g., shapes) and continuous (e.g., position) variables.
Therefore, we set the outcome $\textbf{Y}$ to be a function of all observed and unobserved variables, which can be inferred from the high-dimensional images.

\paragraph{Counterfactual fairness.}
We first show that the popular definition of counterfactual fairness \citep{DBLP:conf/nips/KusnerLRS17} is an incarnation of counterfactual invariance (see Definition \ref{def:counterfactual_invariance}).
Counterfactual fairness corresponds to the commonly posed question, after receiving a consequential decision from an algorithm: ``Would I have gotten the same outcome had I been a different gender, race, or age with all else being equal?''.
Again, we denote the outcome by $\mathbf{Y} \subset \mathbf{V}$ and so-called \emph{protected attributes} by $\mathbf{A} \subseteq \mathbf{V} \setminus \mathbf{Y}$. These attributes typically constitute demographic features such as gender, race, or age protected under anti-discrimination laws (see, e.g., \citet{us_law}).
Collecting all remaining observed covariates into $\textbf{W} \coloneqq \mathbf{V} \setminus \mathbf{Y}$, we define counterfactual fairness as follows, showing that it is a special case of counterfactual invariance.
\begin{definition}[Counterfactual Fairness, Definition 5 by \citet{DBLP:conf/nips/KusnerLRS17}]
\label{def:counterfactual_fairness}
A predictor $\hat{\mathbf{Y}}$ is counterfactually fair if under any context $\mathbf{W} = \mathbf{w}$ and $\mathbf{A} = \mathbf{a}$, it holds $\mathbb{P}_{\hat{\mathbf{Y}}^*_{\mathbf{a}} \mid \mathbf{W}, \mathbf{A}}(\mathbf{y}\mid \mathbf{w}, \mathbf{a}) = \mathbb{P}_{\hat{\mathbf{Y}}^*_{\mathbf{a}'} \mid \mathbf{W}, \mathbf{A}}(\mathbf{y}\mid \mathbf{w}, \mathbf{a})$, for all $\mathbf{y}$ and for any value $\mathbf{a}'$ attainable by $\mathbf{A}$.
\end{definition}
For a concrete example, we consider the commonly used ``Adult'' dataset from the UCI repository \citep{uci_dataset}.
Our goal is to predict whether an individual's income is above a certain threshold $\mathbf{Y} \in \{0,1\}$ based on a collection of (demographic) information including protected attributes such as gender and age.
From a fairness perspective, such a predictor could be used in loan-granting decisions, where one must not discriminate against certain demographic groups. We follow \citet{nabi2018fair,DBLP:conf/aaai/Chiappa19}, where a subset of variables are selected from the dataset and a causal structure is assumed as in Figure \ref{fig:graphs}(b).
We choose gender (considered binary in this dataset) and age (considered continuous) as the protected attributes $\mathbf{A}$. We denote the marital status, level of education, occupation, working hours per week, and work class jointly by $\mathbf{X}$ and combine the remaining observed attributes in $\mathbf{C}$.
With this notation, our goal is to learn a predictor $\hat{\mathbf{Y}}$ that is counterfactually invariant (or fair) in $\mathbf{A}$ with respect to $\mathbf{W} = \mathbf{C} \cup \mathbf{X}$.
We remark that achieving fairness for continuous or even mixed categorical and continuous protected attributes is an ongoing line of research (even for non-causal fairness notions) \citep{Mary19:Fairness,chiappa2021fairness}, but directly supported by \textsc{HSCIC}{}.

\section{Experiments}\label{sec:new_experiments}

We begin our empirical assessment of \textsc{HSCIC}{} on synthetic datasets, where ground truth is known and systematically study how performance is influenced by the dimensionality of the variables. We then move to the dSprites and Adult datasets described in \ref{sec:counterfactual_invariance}.

\begin{wrapfigure}{r}{0.4\textwidth}
\vspace{-5mm}
\centering
   \begin{tikzpicture}[node distance=7mm and 10mm, thick, main/.style = {draw, circle, minimum size=0.8cm}, >={triangle 45}]
    \node[main] (1) {$\scriptstyle \mathbf{A}$};
    \node[main] (2) [right =of 1] {$\scriptstyle \mathbf{X}$};
    \node[main] (3) [right =of 2] {$\scriptstyle \mathbf{Y}$};
    \node[main] (4) [above =of 2] {$\scriptstyle \mathbf{Z}$};
    \draw[->] (1) -- (2);
    \draw[->] (2) -- (3);
    \draw[->] (4) -- (1);
    \draw[->] (4) -- (2);
    \draw[->] (4) -- (3);
    \draw[->] (1) to [out=320, in=220, looseness=1] (3);
  \end{tikzpicture}
  \caption{\textbf{Assumed causal structure for the synthetic experiments} (see Section \ref{sec:synthetic_experiments}). The precise corresponding generative random variables are described in Appendix \ref{app:synthetic}.}\label{fig:graph_experiment_1}
\vspace{-3mm}
\end{wrapfigure}

\subsection{Synthetic experiments}\label{sec:synthetic_experiments}
We simulate different datasets following to the causal graphical structure in Figure \ref{fig:graph_experiment_1}. The datasets are composed of four sets of observed continuous variables (i) the prediction target $\mathbf{Y}$, (ii) the variable(s) we want to be counterfactually invariant in $\mathbf{A}$, (iii) observed covariates that mediate effects from $\mathbf{A}$ on $\mathbf{Y}$, and (iv) observed confounding variables $\mathbf{Z}$.
The goal is to learn a predictor $\hat{\mathbf{Y}}$ that is counterfactually invariant in $\mathbf{A}$ with respect to $\mathbf{W} \coloneqq \mathbf{A} \cup \mathbf{X} \cup \mathbf{Z}$. We consider various artificially generated datasets for this example, which mainly differ in the dimension of the observed variables and their correlations. A detailed description of each dataset is deferred to Appendix \ref{app:synthetic}.

For all synthetic experiments, we train the model using fully connected networks (MLPs). We use the \textbf{MSE} loss $\mathcal{L}_{\textsc{mse}}( \hat{\mathbf{Y}})$ as the predictive loss $\mathcal{L}$ in eq.~\eqref{eq:total_loss} for continuous outcomes $\mathbf{Y}$. We generate $4$k samples from the observational distribution in each setting and use a 80:20 train/test split.
All metrics reported are on the test set.
We perform hyper-parameter tuning for MLP hyperparameters based on a random strategy (see Appendix \ref{app:synthetic} for details).
The $\hscic{\hat{\mathbf{Y}},\mathbf{A}}{\mathbf{Z}}$ term is computed as in eq.~\eqref{eq:estimate_hscic} using a Gaussian kernel with amplitude $1.0$ and length scale $0.1$.
The regularization parameter $\lambda$ for the ridge regression coefficients in eq.~\eqref{eq:cond_measure_estimand} is set to $\lambda = 0.01$.
To estimate \textbf{VCF} we set $d=1000$ and $k=500$.

\paragraph{Model performance with the use of the HSCIC.} We first perform a set of experiments to study the effect of the HSCIC, and to highlight the trade-off between accuracy and counterfactual invariance. For this set of experiments, we generate a dataset as described in Appendix \ref{appendix:synthetic_experiment}. Figure \ref{fig:simulated_exp} (left) shows the \textbf{VCF} for different regularization parameters $\gamma$, as a function of the values of $\mathbf{Z}$. We observe that increasing values of $\gamma$ lead to a consistent decrease of \textbf{VCF}. Figure \ref{fig:simulated_exp} (center) shows the values attained by the \textsc{HSCIC}{} and \textbf{MSE} for increasing parameter $\gamma$, demonstrating the expected trade-off in raw predictive performance and enforcing counterfactual invariance.
Finally, Figure \ref{fig:simulated_exp} (right) highlights the usefulness of \textsc{HSCIC}{} as a measure of counterfactual invariance, being in strong agreement with \textbf{VCF}.

\begin{figure}
  \centering
    \includegraphics[width=0.295\linewidth]{scatterplot_VCF.pdf}
    \includegraphics[width=0.34\linewidth]{plot_mse_hscic.pdf}
    \includegraphics[width=0.34\linewidth]{plot_vcf_hscic.pdf}
  \caption{%
  \textbf{(Left) variance of the counterfactual distributions} for $d=100$ datapoints with lines representing $3^{\mathrm{rd}}$-order polynomial regression and shaded areas being 95\% confidence intervals.
  \textbf{(Center) trade-off between accuracy and counterfactual invariance.} We observe that the \textsc{HSCIC}{} decreases, as the \textbf{MSE} increases. Vertical bars denote standard errors over $15$ different random seeds. \textbf{(Right) Correspondence between the \textsc{HSCIC}{} and the \textbf{VCF}}, for increasing regularizer $\gamma$. Again, vertical bars denote standard errors over $15$ different random seeds.}
  \label{fig:simulated_exp}
\end{figure}

\paragraph{Comparison with baselines.}
We perform a second set of experiments to compare our method against baselines. For this set of experiments we consider two settings, which we refer to as Scenario 1 and Scenario 2. These two settings differ in how the conditioning set $\textbf{Z}$ affects the outcome $\textbf{Y}$, with
Scenario 1 exhibiting higher correlation of $\textbf{Z}$ with both the mediator $\textbf{X}$ and the outcome $\textbf{Y}$. We refer the reader to Appendix \ref{appendix:baselines} for a precise description of Scenario 1 and Scenario 2. We remark that counterfactually invariant training has not received much attention yet. This limits the choice of baselines for experimental comparison. The main work on this topic to date does not provide a suitable benchmark for experimental comparison \citep{DBLP:conf/nips/VeitchDYE21}. In fact, the method by \citet{DBLP:conf/nips/VeitchDYE21} applies to a weaker notion of counterfactual invariance, and it can be used only on specific graphical structures which differ from our set-up. However, since counterfactual fairness is a special case of counterfactual invariance, we can compare to two methods proposed by \citet{DBLP:conf/nips/KusnerLRS17} (in applicable settings). We compare to the \textit{Level 1} (only use non-descendants of $\mathbf{A}$ as inputs to $\hat{\mathbf{Y}}$) and the \textit{Level 3} (assume an additive noise model and in addition to non-descendants, only use the residuals of descendants of $\mathbf{A}$ after regression on $\mathbf{A}$ as inputs to $\hat{\mathbf{Y}}$) approaches of \citet{DBLP:conf/nips/KusnerLRS17}.
We refer to these two baselines as CF1 and CF2 respectively. Note that CF1 satisfies counterfactual invariance perfectly by construction. We compare our method against CF1 and CF2 on both Scenario 1 and Scenario 2. We summarize the results in Table \ref{table:baseline_comparison}.
For a suitable choice of $\gamma$, our method outperforms the baseline CF2 in both \textbf{MSE} and \textbf{VCF} simultaneously.

\paragraph{Multi-dimensional variables.}
We perform a third set of experiments to assess \textsc{HSCIC}{}'s performance in higher dimensions. We consider simulated datasets as described in Appendix \ref{multid_setting}.
In this set of experiments, we independently increase the dimension of $\textbf{Z}$ and $\textbf{A}$ in two different simulations, leaving the rest of the variables unchanged. \Cref{table:multidimensional_results_A,table:multidimensional_results_Z} present the results for different regularization coefficients $\gamma$ and different dimensions of $\textbf{A}$ and $\textbf{Z}$. These results demonstrate that \textsc{HSCIC}{} can handle multi-dimensional variables while maintaining performance, as counterfactual invariance is approached when $\gamma$ increases.

\begin{table}
  \centering
  \caption{\textbf{Performance of the \textsc{HSCIC}{} against baselines} CF1 and CF2 on two synthetic datasets (see Appendix \ref{appendix:baselines}). Notably, for $\gamma=5$ in Scenario 1 and $\gamma=13$ in Scenario 2 we outperform CF2 in \textbf{MSE} and \textbf{VCF} simultaneously.}
  \label{table:baseline_comparison}
  \begin{tabularx}{\columnwidth}{l*6{Z}}
     \toprule
     & \multicolumn{2}{c}{Scenario 1}
     & \multicolumn{2}{c}{Scenario 2} \\
     \cmidrule(lr){2-3} \cmidrule(l){4-5}
     & \textbf{MSE} {\color{gray}$\times 10^3$} &$\textbf{VCF}${\color{gray}$\times 10^3$} & \textbf{MSE}{\color{gray}$\times 10^3$} & $\textbf{VCF}${\color{gray}$\times 10^3$} \\
   \midrule
   $\gamma=0$  & $0.36 \pm 0.5 $ &  $20.33 \pm 1$ & $2 \pm 0.1$ &  $239 \pm 0.4$ \\
   $\gamma=5$   & $16.97 \pm 2$ & $6.84 \pm 1$ & $43 \pm 8$ &  $230 \pm 6$ \\
   $\gamma=10$ &$19.80\pm 2$ & $5.89 \pm 0.6$ & $50 \pm 8$ & $190 \pm 3$ \\
   $\gamma=13$   &$22 \pm 1$ &$5.78 \pm 0.5$ &  $133 \pm 6$ & $157 \pm 7$ \\
   \midrule
   \textrm{CF1}  & $24.44 \pm 3$  & $0$ & $218 \pm 10$  & $0$ \\
   \textrm{CF2}   & $19.50 \pm 2$ & $7.50 \pm 1$ & $137 \pm  4$ & $167 \pm 2$ \\
   \bottomrule

   \end{tabularx}
   \vspace{10pt}
\end{table}

\begin{table}
  \centering
  \caption{\textbf{Results of the \textbf{MSE} ($\times10^5$ for readability), \textsc{HSCIC}{}, \textbf{VCF} for increasing dimension of $\textbf{A}$}, on a synthetic dataset as in Appendix \ref{multid_setting}. In this case, the dimension of the other variables are $\text{dimZ}=1, \text{dimX}=1,\text{dimY}=1$.}
  \label{table:multidimensional_results_A}
  \scriptsize
  \begin{tabularx}{\columnwidth}{@{}Yrrrrrrrrrrrr@{}}
   \toprule
   & \multicolumn{3}{c}{dimA=2}
   & \multicolumn{3}{c}{dimA=5}
   & \multicolumn{3}{c}{dimA=10}
   & \multicolumn{3}{c}{dimA=20} \\
   \cmidrule(lr){2-4} \cmidrule(lr){5-7} \cmidrule(lr){8-10} \cmidrule(l){11-13}
   $\gamma$ & \textbf{MSE} & \textbf{HSCIC} & \textbf{VCF} & \textbf{MSE} & \textbf{HSCIC} & \textbf{VCF} & \textbf{MSE} & \textbf{HSCIC} & \textbf{VCF} & \textbf{MSE} & \textbf{HSCIC} & \textbf{VCF} \\
   \midrule
   $0$   & $35$ & $0.0230$ & $0.025$ & $1.53$ & $0.0177$ & $0.013$ & $10$  & $0.0150$ &  $0.004$& $1.45$& $0.0126$ & $0.0010$ \\
   $\tfrac{1}{2}$  & $86$ & $0.0214$ & $0.022$  & $1580$ & $0.0157$ & $0.010$ & $40$ & $0.0135$& $0.002$ &$33$ & $0.0115$ & $0.0004$ \\
   $1$   & $210$ & $0.0200$ &  $0.220$ & $1450$ & $0.0144$  & $0.008$  & $100$ & $0.0127$ & $0.001$ & $60$& $0.0111$ & $0.0003$ \\
   \bottomrule
   \end{tabularx}%
   \vspace{10pt}
\end{table}

\begin{table}[t]
  \centering
  \caption{\textbf{Results of \textbf{MSE} ($\times 10^5$ for readability), \textsc{HSCIC}{}, \textbf{VCF} for increasing dimension of $\textbf{Z}$}, on a synthetic dataset as in Appendix \ref{multid_setting}. The dimension of the other variables are
  $\text{dimA}=1, \text{dimX}=1,\text{dimY}=1$.}
  \label{table:multidimensional_results_Z}
  \begin{tabularx}{\columnwidth}{@{}Yrrrrrrrrr@{}}
   \toprule
   & \multicolumn{3}{c}{dimZ=2}
   & \multicolumn{3}{c}{dimZ=3}
   & \multicolumn{3}{c}{dimZ=5} \\
   \cmidrule(lr){2-4} \cmidrule(lr){5-7} \cmidrule(l){8-10}
   $\gamma$ & \textbf{MSE} & \textbf{HSCIC} & \textbf{VCF} & \textbf{MSE} & \textbf{HSCIC} & \textbf{VCF} & \textbf{MSE} & \textbf{HSCIC} & \textbf{VCF} \\
   \midrule
   $0$   & $20$ & $0.002024$ & $6.15 \times 10^{-5}$ & $11$ & $0.03407$ & $0.023$ & $53.9$ & $0.033180$ & $0.0062$ \\
   $\tfrac{1}{2}$    & $20$ & $0.002023$ & $1.29 \times 10^{-5}$ & $49$ &$0.03406$ & $0.016$ & $21.9$ & $0.033178$ & $0.0061$ \\
   $1$    &  $20$ & $0.002019$ & $4.40 \times 10^{-6}$ & $57$ & $0.03403$& $0.009$ &$18$ & $0.033176$ & $0.0039$\\
   \bottomrule
   \end{tabularx}%
   \vspace{10pt}
\end{table}

\subsection{High-dimensional image experiments}\label{sec:imagedata}

\begin{figure}
  \centering
  \includegraphics[width=0.32\linewidth]{plot_mse_imagedataset.pdf}
  \includegraphics[width=0.32\linewidth]{plot_hscic_imagedataset.pdf}
  \includegraphics[width=0.32\linewidth]{plot_vcf_imagedataset.pdf}
  \caption{\textbf{Results of \textbf{MSE}, \textsc{HSCIC}{} operator and \textbf{VCF} for the dSprites image dataset experiment.} The \textsc{HSCIC}{} operator decreases steadily with higher values of $\gamma$. Similarly, a necessary increase of \textbf{MSE} can be observed. For both $\gamma=0.5$ and $\gamma=1$ an overall decrease of \textbf{VCF} is observed compared to the not-regularized setting. However, probably due to the high dimensional setting, the analysed metrics do not relevantly vary for the two regularized settings.
  }
  \label{fig:image_exp}
\end{figure}
We perform a set of experiments on a high-dimensional real-world dataset. We consider the image classification task as described in \ref{sec:counterfactual_invariance}. In this set of experiments, we use the dSprites dataset with a causal graph as in Figure \ref{fig:graphs}(c).
The full structural equations are provided in Appendix \ref{app:sec:imageresults}, where we assume a causal graph over the determining factors of the image, and essentially look up the corresponding image in the simulated dataset.
This experiment is particularly challenging due to the mixed categorical and continuous variables in $\mathbf{C}$ (\texttt{shape}, \texttt{y-pos}) and $\mathbf{X}$ (\texttt{color}, \texttt{orientation}), continuous $\mathbf{A}$ (\texttt{x-pos}).
All variables except for scale are assumed to be observed, and all variables jointly with the actual image determine the outcome $\textbf{Y}$.
Our goal is to learn a predictor $\hat{\mathbf{Y}}$ that is counterfactually invariant in the x-position with respect to all other observed variables.
In the chosen causal structure, $\{\texttt{shape}, \texttt{y-pos}\}\in \mathbf{C}$ block all non-directed paths from $\mathbf{A} \cup \mathbf{C} \cup \mathbf{U}$ to $\mathbf{Y}$.
Hence, we seek to achieve $\hat{\mathbf{Y}} \perp \!\!\! \perp \texttt{x-pos} \mid \{\texttt{shape}, \texttt{y-pos}\}$ via the \textsc{HSCIC}{} operator.
To accommodate the mixed input types, $\hat{\mathbf{Y}}$ puts an MLP on top of features extracted from the images via convolutional layers concatenated with features extracted from the remaining inputs via an MLP.

Figure \ref{fig:image_exp} demonstrates that \textsc{HSCIC}{} achieves improved \textbf{VCF} as $\gamma$ increases up to a certain point while affecting \textbf{MSE}, a necessary trade-off.
Adjusting our method to achieve higher stability is an interesting direction for future work.

\subsection{Fairness with continuous protected attributes}

In this experiment we apply our method to the popular UCI Adult dataset \citep{uci_dataset}, described in Section \ref{sec:counterfactual_invariance}. Following previous related work \citep{chiappa2021fairness}, we assume the causal struture depicted in the Appendix \ref{app:sec:adult}. Furthermore, we assume additional latent variables as in Figure \ref{fig:causal_graph_uci}.
Our goal is to learn a classifier for whether

\begin{figure}
\begin{minipage}[t]{.28\textwidth}
\vspace{0pt}
  \includegraphics[width=\linewidth]{kdeplot_VCF_uci.png}
\vspace{0pt}
\end{minipage}
\hspace{-7mm}
\hfill
\begin{minipage}{.4\textwidth}
\hspace{-3.9cm}
  \vspace{-2cm}
  \scriptsize
  \begin{tabular}{@{} l *{3}{c} @{}}
   \toprule
   & \textbf{MSE} & \textbf{HSCIC}& \textbf{VCF}\\
   \midrule
   $\gamma=0$   & $(81.725 \pm  0.0481) \%$ & $0.01067 \pm 3.4 \times 10^{-5}$ & $0.06685 \pm 0.006$\\
   $\gamma=1$   & $(81.547 \pm 0.0745) \%$ & $0.00739 \pm 6.8 \times 10^{-5}$&$0.05156 \pm 0.003$\\
   \bottomrule
  \end{tabular}
\vspace{-1cm}
\end{minipage}
\vspace{-2mm}
  \caption{\textbf{(Left) Probability density mass of the \textbf{VCF} distribution}, for different choices of $\gamma$. We observer less variance and more mass near zero for regularization parameter  $\gamma=1$. \textbf{(Right) Results of accuracy, \textsc{HSCIC}{} and \textbf{VCF}}, showing a considerable improvement of \textsc{HSCIC}{} and \textbf{VCF} for $\gamma=1$.}
  \label{fig:fairness_2}
  \vspace{-3mm}
\end{figure}

 the income of an individual is above a certain threshold, which is counterfactually invariant in both age and gender with respect to a collection of additional demographic attributes.
Here, we use an MLP with binary cross-entropy loss for $\hat{\mathbf{Y}}$.
Since this experiment is based on real data, the true counterfactual distribution cannot be known.
For an evaluation of our method, we thus again follow \citet{chiappa2021fairness} and estimate a possible true underlying structure by inferring the posterior distribution over the unobserved variables using variational autoencoders \citep{KingmaW13}.
Even though this can only provide approximate variances of counterfactual distributions, Figure \ref{fig:fairness_2} (left) shows that \textsc{HSCIC}{} still achieves more counterfactually fair outcome distributions (more mass near zero) than an unconstrained classifier ($\gamma = 0$).
Figure \ref{fig:fairness_2} (right) highlights once more that the \textsc{HSCIC}{} operator reliably decreases also in classification setting, again trading off accuracy.

\section{Discussion and future work}\label{sec:discussion}

We study the problem of learning predictors $\hat{\mathbf{Y}}$ that are ``counterfactually invariant'' in changes of certain covariates.
We put forth a formal definition of counterfactual invariance and describe how it generalizes existing notions.
Next, we provided a novel graphical criterion to characterize counterfactual invariance and reduce it to a conditional independence statement. Our method does not require identifiability of the counterfactual distribution or exclude the possibility of unobserved confounders.
Finally, we propose an efficiently estimable, model-agnostic regularization term to enforce this conditional independence based on kernel mean embeddings of conditional measures, which works for mixed continuous/categorical, multi-dimensional variables. We demonstrate the efficacy of our method in regression and classification tasks involving controlled synthetic scenarios, high-dimensional images, and in a fairness application, where it outperforms existing baselines.

The main limitation of our work, as for all studies in this domain, is the assumption that the causal graph is known. Another limitation is that our methodology is applicable when a graphical criterion is satisfied, i.e., the required variables are observed. From an ethics perspective, the increased robustness of counterfactually invariant (or societal benefits of counterfactually fair) predictors are certainly desirable.
However, this presupposes that all the often untestable assumptions are valid. Overall, causal methodology should not be applied lightly, especially in high-stakes and consequential decisions. A critical analysis of the broader context or systemic factors may hold more promise for societal benefits, than a well-crafted algorithmic predictor.

In light of the limitations mentioned above, an interesting and important direction for future work is to assess the sensitivity of our method to misspecifications of the causal graph or insufficient knowledge of the required blocking set. Finally, we believe our graphical criterion and KMEs-based regularization can also be useful for causal representation learning, where one aims to isolate causally relevant, autonomous factors underlying the data generating process of high-dimensional data.

\subsubsection*{Acknowledgement}
This work was supported by the German Federal Ministry of Education and Research (BMBF): T\"ubingen AI Center, FKZ: 01IS18039B.

\bibliographystyle{apalike}

\section{Introduction and related work}\label{sec:introduction}

Invariance, or equivariance to certain transformations of data, has proven essential in numerous applications of machine learning (ML). Invariance leads to better generalization capabilities of ML models~\citep{arjovsky2019invariant, chen2020simple, bloem2020probabilistic}.
For instance, in image recognition tasks, predictions ought to remain unchanged under scaling, translation, or rotation of the input image. Data augmentation is one of the earliest heuristics developed to promote this kind of invariance that has become indispensable in training successful models, like deep neural networks (DNNs)~\citep{shorten2019survey, NEURIPS2020_44feb009}. State-of-the-art architectures can handle invariance and equivariance requirements by design. Well-known examples include convolutional neural networks (CNNs) for translation invariance of image data~\citep{NIPS2012_c399862d}, group equivariant CNNs for other group transformations~\citep{cohen2016group}, recurrent neural networks (RNNs) and transformers for sequential data~\citep{vaswani2017attention}, DeepSet~\citep{NIPS2017_f22e4747} for sets, and graph neural networks (GNNs) for different types of geometric structures~\citep{battaglia2018relational}.

Many real-world applications in modern ML, however, call for an arguably stronger notion of invariance, called \emph{counterfactual invariance}. This is the case for image classification, algorithmic fairness~\citep{NIPS2016_9d268236, mitchell2021algorithmic}, robustness \citep{buhlmann2020invariance}, and out-of-distribution  generalization \citep{lu2021invariant}. These applications require predictors to exhibit invariance with respect to hypothetical manipulations of the data generating process (DGP)~\citep{peters2016causal,heinze2018invariant,rojas2018invariant,arjovsky2019invariant,buhlmann2020invariance}. In image classification, for instance, we want to build a model that ``would have made the same prediction, if the object position had been different with everything else being equal''. This model ought to be invariant to the DGP of some attributes, such as an object position and light condition. Similarly, in algorithmic fairness \citet{kilbertus2017avoiding,DBLP:conf/nips/KusnerLRS17} introduce notions of interventional and counterfactual fairness. These notions are based on the causal relationships between the observed variables of a model, and they require invariance to the DGP of protected attributes.

Inspired by problems in natural language processing (NLP), \citet{DBLP:conf/nips/VeitchDYE21} study a notion of counterfactual invariance. In their work, \citet{DBLP:conf/nips/VeitchDYE21} provide a method to enforce their notion of counterfactual invariance, based on the maximum mean discrepancy (MMD). However, this method promotes only a \emph{necessary}, but not sufficient observable imprint or proxy of counterfactual invariance during training. Furthermore, their work can only handle discrete random variables, and it applies to specific causal graphs. Motivated by this line of research, we study the following question:
\begin{center}
\emph{How do we learn predictors that are invariant\\under counterfactual changes of certain covariates?}
\end{center}

We propose a new definition of counterfactual invariance that generalizes the notion studied by \citet{DBLP:conf/nips/VeitchDYE21}. Furthermore, we propose a new method that directly enforces counterfactual invariance during the learning process. Unlike existing approaches, the proposed method is \emph{sufficient} for counterfactual invariance, and it works well for both categorical and continuous covariates and outcomes, as well as in multivariate settings

\paragraph{Technical overview.} Our notion of counterfactual invariance inevitably relies on strong assumptions about the DGP, as counterfactuals by their very nature cannot be observed. This is a major challenge, since counterfactual invariance cannot be checked, unless strong prior knowledge of the DGP is available. We circumvent this problem by deriving a connection between counterfactual invariance and conditional independence, under certain observability assumptions in the causal graph. The resulting conditional independence statements can then be enforced via a regularization term based on the flexible Hilbert-Schmidt Conditional Independence Criterion (\textsc{HSCIC}{}) \citep{Park20:CME}.

\section{Preliminaries}\label{sec:preliminaries}

\paragraph{Counterfactual invariance.}
We introduce structural causal models as in \citet{pearlj}.

\begin{definition}[Structural causal model (SCM)]
\label{def:SCM}
A structural causal model is a tuple $(\mathbf{U}, \mathbf{V}, F, \mathbb{P}_{\mathbf{U}})$ such that
\begin{enumerate}
\item $\mathbf{U}$ is a set of background variables that are exogenous to the model;
\item $\mathbf{V}$ is a set of observable (endogenous) variables; \item $F=\{f_{V}\}_{V\in\mathbf{V}}$ is a set of functions from (the domains of) $\mathsf{pa}(V) \cup U_V$ to (the domain of) $V$, where $U_V \subset \mathbf{U}$ and $\mathsf{pa}(V) \subseteq \mathbf{V} \setminus \{V\}$ such that $V = f_V(\mathsf{pa}(V), U_V)$; (iv) $\mathbb{P}_{\mathbf{U}}$ is a probability distribution over the domain of $\mathbf{U}$.
Further, the subsets $\mathsf{pa}(V) \subseteq \mathbf{V} \setminus \{V\}$ are chosen such that the graph $\mathcal{G}$ over $\mathbf{V}$ where the edge $V' \to V$ is in $\mathcal{G}$ if and only if $V' \in \mathsf{pa}(V)$ is a directed acyclic graph (DAG).
\end{enumerate}
\end{definition}
We always denote with $\mathbf{Y} \subset \mathbf{V}$ the outcome (or prediction target), and with $\hat{\mathbf{Y}}$ a predictor for that target. The predictor $\hat{\mathbf{Y}}$ is not strictly part of the SCM, because we get to tune it. However, since it takes inputs from $\mathbf{V}$, we often treat it as an observed variable in the SCM.
Each SCM implies a unique \emph{observational distribution} over $\mathbf{V}$ \citep{pearlj}.
However, SCMs allow us to also define interventional distributions. Given a variable $A \in \mathbf{V}$, an intervention $A \gets a$ amounts to replacing $f_A$ in $F$ with the constant function $A = a$. This yields a new SCM, which induces the \emph{interventional distribution}.
Similarly, we can intervene on multiple variables $\mathbf{V} \supseteq \mathbf{A} \gets \mathbf{a}$.
We then write $\mathbf{Y}^*_{\mathbf{a}}$ for the outcome in the intervened SCM, also called potential outcome. Note that the interventional distribution $\mathbb{P}_{\mathbf{Y}^*_{\mathbf{a}}}(\mathbf{y})$ differs in general from the conditional distribution $\mathbb{P}_{\mathbf{Y}\mid \mathbf{A}}(\mathbf{y} \mid \mathbf{a})$, due to possible unobserved confounding effects.
We can also condition on a set of variables $\mathbf{W} \subseteq \mathbf{V}$ in the (observational distribution of the) original SCM before performing an intervention, which we denote by $\mathbb{P}_{\mathbf{Y}^*_{\mathbf{a}} \mid \mathbf{W}} (\mathbf{y} \mid \mathbf{w})$.
This is a counterfactual distribution: ``Given that we have observed $\mathbf{W} = \mathbf{w}$, what would $\mathbf{Y}$ have been had we set $\mathbf{A} \gets \mathbf{a}$, instead of the values $\mathbf{A}$ had actually taken?''.
Note that the sets $\mathbf{A}$ and $\mathbf{W}$ need not be disjoint. We can now define counterfactual invariance.
\begin{definition}[Counterfactual invariance]\label{def:counterfactual_invariance}
Let $\mathbf{A}$, $\mathbf{W}$ be (not necessarily disjoint) sets of nodes in a given SCM. A predictor $\hat{\mathbf{Y}}$ is counterfactually invariant in $\mathbf{A}$ with respect to $\mathbf{W}$ if $\mathbb{P}_{\hat{\mathbf{Y}}^*_{\mathbf{a}} \mid \mathbf{W}}(\mathbf{y}\mid \mathbf{w}) = \mathbb{P}_{\hat{\mathbf{Y}}^*_{\mathbf{a}'} \mid \mathbf{W}}(\mathbf{y}\mid \mathbf{w})$ almost surely, for all $\mathbf{a}, \mathbf{a}'$ in the domain of $\mathbf{A}$, and all $\mathbf{w}$ in the domain of $\mathbf{W}$.\footnote{With a mild abuse of notation, if $\mathbf{W} = \emptyset$ then the requirement of conditional counterfactual invariance becomes $\mathbb{P}_{\hat{\mathbf{Y}}^*_{\mathbf{a}}}(\mathbf{y}) = \mathbb{P}_{\hat{\mathbf{Y}}^*_{\mathbf{a}'}}(\mathbf{y})$ almost surely, for all $\mathbf{a}, \mathbf{a}'$ in the domain of $\mathbf{A}$.}
\end{definition}
A predictor that is counterfactually invariant can thus be viewed as robust to changes of $\mathbf{A}$, in the sense that the (conditional) post-interventional distribution of $\hat{\mathbf{Y}}$ does not change for different values of the intervention.
In other words, the classification would not change had we gone back and changed the position of the object.
This definition is stronger than previous definitions of counterfactual invariance (see Definition 1.1 by \citet{DBLP:conf/nips/VeitchDYE21}). In fact, the invariance notion by \citet{DBLP:conf/nips/VeitchDYE21} only requires $\mathbb{P}_{\hat{\mathbf{Y}}^*_{\mathbf{a}}}(\mathbf{y}) = \mathbb{P}_{\hat{\mathbf{Y}}^*_{\mathbf{a}'}}(\mathbf{y})$ almost surely, for all $\mathbf{a}, \mathbf{a}'$ in the domain of $\mathbf{A}$.

\paragraph{Kernel mean embeddings and conditional measures.}
Our regularization term heavily relies on kernel mean embeddings (KMEs).
We now highlight the main concepts pertaining KMEs and refer the reader to \citet{smola2007hilbert,learning_kernels_book,prob_stats_book,Muandet17:KME} for more details.
Fix a measurable space $\mathscr{Y}$ with respect to a $\sigma$-algebra $\mathcal{F}_{\mathscr{Y}}$, and consider a probability measure $\mathbb{P}$ on the space $(\mathscr{Y}, \mathcal{F}_{\mathscr{Y}})$. Let $\mathcal{H}$ be a reproducing kernel Hilbert space (RKHS) with a bounded kernel $k_\mathbf{Y} \colon \mathscr{Y} \times \mathscr{Y} \rightarrow \mathbb{R}$, i.e., $k_\mathbf{Y}$ is such that $\sup_{\mathbf{y} \in \mathscr{Y}}k(\mathbf{y},\mathbf{y}) < +\infty$. The kernel mean embedding $\mu_{\mathbb{P}}$ of $\mathbb{P}$ is defined as the expected value of the function $k(\ \cdot \ ,\mathbf{y})$ with respect to $\mathbf{y}$, i.e., $\mu_{\mathbb{P}} \coloneqq \expect{}{k(\ \cdot \  ,\mathbf{y})}$.

The definition of KMEs can be extended to conditional distributions \citep{DBLP:journals/jmlr/FukumizuSG13,DBLP:conf/icml/GrunewalderLGBPP12,DBLP:conf/icml/SongHSF09,DBLP:journals/spm/SongFG13}.
Consider two random variables $\mathbf{Y}$, $\mathbf{Z}$, and denote with $(\Omega_\mathbf{Y}, \mathcal{F}_{\mathbf{Y}})$ and $(\Omega_\mathbf{Z}, \mathcal{F}_{\mathbf{Z}})$ the respective measurable spaces.
These random variables induce a probability measure $\mathbb{P}_{\mathbf{Y}, \mathbf{Z}}$ in the product space $\Omega_\mathbf{Y} \times \Omega_\mathbf{Z}$.
Let $\mathcal{H}_\mathbf{Y}$ be a RKHS with a bounded kernel $k_\mathbf{Y} (\cdot, \cdot )$ on $\Omega_\mathbf{Y}$.
We define the KME of a conditional distribution $\mathbb{P}_{\mathbf{Y}\mid \mathbf{Z}}(\cdot \mid \mathbf{z})$ via $\mu_{\mathbf{Y}\mid \mathbf{Z} = \mathbf{z}} \coloneqq \expect{}{k_\mathbf{Y}(\ \cdot \ ,\mathbf{y})\mid \mathbf{Z} = \mathbf{z}}$.
Here, the expected value is taken over $\mathbf{y}$.
KMEs of conditional measures can be estimated from samples. To illustrate this, consider i.i.d.\ samples $(\mathbf{y}_1, \mathbf{z}_1), \dots, (\mathbf{y}_n, \mathbf{z}_n)$.
Denote with $\hat{K}_\mathbf{Y}$ the kernel matrix with entries $[\hat{K}_\mathbf{Y}]_{i,j} \coloneqq k_\mathbf{Y}(\mathbf{y}_i, \mathbf{y}_j)$.
Furthermore,  let $k_\mathbf{Z}$ be a bounded kernel on $\Omega_\mathbf{Z}$.
Then, $\mu_{\mathbf{Y}\mid \mathbf{Z} = \mathbf{z}}$ can be estimated as
\begin{equation}\label{eq:cond_measure_estimand}
 \hat{\mu}_{\mathbf{Y}\mid \mathbf{Z} = \mathbf{z}} \coloneqq  \sum_{i = 1}^n \hat{w}_{\mathbf{Y}\mid \mathbf{Z}}^{(i)}(\mathbf{z})k_\mathbf{Y}(\cdot, \mathbf{y}_i)\; \text{ with }\; \hat{w}_{\mathbf{Y}\mid \mathbf{Z}}(\cdot) \coloneqq ( \hat{K}_\mathbf{Z} - n \lambda I )^{-1} \left [k_\mathbf{Z}(\cdot,\mathbf{z}_1), \cdots , k_\mathbf{Z}(\cdot,\mathbf{z}_n)\right ]^T
\end{equation}
where, $I$ is the identity matrix and $\lambda$ is a regularization parameter.
Here, $\hat{w}_{\mathbf{X}\mid \mathbf{A}}^{(i)}(\cdot)$, the $i$-th entry of $\hat{w}_{\mathbf{X}\mid \mathbf{A}}(\cdot)$, are the coefficients of kernel ridge regression \citep{DBLP:conf/icml/GrunewalderLGBPP12}.

\section{Counterfactually invariant predictors}\label{sec:ci}
\paragraph{Counterfactual invariance from conditional independence.}
In this section, we establish a simple graphical criterion to express the desired counterfactual invariance as a conditional independence in the observational distribution, rendering it estimable from i.i.d.\ data.
We first repeat the notion of \emph{blocked paths} \citep{pearlj}.
\begin{restatable}{definition}{blockedpath}
\label{def:blocked_path}
Consider a path $\pi$ of causal graph $\mathcal{G}$. A set of nodes $\mathbf{Z}$ blocks $\pi$, if $\pi$ contains a triple of consecutive nodes connected in one of the following ways: $N_i \rightarrow Z \rightarrow N_j $, $N_i \gets Z \rightarrow N_j $, or $N_i \rightarrow M \gets N_j $, with $N_i, N_j \notin \mathbf{Z}$, $Z \in \mathbf{Z}$, and neither $M$ nor any descendent of $M$ is in $\mathbf{Z}$.
\end{restatable}
Using \ref{def:blocked_path}, we give the following criterion for counterfactual invariance.
\begin{restatable}{theorem}{condcounterfactual}
\label{lemma:conditiona_counterfactual}
Let $\mathcal{G}$ be a causal diagram, and consider two (not necessarily disjoint) sets of nodes $\mathbf{A}$, $\mathbf{W}$.
Let $\mathbf{Z}$ be a set of nodes that blocks all non-directed paths from $\mathbf{A}\cup \mathbf{W}$ to $\mathbf{Y}$.
Then, for any SCM compatible with $\mathcal{G}$, any predictor $\hat{\mathbf{Y}}$ that satisfies $\hat{\mathbf{Y}}\perp \!\!\! \perp \mathbf{A}\mid \mathbf{Z}$ is counterfactually invariant in $\mathbf{A}$ with respect to $\mathbf{W}$.
\end{restatable}

\begin{figure}[t]
  \centering
\begin{tikzpicture}[node distance=7mm and 10mm, thick, main/.style = {draw, circle, minimum size=0.8cm}, >={triangle 45}]
    \node[main] (1) {$\scriptstyle \mathbf{A}$};
    \node[main] (2) [right =of 1] {$\scriptstyle \mathbf{X}$};
    \node[main] (3) [right =of 2] {$\scriptstyle \mathbf{Y}$};
    \draw[->] (1) -- (2);
    \draw[->] (2) -- (3);
    \draw[->] (1) to [out=320, in=220, looseness=1] (3);
    \draw[<->, dashed] (1) to [out=45, in=135, looseness=1] (2);
    \node (label) [above =of 1] {\textbf{(a)}};
  \end{tikzpicture}%
  \hfill
  \begin{tikzpicture}[node distance=7mm and 10mm, thick, main/.style = {draw, circle, minimum size=0.8cm}, >={triangle 45}]
    \node[main] (1) {$\scriptstyle \mathbf{A}$};
    \node[main] (2) [right =of 1] {$\scriptstyle \mathbf{X}$};
    \node[main] (3) [right =of 2] {$\scriptstyle \mathbf{Y}$};
    \node[main] (4) [above =of 2] {$\scriptstyle \mathbf{C}$};
    \draw[->] (1) -- (2);
    \draw[->] (2) -- (3);
    \draw[->] (4) -- (2);
    \draw[->] (4) -- (3);
    \draw[->] (1) to [out=320, in=220, looseness=1] (3);
    \node (label) [left =of 4] {\textbf{(b)}};
  \end{tikzpicture}%
  \hfill
    \begin{tikzpicture}[node distance=7mm and 10mm, thick, main/.style = {draw, circle, minimum size=0.8cm}, >={triangle 45}]
    \node[main] (1) {$\textbf{A}$};
    \node[main, dashed] (3) [right =of 1] {$\scriptstyle \textbf{U}$};
    \node[main] (5) [above =of 3] {$\scriptstyle \textbf{C}$};
    \node[main] (4) [right =of 3] {$\scriptstyle \textbf{Y}$};
    \node[main] (6) [above =of 4] {$\scriptstyle \textbf{X}$};
    \draw[->, dashed] (1) -- (3);
    \draw[->] (5) -- (1);
    \draw[->, dashed] (5) -- (3);
    \draw[->] (5) -- (4);
    \draw[->, dashed] (3) -- (4);
    \draw[->] (6) -- (4);
    \draw[->] (1) to [out=320, in=220, looseness=1] (4);
    \node (label) [left =of 5] {\textbf{(c)}};
  \end{tikzpicture}%
  \caption{%
    \textbf{(a) An example for \ref{lemma:conditiona_counterfactual}}, where $\mathbb{P}_{\mathbf{Y}^*_\mathbf{a}}(\mathbf{y})$ is unidentifiable, due to the confounding effect denoted by the dashed double arrow. However, by \ref{lemma:conditiona_counterfactual} any predictor $\hat{\mathbf{Y}}$ such that $\hat{\mathbf{Y}}\perp \!\!\! \perp \mathbf{A}\mid \mathbf{Z}$ with $\mathbf{Z} = \{\mathbf{X}\}$ is counterfactually invariant in $\mathbf{A}$. Hence, the criterion in \ref{lemma:conditiona_counterfactual} is weaker than identifiability.
    \textbf{(b) Assumed causal structure for the Adult dataset}, where $\mathbf{A}$ consists of the protected attributes \emph{gender} and \emph{age}. A predictor is learned for the income $\mathbf{Y}$ that is invariant in $\mathbf{A}$, which corresponds to counterfactual fairness (see \ref{def:counterfactual_fairness}).
    \textbf{(c) Causal structure for the constructed dSprites ground truth}, where $\mathbf{A} = \{\text{\small{Pos.X}}\}$, $\mathbf{U} = \{\text{\small{Scale}}\}$, $\mathbf{C} = \{\text{\small{Shape}}, \text{\small{Pos.Y}}\}$,
    $\mathbf{X} = \{\text{\small{Color}}, \text{\small{Orientation}}\}$, and $\mathbf{Y} = \{\text{\small{Outcome}}\}$.
   $\mathbf{U}$ is unobserved.
    We want a predictor for $\mathbf{Y}$, which is counterfactually invariant in $\mathbf{A}$ with respect to all other observed attributes (see \ref{sec:counterfactual_invariance}).
    }\label{fig:graphs}
\end{figure}

The proof is deferred to Appendix \ref{proof:conditiona_counterfactual}, where we derive the graphical criterion by studying a graphical model similar to the counterfactual graph \citep{DBLP:journals/jmlr/ShpitserP08}.
Our key observation is that the set $\mathbf{Z}$ as in \ref{lemma:conditiona_counterfactual} acts as a valid adjustment set for certain random variables in this joint pre- and post-interventional graphical model.
Using this observation, we derive \ref{theorem:identifiability_app} with similar techniques to covariate adjustment. However, our proof does \emph{not} achieve identification of the counterfactual distributions by simply applying the do-calculus \citep{pearlj}.
Crucially, \emph{the existence of a measurable set $\mathbf{Z}$ as in Theorem \ref{lemma:conditiona_counterfactual} does not imply the identifiability of counterfactual distributions $\mathbb{P}_{\mathbf{Y}^*_{\mathbf{a}}} (\mathbf{y})$} (see Figure \ref{fig:graphs}(a) for a counterexample). In particular,  Theorem \ref{lemma:conditiona_counterfactual} does not rule out hidden confounding in the model. Hence, our method is applicable when (certain parameters of) $\mathbb{P}_{\mathbf{Y}^*_{\mathbf{a}}} (\mathbf{y})$ or $\mathbb{P}_{\mathbf{Y}^*_{\mathbf{a}} \mid \mathbf{W}} (\mathbf{y})$ cannot be learned from observational data.

\paragraph{Counterfactually invariant predictors.}
We now propose our novel regularization term denoted by $\hscic{\hat{\mathbf{Y}},\mathbf{A}}{\mathbf{Z}}$ that promotes counterfactual invariance by encouraging the conditional independence $\hat{\mathbf{Y}} \perp \!\!\! \perp \mathbf{A} \mid \mathbf{Z}$ (see Theorem \ref{lemma:conditiona_counterfactual}). Given a loss function $\mathcal{L}(\hat{\mathbf{Y}})$ for the predictive performance of $\hat{\mathbf{Y}}$, we propose to minimize the following loss
\begin{equation}\label{eq:total_loss}
    \mathcal{L}_{\textsc{total}}(\hat{\mathbf{Y}}) = \mathcal{L}(\hat{\mathbf{Y}}) + \gamma \cdot\textsc{HSCIC}{}(\hat{\mathbf{Y}},\mathbf{A}\mid \mathbf{Z})\:,
\end{equation}
with a regularization parameter $\gamma \ge 0$.
We will now develop and justify the regularization term (see  Corollary~\ref{cor:hscic_independence}) for which $\textsc{HSCIC}(\hat{\mathbf{Y}},\mathbf{A}\mid \mathbf{Z}) = 0$ if and only if $\hat{\mathbf{Y}}\perp \!\!\! \perp \mathbf{A}\mid \mathbf{Z}$.

\paragraph{Hilbert-Schmidt Conditional Independence Criterion (\textsc{HSCIC}{}).}
Consider two random variables $\mathbf{Y}$ and $\mathbf{A}$, and denote with $(\Omega_\mathbf{Y}, \mathcal{F}_{\mathbf{Y}})$ and $(\Omega_\mathbf{A}, \mathcal{F}_{\mathbf{A}})$ the respective measurable spaces. Suppose that we are given two RKHSs $\mathcal{H}_{\mathbf{Y}}$, $\mathcal{H}_{\mathbf{A}}$ over the support of $\mathbf{Y}$ and $\mathbf{A}$ respectively. The tensor product space $\mathcal{H}_{\mathbf{Y}} \otimes \mathcal{H}_{\mathbf{A}}$ is defined as the space of functions of the form $(f \otimes g)(\mathbf{y}, \mathbf{a}) \coloneqq f(\mathbf{y}) g(\mathbf{a})$, for all $f \in \mathcal{H}_{\mathbf{Y}}$ and $g \in \mathcal{H}_{\mathbf{A}}$. The tensor product space yields a natural RKHS structure, with kernel $k$ defined by $k(\mathbf{y} \otimes \mathbf{a}, \mathbf{y}' \otimes \mathbf{a}') \coloneqq k_{\mathbf{Y}}(\mathbf{y}, \mathbf{y}') k_{\mathbf{A}}(\mathbf{a}, \mathbf{a}')$. We refer the reader, i.e., to \citet{DBLP:journals/jmlr/SzaboS17} for more details on tensor product spaces.
With this notation we define:
\begin{definition}[\textsc{HSCIC}{}]
\label{def:definition_hscic}
For (sets of) random variables $\mathbf{Y}$, $\mathbf{A}$, $\mathbf{Z}$ the \textsc{HSCIC}{} between $\mathbf{Y}$ and $\mathbf{A}$ given $\mathbf{Z}$ is defined as the real-valued random variable $\hscic{\mathbf{Y},\mathbf{A}}{\mathbf{Z}} = H_{\mathbf{Y},\mathbf{A}\mid \mathbf{Z}} \circ \mathbf{Z}$. Here, $H_{\mathbf{Y},\mathbf{A}\mid \mathbf{Z}}$ is a real-valued deterministic function, defined as
\[
H_{\mathbf{Y},\mathbf{A}\mid \mathbf{Z}}(\mathbf{z}) \coloneqq \left \| \mu_{\mathbf{Y},\mathbf{A}\mid \mathbf{Z} = \mathbf{z}} - \mu_{\mathbf{Y}\mid \mathbf{Z} = \mathbf{z}}\otimes \mu_{\mathbf{A}\mid \mathbf{Z} = \mathbf{z}} \right \|,
\]
with $\left \| \cdot \right \|$ the metric induced by the inner product of the tensor product space $\mathcal{H}_{\mathbf{X}}\otimes \mathcal{H}_{\mathbf{A}}$.
\end{definition}
The \textsc{HSCIC}{} has the following important property.
\begin{restatable}{theorem}{hscicind}
\label{thm:hscic_independence}
With the notation introduced above, suppose that the kernel $k$ of the tensor product space $\mathcal{H}_{\mathbf{X}}\otimes \mathcal{H}_{\mathbf{A}}$ is characteristic.\footnote{The tensor product kernel $k$ is characteristic if the mapping $\mathbb{P}_{\mathbf{Y}, \mathbf{A}} \mapsto \expect{}{k(\ \cdot \ , \mathbf{y}\otimes \mathbf{a})}$ is injective. Here, the expected value is taken over the pair $(\mathbf{y}, \mathbf{a})$.} Then, $\hscic{\mathbf{Y},\mathbf{A}}{\mathbf{Z}} = 0$ almost surely if and only if $\mathbf{Y}\perp \!\!\! \perp \mathbf{A} \mid \mathbf{Z}$.
\end{restatable}
We provide a proof in Appendix \ref{appendix:hscic_independence}.
The \textsc{HSCIC}{} as a criterion for conditional independence is not a novel concept \citep{Park20:CME}.
However, our definition differs from \citet{Park20:CME}, who define the \textsc{HSCIC}{} via KMEs based on the Bochner conditional expected value.
In contrast, we define the \textsc{HSCIC}{} using the embeddings of conditional measures $\mu_{\mathbf{Y}\mid \mathbf{Z} = \mathbf{z}}$. This allows us to bypass some technical assumptions required by \citet{Park20:CME}, while retaining generality over cross-covariance based methods \citep{DBLP:conf/nips/FukumizuGSS07}.
We refer readers to Appendix \ref{appendix:ind_bochner} and \ref{appendix:ind_cross_covariance} for a comparison with previous approaches.
Combining  Theorem \ref{lemma:conditiona_counterfactual} with Theorem \ref{thm:hscic_independence}, we can now use the \textsc{HSCIC}{} to promote counterfactual invariance.
\begin{corollary}
\label{cor:hscic_independence}
Consider an SCM with causal diagram $\mathcal{G}$ and fix two (not necessarily disjoint) sets of nodes $\mathbf{A}$, $\mathbf{W}$. Let $\mathbf{Z}$ be a set of nodes in $\mathcal{G}$ that blocks all non-directed paths from $\mathbf{A}\cup \mathbf{W}$ to $\mathbf{Y}$. Then, any predictor $\hat{\mathbf{Y}}$ that satisfies $\hscic{\hat{\mathbf{Y}},\mathbf{A}}{\mathbf{Z}} = 0$ almost surely is counterfactually invariant in $\mathbf{A}$ with respect to $\mathbf{W}$.
\end{corollary}

\paragraph{Estimating the \textsc{HSCIC}{} from samples.}
We show how to estimate $\hscic{\mathbf{Y},\mathbf{A}}{\mathbf{Z}}$ from finite samples using eq.~\eqref{eq:cond_measure_estimand}.
Given $n$ samples $\{(\mathbf{y}_i, \mathbf{a}_i, \mathbf{z}_i) \}_{i = 1}^n$, denote by $\hat{K}_{\mathbf{Y}}$ and $\hat{K}_{\mathbf{A}}$ the corresponding kernel matrices for $\mathbf{Y}$ and $\mathbf{A}$. We then estimate the $H_{\mathbf{Y},\mathbf{A}\mid \mathbf{X}}$ as \citep{Park20:CME}
\begin{align}
    \hat{H}^2_{\mathbf{Y},\mathbf{A}\mid \mathbf{Z}}(\cdot) & = \hat{w}^T_{\mathbf{Y},\mathbf{A}\mid \mathbf{Z}}(\cdot) \left ( \hat{K}_{\mathbf{Y}} \odot \hat{K}_{\mathbf{A}}\right ) \hat{w}_{\mathbf{Y},\mathbf{A}\mid \mathbf{Z}}(\cdot) - 2 \left ( \hat{w}^T_{\mathbf{Y} \mid \mathbf{Z}}(\cdot) \hat{K}_{\mathbf{Y}} \hat{w}_{\mathbf{Y},\mathbf{A} \mid \mathbf{X}}(\cdot)  \right )\label{eq:estimate_hscic} \\
    & \cdot \left ( \hat{w}^T_{\mathbf{A} \mid \mathbf{X}}(\cdot) \hat{K}_{\mathbf{A}} \hat{w}_{\mathbf{Y},\mathbf{A} \mid \mathbf{Z}}(\cdot)  \right ) + \left ( \hat{w}^T_{\mathbf{Y} \mid \mathbf{Z}}(\cdot) \hat{K}_{\mathbf{Y}} \hat{w}_{\mathbf{Y} \mid \mathbf{Z}}(\cdot)  \right ) \left ( \hat{w}^T_{\mathbf{A} \mid \mathbf{Z}}(\cdot) \hat{K}_{\mathbf{A}} \hat{w}_{\mathbf{A} \mid \mathbf{Z}}(\cdot)  \right )\;, \nonumber
\end{align}
where $\odot$ is element-wise multiplication.
As mentioned, the vectors $\hat{w}_{\mathbf{Y}\mid \mathbf{Z}}(\cdot)$, $\hat{w}_{\mathbf{A}\mid \mathbf{Z}}(\cdot)$, and $\hat{w}_{\mathbf{Y},\mathbf{A}\mid \mathbf{Z}}(\cdot)$ are found via kernel ridge regression.
Theorem 1 by \citet{DBLP:journals/focm/CaponnettoV07} shows that if the regularization parameters for the kernel ridge coefficients decay
to $0$ sufficiently quickly as $n \to \infty$, then $\hat{H}^2_{\mathbf{Y},\mathbf{A}\mid \mathbf{Z}}(\cdot)$ convergences to $H^2_{\mathbf{Y},\mathbf{A}\mid \mathbf{Z}}(\cdot)$ at a rate of $n^{-c}$, where $c$ is a parameter depending on the complexity of the embeddings and the effective dimensions of $\mathcal{H}_{\mathbf{Y}}$ and $\mathcal{H}_{\mathbf{A}}$.
In practice, computing the \textsc{HSCIC}{} approximation $\hat{H}^2_{\mathbf{Y},\mathbf{A}\mid \mathbf{Z}}(\cdot)$ may be time-consuming.
To speed it up, we can use random Fourier features to approximate the matrices $\hat{K}_{\mathbf{Y}}$ and $\hat{K}_{\mathbf{A}}$ as above \citep{DBLP:conf/nips/RahimiR07,DBLP:conf/icml/AvronKMMVZ17}.
We refer readers to Appendix \ref{appendix:random_fourier_features} for a discussion of the use of random Fourier features.

\paragraph{Measuring counterfactual invariance.}
Besides predictive performance, e.g., mean squared error (\textbf{MSE}) for regression, our key metric of interest is the level of counterfactual invariance achieved by the predictor $\hat{\mathbf{Y}}$.
We quantify the overall counterfactual variance as a single scalar
\begin{equation}\label{eq:vcf}
  \textbf{VCF}(\hat{\mathbf{Y}}) = \mathbb{E}_{\mathbf{w} \sim\mathbb{P}_{\mathbf{W}}}\left[ \var_{\mathbf{a}' \sim \mathbb{P}_{\mathbf{A}}}[\mathbb{E}_{\hat{\mathbf{Y}}^*_{\mathbf{a}'} \mid \mathbf{W}}(\hat{\mathbf{y}} \mid \mathbf{w})]\right]\:.
\end{equation}
For deterministic predictors (i.e., point estimators), which we use in all our experiments, the variance term is zero if and only if counterfactual invariance holds at $\mathbf{w}$ for (almost) all $\mathbf{a}'$ in the support of $\mathbb{P}_{\mathbf{A}}$ (i.e., the prediction remains constant).
Since the variance is non-negative, $\mathbf{VCF}$ (we often drop the argument) is zero if and only if counterfactual invariance holds almost surely.
However, the outer expectation of $\mathbf{VCF}$ may still be sensitive to outliers in its arguments, which is why it can be instructive to also analyze the variances for different conditioning values $\mathbf{Z}=\mathbf{z}$.
To estimate \textbf{VCF} in practice, we pick $d$ datapoints $(\mathbf{w}_i)_{i=1}^d$ from the observed data, and for each compute the counterfactual outcomes $\hat{\mathbf{Y}}^*_{\mathbf{a}'}\mid \mathbf{w}_i$ for $k$ different values of $\mathbf{a}'$.
The inner expectation is simply the predictor output.
We use empirical variances with $k$ examples for each of the $d$ chosen datapoints, and the empirical mean
of the $d$ variances for the outer expectation.
Finally, we note that by \ref{cor:hscic_independence}, \textsc{HSCIC}{} itself also serves as a useful measure of how well counterfactual invariance has been enforced.

\begin{wrapfigure}{R}{0.45\textwidth}
\vspace{-3mm}
\centering
\includegraphics[width=0.35\textwidth]{dsprites_pic.png}
\vspace{1mm}
\caption{\textbf{Examples of images from the dSprites dataset} \citep{dsprites17}. This dataset consists of 2D shapes generated from 6 ground truth independent latent factors.}
\label{fig:dspritesexamples}
\vspace{-5mm}
\end{wrapfigure}

\section{Applications of counterfactual invariance}\label{sec:counterfactual_invariance}
To clarify the benefits and challenges of counterfactual invariance, we now provide brief explanations of potential application areas, which we will subsequently study empirically with real-world data.

\paragraph{Image classification.} Counterfactual invariance serves as a strong notion of robustness in high-dimensional settings, such as image classification: ``Would the truck have been classified correctly had it been winter in this exact situation instead of summer?'' For concrete demonstration, we use the dSprites dataset \citep{dsprites17}, which consists of relatively simple, yet high-dimensional, square black and white images of different shapes (squares, ellipses, etc.), sizes, and orientations (rotation) in different xy-positions (see Figure \ref{fig:dspritesexamples}).
Since this dataset is fully synthetic and labelled (we know all factors for each image), we consider a causal model as depicted in Figure \ref{fig:graphs}(c).
The full structural equations are provided in Appendix \ref{app:sec:imageresults}.
The main goal in this application is to evaluate our method for high-dimensional data, including mixed categorical (e.g., shapes) and continuous (e.g., position) variables.
Therefore, we set the outcome $\textbf{Y}$ to be a function of all observed and unobserved variables, which can be inferred from the high-dimensional images.

\paragraph{Counterfactual fairness.}
We first show that the popular definition of counterfactual fairness \citep{DBLP:conf/nips/KusnerLRS17} is an incarnation of counterfactual invariance (see Definition \ref{def:counterfactual_invariance}).
Counterfactual fairness corresponds to the commonly posed question, after receiving a consequential decision from an algorithm: ``Would I have gotten the same outcome had I been a different gender, race, or age with all else being equal?''.
Again, we denote the outcome by $\mathbf{Y} \subset \mathbf{V}$ and so-called \emph{protected attributes} by $\mathbf{A} \subseteq \mathbf{V} \setminus \mathbf{Y}$. These attributes typically constitute demographic features such as gender, race, or age protected under anti-discrimination laws (see, e.g., \citet{us_law}).
Collecting all remaining observed covariates into $\textbf{W} \coloneqq \mathbf{V} \setminus \mathbf{Y}$, we define counterfactual fairness as follows, showing that it is a special case of counterfactual invariance.
\begin{definition}[Counterfactual Fairness, Definition 5 by \citet{DBLP:conf/nips/KusnerLRS17}]
\label{def:counterfactual_fairness}
A predictor $\hat{\mathbf{Y}}$ is counterfactually fair if under any context $\mathbf{W} = \mathbf{w}$ and $\mathbf{A} = \mathbf{a}$, it holds $\mathbb{P}_{\hat{\mathbf{Y}}^*_{\mathbf{a}} \mid \mathbf{W}, \mathbf{A}}(\mathbf{y}\mid \mathbf{w}, \mathbf{a}) = \mathbb{P}_{\hat{\mathbf{Y}}^*_{\mathbf{a}'} \mid \mathbf{W}, \mathbf{A}}(\mathbf{y}\mid \mathbf{w}, \mathbf{a})$, for all $\mathbf{y}$ and for any value $\mathbf{a}'$ attainable by $\mathbf{A}$.
\end{definition}
For a concrete example, we consider the commonly used ``Adult'' dataset from the UCI repository \citep{uci_dataset}.
Our goal is to predict whether an individual's income is above a certain threshold $\mathbf{Y} \in \{0,1\}$ based on a collection of (demographic) information including protected attributes such as gender and age.
From a fairness perspective, such a predictor could be used in loan-granting decisions, where one must not discriminate against certain demographic groups. We follow \citet{nabi2018fair,DBLP:conf/aaai/Chiappa19}, where a subset of variables are selected from the dataset and a causal structure is assumed as in Figure \ref{fig:graphs}(b).
We choose gender (considered binary in this dataset) and age (considered continuous) as the protected attributes $\mathbf{A}$. We denote the marital status, level of education, occupation, working hours per week, and work class jointly by $\mathbf{X}$ and combine the remaining observed attributes in $\mathbf{C}$.
With this notation, our goal is to learn a predictor $\hat{\mathbf{Y}}$ that is counterfactually invariant (or fair) in $\mathbf{A}$ with respect to $\mathbf{W} = \mathbf{C} \cup \mathbf{X}$.
We remark that achieving fairness for continuous or even mixed categorical and continuous protected attributes is an ongoing line of research (even for non-causal fairness notions) \citep{Mary19:Fairness,chiappa2021fairness}, but directly supported by \textsc{HSCIC}{}.

\section{Experiments}\label{sec:new_experiments}

We begin our empirical assessment of \textsc{HSCIC}{} on synthetic datasets, where ground truth is known and systematically study how performance is influenced by the dimensionality of the variables. We then move to the dSprites and Adult datasets described in \ref{sec:counterfactual_invariance}.

\begin{wrapfigure}{r}{0.4\textwidth}
\vspace{-5mm}
\centering
   \begin{tikzpicture}[node distance=7mm and 10mm, thick, main/.style = {draw, circle, minimum size=0.8cm}, >={triangle 45}]
    \node[main] (1) {$\scriptstyle \mathbf{A}$};
    \node[main] (2) [right =of 1] {$\scriptstyle \mathbf{X}$};
    \node[main] (3) [right =of 2] {$\scriptstyle \mathbf{Y}$};
    \node[main] (4) [above =of 2] {$\scriptstyle \mathbf{Z}$};
    \draw[->] (1) -- (2);
    \draw[->] (2) -- (3);
    \draw[->] (4) -- (1);
    \draw[->] (4) -- (2);
    \draw[->] (4) -- (3);
    \draw[->] (1) to [out=320, in=220, looseness=1] (3);
  \end{tikzpicture}
  \caption{\textbf{Assumed causal structure for the synthetic experiments} (see Section \ref{sec:synthetic_experiments}). The precise corresponding generative random variables are described in Appendix \ref{app:synthetic}.}\label{fig:graph_experiment_1}
\vspace{-3mm}
\end{wrapfigure}

\subsection{Synthetic experiments}\label{sec:synthetic_experiments}
We simulate different datasets following to the causal graphical structure in Figure \ref{fig:graph_experiment_1}. The datasets are composed of four sets of observed continuous variables (i) the prediction target $\mathbf{Y}$, (ii) the variable(s) we want to be counterfactually invariant in $\mathbf{A}$, (iii) observed covariates that mediate effects from $\mathbf{A}$ on $\mathbf{Y}$, and (iv) observed confounding variables $\mathbf{Z}$.
The goal is to learn a predictor $\hat{\mathbf{Y}}$ that is counterfactually invariant in $\mathbf{A}$ with respect to $\mathbf{W} \coloneqq \mathbf{A} \cup \mathbf{X} \cup \mathbf{Z}$. We consider various artificially generated datasets for this example, which mainly differ in the dimension of the observed variables and their correlations. A detailed description of each dataset is deferred to Appendix \ref{app:synthetic}.

For all synthetic experiments, we train the model using fully connected networks (MLPs). We use the \textbf{MSE} loss $\mathcal{L}_{\textsc{mse}}( \hat{\mathbf{Y}})$ as the predictive loss $\mathcal{L}$ in eq.~\eqref{eq:total_loss} for continuous outcomes $\mathbf{Y}$. We generate $4$k samples from the observational distribution in each setting and use a 80:20 train/test split.
All metrics reported are on the test set.
We perform hyper-parameter tuning for MLP hyperparameters based on a random strategy (see Appendix \ref{app:synthetic} for details).
The $\hscic{\hat{\mathbf{Y}},\mathbf{A}}{\mathbf{Z}}$ term is computed as in eq.~\eqref{eq:estimate_hscic} using a Gaussian kernel with amplitude $1.0$ and length scale $0.1$.
The regularization parameter $\lambda$ for the ridge regression coefficients in eq.~\eqref{eq:cond_measure_estimand} is set to $\lambda = 0.01$.
To estimate \textbf{VCF} we set $d=1000$ and $k=500$.

\paragraph{Model performance with the use of the HSCIC.} We first perform a set of experiments to study the effect of the HSCIC, and to highlight the trade-off between accuracy and counterfactual invariance. For this set of experiments, we generate a dataset as described in Appendix \ref{appendix:synthetic_experiment}. Figure \ref{fig:simulated_exp} (left) shows the \textbf{VCF} for different regularization parameters $\gamma$, as a function of the values of $\mathbf{Z}$. We observe that increasing values of $\gamma$ lead to a consistent decrease of \textbf{VCF}. Figure \ref{fig:simulated_exp} (center) shows the values attained by the \textsc{HSCIC}{} and \textbf{MSE} for increasing parameter $\gamma$, demonstrating the expected trade-off in raw predictive performance and enforcing counterfactual invariance.
Finally, Figure \ref{fig:simulated_exp} (right) highlights the usefulness of \textsc{HSCIC}{} as a measure of counterfactual invariance, being in strong agreement with \textbf{VCF}.

\begin{figure}
  \centering
    \includegraphics[width=0.295\linewidth]{scatterplot_VCF.pdf}
    \includegraphics[width=0.34\linewidth]{plot_mse_hscic.pdf}
    \includegraphics[width=0.34\linewidth]{plot_vcf_hscic.pdf}
  \caption{%
  \textbf{(Left) variance of the counterfactual distributions} for $d=100$ datapoints with lines representing $3^{\mathrm{rd}}$-order polynomial regression and shaded areas being 95\% confidence intervals.
  \textbf{(Center) trade-off between accuracy and counterfactual invariance.} We observe that the \textsc{HSCIC}{} decreases, as the \textbf{MSE} increases. Vertical bars denote standard errors over $15$ different random seeds. \textbf{(Right) Correspondence between the \textsc{HSCIC}{} and the \textbf{VCF}}, for increasing regularizer $\gamma$. Again, vertical bars denote standard errors over $15$ different random seeds.}
  \label{fig:simulated_exp}
\end{figure}

\paragraph{Comparison with baselines.}
We perform a second set of experiments to compare our method against baselines. For this set of experiments we consider two settings, which we refer to as Scenario 1 and Scenario 2. These two settings differ in how the conditioning set $\textbf{Z}$ affects the outcome $\textbf{Y}$, with
Scenario 1 exhibiting higher correlation of $\textbf{Z}$ with both the mediator $\textbf{X}$ and the outcome $\textbf{Y}$. We refer the reader to Appendix \ref{appendix:baselines} for a precise description of Scenario 1 and Scenario 2. We remark that counterfactually invariant training has not received much attention yet. This limits the choice of baselines for experimental comparison. The main work on this topic to date does not provide a suitable benchmark for experimental comparison \citep{DBLP:conf/nips/VeitchDYE21}. In fact, the method by \citet{DBLP:conf/nips/VeitchDYE21} applies to a weaker notion of counterfactual invariance, and it can be used only on specific graphical structures which differ from our set-up. However, since counterfactual fairness is a special case of counterfactual invariance, we can compare to two methods proposed by \citet{DBLP:conf/nips/KusnerLRS17} (in applicable settings). We compare to the \textit{Level 1} (only use non-descendants of $\mathbf{A}$ as inputs to $\hat{\mathbf{Y}}$) and the \textit{Level 3} (assume an additive noise model and in addition to non-descendants, only use the residuals of descendants of $\mathbf{A}$ after regression on $\mathbf{A}$ as inputs to $\hat{\mathbf{Y}}$) approaches of \citet{DBLP:conf/nips/KusnerLRS17}.
We refer to these two baselines as CF1 and CF2 respectively. Note that CF1 satisfies counterfactual invariance perfectly by construction. We compare our method against CF1 and CF2 on both Scenario 1 and Scenario 2. We summarize the results in Table \ref{table:baseline_comparison}.
For a suitable choice of $\gamma$, our method outperforms the baseline CF2 in both \textbf{MSE} and \textbf{VCF} simultaneously.

\paragraph{Multi-dimensional variables.}
We perform a third set of experiments to assess \textsc{HSCIC}{}'s performance in higher dimensions. We consider simulated datasets as described in Appendix \ref{multid_setting}.
In this set of experiments, we independently increase the dimension of $\textbf{Z}$ and $\textbf{A}$ in two different simulations, leaving the rest of the variables unchanged. \Cref{table:multidimensional_results_A,table:multidimensional_results_Z} present the results for different regularization coefficients $\gamma$ and different dimensions of $\textbf{A}$ and $\textbf{Z}$. These results demonstrate that \textsc{HSCIC}{} can handle multi-dimensional variables while maintaining performance, as counterfactual invariance is approached when $\gamma$ increases.

\begin{table}
  \centering
  \caption{\textbf{Performance of the \textsc{HSCIC}{} against baselines} CF1 and CF2 on two synthetic datasets (see Appendix \ref{appendix:baselines}). Notably, for $\gamma=5$ in Scenario 1 and $\gamma=13$ in Scenario 2 we outperform CF2 in \textbf{MSE} and \textbf{VCF} simultaneously.}
  \label{table:baseline_comparison}
  \begin{tabularx}{\columnwidth}{l*6{Z}}
     \toprule
     & \multicolumn{2}{c}{Scenario 1}
     & \multicolumn{2}{c}{Scenario 2} \\
     \cmidrule(lr){2-3} \cmidrule(l){4-5}
     & \textbf{MSE} {\color{gray}$\times 10^3$} &$\textbf{VCF}${\color{gray}$\times 10^3$} & \textbf{MSE}{\color{gray}$\times 10^3$} & $\textbf{VCF}${\color{gray}$\times 10^3$} \\
   \midrule
   $\gamma=0$  & $0.36 \pm 0.5 $ &  $20.33 \pm 1$ & $2 \pm 0.1$ &  $239 \pm 0.4$ \\
   $\gamma=5$   & $16.97 \pm 2$ & $6.84 \pm 1$ & $43 \pm 8$ &  $230 \pm 6$ \\
   $\gamma=10$ &$19.80\pm 2$ & $5.89 \pm 0.6$ & $50 \pm 8$ & $190 \pm 3$ \\
   $\gamma=13$   &$22 \pm 1$ &$5.78 \pm 0.5$ &  $133 \pm 6$ & $157 \pm 7$ \\
   \midrule
   \textrm{CF1}  & $24.44 \pm 3$  & $0$ & $218 \pm 10$  & $0$ \\
   \textrm{CF2}   & $19.50 \pm 2$ & $7.50 \pm 1$ & $137 \pm  4$ & $167 \pm 2$ \\
   \bottomrule

   \end{tabularx}
   \vspace{10pt}
\end{table}

\begin{table}
  \centering
  \caption{\textbf{Results of the \textbf{MSE} ($\times10^5$ for readability), \textsc{HSCIC}{}, \textbf{VCF} for increasing dimension of $\textbf{A}$}, on a synthetic dataset as in Appendix \ref{multid_setting}. In this case, the dimension of the other variables are $\text{dimZ}=1, \text{dimX}=1,\text{dimY}=1$.}
  \label{table:multidimensional_results_A}
  \scriptsize
  \begin{tabularx}{\columnwidth}{@{}Yrrrrrrrrrrrr@{}}
   \toprule
   & \multicolumn{3}{c}{dimA=2}
   & \multicolumn{3}{c}{dimA=5}
   & \multicolumn{3}{c}{dimA=10}
   & \multicolumn{3}{c}{dimA=20} \\
   \cmidrule(lr){2-4} \cmidrule(lr){5-7} \cmidrule(lr){8-10} \cmidrule(l){11-13}
   $\gamma$ & \textbf{MSE} & \textbf{HSCIC} & \textbf{VCF} & \textbf{MSE} & \textbf{HSCIC} & \textbf{VCF} & \textbf{MSE} & \textbf{HSCIC} & \textbf{VCF} & \textbf{MSE} & \textbf{HSCIC} & \textbf{VCF} \\
   \midrule
   $0$   & $35$ & $0.0230$ & $0.025$ & $1.53$ & $0.0177$ & $0.013$ & $10$  & $0.0150$ &  $0.004$& $1.45$& $0.0126$ & $0.0010$ \\
   $\tfrac{1}{2}$  & $86$ & $0.0214$ & $0.022$  & $1580$ & $0.0157$ & $0.010$ & $40$ & $0.0135$& $0.002$ &$33$ & $0.0115$ & $0.0004$ \\
   $1$   & $210$ & $0.0200$ &  $0.220$ & $1450$ & $0.0144$  & $0.008$  & $100$ & $0.0127$ & $0.001$ & $60$& $0.0111$ & $0.0003$ \\
   \bottomrule
   \end{tabularx}%
   \vspace{10pt}
\end{table}

\begin{table}[t]
  \centering
  \caption{\textbf{Results of \textbf{MSE} ($\times 10^5$ for readability), \textsc{HSCIC}{}, \textbf{VCF} for increasing dimension of $\textbf{Z}$}, on a synthetic dataset as in Appendix \ref{multid_setting}. The dimension of the other variables are
  $\text{dimA}=1, \text{dimX}=1,\text{dimY}=1$.}
  \label{table:multidimensional_results_Z}
  \begin{tabularx}{\columnwidth}{@{}Yrrrrrrrrr@{}}
   \toprule
   & \multicolumn{3}{c}{dimZ=2}
   & \multicolumn{3}{c}{dimZ=3}
   & \multicolumn{3}{c}{dimZ=5} \\
   \cmidrule(lr){2-4} \cmidrule(lr){5-7} \cmidrule(l){8-10}
   $\gamma$ & \textbf{MSE} & \textbf{HSCIC} & \textbf{VCF} & \textbf{MSE} & \textbf{HSCIC} & \textbf{VCF} & \textbf{MSE} & \textbf{HSCIC} & \textbf{VCF} \\
   \midrule
   $0$   & $20$ & $0.002024$ & $6.15 \times 10^{-5}$ & $11$ & $0.03407$ & $0.023$ & $53.9$ & $0.033180$ & $0.0062$ \\
   $\tfrac{1}{2}$    & $20$ & $0.002023$ & $1.29 \times 10^{-5}$ & $49$ &$0.03406$ & $0.016$ & $21.9$ & $0.033178$ & $0.0061$ \\
   $1$    &  $20$ & $0.002019$ & $4.40 \times 10^{-6}$ & $57$ & $0.03403$& $0.009$ &$18$ & $0.033176$ & $0.0039$\\
   \bottomrule
   \end{tabularx}%
   \vspace{10pt}
\end{table}

\subsection{High-dimensional image experiments}\label{sec:imagedata}

\begin{figure}
  \centering
  \includegraphics[width=0.32\linewidth]{plot_mse_imagedataset.pdf}
  \includegraphics[width=0.32\linewidth]{plot_hscic_imagedataset.pdf}
  \includegraphics[width=0.32\linewidth]{plot_vcf_imagedataset.pdf}
  \caption{\textbf{Results of \textbf{MSE}, \textsc{HSCIC}{} operator and \textbf{VCF} for the dSprites image dataset experiment.} The \textsc{HSCIC}{} operator decreases steadily with higher values of $\gamma$. Similarly, a necessary increase of \textbf{MSE} can be observed. For both $\gamma=0.5$ and $\gamma=1$ an overall decrease of \textbf{VCF} is observed compared to the not-regularized setting. However, probably due to the high dimensional setting, the analysed metrics do not relevantly vary for the two regularized settings.
  }
  \label{fig:image_exp}
\end{figure}
We perform a set of experiments on a high-dimensional real-world dataset. We consider the image classification task as described in \ref{sec:counterfactual_invariance}. In this set of experiments, we use the dSprites dataset with a causal graph as in Figure \ref{fig:graphs}(c).
The full structural equations are provided in Appendix \ref{app:sec:imageresults}, where we assume a causal graph over the determining factors of the image, and essentially look up the corresponding image in the simulated dataset.
This experiment is particularly challenging due to the mixed categorical and continuous variables in $\mathbf{C}$ (\texttt{shape}, \texttt{y-pos}) and $\mathbf{X}$ (\texttt{color}, \texttt{orientation}), continuous $\mathbf{A}$ (\texttt{x-pos}).
All variables except for scale are assumed to be observed, and all variables jointly with the actual image determine the outcome $\textbf{Y}$.
Our goal is to learn a predictor $\hat{\mathbf{Y}}$ that is counterfactually invariant in the x-position with respect to all other observed variables.
In the chosen causal structure, $\{\texttt{shape}, \texttt{y-pos}\}\in \mathbf{C}$ block all non-directed paths from $\mathbf{A} \cup \mathbf{C} \cup \mathbf{U}$ to $\mathbf{Y}$.
Hence, we seek to achieve $\hat{\mathbf{Y}} \perp \!\!\! \perp \texttt{x-pos} \mid \{\texttt{shape}, \texttt{y-pos}\}$ via the \textsc{HSCIC}{} operator.
To accommodate the mixed input types, $\hat{\mathbf{Y}}$ puts an MLP on top of features extracted from the images via convolutional layers concatenated with features extracted from the remaining inputs via an MLP.

Figure \ref{fig:image_exp} demonstrates that \textsc{HSCIC}{} achieves improved \textbf{VCF} as $\gamma$ increases up to a certain point while affecting \textbf{MSE}, a necessary trade-off.
Adjusting our method to achieve higher stability is an interesting direction for future work.

\subsection{Fairness with continuous protected attributes}

In this experiment we apply our method to the popular UCI Adult dataset \citep{uci_dataset}, described in Section \ref{sec:counterfactual_invariance}. Following previous related work \citep{chiappa2021fairness}, we assume the causal struture depicted in the Appendix \ref{app:sec:adult}. Furthermore, we assume additional latent variables as in Figure \ref{fig:causal_graph_uci}.
Our goal is to learn a classifier for whether

\begin{figure}
\begin{minipage}[t]{.28\textwidth}
\vspace{0pt}
  \includegraphics[width=\linewidth]{kdeplot_VCF_uci.png}
\vspace{0pt}
\end{minipage}
\hspace{-7mm}
\hfill
\begin{minipage}{.4\textwidth}
\hspace{-3.9cm}
  \vspace{-2cm}
  \scriptsize
  \begin{tabular}{@{} l *{3}{c} @{}}
   \toprule
   & \textbf{MSE} & \textbf{HSCIC}& \textbf{VCF}\\
   \midrule
   $\gamma=0$   & $(81.725 \pm  0.0481) \%$ & $0.01067 \pm 3.4 \times 10^{-5}$ & $0.06685 \pm 0.006$\\
   $\gamma=1$   & $(81.547 \pm 0.0745) \%$ & $0.00739 \pm 6.8 \times 10^{-5}$&$0.05156 \pm 0.003$\\
   \bottomrule
  \end{tabular}
\vspace{-1cm}
\end{minipage}
\vspace{-2mm}
  \caption{\textbf{(Left) Probability density mass of the \textbf{VCF} distribution}, for different choices of $\gamma$. We observer less variance and more mass near zero for regularization parameter  $\gamma=1$. \textbf{(Right) Results of accuracy, \textsc{HSCIC}{} and \textbf{VCF}}, showing a considerable improvement of \textsc{HSCIC}{} and \textbf{VCF} for $\gamma=1$.}
  \label{fig:fairness_2}
  \vspace{-3mm}
\end{figure}

 the income of an individual is above a certain threshold, which is counterfactually invariant in both age and gender with respect to a collection of additional demographic attributes.
Here, we use an MLP with binary cross-entropy loss for $\hat{\mathbf{Y}}$.
Since this experiment is based on real data, the true counterfactual distribution cannot be known.
For an evaluation of our method, we thus again follow \citet{chiappa2021fairness} and estimate a possible true underlying structure by inferring the posterior distribution over the unobserved variables using variational autoencoders \citep{KingmaW13}.
Even though this can only provide approximate variances of counterfactual distributions, Figure \ref{fig:fairness_2} (left) shows that \textsc{HSCIC}{} still achieves more counterfactually fair outcome distributions (more mass near zero) than an unconstrained classifier ($\gamma = 0$).
Figure \ref{fig:fairness_2} (right) highlights once more that the \textsc{HSCIC}{} operator reliably decreases also in classification setting, again trading off accuracy.

\section{Discussion and future work}\label{sec:discussion}

We study the problem of learning predictors $\hat{\mathbf{Y}}$ that are ``counterfactually invariant'' in changes of certain covariates.
We put forth a formal definition of counterfactual invariance and describe how it generalizes existing notions.
Next, we provided a novel graphical criterion to characterize counterfactual invariance and reduce it to a conditional independence statement. Our method does not require identifiability of the counterfactual distribution or exclude the possibility of unobserved confounders.
Finally, we propose an efficiently estimable, model-agnostic regularization term to enforce this conditional independence based on kernel mean embeddings of conditional measures, which works for mixed continuous/categorical, multi-dimensional variables. We demonstrate the efficacy of our method in regression and classification tasks involving controlled synthetic scenarios, high-dimensional images, and in a fairness application, where it outperforms existing baselines.

The main limitation of our work, as for all studies in this domain, is the assumption that the causal graph is known. Another limitation is that our methodology is applicable when a graphical criterion is satisfied, i.e., the required variables are observed. From an ethics perspective, the increased robustness of counterfactually invariant (or societal benefits of counterfactually fair) predictors are certainly desirable.
However, this presupposes that all the often untestable assumptions are valid. Overall, causal methodology should not be applied lightly, especially in high-stakes and consequential decisions. A critical analysis of the broader context or systemic factors may hold more promise for societal benefits, than a well-crafted algorithmic predictor.

In light of the limitations mentioned above, an interesting and important direction for future work is to assess the sensitivity of our method to misspecifications of the causal graph or insufficient knowledge of the required blocking set. Finally, we believe our graphical criterion and KMEs-based regularization can also be useful for causal representation learning, where one aims to isolate causally relevant, autonomous factors underlying the data generating process of high-dimensional data.

\subsubsection*{Acknowledgement}
This work was supported by the German Federal Ministry of Education and Research (BMBF): T\"ubingen AI Center, FKZ: 01IS18039B.

\bibliographystyle{apalike}

redpajama_set_name: RedPajamaArXiv