# Verified Training for Counterfactual Explanation Robustness under Data Shift

Anna P. Meyer , Yuhao Zhang∗, Aws Albarghouthi & Loris D’Antoni
Department of Computer Sciences
University of Wisconsin - Madison
Madison, WI, USA
{annameyer, yuhaoz, aws<EMAIL_ADDRESS>
Equal contribution

###### Abstract

Counterfactual explanations (CEs) enhance the interpretability of machine
learning models by describing what changes to an input are necessary to change
its prediction to a desired class. These explanations are commonly used to
guide users’ actions, e.g., by describing how a user whose loan application
was denied can be approved for a loan in the future. Existing approaches
generate CEs by focusing on a single, fixed model, and do not provide any
formal guarantees on the CEs’ future validity. When models are updated
periodically to account for data shift, if the generated CEs are not _robust_
to the shifts, users’ actions may no longer have the desired impacts on their
predictions. This paper introduces VeriTraCER, an approach that _jointly_
trains a classifier and an explainer to explicitly consider the robustness of
the generated CEs to small model shifts. VeriTraCER optimizes over a carefully
designed loss function that ensures the _verifiable_ robustness of CEs to
local model updates, thus providing deterministic guarantees to CE validity.
Our empirical evaluation demonstrates that VeriTraCER generates CEs that (1)
are verifiably robust to small model updates and (2) display competitive
robustness to state-of-the-art approaches in handling empirical model updates
including random initialization, leave-one-out, and distribution shifts.

##  1 Introduction

Machine learning models are increasingly used to support decision-making in
sectors such as banking, education, social services, and criminal justice. Due
to the high stakes of these decision-making settings, and the fact that model
internals are often both proprietary and too complex for humans to understand
directly, laws such as the GDPR (European Commission, 2018) and ECOA (Consumer
Financial Protection Bureau, 2018) require that explanations be offered to
users who are subject to these models’ predictions. Explanations commonly take
the form of counterfactual explanations (CEs), which serve as guidelines for
how an input can change in order to receive a different decision in the
future. For instance, an individual who is denied a loan may receive a CE that
says their application would have been accepted had their salary been $$5000$
higher. If the applicant wishes to obtain the loan, they can work to increase
their salary and then reapply.

However, machine learning models must be periodically updated to account for
new data and avoid declining performance due to distribution shift. Even if
the CE is valid (i.e., produces the desired prediction) at the time it is
generated, there is no guarantee that it will remain valid after one or more
routine model updates. So, the individual who successfully raises their salary
by the requested $$5000$ may still be rejected if they reapply for a loan 6
months later, because the model internals shifted.

To preserve the validity of CEs across model shifts, we want to generate CEs
that are _robust_ to small model changes. Existing work aims to generate
robust CEs by increasing the distance from the original input to its CE
(Hamman et al., 2023), finding CEs in areas of low Lipschitz constants (Black
et al., 2021), or using a minmax objective in the CE generation process
(Upadhyay et al., 2021). While all of these methods yield improved robustness
over standard training, they fail to provide _formal guarantees_ on the CEs’
robustness, they require solving _expensive_ optimization problems to generate
each CE, and they generate CEs with respect to a _fixed_ model, even if
solving the same optimization problem for a different model – such as one that
may be adopted in the future – would yield a better and _still valid_ CE.

Our approach, VeriTraCER, fundamentally reframes the problem of generating
robust CEs in two ways: first, we consider the _multiplicity set_ of similar
models, noting that this set likely contains models that may be adopted in the
future when accounting for slight data shifts. Second, we adapt existing work,
CounterNet, which considers the model training and CE generation processes as
a single pipeline (Guo et al., 2023). Specifically, VeriTraCER uses the
existing conception of training a model that jointly performs classification
and CE generation, but promotes CE robustness by optimizing CE generation over
a multiplicity set of classifiers, rather than a single classifier. A key part
of how VeriTraCER obtains robustness is to use verified training (Gowal et
al., 2018; Zhang et al., 2020) to _deterministically certify_ when the
predictions and the CEs our model generates are robust to small model changes,
as indicated by an $l_{p}$-bound on the classifier’s parameters, subject to
restrictions that hold the original prediction constant. We develop a new
variation of verified training, Simul-CROWN, that allows us to obtain tighter
bounds than existing approaches. We show that CEs produced by VeriTraCER are
not only certifiably robust to small model shifts, but are also robust to
other empirical forms of model updates—e.g., training with different random
seeds or different data subsets. An added benefit is that CE generation using
VeriTraCER is very fast (equivalent to inference speed), with only a modest
increase in training time.

In summary, we make the following key contributions: (1) We propose robust CE
generation over a multiplicity set of models, (2) we develop a loss function
to jointly train an accurate model and a robust and valid CE generator, (3) we
design a new verified training algorithm, Simul-CROWN, to soundly
overapproximate our loss function during training, and (4) we show that our
technique, Verified Training for CE Robustness (VeriTraCER), achieves high CE
robustness, both to $l_{p}$-bounded model shifts and real-world model updates.

##  2 Related Work

##### Explainable AI

Explanations for model predictions can come from insights about the model
(linear regression coefficients, decision tree rules, or feature activations),
or more commonly for black-box models, from a post-hoc technique. Post-hoc
explanation techniques are typically either feature-based (e.g., Lundberg &
Lee (2017); Ribeiro et al. (2016); Simonyan et al. (2013); Smilkov et al.
(2017)) and aim to describe which input features are relevant, or
counterfactual-based (e.g., Karimi et al. (2020); Looveren & Klaise (2019);
Poyiadzi et al. (2020); Ustun et al. (2018); Wachter et al. (2018)) and aim to
describe how the original instance needs to change to get a different
prediction. In this work, we focus on CEs due to the fact that they provide
direct guidance to users. CEs are typically found post-hoc, based on an
optimization problem over a fixed model. An exception is the CounterNet
training procedure (Guo et al., 2023), which jointly trains a model and CE
generator.

##### Robust explanations

“Explanation robustness” refers to multiple phenomena in the literature; the
focus can be robustness with respect to the input (i.e., whether perturbing
the input yields a similar explanation) or with respect to model changes
(i.e., whether changing the model yields a similar prediction and explanation
for a fixed input). In this work, we focus on the latter definition.

Various works have explored robust CE generation for specific model classes,
such as for tree-based ensembles (Dutta et al., 2022; Forel et al., 2022) or
(locally) linear models (Bui et al., 2022; Nguyen et al., 2023). Other work on
more complex model classes has shown that it is possible to adversarially
change a model to keep identical predictions, but drastically change the
associated feature-based explanation (Anders et al., 2020; Heo et al., 2019;
Slack et al., 2020). Similarly, a given CE can be invalidated by another
equivalently-accurate model in a neural network setting (Hamman et al., 2023).
But it is still possible to improve (average) explanation robustness to model
shift, such as by increasing the local smoothness of the model (Black et al.,
2021; Dombrowski et al., 2022; Meyer et al., 2023; Srinivas et al., 2022). For
CEs, increasing the distance to the CE is also commonly thought to increase
robustness to model shift (Jiang et al., 2023; Pawelczyk et al., 2020).
However, this heuristic does not always work for deep models (Black et al.,
2021); instead, the agreement of points in an epsilon-ball around the CE is
important (Hamman et al., 2023). However, all of these papers consider CE
generation with respect to a fixed model – by contrast, our process jointly
trains a model and CE generator in order to obtain robustness. Also, most
existing methods only evaluate CE robustness based on empirical model shift.
An exception is Jiang et al. (2023), who consider CE robustness to a set of
models (via $l_{p}$-norm-bounded changes to the parameters of a fixed model).
We consider the same $l_{p}$-norm bounded setting, but also evaluate our
approach on empirical model shifts.

##  3 Problem Definition

Let $f$ be a neural net classifier that learns its parameters $\theta_{f}$
with a loss function $\mathcal{L}$ and a training set
$D=\\{(\mathbf{x}_{i},y_{i})\\}_{i=1}^{n}$. We will assume that $f$ is a
binary classifier, i.e., assigns $\mathbf{x}$ a label $y\in\\{0,1\\}$. A
distribution shift may occur and yield a new training set
$D^{\prime}=\\{\mathbf{x}_{i},y_{i})\\}_{i=1}^{k}$. We _finetune_ $f$ (i.e.,
train for an additional small number of epochs) or retrain from scratch in the
presence of $\ell_{2}$ regularization to yield a shifted model
$f_{\mathrm{m}}$. Given $f$ and $f_{\mathrm{m}}$ with parameters $\theta_{f}$
and $\theta_{f_{\mathrm{m}}}$, respectively, we measure their distance
$\|\theta_{f}-\theta_{f_{\mathrm{m}}}\|_{p}$ as the largest $\ell_{p}$
parameter distance among all layers (both weight matrices and bias vectors) in
$f$ and $f_{\mathrm{m}}$. We make the assumption that the finetuning (or
retraining with regularization) process will yield a model $f_{\mathrm{m}}$
such that $\|\theta_{f}-\theta_{f_{\mathrm{m}}}\|_{p}\leq\delta$ for some
small $\delta$. For simplicity, we use a single $\delta$ for the $l_{p}$ bound
across all layers. However, our approach easily extends to using different
$\delta$ values for different layers.

##### Counterfactual explanations

A counterfactual explanation (CE) generator $g$ returns a CE
$\mathbf{x}^{\prime}$ with respect to a function $f$ and an original input
$\mathbf{x}$, i.e., $\mathbf{x}^{\prime}=g(f,\mathbf{x})$. We say that
$\mathbf{x}^{\prime}$ is _valid_ if it satisfies $f(\mathbf{x}^{\prime})\neq
f(\mathbf{x})$. The distance between $\mathbf{x}$ and $\mathbf{x}^{\prime}$ is
called _proximity_ and should be minimized subject to the validity constraint.
Additional constraints (e.g., not changing immutable features like race or
gender) can be placed on the CE generation process, e.g., by using a custom
distance metric that heavily penalizes changing those features. Our approach
is compatible with that type of modification, but we assume an $l_{1}$-norm
distance metric for simplicity.

##### Multiplicity set of $f$

Even though $f$ is the result of an optimization process, it likely is not the
only model that performs well on the training data due to model multiplicity
(D’Amour et al., 2022; Marx et al., 2020). Furthermore, the finetuning or
retraining process will yield a distinct model $f_{\mathrm{m}}$ that is close
to $f$. We adopt a common assumption from the literature (Jiang et al., 2023;
Upadhyay et al., 2021) that the parameter distance between $f$ and
$f_{\mathrm{m}}$ will be bounded by $l_{p}$ norms. We define the _multiplicity
set_ with respect to a model $f$ and an input $\mathbf{x}$ as follows.

######  Definition 3.1 (Multiplicity set).

Given a model $f$ with parameters $\theta_{f}$, an input $\mathbf{x}$, an
$l_{p}$ norm, and a bound $\delta$, we define the $\delta$-robust multiplicity
set as $\mathcal{M}_{f,\mathbf{x}}=\\{f_{\mathrm{m}}\mid
f(\mathbf{x})=f_{\mathrm{m}}(\mathbf{x})\wedge\|\theta_{f}-\theta_{f_{\mathrm{m}}}\|_{p}\leq\delta\\}$.

Intuitively, the multiplicity set contains models that have comparable
performance to $f$ and have a similar set of model weights. Crucially, it also
is intended to contain models $f_{m}$ that may be the result of model updates
due to distribution shift. Note that our definition limits us to models with
the same prediction on $\mathbf{x}$. As we are primarily interested in CEs, it
no longer makes sense to consider adopting the CE if the prediction changes.

##### CE robustness

_Data shift necessitates that models change over time_ , thus creating the
risk that CEs generated by a model will be invalidated in the future. We
define _CE robustness_ as follows.

######  Definition 3.2 ($\mathcal{M}_{f,\mathbf{x}}$-robustness of a CE).

Given a binary classifier $f$ and an input $\mathbf{x}$, we say that a CE
$\mathbf{x}^{\prime}$ is _robust_ if it is a valid CE for all models in the
multiplicity set, i.e., if $\forall
f_{\mathrm{m}}\in\mathcal{M}_{f,\mathbf{x}}$,
$f_{\mathrm{m}}(\mathbf{x}^{\prime})=1-f(\mathbf{x})$.

##### Goal

Our goal is to devise a training algorithm that yields a model $f$ and a CE
generator $g$ such that the CEs generated by $g$ will be
$\mathcal{M}_{f,\mathbf{x}}$-robust (i.e., robust to small changes in the
model $f$). In other words, we want to maximize the number of robust CEs that
we generate on some dataset (e.g., on the training dataset during training, or
on a test or validation dataset post-training). Note that our goal differs
from all existing works, which propose algorithms to generates robust CEs
according to a fixed model $f_{\mathrm{fix}}$. The following example
illustrates how focusing on the model selection process—rather than just on
the CE generation technique—can yield better and more robust CEs.

(a) $f_{0}=x_{1}-x_{2}-2$. The optimal robust CE of $\mathbf{x}=(4,1)$ is
$(0,0)$.

(b) $f_{1}=x_{1}-x_{2}-1.5$. There exists no robust CE of $\mathbf{x}=(4,1)$.

(c) $f_{2}=x_{1}-x_{2}-2.5$. The opt. robust CE of $\mathbf{x}=(4,1)$ is
$(0,0.5)$.

Figure 1: Plots of three linear models and their multiplicity sets. The black
line shows the original linear model $f_{i}$. The green and red regions
contain all samples that are robust under the multiplicity set
$\mathcal{M}_{f_{i},\mathbf{x}}$, receiving predictions $+$ and $-$,
respectively. The yellow region corresponds to all samples that are not robust
under the multiplicity set $\mathcal{M}_{f_{i},\mathbf{x}}$.

######  Example 3.1.

Suppose we have an input $\mathbf{x}=(4,1)$ and the three equally-accurate
linear models $f_{0},f_{1},f_{2}$ shown in Figure 1. Considering that these
models are equally accurate, we can use CE robustness to small model shifts as
a secondary criteria for choosing the best model. To do so, we consider the
multiplicity set $\mathcal{M}_{f_{i},\mathbf{x}}$ of $f_{i}$ containing linear
models with a bound $\delta=2$ on the $l_{\infty}$ norm. In Figure 1, the
yellow regions represent all samples that are not robust to
$\mathcal{M}_{f_{i},\mathbf{x}}$. The optimal robust CE of $\mathbf{x}$ on
$f_{0}$ is $\mathbf{x}^{\prime}_{0}=(0,0)$. However, $\mathbf{x}$ has no
robust CE on $f_{1}$, and its optimal robust CE on $f_{2}$ is
$\mathbf{x}^{\prime}_{2}=(0,0.5)$. In other words, if
$f_{\mathrm{fix}}=f_{1}$, then no CE generation algorithm $g$ will be able to
find a robust counterfactual. But if $f_{\mathrm{fix}}=f_{2}$, rather than
$f_{0}$, we can find a better robust CE (since $\mathbf{x}^{\prime}_{2}$ is
closer to $\mathbf{x}$ than $\mathbf{x}^{\prime}_{0}$ is). So, if we train $f$
and $g$ jointly, we can try to end up with a model closer to $f_{2}$: that is,
a model that performs well and also yields robust CEs that are of high-quality
(i.e., close to the original sample).

##  4 Approach

Our approach, VeriTraCER, is a training algorithm that takes a training set
$D$ as input and outputs a model $f$ and a CE generator $g$ such that the CEs
generated by $g$ have a high rate of $\mathcal{M}_{f,\mathbf{x}}$-robustness.
Devising such a training algorithm requires us to solve the following two
challenges. First, training a model to produce robust explanations requires
reasoning in tandem about the accuracy of the model $f$ and the validity,
quality, and robustness of the CEs produced by the CE generator $g$ (Section
4.1). Second, we need to be able to _soundly overapproximate_ our robust loss
function, as using gradient descent to approximate the loss will be both
inefficient and may overlook some
$f_{\mathrm{m}}\in\mathcal{M}_{f,\mathbf{x}}$ (Section 4.2).

###  4.1 Training the model $f$ and the CE generator $g$ in tandem

To intertwine model training and robust CE generation, we build off of
CounterNet (Guo et al., 2023), which jointly learns a model $f$ and a (non-
robust) CE generator $g$. Intuitively, we have two loss functions,
$\mathcal{L}_{f}$ and $\mathcal{L}_{g}$, that are used to optimize $f$ and
$g$, respectively. In each training epoch, there are three steps. First, we
generate CEs $\mathbf{x}^{\prime}_{i}=g(f,\mathbf{x}_{i})$ according to the
current weights of $f$ and $g$. Next, we optimize the parameters of $f$ using
gradient descent (WRT $f$) on $\mathcal{L}_{f}$, and finally, we optimize the
parameters of $g$ using gradient descent (WRT $g$) on $\mathcal{L}_{g}$.

We denote the output of model $f$ after the final sigmoid layer but before
discretization as $\bar{f}$. The function $\bar{f}(\mathbf{x})\in[0,1]$ is
frequently employed in loss functions for improved optimization compared to
the hard label $f(\mathbf{x})\in\\{0,1\\}$. We define helper loss functions
$\mathcal{L}_{\mathrm{A}}$, $\mathcal{L}_{\mathrm{V}}$, and
$\mathcal{L}_{\mathrm{Q}}$ to promote accuracy, CE validity, and CE quality
(i.e., proximity to the original sample $\mathbf{x}_{i}$), respectively. All
three losses were used in the original CounterNet loss, and can be any binary
loss function, e.g., MSE. $\mathcal{L}_{\mathrm{R}}$ is our proposed RobustCE
loss, and aims to measure the robustness of a CE $\mathbf{x}^{\prime}$ with
respect to $\mathcal{M}_{f,\mathbf{x}}$ (as defined by hyperparameters
$\delta$ and $p$, which are omitted in the equations below for brevity). This
loss function should capture the worst classifier
$f_{\mathrm{m}}\in\mathcal{M}_{f,\mathbf{x}}$, i.e., the model that does most
poorly on the CE $\mathbf{x}^{\prime}$. We define
$\mathcal{L}_{\mathrm{R}}(\mathbf{x},\mathbf{x}^{\prime},\theta_{f})=\max_{f_{\mathrm{m}}\in\mathcal{M}_{f,\mathbf{x}}}\mathcal{L}_{\mathrm{MSE}}(\bar{f_{\mathrm{m}}}(\mathbf{x}^{\prime}),1-y)$.
With those definitions in mind, we define $\mathcal{L}_{f}$ and
$\mathcal{L}_{g}$ as follows:

$\mathcal{L}_{f}(\mathbf{x}_{i},y_{i},\mathbf{x}^{\prime}_{i},\theta_{f})=\lambda_{1}\mathcal{L}_{\mathrm{A}}(\bar{f}(\mathbf{x}_{i}),y_{i})+\lambda_{2}\mathcal{L}_{\mathrm{R}}(\mathbf{x}_{i},\mathbf{x}^{\prime}_{i},\theta_{f})$
(1)
$\mathcal{L}_{g}(\mathbf{x}_{i},y_{i},g(f,\mathbf{x}_{i}),\theta_{f})=\lambda_{3}\mathcal{L}_{\mathrm{Q}}(\mathbf{x}_{i},g(f,\mathbf{x}_{i}))+\lambda_{4}\mathcal{L}_{\mathrm{V}}(\bar{f}(g(f,\mathbf{x}_{i})),1-y_{i})+\lambda_{2}\mathcal{L}_{\mathrm{R}}(\mathbf{x}_{i},g(f,\mathbf{x}_{i}),\theta_{f})$
(2)

###  4.2 Computing the RobustCE loss

Exactly minimizing $\mathcal{L}_{\mathrm{R}}$ is expensive: we either need to
consider infinitely many reasonable models in $\mathcal{M}_{f,\mathbf{x}}$ or
we need to rely on expensive gradient descent to approximate the model
$f_{\mathrm{m}}$ with no guarantee that it is the worst-case model. To address
this, we use abstract interpretation (Cousot & Cousot, 1977) to efficiently
compute an upper bound
$\mathcal{L}_{\mathrm{R}}^{\sharp}(\mathbf{x},\mathbf{x}^{\prime},\theta_{f})$
on the RobustCE loss. This upper bound, when minimized and sufficiently tight,
ensures a reduction in the RobustCE loss simultaneously. The challenge of
overapproximation is how to compute the upper bound
$\mathcal{L}_{\mathrm{R}}^{\sharp}(\mathbf{x},\mathbf{x}^{\prime},\theta_{f})$
as tightly as possible: overly loose bounds may hamper prediction accuracy, CE
validity, and CE quality, without meaningfully enhancing CE robustness.

In Section 4.2.1, we show how to compute $\mathcal{L}_{\mathrm{R}}^{\sharp}$
using two existing abstract interpretation techniques, IBP (Gowal et al.,
2018) and CROWN-IBP (Zhang et al., 2020). Then, in Section 4.2.2 we introduce
a new technique that yields tighter upper bounds on the robust loss while
maintaining efficiency similar to CROWN-IBP.

####  4.2.1 Overapproximating the RobustCE loss using existing techniques

Interval bound propagation (IBP) allows us to evaluate a function on an
infinite set of inputs represented as a hyperrectangle in $\mathbb{R}^{n}$. We
will use an interval $\mathbf{z}^{\sharp}=[\mathbf{z}^{L},\mathbf{z}^{U}]$,
where $\mathbf{z}^{L},\mathbf{z}^{U}\in\mathbb{R}^{n}$ and $\forall 1\leq
i\leq n,\ \mathbf{z}^{L}_{i}\leq\mathbf{z}^{U}_{i}$, to denote the set of all
$n$-dimensional vectors whose $i$-th element is between $\mathbf{z}^{L}_{i}$
and $\mathbf{z}^{U}_{i}$, inclusive. Previous work (Gowal et al., 2018; Zhang
et al., 2021) has used IBP to overapproximate the worst-case loss of test-time
adversarial attacks. These works look at perturbations of an input
$\mathbf{x}$, i.e., an infinite set of possible $\mathbf{x}$, as evaluated on
a fixed model. By contrast, we want to use IBP to overapproximate the
infinitely many models in $\mathcal{M}_{f,\mathbf{x}}$.

To use IBP, we will relax the definition of $\mathcal{M}_{f,\mathbf{x}}$ to
$\mathcal{M}_{f}$ by removing the requirement that the classification of
$\mathbf{x}$ remains constant, i.e.,
$\mathcal{M}_{f}=\\{f_{\mathrm{m}}\mid\|\theta_{f}-\theta_{{f_{\mathrm{m}}}}\|_{p}\leq\delta\\}$.
After relaxation, the components of $\theta_{f}$ (the weights and biases in
linear layers) can be overapproximated by intervals. Therefore, we apply IBP
to overapproximate the RobustCE loss under the multiplicity $\mathcal{M}_{f}$
using interval arithmetic.

Figure 2: Our approach Simul-CROWN achieves a tighter overapproximation than
IBP and CROWN-IBP because the latter techniques include portions of the red
region where the CE is not robust.

CROWN-IBP (Zhang et al., 2020) achieves tighter bounds than IBP by
incorporating a backward tightening technique employed in $\alpha$-CROWN
(Zhang et al., 2018; Singh et al., 2019). In the evaluation of a function
$f_{\theta^{\sharp}}$ whose parameters $\theta^{\sharp}$ are intervals, CROWN-
IBP finds
$\boldsymbol{\alpha}^{l},\boldsymbol{\alpha}^{u},\beta^{l},\beta^{u}$ such
that
$\sigma(\boldsymbol{\alpha}^{l}\theta+\beta^{l})\leq\bar{f}_{\theta}(\mathbf{x})\leq\sigma(\boldsymbol{\alpha}^{u}\theta+\beta^{u})$
for all $\theta\in\theta^{\sharp}$. Then, these equations
($\boldsymbol{\alpha}^{l}\theta+\beta^{l}$ and
$\boldsymbol{\alpha}^{u}\theta+\beta^{u}$) serve to bound $\bar{f}$, allowing
us to do IBP with tighter – yet still sound – bounds. Similar to the adaption
of IBP, we need to relax the definition of $\mathcal{M}_{f,\mathbf{x}}$ to
$\mathcal{M}_{f}$ for CROWN-IBP because it cannot deal with the additional
requirement of prediction robustness that $\mathcal{M}_{f,\mathbf{x}}$
requires. After this relaxation, the parameters in all layers (i.e., the
weight matrices and bias vectors), can be overapproximated into intervals.
When applying CROWN-IBP to overapproximate the RobustCE loss, the interval
$\theta^{\sharp}$ becomes the intervals of parameters in all layers, and the
$\boldsymbol{\alpha}^{l},\boldsymbol{\alpha}^{u},\beta^{l},\beta^{u}$ are
computed based on a specific input $\mathbf{x}$ and the lower and upper bounds
of parameters, i.e., $\theta^{L}$ and $\theta^{U}$.

######  Theorem 4.1 (Soundness and Tightness).

For any $\mathbf{x}$, $\mathbf{x}^{\prime}$, and $f$, the CROWN-IBP-
overapproximated loss
$\mathcal{L}_{\mathrm{R}}^{\sharp\mathrm{CIBP}}(\mathbf{x},\mathbf{x}^{\prime},\theta_{f})$
is an upper bound of the RobustCE loss
$\mathcal{L}_{\mathrm{R}}(\mathbf{x},\mathbf{x}^{\prime},\theta_{f})$ and a
lower bound of the IBP-overapproximated loss
$\mathcal{L}_{\mathrm{R}}^{\sharp\mathrm{IBP}}(\mathbf{x},\mathbf{x}^{\prime},\theta_{f})$.
Formally,

$\mathcal{L}_{\mathrm{R}}^{\sharp\mathrm{IBP}}(\mathbf{x},\mathbf{x}^{\prime},\theta_{f})\geq\mathcal{L}_{\mathrm{R}}^{\sharp\mathrm{CIBP}}(\mathbf{x},\mathbf{x}^{\prime},\theta_{f})\geq\mathcal{L}_{\mathrm{R}}(\mathbf{x},\mathbf{x}^{\prime},\theta_{f})$

######  Example 4.1.

Consider an interval linear layer with weight matrix
$W^{\sharp}=([-1,3],[-3,1])$ and bias scalar $b^{\sharp}=[-4,0]$. For an input
$\mathbf{x}=(-4,-1)$, CROWN-IBP concatenates $W^{\sharp}$ and $b^{\sharp}$ to
form $\theta^{\sharp}=([-1,3],[-3,1],[-4,0])$ and computes
$\boldsymbol{\alpha}^{l}=\boldsymbol{\alpha}^{u}=(\mathbf{x};1)=(-4,1,1)$
(note that we append 1 to $\mathbf{x}$ as the multiplicative factor for
$b^{\sharp}$) and $\beta^{l}=\beta^{u}=0$. Then we can symbolically compute
the output based on $\boldsymbol{\alpha}^{l}{}^{T}\theta+\beta^{l}$ (or
$\boldsymbol{\alpha}^{u}{}^{T}\theta+\beta^{u}$ because the two are the same),
i.e., $(-4,-1,1)^{\top}([-1,3],[-3,1],[-4,0])$. To compute the concrete upper
bounds, we apply interval arithmetic to the vector product and get $[-17,7]$.
Note that IBP and CROWN-IBP yield identical lower and upper bounds for this
linear model because it has only one layer.

However, the bounds computed by CROWN-IBP are still loose due to the
relaxation from $\mathcal{M}_{f,\mathbf{x}}$ to $\mathcal{M}_{f}$ as
illustrated in the following example.

######  Example 4.2.

Consider the same interval linear model as in Example 4.1. Note that the
midpoint of this interval network corresponds to our “true” $f$, i.e.,
$W=(1,-1)$ and $b=-2$. Let $\bar{f}(\mathbf{x})=\sigma(W\mathbf{x}+b)$. For
the input $\mathbf{x}=(4,1)$ and its CE $\mathbf{x}^{\prime}=(-4,-1)$, we have
$\bar{f}(\mathbf{x})>0.5$ and $\bar{f}(\mathbf{x}^{\prime})<0.5$. IBP and
CROWN-IBP are unable to prove that the CE $\mathbf{x}^{\prime}$ is robust even
though $\mathcal{M}_{f,\mathbf{x}}$ does not intersect with the “unsafe” red
region shown in Figure 2. To see this limitation, let us consider a model
$f_{\mathrm{m}}$ with $W^{\prime}=(-1,-3)\in W^{\sharp}$ and $b^{\prime}=0\in
b^{\sharp}$ that is in $\mathcal{M}_{f}$ but outside
$\mathcal{M}_{f,\mathbf{x}}$ (the output of $\bar{f_{\mathrm{m}}}(\mathbf{x})$
is $\sigma(-7)<0.5$). $f_{\mathrm{m}}$ has an output of
$\sigma(W^{\prime}\mathbf{x}^{\prime}+b^{\prime})=\sigma(7)>0.5$ on the CE
$\mathbf{x}^{\prime}$ meaning that the CE is invalid on $f_{\mathrm{m}}$. This
limitation arises from the relaxation from $\mathcal{M}_{f,\mathbf{x}}$ to
$\mathcal{M}_{f}$. Our approach (Section 4.2.2) overcomes this limitation by
overapproximating $\mathcal{M}_{f,\mathbf{x}}$ instead of its relaxation
$\mathcal{M}_{f}$ and is able to prove that $\mathbf{x}^{\prime}$ is
$\mathcal{M}_{f,\mathbf{x}}$-robust.

####  4.2.2 Overapproximation the RobustCE Loss using Simul-CROWN

The primary challenge faced by IBP and CROWN-IBP in handling
$\mathcal{M}_{f,\mathbf{x}}$ is their inability to _simultaneously_ reason
about the overapproximation of $f(\mathbf{x})$ and $f(\mathbf{x}^{\prime})$.
In this section, we present our approach, Simul-CROWN, which addresses this
challenge to achieve a more precise overapproximation of the RobustCE loss.

To simultaneously overapproximate $f(\mathbf{x})$ and
$f(\mathbf{x}^{\prime})$, we first reframe the definition of the RobustCE
loss,
$\mathcal{L}_{\mathrm{R}}=\max_{f_{\mathrm{m}}\in\mathcal{M}_{f,\mathbf{x}}}\mathcal{L}_{\mathrm{MSE}}(\bar{f_{\mathrm{m}}}(\mathbf{x}^{\prime}),1-\hat{y})$.
Namely, if $\hat{y}=1$, $\mathcal{L}_{\mathrm{R}}$ is maximized when the
worst-case classifier is
$[\mathrm{arg}\max_{f_{\mathrm{m}}\in\mathcal{M}_{f}}\bar{f_{\mathrm{m}}}(\mathbf{x}^{\prime})\quad\textrm{s.t.
}\bar{f_{\mathrm{m}}}(\mathbf{x})\geq 0]$. However, if $\hat{y}=0$,
$\mathcal{L}_{\mathrm{R}}$ is maximized when the worst-case classifier is
$[\mathrm{arg}\max_{f_{\mathrm{m}}\in\mathcal{M}_{f}}\bar{f_{\mathrm{m}}}(\mathbf{x}^{\prime})\quad\textrm{s.t.
}\bar{f_{\mathrm{m}}}(\mathbf{x})\geq 0]$. Using this insight, we apply CROWN-
IBP to obtain lower and upper bounds in $\mathcal{M}_{f}$ for each training
instance $\mathbf{x}_{i}$. For all parameters $\theta\in\theta^{\sharp}$, we
have
$\sigma(\boldsymbol{\alpha}^{l}_{i}\theta+\beta^{l}_{i})\leq\bar{f}_{\theta}(\mathbf{x})\leq\sigma(\boldsymbol{\alpha}^{u}_{i}\theta+\beta^{u}_{i})$
and
$\sigma(\boldsymbol{\mu}^{l}_{i}\theta+\nu^{l}_{i})\leq\bar{f}_{\theta}(\mathbf{x})\leq\sigma(\boldsymbol{\mu}^{u}_{i}\theta+\nu^{u}_{i})$.
We compute $\boldsymbol{\alpha}_{i}$ and $\beta_{i}$ based on $\mathbf{x}_{i}$
and $\theta^{\sharp}$, while $\boldsymbol{\mu}_{i}$ and $\nu_{i}$ are computed
based on $\mathbf{x}^{\prime}_{i}$ and $\theta^{\sharp}$.

######  Theorem 4.2 (Overapproximation by Simul-CROWN).

Conditional on the value of $\hat{y}$, the optimal value of
$\mathcal{L}_{\mathrm{R}}$ can be upper bounded by the solution to one of the
following cases,

$\displaystyle\begin{cases}t_{\hat{y}=1}=\underset{\theta\in\theta^{\sharp}}{\max}~{}~{}\boldsymbol{\mu}^{u}_{i}\theta+\nu^{u}_{i}\quad\textrm{s.t.
}\boldsymbol{\alpha}^{u}_{i}\theta+\beta^{u}_{i}\geq 0,&\\\
t_{\hat{y}=0}=\underset{\theta\in\theta^{\sharp}}{\min}~{}~{}~{}\boldsymbol{\mu}^{l}_{i}\theta+\nu^{l}_{i}~{}\quad\textrm{s.t.
}\boldsymbol{\alpha}^{l}_{i}\theta+\beta^{l}_{i}\leq 0,&\end{cases}$ (3)

That is, if $\hat{y}=1$, we have
$\mathcal{L}_{\mathrm{R}}\leq\mathcal{L}_{\mathrm{MSE}}(\sigma(t_{\hat{y}=1}),0)$,
and if $\hat{y}=0$,
$\mathcal{L}_{\mathrm{R}}\leq\mathcal{L}_{\mathrm{MSE}}(\sigma(t_{\hat{y}=0}),1)$.

We can use a solver to address each case in LABEL:eq:sim_overapprox by
encoding it into a linear programming (LP) problem. However, solving an LP for
each training batch is time-consuming and breaks the gradient information
required to optimize the RobustCE loss. Appendix A introduces a greedy
algorithm (Algorithm 2) that solves LABEL:eq:sim_overapprox in $O(n\log n)$
time, where $n$ is the size of parameters, and can be implemented in PyTorch
preserving the gradient information.

######  Theorem 4.3 (Soundness and Tightness).

For any $\mathbf{x}$, $\mathbf{x}^{\prime}$, and $f$, the overapproximated
loss
$\mathcal{L}_{\mathrm{R}}^{\sharp\mathrm{Ours}}(\mathbf{x},\mathbf{x}^{\prime},\theta_{f})$
is an upper bound of the RobustCE loss
$\mathcal{L}_{\mathrm{R}}(\mathbf{x},\mathbf{x}^{\prime},\theta_{f})$ and a
lower bound of the CROWN-IBP-overapproximated loss
$\mathcal{L}_{\mathrm{R}}^{\sharp\mathrm{CIBP}}(\mathbf{x},\mathbf{x}^{\prime},\theta_{f})$.
Formally,

$\mathcal{L}_{\mathrm{R}}^{\sharp\mathrm{CIBP}}(\mathbf{x},\mathbf{x}^{\prime},\theta_{f})\geq\mathcal{L}_{\mathrm{R}}^{\sharp\mathrm{Ours}}(\mathbf{x},\mathbf{x}^{\prime},\theta_{f})\geq\mathcal{L}_{\mathrm{R}}(\mathbf{x},\mathbf{x}^{\prime},\theta_{f})$

Note that
$\mathcal{L}_{\mathrm{R}}^{\sharp\mathrm{Ours}}(\mathbf{x},\mathbf{x}^{\prime},\theta_{f})=\mathcal{L}_{\mathrm{R}}^{\sharp\mathrm{CIBP}}(\mathbf{x},\mathbf{x}^{\prime},\theta_{f})$
when for all $1\leq i\leq n$,
$\boldsymbol{\mu}^{u}_{i}\boldsymbol{\alpha}^{u}_{i}>0$ if $\hat{y}=1$ (or
$\boldsymbol{\mu}^{l}_{i}\boldsymbol{\alpha}^{l}_{i}>0$ if $\hat{y}=0$).

The following example shows that Simul-CROWN is tighter than the CROWN-IBP-
overapproximated loss.

######  Example 4.3.

Consider the same linear model given in Example 4.2, for the original input
$\mathbf{x}=(4,1)$ and the corresponding CE $\mathbf{x}^{\prime}=(-4,-1)$, we
have the following optimization problem according to LABEL:eq:sim_overapprox.

$\displaystyle\underset{\theta\in\theta^{\sharp}}{\max}-4z_{1}-z_{2}+z_{3}\quad\textrm{s.t.
}4z_{1}+z_{2}+z_{3}\geq 0,\text{ where
}\theta^{L}=[-1,-3,-4],\theta^{U}=[3,1,0]$

According to Algorithm 2, the optimal value $0$ is obtained when
$\theta=[-\frac{1}{4},1,0]$. In other words, the output of the worst-case
classifier in $\mathcal{M}_{f,\mathbf{x}}$ is $\sigma(0)$ computed by Simul-
CROWN. This upper bound is tighter than the ones computed by IBP and CROWN-IBP
($\sigma(7)$ in Example 4.2).

##  5 Experimental Evaluation

In this section, we present our evaluation of VeriTraCER. First, we describe
the experiment setup (Section 5.1), and next we show results (Section 5.2) on
$\delta$-robustness, cross-model validity, and CE quality.

###  5.1 Experimental setup

##### Datasets

We perform our evaluation on two _real-world distribution shift_ datasets,
Cardiotocography (CTG) (Campos & Bernardes, 2010) and WHO (Rajarshi, 2017), as
well as three other datasets commonly used in the robust counterfactual
literature, Home Equity Line of Credit (HELOC) (FICO, 2018), Taiwanese Credit
(TC) (Yeh & Lien, 2009), and Open University Learning Analytics (OULA)
(Kuzilek et al., 2017). CTG aims to predict whether a fetal cardiotocogram is
healthy, WHO aims to predict whether a country’s life expectancy is above the
median, HELOC aims to predict the likelihood that home buyers will repay a
loan within 2 years, TC aims to predict default on their credit card payments,
and OULA aims to predict whether a student will pass an online class.
Additional details including dataset composition and preprocessing steps such
as binarization and distribution shift construction are available in the
appendix.

##### Metrics

Given an $\mathbf{x}$ and a corresponding CE $\mathbf{x}^{\prime}$,
$\delta$-robustness is satisfied when we can certify that
$\mathbf{x}^{\prime}$ obeys
$f_{\mathrm{m}}(\mathbf{x}^{\prime})=f(\mathbf{x}^{\prime})$ for all
$f_{\mathrm{m}}\in\mathcal{M}_{f,\mathbf{x}}=\\{f_{\mathrm{m}}\mid
f(\mathbf{x})=f_{\mathrm{m}}(\mathbf{x})\land\|\theta_{f_{\mathrm{m}}}-\theta_{f}\|_{p}\leq\delta\\}$
using Simul-CROWN. We dynamically determine distinct $\delta_{i}$ values for
each layer (see the appendix for details); for simplicity, we will continue to
write $\delta$-robustness as if $\delta$ is a constant. By contrast,
$\mathbf{x}^{\prime}$ exhibits cross-model validity with respect to two models
$f$ and $f_{\mathrm{m}}$ when
$f(\mathbf{x}^{\prime})=f_{\mathrm{m}}(\mathbf{x}^{\prime})$. In practice, we
compute empirical cross-model validity across sets of models trained with one
of three variations: different random initialization (RI) for training $f$ and
$f_{\mathrm{m}}$, different training datasets constructed by randomly removing
1% of the data (LOO) for $f$ and $f_{\mathrm{m}}$, or different training
datasets under data shift (DS) for $f$ and $f_{\mathrm{m}}$. We refer to the
_$\delta$ -robustness rate_ as the fraction of data points in a test set that
exhibit $\delta$-robustness, and similarly for _cross-model validity rate_. As
a secondary metric, we consider how feasible the counterfactual is to
implement through its proximity, i.e., the $l_{1}$ distance between
$\mathbf{x}$ and $\mathbf{x}^{\prime}$. We consider two other distance metrics
in the appendix.

##### Baseline comparisons

We compare our approach with the following existing robust counterfactual
generation methods: Counternet (CN) (Guo et al., 2023) without our
robustifying modifications; ROAR (Upadhyay et al., 2021), which finds robust
counterfactuals by using adversarial training in the counterfactual generation
process; and SNS (Black et al., 2021), which finds CEs in regions with a low
Lipschitz constant. For ROAR and SNS, we train NNs with the same architecture
as the predictor part of CounterNet; we optimize this model using a loss
function that solely prioritizes accuracy. We implement all three techniques
(IBP, CROWN-IBP, Simul-CROWN) from Section 4.2. We include results for Simul-
CROWN here and put results for the others in the ablation study in the
appendix.

##### Experimental procedure

To evaluate $\delta$-robustness, we select layer-specific $\delta_{i}$ (see
appendix for details) and perform training and all evaluation for that
setting. To evaluate cross-model validity for RI and LOO, we train 10 models,
generate CEs for each, and report the average fraction of these CEs that
remain valid across the other 9 models. To evaluate cross-model validity for
DS, we train a model on the original data, then finetune for a small number
(typically 20) additional epochs on the new data to obtain the shifted model.

Table 1: Fraction of samples that are pair-wise cross-model validity across 10 model trained with different random initializations (RI) or different segments of data randomly removed (leave-one-out, or LOO). Standard deviations are in parenthesis. Best result is in bold and second-best result is underlined. | Random Initialization | Leave-One-Out
---|---|---
| HELOC | TC | OULA | HELOC | TC | OULA
VeriTraCER |  0.93 (0.03) |  0.88 (0.14) |  0.94 (0.05) | 0.86 (0.14) |  0.97 (0.01) |  0.96 (0.05)
CN |  0.92 (0.07) | 0.82 (0.26) | 0.92 (0.07) |  0.91 (0.07) | 0.96 (0.02) |  0.97 (0.01)
ROAR | 0.88 (0.02) | 0.44 (0.25) |  0.93 (0.05) | 0.88 (0.02) | 0.88 (0.01) | 0.95 (0.04)
SNS | 0.90 (0.04) |  0.98 (0.02) | 0.80 (0.11) |  0.94 (0.05) |  1.00 (0.00) | 0.90 (0.06)
Table 2: Fraction of samples whose CEs are valid and robust after finetuning with a distribution shift. Standard deviation over 10 trials in parentheses. | CTG | WHO
---|---|---
VeriTraCER |  0.987 (0.010) |  0.995 (0.027)
CN |  0.981 (0.019) |  0.967 (0.030)
ROAR | 0.570 (0.104) | 0.884 (0.078)
SNS | 0.407 (0.024) | 0.846 (0.050)

###  5.2 Results

#####  $\delta$-robustness

VeriTraCER exhibits high levels of $\delta$-robustness: 70.24% (±17.41%) of
test samples exhibit $\delta$-robustness for HELOC, 81.08% (±4.51%) for TC,
and 96.96% (±0.96%) for OULA. For the same value of $\delta$, CN has a
$\delta$-robustness for 26.29% (±18.29%) of HELOC samples, but 0% for TC and
OULA samples. Likewise, the CEs generated by ROAR and SNS for models trained
in a standard way have 0% $\delta$-robustness on all three datasets.

##### Random initialization and leave-one-out results

Table 1 shows what fraction of counterfactuals, generated for a particular
model, remain valid for a model trained (a) using a different random seed, or
(b) with a different 1% of the training data removed. We note that most CE and
dataset combinations have high ($>90$%) cross-model validity rates; however,
VeriTraCER has the highest or second-highest cross-model validity rate for
most settings.

##### Real-world distribution shifts

Table 2 shows the robustness of CEs to real-world distribution shifts. We see
that VeriTraCER significantly outperforms the non-CounterNet baseline methods,
especially on CTG.

##### Counterfactual Quality

Table 3: Average proximity for valid CEs. (*) indicates that the CE technique does not correctly account for categorical features. | CTG | WHO | HELOC | TC | OULA
---|---|---|---|---|---
VeriTraCER | 0.250 | 0.216 | 0.099 | 0.244 | 0.179
CN | 0.220 | 0.164 | 0.107 | 0.234 | 0.169
ROAR | 0.189 | 0.031 | 0.032 |  0.018* |  0.015*
SNS | 0.038 | 0.031 | 0.041 |  0.015* |  0.013*

A weakness of VeriTraCER is that the generated CEs are a larger distance from
the original CEs. The data on proximity (the $l_{1}$-distance between an input
$\mathbf{x}$ and its CE $\mathbf{x}^{\prime}$ is) is summarized in Table 3.
Note that the two datasets where VeriTraCER performs worst relative to the
other techniques (TC and OULA) have categorical features, which ROAR does not
modify, and SNS does not handle properly (i.e., by breaking one-hot
encodings).

##  6 Conclusions

We have presented VeriTraCER, a training algorithm that jointly produces a
model $f$ and a CE generator $g$ such that the CEs generated by $g$ will be
robust to small changes in the weights of $f$. We do this by minimizing an
upper bound on our robust loss, i.e., we minimize the loss on the _worst-case_
model in the multiplicity set $\mathcal{M}_{f,\mathbf{x}}$. We provide a
refinement of interval-bound propagation, Simul-CROWN, allowing the over-
approximation to be tighter than other state-of-the-art approaches. Our
approach is able to to find CEs that are _certifiably robust_ at high rates
(typically over 90%). This carries over to high robustness for empirical model
updates, such as retraining with a different random seed. In particular, we
outperform state-of-the-art approaches on finetuning after real-world
distribution shifts. The tradeoff is that VeriTraCER generally yields CEs that
are a larger distance form the original sample. In some settings, the
additional robustness as well as fast CE generation time may be worth the
additional recourse costs, but future work should aim to provide similar
deterministic guarantees and empirical performance with smaller CE distances.

## Reproducibility Statement

Our code is available at https://github.com/ForeverZyh/robust_cfx.

## References

  * Anders et al. (2020) Christopher Anders, Plamen Pasliev, Ann-Kathrin Dombrowski, Klaus-Robert Müller, and Pan Kessel.  Fairwashing explanations with off-manifold detergent.  In Hal Daumé III and Aarti Singh (eds.), _Proceedings of the 37th International Conference on Machine Learning_ , volume 119 of _Proceedings of Machine Learning Research_ , pp. 314–323. PMLR, 13–18 Jul 2020.  URL https://proceedings.mlr.press/v119/anders20a.html.
  * Black et al. (2021) Emily Black, Zifan Wang, Matt Fredrikson, and Anupam Datta.  Consistent counterfactuals for deep models.  _arXiv preprint arXiv:2110.03109_ , 2021.
  * Bui et al. (2022) Ngoc Bui, Duy Nguyen, and Viet Anh Nguyen.  Counterfactual Plans under Distributional Ambiguity.  In _International Conference on Learning Representations_ , 2022.  URL https://openreview.net/forum?id=noaG7SrPVK0.
  * Campos & Bernardes (2010) D. Campos and J. Bernardes.  Cardiotocography.  UCI Machine Learning Repository, 2010.  DOI: https://doi.org/10.24432/C51S4N.
  * Consumer Financial Protection Bureau (2018) Consumer Financial Protection Bureau.  12 CFR Part 1002 - Equal Credit Opportunity Act (Regulation B).  Technical report, Consumer Financial Protection Bureau, January 2018.  URL https://www.consumerfinance.gov/rules-policy/regulations/1002/.
  * Cousot & Cousot (1977) Patrick Cousot and Radhia Cousot.  Abstract interpretation: A unified lattice model for static analysis of programs by construction or approximation of fixpoints.  In _Conference Record of the 4th ACM Symposium on Principles of Programming Languages, POPL 1977, Los Angeles, California, USA, January 1977_ , pp. 238–252, 1977.  doi: 10.1145/512950.512973.  URL https://doi.org/10.1145/512950.512973.
  * D’Amour et al. (2022) Alexander D’Amour, Katherine Heller, Dan Moldovan, Ben Adlam, Babak Alipanahi, Alex Beutel, Christina Chen, Jonathan Deaton, Jacob Eisenstein, Matthew D. Hoffman, Farhad Hormozdiari, Neil Houlsby, Shaobo Hou, Ghassen Jerfel, Alan Karthikesalingam, Mario Lucic, Yian Ma, Cory McLean, Diana Mincu, Akinori Mitani, Andrea Montanari, Zachary Nado, Vivek Natarajan, Christopher Nielson, Thomas F. Osborne, Rajiv Raman, Kim Ramasamy, Rory Sayres, Jessica Schrouff, Martin Seneviratne, Shannon Sequeira, Harini Suresh, Victor Veitch, Max Vladymyrov, Xuezhi Wang, Kellie Webster, Steve Yadlowsky, Taedong Yun, Xiaohua Zhai, and D. Sculley.  Underspecification presents challenges for credibility in modern machine learning.  _J. Mach. Learn. Res._ , 23(1), jan 2022.  ISSN 1532-4435.
  * Dombrowski et al. (2022) Ann-Kathrin Dombrowski, Christopher J Anders, Klaus-Robert Müller, and Pan Kessel.  Towards robust explanations for deep neural networks.  _Pattern Recognition_ , 121:108194, 2022.
  * Dutta et al. (2022) Sanghamitra Dutta, Jason Long, Saumitra Mishra, Cecilia Tilli, and Daniele Magazzeni.  Robust Counterfactual Explanations for Tree-Based Ensembles.  In _International Conference on Machine Learning_ , pp. 5742–5756. PMLR, 2022.
  * European Commission (2018) European Commission.  General Data Protection Regulation (GDPR), 2018.  URL https://gdpr.eu/tag/gdpr/.
  * FICO (2018) FICO.  Explainable Machine Learning Challenge, 2018.  URL https://community.fico.com/s/explainable-machine-learning-challenge?tabset-158d9=d157e.
  * Forel et al. (2022) Alexandre Forel, Axel Parmentier, and Thibaut Vidal.  Robust Counterfactual Explanations for Random Forests.  _arXiv preprint arXiv:2205.14116_ , 2022.
  * Gowal et al. (2018) Sven Gowal, Krishnamurthy Dvijotham, Robert Stanforth, Rudy Bunel, Chongli Qin, Jonathan Uesato, Relja Arandjelovic, Timothy A. Mann, and Pushmeet Kohli.  On the effectiveness of interval bound propagation for training verifiably robust models.  _CoRR_ , abs/1810.12715, 2018.  URL http://arxiv.org/abs/1810.12715.
  * Guo et al. (2023) Hangzhi Guo, Thanh Hong Nguyen, and Amulya Yadav.  Counternet: End-to-end training of prediction aware counterfactual explanations.  In Ambuj Singh, Yizhou Sun, Leman Akoglu, Dimitrios Gunopulos, Xifeng Yan, Ravi Kumar, Fatma Ozcan, and Jieping Ye (eds.), _Proceedings of the 29th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, KDD 2023, Long Beach, CA, USA, August 6-10, 2023_ , pp. 577–589. ACM, 2023.  doi: 10.1145/3580305.3599290.  URL https://doi.org/10.1145/3580305.3599290.
  * Hamman et al. (2023) Faisal Hamman, Erfaun Noorani, Saumitra Mishra, Daniele Magazzeni, and Sanghamitra Dutta.  Robust counterfactual explanations for neural networks with probabilistic guarantees.  In Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett (eds.), _Proceedings of the 40th International Conference on Machine Learning_ , volume 202 of _Proceedings of Machine Learning Research_ , pp. 12351–12367. PMLR, 23–29 Jul 2023.  URL https://proceedings.mlr.press/v202/hamman23a.html.
  * Heo et al. (2019) Juyeon Heo, Sunghwan Joo, and Taesup Moon.  _Fooling Neural Network Interpretations via Adversarial Model Manipulation_.  Curran Associates Inc., Red Hook, NY, USA, 2019.
  * Jiang et al. (2023) Junqi Jiang, Francesco Leofante, Antonio Rago, and Francesca Toni.  Formalising the robustness of counterfactual explanations for neural networks.  In Brian Williams, Yiling Chen, and Jennifer Neville (eds.), _Thirty-Seventh AAAI Conference on Artificial Intelligence, AAAI 2023, Thirty-Fifth Conference on Innovative Applications of Artificial Intelligence, IAAI 2023, Thirteenth Symposium on Educational Advances in Artificial Intelligence, EAAI 2023, Washington, DC, USA, February 7-14, 2023_ , pp. 14901–14909. AAAI Press, 2023.  doi: 10.1609/aaai.v37i12.26740.  URL https://doi.org/10.1609/aaai.v37i12.26740.
  * Karimi et al. (2020) Amir-Hossein Karimi, Gilles Barthe, Borja Balle, and Isabel Valera.  Model-agnostic counterfactual explanations for consequential decisions.  In Silvia Chiappa and Roberto Calandra (eds.), _Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics_ , volume 108 of _Proceedings of Machine Learning Research_ , pp. 895–905. PMLR, 26–28 Aug 2020.  URL https://proceedings.mlr.press/v108/karimi20a.html.
  * Kuzilek et al. (2017) Jakub Kuzilek, Martin Hlosta, and Zdenek Zdrahal.  Open university learning analytics dataset.  _Scientific Data_ , 4:170171, 2017.  doi: 10.1038/sdata.2017.171.
  * Looveren & Klaise (2019) Arnaud Van Looveren and Janis Klaise.  Interpretable Counterfactual Explanations Guided by Prototypes.  _ArXiv_ , abs/1907.02584, 2019.
  * Lundberg & Lee (2017) Scott M Lundberg and Su-In Lee.  A Unified Approach to Interpreting Model Predictions.  In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (eds.), _Advances in Neural Information Processing Systems_ , volume 30. Curran Associates, Inc., 2017.  URL https://proceedings.neurips.cc/paper/2017/file/8a20a8621978632d76c43dfd28b67767-Paper.pdf.
  * Marx et al. (2020) Charles Marx, Flavio Calmon, and Berk Ustun.  Predictive multiplicity in classification.  In _International Conference on Machine Learning_ , pp. 6765–6774. PMLR, 2020.
  * Meyer et al. (2023) Anna P. Meyer, Dan Ley, Suraj Srinivas, and Himabindu Lakkaraju.  On Minimizing the Impact of Dataset Shifts on Actionable Explanations.  In _The 39th Conference on Uncertainty in Artificial Intelligence_ , 2023.  URL https://openreview.net/forum?id=mtd904kJUs.
  * Nguyen et al. (2023) Duy Nguyen, Ngoc Bui, and Viet Anh Nguyen.  Distributionally robust recourse action.  In _The Eleventh International Conference on Learning Representations_ , 2023.  URL https://openreview.net/forum?id=E3ip6qBLF7.
  * Pawelczyk et al. (2020) Martin Pawelczyk, Klaus Broelemann, and Gjergji Kasneci.  On counterfactual explanations under predictive multiplicity.  In _Conference on Uncertainty in Artificial Intelligence_ , pp. 809–818. PMLR, 2020.
  * Poyiadzi et al. (2020) Rafael Poyiadzi, Kacper Sokol, Raul Santos-Rodriguez, Tijl De Bie, and Peter Flach.  FACE: Feasible and Actionable Counterfactual Explanations.  In _Proceedings of the AAAI/ACM Conference on AI, Ethics, and Society_ , AIES ’20, pp. 344–350, New York, NY, USA, 2020. Association for Computing Machinery.  ISBN 9781450371100.  doi: 10.1145/3375627.3375850.  URL https://doi.org/10.1145/3375627.3375850.
  * Rajarshi (2017) Kumar Rajarshi.  Life Expectancy (WHO), 2017.  URL https://www.kaggle.com/datasets/kumarajarshi/life-expectancy-who.
  * Ribeiro et al. (2016) Marco Túlio Ribeiro, Sameer Singh, and Carlos Guestrin.  ”Why Should I Trust You?”: Explaining the Predictions of Any Classifier.  In Balaji Krishnapuram, Mohak Shah, Alexander J. Smola, Charu C. Aggarwal, Dou Shen, and Rajeev Rastogi (eds.), _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, San Francisco, CA, USA, August 13-17, 2016_ , pp. 1135–1144. ACM, 2016.  ISBN 978-1-4503-4232-2.  doi: 10.1145/2939672.2939778.
  * Simonyan et al. (2013) Karen Simonyan, Andrea Vedaldi, and Andrew Zisserman.  Deep Inside Convolutional Networks: Visualising Image Classification Models and Saliency Maps.  _CoRR_ , abs/1312.6034, 2013.
  * Singh et al. (2019) Gagandeep Singh, Timon Gehr, Markus Püschel, and Martin T. Vechev.  An abstract domain for certifying neural networks.  _Proc. ACM Program. Lang._ , 3(POPL):41:1–41:30, 2019.  doi: 10.1145/3290354.  URL https://doi.org/10.1145/3290354.
  * Slack et al. (2020) Dylan Slack, Sophie Hilgard, Emily Jia, Sameer Singh, and Himabindu Lakkaraju.  Fooling LIME and SHAP: Adversarial attacks on post hoc explanation methods.  In _Proceedings of the AAAI/ACM Conference on AI, Ethics, and Society_ , AIES ’20, pp. 180–186, New York, NY, USA, 2020. Association for Computing Machinery.  ISBN 9781450371100.  doi: 10.1145/3375627.3375830.  URL https://doi.org/10.1145/3375627.3375830.
  * Smilkov et al. (2017) Daniel Smilkov, Nikhil Thorat, Been Kim, Fernanda Viégas, and Martin Wattenberg.  Smoothgrad: removing noise by adding noise.  _arXiv preprint arXiv:1706.03825_ , 2017.
  * Srinivas et al. (2022) Suraj Srinivas, Kyle Matoba, Himabindu Lakkaraju, and François Fleuret.  Efficient Training of Low-Curvature Neural Networks.  In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho (eds.), _Advances in Neural Information Processing Systems_ , 2022.  URL https://openreview.net/forum?id=2B2xIJ299rx.
  * Upadhyay et al. (2021) Sohini Upadhyay, Shalmali Joshi, and Himabindu Lakkaraju.  Towards robust and reliable algorithmic recourse.  In Marc’Aurelio Ranzato, Alina Beygelzimer, Yann N. Dauphin, Percy Liang, and Jennifer Wortman Vaughan (eds.), _Advances in Neural Information Processing Systems 34: Annual Conference on Neural Information Processing Systems 2021, NeurIPS 2021, December 6-14, 2021, virtual_ , pp. 16926–16937, 2021.  URL https://proceedings.neurips.cc/paper/2021/hash/8ccfb1140664a5fa63177fb6e07352f0-Abstract.html.
  * Ustun et al. (2018) Berk Ustun, Alexander Spangher, and Yang Liu.  Actionable recourse in linear classification.  _arXiv preprint arXiv:1809.06514_ , 2018.
  * Wachter et al. (2018) Sandra Wachter, Brent Mittelstadt, and Chris Russell.  Counterfactual Explanations Without Opening the Black Box: Automated Decisions and the GDPR.  _Harvard Journal of Law & Technology_, 31:841–887, 04 2018.  doi: 10.2139/ssrn.3063289.
  * Yeh & Lien (2009) I-Cheng Yeh and Che-hui Lien.  The comparisons of data mining techniques for the predictive accuracy of probability of default of credit card clients.  _Expert Systems with Applications_ , 36(2, Part 1):2473–2480, 2009.  ISSN 0957-4174.  doi: https://doi.org/10.1016/j.eswa.2007.12.020.  URL https://www.sciencedirect.com/science/article/pii/S0957417407006719.
  * Zhang et al. (2018) Huan Zhang, Tsui-Wei Weng, Pin-Yu Chen, Cho-Jui Hsieh, and Luca Daniel.  Efficient neural network robustness certification with general activation functions.  In Samy Bengio, Hanna M. Wallach, Hugo Larochelle, Kristen Grauman, Nicolò Cesa-Bianchi, and Roman Garnett (eds.), _Advances in Neural Information Processing Systems 31: Annual Conference on Neural Information Processing Systems 2018, NeurIPS 2018, December 3-8, 2018, Montréal, Canada_ , pp. 4944–4953, 2018.  URL https://proceedings.neurips.cc/paper/2018/hash/d04863f100d59b3eb688a11f95b0ae60-Abstract.html.
  * Zhang et al. (2020) Huan Zhang, Hongge Chen, Chaowei Xiao, Sven Gowal, Robert Stanforth, Bo Li, Duane S. Boning, and Cho-Jui Hsieh.  Towards stable and efficient training of verifiably robust neural networks.  In _8th International Conference on Learning Representations, ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020_. OpenReview.net, 2020.  URL https://openreview.net/forum?id=Skxuk1rFwB.
  * Zhang et al. (2021) Yuhao Zhang, Aws Albarghouthi, and Loris D’Antoni.  Certified robustness to programmable transformations in lstms.  In Marie-Francine Moens, Xuanjing Huang, Lucia Specia, and Scott Wen-tau Yih (eds.), _Proceedings of the 2021 Conference on Empirical Methods in Natural Language Processing, EMNLP 2021, Virtual Event / Punta Cana, Dominican Republic, 7-11 November, 2021_ , pp. 1068–1083. Association for Computational Linguistics, 2021.  doi: 10.18653/V1/2021.EMNLP-MAIN.82.  URL https://doi.org/10.18653/v1/2021.emnlp-main.82.

##  Appendix A Algorithms

Algorithm 1 gives an overview of our robust training approach. We concurrently
train a model $f$ and a CE generator $g$. Alternatively, it is possible to use
a fixed CE generator and to train $f$ to be a model whose CEs – with respect
to the fixed $g$ – are likely to be robust. In that case, line 1 of the
algorithm only initializes $f$, and line 5 of the algorithm is skipped.

Algorithm 1 Joint Training of $f$ and $g$

0: Training set $D={(\mathbf{x}_{i},y_{i})}_{i=1}^{n}$, Hyper-parameters
$\delta,p$ in Definition 3.1

0: Classifier $f$ and robust CE generator $g$

1: Initialize $f$ and $g$

2: for $e=1$ to $\mathrm{maxepoch}$ do

3:
$\mathbf{x}^{\prime}=g(f,\mathbf{x}_{1}),g(f,\mathbf{x}_{2}),\ldots,g(f,\mathbf{x}_{n})$

4: Optimize $\theta_{f}$ via
$\nabla_{\theta_{f}}\sum_{i=1}^{n}\mathcal{L}_{f}(\mathbf{x}_{i},y_{i},\mathbf{x}^{\prime}_{i},\theta_{f},\delta,p)$

5: Optimize $\theta_{g}$ via
$\nabla_{\theta_{g}}\sum_{i=1}^{n}\mathcal{L}_{g}(\mathbf{x}_{i},y_{i},g(f,\mathbf{x}_{i}),\theta_{f},\delta,p)$

6: end for

7: return $f,g$

Algorithm 2 shows our algorithm to solve LABEL:eq:sim_overapprox. Notably,
this algorithm – when implemented in PyTorch – preserves gradient information.
If we were to naïvely solve a linear programming (LP) problem to optimize the
RobustCE loss within each training batch, the gradient information would be
lost and thus the optimization process would not work.

Algorithm 2 Solving LABEL:eq:sim_overapprox when $\hat{y}=1$

0: Bounded parameters $\theta^{\sharp}=(\theta^{L},\theta^{U})$, and
coefficients from CROWN-IBP
$\boldsymbol{\mu}^{u},\nu^{u},\boldsymbol{\alpha}^{u},\beta^{u}$

0:
$\underset{\theta\in\theta^{\sharp}}{\max}~{}~{}\boldsymbol{\mu}^{u}\theta+\nu^{u}\quad\textrm{s.t.
}\boldsymbol{\alpha}^{u}\theta+\beta^{u}\geq 0$

1: for $i=1$ to $n$ do

2: if
$\boldsymbol{\alpha}^{u}_{i}>0\vee(\boldsymbol{\alpha}^{u}_{i}=0\wedge\boldsymbol{\mu}^{u}_{i}>0)$
then

3: $\theta_{i}=\theta_{i}^{U}$

4: else

5: $\theta_{i}=\theta_{i}^{L}$

6: end if

7: end for

8: $s\leftarrow\boldsymbol{\alpha}^{u}\theta+\beta^{u}$,
$s^{\prime}\leftarrow\boldsymbol{\mu}^{u}\theta+\nu^{u}$

9: if $s<0$ then

10: return $-\infty$ // Constraint not satisfied

11: end if

12: $I\leftarrow\\{1\leq i\leq
n\mid\boldsymbol{\mu}^{u}_{i}\boldsymbol{\alpha}^{u}_{i}<0\\}$

13: Sort the indices list $I$ descendingly by
$-\frac{\boldsymbol{\mu}^{u}_{i}}{\boldsymbol{\alpha}^{u}_{i}}$

14: for $i\in I$ do

15: $\delta\leftarrow|\boldsymbol{\alpha}^{u}_{i}|(\theta^{U}-\theta^{L})$,
$\delta^{\prime}\leftarrow|\boldsymbol{\mu}^{u}_{i}|(\theta^{U}-\theta^{L})$

16: if $\delta>s$ then

17: $s^{\prime}\leftarrow s^{\prime}+\delta^{\prime}\frac{s}{\delta}$

18: break

19: end if

20: $s\leftarrow s-\delta$, $s^{\prime}\leftarrow s^{\prime}+\delta^{\prime}$

21: end for

22: return $s^{\prime}$

##  Appendix B Proofs

We provide proofs of Theorems 4.1 and 4.3.

###### Proof of Theorem 4.1.

First, we will show that the IBP loss provides a sound upper bound on
$\mathcal{L}_{\mathrm{R}}$. By the soundness of IBP, we have

$\mathcal{L}_{\mathrm{R}}^{\sharp\mathrm{IBP}}(\mathbf{x},\mathbf{x}^{\prime},\theta_{f})\geq\max_{f_{\mathrm{m}}\in\mathcal{M}_{f}}\mathcal{L}_{\mathrm{MSE}}(\bar{f_{\mathrm{m}}}(\mathbf{x}^{\prime}),1-\hat{y})$

Next, note that $\mathcal{M}_{f,\mathbf{x}}\subseteq\mathcal{M}_{f}$ ensures
the overapproximation can still capture the worst-case classifier. Formally,

$\max_{f_{\mathrm{m}}\in\mathcal{M}_{f}}\mathcal{L}_{\mathrm{MSE}}(\bar{f_{\mathrm{m}}}(\mathbf{x}^{\prime}),1-\hat{y})\geq\mathcal{L}_{\mathrm{R}}(\mathbf{x},\mathbf{x}^{\prime},\theta_{f})$

The proof of CROWN-IBP loss is tighter than the IBP loss is given in Zhang et
al. (2020). ∎

###### Proof of Theorem 4.3.

Without loss of generality, we consider the case $\hat{y}=1$. Simul-CROWN
applies CROWN-IBP to obtain lower and upper bounds in $\mathcal{M}_{f}$. For
all parameters $\theta\in\theta^{\sharp}$, we have
$\sigma(\boldsymbol{\alpha}^{l}\theta+\beta^{l})\leq\bar{f}_{\theta}(\mathbf{x})\leq\sigma(\boldsymbol{\alpha}^{u}\theta+\beta^{u})$
and
$\sigma(\boldsymbol{\mu}^{l}\theta+\nu^{l})\leq\bar{f}_{\theta}(\mathbf{x})\leq\sigma(\boldsymbol{\mu}^{u}\theta+\nu^{u})$.
The coefficients
$\boldsymbol{\alpha}^{l},\beta^{l},\boldsymbol{\alpha}^{u},\beta^{u},\boldsymbol{\mu}^{l},\nu^{l},\boldsymbol{\mu}^{u},\nu^{u}$
in Simul-CROWN are the same as the ones in CROWN-IBP.

Simul-CROWN computes the upper bound of $\mathcal{L}_{\mathrm{R}}$ as
$\mathcal{L}_{\mathrm{R}}^{\sharp\mathrm{Ours}}(\mathbf{x},\mathbf{x}^{\prime},\theta_{f})=\mathcal{L}_{\mathrm{MSE}}(\sigma(t_{\emph{Ours}}),0)$,
where $t_{\emph{Ours}}$ is obtained in the following optimization problem,

$\displaystyle
t_{\emph{Ours}}=\underset{\theta\in\theta^{\sharp}}{\max}~{}~{}\boldsymbol{\mu}^{u}\theta+\nu^{u}\quad\textrm{s.t.
}\boldsymbol{\alpha}^{u}\theta+\beta^{u}\geq 0.$ (4)

CROWN-IBP computes the upper bound of $\mathcal{L}_{\mathrm{R}}$ as
$\mathcal{L}_{\mathrm{R}}^{\sharp\mathrm{CIBP}}(\mathbf{x},\mathbf{x}^{\prime},\theta_{f})=\mathcal{L}_{\mathrm{MSE}}(\sigma(t_{\emph{CIBP}}),0)$,
where $t_{\emph{CIBP}}$ is obtained in the following optimization problem,

$\displaystyle
t_{\emph{CIBP}}=\underset{\theta\in\theta^{\sharp}}{\max}~{}~{}\boldsymbol{\mu}^{u}\theta+\nu^{u}.$
(5)

As Equation 4 has an additional constraint than Equation 5, we have
$t_{\emph{CIBP}}\geq t_{\emph{Ours}}$. Therefore, we have
$\mathcal{L}_{\mathrm{MSE}}(\sigma(t_{\emph{CIBP}}),0)\geq\mathcal{L}_{\mathrm{MSE}}(\sigma(t_{\emph{Ours}}),0)$,
which leads to

$\mathcal{L}_{\mathrm{R}}^{\sharp\mathrm{CIBP}}(\mathbf{x},\mathbf{x}^{\prime},\theta_{f})\geq\mathcal{L}_{\mathrm{R}}^{\sharp\mathrm{Ours}}(\mathbf{x},\mathbf{x}^{\prime},\theta_{f}).$

∎

##  Appendix C Additional experimental results

###  C.1 Experimental setup

##### Details on datasets

Details on dataset size and feature composition (number of features, as well
as continuous/categorical breakdown) are summarized in Table 4. We create our
own distribution shift for CTG and WHO as follows. CTG originally has three
classes (normal, suspect, and pathological). We create a distribution shift by
first training on only the “normal” and “suspect” samples. Then, we add in the
”pathological” samples (but assign them the same label as the “suspect”
samples). WHO contains data about countries’ life expectancies across several
years (2000-2015). We create a temporal data shift by using pre-2012 data as
the original dataset, and all data as the shifted dataset. When preprocessing
WHO, we also binarize the outcome variable. Instead of predicting the exact
life expectancy in years, we predict whether or not it is greater than the
pre-2012 median life expectancy.

Table 4: Dataset composition Dataset | Size (orig.) | Size (shifted) | # cont. feat. | # cat. feat.
---|---|---|---|---
HELOC | 9871 | - | 22 | 0
TC | 30000 | - | 14 | 9
OULA | 32593 | - | 22 | 8
CTG | 1950 | 2126 | 22 | 0
WHO | 2196 | 2928 | 18 | 0

##### Choosing $\delta$

We dynamically determine a distinct $\delta_{i}$ value for the $i$-th layer
with parameter $\theta_{i}$ based on its $l_{p}$ norm, denoted as
$\delta_{i}=\kappa\|\theta_{i}\|_{p}$, where $\kappa\in[0,1]$ is the ratio of
$\delta_{i}$ to the layer norm. This design is motivated by two
considerations. Firstly, parameters across different layers have different
$l_{p}$ norms due to layer size and depth differences. Consequently, a fixed
$\delta$ may be overly restrictive for some layers and lenient for others.
Secondly, using a single fixed $\delta$ could potentially inflate the norm of
each layer artificially to enhance $\delta$-Robustness because larger layer
norms make the $\delta$-robust multiplicity set relatively smaller under a
fixed $\delta$.

###  C.2 Class-level data on CTG

We include additional analysis on class-level performance for CTG. Recall that
the original dataset has 3 classes: normal (N), suspect (S), and pathological
(P). We train the original model only on N and S data, and then finetune a
shifted model with the addition of P data (assigned to the same class as S).
We perform this analysis because the dataset is very imbalanced: around 78%
belongs to class N. That is, we want to ensure that the finetuning process
actually adapts the model and CE generator to the shifted data (rather than
just yielding an overall high performance at the expense of the minority
groups S and P).

Before finetuning, all CEs are valid regardless of class. After finetuning, we
find that overall, 98.7% of CEs generated for CTG remain valid. For samples
with class=S (i.e., the minority samples that are represented in the training
data), 95.9% of CEs are valid after finetuning. And for samples with class=P
(i.e., only represented in the finetuning data), 96.8% of CEs remain valid.
So, while robustness to finetuning is slightly lower for the minority classes,
robustness is very high for all groups.

###  C.3 Ablation Studies

##### Other distance metrics

In addition to proximity (see Section 5.1 for a definition), we consider
_sparsity_ and _distance to the data manifold (DDM)_ as CE quality metrics.
Sparsity is defined as the fraction of features that change, i.e., given
$\mathbf{x}$ and its $\mathbf{x}^{\prime}$,
$\textsc{SPARS}(\mathbf{x},\mathbf{x}^{\prime})=\frac{1}{d}\sum_{i}^{d}\mathbb{1}[\mathbf{x}_{i}\neq\mathbf{x}^{\prime}_{i}]$.
We approximate DDM by taking the distance to the nearest point in the training
data (by using $k$-nearest neighbors with $k=1$). For all three metrics, a
lower score indicates smaller distance, which indicates that the CE is of
higher quality because it will be easier for the user to obtain. Table 5
contains data on sparsity and DDM; see Table 3 in the main paper for proximity
data.

Table 5: Average sparsity and distance to the data manifold (DDM) for valid CE. (*) indicates that the CE technique does not correctly account for categorical features. |  | CTG | WHO | HELOC | TC | OULA
---|---|---|---|---|---|---
Proximity | VeriTraCER | 0.780 | 0.758 | 0.817 | 0.758 | 0.424
CN | 0.765 | 0.718 | 0.791 | 0.751 | 0.403
ROAR | 0.900 | 0.326 | 0.286 | 0.166* | 0.399*
SNS | 0.585 | 0.338 | 0.541 | 0.210* | 0.292*
DDM | VeriTraCER | 0.130 | 0.103 | 0.050 | 0.123 | 0.051
CN | 0.109 | 0.064 | 0.049 | 0.108 | 0.039
ROAR | 0.072 | 0.042 | 0.048 | 0.020* | 0.028*
SNS | 0.052 | 0.040 | 0.054 | 0.015* | 0.024*

##### Relation of our approach’s effectiveness and the tightness of the bound.

We compare the three different techniques from Section 4.2. Using an looser
bound such as IBP or CROWN-IBP in place of Simul-CROWN is advantageous for
computation speed, however, the overly loose bounds may hurt the
effectiveness. For instance, we expect to see lower model accuracy and higher
CE distance metrics with IBP and CROWN-IBP than with Simul-CROWN. However,
Table 6 shows that three techniques perform largely similarly. Note that the
$\delta$-Robustness rate is computed post-training using Simul-CROWN as it
achieves tighter over-approximation than the other two approaches.

Table 6: Fraction of samples that exhibit $\delta$-robustness and empirical random initialization (RI) robustness, along with test accuracy and CE quality metrics for IBP, CROWN-IBP, and VeriTraCER. Prox. is proxmity, Spars. is sparsity, and Man. is distance to the data manifold. |  |  $\Delta$-Rob. | RI | Acc. | Prox. | Spars. | Man.
---|---|---|---|---|---|---|---
HELOC | IBP | 0.746 | 0.920 (0.079) | 0.742 | 0.10 (0.03) | 0.82 (0.08) | 0.05 (0.01)
C-IBP | 0.721 | 0.918 (0.064) | 0.742 | 0.10 (0.03) | 0.82 (0.08) | 0.05 (0.01)
VeriTraCER | 0.788 | 0.932 (0.028) | 0.738 | 0.10 (0.03) | 0.82 (0.08) | 0.05 (0.01)
TC | IBP | 0.902 | 0.862 (0.187) | 0.812 | 0.22 (0.06) | 0.74 (0.09) | 0.10 (0.03)
C-IBP | 0.872 | 0.900 (0.120) | 0.813 | 0.22 (0.06) | 0.75 (0.09) | 0.10 (0.03)
VeriTraCER | 0.813 | 0.884 (0.141) | 0.806 | 0.24 (0.05) | 0.76 (0.09) | 0.12 (0.03)
OULA | IBP | 0.880 | 0.922 (0.079) | 0.929 | 0.19 (0.04) | 0.45 (0.12) | 0.06 (0.02)
C-IBP | 0.968 | 0.941 (0.055) | 0.930 | 0.18 (0.04) | 0.43 (0.12) | 0.05 (0.02)
VeriTraCER | 0.970 | 0.941 (0.054) | 0.929 | 0.18 (0.04) | 0.42 (0.12) | 0.05 (0.02)

###  C.4 Time complexity

For VeriTraCER, it takes 8.24 (±0.14) s. to train one epoch on HELOC, 126.61
(±1.92) s. to train one epoch on OULA, and 27.66 (±0.36) s. to train one epoch
on TC. By contrast, training one epoch with standard training on the same
hardware takes 0.55 (±0.01) s. for HELOC, 5.39 (±0.06) s. for OULA, and 1.80
(±0.01) s. for TC. We train each dataset for 100 epochs, so with VeriTraCER,
it takes on the order of 14 minutes for HELOC, 3 hours for OULA, 45 minutes
for TC. With standard training, those times are reduced to around 1 minute, 9
minutes, and 3 minutes, respectively.

Table 7 summarizes the time required to generate one CE, given a trained
model. (CounterNet is not included, as it takes the same amount of time as
VeriTraCER.) We see that VeriTraCER can generate CEs near-instantaneously,
while ROAR and SNS must solve complex optimization problems for each instance.

Table 7: Average time, in seconds, to generate one counterfactual given a trained model. The standard deviation for VeriTraCER is negligible ($<10^{-4}$) for all datasets. | HELOC | OULA | TC
---|---|---|---
VeriTraCER |  $<$0.001 |  $<$0.001 |  $<$0.001
ROAR | 3.25 ±0.72 | 2.75 ±0.52 | 3.33 ±0.40
SNS | 5.25 ±0.21 | 5.92 ±0.21 | 5.89 ±0.20


# Entangled Photons Enabled Time- and Frequency-Resolved Coherent Raman
Spectroscopy in Condensed Phase Molecules

Zhedong Zhang<EMAIL_ADDRESS>City University of Hong Kong, Kowloon
Tong, Hong Kong SAR  Tao Peng  Texas A$\&$M University, College Station, Texas
77843, United States  Xiaoyu Nie  Xi’an Jiaotong University, Xi’an, Shaanxi
710049, China  Girish S. Agarwal  Texas A$\&$M University, College Station,
Texas 77843, United States  Marlan O. Scully  Texas A$\&$M University, College
Station, Texas 77843, United States  Baylor University, Waco, Texas 76704,
United States  Princeton University, New Jersey 08544, United States

###### Abstract

We develop an ultrafast frequency-resolved Raman spectroscopy with entangled
photons for polyatomic molecules in condensed phases, to probe the electronic
and vibrational coherences. Using quantum correlation between the photons, the
signal shows the capability of both temporal and spectral resolutions that are
not accessible by either classical pulses or the fields without entanglement.
We develop a microscopic theory for this Raman spectroscopy, revealing the
electronic coherence dynamics which often shows a rapid decay within $\sim
50$fs. The heterodyne-detected Raman signal is further developed to capture
the phases of electronic coherence and emission in real-time domain.

Introduction.–Over last few decades quantum entanglement has emerged to have
considerable applications in testing fundamental laws of physics [1, 2, 3, 4,
5, 6, 7], as well as revolutionizing photonic applications, e.g., quantum
imaging [8], metrology [9, 10], sensing [11, 12], and quantum computing [13,
14].

The field of spectroscopy especially the nonlinear spectroscopy can benefit
enormously from the quantum states of light. The role of quantum entanglement
in the context of two-photon spectroscopy was realized rather early. Fei et
al. demonstrated entanglement-induced transparency in two-photon processes
[15]. The occurrence of simultaneous excitation two independent atoms due to
entangled photons have been demonstrated in recent experiments [16, 17]. In
general the entangled photons have demonstrated incredible power in developing
new spectroscopic techniques for complex materials [18, 19, 20, 21]. Thanks to
the technical advance in pulse shaping, the two-photon absorption with photon
entanglement has been demonstrated in various experiments, incredibly
suppressing the intermediate state relaxation that yields the optimal
efficiency of populating higher excited states [22, 23, 24, 25, 26]. The
entangled states of photons provide extra control knobs for selective
molecular relaxation and radiative processes. The related signals may be
extended into sub-picosecond time-resolved regimes, resulting in an ultrafast
optical technique to probe the dynamics of various excitations. Recent work on
stimulated Raman and pump-probe spectroscopies showed the enhanced resolution
via photon entanglement through the model calculations [27, 28, 29].
Nevertheless, quantum states of light may be powerful for exploiting the
ultrafast electronic processes in semiconductor materials. This has been
indicated from the semiconductor quantum-light sources taking the advantage of
generating high photon flux and the scalability [30, 31]. More delicate
information about the structure and dynamics in metal nanostructures and low-
dimensional materials can be read out from the photon correlation [32, 33, 34,
37, 36, 35]. Subsequent theoretical proposals aiming at ultrafast multi-photon
coincidence counting spectrum revealed the highly-resolved excitonic
relaxation pathways and charge fluctuation through manipulating the temporal
and spectral gates, not accessible by classical pulses [37, 38, 39]. It turns
out the unusual band spectrum of entangled photons responsible for the
superresolution nature beyond Heisenberg conjugation renders the quantum-light
spectroscopy as a new promising route towards novel applications in ultrafast
and ultra-sensitive tomography [18, 40, 41, 42, 43, 44].

In this Letter, we demonstrate new capabilities brought out by the use of
entangled photons in femtosecond time-resolved coherent Raman spectroscopy, to
monitor the electronic and vibrational coherences. This leads to the QFRS
(Quantum femtosecond Raman spectroscopy). Specifically, we show the time-
frequency resolutions, obtained by the photon entanglement. Explicit results
are presented for a quantum extension of FAST CARS (Femtosecond adaptive
spectroscopic technique for coherent anti-Stokes Raman scattering) [45, 46,
47]. Extensive works have been done in the population relaxation with
entangled-photon spectroscopy [28, 48]. The ultrafast dynamics of electronic
coherence fundamentally important for reaction kinetics and internal
conversions, however, is a challenging task due to its rapid decay [49, 50,
51, 52, 53, 54]. We show how the entangled-photon spectroscopy could be useful
here.

Figure 1: (a) Schematic of entangled twin photons as ultrafast probe for
molecules, where the parametric down conversion through a beta baruim borate
(BBO) crystal and multi-photon detection are presented. (b) Quantum heterodyne
scheme of detection where a local oscillator interferes with emitted photons
in s arm; Detector 1 records the interference term whereas detector 2 records
the idler photon number, yielding the joint detection. (c) Level scheme of
microscopic model in Quantum FAST CARS. (d) Level scheme of microscopic model
in QFRS for electronically excited states. (e) Quantum FAST CARS signal from
Eq.(3), taking the 4 Raman-active modes $\text{A}_{1}$, E and $\text{T}_{2}$
in methane ($\mathrm{CH}_{4}$).

Quantum FAST CARS.–To start off, we consider the FAST CARS with entangled
photons, where the molecules are driven by a pair of classical pulses. After a
time delay, the beam in s arm of the entangled twin photons serves as a Raman
probe off-resonantly interacting with molecules, whereas the idler beam
propagates freely and provides a reference, as depicted in Fig.1(a). Referring
to the level scheme in Fig.1(c), the field-molecule interaction of Raman
process is

$\begin{split}V(t)=\sum_{j=1}^{N}\sum_{b}\alpha_{bg}^{(j)}|b\rangle\langle
g|_{j}E_{s}(t)E_{s}^{\dagger}(t)e^{i(\omega_{b}-\omega_{g})t}+\text{h.c.}\end{split}$
(1)

where $\alpha_{bg}^{(j)}$ is the Raman polarizability and the field $E_{s}(t)$
contains multiple frequency modes [55]. $N$ denotes the number of molecules.
$|\Psi\rangle=\frac{1}{\sqrt{{\cal
N}}}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\text{d}\omega_{s}\text{d}\omega_{i}\
\Phi(\omega_{s},\omega_{i})|1_{\omega_{s}},1_{\omega_{i}}\rangle$ and

$\begin{split}\Phi(\omega_{s},\omega_{i})=E_{0}(\omega_{s}+\omega_{i})\text{sinc}\left[\frac{\Delta
k(\omega_{s},\omega_{i})L}{2}\right]e^{i\Delta
k(\omega_{s},\omega_{i})L/2}\end{split}$ (2)

is the wave function for the entangled twin photons with $\Delta
k(\omega_{s},\omega_{i})L=(\omega_{s}-\frac{\omega_{0}}{2})T_{s}+(\omega_{i}-\frac{\omega_{0}}{2})T_{i}$
where $T_{s}(T_{i})$ is the time delay between the photons in $s$ (reference)
arm and the SPDC pump field, due to the group velocity dispersion in the
nonlinear crystal. ${\cal N}$ is a normalization factor.

A joint detection of the spectral-resolved transmissions in two arms gives the
intensity-correlated signal $S(\omega,\omega_{i};T)=\langle
E_{s}^{\dagger}(\omega)E_{i}^{\dagger}(\omega_{i})E_{i}(\omega_{i})E_{s}(\omega)\rangle_{\rho}$
where $E_{s}(\omega)$ and $E_{i}(\omega_{i})$ denote the Fourier component of
the electric fields in s and reference arms, respectively. $T$ is the time
delay of $s$-arm photons relative to the pump pulse. Using the perturbation
expansion against Eq.(1), we find the 8-point field correlation function
playing a key role in understanding the unusual spectroscopic properties of
Raman signals with photon entanglement. A lengthy algebra leads to the Quantum
FAST CARS signal

$\begin{split}S_{\text{QV}}(\omega,\omega_{i};T)=&\frac{N(N-1)}{8\pi^{2}{\cal
N}}|{\cal E}_{s,\omega}|^{6}|{\cal
E}_{i,\omega_{i}}|^{2}\bigg{|}\sum_{b}\alpha_{bg}^{*}\rho_{bg}(T)\\\\[4.26773pt]
&\qquad\qquad\times\Phi(\omega-\omega_{bg}-i\gamma_{bg},\omega_{i})\bigg{|}^{2}\end{split}$
(3)

with $\rho_{bg}(t)=\rho_{bg}e^{-(i\omega_{bg}+\gamma_{bg})t}$ where
$\gamma_{bg}^{-1}$ quantifies the dephasing of the vibrational coherence.
${\cal E}_{s(i),\omega}=\sqrt{2\pi\omega_{s(i)}/\text{V}}$ and V is bulk
volume. Besides the enhanced signal intensity via coherence, Eq.(3) indicates
the spectral resolution governed by photon entanglement that contains new
control knobs such as photon arrival time and bandwidth.

Figure 2: (First column) (a) Intensity-correlated QFRS for time-evolving
electronic coherence versus the delay $T$ between entangled photons and
resonant pump pulse; (b) Same as (a) but for heterodyne detection. (c,d) Same
as (a,b) but using classical pulse as probe having bandwidth
$\sigma_{0}=0.82$eV. (e,f) Same as (a,b) but using uncorrelated separable
state of photons with the bandwidth $\sigma_{0}=0.82$eV
($\sigma_{0}^{-1}=5$fs). Parameters are taken from 4-oriented
amino-4’nitrostilbene [56], i.e., $\omega_{e}=7.1$eV, $\omega_{e_{1}}=5.3$eV,
$\omega_{e_{2}}=5.7$eV, $v_{h}=0.26$eV, $F_{1}=2.2,\ F_{2}=1.3$,
$D_{1}^{-1/2}=30$fs, $D_{2}^{-1/2}=20$fs and $T_{s}=T_{i}=30$fs.

For a neat understanding of the signal, we focus on the optimal photon
entanglement, i.e., $T_{s}=T_{i}$. The two-photon amplitude in Eq.(2)
indicates a Raman resonance $\omega-\omega_{\text{pr}}=\omega_{bg}$ defining
the probe frequency $\omega_{\text{pr}}=\omega_{0}-\omega_{i}$. This leads to
the spectral resolution $T_{s}^{-1}$, which reveals a considerable enhancement
not accessible by classically shaped pulses. To make this clear, we calculate
the Raman signals with two incoherent states of light: (i) Fock state of
photons that is fully separable and (ii) the psuedo-thermal light
$\rho_{\text{pht}}=\frac{1}{{\cal
N}}\iint\text{d}v\text{d}v_{i}|\Phi(v,v_{i})|^{2}|1_{s,v};1_{i,v_{i}}\rangle\langle
1_{s,v};1_{i,v_{i}}|$ having classical correlation which has been used in
information processing. For (i), the spectral resolution is governed by the
single-photon amplitude $\Phi_{s}(\omega_{s}-\omega_{0})\propto
e^{-(\omega_{s}-\omega_{0})^{2}/2\sigma_{0}^{2}}$ having the identical time
duration $\sigma_{0}^{-1}$ as the pump pulse, as seen from Eq.(16) in
Ref.[55]. In CARS experiments, the conjugated time and frequency resolutions
results in the challenge to distinguish the vibrational modes once
$|\omega_{b_{1},b_{2}}|\ll\sigma_{0}$. Using photon entanglement, the spectral
resolution can be controlled independently and can be therefore greatly
enhanced such that $T_{s}^{-1}<|\omega_{b_{1},b_{2}}|\ll\sigma_{0}$, evident
by Fig.1(e) that keeps high temporal resolution. For (ii) as entailed by
Eq.(15) in Ref.[55], the signal shows a high spectral resolution $T_{s}^{-1}$,
whereas the temporal resolution is no longer there. This is expected from the
random phase between photon pairs that results in arbitrary arrival time of
photons.

Intensity-correlated QFRS.–To monitor the ultrafast dynamics of electronic
coherence in excited states of molecules, the probe pulse has to act after a
delay $T$ following a resonant photoexcitation or a propagation of the excited
states, depicted in Fig.1(d). The Raman signal is produced when the s-arm
photon scatters off the excited-state coherence resulting in the electronic
polarizability $\alpha(t)$ describing the Raman transitions between excited
states. Proceeding via the algebra similar as before, we find the QFRS signal
[55]

$\begin{split}S_{\text{QE}}(\omega,\omega_{i};T)=\frac{N(N-1)}{32\pi^{4}{\cal
N}}|{\cal E}_{s,\omega}|^{6}|{\cal
E}_{i,\omega_{i}}|^{2}\bigg{|}\sum_{e^{\prime}}\alpha_{e,e^{\prime}}^{*}f_{e,e^{\prime}}(T)\bigg{|}^{2}\end{split}$
(4)

where the Raman line-shape function is defined as

$\begin{split}f_{e,e^{\prime}}(T)&=\int_{-\infty}^{\infty}\text{d}\tau\int_{-\infty}^{\infty}\text{d}\omega^{\prime}\rho_{e,e^{\prime}}(\tau)e^{i(\omega-\omega^{\prime})(\tau-T)}\Phi(\omega^{\prime},\omega_{i}).\end{split}$
(5)

Eq. (5) indicates the role of the photon entanglement whose unusual band
properties may provide versatile tool for controlling the ultrafast electron
dynamics in molecules.

To elaborate this, the electronically excited-state dynamics coupled to
nuclear motions has to be considered. This may facilitate a quantitative
understanding of the inhomogeneous dephasing in realistic materials. We adopt
the Fröhlich-Holstein model for molecules

$\begin{split}H=H_{\text{vib}}|g\rangle\langle
g|+\sum_{i}\left(\omega_{e_{i},g}+\lambda_{s}^{(i),2}v_{s}+H_{\text{vib}}^{i}\right)|e_{i}\rangle\langle
e_{i}|\end{split}$ (6)

and $H_{\text{vib}}=\sum_{s}v_{s}b_{s}^{\dagger}b_{s}$ and
$H_{\text{vib}}^{i}=\sum_{s}\left[v_{s}b_{s}^{\dagger}b_{s}-\lambda_{s}^{(i)}v_{s}(b_{s}+b_{s}^{\dagger})\right]$,
where $b_{s}$ is the annihilation operator of the $s$-th vibration having the
frequency $v_{s}$. $\lambda_{s}^{(i)}$ quantifies the vibronic coupling for
the $i$-th electronic state. Propagating the nuclear wave packets and using
Eq.(4) and (5), the intensity-correlated QFRS signal can be found

$\begin{split}&S_{\text{QE}}(\omega,\omega_{i};T)=\frac{N(N-1)}{8\pi^{2}{\cal
N}}|{\cal E}_{s,\omega}|^{6}|{\cal E}_{i,\omega_{i}}|^{2}\\\\[4.26773pt]
&\qquad\qquad\times\bigg{|}\sum_{e_{1}}\sum_{n=0}^{\infty}\alpha_{e,e_{1}}^{*}g_{n,e_{1}}(\omega,T)\rho_{e,e_{1}}^{(n)}(T)\bigg{|}^{2}\end{split}$
(7)

where $g_{n,e_{j}}(\omega,T)=\int_{-\infty}^{\infty}\text{d}\tau\
e^{i(\omega-\tilde{\omega}_{e,e_{j}}-nv_{h})\tau}e^{-D_{j}(\tau^{2}+2T\tau)}\tilde{\Phi}(\tau,\omega_{i})$
and
$\tilde{\Phi}(\tau,\omega_{i})\equiv\frac{1}{2\pi}\int_{-\infty}^{\infty}\Phi(\omega,\omega_{i})e^{-i\omega\tau}\text{d}\omega$
[55].
$\rho_{e,e_{j}}^{(n)}(t)/\rho_{e,e_{j}}(0)=S_{n,j}e^{-i(\tilde{\omega}_{e,e_{j}}+nv_{h})t-D_{j}t^{2}}$
is the vibronic coherence with $n$ harmonics of the high-frequency vibrations
and $S_{n,j}=e^{-F_{j}}F_{j}^{n}/n!$ is the Franck-Condon factor.
$D_{j}=m(\lambda_{l}^{(j)}-\lambda_{l}^{(e)})^{2}$ and
$F_{j}=(\lambda_{h}^{(j)}-\lambda_{h}^{(e)})^{2}$ quantify the coupling of
excitons to low- and high-frequency vibrations, respectively [57]. The
superposition of vibronic coherences yields the electronic coherence, i.e.,
$\rho_{e,e_{j}}(t)=\sum_{n=0}^{\infty}\rho_{e,e_{j}}^{(n)}(t)$.

Notably, for the timescale $T<1/\sqrt{D_{j}}$, we can find the line-shape
function for arbitrary two-photon amplitude, namely,
$g_{n,e_{j}}(\omega,T)\simeq\Phi(\omega-\tilde{\omega}_{e,e_{j}}-nv_{h}+2iD_{j}T,\omega_{i})$.
Eq.(7) indicates the multiple Raman resonance
$\omega-\omega_{\text{pr}}=\tilde{\omega}_{e,e_{1}}+nv_{h}$ in which the
intensity is governed by the Franck-Condon factor quantifying the vibronic
interaction. The spectral resolution is $T_{s}^{-1}$ that may be greatly
enhanced, as dictated by the two-photon amplitude. The Raman signals in Eq.(7)
are highly time- and frequency-resolved, revealing the ultrafast coherent
dynamics of the electronically excited states.

Figure 3: Heterodyne-detected QFRS and dynamics of vibronic coherence. (Left
column) (a) Heterodyne-detected QFRS with the phase associated with
$\Phi(\bar{\omega},\omega_{i})$ being $-\pi/2$ (a certain phase of local
oscillator). (b, up) Red line is for the 1D slice of (a) at
$\omega-\omega_{\text{pr}}=\omega_{e,e_{1}}+2v_{h}=2.32$eV; Black line is for
the imaginary part of $n=2$ harmonics of the vibronic coherence
$\rho_{e,e_{1}}$, i.e., $\text{Im}[\rho_{e,e_{1}}^{(2)}]$. (b, down) Same as
(b, up) but for $\omega-\omega_{\text{pr}}=\omega_{e,e_{2}}+v_{h}=1.66$eV and
imaginary part of $n=1$ harmonics of the vibronic coherence $\rho_{e,e_{2}}$.
(c) Imaginary part of time-evolving vibronic coherence, i.e.,
$\text{Im}[\rho_{e,e_{1}}]$ (black) and $\text{Im}[\rho_{e,e_{2}}]$ (red).
(Right column) (d) Same as (a) but for the phase associated with
$\Phi(\bar{\omega},\omega_{i})$ being $-\pi$. (e) Same as (b) but for 1D
slices of (d) (red lines) and $\text{Re}[\rho_{e,e_{1}}^{(2)}],\
\text{Re}[\rho_{e,e_{2}}^{(1)}]$ (black lines). (f) Same as (c) but for
$\text{Re}[\rho_{e,e_{1}}]$ and $\text{Re}[\rho_{e,e_{2}}]$. Parameters are
the same as Fig. 2.

Fig. 2(a) illustrates the intensity-correlated QFRS, where the electronic
coherence $\rho_{e,e_{i}};i=1,2$ is monitored, yielding the Raman resonance
corresponding to the two electronically excited states at
$\omega_{e_{1}}=5.3$eV and $\omega_{e_{2}}=5.7$eV. Using the short pulse, we
can readily resolves the vibronic couplings. In particular, the signal
contains the Raman resonance
$\omega-\omega_{\text{pr}}=\omega_{e,e_{1}}+nv_{h}\ (n=1,2,3,4)$ and
$\omega_{e,e_{2}}+mv_{h}\ (m=0,1,2)$, where the peaks at
$\omega_{e,e_{1}}+nv_{h}$ decay slower and dominate at longer time delay $T$.
This is because of the weaker influence from the low-frequency vibrations
leading to slower dephasing of electronic coherence. It is worth noting that
only a few harmonics of the vibrational modes can be seen from the Raman
signal for electronic coherence. This results from the Franck-Condon effect
yielding a Poissonian distribution of the vibronic energy seen from
$S_{\text{QE}}(\omega,\omega_{i};T=0)$. Fig. 2(c) illustrates the intensity-
correlated QFRS with a classical probe pulse. The broadband pulse has a high
temporal resolution, but the individual electronically excited states, as well
as their coupling to vibrations are poorly resolved. This has been a long
existing bottleneck in ultrafast Raman spectroscopy.

We further examine the Fock state of photons that is fully separable with
neither correlation nor entanglement, in Fig. 2(e). It shows that the fully
separable Fock photon state yields as poor resolutions as the Raman signals
with classical probe. This is because the field correlation function is
factorized into a product of field intensities.

It is worth noting that the intensity-correlated QFRS may smear out the phase
information of the excited states, as seen from the spectral-resolved
condition $T_{s}^{-1}<|\omega_{e_{1},e_{2}}|$ so that
$|\rho_{e,e^{\prime}}|^{2}$ dominates the signal. A new spectroscopic signal
is therefore essential to be exploited, as will be entailed next.

Heterodyne-detected QFRS.–We let the s-arm photons serve as a local oscillator
interfering with the emitted photons, as depicted in Fig.1(b). The field
directed into detector 1 carries the superposition state
$|e\rangle+e^{i\phi}|e^{\prime}\rangle$ assuming 50/50 beam splitters. The
joint detection of photons with anti-Stokes shift yields the measurement
$\propto\text{Re}(e^{-i\phi}\rho_{e,e^{\prime}})$. In such a spirit, we find
the heterodyne-detected signal explicitly [55]

$\begin{split}S_{\text{QEHD}}(\omega,\omega_{i};T)=\frac{N|{\cal
E}_{s,\omega}|^{2}|{\cal E}_{i,\omega_{i}}|^{2}}{\pi{\cal
N}}\text{Im}\bigg{[}\sum_{e^{\prime}}\alpha_{e,e^{\prime}}^{*}\Phi^{*}(\bar{\omega},\omega_{i})f_{e,e^{\prime}}(T)\bigg{]}\end{split}$
(8)

with the same definition of Raman line-shape function and electronic
polarizability as in Eq.(4) and (5).

For the molecules in condensed phases, we incorporate into the signal the
nuclear dynamics interacting with electronically excited states. The
heterodyne-detected QFRS reads

$\begin{split}S_{\text{QEHD}}(\omega,\omega_{i};T)=&\frac{2N}{{\cal N}}|{\cal
E}_{s,\omega}|^{2}|{\cal
E}_{i,\omega_{i}}|^{2}\sum_{e_{1}}\sum_{n=0}^{\infty}\text{Im}\Big{[}\alpha_{e,e_{1}}^{*}\rho_{e,e_{1}}^{(n)}(T)\\\\[2.84544pt]
&\qquad\qquad\times\Phi^{*}(\bar{\omega},\omega_{i})g_{n,e_{1}}(\omega,T)\Big{]}\end{split}$
(9)

where $g_{n,e_{j}}(\omega,T)$ is of the identical definition as in Eq.(7).
From Eq.(9), $S_{\text{HD}}(\omega,\omega_{i};T)$ visualizes the full
transient coherence dynamics including phase and intensity, via the amplitude
$\Phi(\bar{\omega},\omega_{i})$ whose phase varies with the optical path.

To track the fast coherence dynamics in further, we plot the heterodyne-
detected QFRS in Fig. 2(b), which unambiguously reveals the oscillations. The
heterodyne signal predicts the same Raman resonance as the signal in Fig.
2(a). The advantage of photon entanglement is shown in Fig. 2(d) which
produces a much worse spectral resolution with a classical broadband probe,
despite of the temporal resolution. This is extensively supported from the
comparison to the heterodyne signal using Fock state of photons, as depicted
in Fig. 2(e) where neither correlation nor entanglement has been involved.

Fig. 3 shows the heterodyne-detected QFRS with different phases of the local
oscillator field that leads to a global phase in
$\Phi(\bar{\omega},\omega_{i})$. It can be seen from Fig. 3(a) and 3(d) the
fast oscillations that selectively track the real-time dynamics of the
coherence phase (imaginary and real parts), given the high spectral resolution
revealing the vibronic states. As time evolves, Fig. 3(b) and 3(e) illustrate
the Fourier components of the vibronic coherence (imaginary and real parts)
resolved from the oscillations with various frequencies, e.g.,
$\omega_{e,e_{1}}+2v_{h}$ and $\omega_{e,e_{2}}+v_{h}$ as a result of the
vibronic coupling. This leads to the global beating feature of the vibronic
coherence, depicted in Fig. 3(c) and 3(f) revealing the phase information.

Experimental feasibility.–Since the experiments using electronically ground
states are easier and have been accomplished, we will focus on the feasibility
of Quantum FAST CARS. The vibrational coherence is scattered off by the s
photons entangled with the idler photons having a delay from the Stokes pulse.
Dramatically different from the FAST CARS, the probe pulse using entangled
twin photons can be broadband, rather than the narrow-band one. Typically, the
time duration of fast probe is taken to be $\sigma_{0}^{-1}=35$fs and the
delay controlled by the crystal length is $T_{s}=1$ps yielding the Raman
signal with frequency resolution of 33cm-1 sufficient for a variety of organic
molecules. This allows to monitor the fast-evolving vibrational coherence,
provided the same spectral resolution as the FAST CARS.

We employ the methane ($\mathrm{CH}_{4}$) which is a ubiquitous molecule with
well characterized spectral lines. The Raman-active transition
$\mathrm{A}_{1}$ at 2914$\mathrm{cm}^{-1}$ of methane, for instance, can be
probed with focusing the femtosecond pump and Stokes pulses followed by s
photons. The s photon pulse intensity is of $10^{-9}$ magnitude lower than the
classical pulse. In a typical CARS scheme the electronic transition dipole is
$\sim 1$ debye with pump detunings on the order of $10^{5}$
$\mathrm{cm}^{-1}$, the quantum light signal is only $\sim
10^{-1}/\mathrm{s}$. Nevertheless, approaching resonance by reducing detuning
with 2 orders, may give an extra 4 orders enhancement of the signal without
laser damage. One can then expect the quantum FAST CARS signal at a level of
$10^{3}/\mathrm{s}$, which is sufficient for the detection. We note here that
the quantum FAST CARS scheme possesses an inherent advantage due to the multi-
photon coincidence counting. This plays an essential role in increasing the
signal-to-noise ratio for quantum FAST CARS.

Summary.–We developed a new time- and frequency-resolved Raman spectroscopy
using entangled photons and interferometry, as a quantum extension of FAST
CARS. Substantially, QFRS was further proposed to probe the electronic
coherence. In contrast to existing paradigms [18, 24, 48], QFRS is sensitive
to the electronic coherence only and is population-free, making it uniquely
suitable for monitoring the electronically excited state dynamics during a
short timescale. We developed a microscopic theory for the signals involving
the intensity- and heterodyne-detected modes. These demonstrated that QFRS can
measure the time-varying energy gap between the vibronic states with a high
spectral resolution. The rapid oscillation and decay of electronic coherence
can be fully visible in the time-resolved spectrum. Such temporal- and
spectral-resolved nature is not attainable in conventional Raman spectra. The
QFRS, as a new coherent Raman technique, offers new perspective for the
ultrafast dynamics in photophysical systems including low-dimensional
semiconductor materials, exciton polaritons and nano-plasmonics.

Z.D.Z. and T.P. contributed equally to this work.

We gratefully thank Renbao Liu from CUHK for the instructive discussions and
suggestions. Z.D.Z. gratefully acknowledges the support of ARPC-CityU new
research initiative/infrastructure support from central (No. 9610505). T.P.,
G.S.A. and M.O.S. gratefully acknowledge the support of Air Force Office of
Research (Grant No. FA9550-20-1-0366), Office of Naval Research (Grant No.
N00014-20-1-2184), Robert A. Welch Foundation (Grant No. A-1261; A-1943), and
National Science Foundation (Grant No. PHY-2013771).

## References

  * [1] C.-K. Hong, Z.-Y. Ou and L. Mandel, Phys. Rev. Lett. 59, 2044-2046 (1987)
  * [2] P. G. Kwiat, K. Mattle, H. Weinfurter, A. Zeilinger, A. V. Sergienko and Y. Shih, Phys. Rev. Lett. 75, 4337-4342 (1995)
  * [3] Y. Shih and C. O. Alley, Phys. Rev. Lett. 61, 2921-2924 (1988)
  * [4] J. D. Franson, Phys. Rev. Lett. 62, 2205-2208 (1989)
  * [5] J. C. Howell, R. S. Bennink, S. J. Bentley and R. W. Boyd, Phys. Rev. Lett. 92, 210403-210406 (2004)
  * [6] G. Weihs, T. Jennewein, C. Simon, H. Weinfurter and A. Zeilinger, Phys. Rev. Lett. 81, 5039-5043 (1998)
  * [7] Y.-H. Kim, R. Yu, S. P. Kulik, Y. Shih and M. O. Scully, Phys. Rev. Lett. 84, 1-5 (2000)
  * [8] T. B. Pittman, Y. Shih, D. Strekalov and A. V. Sergienko, Phys. Rev. A 52, 3429-3432(R) (1995)
  * [9] M. D’Angelo, M. V. Chekhova and Y. Shih, Phys. Rev. Lett. 87, 013602-013605 (2001)
  * [10] M. W. Mitchell, J. S. Lundeen and A. M. Steinberg, Nature 429, 161-164 (2004)
  * [11] C. L. Degen, F. Reinhard and P. Cappellaro, Rev. Mod. Phys. 89, 035002-035040 (2017)
  * [12] S. Pirandola, B. R. Bardhan, T. Gehring, C. Weedbrook and S. Lloyd, Nat. Photonics 12, 724-733 (2018)
  * [13] P. Walther, et al., Nature 434, 169-176 (2005)
  * [14] J. L. O’Brien, Science 318, 1567-1570 (2007)
  * [15] H.-B. Fei, B. M. Jost, S. Popescu, B. E. A. Saleh and M. C. Teich, Phys. Rev. Lett. 78, 1679-1682 (1997)
  * [16] A. Muthukrishnan, G. S. Agarwal and M. O. Scully, Phys. Rev. Lett. 93, 093002-093005 (2004)
  * [17] W. Ren, et al., Phys. Rev. Lett. 125, 133601-133606 (2020)
  * [18] K. E. Dorfman, F. Schlawin and S. Mukamel, Rev. Mod. Phys. 88, 045008-045074 (2016)
  * [19] S. Mukamel, et al., J. Phys. B.: At. Mol. Opt. Phys. 53, 072002-072046 (2020)
  * [20] M. G. Raymer, A. H. Marcus, J. R. Widom and D. L. P. Vitullo, J. Phys. Chem. B 117, 15559-15575 (2013)
  * [21] H. Oka, J. Chem. Phys. 152, 044106-044111 (2020)
  * [22] F. Schlawin, K. E. Dorfman, B. P. Fingerhut and S. Mukamel, Nat. Commun. 4, 1782-1788 (2013)
  * [23] D.-lk Lee and T. Goodson III, J. Phys. Chem. B 110, 25582-25585 (2006)
  * [24] R. de J. Leon-Montiel, J. Svozilik, J. P. Torres and A. B. U’Ren, Phys. Rev. Lett. 123, 023601-023605 (2019)
  * [25] O. Varnavski, B. Pinsky and T. Goodson III, J. Phys. Chem. Lett. 8, 388-393 (2017)
  * [26] J. P. Villabona-Monsalve, O. Varnavski, B. A. Palfey and T. Goodson III, J. Am. Chem. Soc. 140, 14562-14566 (2018)
  * [27] F. Schlawin, K. E. Dorfman and S. Mukamel, J. Chem. Phys. 154, 104116-104121 (2021)
  * [28] F. Schlawin, K. E. Dorfman and S. Mukamel, Phys. Rev. A 93, 023807-023814 (2016)
  * [29] A. A. Svidzinsky, G. S. Agarwal, A. Classen, A. V. Sokolov, A. M. Zheltikov, M. S. Zubairy and M. O. Scully, to be published
  * [30] Y. Chen, M. Zopf, R. Keil, F. Ding and O. G. Schmidt, Nat. Commun. 9, 2994-3000 (2018)
  * [31] A. J. Shields, Nat. Photonics 1, 215-223 (2007)
  * [32] K. Wu, J. Chen, J. R. McBride and T. Lian, Science 349, 632-636 (2015)
  * [33] G. Tagliabue, et al., Nat. Mater. 19, 1312-1318 (2020)
  * [34] A. M. Jones, et al., Nat. Phys. 12, 323-327 (2016)
  * [35] Z. D. Zhang, K. Wang, Z. H. Yi, M. S. Zubairy, M. O. Scully and S. Mukamel, J. Phys. Chem. Lett. 10, 4448-4454 (2019)
  * [36] D. W. Wang, R. B. Liu, S.-Y. Zhu and M. O. Scully, Phys. Rev. Lett. 114, 043602-043606 (2015)
  * [37] Z. D. Zhang, P. Saurabh, K. E. Dorfman, A. Debnath and S. Mukamel, J. Chem. Phys. 148, 074302-074314 (2018)
  * [38] S. Asban, K. E. Dorfman and S. Mukamel, Proc. Natl. Sci. Acad. U.S.A. 116, 11673-11678 (2019)
  * [39] L. Ye, J. R. Rouxel, D. Cho and S. Mukamel, Proc. Natl. Sci. Acad. U.S.A. 116, 395-400 (2019)
  * [40] S.-N. Zhu, Y.-Y. Zhu, Y.-Q. Qin, H.-F. Wang, C.-Z. Ge and N.-B. Ming, Phys. Rev. Lett. 78, 2752-2755 (1997)
  * [41] A. Halpin, et al., Nat. Chem. 6, 196-201 (2014)
  * [42] K. W. Stone, K. Gundogdu, D. B. Turner, X. Li, S. T. Cundiff, K. A. Nelson, Science 324, 1169-1173 (2009)
  * [43] B. E. A. Saleh, B. M. Jost, H.-B. Fei and M. C. Teich, Phys. Rev. Lett. 80, 3483-3486 (1998)
  * [44] P. Wang, C. Chen and R. B. Liu, Chin. Phys. Lett. 38, 010301-010307 (2021)
  * [45] M. O. Scully, et al., Proc. Natl. Sci. Acad. U.S.A. 99, 10994-11001 (2002)
  * [46] G. I. Petrov, et al., Proc. Natl. Sci. Acad. U.S.A. 104, 7776-7779 (2007)
  * [47] D. Pestov, et al., Science 316, 265-268 (2007)
  * [48] K. E. Dorfman, F. Schlawin and S. Mukamel, J. Phys. Chem. Lett. 5, 2843-2849 (2014)
  * [49] A. J. Musser, et al., Nat. Phys. 11, 352-357 (2015)
  * [50] M. Seel and W. Domcke, J. Chem. Phys. 95, 7806-7822 (1991)
  * [51] Z. D. Zhang and J. Wang, Sci. Rep. 6, 37629-37637 (2016)
  * [52] J. A. Campos-Gonzalez-Angulo, R. F. Ribeiro and J. Yuen-Zhou, Nat. Commun. 10, 4685-4692 (2019)
  * [53] G. S. Engel, et al., Nature 446, 782-786 (2007)
  * [54] M. Kowalewski, K. Bennett, K. E. Dorfman and S. Mukamel, Phys. Rev. Lett. 115, 193003-193007 (2015)
  * [55] See Supplemental Material for a detailed description of the coherent Raman signals with entangled photons.
  * [56] J. R. Rouxel, M. Kowalewski, K. Bennett and S. Mukamel, Phys. Rev. Lett. 120, 243902-243906 (2018)
  * [57] P. F. Barbara, T. J. Meyer and M. A. Ratner, J. Phys. Chem. 100, 13148-13168 (1996)