# Accurate estimation of feature importance faithfulness for tree models

Mateusz Gajewski<EMAIL_ADDRESS>IDEAS NCBR
o Faculty of Computing and Telecommunications
Poznan University of Technology
Poznan, Poland  Adam Karczmarz<EMAIL_ADDRESS>Faculty of
Mathematics, Informatics and Mechanics
University of Warsaw
Warsaw, Poland
o IDEAS NCBR
o Mateusz Rapicki<EMAIL_ADDRESS>Faculty of Mathematics,
Informatics and Mechanics
University of Warsaw
Warsaw, Poland
o Piotr Sankowski<EMAIL_ADDRESS>Faculty of Mathematics,
Informatics and Mechanics
University of Warsaw
Warsaw, Poland
o IDEAS NCBR
o

###### Abstract

In this paper, we consider a perturbation-based metric of predictive
faithfulness of feature rankings (or attributions) that we call _PGI squared_.
When applied to decision tree-based regression models, the metric can be
computed accurately and efficiently for arbitrary independent feature
perturbation distributions. In particular, the computation does not involve
Monte Carlo sampling that has been typically used for computing similar
metrics and which is inherently prone to inaccuracies.

Moreover, we propose a method of ranking features by their importance for the
tree model’s predictions based on PGI squared. Our experiments indicate that
in some respects, the method may identify the globally important features
better than the state-of-the-art SHAP explainer.

##  1 Introduction

One of the key challenges in deploying modern machine learning models in such
areas as medical diagnosis lies in the ability to indicate why a certain
prediction has been made. Such an indication may be of critical importance
when a human decides whether the prediction can be relied on. This is one of
the reasons various aspects of explainability of machine learning models have
been the subject of extensive research lately (see, e.g., [3]).

For some basic types of models (e.g., single decision trees), the rationale
behind a prediction is easy to understand by a human. However, predictions of
more complex models (that offer much better accuracy, e.g., based on neural
networks or decision tree ensembles) are also much more difficult to
interpret. Accurate and concise explanations understandable to humans might
not always exist. In such cases, it is still beneficial to have methods
_giving a flavor_ of what factors might have influenced the prediction the
most.

Local feature attribution111It should be contrasted with _global_ feature
attribution whose aim is to measure the feature’s total impact on the model.
constitutes one of such general approaches. For a fixed input $x$ and the
model’s prediction $f(x)$, each feature is assigned a weight “measuring” the
feature’s impact on the prediction. Ideally, a larger (absolute) weight should
correspond to a larger importance of the feature for the prediction.

Many local attribution methods have been proposed so far, e.g., [12, 14, 15].
They are often completely model-agnostic, that is, they can be defined for and
used with any model which is only accessed in a black-box way for computing
the attribution. In some of these cases, focusing on a particular model’s
architecture can lead to much more efficient and accurate attribution
computation algorithms (e.g., [10] for tree ensemble models).

With the large body of different attribution methods, it is not clear which
one should one use, especially since there is not a single objective measure
of their reliability. To help deal with this issue, several notions of
explanation quality have been proposed in the literature, such as faithfulness
(fidelity) [11, 17], stability [2], or fairness [4, 8]. [1] developed an open-
source benchmark OpenXAI automatically computing variants of these measures
for a number of proposed attribution methods.

Out of these notions, our focus in this paper is on faithfulness.

Roughly speaking, an attribution method deserves to be called faithful if the
features deemed important by the attribution are truly important for the
decision-making process of the model. The _perturbation-based_ methods
constitute one popular class of approaches to measuring faithfulness.
Intuitively, perturbing features deemed impactful should generally lead to a
significant change in the prediction. On the contrary, manipulating
unimportant features should not make a big difference.

One concrete cleanly-defined perturbation-based faithfulness metric for
regression problems is the _prediction gap on important feature perturbation_
, PGI in short. PGI is a faithfulness measure of choice in the OpenXAI
benchmark [1]. Having fixed some subset of important features $S$ derived from
the obtained attribution, the _prediction gap_ wrt. $S$ is defined as the
expected absolute change in the prediction when an independent Gaussian noise
is applied to each coordinate $x_{j}$, $j\in S$. PGI can be defined as either
the prediction gap wrt. a fixed set of important features (e.g., $k$ top-
scoring features for some fixed $k$ [8]), or an average prediction gap over
many important feature sets (as in [1]). In particular, the PGI metric depends
_only_ on the ordering of features by importance induced by the attribution.
Other perturbation-based methods include e.g., PGU [1], comprehensiveness,
sufficiency [7], remove-and-retrain [9], deletion [13].

In practice, PGI and other methods using random perturbations are computed
using the Monte Carlo method which is inherently prone to inaccuracies, unless
a very large number of samples is used.

It is interesting to ask whether PGI or similar random perturbation-based
faithfulness matrices can be computed more accurately, say, with the
approximation error caused merely by the floating-point arithmetic’s rounding
errors.

###  1.1 Our contribution

####  1.1.1 Exact computation

First, we consider random perturbation-based measurement of faithfulness from
the point of view of efficient and accurate computation. Of course, it is
unrealistic to assume that such a quantity can be computed exactly (or even
beyond the Monte Carlo method) given only black-box access to the model and
the perturbation distribution. Hence, we focus on a concrete model
architecture: tree ensemble models. Tree ensemble models remain a popular
choice among practitioners since they are robust, easy to tune, and fast to
train.

Unfortunately, computing the _prediction gap_ , as defined in [1] (wrt. a
fixed set $S$ of perturbed features, see Section 2 for a formal definition),
does not look very tractable given the rich structure that tree ensemble
models offer. This is because therein, the expectation is taken over the
_absolute value_ of $f(x)-f(x^{\prime})$, where $x^{\prime}$ is $x$ with
important features perturbed. The absolute value is not very mathematically
convenient to work with when combined with expectations. On the other hand,
dropping the absolute value here would not make much sense: a small random
perturbation (positive or negative) of a feature could in principle lead to a
large positive or negative prediction difference. These, in turn, could cancel
out in expectation, making the metric close to zero even though the feature is
visibly of high importance.

To deal with the tractability problem of measuring PGI on tree models, we
propose a similar metric, $\text{PGI}^{2}$, where instead of taking the
absolute value, we compute the expectation of the square prediction gap
$(f(x)-f(x^{\prime}))^{2}$. This eliminates the cancellation problem outlined
above. Of course, introducing the square can alter the evaluation of the
feature’s importance wrt. PGI as defined previously. However, intuitively, it
amplifies the big prediction differences, and at the same time marginalizes
the small, which is a generally good property if we seek concise explanations.

We prove that for tree ensemble models with $n$ nodes in total, the squared
prediction gap can be computed in $O(n^{2})$ time exactly for any fixed choice
of subset of perturbed important features, assuming the cumulative
distribution function of feature-wise random perturbations can be evaluated in
constant time. That is, the obtained running time bound does not require that
the random perturbations are Gaussian or their distributions are equal across
features.

The OpenXAI benchmark calculates PGI by taking the average expected prediction
gap obtained when perturbing $1,2,\ldots,d$ most important features (where $d$
denotes the total number of features; see [1, Appendix A]). In such a case,
$d$ different expectations of absolute prediction gap are estimated, each of
them via Monte Carlo sampling. In our case, $\text{PGI}^{2}$ is defined as the
average expected squared prediction gap when perturbing $1,2,\ldots,d$ most
important features. Using our prediction gap algorithm, $\text{PGI}^{2}$ is
computed exactly (assuming infinite precision arithmetic; in reality, up to
the precision error incurred by the usage of floating-point arithmetic) in
$O(n^{2}d)$ time.

##### Experiments.

We accompany our theoretical development with an experimental comparison of
the results obtained using the exact algorithm and Monte Carlo sampling.

The experiments confirm that Monte Carlo simulation gives noticeably less
accurate results within the allocated time (iteration) budget comparable to
the running time of our algorithm. Specifically, in our experiments, the
normalized mean absolute error (NMAE) is typically between $0.02$ and $0.3$.
The error measured this way generally decreases when increasing the size of
the model which follows from the quadratic complexity of our algorithm. It is
also decreasing in the magnitude of perturbations used: sampling is generally
better at estimating larger quantities, whereas tiny perturbations lead to
tiny prediction gaps.

####  1.1.2 $\text{PG}^{2}$-based greedy feature ordering.

We also investigate the possibility of using our exact $\text{PG}^{2}$
algorithm as a base of a feature importance ranking algorithm. We stress that
finding a feature ordering optimizing the $\text{PGI}^{2}$ metric is a highly
non-trivial task, as the metric depends on the entire ordering of the
features. Checking all the features’ permutations is infeasible in most cases.

We consider constructing the feature ranking using a _greedy_ $\text{PG}^{2}$
heuristic: the $i$-th most important feature $f$ is chosen so that it
optimizes the prediction gap together with the already chosen $i-1$ most
important features plus $f$.

We next compare the greedy $\text{PG}^{2}$ ranking with the ranking produced
using the popular SHAP feature attribution method for tree ensembles [10].
First of all, our experiments confirm that the greedy $\text{PG}^{2}$ ranking
yields, on average, better $\text{PGI}^{2}$ scores than SHAP.

While that such a property holds is non-obvious, it is perhaps not very
surprising either. This is why we also compare our ranking method with SHAP in
terms of average _conciseness_. Specifically, we consider three natural
metrics of averaged feature importance across the test data set, each of them
based solely on the feature ordering and clearly promoting the most important
features. For each of these metrics, we compare the entropy of the obtained
averaged importances. Intuitively, a smaller entropy indicates that only few
features are identified as critically important to the model.

Our experiments show that the greedy $\text{PG}^{2}$ ranking used with
perturbations of sufficiently large standard deviation can lead to global
importance scores with smaller entropy compared to SHAP.

##  2 Preliminaries

Let us denote by $f:\mathbb{R}^{d}\to\mathbb{R}$ the output function of the
considered regression model. We sometimes also use $f$ to refer to the model
itself. The input $x\in\mathbb{R}^{d}$ to $f$ is called a feature vector, and
we denote by $x_{i}$ the value of the $i$-th feature of $x$. Thus, we identify
the set of features with $[d]$.

##### Tree models.

When talking about decision trees, we assume them to be binary and based on
single-value splits. That is, each non-leaf node $v$ of a decision tree
$\mathcal{T}$ has precisely two children $a_{v},b_{v}$. It is also assigned a
feature $q_{v}\in[d]$ and a threshold value $t_{v}\in\mathbb{R}$. Each leaf
node $l\in\mathcal{T}$ is in turn assigned a value $y_{v}\in\mathbb{R}$. We
denote by $\mathcal{L}(\mathcal{T})$ the set of leaves of the tree
$\mathcal{T}$.

The output $f_{\mathcal{T}}(x)$ of the tree $\mathcal{T}$ is computed by
following a root-leaf path in $\mathcal{T}$: at a non-leaf node
$v\in\mathcal{T}$, we descend either to the child $a_{v}$ if
$x_{q_{v}}<t_{v}$, or to $b_{v}$ otherwise. When a leaf $l$ is eventually
reached, its value $y_{l}$ is returned. We denote by $n_{\mathcal{T}}$ the
number of nodes in $\mathcal{T}$.

When considering tree ensemble models $(\mathcal{T})_{i=1}^{m}$, the output
$f(x)$ of the model is simply the sum of outputs $f_{\mathcal{T}_{i}}(x)$ of
its $m$ individual trees. We generally use
$n=\sum_{i=1}^{m}n_{\mathcal{T}_{i}}$ to refer to the _total size_ of the tree
ensemble model.

##### Prediction gap on important feature perturbation.

Let a _ranking_ $\pi\in[d]\to[d]$ be a permutation of features giving their
importance from highest to lowest. Denote by $\pi[1..\ell]$ the set
$\\{\pi(1),\ldots,\pi(\ell)\\}$, that is, the $\ell$ top ranking features in
$\pi$. For any $x\in\mathbb{R}^{d}$ and $S\subseteq[d]$, let
$\text{perturb}(x,S)$ be some distribution on feature vectors
$x^{\prime}\in\mathbb{R}^{d}$ satisfying $x^{\prime}_{i}=x_{i}$ for all
$i\notin S$. AKS+ [1] define:

$\mathrm{PGI}(x,\pi,k):=\mathbb{E}_{x^{\prime}\sim\text{perturb}(x,\pi[1..k])}\left[|f(x^{\prime})-f(x)|\right].$

Based on that, they define an averaged metric:

$\mathrm{PGI}(x,\pi):=\frac{1}{d}\sum_{k=1}^{d}\mathrm{PGI}(x,\pi,k).$

We analogously define the $\mathrm{PGI}^{2}$ metric:

$\displaystyle\mathrm{PGI}^{2}(x,\pi,k)$
$\displaystyle:=\mathbb{E}_{x^{\prime}\sim\text{perturb}(x,\pi[1..k])}\left[\left(f(x^{\prime})-f(x)\right)^{2}\right],$
$\displaystyle\mathrm{PGI}^{2}(x,\pi)$
$\displaystyle:=\frac{1}{d}\sum_{k=1}^{d}\mathrm{PGI}^{2}(x,\pi,k).$

##  3 Computing $\mathrm{PGI}^{2}$ for tree ensemble models

Consider distributions $\text{perturb}(x,S)$ such that
${x^{\prime}\sim\text{perturb}(x,S)}$ satisfies:

  1. 1.

for every $i\in S$, $x^{\prime}_{i}=x_{i}+\delta_{i}$, where
$\delta_{i}\sim\mathcal{D}_{i}$, and

  2. 2.

all $\delta_{j}$ ($j\in S$) are independent random variables.

Moreover, we assume that for all $i\in S$, the c.d.f. $F_{\mathcal{D}_{i}}$ of
$\mathcal{D}_{i}$ can be evaluated in $O(1)$ time.

In the following, we show how to compute

$\mathrm{PG}^{2}(x,S):=\mathbb{E}_{x^{\prime}\sim\text{perturb}(x,S)}\left[(f(x^{\prime})-f(x))^{2}\right]$
(1)

under the assumptions made in the case when $f$ is a tree ensemble model
$(\mathcal{T})_{i=1}^{m}$. Note that the metric $\mathrm{PGI}^{2}(x,\pi)$ for
a given ranking $\pi$ can be computed by running the algorithm $d$ times for
different subsets $S$ of the form $\pi[1..k]$.

Put $c=f(x)$. Let $x^{\prime}\sim\text{perturb}(x,S)$. Recall that
${f(x^{\prime})=\sum_{i=1}^{m}f_{\mathcal{T}_{i}}(x^{\prime})}$. For every
leaf node $v$ of some $\mathcal{T}_{i}$, consider a random indicator variable
$X_{v}$ such that $X_{v}=1$ iff $\mathcal{T}_{i}$ evaluates to the value
$y_{v}$ of leaf $v$ given $x^{\prime}$ as input. In particular, for any
$i=1,\ldots,m$, there is exactly one $v\in\mathcal{L}(\mathcal{T}_{i})$ such
that $X_{v}=1$. Set
$\mathcal{L}=\bigcup_{i=1}^{m}\mathcal{L}(\mathcal{T}_{i})$. Then:

$f(x^{\prime})=\sum_{v\in\mathcal{L}}X_{v}\cdot y_{v}.$

By linearity of expectation, we obtain:

$\displaystyle\mathrm{PG}^{2}(x,S)$
$\displaystyle=\mathbb{E}\left[\left(\sum_{v\in\mathcal{L}}X_{v}y_{v}-c\right)^{2}\right]$
$\displaystyle=c^{2}+\mathbb{E}\left[\left(\sum_{u\in\mathcal{L}}X_{u}y_{u}\right)\left(\sum_{v\in\mathcal{L}}X_{v}y_{v}-2c\right)\right]$
$\displaystyle=c^{2}+\sum_{\begin{subarray}{c}u,v\in\mathcal{L}\\\ u\neq
v\end{subarray}}y_{u}y_{v}\Pr[X_{u}=1\land X_{v}=1]$
$\displaystyle\qquad\qquad+\sum_{u\in\mathcal{L}}\Pr[X_{u}=1]\cdot
y_{u}\cdot(y_{u}-2c).$

The above formula reduces computing $\mathrm{PG}^{2}(x,S)$ to calculating
probabilities of the form $\Pr[X_{u}=1\land X_{v}=1]$ (where $u\neq v$) or of
the form $\Pr[X_{u}=1]$. We focus on the former task as dealing with the
latter is analogous but easier.

For a node $w$, let $Q_{w}$ be the set of features appearing on the root-$w$
path (in the tree whose $w$ is the node of). Observe that for
$w\in\mathcal{L}$, $X_{w}=1$ holds iff for each feature $q\in Q_{w}$,
$x^{\prime}_{q}$ falls into a certain interval $I_{w,q}$. Namely, if the root-
to-$w$ path consists of nodes $w_{1},\ldots,w_{k}=w$, then $I_{w,q}$ can be
constructed as follows. Start with the interval $(-\infty,\infty)$, and follow
the path downwards. For each encountered node $w_{i}$ splitting on $q$,
intersect the current $I_{w,q}$ with either:

  * •

$(-\infty,t_{w_{i}})$ if $w_{i+1}=a_{w_{i}}$, or

  * •

$[t_{w_{i}},\infty)$ if $w_{i+1}=b_{w_{i}}$, or

  * •

$\emptyset$ otherwise.

It follows that $X_{u}=1\land X_{v}=1$ iff for all $q\in Q_{u}\cup Q_{v}$,
$x^{\prime}_{q}\in I_{u,q}\cap I_{v,q}$. Hence, if we denote by
$l_{u,v,q},r_{u,v,q}$ the respective endpoints of $I_{u,q}\cap I_{v,q}$, by
the independence of all perturbations $\delta_{j}$, we obtain that
$\Pr[X_{u}=1\land X_{v}=1]$ equals $\Pi(u,v)$, where

$\Pi(u,v)=\left(\prod_{q\in(Q_{u}\cup Q_{v})\setminus S}[x_{q}\in I_{u,q}\cap
I_{v,q}]\right)\cdot\left(\prod_{q\in(Q_{u}\cup Q_{v})\cap
S}\left(F_{\mathcal{D}_{q}}(r_{u,v,q}-x_{q})-F_{\mathcal{D}_{q}}(l_{u,v,q}-x_{q})\right)\right).$
(2)

Implemented naively, computing a single value $\Pi(u,v)$ takes
$O(D_{u}+D_{v})$ time, where $D_{u},D_{v}$ are the depths of $u$ and $v$ in
their respective trees. Indeed, note that by traversing the root-to-$u$ and
root-to-$v$ paths, we can find all the relevant intervals $I_{u,q},I_{v,q}$
that the product $\Pi(u,v)$ depends on. Moreover, $\Pi(u,v)$ has at most
$D_{u}+D_{v}$ factors, each of which can be evaluated in $O(1)$ time.

Computing the values $\Pi(u,v)$ through all $\Theta(n^{2})$ pairs
$u,v\in\mathcal{L}$ could take $O(n^{2}D)$ time, where $D$ is a global bound
on the depth of the trees in the ensemble.

However, the computation can be optimized to run in $O(n^{2})$ time by
computing the probabilities in an adequate order. Namely, note that, say, the
value $\Pi(a_{u},v)$ can be computed based on $\Pi(u,v)$ in $O(1)$ time. This
is because the interval $I_{a_{u},q}$ can differ from the corresponding
interval $I_{u,q}$ only for the feature $q=q_{u}$, and it can be obtained by
splitting $I_{u,q}$ with the threshold $t_{u}$. Hence, only one factor of the
product (2) has to be added or replaced in order to obtain $\Pi(a_{u},v)$ from
$\Pi(u,v)$. An analogous trick applies to computing $\Pi(b_{u},v)$,
$\Pi(u,a_{v})$, $\Pi(u,b_{v})$ from $\Pi(u,v)$. Consequently, by iterating
through nodes $u$ in pre-order, and then, once ${u\in\mathcal{L}}$ is fixed,
through nodes $v$ in pre-order (the order of processing trees does not
matter), all the required values $\Pi(u,v)$ can be computed in $O(n^{2})$
total time.

######  Theorem 3.1.

Let $f$ be a tree ensemble model whose trees have $n$ nodes in total. Let
$x\in\mathbb{R}^{d}$. Then for any $S\subseteq[d]$, $\mathrm{PG}^{2}(x,S)$ can
be computed in $O(n^{2})$ time. In particular, given a feature ranking $\pi$,
for any $k\in[d]$, $\mathrm{PGI}^{2}(x,\pi,k)$ can be computed in $O(n^{2})$
time, whereas $\mathrm{PGI}^{2}(x,\pi)$ can be computed in $O(n^{2}d)$ time.

##  4 Experiments

###  4.1 Datasets

Following the previous work (e.g., [1, 8]), in the experiments we assume that
the (real) perturbations used for computing the prediction gap (either
exactly, or approximately) come from the same normal distribution
$\mathcal{D}=\mathcal{N}(0,\sigma^{2})$. This is why we need to assume that
all the features have numerical values.

Using the same perturbation distribution for all the features requires the
features be standardised by subtracting the mean and dividing by the standard
deviation. This ensures that if a feature is perturbed by a noise variable
drawn from $\mathcal{D}$, the perturbation will have a similar effect. More
specifically, it prevents a perturbation from $\mathcal{D}$ being relatively
significant for one feature and insignificant for another.

In our experiments we used two datasets:

  1. 1.

Red Wine Quality Dataset [5] The dataset contains 11 features related to
physico-chemical attributes of wine, all numerical and continuous. The task is
to predict the score of a wine, which is an integer between 1 and 10. The task
can be considered as both a regression and a classification task. The dataset
contains $1\ 599$ examples.

  2. 2.

California Housing [16] The dataset contains information from the 1990
Californian census on some characteristics of a given Californian county and
the median house value in that county. There are 8 numerical characteristics,
such as length or number of bedrooms, and one categorical characteristic -
proximity to the ocean. For the reasons outlined before, we decided to drop
this feature and use a modified dataset. The task is to predict the median
value of the house. The dataset contains $20\ 640$ examples.

###  4.2 Models

All tree ensembles are implemented in the XGBoost library [6]. For each of the
two datasets, three models of the following types were trained:

  * •

Single tree – a single decision tree with `max_depth` set to 4.

  * •

Bigger – gradient boosted trees with `max_depth` set to 4 and number of trees
limited to 40.

In each case, the data set was split 80:20 into training and test sets. In
addition, for each model and its constraint, a grid search was performed over
selected hyperparameters. The details of the training process can be found in
Appendix B.

###  4.3 Comparing squared prediction gap algorithms

In this section we compare our algorithm computing the squared prediction gap
(as defined in Equation (1)) for tree ensembles described in Section 3 and the
Monte Carlo sampling-based method. Specifically, the latter algorithm, given
$x$ and the set $S\subseteq[d]$ or perturbed features, evaluates
$f(x^{\prime})$ for $i$ randomly sampled $x^{\prime}\sim\text{perturb}(x,S)$
(see Section 3) and records the average value of $(f(x^{\prime})-f(x))^{2}$.
We will call $i$ the number of _iterations_ in the following.

The subject of the comparison was the accuracy of the results produced by the
sampling algorithm. Since the accuracy of the sampling algorithm increases
with the number of iterations $i$, to make the comparison fair, the number $i$
was chosen so that the wall time of both algorithms was close. Both algorithms
were implemented in Python and all the arithmetic operations in the algorithms
were carried out using the `numpy.float32` type.

The computations were carried out on a FormatServer THOR E221 (Supermicro)
server equipped with two AMD EPYC 7702 64-Core processors and 512 GB of RAM
with operation system Ubuntu 22.04.1 LTS.

The comparison between the algorithms was carried out in different settings.
The variable parameters were:

  * •

Model type $m$ trained on a selected dataset. In our case, there were four
models in total (Section 4.2) trained on two datasets (Section 4.1).

  * •

The standard deviation $\sigma$ of a Gaussian used to perturb a feature. In
our experiment, the evaluated values were $\\{0.01,0,03,0.1,0.3\\}$.

To carry out the experiment having fixed the above parameters, we used the
following procedure. For each number of iterations
${i\in\\{100,500,1000,1500,2000,3000,4000,5000,6000\\}}$, we recorded the
average differences in results and computation times over $N=20\ 000$ samples.
These samples varied in the size of the perturbed feature set $S$ (which could
be, e.g., $1,\ldots,11$ for the Wine Quality dataset), and we ensured that all
subset sizes were equally represented in the samples. Each sample involved a
randomly selected point from a test dataset and a set of features of the
desired size. Both algorithms were then run with appropriate parameters on
each prepared sample. We compared the running times and results of the two
algorithms. The error was recorded as the Normalised Mean Absolute Error
(NMAE), defined as:

$\text{NMAE}(y,\hat{y})=\frac{\sum_{j=1}^{N}\mid
y_{j}-\hat{y_{j}}\mid}{\sum_{j=1}^{N}{\mid y_{j}\mid}}$

where $y_{j}$ is the correct value of the prediction gap calculated by the
exact algorithm, $\hat{y_{j}}$ is the prediction gap calculated by the
sampling algorithm. NMAE was used as opposed to the relative error to suppress
the effect of relatively large floating-point errors present when operating on
very small numbers.

The results obtained in the experiment – the error as a function of the number
of iterations – are depicted using the “Mean” curve for $\sigma=0.3$ and
Bigger model for the Red Wine Quality dataset in Figure 1 and for
$\sigma=0.03$ and Bigger model for the Californian Housing dataset in Figure
2. Analogous figures for other settings $(m,\sigma)$ can be found in Appendix
B.

Additionally, the figures contain separate error curves for some selected
sizes $|S|$ of the perturbed features set. As can be observed for a given
iteration, the NMAE error behaves similarly for all subset sizes. Overall, the
difference between the exact and sampling algorithms decreases as the number
of iterations in the sampling algorithm increases. This is indeed anticipated
since Monte Carlo sampling yields better accuracy with iteration increase; in
fact the relative error decreases proportionally to $1/\sqrt{i}$.

Figure 1: Results of precision comparison described in 4.3 for $\sigma=0.3$
and for Bigger model for Red Wine Quality dataset. Figure 2: Results of
precision comparison described in 4.3 for $\sigma=0.03$ and for Bigger model
for California housing dataset.

For the fair comparison, we chose the number of iterations (out of $9$
considered as explained above) in the sampling algorithm that came closest to
the running time of the exact algorithm for each setting. The chosen number of
iterations and the associated NMAE are shown in Tables 1, 2, 3 4.

Table 1: Number of iterations needed to perform to match execution time of exact algorithm on single tree model for Wine Quality Dataset described in 4.3. $\sigma$ | Number of iterations | Relative error
---|---|---
0.01 | 100 | 0.08
0.03 | 100 | 0.14
0.1 | 100 | 0.15
0.3 | 1000 | 0.05
Table 2: Number of iterations needed to perform to match execution time of exact algorithm on Bigger model for Wine Quality Dataset described in 4.3. $\sigma$ | Number of iterations | Relative error
---|---|---
0.01 | 100 | 0.08
0.03 | 100 | 0.11
0.1 | 100 | 0.12
0.3 | 1000 | 0.04
Table 3: Number of iterations needed to perform to match execution time of exact algorithm on single tree model for Californian Housing dataset described in 4.3. $\sigma$ | Number of iterations | Relative error
---|---|---
0.01 | 100 | 0.13
0.03 | 100 | 0.15
0.1 | 100 | 0.15
0.3 | 1000 | 0.04
Table 4: Number of iterations needed to perform to match execution time of exact algorithm on Bigger tree model for Californian Housing dataset described in 4.3. $\sigma$ | Number of iterations | Relative error
---|---|---
0.01 | 500 | 0.3
0.03 | 500 | 0.14
0.1 | 1000 | 0.06
0.3 | 2000 | 0.03

We conclude that the “fair” comparison confirms that the Monte Carlo
simulation gives noticeably less accurate results than our exact method. We
observe that the error generally decreases when increasing the size of the
model. This is justified by the fact that our algorithm’s running time
increases quadratically with $n$, whereas a single iteration of the Monte
Carlo method runs in time at most linear in $n$. As a result, the bigger the
model size is, the more iterations the Monte Carlo method can perform within
the same time budget. Moreover, the error decreases with the standard
deviation of the distribution used to obtain perturbations. This is also
expected: for smaller perturbations, the floating-point rounding errors are of
higher significance.

Table 5: Number of iterations needed to perform to match execution time of exact algorithm on Bigger Californian Housing model described in 4.3. $\sigma$ | Number of iterations | Relative error
---|---|---
0.01 | 500 | 0.30
0.03 | 500 | 0.14
0.1 | 1000 | 0.06
0.3 | 2000 | 0.03

##  5 Feature ranking using $\mathrm{PG}^{2}$

We propose a method for ranking the features of a model using the squared
prediction gap called _greedy_ $\text{PG}^{2}$. Let $x\in\mathbb{R}^{d}$.
Suppose the feature perturbations $\delta_{i}$ all come from
$\mathcal{D}=\mathcal{N}(0,\sigma)$. The ranking
$\pi_{\mathrm{PG}}^{\sigma}(x)$ is constructed in the following inductive
manner.

  1. 1.

For every $i\in[d]$, compute $p_{i}=\mathrm{PG}^{2}(x,\\{i\\})$. Set the most
important feature $\pi^{\sigma}_{\mathrm{PG}}(x)(1)$ to be the feature $i$ for
which $p_{i}$ is the highest.

  2. 2.

Assume $\pi^{\sigma}_{\mathrm{PG}}(x)(1),...,\pi^{\sigma}_{\mathrm{PG}}(x)(l)$
are already computed. Let
${S=\\{\pi^{\sigma}_{\mathrm{PG}}(x)(1),...,\pi^{\sigma}_{\mathrm{PG}}(x)(l)\\}}$.
For every $i\in[d]\setminus S$, compute
$p_{i}:=\mathrm{PG}^{2}(x,S\cup\\{i\\})$. Let the next important feature
$\pi^{\sigma}_{\mathrm{PG}}(x)(l+1)$ be the feature $i$ for which $p_{i}$ is
the highest.

In the following, we compare the feature rankings obtained using the above
method with the rankings obtained using the popular SHAP attribution framework
[10].

For a given $x\in\mathbb{R}^{d}$, SHAP produces an attribution vector
$(\phi)_{i=1}^{d}$ whose elements sum up to $f(x)-\mathbb{E}[f(x)]$ but can be
both positive and negative [12]. We define the ranking $\pi_{\text{SHAP}}(x)$
to be obtained by sorting the features $[d]$ by the _absolute value_ of
$\phi_{i}$, from highest to lowest. This is justified by the fact that the
predictions $f(x)$ can be smaller than the expected value; in such a case, all
the attributions $\phi_{i}$ might be negative.

###  5.1 Global comparison vs SHAP wrt. $\mathrm{PGI}^{2}$

We perform a _global_ comparison of $\pi^{\sigma}_{\mathrm{PG}}$ and
$\pi_{\text{SHAP}}$ wrt. the $\mathrm{PGI}^{2}$ metric, that is, we check
which ranking method gives higher $\mathrm{PGI}^{2}$ scores on average. More
formally, fix a dataset/model combination $m$ and a perturbation distribution
$\mathcal{D}^{\prime}=\mathcal{N}(0,\sigma^{\prime})$, where potentially
$\sigma^{\prime}\neq\sigma$, as we did in Section 4.3. Let
$X\subseteq{\mathbb{R}}^{d}$ be the corresponding dataset used in $m$. Then,
define the average $\mathrm{PGI}^{2}$ scores as

$\displaystyle\overline{{\mathrm{PGI}}^{2}}(\pi_{\text{SHAP}})$
$\displaystyle=\frac{1}{|X|}\sum_{x\in
X}\mathrm{PGI}^{2}(x,\pi_{\text{SHAP}}(x)),$
$\displaystyle\overline{{\mathrm{PGI}}^{2}}(\pi_{\mathrm{PGI}}^{\sigma})$
$\displaystyle=\frac{1}{|X|}\sum_{x\in
X}\mathrm{PGI}^{2}(x,\pi_{\mathrm{PGI}}^{\sigma}(x)).$

Figures 3, 4, 5 and 6 present the comparison of the average scores
$\overline{{\mathrm{PGI}}^{2}}(\pi_{\text{SHAP}})$ and
$\overline{{\mathrm{PGI}}^{2}}(\pi_{\mathrm{PGI}}^{\sigma})$ for $\sigma=0.3$
as a function of $\sigma^{\prime}$. We observe that even if the parameter
$\sigma$ used when computing the ranking $\pi_{\mathrm{PG}}^{\sigma}$ differs
from the parameter $\sigma^{\prime}$ used for computing the $\mathrm{PGI}^{2}$
metrics, $\overline{{\mathrm{PGI}}^{2}}(\pi_{\mathrm{PGI}}^{\sigma})$ is
consistently better than $\overline{{\mathrm{PGI}}^{2}}(\pi_{\text{SHAP}})$
for the Bigger models. For Single tree models, this is the case only for
sufficiently large $\sigma^{\prime}$, but in particular for
$\sigma^{\prime}=\sigma$. The figures presenting the results for other choices
of $\sigma$ can be found in Appendix B.2.

Figure 3: Single tree model for the Wine Quality dataset. Figure 4: Bigger
model for the Wine Quality dataset. Figure 5: Single tree model for the
Housing dataset. Figure 6: Bigger model for the Housing dataset

We conclude that the greedy optimization of the $\mathrm{PGI}^{2}(x,\pi)$
metric performs reasonably well, even though its correctness could only be
guaranteed if it optimized the $\mathrm{PGI}^{2}(x,\pi,1)$ metric (see Section
2) instead.

###  5.2 Conciseness comparison vs SHAP

We also compare the global conciseness of the feature rankings
$\pi_{\text{SHAP}}$ and $\pi_{\mathrm{PG}}^{\sigma}$. Intuitively, we would
like to declare a feature ranking method more concise if it tends to produce
rankings that identify the most important features of the model more clearly.
For example, the information that two or three features of the model are
identified as much more important than the others seems to say more about the
model’s behavior in general than saying that all the features are equally
important on average.

Our conciseness evaluation is performed using four different methods of
converting a feature ranking $\pi(x)$ to a score vector
$r^{\pi}(x):[d]\to\mathbb{R}_{\geq 0}$: _geom_ , _top1_ , _top2_ , _top3_ :

  * •

_geom_ sets $(r^{\pi}(x))_{i}=1/2^{k}$ if $\pi(k)=i$, that is, the scores
increase geometrically with the importance.

  * •

_top $k$_, for $k\in\\{1,2,3\\}$, sets $(r^{\pi}(x))_{i}=1$ iff
${i\in\\{\pi(1),\ldots,\pi(k)\\}}$. In other words, the top $k$ features get a
unit score, and all others get a zero score.

Having fixed (as before) the model-dataset combination $m$, $\sigma$, and a
scoring method, we compute the aggregate feature scores $R_{i}$ through all
data points $x$ from the test dataset $X$:

$R_{i}=\sum_{x\in X}(r^{\pi}(x))_{i}.$

Finally, the conciseness $C_{\pi}$ is computed as the entropy of the
normalized vector $R_{1},\ldots,R_{d}$:

$C_{\pi}=\sum_{i=1}^{d}\frac{R_{i}}{\sum_{j}R_{j}}\log\left(\frac{R_{i}}{\sum_{j}R_{j}}\right).$

Intuitively, higher entropy means more uncertainty in the results.

Tables 6,7, 8, and 9 present the results of our comparison for
$\sigma\in\\{0.1,0.3\\}$ for Wine Quality models, and additionally for
$\sigma=1$ for the Housing models. We observe that for the “single tree”
models, the entropy that SHAP ranking yields is slightly smaller, for most of
the scoring methods. However, there is not much uncertainty in the produced
rankings since the models base their predictions only on a few features.
Indeed, in the case of “single tree” models for both Wine Quality and
Californian Housing datasets, only three features were present in trees.

For the bigger models consisting of numerous trees, all the features are used.
In this more complex case, where the top features are not identified
immediately, the greedy $\mathrm{PG}$-based ranking yields more concise
rankings for the _geom_ , _top2_ and _top3_ metrics for the Wine Quality
dataset. For the Housing dataset, however, the obtained entropies for
$\pi_{\text{SHAP}}$ are also consistently smaller than for the
$\pi_{\mathrm{PG}}^{\sigma}$ rankings for $\sigma\in\\{0.1,0.3\\}$. Only for
$\sigma=1$, the greedy $\mathrm{PG}^{2}$ method looks generally more concise.

Table 6: Entropy comparison between SHAP rankings and PG rankings for Bigger model for Wine Quality Dataset. score | SHAP |  PG $\sigma=0.1$ |  PG $\sigma=0.3$ |  PG $\sigma=1.0$
---|---|---|---|---
_geom_ | 3.245 | 3.285 | 3.157 | 2.740
_top1_ | 2.643 | 3.059 | 2.792 | 1.730
_top2_ | 3.120 | 3.232 | 3.019 | 2.402
_top3_ | 3.305 | 3.300 | 3.189 | 2.846
Table 7: Entropy comparison between SHAP rankings and PG rankings for Bigger model for Californian Housing Dataset. score | SHAP |  PG $\sigma=0.1$ |  PG $\sigma=0.3$ |  PG $\sigma=1.0$
---|---|---|---|---
_geom_ | 2.570 | 2.774 | 2.624 | 2.383
_top1_ | 1.748 | 2.324 | 1.979 | 0.534
_top2_ | 2.394 | 2.651 | 2.377 | 2.176
_top3_ | 2.634 | 2.834 | 2.676 | 2.622
Table 8: Entropy comparison between SHAP rankings and PG rankings for Single tree model for Wine Quality Dataset. score | SHAP | PG, $\sigma=0.1$ | PG, $\sigma=0.3$ | PG, $\sigma=1.0$
---|---|---|---|---
_geom_ | 2.644 | 2.463 | 2.328 | 2.316
_top1_ | 1.725 | 1.482 | 1.372 | 1.305
_top2_ | 2.112 | 1.988 | 1.654 | 1.585
_top3_ | 2.432 | 2.207 | 2.073 | 2.012
Table 9: Entropy comparison between SHAP rankings and PG rankings for Single tree model for Californian Housing Dataset. score | SHAP |  PG $\sigma=0.1$ |  PG $\sigma=0.3$ |  PG $\sigma=1.0$
---|---|---|---|---
_geom_ | 2.094 | 2.197 | 2.175 | 2.138
_top1_ | 1.077 | 0.884 | 0.722 | 0.031
_top2_ | 1.549 | 1.580 | 1.602 | 1.633
_top3_ | 1.585 | 1.998 | 2.020 | 2.066

## References

  * AKS+ [22] Chirag Agarwal, Satyapriya Krishna, Eshika Saxena, Martin Pawelczyk, Nari Johnson, Isha Puri, Marinka Zitnik, and Himabindu Lakkaraju.  Openxai: Towards a transparent evaluation of model explanations.  Advances in Neural Information Processing Systems, 35:15784–15799, 2022.
  * AMJ [18] David Alvarez-Melis and Tommi S Jaakkola.  On the robustness of interpretability methods.  arXiv preprint arXiv:1806.08049, 2018.
  * BH [21] Nadia Burkart and Marco F. Huber.  A survey on the explainability of supervised machine learning.  J. Artif. Intell. Res., 70:245–317, 2021.
  * BZH+ [22] Aparna Balagopalan, Haoran Zhang, Kimia Hamidieh, Thomas Hartvigsen, Frank Rudzicz, and Marzyeh Ghassemi.  The road to explainability is paved with bias: Measuring the fairness of explanations.  In Proceedings of the 2022 ACM Conference on Fairness, Accountability, and Transparency, pages 1194–1206, 2022.
  * CCA+ [09] Paulo Cortez, A. Cerdeira, F. Almeida, T. Matos, and J. Reis.  Wine Quality.  UCI Machine Learning Repository, 2009.  DOI: https://doi.org/10.24432/C56S3T.
  * CG [16] Tianqi Chen and Carlos Guestrin.  Xgboost: A scalable tree boosting system.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’16. ACM, August 2016.
  * DJR+ [20] Jay DeYoung, Sarthak Jain, Nazneen Fatema Rajani, Eric Lehman, Caiming Xiong, Richard Socher, and Byron C. Wallace.  ERASER: A benchmark to evaluate rationalized NLP models.  In Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics, ACL 2020, Online, July 5-10, 2020, pages 4443–4458. Association for Computational Linguistics, 2020.
  * DUA+ [22] Jessica Dai, Sohini Upadhyay, Ulrich Aivodji, Stephen H Bach, and Himabindu Lakkaraju.  Fairness via explanation quality: Evaluating disparities in the quality of post hoc explanations.  In Proceedings of the 2022 AAAI/ACM Conference on AI, Ethics, and Society, pages 203–214, 2022.
  * HEKK [19] Sara Hooker, Dumitru Erhan, Pieter-Jan Kindermans, and Been Kim.  A benchmark for interpretability methods in deep neural networks.  In Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada, pages 9734–9745, 2019.
  * LEC+ [20] Scott M. Lundberg, Gabriel G. Erion, Hugh Chen, Alex J. DeGrave, Jordan M. Prutkin, Bala Nair, Ronit Katz, Jonathan Himmelfarb, Nisha Bansal, and Su-In Lee.  From local explanations to global understanding with explainable AI for trees.  Nat. Mach. Intell., 2(1):56–67, 2020.
  * LKWN [21] Yang Liu, Sujay Khandagale, Colin White, and Willie Neiswanger.  Synthetic benchmarks for scientific research in explainable machine learning.  In Proceedings of the Neural Information Processing Systems Track on Datasets and Benchmarks 1, NeurIPS Datasets and Benchmarks 2021, December 2021, virtual, 2021.
  * LL [17] Scott M Lundberg and Su-In Lee.  A unified approach to interpreting model predictions.  In Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
  * PDS [18] Vitali Petsiuk, Abir Das, and Kate Saenko.  RISE: randomized input sampling for explanation of black-box models.  In British Machine Vision Conference 2018, BMVC 2018, Newcastle, UK, September 3-6, 2018, page 151. BMVA Press, 2018.
  * PMT [18] Gregory Plumb, Denali Molitor, and Ameet Talwalkar.  Model agnostic supervised local explanations.  In Advances in Neural Information Processing Systems 31: Annual Conference on Neural Information Processing Systems 2018, NeurIPS 2018, December 3-8, 2018, Montréal, Canada, pages 2520–2529, 2018.
  * RSG [16] Marco Túlio Ribeiro, Sameer Singh, and Carlos Guestrin.  ”Why should I trust you?”: Explaining the predictions of any classifier.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, San Francisco, CA, USA, August 13-17, 2016, pages 1135–1144. ACM, 2016.
  * Tor [23] Luis Torgo.  California housing dataset.  https://www.dcc.fc.up.pt/~ltorgo/Regression/cal_housing.html, 2023.  Accessed on: October 10, 2023.
  * YHS+ [19] Chih-Kuan Yeh, Cheng-Yu Hsieh, Arun Suggala, David I Inouye, and Pradeep K Ravikumar.  On the (in) fidelity and sensitivity of explanations.  Advances in Neural Information Processing Systems, 32, 2019.

##  Appendix A Models training

The grid search was performed by using class $GridSearchCV$ from library
XGBoost [6] with numbers of folds in cross-validation set to $5$. Additionally
to parameters described in 4.2, the grid search included also parameters with
the following possible values:

  * •

$eta$ with possible values in $\\{0.01,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9\\}$

  * •

$subsample$ with possible values in
${0.01,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9}$

##  Appendix B Additional experiments results

###  B.1 Results from precision experiments

Figure 7: Single tree model for Wine Quality dataset

Figure 8: Bigger model for Wine Quality dataset

Figure 9: Single tree model for Californian Housing dataset

Figure 10: Bigger model for Californian Housing dataset

###  B.2 Results from PGI comparison experiment

Figure 11: Single tree model for the Wine Quality dataset.

Figure 12: Bigger model for the Wine Quality dataset.

Figure 13: Single tree model for Californian Housing dataset

Figure 14: Bigger model for Californian Housing dataset


# On the Roughness of Two-Dimensional Materials

Pengjie Shi, Zhiping Xu  Corresponding author<EMAIL_ADDRESS>(Z.X.)
Applied Mechanics Laboratory, Department of Engineering Mechanics, Tsinghua
University, Beijing 100084, China

###### Abstract

This study examines the roughness of mechanically cleaved edges in 2D crystals
and glasses using molecular dynamics simulations with chemically accurate
machine-learning force fields. Our results show that ultra-flat armchair and
zigzag edges can be achieved in graphene by aligning the loading direction
with specific lattice orientations. Deviations from these orientations create
kinks between the atomically smooth armchair and zigzag segments, with
increased irregularities when dynamic effects are considered. Fracture
mechanics analysis highlights the kinetic and dynamic factors contributing to
crack deflection and edge roughening. In three-atom-thick 2D silica crystals,
the relationship between edge morphologies and cleavage conditions is modified
by their bilayer structure and sublattice asymmetry. In 2D silica glasses,
this correlation is further disrupted by topological disorder. These insights
are crucial for minimizing edge roughness in 2D materials, which is essential
for their performance in mechanical and electronic applications.

## Introduction

The roughness of a surface or an edge is an important measure of its quality.
Pursuing the lowest roughness is one of the eternal interests in materials
science and engineering, which finds its significance in a wide spectrum of
applications from tribology to microelectronics. Atomistically smooth surfaces
and edges offer structural superlubricity by harnessing their
incommensurability [1], templates for epitaxial growth of high-quality
semiconductors [2], and electrical contacts approaching the quantum limit [3,
4].

The roughness of materials is shaped by their growth processes and fabrication
techniques. The physics behind the roughening or smoothing processes can be
understood from thermodynamic or kinetic arguments. From the thermodynamic
perspective, surfaces or edges with the minimum energy densities are preferred
during its growth process, provided the atoms can diffuse on them and reside
on to the most energy-favorable sites. However, the surface tension of finite-
sized crystals typically induces faceting and the formation of polyhedrons
governed by Wulff’s construction principles, resulting in steps and kinks that
are thermodynamically or kinetically controlled [5]. In post-growth material
processing, mechanical polishing and physical or (electro)chemical etching are
widely used to create smooth surfaces and edges down to the Ångstrom level
[6].

Recent findings indicate that even minor edge defects in 2D crystals can
significantly degrade material performance, thereby limiting their promising
applications in elastic strain engineering (ESE) [7, 8]. Consequently, kinetic
control of edge roughness has garnered considerable attention. Following
principles of fracture mechanics, controlled cleavage produces atomically
smooth edges in 2D crystals [9]. However, the underlying physics of achieving
minimum roughness remains to be clarified. It has been argued in the fracture
mechanics community that a substance is considered truly brittle when its
fracture toughness is twice the specific surface free energy of the solid,
that is, $G_{\rm c}=2\gamma$. A ductile crystal fails by plastic shear,
through the motion of dislocations, for example, and the measured values of
$G_{\rm c}$ can rise to as much as $10^{5}\gamma$ [10]. However, the value of
$G_{\rm c}$ can be significantly larger than $2\gamma$ for the non-equilibrium
nature of fracture, even for extremely brittle crystals that fail to be bond-
breaking where a non-equilibrium measure $\Gamma$ is needed [11, 12]. Kinetic
control plays a crucial role.

Notably, for crystalline surfaces or edges with well-defined orientations, the
atomic flatness at a length scale much larger than the lattice constants is
not assured. For a 2D crystal in the honeycomb lattice, a ‘chiral’ edge is
composed by armchair and zigzag segments (Fig. 1a). The edge roughness is
defined by the numbers of consecutive segments $n_{\rm A}$ and $n_{\rm Z}$ and
can be significant if they are large, while the orientation depends only on
the ratio $n_{\rm A}/n_{\rm Z}$ (Fig. 1b). In contrast, in glassy materials
with a random bonding network, the minimum roughness can be lower than the
length of a single bond (Fig. 1c). However, the formation of local defects,
such as voids, can increase the roughness (Fig. 1d). For these reasons, it is
of practical interest to explore the lower bounds of or the minimum roughness
in crystalline and glassy materials, which can guide the preparation of
atomistically smooth surfaces with certain orientations.

In this study, we performed full-atom molecular dynamics (MD) simulations to
investigate the edge roughening mechanisms in 2D materials, motivated by the
simplicity of atomic-level structures and the recent surge in interest. By
resolving the evolution of atomic-level structures during the cleavage process
at the first-principles level, we elucidated the thermodynamic and kinetic
factors that govern the morphology and roughness of the resulting edges. This
research provides valuable insights into the fracture mechanics of 2D
structures and advances material fabrication techniques aimed at achieving
extreme or specified flatness.

## Results

Under stretch, shear, or tearing, graphene edges can be mechanically cleaved
along specific macroscopic directions. Considering the cleavage process as a
problem of fracture, the edge roughness is controlled by the local kinetics of
bond distortion and breakage, which includes crack deflection by mode mixing
or microstructural heterogeneity, and maybe crack branching triggered by
dynamical stability.

To explore the roughening process of 2D crystals at the atomic scale, MD
simulations are performed using neural-network force field for fracture
(NN-F3) trained by density functional theory (DFT) calculations, which allows
the use of large-scale models that are required for fracture mechanics
studies. 2D models of $50$ nm $\times$ $15$ nm and $100$ nm $\times$ $30$ nm
are constructed for graphene and silica, respectively (Fig. 2a-d). The model
sizes are comparable in the number of unit cells since the lattice constants
of 2D silica are approximately twice that of graphene. Edge pre-cracks are
created perpendicular to the loading directions by assuming the mode-I
condition. For crystalline graphene, we find that the edge morphology evolves
during propagation, which is characterized by the sequence ($n_{{\rm
A}i}$,$n_{{\rm Z}i}$). The stabilization of the ($n_{{\rm A}i}$,$n_{{\rm
Z}i}$) patterns indicates that the size or boundary effects can be precluded
in the discussion, which allows us to assess the kinetic control of edge
roughness (Fig. 2e,f).

The simulation results show that, under uniaxial tension along the close-to-
armchair direction ($\theta_{\mathrm{Z}}<9.82^{\circ}$), only zigzag edges are
cleaved. Conversely, for tension along the close-to-zigzag directions
($\theta_{\mathrm{Z}}>25.87^{\circ}$), only armchair edges cleave. The cleaved
edges are thus atomistically smooth with a roughness of $0$. Beyond these two
boundaries, cracks deflect from the zigzag or armchair edges via nucleating
lattice kinks due to the rise of a mode-II component of the stress field that
twists the crack tip (Fig. 2f). The strength of this twist is quantified by a
mode-mixing factor, $M=K_{\rm II}/K_{\rm I}$, and the critical values $M_{\rm
c}$ can be identified from the simulation results (Fig. 2e).

As the elastic driving force $G$ exceeds the fracture resistance $G_{\rm c}$
of the lattice, cracks start to grow. In this condition, the relation between
the direction of crack edge, $\theta_{\rm Edge}$, and $\theta_{\mathrm{Z}}$
can be determined from the maximum energy release rate (MERR) criterion as
[12]

$\theta_{\rm Edge}=\arg\min\frac{G_{\rm c}\left(\theta_{\rm
Edge}\right)}{\cos^{2}\left(\theta_{\mathrm{Z}}-\theta_{\rm Edge}\right)},$
(1)

where $G_{\rm c}\left(\theta_{\rm Edge}\right)=2\Gamma_{\rm
A}\left[\sin\left(\theta_{\rm
Edge}\right)+A_{G}\sin\left(30^{\circ}-\theta_{\rm Edge}\right)\right]$ is the
orientation-dependent fracture toughness. The anisotropic factor of
$A_{G}=\Gamma_{\rm Z}/\Gamma_{\rm A}=0.96$ is calculated by the non-
equilibrium edge energy density $\Gamma$ as a practice measure of the fracture
roughness ($G_{\rm c}=2\Gamma$) [12]. For
$\theta_{\mathrm{Z}}\in[9.82^{\circ},20.63^{\circ}]$, $n_{{\rm Z}i}=1$, with
$n_{{\rm A}i}\geq 1$. For
$\theta_{\mathrm{Z}}\in[20.63^{\circ},25.87^{\circ}]$, $n_{{\rm A}i}=1$, with
$n_{{\rm Z}i}\geq 1$ (Fig. 3a,c). The predicted boundaries $\theta_{\rm
Z}=9.9^{\circ}$ and $\theta_{\rm Z}=25.0^{\circ}$ by scanning the solutions of
the MERR criterion (Eq. 1) agree well with the simulation results
($\theta_{\rm Z}=9.82^{\circ}$ and $\theta_{\rm Z}=25.87^{\circ}$).
Interestingly, the cracks proceed along the deflected direction only for a
single step before the mode-mixing factor is effectively reduced, before
$M_{\rm c}$. These results with $n_{{\rm A}i}=1$ or $n_{{\rm Z}i}=1$ along the
roughened edges imply a one-to-one correspondence between crack direction and
edge roughness. Specifically at $\theta_{\mathrm{Z}}=20.63^{\circ}$, we have
both $n_{{\rm A}i}=1$ and $n_{{\rm A}i}=1$. The crack advances along the
direction of $\theta_{\rm Edge}=19.11^{\circ}$ in this condition, resulting in
a local minimum of the edge roughness in the range of
$\theta_{\mathrm{Z}}\in[9.82^{\circ},25.87^{\circ}]$.

MD simulations of 2D crystalline silica reveal trends in the sequence
($n_{{\rm A}i}$,$n_{{\rm Z}i}$) (Figs. 2g,h and 3b,c). However, the three-
atom-thick structure results in more complex edge morphologies (Fig. 1c). At
the tip of a propagating crack, the two Si-O bonds experiencing the greatest
forces can be cleaved in both the upper and lower layers. This asymmetry in
bond cleavage and the resulting lattice distortion contribute to rougher edges
compared to the atomistically sharp edges in graphene (Fig. 3e-h). The bilayer
structure and the presence of dangling bonds in broken Si-O-Si units further
exacerbate atomic-level structural irregularities. Consequently, the limits
$n_{{\rm A}i}=1$ or $n_{{\rm Z}i}=1$ are rarely achieved, leading to perturbed
directional dependence of edge roughness with varying loading directions.
Quantitatively, the anisotropy in fracture toughness for 2D crystalline silica
is more pronounced than in graphene, with $\Gamma_{\mathrm{A}}=1.004$ eV/Å,
$\Gamma_{\mathrm{Z}}=0.906$ eV/Å, and $A_{G}=0.90$ from our NN-F3
calculations, resulting in a shift of the low (high) roughness regions. The
predicted boundaries ($\theta_{\rm Z}=13.8^{\circ}$ and $28.1^{\circ}$) from
the MERR criterion (Eq. 1) agree well with the simulation results (Fig. 3c).

To turn on the inertial effects, we perform MD simulations with loading rates
in the range of $10^{-3}$-$10^{-5}$ ps-1 (Fig. 4a-f). The roughness of cleaved
edges is modified in the dynamic fracture process. An increase in the strain
rate leads to a higher line density of kinks along the edges, while the
overall direction of these edges remains almost unchanged. The conclusion
holds for 2D crystals of both graphene and silica, aligning with the common
knowledge that the inertial effect induces a mixed-mode correction at the
crack tip [13]. However, our discussion on the dynamics will not be elaborated
since an in-depth discussion needs much larger models where crack
instabilities can be triggered and observed.

Compared to 2D crystals, the disordered bonding network in 2D glassy silica
leads to amplified irregularities along cleaved edges. The crack paths are
shaped by lattice kinks and the disorder trapping effect, often guided by
voids forming ahead of the crack tip (Fig. 2i,j) [14]. These voids are
nucleated due to stress concentration at the crack tip ($\sim r^{{-1/2}}$) and
heterogeneity in the material strength due to the presence of lattice defects
or topological ones (e.g., pentagons and heptagons). Compared to the 2D silica
crystals, crack deflection in the glassy structures is less sharp and often
strongly distorted. This irregularity in the edge profiles indicates a
statistical nature of material strength. The roughness of 2D glasses increases
with the level of disorder, quantified by the standard deviation of ring
distributions,
$s=\sqrt{E\left[P^{2}\left(n\right)\right]-E^{2}\left[P\left(n\right)\right]}$.
Here, $n$ is the number of atoms in the minimal ring, and $E$ is the mean of
its distribution $P\left(n\right)$. For 2D silica, a transition from the
crystalline to glassy states occurs at $s\approx s_{\rm c}=0.45$, beyond which
the alignment between the loading direction and lattice orientations becomes
irrelevant. At low $s$, the crack velocity approaches a constant, indicating
that the effect of lattice trapping is negligible. However, for $s>s_{\rm c}$,
cracks can be trapped by the disorder. As disorder increases, the
instantaneous crack speed exhibits a distinct distribution pattern. At higher
disorder levels, the cumulative probability density function (CDF) of the
crack velocity converges to an exponential distribution (Fig. 4g,h),
indicating a less correlated nature.

## Discussion and Conclusion

The edge roughness (_Ra_ , the arithmetic average of absolute height
deviations measured from the mean line, Fig. 3) in the 2D crystals and glasses
depends on the loading direction and the lattice orientation. Understanding
these dependencies allows for better control of roughness during mechanical
cleavage. Firstly, selecting the direction of uniaxial tension is crucial in
determining the edge roughness of 2D crystals and glasses with low disorder.
By aligning the loading direction with low-energy edge motifs (zigzag or
armchair), one can achieve nearly atomically flat edges within a narrow range.
Outside this range, low roughness edges can be formed by lattice kinks joining
the armchair and zigzag segments, with roughness decreasing as kink density
increases. Secondly, dynamic fracture processes such as crack deflection and
bifurcation should be avoided at high strain rates. At low loading rates,
minimizing inertial effects ensures conditions for achieving minimum edge
roughness. Elastic anisotropy can be leveraged to surpass the fracture
resistance anisotropy inherent in the underlying crystalline structures. For
example, creating wrinkles in 2D crystals forces cracks to advance in the
direction with the highest stress intensity factor (SIF), perpendicular to the
wrinkling patterns [9]. Consequently, edges with specific orientations can be
produced. Notably, the evidence of $n_{\rm A}=1$ or $n_{\rm Z}=1$ under quasi-
static loading conditions links the macroscopic crack orientation to the
atomic-level edge structures and their associated roughness.

Our discussion extends to existing experimental facts of edge morphologies and
roughness reported for 2D materials. In graphene, edge morphologies and
roughness are closely tied to processing techniques such as peeling [15],
chemical vapor deposition (CVD) [16], and dry or gas-phase etching [17, 18].
Without controlling the loading conditions, the peeling method yields a
roughness of approximately $4.8$ nm. In contrast, the dry etching method can
achieve roughness as low as $1.5$ nm. Additionally, edge morphology can be
further disrupted by chemical asymmetry, as evidenced by the experimental
observation of branched fractured edges in hexagonal boron nitride (h-BN)
[19]. By controlling loading directions, loading rates, or elastic anisotropy,
atomically smooth edges can be achieved. Regarding edge generation experiments
for 2D glass, no relevant literature has been found to the best of our
knowledge. However, our results suggest that a low degree of disorder is
favored to suppress edge roughness.

In summary, we elucidate the origins of edge roughening during the mechanical
cleavage of 2D crystals and glasses. Chemically accurate neural network force
fields enable quantitative studies of atomistic fracture processes. The
anisotropy in fracture toughness of 2D crystals and the atomic-scale
heterogeneity in 2D glasses lead to localized deflection of crack tips.
Minimizing edge roughness can be achieved by carefully selecting loading
directions and mitigating dynamic effects, particularly relevant for
tribological and microelectronics applications. Although our findings apply to
3D solids in general, the presence of kinks along terraces of cleaved surfaces
introduces additional complexity. Notably, the simulation results also provide
valuable insights into fracture mechanics by delivering direct data on crack
deflection and bifurcation [11, 20, 21].

## Models and Methods

### Atomic-level structural models

To study the mechanical cleavage of 2D crystals and glasses, we choose
graphene and 2D silica (SiO2) as our model systems, which share the same
honeycomb lattice but different chemistry. A precrack is created along the $x$
direction. Displacement-controlled uniaxial tensile tests are then performed
along the $y$ direction. The lattice orientation is measured by an angle
$\theta$ between the $x$ direction and the nearest zigzag motif (Fig. 2a). The
cleaved edges are composed of armchair and zigzag segments. A sequence of
$\left(n_{\rm A1},n_{\rm Z1},n_{\rm A2},n_{\rm Z2},\cdots\right)$ is used to
describe the edges, where $n_{\rm A}$ and $n_{\rm Z}$ are the numbers of
consecutive armchair and zigzag structural units, respectively.

The atomic-level structures of 2D silica glasses are constructed by adapting a
dual-switch procedure [22], where the topological disorder is generated by
$10^{4}$ Monte Carlo steps with a fictitious temperature of $10^{-4}$ [23].
The product structures follow the Aboav-Weaire law with a model parameter of
$\alpha=1/3$ and agree with the experimental evidence, where small rings tend
to sit around larger ones [22].

### MD simulations

To capture the non-equilibrium nature of fracture processes, we use the
chemically-accurate neural-network force field for fracture (NN-F3) in our
molecular dynamics (MD) simulations [19, 12, 14]. In the workflow to construct
NN-F3 models, we pre-sample the 3D space of basal-plane strain states and
introduce an active-learning workflow to explore the bond-breakage and
(re)formation events during cracking using the end-to-end, symmetry-preserving
DeepPot-SE framework [24] and DeePMD-kit [25]. The pre-sampling technique
boosts the efficiency in exploring the high-dimensional space of strain
states, and the active-learning process ensures the accuracy of predictions by
enriching the training dataset. The technical details, scripts, and datasets
can be found in our recent reports [19, 12, 14]. The accuracy of NN-F3 is
confirmed to offer a root mean square error (RMSE) of $\leq 2$ meV/atom for
the energy per atom, $\leq 50$ meV/Å for atomic forces, and $\leq 30$ meV/Å2
for stress, respectively.

The generated NN-F3 are sufficiently accurate for the study of mechanical
cleavage. For a typical fracture problem of 2D fracture, the local processes
at the crack-tip and the far-field stress are concerned. As a result, the
cubic scaling of first-principles (e.g., density functional theory or DFT)
calculations hinders their applicability in simulating the fracture problem
where a significant size effect exists. In comparison, although the NN-F3
simulations are $\sim 100$ times slower than the empirical force fields (FFs,
e.g., Stillinger-Weber) for the models studied in this work, the computational
cost is much lower than DFT calculations, and the same linear scaling to the
number of atoms is maintained.

Quasi-static MD simulations are conducted to extract the stress-strain
relation and fracture patterns. The inertial effects are explored by MD
simulations at specific loading rates. The calculation of stress intensity
factors (SIFs) follows Refs. [26, 27], where the displacement field in an
annular region around the crack tip (see the inset of Fig. 2c) is fitted to
the Williams power expansion [28].

###### Acknowledgements.

This study was supported by the National Natural Science Foundation of China
through grants 12425201 and 52090032. The computation was performed on the
Explorer 1000 cluster system of the Tsinghua National Laboratory for
Information Science and Technology.

## Figures and Figure Captions

Figure 1:  (a) Crystalline edges of a honeycomb lattice, which are composed of
armchair (blue) and zigzag (red) segments. It should be noted that the same
edge orientation does not necessarily correspond to the same edge roughness
defined by $Ra$, the arithmetic average of absolute height deviations measured
from the mean line (b). (c) Typical edges of an amorphous lattice with
‘minimum’ roughness, displaying an irregular ring distribution. (d) Rough
edges of 2D glasses with irregular patterns including those created by voids
nucleated during fracture.  Figure 2:  (a) Uniaxial tension setup in molecular
dynamics (MD) simulations. (b-d) Typical crack-tip configurations in graphene
(b), 2D silica crystals (c) and glasses (d). (e-j) Mode-mixing factors
($K_{\rm II}/K_{\rm I}$, panels e, g, i) and edge morphologies (f, h, j) for
graphene (e-f), 2D silica crystals (g-h) and glasses (i-j). The dash lines
mark the range of critical mixing-mode factor $M_{\rm c}$, where the lower and
upper bounds are specifically annotated for graphene.  Figure 3:  (a-b)
Statistics of armchair and zigzag edge segments ($m$, $n$) the cleaved edges
in graphene (a) and 2D silica crystals (b), where the samples are cleaved
along direct orientation $\theta_{\rm Z}$. (c) Dependence of edge roughness
($Ra$) on the loading direction for graphene and 2D silica crystals. (d)
Dependence of edge roughness on the degree of disorder of 2D silica glasses,
defined as the standard deviation of the ring distribution,
$s=\sqrt{E\left[P^{2}\left(n\right)\right]-E^{2}\left[P\left(n\right)\right]}$.
(e,f) Edge morphologies of 2D silica crystals exhibiting irregularities due to
their bilayer structure and sublattice asymmetry.  Figure 4:  (a-f) Edge
morphologies of graphene (a-c) and 2D silica crystals (d-f) cleaved under
dynamic fracture. The loading rates are $10^{-3}$ (a,d), $10^{-4}$ (b,e), and
$10^{-5}$ (c,f), respectively. (g,h) Cumulative probability density function
(CDF) of the instantaneous crack velocities at different levels of disorder,
$s$ (g), and the deviation from an exponential distribution, $d$ (h).

## References

  * Hod  _et al._ [2018] O. Hod, E. Meyer, Q. Zheng, and M. Urbakh, Structural superlubricity and ultralow friction across the length scales, Nature 563, 485 (2018).
  * Zhao  _et al._ [2024] J. Zhao, P. Ji, Y. Li, R. Li, K. Zhang, H. Tian, K. Yu, B. Bian, L. Hao, X. Xiao, _et al._ , Ultrahigh-mobility semiconducting epitaxial graphene on silicon carbide, Nature 625, 60 (2024).
  * Wang  _et al._ [2011] X. Wang, Y. Ouyang, L. Jiao, H. Wang, L. Xie, J. Wu, J. Guo, and H. Dai, Graphene nanoribbons with smooth edges behave as quantum wires, Nat. Nanotechnol. 6, 563 (2011).
  * Li  _et al._ [2023] W. Li, X. Gong, Z. Yu, L. Ma, W. Sun, S. Gao, Ç. Köroğlu, W. Wang, L. Liu, T. Li, _et al._ , Approaching the quantum limit in two-dimensional semiconductor contacts, Nature 613, 274 (2023).
  * Pimpinelli and Villain [1999] A. Pimpinelli and J. Villain,  _Physics of Crystal Growth_ (Cambridge University Press, 1999).
  * Zhang  _et al._ [2023] Y. Zhang, Y. Zhang, K. Gu, L. Zhang, Y. Zhu, D. Liu, and H. Deng, A general and ultrafast polishing method with truly atomic roughness, J. Phys. Chem. Lett. 14, 9441 (2023).
  * Cao  _et al._ [2020] K. Cao, S. Feng, Y. Han, L. Gao, T. Hue Ly, Z. Xu, and Y. Lu, Elastic straining of free-standing monolayer graphene, Nat. Commun. 11, 284 (2020).
  * Han  _et al._ [2020] Y. Han, S. Feng, K. Cao, Y. Wang, L. Gao, Z. Xu, and Y. Lu, Large elastic deformation and defect tolerance of hexagonal boron nitride monolayers, Cell Rep. Phys. Sci. 1 (2020).
  * Feng  _et al._ [2023] S. Feng, W. Wang, S. Wang, X. Cui, Y. Yang, F. Xu, L. Liu, and Z. Xu, Controlling and visualizing fracture of 2D crystals by wrinkling, J. Mech. Phys. Solids. 174, 105253 (2023).
  * Kelly  _et al._ [1967] A. Kelly, W. Tyson, and A. Cottrell, Ductile and brittle crystals, Phil. Mag. 15, 567 (1967).
  * Lawn [1993] B. Lawn,  _Fracture of Brittle Solids_ (Cambridge University Press, 1993).
  * Shi  _et al._ [2024] P. Shi, S. Feng, and Z. Xu, Non-equilibrium nature of fracture determines the crack paths, Extreme Mech. Lett. 68, 102151 (2024).
  * Freund [1998] L. B. Freund,  _Dynamic Fracture Mechanics_ (Cambridge University Press, 1998).
  * Shi and Xu [2024a] P. Shi and Z. Xu, Strength of 2D glasses explored by machine-learning force fields, arXiv preprint arXiv:2407.13507 (2024a).
  * Qu  _et al._ [2022] C. Qu, D. Shi, L. Chen, Z. Wu, J. Wang, S. Shi, E. Gao, Z. Xu, and Q. Zheng, Anisotropic fracture of graphene revealed by surface steps on graphite, Phys. Rev. Lett. 129, 026101 (2022).
  * Tian  _et al._ [2011] J. Tian, H. Cao, W. Wu, Q. Yu, and Y. P. Chen, Direct imaging of graphene edges: Atomic structure and electronic scattering, Nano Lett. 11, 3663 (2011).
  * Zhao  _et al._ [2021] W. Zhao, H. Li, X. Xiao, Y. Jiang, K. Watanabe, T. Taniguchi, A. Zettl, and F. Wang, Nanoimaging of low-loss plasmonic waveguide modes in a graphene nanoribbon, Nano Lett. 21, 3106 (2021).
  * Wang and Dai [2010] X. Wang and H. Dai, Etching and narrowing of graphene from the edges, Nat. Chem. 2, 661 (2010).
  * Shi and Xu [2024b] P. Shi and Z. Xu, Exploring fracture of h-BN and graphene by neural network force fields, J. Phys. Condens. Matter 36, 415401 (2024b).
  * Liu and Wei [2021] Z.-E. Liu and Y. Wei, An analytical solution to the stress fields of kinked cracks, J. Mech. Phys. Solids 156, 104619 (2021).
  * Liu and Wei [2023] Z.-E. Liu and Y. Wei, A semi-analytical solution to the stress intensity factors of branched cracks, J. Mech. Phys. Solids 179, 105351 (2023).
  * Ebrahem  _et al._ [2020] F. Ebrahem, F. Bamer, and B. Markert, Vitreous 2D silica under tension: From brittle to ductile behaviour, Mater. Sci. Eng. A 780, 139189 (2020).
  * Font-Clos  _et al._ [2022] F. Font-Clos, M. Zanchi, S. Hiemer, S. Bonfanti, R. Guerra, M. Zaiser, and S. Zapperi, Predicting the failure of two-dimensional silica glasses, Nat. Commun. 13, 2820 (2022).
  * Zhang  _et al._ [2018] L. Zhang, J. Han, H. Wang, W. Saidi, R. Car, _et al._ , End-to-end symmetry preserving inter-atomic potential energy model for finite and extended systems, Adv. Neural Inf. Process. Syst. 31 (2018).
  * Wang  _et al._ [2018] H. Wang, L. Zhang, J. Han, and W. E, DeePMD-kit: A deep learning package for many-body potential energy representation and molecular dynamics, Comput. Phys. Commun. 228, 178 (2018).
  * Ayatollahi and Nejati [2011] M. Ayatollahi and M. Nejati, An over-deterministic method for calculation of coefficients of crack tip asymptotic field from finite element analysis, Fatigue $\&$ Fract. Eng. Mater. Struct. 34, 159 (2011).
  * Wilson  _et al._ [2019] M. A. Wilson, S. J. Grutzik, and M. Chandross, Continuum stress intensity factors from atomistic fracture simulations, Comput. Methods Appl. Mech. Eng. 354, 732 (2019).
  * Williams [1957] M. L. Williams, On the stress distribution at the base of a stationary crack, J. Appl. Mech. 24, 109 (1957).