# FASTER-CE: Fast, Sparse, Transparent, and Robust Counterfactual Explanations

Shubham Sharma
University of Texas at Austin
<EMAIL_ADDRESS>
&Alan H. Gee
Amira Learning
<EMAIL_ADDRESS>
&Jette Henderson
CognitiveScale
<EMAIL_ADDRESS>
&Joydeep Ghosh
University of Texas at Austin
<EMAIL_ADDRESS>

###### Abstract

Counterfactual explanations have substantially increased in popularity in the
past few years as a useful human-centric way of understanding individual
black-box model predictions. While several properties desired of high-quality
counterfactuals have been identified in the literature, three crucial
concerns: the speed of explanation generation, robustness/sensitivity and
succinctness of explanations (sparsity) have been relatively unexplored. In
this paper, we present FASTER-CE: a novel set of algorithms to generate fast,
sparse, and robust counterfactual explanations. The key idea is to efficiently
find promising search directions for counterfactuals in a latent space that is
specified via an autoencoder. These directions are determined based on
gradients with respect to each of the original input features as well as of
the target, as estimated in the latent space. The ability to quickly examine
combinations of the most promising gradient directions as well as to
incorporate additional user-defined constraints allows us to generate multiple
counterfactual explanations that are sparse, realistic, and robust to input
manipulations. Through experiments on three datasets of varied complexities,
we show that FASTER-CE is not only much faster than other state of the art
methods for generating multiple explanations but also is significantly
superior when considering a larger set of desirable (and often conflicting)
properties. Specifically we present results across multiple performance
metrics: sparsity, proximity, validity, speed of generation, and the
robustness of explanations, to highlight the capabilities of the FASTER-CE
family.

##  1 Introduction

Figure 1: The FASTER-CE method: An autoencoder is trained on the input
dataset. Then, the latent space is extensively sampled to generate random
latent space vectors that correspond to random decoder outputs, and for every
random perturbation, the value of every feature from the decoder output, and
the black-box prediction for the decoder output are noted. For every feature
in the dataset and for the black-box prediction, a linear model is trained
with the latent space vectors as inputs to the linear model and the feature or
the black-box prediction as the output. After finding these linear hyperplanes
associated with each feature and the black-box prediction, user-defined
constraints are incorporated and a combination of QR factorization with
respect to the normal vectors of the hyperplanes and finding interesections
between hyperplanes gives us a direction to move to find candidate
explanations. Multiple candidates are found for the test point for which we
need an explanation and checked for validity and post-process, and a set of
final explanations is returned.

As machine learning models get deployed in high-stake applications such as
finance and healthcare, providing human-centric explanations for the decisions
made by such models is becoming increasingly important. Not surprisingly, a
wide variety of explanation types with associated algorithms have been
proposed and deployed recently [towards], including feature attribution
methods, prototypes and counterfactuals. Among these, counterfactual
explanations have found increased traction as a useful, human-centric way of
providing individual-level explanations for the decision of a black-box model.
Originally introduced in [wachter2017counterfactual], counterfactual
explanations attempt to answer the question: "What would have to change in the
input to change the outcome of the model?". Such an explanation can be used as
a means of actionable recourse for an end-user subjected to a model decision.
The original definition [wachter2017counterfactual] focussed on finding a
valid counterfactual as close to the input as possible, where "valid" means
that the class label given by the model is different for the counterfactual.
However it was soon apparent that useful counterfactuals exhibit several other
desirable characteristics as well [verma2020counterfactual,
mothilal2020explaining, sharma2020certifai]. The sparsity of an explanation,
i.e., the number of attributes/features that change to obtain the
counterfactual, is one such property. For an end-user, changing just one
attribute may be easier and more understandable compared to changing multiple
features in order to receive a different outcome. However, including a
sparsity constraint in existing methods that provide black-box counterfactual
explanations as an objective is fairly non-trivial, since many of these
methods rely on gradient descent [wachter2017counterfactual,
mothilal2020explaining, balasubramanian2020latent] and an L0 penalty is non-
differentiable, so surrogates may be needed. Another key concern is the
robustness of counterfactual explanations. Recently, it has been shown that
counterfactual explanations are vulnerable to input perturbations as well as
adversarial attacks [slack2021counterfactual]. The realism and actionability
of an explanation is also desirable. Counterfactuals that are "in-
distribution" are preferred, so are those where the indicated change is
feasible. For example, a counterfactual that states that one needs to decrease
her age to get a more preferred outcome, indicates an infeasible change and
hence is not an actionable means of recourse for that end-user.

However, all these desirable properties can be at odds with each other. For
example, the nearest point with a different outcome may be in an infeasible
change direction. Moreover, given multiple counterfactuals for a given input,
the most actionable choice for recourse can be subjective and thus should be
an end-user’s choice. We have also come across applications such as automated
approval for medical procedures, where many decisions are being made each hour
and explanations are needed for each one. In such cases, the speed of
explanation generation becomes crucial because of the (near) real-time
business requirements, even though it may cause one to miss some
counterfactuals that could have been more actionable.

In this paper, we propose FASTER-CE: FAst, Sparse, TransparEnt, and Robust
Counterfactual Explanations. The FASTER-CE method is shown in Figure 1. First,
an autoencoder is trained on the input dataset. Then, the latent space of this
autoencoder is sampled from to generate a set of latent vectors
$\\{\tilde{\textbf{z}}\\}$ and corresponding decoder outputs
$\\{\textbf{x'}\\}$. Then, we train a set of linear models mapping the latent
vectors to each individual input feature now considered as the target, and one
more linear models that uses the black-box predictions $y^{\prime}$ of the
decoder output as target. Furthermore, using QR factorization to incorporate
user-defined constraints and by finding the intersection of hyperplanes
defined by the linear models, we identify promising direction(s) of
counterfactual search and use line search to generate a set of candidate
counterfactual explanations for a given test data point $\textbf{x}_{test}$.
The validity of these explanations is checked, post-processing is carried out
if need be to induce sparsity, and the set of valid counterfactual
explanations is returned.

There are multiple advantages of using this process to generate counterfactual
explanations. A reasonably accurate autoencoder provides a lower dimensional
representation of the data in which the search for counterfactuals is quicker,
more realistic and robust. Moreover, training the autoencoder and the linear
models is a one-time process. Finding counterfactual explanations for any
future input then is simply based on a projection of its latent representation
onto the intersection of hyperplanes and moving along specific directions.
This is significantly faster than performing gradient descent
[mothilal2020explaining], using an evolutionary algorithm
[sharma2020certifai], or the many existing methods that require an
optimization procedure for every test point. Through QR factorization, we are
able to generate the sparsest explanations, i.e., where only one feature
changes (when feasible). We can also generate multiple explanations
efficiently with multiple features changing, if required. Finally, the use of
linear models makes it easier to specify a margin around the boundary in which
an explanation should not exist, thereby allowing us to produce explanations
that have more robust recourse, as defined as in [slack2021counterfactual].

We propose three different variants of the FASTER-CE algorithm to generate
explanations reflecting different tradeoffs involving proximity, sparsity, and
realism. We theoretically motivate how our method produces more robust
explanations. Experiments on three different datasets demonstrate the
effectiveness of our approach. We are able to generate multiple explanations
with different objectives much faster than existing methods. We also discuss
the trade-offs between these objectives experimentally. To the best of our
knowledge, this is the first work providing a set of algorithms to generate
multiple explanations along key the desirable axes: robustness, sparsity,
realism, validity, proximity, and speed.

##  2 Related Work

Counterfactual explanations were originally defined and motivated as a useful
means of explaining artificial intelligence models in
[wachter2017counterfactual]. Gradient descent in the input space is used in
[wachter2017counterfactual] to generate explanations. Subsequently, numerous
other methods [mothilal2020explaining, sharma2020certifai,
balakrishnan2021towards, barr2021counterfactual, poyiadzi2020face,
guidotti2019factual, looveren2021interpretable, joshi2018xgems, keane2021if,
pawelczyk2020learning, schleich2021geco] have been proposed to generate these
explanations. A recent survey [stepin2021survey] provides details on many such
algorithms. However, several challenges exist towards the generation of these
explanations [sokol2019counterfactual, barocas2020hidden,
artelt2019computation]. For brevity, we only mention the methods most relevant
to our work in this section.

CERTIFAI [sharma2020certifai] provides multiple desirable objectives including
proximity, validity, sparsity, and realism of explanations in the generation
of counterfactuals. The method uses a genetic algorithm to generate candidate
counterfactuals. The use of an evolutionary algorithm allows one to more
easily incorporate a variety of objectives or domain constraints compared to
gradient based optimization methods, but it can be quite slow, specially for
large, high-dimensional data. DiCE is notable for incorporating a diversity
constraint to provide multiple counterfactual explanations, but is also slow,
and using gradient descent results in explanations that might not be robust
[slack2021counterfactual]. Latent-CF [balasubramanian2020latent] and
Sharpshooter [barr2021counterfactual] are two recent approaches that utilize a
latent space (like FASTER-CE) to find counterfactual explanations. Both
methods demonstrate the usefulness of representation learning. However, the
gradient descent based Latent-CF is vulnerable to providing explanations that
are not robust. Sharpshooter is a fast algorithm, but does not take into
account the sparsity or robustness of explanations. Finally,
balakrishnan2021towards use latent space perturbations to analyze the bias in
face recognition algorithms. The method uses such perturbations to generate
images for bias analysis, and not towards counterfactual explanations.

##  3 Theory

Given an input $\vec{x}$, a classifier $\vec{f}$, and a distance metric $d$, a
counterfactual explanation $\vec{c}$ can be found by solving the optimization
problem:

$\min_{\vec{c}}d(\vec{x},\vec{c})~{}~{}~{}\textrm{s.t.}\;\vec{f}(\vec{x})\neq\vec{f}(\vec{c})$
(1)

However, the solution to such an optimization problem is non-trivial,
especially when the classifier is assumed to be a black-box. Moreover, multi-
objective extensions for this formulation (where the objectives are beyond
proximity and validity) are needed to generate sparse, realistic, and robust
explanations. In this section, we describe our method to generate
counterfactual explanations. FASTER-CE utilizes the same underlying method and
provides three algorithms, each of which has different primary objectives: 1)
generating counterfactual explanations with only distance to the input as
consideration (where distance is defined based on Euclidean distance for
continuous features and simple matching distance for categorical features), 2)
generating sparse counterfactual explanations with sparsity as the main
objective and distance as a secondary objective, and 3) generating realistic
counterfactual explanations with sparsity and distance as secondary
objectives. A user can choose a specific version of FASTER-CE that best
matches her objectives, or simply deploy all three algorithms to generate
multiple explanations, since our proposed method is extremely quick.
Experimentally, we show that running all three algorithms to generate multiple
diverse explanations is still much faster than existing methods that generate
just one explanation.

Consider an autoencoder with an encoder $f_{\phi}$ and a decoder $g_{\theta}$
(Figure 1) that maps input x into a latent space vector z using the encoder.
The decoder maps the latent space vector to a reconstructed input x’. The
autoencoder is trained to minimize the distance (d) based loss =
$d(\textbf{x},\;g_{\theta}(f_{\phi}(\textbf{x})))$, The distance for
continuous features is the normalized Euclidean distance, and for categorical
features is one minus the cosine similarity between one-hot encoded vectors.

We first randomly sample the latent space of the trained autoencoder
extensively, generate decoder outputs for each of these random samples, and
query the black-box model for these decoder outputs. We then calculate the
principal axes of variation for the attributes and for the model output.
Formally, let there be a list of $n$ attributes (features in the input
dataset). We generate a training set $T_{z}$ =
$\\{z_{i},x^{\prime}_{i}\\}_{i=1}^{k}$ , where $x^{\prime}_{i}$ is the decoder
output for latent vector $z_{i}$. The value for attribute $j$ may be
continuous i.e. $\in$ [0, 1] or binary i.e. $\in$ 0, 1. We produce $T_{z}$ as
follows: First, we sample a large number of values of $z_{i}$ from the latent
space of the autoencoder. We obtain the corresponding $x^{\prime}_{i}$ by
passing these perturbations through the decoder. For each feature $j$, we use
$T_{z}$ to find a $T-1$ dimensional hyperplane $h_{j}=(n_{j},b_{j})$, where
$T$ is the dimension of the $T_{z}$, $n_{j}$ is the normal vector and $b_{j}$
is the offset. To find these hyperplanes, for continuous features, we train a
lasso regression model. For categorical features we train a support vector
machine (SVM) classifier with a linear kernel. Simultaneously, we find the
prediction of the black-box model $y^{\prime}$ for each of the decoder outputs
corresponding to every randomly generated latent space vector, and train an
SVM to find the prediction hyperplane
$h_{y^{\prime}}=(n_{y^{\prime}},b_{y^{\prime}})$. The hyperplane $h_{j}$
(approximately) represents a contour of constant values of feature $j$, and
the normal vector $n_{j}$ specifies the gradient direction for the feature.
Similarly, the hyperplane $h_{y^{\prime}}$ specifies (approximately) a contour
of constant values for predictions $y^{\prime}$, and the normal vector
$n_{y^{\prime}}$ specifies the gradient direction for the prediction.

The feature and prediction hyperplanes are then used to find counterfactual
explanations. The method to find counterfactual explanations for a point
$x_{test}$ without any consideration for sparsity or user-constraints only
depends on the prediction hyperplane and is shown in Algorithm 1. The test
point $x_{test}$ is fed into the encoder to generate the latent space vector
$z_{test}$. Then, $z_{test}$ is projected onto the hyperplane $h_{y^{\prime}}$
to get $z_{proj}$. Finally, we use line-search in the direction
$n_{y^{\prime}}$ or $-n_{y^{\prime}}$ (the direction opposite to the direction
of projection) in small incremental steps and query the black-box model until
a valid counterfactual explanation is found. Linearly increasing increments
are used, but multiplicative updates can also be considered. Choosing between
$n_{y^{\prime}}$ and $-n_{y^{\prime}}$ depends on whether the prediction of
the black-box for the test input point matches the prediction of the linear
prediction hyperplane. Specifically, if the predictions match, then the
decision of the linear hyperplane aligns with the decision of the black-box,
and the counterfactual search direction will be in the direction of
projection. Otherwise, we move in a direction opposite to the direction of
projection until a valid counterfactual explanation is found. Even if the
global linear model is unable to approximate well for a highly non-linear
black-box model, the counterfactual explanation is likely to remain valid for
the original black-box.

Algorithm 1 FASTER-CE1: Counterfactual Explanation

Inputs: Test point $x_{test}$. Black-box model $F(x)$, Learned Encoder
$f_{\phi}$, Learned Decoder $g_{\theta}$, Learned prediction hyperplane
$h_{y^{\prime}}=(n_{y^{\prime}},b_{y^{\prime}})$

Hyperparameters: $\epsilon$ the incremental step value , $\delta\epsilon$,
$\epsilon_{max}$

$z_{test}=f_{\phi}(x_{test})$

$z_{proj}$ = Projection of $z_{test}$ on hyperplane $h_{y^{\prime}}$

Initialize $\epsilon$ = 0, $c(\epsilon)$ = $\\{\\}$ the variable that stores
the counterfactual explanations

if prediction($n_{y^{\prime}}.f_{\phi}(x_{test})+b_{y^{\prime}})=F(x_{test})$
then

$s=sign(\frac{z_{proj}}{||z_{proj}||})$ Signed direction vector of projection,
direction of counterfactual is same

else

$s=-sign(\frac{z_{proj}}{||z_{proj}||})$ If prediction through prediction
hyperplane and prediction of input do not match, move in opposite direction to
projection

end if

while $c(\epsilon)$ is not empty AND $\epsilon$ $\leq$ $\epsilon_{max}$ do

$z_{proj}=z_{proj}+\epsilon.s.\frac{n_{y^{\prime}}}{||n_{y^{\prime}}||}$

Append $c(\epsilon)$ with $g_{\theta}(z_{proj})$

if $F(c(\epsilon))=F(x_{test})$ then

$c(\epsilon)$ = $\\{\\}$

end if

$\epsilon$ = $\epsilon$ \+ $\delta\epsilon$

end while

###  3.1 Sparse Counterfactual Explanations

Algorithm 2 shows the method to produce sparse explanations using the feature
hyperplanes along with the prediction hyperplane. For simplicity, we show an
example where a sparse explanation is obtained by changing only feature $a$.
We project the test point’s latent space vector $z_{test}$ onto the
intersection between the prediction hyperplane $h_{y^{\prime}}$ and the
feature hyperplane $h_{a}$. To do this, we iteratively project $z_{test}$ onto
each of these two hyperplanes till the intersection is found. Then, we find
new direction vectors $c_{j}$ for every feature $j$ such that they are
orthogonal to each other (called orthogonalization). Then, to only change
feature $a$, we start from the intersection of the hyperplanes
$h_{y^{\prime}}$ and $h_{a}$ and incrementally move in direction $c_{a}$ or
$-c_{a}$ (dependent on the test points black-box prediction and prediction
through the prediction hyperplane), repeatedly querying the black-box model
until a valid explanation is found.

Even though we ensure sparsity in the latent space, the decoded output may
still have more than one feature change. Such an explanation would be more
realistic since the autoencoder has learned the distribution of the dataset.
However, if sparsity is an over-riding objective i.e. only one feature is
required to be changed, we post-process to try to induce it. Specifically, if
more than one feature has changed between the input and counterfactual (i.e.
the L0 distance between $c(\epsilon)$ and $x_{test}$ is greater than 1), we
set all features except for feature $a$ (the feature required to be sparse)
between the input and counterfactual to be equal, and check if this is a valid
counterfactual. If it is, this is the final sparse counterfactual, otherwise
we perturb $c(\epsilon)$ randomly (in the input space) and check if we can
find a valid sparse explanation. This post-processing comes at the cost of
speed and potentially the realism of the counterfactual explanation. We
discuss this trade-off in the experiments section.

Algorithm 2 FASTER-CE2: Sparse Counterfactual Explanations

Inputs: Test point $x_{test}$. Black-box model $F(x)$, Learned Encoder
$f_{\phi}$, Learned Decoder $g_{\theta}$, Learned prediction hyperplane
$h_{y^{\prime}}=(n_{y^{\prime}},b_{y^{\prime}})$, Learned feature hyperplanes
$h_{j}=(n_{j},b_{j})$

Hyperparameters: $\epsilon$ the incremental step value , $\delta\epsilon$,
$\epsilon_{max}$

$z_{test}=f_{\phi}(x_{test})$

$z_{proj}$ = Projection of $z_{test}$ on intersection of hyperplanes
$h_{y^{\prime}}$ and $h_{a}$ // $a$ is one specific feature

Start orthogonalization of feature hyperplane direction vectors: Q,R <\- QR
Factorization of matrix $[n_{1},n_{2},...,n_{j},...n_{n}]$

for $i$ in $\\{1,...,n\\}$ do

$c_{i}=n_{i}$

for $k$ in $\\{1,...,n\\}$ do

if $i\neq k$ then

$c_{i}=c_{i}-\frac{Q_{k}.<Q_{k},c_{k}>}{<Q_{k}.Q_{k}>}$

end if

end for

end for

Initialize $\epsilon$ = 0, $c(\epsilon)$ = $\\{\\}$ the variable that stores
the counterfactual explanations

if prediction($n_{y^{\prime}}.f_{\phi}(x_{test})+b_{y^{\prime}})=F(x_{test})$
then

$s=sign(\frac{z_{proj}}{||z_{proj}||})$ Signed direction vector of projection,
direction of counterfactual is same

else

$s=-sign(\frac{z_{proj}}{||z_{proj}||})$ If prediction through prediction
hyperplane and prediction of input do not match, move in opposite direction to
projection

end if

while $c(\epsilon)$ is not empty do

$z_{proj}=z_{proj}+\epsilon.s.\frac{c_{a}}{||c_{a}||}$

Append $c(\epsilon)$ with $g_{\theta}(z_{proj})$

if $F(c(\epsilon))=F(x_{test})$ then

$c(\epsilon)$ = $\\{\\}$

end if

$\epsilon$ = $\epsilon$ \+ $\delta\epsilon$

end while

Post-processing

if $||c(\epsilon)-x_{test}||_{0}>1$ then

$c(\epsilon)=x_{test}$ for all features except feature $a$

if $F(c(\epsilon))=F(x_{test})$ then

Perturb feature $a$ till $F(c(\epsilon))\neq F(x_{test})$

if $F(c(\epsilon))=F(x_{test})$ then

No valid counterfactual found

else

Return $c(\epsilon)$

end if

end if

end if

###  3.2 Realistic and Multiple Counterfactual Explanations

One major concern with Algorithm 2 is that a valid counterfactual that
involves changing just one feature, may not exist. Moreover, if highly
correlated features are not changed as well, then the explanation may not be
realistic. Since the autoencoder’s latent space learns the distribution of the
data, using perturbations from it can help in producing in-distribution
explanations and alleviate this issue.

To provide realistic explanations, user-defined constraints can be
incorporated and QR factorization is carried out based on these constraints.
For example, if a user cannot substantially change income level without also
changing education status, the user can specify this as a constraint and the
FASTER-CE algorithm considers orthogonlization via QR factorization of the
normals of the hyperplanes only with respect to features that cannot change.
The intersection of hyperplanes is found for features that can change, and the
direction vector $c_{j}$ is only orthogonal to features that cannot change. We
show this in Algorithm 3, detailed in the supplementary material. Since
FASTER-CE is extremely quick (as shown in the experiments), running the method
repeatedly not only with the three proposed algorithms but also with different
sets of user constraints using Algorithm 3 can generate multiple explanations
with different feature changes and this is still faster than running most
other popular methods to generate just one explanation.

###  3.3 Robustness of Counterfactual Explanations

A notable definition of robustness in the counterfactual explanation
literature [dominguez2021adversarial] states that counterfactual explanations
(and hence recourse recommendations) should remain valid (i.e. lead to the
same outcome) for all plausible individuals (adequately) similar to the
individual seeking recourse. It has been shown in [dominguez2021adversarial]
that for a linear model with a set of actionable features that are unbounded,
the existence of robust recourse is guaranteed, and for a linear model, the
problem is equivalent to the standard robustness problem for a classifier:
$h^{\prime}(x)=\,<w,x>\,\geq b+||w^{*}||\epsilon$
Any recourse recommendation or counterfactual explanation would be
adversarially robust by a level of $\epsilon$ i.e. providing more robust
explanations is equivalent to having a higher value of $\epsilon$. Since
FASTER-CE1 uses a linear hyperplane to find counterfactual explanations and
does a line search in the direction of the normal to the prediction
hyperplane, we can generate a more robust means of recourse by making
incremental steps in the line-search beyond a counterfactual explanation that
may have been optimized based on proximity, sparsity, or realism. Hence, more
robust explanations would come at the cost of the distance between the input
and a valid counterfactual explanation: the closer the input and explanation
are, the closer the explanation is to the decision boundary. In the
experiments, we discuss the trade-off between producing more robust
explanations and the proximity objective.

##  4 Experiments and Results

Experiments are performed on three different datasets of varying complexity:
UCI Adult [kohavi1996scaling], German Credit [Dua:2019], and the Lending Club
datasets [bworld]. The UCI Adult dataset contains 48842 instances and 14
features, and the model is trained to predict if an individual’s income is
greater than or less than 50k. The German Credit dataset has 1000 instances
and 20 features, and the model is trained to predict an individual’s credit
risk (good or bad). The Lending Club dataset has 10000 instances and 55
features, and the model is trained to predict loan decisions based on whether
an individual will pay back their loan. Each dataset is used to train a multi-
layered perceptron neural network which is treated as the black-box model. We
consider a neural network as the black-box model since it achieves state of
the art accuracies across datasets, and to validate that the proposed approach
that approximates these non-linear models linearly can find good
counterfactual explanations. We use a standard autoencoder instead of using a
variational autoencoder (VAE) since a VAE assumes a gaussian distribution in
the latent space, and we did not want any such assumptions for tabular
datasets. An autoencoder is trained on each of the three datasets with a 70,
15, 15 split between training set, validation set, and test set. The details
of the architectures and performances of the neural networks and autoencoders
for each dataset are provided in the supplementary material. This section
first reports example counterfactual explanations for each of the FASTER-CE
algorithms. Thereafter, we compare our method in terms of various metrics for
counterfactual explanations to some state of the art methods and discuss the
trade-offs between these metrics. We then show the effect of producing more
robust explanations using our algorithm. All experimental results are averaged
across five runs and error bars are also provided.

###  4.1 Example Counterfactual Explanations

We show example counterfactual explanations for the UCI Adult dataset in Table
2 for a test point for all three algorithms. We only show the features that
changed between the input and counterfactual explanation. FASTER-CE2 is used
three times to generate three different types of sparse explanations
corresponding to three different features. FASTER-CE3 is used with the
constraint that the education and occupation features are the two features
that can change.

The results are shown in Table 2. The test point (input) received a negative
decision from the black-box model. When using FASTERCE-1, which only accounts
for distance from the prediction hyperplane, the occupation and capital gain
features change. When we induce sparsity with respect to a feature using
FASTER-CE2, only one feature changes. In the table, we show three examples of
just one feature changing. When we include a user-defined constraint that both
education and occupation can change, both features change in the resultant
counterfactual. For the person subject to the negative decision, any of these
decisions could be valuable. Example results for the German Credit dataset and
the Lending Club dataset and are provided in the supplementary material.

Table 1: Counterfactual Explanations for the UCI-Adult dataset for all three algorithms. Any feature that has changed with respect to the input is marked in bold Features | Education | Occupation | Capital Gain | Hours-per-week | Prediction
---|---|---|---|---|---
Input | Assoc-acdm | Professional | 0 | 40 | <=50k
FASTER-CE1 | Assoc-acdm | White-collar | 33704 | 40 | >50k
FASTER-CE2 | Doctorate | Professional | 0 | 40 | >50k
FASTER-CE2 | Assoc-acdm | Professional | 99985 | 40 | >50k
FASTER-CE2 | Assoc-acdm | Professional | 0 | 76 | >50k
FASTER-CE3 | Doctorate | White-collar | 0 | 40 | >50k

###  4.2 FASTER-CE performance and comparison to other methods

We compare the three proposed algorithms to three other popular methods that
provide counterfactual explanations for all three datasets. All evaluations
are carried out on the test sets (randomly chosen) for each of the three
datasets. We compare to DiCE [mothilal2020explaining], CERTIFAI
[sharma2020certifai], Latent-CF [balasubramanian2020latent], and
[looveren2021interpretable] (we call this method ICEGP). While CERTIFAI, DiCE,
and ICEGP are three of the most widely used techniques to generate
counterfactual explanations, Latent-CF is a recent method that uses gradient
descent in the latent space to find counterfactual explanations. Latent-CF
requires access to predicted class probabilities and not just the class label,
hence is not truly a black-box method.

Figure 2: Comparison of methods to FASTER-CE for the UCI Adult dataset, German
Credit Dataset, and Lending Club dataset. Validity and sparsity are in
percentages. Higher validity is better; sparsity value is lower if less
features change and hence a lower value is better; the measure for proximity
is just the distance between the input and counterfactual explanation and
hence a lower proximity is better, and lower runtime per explanation is
better. Time is measured in seconds.

We use validity, sparsity, proximity, and efficiency to evaluate the
counterfactual explanations found. For validity, we report the percentage of
counterfactual explanations that are valid. For sparsity, we compute the
average number of features that change between the input and counterfactual
explanation for all test points (hence, lower is better). For proximity, the
distance between the input and counterfactual explanation is measured by
taking the normalized Euclidean distance (normalized by the median absolute
deviation) for continuous features and simple matching distance for
categorical features, as in [sharma2020certifai]. We report the average
distance between the input and counterfactual explanations for all input
(test) points. Efficiency is reported via two different types of runtime:
total runtime for the algorithm across all test points when the algorithm is
run from start to end, and average time per counterfactual explanation
generation when the back-bone for any algorithm is already setup. For example,
for FASTER-CE, training the autoencoder and training the linear models is a
fixed cost with respect to time, and hence in deployment, only the three
algorithms as defined in the theory are run. Total runtime results are in the
supplementary material. For FASTER-CE2, we average the results that we get by
trying to induce sparsity across all features. For FASTER-CE3, we induce the
same feature constraints across the test set based on characteristics of the
dataset (e.g. gender or race cannot be changed for the UCI Adult dataset). The
results for the comparison for all three datasets are shown in Figure 3 for
the UCI Adult, the German credit, and the Lending Club datasets.

Table 2: Robustness of recourse (in percentage) measured in terms of counterfactual explanations that still have the same prediction after perturbing the input for different values of $\delta\epsilon$ and proximity of counterfactual explanations for the lending club dataset. Metric | Robustness | Proximity
---|---|---
FASTER-CE1 $\delta\epsilon=0.05$ | 86.59 $\pm$ 2.54 | 1.29 $\pm$ 0.22
FASTER-CE1 $\delta\epsilon=0.1$ | 88.70 $\pm$ 3.19 | 1.38 $\pm$ 0.18
FASTER-CE1 $\delta\epsilon=0.3$ | 100 $\pm$ 0 | 4.08 $\pm$ 0.34
FASTER-CE1 $\delta\epsilon=0.5$ | 100 $\pm$ 0 | 6.41 $\pm$ 0.28
FASTER-CE1 $\delta\epsilon=1$ | 100 $\pm$ 0 | 13.43 $\pm$ 0.67
DiCE | 93.41 $\pm$ 2.55 | 3.12 $\pm$ 0.27
CERTIFAI | 89.72 $\pm$ 2.93 | 1.88 $\pm$ 0.36
Latent-CF | 90.51 $\pm$ 3.21 | 2.24 $\pm$ 0.20

The validity of explanations is 100% for all methods except FASTER-CE2. This
is because the severe restriction of changing just one feature may not yield a
valid counterfactual explanation. This is also why the time taken for FASTER-
CE2 is higher: trying to find a valid explanation takes longer. However,
explanations produced by FASTER-CE2 are more sparse compared to any of the
other methods across all three datasets.

CERTIFAI performs well in terms of proximity for datasets with lesser number
of features (UCI Adult and German Credit) and FASTER-CE1 performs comparably.
However, with a larger dimensional dataset (Lending Club), FASTER-CE1 performs
better. This is because CERTIFAI uses a genetic algorithm that does not scale
well as the dimensionality of the dataset grows when the number of generations
for the algorithm is restricted. Our method relies on projections from the
lower dimensional latent space, and hence is able to find explanations that
are closer much faster.

Across all datasets, FASTER-CE is much faster for runtime per explanation,
which is crucial for explainability methods in deployment. Across datasets,
FASTER-CE2 and FASTER-CE3 are able to find the intersection between
hyperplanes with a maximum of 20 iterations. FASTER-CE3 is able to provide a
reasonable trade-off between sparsity, validity, and the time taken to
generate explanations. The choice of which algorithm to use is up to the end-
user based on explanation needs. Since the time taken for generating an
individual explanation across all FASTER-CE algorithms is still significantly
less than producing one explanation for other methods for a larger dataset, a
user can run all algorithms and choose the means of recourse that is best for
them.

Robust Explanations. To examine the robustness of FASTER-CE based explanations
we randomly perturb the features in the input that changed between the input
and counterfactual explanation found using FASTER-CE1, and find the
corresponding counterfactual explanation using FASTER-CE1. We report the
proximity and robustness of recourse (measured by checking if the recourse
remains valid) for the counterfactual explanation (a robust recourse would
keep the counterfactual explanation valid after perturbation). We do this for
different values of $\delta\epsilon$, and report the results for the lending
club dataset, and compare to DiCE, CERTIFAI and Latent-CF. The results are in
Table 5. FASTER-CE is able to produce explanations that are more robust to
input manipulations. This comes at the cost of proximity: producing more
robust explanations results in the input and counterfactual explanation being
further away. However, the proximity of explanations is not much compromised
and in practice one may prefer the robust explanation version for this reason.

##  5 Conclusion and Future Work

This paper introduces FASTER-CE, a novel set of algorithms to generate
counterfactual explanations for a black-box model based on multiple quality
objectives. This work was motivated by actual business applications that
demanded much faster generation of explanations that were likely actionable
based on multiple criteria - a need that was not met by any existing
algorithm. FASTER-CE can produce sparse, realistic, and robust explanations
efficiently. For future work, we would like to examine the tradeoff in
building multiple (locally) linear models for different points as opposed to
training one global linear model per target variable. This will take more
time, but may produce even better counterfactuals. Further, we would be
interested in studying the effect of adding a realism constraint beyond
providing user-defined constraints, such as including an objective that tries
to optimize the distance between a counterfactual explanation and class
prototypes in the z-space. It would also be interesting to investigate the
variability in explanations with drift in data. While the total runtime for
FASTER-CE is still less than existing methods (which also do not consider
drift) and hence it would still be faster for new data, we would like to
extend the training procedure to account for such streaming data.

##  Appendix A Supplementary Material

Additional experiments and results, along with a discussion on the method is
provided in this section.

###  A.1 Realistic and Multiple Counterfactual Explanations

One major concern with Algorithm 2 is that a valid counterfactual that
involves changing just one feature, may not exist. Moreover, if highly
correlated features are not changed as well, then the explanation may not be
realistic. Since the autoencoder’s latent space learns the distribution of the
data, using perturbations from it can help in producing in-distribution
explanations and alleviate this issue.

To provide realistic explanations, user-defined constraints can be
incorporated and QR factorization is carried out based on these constraints.
For example, if a user cannot substantially change income level without also
changing education status, the user can specify this as a constraint and the
FASTER-CE algorithm considers orthogonlization via QR factorization of the
normals of the hyperplanes only with respect to features that cannot change.
The intersection of hyperplanes is found for features that can change, and the
direction vector $c_{j}$ is only orthogonal to features that cannot change. We
show this in Algorithm 3 for a case when the constraint is based on two
features, but the method can easily be extended for any such constraints.
Moreover, if the values associated with any feature need to be within a range,
the corresponding latent space vectors can be generated by random sampling
such that latent vectors that have decoder outputs beyond this range are not
considered to train the linear model. If user-defined constraints are not
available, this algorithm can be run repeatedly for different combinations of
feature changes to generate multiple diverse counterfactual explanations.
Since FASTER-CE is extremely quick, running the method repeatedly is still
faster than running other existing methods to generate one explanation (as
shown in the experiments).

Algorithm 3 FASTER-CE3: Realistic Counterfactual Explanations

Inputs: Test point $x_{test}$. Black-box model $F(x)$, Learned Encoder
$f_{\phi}$, Learned Decoder $g_{\theta}$, Learned prediction hyperplane
$h_{y^{\prime}}=(n_{y^{\prime}},b_{y^{\prime}})$, Learned feature hyperplanes
$h_{j}=(n_{j},b_{j})$, User-defined constraints: Only change features $1$ and
$2$

Follow Algorithm 2 till projection step

Start orthogonalization of feature hyperplane direction vectors: Q,R <\- QR
Factorization of matrix $[n_{1},n_{2},...,n_{j},...n_{n}]$

for $i$ in $\\{1,...,n\\}$ do

$c_{i}=n_{i}$ Adding user-defined constraints

for $k$ in $\\{1,...,n\\}$ do

if ($i=1$ AND $k=2$) OR ($i=2$ AND $k=1$) then

$c_{i}=n_{i}$

else if $i\neq k$ then

$c_{i}=c_{i}-\frac{Q_{k}.<Q_{k},c_{k}>}{<Q_{k}.Q_{k}>}$

end if

end for

end for

Initialize $\epsilon$ = 0, $c(\epsilon)$ = $\\{\\}$ the variable that stores
the counterfactual explanations

if prediction($n_{y^{\prime}}.f_{\phi}(x_{test})+b_{y^{\prime}})=F(x_{test})$
then

$s=sign(\frac{z_{proj}}{||z_{proj}||})$ Signed direction vector of projection,
direction of counterfactual is same

else

$s=-sign(\frac{z_{proj}}{||z_{proj}||})$ If prediction through prediction
hyperplane and prediction of input do not match, move in opposite direction to
projection

end if

while $c(\epsilon)$ is not empty do

$z_{proj}=z_{proj}+\epsilon.s.\frac{c_{j}}{||c_{j}||}$

Append $c(\epsilon)$ with $g_{\theta}(z_{proj})$

if $F(c(\epsilon))=F(x_{test})$ then

$c(\epsilon)$ = $\\{\\}$

end if

$\epsilon$ = $\epsilon$ \+ $\delta\epsilon$

end while

###  A.2 Autoencoder and Black-Box details

The details of the architectures and the accuracies of the neural networks and
autoencoders for each dataset are shown in Table 1. For every dataset, the
categorical features are one-hot encoded and the continuous features are
normalized. The architectures for both the neural networks and the
autoencoders for each dataset have been experimented upon and the values
associated with the best accuracy and best loss (lower is better),
respectively, are reported.

###  A.3 Example Counterfactual Explanations

Results similar to the UCI Adult dataset are observed for the German Credit
dataset and the Lending Club dataset, as shown in Table 4 and Table 5
respectively. Each of these explanations are valid counterfactuals, and can be
used as a means of recourse. Moreover, each of the explanations is close to
the input, and has different characteristics.

Table 3: Details on the autoencoder and black-box neural network models for the three datasets. NN is the neural network, AE is the autoencoder. Each number in the encoder architecture corresponds to the size of the layer (Decoder Architecture is just the reverse) AEloss: Reconstruction loss of the autoencoder (lower is better) Characteristics | NN layers | NN layer size | NN Accuracy | Encoder arch. | AE loss
---|---|---|---|---|---
UCI Adult | 2 | 30 | 83.1 | (42,30,20,10) | 0.14
German Credit | 2 | 40 | 75.74 | (58,40,30,10) | 0.09
Lending Club | 5 | 60 | 83.92 | (168,128,64,32,16) | 0.19
Table 4: Counterfactual Explanations for the German Credit dataset for all three algorithms. Any feature that has changed with respect to the input is marked in bold. For ease of understanding, categorical feature values have been changed in the . For status of checking account: Low - A11, Medium - A12, and High - A13 in the dataset. For payment status of previous credit: All credits paid back duly - A31, Delay - A33 Features | Status | Duration | Credit amount | Previous credit | Prediction
---|---|---|---|---|---
Input | Low | 36 | 6887 | Delay | Bad
FASTER-CE1 | Medium | 24 | 6887 | Paid back duly | Good
FASTER-CE2 | High | 36 | 6887 | Delay | Good
FASTER-CE2 | Low | 36 | 6887 | Paid back duly | Good
FASTER-CE3 | Medium | 16 | 6887 | Delay | Good
Table 5: Counterfactual Explanations for the Lending Club dataset for all three algorithms. Any feature that has changed with respect to the input is marked in bold Features | Interest Rate | FICO | Debt-To-Income Ratio | Prediction
---|---|---|---|---
Input | 0.1071 | 707 | 14.29 | Fully Paid
FASTER-CE1 | 0.1071 | 682 | 15.13 | Charged Off
FASTER-CE2 | 0.193 | 707 | 14.29 | Charged Off
FASTER-CE2 | 0.1071 | 614 | 14.29 | Charged Off
FASTER-CE3 | 0.1071 | 682 | 15.13 | Charged Off

###  A.4 Total runtime

The total runtime (Figure 3 across methods is comparable for lower dimensional
datasets. However, the total time for the higher dimensional dataset increases
drastically for other methods, as does the individual average time as shown in
the main paper.

Figure 3: Comparison of methods to FASTER-CE for the UCI Adult dataset, German
Credit Dataset, and Lending Club dataset. Time is measured in seconds.

###  A.5 Discussion

As shown in the experiments above, the three proposed algorithms can generate
a diverse set of explanations efficiently and these explanations are valid,
close to the input, sparse, realistic, and robust. However, there are trade-
offs among all these objectives, as demonstrated by the results in Figure 3
and Table 5. $\delta\epsilon$ controls the trade-off between not only
robustness and proximity, but also runtime. For the results in Figure 3, we
use $\delta\epsilon=0.1$ across datasets, and find that a valid counterfactual
explanation is returned in not more than 20 iterations in most cases for
FASTER-CE1 and FASTER-CE3. A larger value of $\delta\epsilon$ might lead to
finding a valid explanation quickly that is also more robust, compared to
using a smaller value, but the counterfactual explanation will be further away
from the input. Furthermore, as shown by FASTER-CE2, sparsity as a concrete
objective can lead to producing invalid counterfactuals. The runtime of the
algorithm to produce sparse explanations is also more than the other FASTER-CE
algorithms because of the post-processing step. Experimentally, we have found
that a dataset with more categorical features that have fewer categories (eg.
a binary protected attribute for gender) lead to more invalid explanations.
Moreover, the proximity of these explanations on average are also worse than
when the sparsity constraint is not considered. While we generate individually
realistic explanations by including user constraints, evaluating the realism
of an explanation quantitatively can be done at a global level by comparing
the distribution of the counterfactual explanations to the input dataset. We
refrain from doing this evaluation because a counterfactual explanation is
local in nature, and how realistic it is should be determined by an end-user
subject to the decision of the model.

We use the definition of robustness of recourse from
[dominguez2021adversarial], which defines robustness in terms of the validity
of recourse. Given this definition, we show that FASTER-CE algorithms can
easily provide more robust recourse than existing methods, and to the best of
our knowledge, FASTER-CE is the only existing counterfactual explanation
method that tries to provide robust recourse with respect to the definition in
[dominguez2021adversarial]. However, there can be other possible notions of
robustness. For example, similar individuals should achieve similar recourse
not only in terms of validity, but also proximity, sparsity, and realism.
Investigating robustness of recourse in terms of all these factors is an
important and open problem.

It is also interesting to note that this method may potentially lend itself
towards understanding the bias of a black-box model. Specifically, if the
linear hyperplane associated with a protected attribute like gender has a
normal vector that is perpendicular to the normal vector of the prediction
hyperplane, this means that the prediction is not impacted directly by the
protected attribute. Furthermore, if the normal vector of the protected
attribute hyperplane is orthogonal to all other features, this protected
attribute also doesn’t effect the values of other features in the dataset.
However, given that we are training an autoencoder to learn the distribution
of data, and then approximating the black-box with a linear hyperplane, such
an analysis needs to be approached with caution. Even though these
approximations don’t significantly effect the generation of a counterfactual
explanation owing to the conditional constraints to find the direction of
constant search and constant querying of the black-box model in our
algorithms, the fairness of a model might still differ between a linear
approximation of the model and the original black-box.


Control laws for continuous-time dynamical systems are most often implemented via digital controllers using a
sample-and-hold technique. Numerical discretization of the continuous system is an integral part of subsequent
analysis. Feedback linearizability of such sampled systems is dependent upon the choice of discretization map
or technique. In this article, for feedback linearizable continuous-time systems, we utilize the idea of retraction maps to construct discretizations that are feedback linearizable as well. We also propose a method to functionally compose discretizations to obtain higher-order integrators that are feedback linearizable.

§ INTRODUCTION

Digital controllers facilitate the implementation of continuous-time control systems via discretization. For non-autonomous systems i.e., for systems with inputs this is done via (a) sample and hold technique where the control input is held constant between two sampling intervals and (b) a discretization scheme that solves the evolution of the continuous-time dynamical systems numerically. Different numerical schemes result in different discretizations of the continuous time systems. On Euclidean spaces i.e., for systems evolving on $\R[n]$, some of the common numerical integration schemes are Euler Integrations methods, Runge-kutta-based methods, Simpson 1/3 rule, etc. <cit.>. While these schemes perform well for systems evolving in euclidean spaces, when implemented for systems evolving on general manifolds, they do not guarantee that the system states stay on the manifold. In order to maintain the non-euclidean structure of the underlying manifold one would like to construct integrators that respect the underlying geometry of the continuous-time dynamical system. Such integrators are called geometric integrators and these result in more accurate long-term behavior. A summary of geometric integrator schemes is given in <cit.>. Retraction maps  are a generalization of euclidean discretizations on non-euclidean manifolds (see <cit.>). Retraction maps allow us to construct geometric discretizations that guarantee the system states stay on the manifold.

Feedback linearization allows us to transform a nonlinear control system into a linear system via a coordinate transformation and invertible control feedback. This allows us to utilize the control design methods such as pole placement etc  available for linear systems to synthesize controls for the nonlinear system. A study of feedback linearization for continuous time systems is provided in <cit.> and references therein. A discrete-time equivalent of feedback linearization is studied in <cit.>. Sampling, in general, does not preserve the feedback linearization, i.e., given a feedback linearizable continuous time system, under sample and hold method and a particular choice of discretization the resulting discrete-time system need not be feedback linearizable (in discrete time) in general. In <cit.>, with help of an example shows that a continuous system, under one discretization, may result in a feedback linearizable discrete-time system while it may not do so under some other discretization. Under exact discretization methods, the feedback linearizability of the sampled time systems has been studied in <cit.>. Since feedback linearization allows us to utilize the advantages of  linear control theory, it is of interest to find (numerical) discretizations that are feedback linearizable.

Contribution: In this article, given a feedback linearization we utilize retraction maps to construct discretizations  that are feedback linearizable. We also provide a way to compose these discretizations to generate symmetric discretizations that are accurate upto second order while maintaining feedback linearizability. However, this requires multi-rate sampling.

Organization: The article is organized as follows: in Section <ref> and <ref>, we provide a brief introduction on retraction maps and continuous time feedback linearization respectively. The retraction maps are defined for autonomous systems, we extend this notion to nonautonomous systems. In section <ref> we provide our main result where we utilize the retractions maps to construct first-order discretizations  that are feedback linearizable and in Section <ref> we present a way to compose first-order discretizations to construct higher-order feedback linearizable schemes. This is obtained via multi-rate sampling schemes. We demonstrate these results on a simple example in Section <ref>.

§ RETRACTION AND  DISCRETIZATION MAPS

Let $M$ be an $n$ dimensional manifold and $TM$ be the associated tangent bundle. Let $TM\ni(x,v)\mapsto \tau_M(x,v)\coloneqq x$ be the canonical projection onto the manifold.

Let $U\subset TM$ be an open neighborhood of the zero section of the tangent bundle $TM$. $U\ni(x,v)\mapsto R(x,y)\in M\times M$ is a discretization map if it satisfies

* $(x,0_x)\mapsto R(x,0_x) = (x,x)$

* $T_{(x,0_x)}R^2-T_{(x,0_x)}R^1={\rm Id}_{T_xM} \colon T_{(x,0_x)}T_xM\simeq T_xM \lra T_xM$ is equal to the identity map on $T_xM$ for any $x$ in $M$.

As a consequence is easy to show that  any discretization map $R$ is a local diffeomorphism.

Let $R$ be a discretization on $M$. Consider the inversion map $(x,y)\ni M\times M \mapsto \I_M(x,y)\coloneqq (y,x)\in M\times M$. The adjoint of $R$ is defined by $U\ni(x,v)\mapsto R^*(x,v) \coloneqq \I_M(R(x,-v))$.

A discretization is called symmetric if $R= R^*$.

Given $X\in\mathfrak{X}$ a vector field on $M$ and a discretization map $R$ we have the following discretization map and a fixed time discretization map $t\mapsto(t-\alpha h, t+(1-\alpha)h)$, $\alpha\in[0,1]$.

The discretization of $X$ defined by

\[
R^{-1}(x_k, x_{k+1})=hX\underbrace{( \tau_M( R^{-1}(x_k, x_{k+1})) )}_{\in M}
\]

is  a first-order discretization of $X$ and second-order if $R$ is symmetric.

Let $M$ and $N$ be $n$ dimensional manifold and $M\ni x\mapsto \phi(x)\eqqcolon y\in N$ be a diffeomorphism. For a given discretization $R$ on $M$, $R_\phi := (\phi\times\phi)\circ R\circ T\phi^{-1}$ is a  discretization on $N$ (see Figure <ref>).

For any given $y\in N$ we have that

\begin{eqnarray*}
R_{\phi}(y,0_y)&=&\left( (\phi\times \phi)\circ  R\circ T\phi^{-1}\right)(y,0_y)\\
&=&\left((\phi\times \phi)\circ  R\circ T\phi^{-1}\right)(\phi(x),0_{\phi(x)})
\\
&=&(\phi\times \phi)^{-1}R(x, 0_x)\\
&=&(\phi\times \phi)^{-1}(x,x) = (y,y)
\end{eqnarray*}

this proves the first condition.

Now, given a vector $u_y\in T_yN$ we

\begin{eqnarray*}
& &\left(T_{(y,0_y)}R^2_{\phi}-T_{(y,0_y)}R^1_{\phi}\right)(y,u_y)\\
(\phi\circ R^1\circ T\phi^{-1})(y,su_y)\\&~-&	(\phi\circ R^2\circ T\phi^{-1})(y,su_y)
R^1(t (T\phi^{-1})(y,u_y)) \\~-	 R^2(t(T\phi^{-1})(y,u_y))
&=&T_{y} \phi((T\phi^{-1})(y,u_y) ) =(y,u_y)
\end{eqnarray*}

using the linearity of the map $T_y\phi$ and that $R$ is a discretization map.

Using the inversion map $\I_M$ one can easily show that $R^*_\phi = (\phi\times\phi)\circ R^*\circ T\phi^{-1}$ is the adjoint discretization of $R_\phi$. Further, $R_\phi$ is symmetric if $R$ is symmetric. From definition of $R$ and $R_\phi$, Figure <ref> commutes.

(m) [matrix of math nodes,row sep=3em,column sep=4em,minimum width=2em]

TM     TN

M×M     N×N
(m-1-1) edge node [left] $R$ (m-2-1)
edge node [above] $T\phi$ (m-1-2)
(m-2-1.east|-m-2-2) edge node [below] $\phi\times\phi$
(m-1-2) edge node [right] $R_{\phi}$ (m-2-2);

$R$ and $R_\phi$ commutes as shown above

§ CONTINUOUS TIME CONTROL SYSTEM

Let $M$ be an $m$ dimensional manifold and $U\subset\R[n]$ be open. For each $u\in U$ let $X(\cdot,u)\in \mathfrak{X}(M)$ be a vector field on $M$. Then for a fixed $T>0$, a continuous-time dynamical system on $M$ is given by

\begin{equation}
\label{ctssys}
\odv{}{t}x(t) = X(x(t),u(t)) ~ \forall~ t\in[0,T],
\end{equation}

with $t\mapsto x(t)\in M$ and $t\mapsto u(t)$ for all $t\in[0,T]$. A point $(x_0,u_0)\in M\times U$ is said to be an equilibrium point of (<ref>) if $X(x_0,u_0) = 0$.

§.§ Feedback Linearization of Continuous Time systems

Let $M$ and $N$ be two $n$-dimensional manifolds and $\phi:M\lra N$ be a diffeomorphism. Let $X\in \mathfrak{X}(M)$ be a vector field on $M$. Then $X_\phi \coloneqq T\phi\circ X\circ\phi^{-1}$ is a vector field on $N$. Further for the dynamical system

\begin{equation}
\label{ctssys1}
\odv{}{t}y(t) = X_\phi(y(t),u(t))~\forall t\in[0,T]
\end{equation}

with $y(0) = \phi(x(0))$ satisfy $y(t) = \phi(x(t))$  where $x(t)$ is a solution of (<ref>) with

Let $\Open{x_0}\ni x_0$ and $\Open{u_0}\ni u_0$ be open neighorhoods around $x_0$ and $u_0$ of $M$ and $U$, respectively. Let $\Open{x_0}\ni x\mapsto\phi(x)\coloneqq y\in N:= \R[n]$ be a diffeomorphism to its image and  $\Open{x_0}\times\Open{u_0}\ni (x,u)\mapsto\psi(x,u)\coloneqq v\in\R[m]$ such that for each fixed $x$, $\psi(x,\cdot):U\lra\R[n]$ is invertible.  A given continuous time system (<ref>) is said to be (locally)  feedback linearizable around $(x_0,u_0)$ on $\Open{x_0}\times\Open {y_0}$ if there exists matrices $A\in\R[n\times n]$ and $B\in\R[n\times m]$ such that $X_\phi (y,u) = Ay+Bv$
with $v = \psi(\phi^{-1}y,u)$. The feedback linearized dynamical system is given by

\begin{equation}
\label{ctslin}
\odv{}{t}y(t) = Ay(t)+Bv(t) ~\forall t\in [0,T].
\end{equation}

For background on feedback linearization, we refer the reader to <cit.> and references therein.

§ NUMERICAL DISCRETIZATION OF CONTINUOUS-TIME SYSTEMS

Continuous time control systems are implemented via digital controllers using sample and hold method where the control input $u$ is held constant at a fixed value between two successive samples i.e., $u(t) = u_k~\forall~t\in[t_k,t_k+h[~\forall k\in\N$ where $h$ is the fixed sampling period. Further, since analytical solutions for (<ref>) are often not available in closed form, the solutions are to be approximated numerically.

Let $U\subset\R[m]$ be open and for each $u\in U$, $X(\cdot,u)\in\mathfrak{X}(M)$ is a vector field on $x$. Let $R$ be a discretization map on $M$ then using Proposition <ref> a discretization of $X(\cdot,u)$ is defined by

\begin{equation}
\label{disc_sys_1}
R^{-1}(x_k,x_{k+1}) = hX(\tau_M(R^{-1}(x_k,x_{k+1})),u_k)
\end{equation}

where the control input $u$ is held constant over the interval $[t_k,t_{k+1}[$ i.e., $u(t) = u_k~\forall~ t\in[t_k,T_{k+1}]$.

A different choice for $R$ leads to different numerical discretization schemes. For example, on Euclidean spaces ($M=\R[m]$),  $R(x,v) = (x,x+v)$ results in the Explicit Euler discretization scheme  $x_{k+1}= x_k+hX(x_k,u_k)$, while $R^*(x,v) = (x-v,x)$ defines the Implicit Euler Discretization with $x_{k+1}-hX(x_{k+1},u_k) = x_k$. Solving (<ref>) implicitly for $x_{k+1}$, the sampled discrete-time system can be explicitly written as

\begin{equation}
\label{disc_sys}
x_{k+1} = F(x_k,u_k;h)
\end{equation}

where $x_k\in M$ and $u_k\in U$ for all $k\in\N$ and $M\times U\ni(x,u)\mapsto f(x,u)\in M$ is a smooth map (if $F$ is not well defined on entire $M$ one may very well work with a local definition of $F$, by replacing $M$ with an open neighborhood around $x_0$ in $M$). From the properties of retraction maps, one can show that at equilibrium point $(x_0,u_0)$ one has $F(x_0,u_0;h)= x_0$.

§.§ Feedback Linearization of Discrete time system

The idea of feedback linearization can be extended to discrete-time systems as follows. Consider the discrete-time system given by (<ref>).

Let $\Open{x_0}\ni x_0$ and $\Open{u_0}\ni u_0$ be open neighborhoods around $x_0$ and $u_0$. Let $\Open{x_0}\ni x\mapsto y \eqqcolon \phi(x) \in N\coloneqq \R[n]$ be a diffeomorphism onto its image, and $\Open{x_0}\times\Open{u_0}\ni(x,u)\mapsto v\eqqcolon\psi(x,u)\in\R[m]$ be such that for each $x$, $\psi(x,\cdot)$ is locally invertible.  The discrete-time system (<ref>) is said to be feedback linearizable if there exist matrices $A_h$ and $B_h$ such that

\begin{equation*}
\phi(F_h(x,u)) = A_h\phi(x) +B_h\psi(x,u) = A_hy +B_hv
\end{equation*}

The discrete-time system (<ref>) is linearized to

\begin{equation}
\label{disc_lin}
y_{k+1} = A_hy_k+B_hv_k
\end{equation}

The feedback linearizability of discrete-time systems has been dealt with in great detail in <cit.>. For sampled time continuous time system the feedback linearizability is in general not preserved, i.e., a feedback linearizable continuous time system when implemented with sample and hold may not result in a feedback linearizable discrete-time system. The linearizability is not only dependent upon the underlying continuous-time system but also on the choice of discretization (see <cit.>). Using this as our motivation we are interested in the following problem – given a feedback linearizable continuous time system (<ref>) is it possible to construct a numerical discretization (<ref>)
that is feedback linearizable?

§.§ Constructing feedback linearizable discretization maps

Consider the continuous time system given by (<ref>). Let $\phi$ and $\psi$ be as in Definition (<ref>). Suppose (<ref>) is feedback linearizable to (<ref>). Keeping $v(t) = v_k$ for all $t\in[t_k,t_{k+1}[$, let $R$ be a discretization for (<ref>) such that it preserves the linearity of (<ref>) i.e. it results in a discrete system

\begin{equation*}
y_{k+1} = A_hy_k+B_hv_k
\end{equation*}

where $A_h$ and $B_h$ are fixed matrices of appropriate order. Given a discretization map $R$ for (<ref>), using Proposition <ref> one can construct a discretization map

\begin{equation}
\label{4.7}
R_{\phi^{-1}} = (\phi\times\phi)^{-1}\circ R\circ T\phi
\end{equation}

and a discretization scheme

\begin{equation}
\label{disc_sys_2}
R^{-1}_{\phi^{-1}}(x_k,x_{k+1}) = hX(\tau_M(R_{\phi^{-1}}^{-1}(x_k,x_{k+1}),u_k))
\end{equation}

Let $R$ be a discretization map  for the feedback linearized continuous system (<ref>) preserving linearity, then the nonlinear system (<ref>) has a discretization given by (<ref>) that is feedback linearizable in the discrete-time domain.

Define $y_k = \phi(x_k)$ and $\psi(x_k,u_k) = v_k$ for all $k$. From (<ref>) we have

\begin{eqnarray*}
R_{\phi^{-1}}^{-1}(x_k,x_{k+1}) &=& (T\phi^{-1}\circ R^{-1}\circ(\phi\times\phi))(x_k,x_{k+1})\\
&=& (T\phi^{-1}\circ R^{-1})(y_k,y_{k+1}),
\end{eqnarray*}

\begin{eqnarray*}
&=&X\big(\tau_M((T\phi^{-1}\circ R^{-1})(y_k,y_{k+1})),u_k)\big)\\
&=&X_\phi\big(\tau_N(R^{-1}(y_k,y_{k+1})), u_k)\big)
\end{eqnarray*}

From (<ref>) we have $(x_k,x_{k+1}) = R_{\phi^{-1}}(hX(\tau_M(R_{\phi^{-1}}^{-1}(x_k,x_{k+1})))$ and therefore we have

\begin{eqnarray*}
&&(\phi\times\phi)(x_k,x_{k+1})\\&=& (\phi\times\phi)\circ R_{\phi^{-1}}(hX(\tau_M(R_{\phi^{-1}}^{-1}(x_k,x_{k+1})),u_k)\\
\end{eqnarray*}

and therefore we have

\begin{eqnarray*}
(y_k,y_{k+1}) = R(hX_\phi(\tau_N(R^{-1}(y_k,y_{k+1})),u_k)).
\end{eqnarray*}

But since $R$ preserves linearity and  $X_\phi (y,u) = Ay+Bv$
with $v = \psi(\phi^{-1}y,u)$, we have

\begin{eqnarray*}
y_{k+1} &=& A_hy_k +B_h\psi(x_k,u_k)\\
&=& A_hy_k+B_hv_k
\end{eqnarray*}

This concludes the fact that discretization given by (<ref>) is feedback linearizable under the coordinate change $x\mapsto\phi(x) = y$ and a modified control input $(x,u)\mapsto \psi(x,u) = v$.

It is important to note that independent of the order of $R$ one can  ensure an accuracy of $R_{\phi}$ up to the first order. This is due to the fact that while implementing (<ref>) via the sample and hold, the control input $u$ is to be held constant on the interval $[t_k,t_{k+1}[$. This is in general not possible while simultaneously keeping the linearized control input $v$ constant over $[t_k,t_{k+1}[$ as $v(t) =\psi(x(t),u(t))$. Instead of employing the exact control input $u(t)$ over the interval, we apply the control $u_k$ satisfying $v_k = \psi(x_k,u_k)$ for all $t$ where $x_k$ is the state sampled at $t = t_k$.

§.§ Linearizability of Adjoint discretization

Given a discretization map, $R_{\phi^{-1}}$ one can construct an adjoint discretization $R^*_{\phi^{-1}}$ as given by the Definition <ref>. From proposition <ref> we have

\begin{equation}
\label{disc_len_232}
R^*_{\phi^{-1}} = (\phi\times\phi)^{-1}\circ R^*\circ T\phi.
\end{equation}

Using the definition of $R^*$ and the inversion map, $R^*$ induces a following discretization scheme

\begin{equation*}
(R^*)^{-1}(y_{k+1},y_k) = -hX_\phi(\tau_N((R^*)^{-1}(y_{k+1},y_k),v  _k)).
\end{equation*}

Suppose $R$ is such that it preserves the linearity of (<ref>) for $-h$ as well i.e.,  the following discretization

\begin{eqnarray*}
R^{-1}(y_k,y_{k+1}) = -hX_\phi(\tau_N(R(y_k,y_{k+1})),u_k)
\end{eqnarray*}

results in a linear discrete-time system of the form

\begin{equation}
y_{k+1} = A_{-h}y_k + B_{-h}v_k
\end{equation}

then we have the following result.

Let $R$ be a discretization of (<ref>) preserving linearity for $h$ as well as $-h$. Let $R^*$ be the adjoint of $R$, then $R^*_{\phi^{-1}}$ given by (<ref>) results in a discretization

\begin{equation}
x_{k+1} = F^*(x_{k},u_{k};h)
\end{equation}

i.e., feedback is linearizable. Moreover, the linearizing coordinate is given by $x\mapsto\phi(x):=y$ and the linearized system is given by

\begin{equation}
\label{disc_sys_lin2}
y_k = A_{-h}y_{k+1}+B_{-h}v_k
\end{equation}

with $v_k = \psi(x_{k+1},u_k)$.

The proof of the above theorem follows a similar process to that of Theorem <ref> and is hence omitted. Moreover, the control input $u_k
$ is to be calculated implicitly from the control input $v_k$. Similar to $R_{\phi^{-1}}$, $R_{\phi^{-1}}^*$ is also accurate upto first order.

§ CONSTRUCTING HIGHER ORDER DISCRETIZATIONS

Consider the continuous time system (<ref>), then for a given discretization (<ref>) the $n$-step global error is defined as
$$e(k) \coloneqq x(t_k)-x_k,$$
where $x(t_k)$ is the exact solution of (<ref>) evaluated at $t= hk$. The one-step truncated error at $t_k$ is

\begin{equation*}
\tilde{x}_k = \left(F(x(t_k),u_k))-x(t_k)\right)/{h}
\end{equation*}

For an $r$-order discretizations, the one-step truncation error is bounded above by $\norm{\tilde{x}_k}\leq K\norm{h}^{r}$ <cit.>. For first-order methods, the error varies linearly with the stepsize, therefore one requires a smaller stepsize to have better accuracy. For instance, in order to have an accuracy of an order of $10^{-4}$ the stepsize $h$ is to be of the order of $10^{-4}$, whereas for second-order methods a stepsize of an order of $10^{-2}$ shall suffice. The discretizations $R$ in Definition <ref> are in general first order. However, if $R$ is symmetric then it is accurate up to second accurate. This serves as our motivation to construct symmetric discretization.

§.§ Symmetric Discretizations

Let $M$ be an $n$ dimensional manifold and $X\in\mathfrak{X}(M)$ be a vector field on $M$. Let $R$ be a discretization map on $M$ and $R^*$ be its associated adjoint. Define by composition a discretization scheme as follows :

\begin{equation}
\label{symm_int}
\begin{split}
R^{-1}(x_k,x_{k+1/2}) &= \frac{h}{2}X(\tau_M(R^{-1}(x_k,x_{k+1/2})))\\
(R^*)^{-1}(x_{k+1/2},x_{k+1}) &= \frac{h}{2}X(\tau_M((R^*)^{-1}(x_{k+1/2},x_{k+1})))
\end{split}
\end{equation}

In the above equation, $x_{k+1/2}\in M$ is to be taken as an intermediate point and is to be solved implicitly to get a discrete system of type (<ref>). For any $x,y\in M$ the following is true

\begin{equation}
\label{5.2}
(R^*)^{-1}(x,y) =  I_{TM}\circ R^{-1}(y,x)
\end{equation}

\begin{equation}
\label{5.3}
\end{equation}

where $ TM\ni(x,v_x)\mapsto I_{TM}(x,v_x)=(x,-v_x)\in TM$.

Discretization given by (<ref>) is symmetric.

Replacing $(x_k,x_{k+1})$ and $h$ with $(x_{k+1},x_k)$ and $-h$ in (<ref>) we have

\begin{equation*}
\begin{split}
R^{-1}(x_{k+1},x_{k+1/2}) &= -\frac{h}{2}X(\tau_M(R^{-1}(x_{k+1},x_{k+1/2})))\\
(R^*)^{-1}(x_{k+1/2},x_k) &= -\frac{h}{2}X(\tau_M((R^*)^{-1}(x_{k+1/2},x_{k})))
\end{split}
\end{equation*}

Using (<ref>) and (<ref>) we get

\begin{align*}
(R^*)^{-1}(x_{k+1/2},x_{k+1})&= \frac{h}{2}X(\tau_M((R^*)^{-1}(x_{k+1/2},x_{k+1}))),\\
R^{-1}(x_k,x_{k+1/2})&= \frac{h}{2}X(\tau_M(R^{-1}(x_k,x_{k+1/2})))
\end{align*}

i.e., (<ref>), thereby proving (<ref>) is symmetric.

Since (<ref>) is symmetric and therefore is second order.

For nonautonomous systems, the control input $u_k$ is held constant between $t \in [t_k,t_{k+1}[$, (<ref>) is then modified as

\begin{equation}
\label{symm_int_na}
\begin{split}
R^{-1}(x_k,x_{k+1/2}) &= \frac{h}{2}X(\tau_M(R^{-1}(x_k,x_{k+1/2})),u_k)\\
(R^*)^{-1}(x_{k+1/2},x_{k+1})\\ = &\frac{h}{2}X(\tau_M((R^*)^{-1}(x_{k+1/2},x_{k+1})),u_k)
\end{split}
\end{equation}

Under closed-loop performance i.e., applying a feedback control $u_k =u(x_k)$ (<ref>) loses its symmetric nature. This can be overcome by employing multirate sampling methods.

§.§ Multirate Sampling

Consider a continuous time system given by (<ref>). Let $h$ be the sampling time interval i.e., $x_k= x(t_k)$ and $t_{k+1} = t_k+h$. For a fixed $N\in \{1,2,\ldots,n\}$, and for each $i\in\{1,2,\ldots N\}$ let $(x,u)\mapsto F_i(x,u)\eqqcolon F_i^u(x)$ be discretizations of (<ref>). The $N^{th}$ step evolution is then given by

\begin{equation}
x_{k+N} =F_N^{u_{k+N-1}}\circ\ldots\circ F_2^{u_{k+1}}\circ F_1^{u_k}(x_k)
\end{equation}

Sampling states $x_k$ at a rate $N$ times that of control input $u_k$ we get a multistep discretization given by

\begin{equation}
x_{k+N} = \bar F(x_k,u_k,\ldots u_{k+N-1})
\end{equation}

The control input $u_k,\ldots,u_{k+N-1}$ are to be computed a priori at $t_k$ are functions of the state $x_{Nm}$, $Nm<k$.

Setting $N=2$ and $F_1$ and $F_2$ as $F$ and $F^*$ from (<ref>) and (<ref>) respectively. Under multirate sampling, the discrete system generated by (<ref>) is given by

\begin{eqnarray*}
x_{k+1/2} &=& F(x_k,u_k;h/2)\\
x_{k+1/2} &=& F(x_{k+1},u_{k+1};-h/2).
\end{eqnarray*}

Let $u_k = u(x_k)$ be a closed-loop control input for discretization (<ref>). Setting $u_{k+1} = u(x_{k+1})$ renders (<ref>) symmetric and the discretization is given by

\begin{equation}
\label{symm_int_Na1}
F(x_k,u(x_k);h/2) =F(x_{k+1},u(x_{k+1});-h/2)
\end{equation}

which is symmetric and therefore is accurate up to second order. Corresponding continuous time control input is

\begin{eqnarray*}
u(t) = \begin{cases}
&u_k,~ \quad t\in[t_k,t_{k+1/2}[\\
&u_{k+1},~\quad t\in[t_{k+1/2},t_{k+1}[
\end{cases}
\end{eqnarray*}

where $t_{k+1/2} = t_k+\frac{h}{2}$ and $t_{k+1}= t_k+h$.

Consider the continuous time system given by (<ref>). Let $R_\phi$ be its discretization as given by (<ref>) and (<ref>) be its associated discretization. Then one can construct a symmetric discretization given by (<ref>), the resulting discrete system given by

\begin{eqnarray}
\label{symm_disc_sys}
F(x_k,u_k;h/2) = F(x_{k+1},u_{k+1};-h/2)
\end{eqnarray}

is symmetric and is of second order. Moreover, (<ref>) is feedback linearizable under coordinates $x\mapsto \phi(x)\eqqcolon y$ and the modified control input is given by $(x,u)\mapsto \psi(x,u)\eqqcolon v$. The linearized system is given by

\begin{equation}
\label{disc_symm_lin}
A_{h'}y_{k} + B_{h'}v_k =  A_{-h'}y_{k+1} + B_{-h'}v_{k+1}
\end{equation}

where $h' = h/2$.

(<ref>) being symmetric is trivial from the definition itself. Since $\phi$ is a local diffeomorphism we have

\begin{eqnarray*}
\phi(F(x_k,u_k;h')) &=& \phi(F(x_{k+1},u_{k+1});-h')\\
A_{h'}\phi(x_k)+B_{h'}\psi(x_k,u_k) &=& A_{-h'}\phi(x_{k+1})+\\ &~&~B_{-h'}\psi(x_{k+1},u_{k+1})\\
A_{h'}y_{k}+B_{h'}v_k &=&  A_{-h'}y_{k+1} + B_{-h'}v_{k+1},
\end{eqnarray*}

thereby completing the proof.

The control input $v_k$, $v_{k+1}$ can be computed apriori at $t=t_k$ from $z_k$. The control input $u_k$ and $u_{k+1}$ are than computed implicitly solving $\psi(x_k,u_k)= v_k$, with $x_k =\phi^{-1}(z_k)$.

Theorem <ref> is different from the result in <cit.> in the sense that here the rate of multi-sampling is fixed apriori while <cit.> the order of discretization is chosen so that the resulting scheme is feedback linearizable.

§ EXAMPLE

In order to demonstrate the ideas discussed we consider the following example. Consider the following dynamical system evolving on $M=\R[2]$ and $U=\R$.

\begin{equation}
\label{ex1}
\odv{}{t}\pmat{x_1(t)\\x_2(t)} = \pmat{a\sin(x_2(t))\\-(x_1(t))^2+u(t)}
\end{equation}

where $a\in\R\setminus\{0\}$ is fixed and given. Define $\phi(x_1,x_2) =(x_1,a\sin(x_2))\coloneqq (y_1,y_2)$ and $\psi(x_1,x_2,u) =(-(x_1)^2+u)a\cos(x_2) \coloneqq v$. The linearized system is then given by

\begin{equation}
\label{ex_lin}
\odv{}{t}\pmat{y_1(t)\\y_2(t)} = \pmat{y_2(t)\\v(t)}.
\end{equation}

Choosing the Explicit Euler discretization i.e., $\R[2]\times\R[2]\ni(y,w)\mapsto R(y,w) = (y,y+w))\in\R[2]\times\R[2]$, the discrete system is given by

\begin{equation}
\label{exln_dis}
\begin{split}
y_{1,k+1} &= y_{1,k}+hy_{2,k},\\
y_{2,k+1} &= y_{2,k}+hv_k.
\end{split}
\end{equation}

Lifting $R$ via $\phi^{-1}$ to get a discretization for (<ref>), $R_{\phi^{-1}}$ induces the following discretization scheme

\begin{equation}
\tag{EES}
\label{ex_dis_el}
\begin{split}
x_{1,k+1} &= x_{1,k}+ha\sin(x_{2,k})\\
x_{2,k+1} &= \arcsin(\sin(x_{2,k})+h(-x_{1,k}^2+u_k)\cos(x_{2,k})).
\end{split}
\end{equation}

which can be compactly written as $x_{k+1}= F(x_k,u_k;h)$ with $F((x_1,x_2),u;h) = (x_1+ha\sin(x_2),\arcsin(\sin(x_2)+h(-x_1^2+u)\cos(x_2)))$. It is easy to see that (<ref>) is feedback linearizable for around the equilibrium point $(0,0,0)$.

The associated adjoint scheme $R^*(y,w) = (y-w,w)$ defines the Implicit Euler Discretization. Lifting $R^*$, $R_{\phi^{-1}}^*$ induces the following discretization scheme

\begin{equation}
\tag{IES}
\label{IES1}
x_{k} = F(x_{k+1},u_k;-h)
\end{equation}

A symmetric integrator can be defined for Example (<ref>) by composing $R_{\phi^{-1}}$ and $R_{\phi^{-1}}^*$. Using multi-rate sampling with states sampled at a
rate twice that of control. The symmetric integrator for (<ref>) is given by

\begin{equation}
\tag{SES}
\label{ex_dis_sym}
F(x_{k+1},u_{k+1};-h) = F(x_{k},u_k;h).
\end{equation}

One can check that (<ref>) is also feedback linearizable and the linearized system is given by

\begin{equation}
A_{-h'}y_{k+1} +B_{-h'}v_{k+1} = A_{h'}y_{k}+B_{h'}v_k
\end{equation}

with $A_{h'} = \pmat{1&h'\\0&1}$ and $B_{h} = \pmat{0\\h'}$ for $h'= h/2$ and $v_k = \psi(x_k,u_k)$.

The three schemes were simulated under the following parameters: For all three schemes, we had $a=1$, and the initial condition was chosen as $x(0) = (0.25,\pi/6)$. For $K=\bmat{-10&-10}$, the control schemes were chosen as in Table <ref>.

Discretization    Associated Control
(<ref>)     $v_k = Ky_k$
(<ref>)     $v_k = Ky_{k+1}$
(<ref>)     $v_k = Ky_k,~~ v_{k+1} = Ky_{k+1}$

Control input for various discretization schemes. ($K = \bmat{-10&-10}$)

The schemes were simulated for various step sizes and the error was compared with the standard ODE solver (ODE45) available in MATLAB. The simulation was run for $t\in[0,5]$ and the trajectories for the discrete-time system (<ref>) and continuous-time system (<ref>) are plotted in Figure <ref>. The corresponding error is plotted in Figure <ref>. The control input is plotted in <ref>. Similarly, the system trajectory  for (<ref>) is plotted in Figure <ref>, and the control input and error in Figures<ref> and <ref> respectively. To show the multi-rate sampling,  A zoomed-in version of the control signal around $t_k=0$ is also plotted in figure <ref>. It can be seen that the absolute error for (<ref>) is significantly smaller than that of (<ref>). In Figure <ref>, we compare the (percentage) relative error $100*\norm{e(t_k)}/\norm{x(t_k)}$, for (<ref>), (<ref>) and (<ref>). While percentage error increases for (<ref>) and (<ref>), the absolute error is actually quite small (the increase is partly because of the precision errors when $x$ becomes small). For various stepsizes,  Table <ref> records the order of the error magnitude for the various discretization. It can be seen that (<ref>) outperform (<ref>) and (<ref>) and the error is proportional to $h^2$.

3l|Order of error magnitude
1|l|Stepsize     1l|(<ref>)     1l|(<ref>)     (<ref>)
1|l|$h = 10^{-1} $     1l|$10^{-1}$     1l|$10^{-1}$     $10^{-2}$
1|l|$h= 10^{-2}$      1l|$10^{-3}$     1l|$10^{-3}$     $10^{-4}$
1|l|$h=10^{-3}$       1l|$10^{-4}$     1l|$10^{-4}$     $10^{-6}$
1|l|$h=10^{-5}$       1l|$10^{-5}$     1l|$10^{-4}$     $10^{-8}$

Order of error magnitude for various step sizes.

at =(0,-0.25),
xlabel=$t_k$ (s),
ylabel= ,
xmin= 0, xmax=5,
ymin= -0.4, ymax=0.6,
legend pos=north east,
no markers,
minor tick num =5,
minor tick style=draw=none,
minor grid style=thin,color=black!10,
major grid style=thin,color=black!10,
tick align=outside,
axis x line*=none,
axis y line*=none,

[color=blue,smooth,thick] table [x=t,y=X1,col sep=comma]data_ees.csv;
[color=red,smooth,thick] table [x=t,y=X2,col sep=comma]data_ees.csv;
[color=green,smooth,dashed,thick] table [x=t,y=x1,col sep=comma]data_ees.csv;
[color=black,smooth,dashed,thick] table [x=t,y=x2,col sep=comma]data_ees.csv;

System State $x_k$ for (<ref>) plotted against exact discretization (ODE45) $x(t_k)$ for stepsize $h=10^{-2}$ and $t_k\in[0,5]$.

at =(0,-0.25),
xlabel=$t_k$ (s),
ylabel= ,
xmin= 0, xmax=5,
ymin= -10, ymax=2.0,
legend pos=south east,
no markers,
minor tick num =5,
minor tick style=draw=none,
minor grid style=thin,color=black!10,
major grid style=thin,color=black!10,
tick align=outside,
axis x line*=none,
axis y line*=none,

[color=blue,smooth,thick] table [x=t,y=U1,col sep=comma]data_ees.csv;
[color=red,smooth,dashed,thick] table [x=t,y=u1,col sep=comma]data_ees.csv;

Control input $u_k$ for (<ref>) plotted against exact discretization (ODE45) $u(t_k)$ for stepsize $h=10^{-2}$ and $t_k\in[0,5]$.

at =(0,-0.25),
xlabel=$t_k$ (s),
xmin= 0, xmax=5,
ymin= 0, ymax=0.006,
legend pos=south east,
no markers,
minor tick num =5,
minor tick style=draw=none,
minor grid style=thin,color=black!10,
major grid style=thin,color=black!10,
tick align=outside,
axis x line*=none,
axis y line*=none,

[color=blue!80,smooth,thick] table [x=t,y=hn,col sep=comma]data_ees.csv;

Magnitude of global error for (<ref>) for a stepsize $h=10^{-2}$ and $t_k\in[0,5]$.

at =(0,-0.25),
xlabel=$t_k$ (s),
ylabel= ,
xmin= 0, xmax=5,
ymin= -0.4, ymax=0.6,
legend pos=north east,
no markers,
minor tick num =5,
minor tick style=draw=none,
minor grid style=thin,color=black!10,
major grid style=thin,color=black!10,
tick align=outside,
axis x line*=none,
axis y line*=none,

[color=blue,smooth,thick] table [x=t,y=X1,col sep=comma]data_sym.csv;
[color=red,smooth,thick] table [x=t,y=X2,col sep=comma]data_sym.csv;
[color=green,smooth,dashed,thick] table [x=t,y=x1,col sep=comma]data_sym.csv;
[color=black,smooth,dashed,thick] table [x=t,y=x2,col sep=comma]data_sym.csv;

System State $x_k$ and for (<ref>) plotted against exact discretization (ODE45) $x(t_k)$ for stepsize $h=10^{-2}$ and $t_k\in[0,5]$.

at =(0,-0.25),
xlabel=$t_k$ (s),
ylabel= ,
xmin= 0, xmax=5,
ymin= -10, ymax=2.0,
legend pos=south east,
no markers,
minor tick num =5,
minor tick style=draw=none,
minor grid style=thin,color=black!10,
major grid style=thin,color=black!10,
tick align=outside,
axis x line*=none,
axis y line*=none,

[color=blue,smooth,thick] table [x=t,y=U,col sep=comma]data_sym_cont.csv;
[color=red,smooth,dashed,thick] table [x=t,y=u,col sep=comma]data_sym_cont.csv;

Control input $u_k$ for (<ref>) plotted against exact discretization (ODE45) $u(t_k)$ for stepsize $h=10^{-2}$ and $t_k\in[0,5]$.

at =(0,-0.25),
xlabel=$t_k$ (s),
xmin= 0, xmax=5,
ymin= 0, ymax=0.0006,
legend pos=south east,
no markers,
minor tick num =5,
minor tick style=draw=none,
minor grid style=thin,color=black!10,
major grid style=thin,color=black!10,
tick align=outside,
axis x line*=none,
axis y line*=none,

[color=blue!80,smooth,thick] table [x=t,y=hn,col sep=comma]data_sym.csv;

Magnitude of global error for (<ref>) for a stepsize $h=10^{-2}$ and $t_k\in[0,5]$.

at =(0,-0.25),
xlabel=$t_k$ (s),
ylabel= ,
xmin= 0, xmax=0.05,
ymin= -10, ymax=-4,
legend pos=south east,
no markers,
minor tick num =5,
minor tick style=draw=none,
minor grid style=thin,color=black!10,
major grid style=thin,color=black!10,
tick align=outside,
axis x line*=none,
axis y line*=none,

[color=blue,const plot,thick] table [x=t1,y=U1,col sep=comma]data_sym_cont.csv;
[color=red,const plot,dashed,thick] table [x=t1,y=u1,col sep=comma]data_sym_cont.csv;

Zoomed-in Control input $u_k$ for (<ref>)  stepsize $h=10^{-2}$ and $t_k\in[0,0.05]$.

at =(0,-0.25),
xlabel=$t_k$ (s),
xmin= 0, xmax=5,
ymin= 0, ymax=2.25,
legend pos=south east,
no markers,
minor tick num =5,
minor tick style=draw=none,
minor grid style=thin,color=black!10,
major grid style=thin,color=black!10,
tick align=outside,
axis x line*=none,
axis y line*=none,

[color=blue!80,smooth,thick] table [x=t,y=hr,col sep=comma]data_ees.csv;
[color=red!80,smooth,thick] table [x=t,y=hr,col sep=comma]data_ies.csv;
[color=green!80,smooth,thick] table [x=t,y=hr,col sep=comma]data_sym.csv;

Comparing (percentage) relative error for (<ref>), (<ref>) and (<ref>) for $h=10^{-2}$

§ CONCLUSIONS

In this article, we have utilized the idea of retraction maps and their lifts under diffeomorphism to construct feedback linearizable discretization. Given a continuous-time feedback linearizable system, we show that one can build first-order discretization that preserves feedback linearizability. This is done by lifting a discretization of the linearized continuous time system. We have also shown a way to functionally compose two first-order discretizations to design second-order discretizations that are feedback linearizable. However, this comes at the cost of multi-rate sampling.

§ ACKNOWLEDGMENT

We would like to thank Debasish Chatterjee
(Professor, Systems and Control Engineering, IIT Bombay) for his ideas and  discussions which were crucial in the development of this work.