# Efficient Low-Latency Dynamic Licensing for Deep Neural Network Deployment
on Edge Devices

Toan Pham Van<EMAIL_ADDRESS>R&D Lab, Sun* Inc.HanoiVietnam ,
Ngoc N. Tran<EMAIL_ADDRESS>R&D Lab, Sun*
Inc.HanoiVietnam ,  Hoang Pham Minh<EMAIL_ADDRESS>R&D Lab,
Sun* Inc.HanoiVietnam ,  Tam Nguyen Minh<EMAIL_ADDRESS>R&D Lab, Sun* Inc.HanoiVietnam and  Thanh Ta Minh<EMAIL_ADDRESS>Le Quy
Don Technical UniversityHanoiVietnam

(2020)

###### Abstract.

Along with the rapid development in the field of artificial intelligence (AI),
especially deep learning, deep neural network (DNN) applications are becoming
more and more popular in reality. To be able to withstand the heavy load from
mainstream users, deployment techniques are essential in bringing neural
network models from research to production. Among the two popular computing
topologies for deploying neural network models in production are cloud-
computing and edge-computing. Recent advances in communication technologies,
along with the great increase in the number of mobile devices, has made edge-
computing gradually become an inevitable trend. In this paper, we propose an
architecture to solve deploying and processing deep neural networks on edge-
devices by leveraging their synergy with the cloud and the access-control
mechanisms of the database. Adopting this architecture allows low-latency DNN
model updates on devices. At the same time, with only one model deployed, we
can easily make different versions of it by setting access permissions on the
model weights. This method allows for dynamic model licensing, which benefits
commercial applications.

Edge Computing, Artificial Intelligence, Dynamic Licensing, Deployment
Architecture

††journalyear: 2020††copyright: acmcopyright††conference: 2020 The 3rd
International Conference on Computational Intelligence and Intelligent
Systems; November 13–15, 2020; Tokyo, Japan††booktitle: 2020 The 3rd
International Conference on Computational Intelligence and Intelligent Systems
(CIIS 2020), November 13–15, 2020, Tokyo, Japan††price: 15.00††doi:
10.1145/3440840.3440860††isbn: 978-1-4503-8808-5/20/11††ccs: Software and its
engineering Software version control††ccs: Computer systems organization
Embedded software

##  1\. Introduction

Deep neural networks (DNNs) have become essential in real-world applications,
and thus started the quest for efficient DNN deployment architectures. There
are two popular computing topologies for this task: cloud computing and edge
computing. Cloud computing (Mell et al., 2011), a technology that has a long
history and numerous achievements, processes everything in the cloud, making
it easy for system upgrades, source code updates, and flexible scaling.
Deploying DNN models to the cloud would also take advantage of the powerful
computing power of these systems; and having only one centralized version for
all users would simplify weight updates and model performance monitoring.
However, cloud computing also has its limitations in terms of waste of system
resources, as it cannot make use of the computing power of edge devices. As a
result, the overall cost of maintaining a cloud system is very high. According
to the prediction of Ericsson, in 2024, 45% of global internet data will be
generated by the Internet-of-Things (IoT) devices (Report, 2019). Transferring
a great amount of data back-and-forth between edge devices and the cloud is
intractable as it causes excessive strain to the network infrastructure.
Furthermore, with the development of hardware technology, edge-devices have
more and more powerful computing capabilities. For that reason, moving from
cloud to edge is historically inevitable.

Instead of cloud computing, pushing DNNs to the edge of the network enhances
both efficiency and personalization. The AI model is still trained on the
cloud to utilize its great computing power. Afterward, the trained model gets
compressed to work with more computationally-limited hardware, then deployed
to the edge devices. This setting allows model updates to be distributed
quickly, simultaneously on all devices, and while trying to save as much
resources as possible. Besides, publishing a model with different licenses is
also an interesting issue for commercial applications, where a model can be
shipped with different performances depending on the user’s license. In this
paper, we propose a deployment architecture to deploy the DNN model to edge-
devices with four advantages as follows:

  * •

Efficient: In our architecture, we use an in-cloud database to store the DNN
model weights with their corresponding layer names and indices. This design
choice helps us centrally manage different versions of weights and easily
pushes or updates the changes of DNN model. As a result, the deployment and
update process becomes more efficient.

  * •

Low-Latency Update: Downloading all the weights of the DNN model can take a
long time on edge devices. This operation is only needed the first time when
we load model weights onto the edge devices. However, for later updates, we
might not need to download all the model weights but only the modified ones.
Based on such analysis, we propose an algorithm to download only modified
weights to achieve low-latency updates.

  * •

Version Management: Managing versions is a common task when deploying a DNN
model in production. Our method provides a solution to track changes with
commit history, update new versions, and rollback to an older version, similar
to Version Control Systems.

  * •

Dynamic Licensing: This feature is especially meaningful for commercial
applications. Our method allows the model owner to create unlimited licenses
of their model with different accuracies but while requires only one set of
weights saved on the database.

The rest of the paper is organized as follows. Section 2 provides a brief
survey of related works. Next, in Section 3, we describe our proposed method.
Experiment settings, including dataset information and data processing
methods, are described in Section 4; while experimental results are detailed
in Section 5. Finally, Section 6 presents our conclusions.

##  2\. Related Works

###  2.1. Edge-based vs Cloud-based AI

Currently, model inferences are mostly performed in the cloud; but as the
diversity of DNN applications grows, other alternatives to the centralized
training and inference strategy are required to lessen the burden on network
infrastructures (Li et al., 2018). To that, we opt for Edge Computing (Shi and
Dustdar, 2016): a distributed computing paradigm where software-defined
networks are built to decentralize data and provide results expected to be the
same as which of cloud computing (Shi et al., 2020). Solving these problems,
edge DNN aims to process DNN models directly on edge devices. However, to
utilize their innate computing power, edge computing faces more resource
allocation problems, due to the inherent difference in hardware architecture,
the need to sample inputs from built-in peripherals, bandwidth constraints and
more (Shi et al., 2016, 2020). To offset these problems, deployment on edge-
devices requires DNN models to be minimal enough for fast inference and low-
latency updates. As a result, it requires various powerful optimizations to
achieve the required system efficiency. The difference between cloud-based and
edge-based for DNN application is demonstrated in Figure 1. The left figure
shows the cloud-based deployment with training, and the right figure shows the
edge-based architecture where inference has been offloaded to the edge
devices.

Figure 1. Cloud-based DNN (left) vs Edge-based DNN architecture (right)

###  2.2. Deep Neural Network (DNN)

A deep neural network is, in essence, a composition of more than two non-
linear, simpler functions, which can approximate any arbitrarily complex
function (Goodfellow et al., 2016) (Universal Approximate Theorem (Csáji,
2001)). Mathematically, a DNN is defined as:

(1) $f_{\theta(x)}=f_{\theta_{n}}\circ f_{\theta_{2}}\circ f_{\theta_{1}}(x),$

where

$\theta=\left\\{\theta_{1},\theta_{2},\ldots,\theta_{n}\right\\}$

is its parameterization.

The most common task is classification, where the model would output whether a
data is of a particular class. To build a DNN for the job, we first choose a
model architecture that can find meaningful patterns to output. Subsequently,
all that is left to do is to learn the parameters of that network through
training.

To train a DNN model, a scalar objective function must be defined to quantify
the quality of the model; and thus optimizing this function will improve the
model’s performance as a result. In the case of classification, we would
minimize MLE as our loss function. Optimizing this loss function using
gradient-based methods (Zhang, 2019) is a common practice in deep learning
algorithms. This “allows the information from the cost to flow backward
through the network to compute the gradient” (Goodfellow et al., 2016). The
optimal parameters (weights) are approximated via the iteration method.

###  2.3. Model Compression Methods

Modern DNNs are expensive in both computation and memory storage, which leads
to difficulties in deployment on mobile devices. Consequently, model
compression – compressing large, complex models into a lighter, simpler one
without significant loss in accuracy – has become indispensable to the DNNs
deployment.

####  2.3.1. Model Pruning

An enormously trained model in deep learning contains a large amount of
redundancy (Cheng et al., 2017) in the form of unimportant weights that have
little contribution to the final output. Pruning is a method of model
compression, lightening the architecture by cutting off those unimportant
connections or weights, trading a minor loss in quality for performance (Cheng
et al., 2017). Methods of identifying uninformative weights are varied: from
using second-order derivatives information (Cun et al., 1990), similarity in
neurons (Srinivas and Babu, 2015), to the state-of-the-art Movement Pruning
(Sanh et al., 2020) which use the change’s magnitude in weight updating. Among
these above complex methods, one stands out with its simplicity and
effectiveness: Magnitude Pruning (Han et al., 2016) – eliminating weights
which have magnitudes less than some threshold, followed by fine-tuning to
achieve the highest possible accuracy on the (significantly) pruned model.

####  2.3.2. Quantization

Quantization is the method of reducing the representation accuracy of the
network’s weights by storing them with lower numbers of bits. This leads to
faster computation, decreases in-memory storage with a trade-off of lower
accuracy. 8-bit weight quantization (Vanhoucke et al., 2011) or 16-bit fixed-
point with stochastic rounding representation (Gupta et al., 2015) are well-
known examples of this method, for their effectiveness heavily outweighs their
loss of accuracy (Cheng et al., 2017).

####  2.3.3. Weight Sharing

Deep Compression (Han et al., 2016) also introduces weight-sharing that works
well with quantization. This method can be applied after quantization. It
divides the weights into $k$ clusters, where weights in one cluster are closed
in magnitude, which afterward will be set to be identical. Consequently, a
reduction in storage and computation occurs since one only needs to store a
sparse cluster-index matrix and a hashtable for quantized value-lookup instead
of a dense matrix.

###  2.4. Database and Query

####  2.4.1. Database

A database is a collection of data, typically describing the activities of and
among related entities (Ramakrishnan and Gehrke, 2000). For different problem
requirements and types of information stored, various types of databases are
available for the task: relational database, NoSQL database, graph database,
etc. (TutorialsPoint, [n.d.]). The most mature and widely used database
systems in production today are relational database management systems, which
can be found in most applications, such as e-commerce, social networks,
retails, etc.

DNNs consist of neurons organized into layers of neurons and the connections
between them. This hierarchical architecture encourages us to represent a DNN
with a relational database. Storing a DNN in the database would make it easy
to update weights, ensuring weight constraints, and manage weight access
permissions.

####  2.4.2. GraphQL

Querying from a database fast and efficiently is impossible with the
traditional RESTful API (Gao et al., 2011), so we opted for a better choice
for the job: GraphQL (Taelman et al., 2018). For all its popularity’s worth,
one aspect where REST API falls short is its strict inflexible specifications,
where custom requests for various information would require a lot of queries
to different API endpoints. Not only this would require a complicated control
flow to support many types of requests, but also combined latency of all these
queries introduces a noticeable lag. Moreover, if one endpoint were to fail,
the request would not be able to go through later endpoints – this
bottlenecking phenomenon is common with REST API. However, these problems do
not occur with GraphQL, being a client-oriented query language with a flexible
data structure that can adapt to any type of client’s demands. GraphQL
supports operations similar to REST API, but with a hierarchical structure
that is client-friendly.

##  3\. Proposed Method

###  3.1. Our Architecture

####  3.1.1. Combining Edge-Devices with weights storage database

The traditional Edge-AI architecture (mentioned in subsection 2.1) only has
two main partitions: edge devices and the cloud. The cloud in this computing
topology has the roles of training and storing DNN models, managing model
versions, storing user data, and model licenses. That is a huge amount of
tasks to manage; especially with the model licensing, where the server needs
to store a large number of model versions. Imagine we have 10 versions of the
model, each of which has 10 different licenses: the cloud would have to store
100 copies of the model to accommodate that. This leads to a waste of system
resources. Therefore, we propose our new pipeline as demonstrated in Figure 2.
Our proposed architecture as the same as an additional module of traditional
edge computing in deployment the DNN models. Instead of store all weight and
DNN topology in a server, we store the model weights in a database placed on a
different cloud. By splitting the traditional single unified cloud into two
with their dedicated tasks, we both lessen the workload and let each server do
one task and do it well. Moreover, we only store incremental changes across
model versions and only license-specific indices of the production model
version through our permission management system. Our version management is
further elaborated in subsection 3.4.

Figure 2. Our architecture with weight storage in database

####  3.1.2. Updating Weight Versions

Our system enables a flexible database weight update when there is a newer
version available. Our procedure is as follows: we check the existence of each
weight variable of the new version and whether its value changed. If a
parameter is not yet in the table or its value is different from the old one,
a new entry in the Weight table will be created to store it. Otherwise, no
update in the database will be made. This simple strategy allows significantly
more time and storage efficiency in weight-updating, compared to the
traditional method of storing parameters of each version in a separate table.

In the other direction, when the database has a newer version than which on
the edge device, the user may want to update his version to the newest one
(possibly as a response to a push notification). In this situation, the device
starts the update by sending our storage server the current model version it
has. Then, the server responds with values and indices of the weights that are
either newly created or updated. This process is very similar to the previous
one, but the benefits of not fetching unchanged weights to reduce extra
bandwidth are more noticeable: this be analyzed in more detail later.

###  3.2. Model Compression Strategy

To deploy machine learning models to mobile devices, we lessen the burden of
computation and storage by compressing our models before storing them in the
database.

Figure 3. Model compression pipeline

Figure 3 describes the pipeline of our model compression: Firstly, we apply
magnitude pruning on our models to cut off uninformative weights.
Subsequently, we fine-tune models to obtain optimal weights on the pruned
architecture. Next, quantization (converting weights from 64-bit to 8-bit
representation) followed by weight sharing is used to further compress the
models. Despite its simplicity, this pipeline results in highly effective
model compression with acceptable performance loss.

###  3.3. Weight Database Storage

After compressing models using our compression pipeline, we store models’
weights into Postgres — our database of choice. In order to connect to the
Postgres database, we opt for GraphQL API instead of REST API due to its
advanced features that allow fast and precise data access. Particularly,
Hasura Engine, a GraphQL server, is utilized connect to Postgres in realtime.
We use Hasura for querying data from the database, and Django for adding data
to the database. To harmonize with GraphQL syntaxes, we design our model
weight database with tables for Model, Layer, Weights, Version, and Accuracy,
as illustrated in Figure 4.

Figure 4. Weight Database Schema

After reducing the model’s size with compression, we save the optimized
weights in our database along with their metadata: the weight’s flattened
index, weight value, date added, corresponding layer and version table
entries’ foreign keys. Weight arrays extracted from the model instance with
common deep learning frameworks are in the same format, being a dictionary
with layer names as keys and weight arrays as values, greatly streamlining our
storage procedure. After flattening the weight arrays, we save the weights in
the aforementioned format. Thanks to model pruning, we only need to save the
nonzero weights of the resulting sparse matrices, greatly saving our
database’s space. This process is time-efficient thanks to Django Object-
Relational Mapper’s (ORM) bulk insert mechanism.

Conversely, to load weight values from the database, we first build an empty
(zeroed) model with its architecture layer-by-layer with their corresponding
Layer’s table entries. Then, we load the weight values and place them in the
layers individually in the locations indicated by their flattened indices.

###  3.4. Weight Version Control

During the design process, we might come up with with various versions for one
model. For example, during model selection, one may try cross-validation to
choose over a range of possible hyperparameters. Changing a hyperparameter set
would change every weight in that model, and the new version’s record in the
database would be a completely new entry. Another possibility is fine-tuning
models while freezing layers to incorporate minor improvements: in this case
only a subset of weights would change. In this case, the new version only
stores the indices and values of the changed weights.

Model versions and their parameters are stored in the Version and Weight
tables respectively. In the Version database, there is a status field of the
Boolean type which determines whether the major version is used for
production. These major versions are completely different from each other, so
only one is in use at a time. Then, for minor versions of the production
model, each parameter component is linked with its corresponding version ID
for efficient weight updates. Specifically, when a model gets an incremental
update, only the changed weights are stored in new entries, while the others
remain unchanged in the Weight database. This improves model retrieving since
updating a new version only requires a fraction of the weights to be
retrieved.

Accuracy table is connected with the Version, Layer, Model tables, and most
importantly, the Weight table. In Accuracy, for a particular version accuracy
value, we save a set of weight value ranges to be used that would give us the
desired performance. The idea is that we artificially cause worse model
performance by not using weights in some particular magnitude ranges. With
information from the requesting user’s license, a version with the
corresponding accuracy will be shipped.

###  3.5. Weight Licensing

Given that machine learning models are our product, we need to have a customer
licensing mechanism for proprietary reasons. A tangible example would be as
follows: we offer a 3-layer-perceptron model with 98% accuracy through a
freemium business model. Using a permission control mechanism, we do not let
free-tier users to access weights of the first layers with magnitudes between
$0.5$ and $0.8$; instead, these weights would be set to 0 (similar to the
pruning process). The model’s accuracy then drops to 70% in our experiment.

The aforementioned Accuracy table contains masks over our deployed machine
learning model’s weights. The more weights the mask hides, the worst accuracy
a version has. To get weight indices corresponding to some specific desired
accuracy, all current deployed model’s weight values are divided into equal-
sized intervals based on their magnitude. We perform gradual magnitude pruning
on the model, results in a corresponding gradual accuracy reduction. The
pruning process terminates when the required accuracy is observed. Algorithm 1
describes our method of retrieving pruned models based on desired accuracy.

divide weight range into $k$ smaller equal-sized intervals

initialize a list of cut-off intervals

for all intervals do

for all model’s layers do

cut off weights that have values in that interval

append interval into cut-off interval list

if accuracy of pruned model is close to the target then

break the pruning process

end if

end for

end for

return uncut interval lists

Algorithm 1 Pruning model based on accuracy

For customers with licenses that fall into our predefined tiers, we use Static
Licensing: we access the Accuracy table, which is directly linked to the
Weights table and stores the weight ranges for every layer and their
corresponding accuracies. These weight range-accuracy relationships are
evaluated beforehand; and at deployment only lookups are needed to ship the
appropriate set of weights according to the user’s subscription tier. On the
other hand, if the client requires custom performance tier, we turn to Dynamic
Licensing, which evaluates these ranges on-demand.

##  4\. Main Features of Our Method

###  4.1. Efficient Deployment

As mentioned in subsubsection 3.1.1, we break the work that is traditionally
all processed in one cloud server into two highly coherent tasks. Even with
decent scaling solutions like Docker Swarm or Kubernetes, if the basic worker
instance has to work with both training and shipping new weights, it will be a
lot less efficient than having separate instances doing their respective
tasks.

###  4.2. Version Management

Our design stores weight incrementally, instead of a completely new whole
model every new version. With this, not only we make the previously analyzed
Efficient Deployment possible, we save storage by only keeping changes and not
the unchanged weights. Moreover, this allows for skipping intermediate
patches. Given that every weight entry in our database is stored along with
the version it was last updated, the customer can query for all new weights
throughout his missed updates in one go, instead of having to download
individual version updates and apply them one-by-one gradually.

###  4.3. Low-Latency Update

This approach’s contribution to fast evaluation is manifold: first, the
customer’s model evaluation is carried out within their device. This is
crucial for time-critical operations, such as self-driving cars, where
milliseconds matter. This is also important to always-on systems, such as
health gadgets with real-time illness detection – a model evaluation server
experiencing difficulty at the worst possible time may cost lives. Moreover,
the data for evaluation stays on the customer’s device, which is a huge plus
for the privacy-conscious. For us, this also offloads work to the client’s
device, which is getting more and more capable as technology develops. Second,
our setup allows fast weight updates from our database server. This would
greatly benefit interdisciplinarity that value timely updates of even
minuscule incremental performance, such as stock trading. For large hedge
funds, even just a 0.1% percent increase in accuracy shipped 1 millisecond
earlier could mean millions of dollars.

###  4.4. Dynamic and Static Licensing

Most businesses offer their service tiers in terms of speed and/or
availability, but here we introduce a new factor: accuracy. Most trial users
would want a not-too-accurate but fast free-tier service since their
applications would not need more than that, while their patience is limited.
Moreover, our design allows dynamic licensing which lets the clients custom
their version’s performance beyond our predefined tiers. This allows
flexibility to match whatever the most budget-conscious client needs.

##  5\. Experiments

###  5.1. Experimental Settings

For our experiments, we use the following technologies:

  * •

Django Framework: database connecter for inserting and updating data into
database with Django ORM (Forcier et al., 2008).

  * •

Keras (Team, [n.d.]), TensorFlow (Google, [n.d.]): deep learning frameworks
for building DNN architectures.

  * •

PostgreSQL (PostgreSQL, [n.d.]): a relational database backend,

  * •

Hasura (Hasura, [n.d.]): an instant GraphQL engine with built-in authorization
for querying weights,

  * •

Docker (Docker, [n.d.]): machine virtualization for simulating real case
studies in this paper.

To verify our system’s effectiveness, we measure the amount of space needed to
store weights in our database with gradual additions of various optimizations.
The results of these experiments are listed in Table 1.

###  5.2. Experimental Results

Table 1. The cost of memory storage No. of params | Full params | Pruning 80% |  Pruning 80% + Quantization
---|---|---|---
109386 | 13MB | 2.92MB |  2.34MB
101770 | 12MB | 2.65MB |  2.09MB

Table 1 shows our results in saving weights of models with various
optimizations. With more than 100,000 weights, saving all weights in Postgres
takes 13MB of space. Through pruning the unimportant 80% of the weights and
not storing their unused bits with a more concise representation, the storage
size reduces to 2.09MB.

##  6\. Conclusion

In this paper, we successfully design an architecture utilizing an In-cloud
database for efficient deployment and centralization management. The merits of
our procedure are additionally represented by low-latency update (which only
downloads modified weights on to edge devices), flexible version managed
(allows track on history commitment changes, new version update, old version
rollback), dynamic and static licensing (provides both ready-to-use and
custom-accuracy models). With this novel deployment architecture, we hope our
contributions will be beneficial to numerous DNN applications.

###### Acknowledgements.

This work is partially supported by Sun-Asterisk Inc. We would like to thank
our colleagues at Sun-Asterisk Inc for their advice and expertise. Without
their support, this experiment would not have been accomplished.

## References

  * (1)
  * Cheng et al. (2017) Yu Cheng, Duo Wang, Pan Zhou, and Tao Zhang. 2017\.  A Survey of Model Compression and Acceleration for Deep Neural Networks.  arXiv:1710.09282 [cs.LG]
  * Csáji (2001) Balázs Csanád Csáji. 2001\.  Approximation with Artificial Neural Networks. In _Faculty Sci. Etvs Lornd Univ. Hungary_.
  * Cun et al. (1990) Yann Le Cun, John S. Denker, and Sara A. Solla. 1990\.  _Optimal Brain Damage_.  Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 598–605.
  * Docker ([n.d.]) Docker. [n.d.].  Docker Documentation.  https://www.docker.com/
  * Forcier et al. (2008) Jeff Forcier, Paul Bissex, and Wesley J Chun. 2008\.  _Python web development with Django_.  Addison-Wesley Professional.
  * Gao et al. (2011) Lei Gao, Chunhong Zhang, and Li Sun. 2011.  RESTful Web of Things API in sharing sensor data. In _2011 International Conference on Internet Technology and Applications_. IEEE, 1–4.
  * Goodfellow et al. (2016) Ian Goodfellow, Yoshua Bengio, and Aaron Courville. 2016\.  _Deep Learning_.  The MIT Press.
  * Google ([n.d.]) Google. [n.d.].  Tensorflow Documentation.  https://www.tensorflow.org/
  * Gupta et al. (2015) Suyog Gupta, Ankur Agrawal, Kailash Gopalakrishnan, and Pritish Narayanan. 2015.  Deep Learning with Limited Numerical Precision.  arXiv:1502.02551 [cs.LG]
  * Han et al. (2016) Song Han, Huizi Mao, and W. Dally. 2016.  Deep Compression: Compressing Deep Neural Network with Pruning, Trained Quantization and Huffman Coding.  _CoRR_ abs/1510.00149 (2016).
  * Hasura ([n.d.]) Hasura. [n.d.].  Hasura Documentation.  https://hasura.io/
  * Li et al. (2018) He Li, Kaoru Ota, and Mianxiong Dong. 2018.  Learning IoT in edge: Deep learning for the Internet of Things with edge computing.  _IEEE network_ 32, 1 (2018), 96–101.
  * Mell et al. (2011) Peter Mell, Tim Grance, et al. 2011\.  The NIST definition of cloud computing.  (2011).
  * PostgreSQL ([n.d.]) PostgreSQL. [n.d.].  PostgreSQL Documentation.  https://www.postgresql.org/
  * Ramakrishnan and Gehrke (2000) Raghu Ramakrishnan and Johannes Gehrke. 2000.  _Database management systems_.  McGraw-Hill.
  * Report (2019) Ericsson Mobility Report. 2019\.  NB-IoT and Cat-M technologies will account for 45 percent of cellular IoT connections in 2024.  https://iotbusinessnews.com/2019/06/11/50915-nb-iot-and-cat-m-technologies-will-account-for-45-percent-of-cellular-iot-connections-in-2024
  * Sanh et al. (2020) Victor Sanh, Thomas Wolf, and Alexander M. Rush. 2020.  Movement Pruning: Adaptive Sparsity by Fine-Tuning.  arXiv:2005.07683 [cs.CL]
  * Shi et al. (2016) Weisong Shi, Jie Cao, Quan Zhang, Youhuizi Li, and Lanyu Xu. 2016.  Edge computing: Vision and challenges.  _IEEE internet of things journal_ 3, 5 (2016), 637–646.
  * Shi and Dustdar (2016) Weisong Shi and Schahram Dustdar. 2016.  The promise of edge computing.  _Computer_ 49, 5 (2016), 78–81.
  * Shi et al. (2020) Yuanming Shi, Kai Yang, Tao Jiang, Jun Zhang, and Khaled B Letaief. 2020.  Communication-efficient edge AI: Algorithms and systems.  _arXiv preprint arXiv:2002.09668_ (2020).
  * Srinivas and Babu (2015) Suraj Srinivas and R. Venkatesh Babu. 2015.  Data-free parameter pruning for Deep Neural Networks.  arXiv:1507.06149 [cs.CV]
  * Taelman et al. (2018) Ruben Taelman, Miel Vander Sande, and Ruben Verborgh. 2018\.  GraphQL-LD: linked data querying with GraphQL. In _ISWC2018, the 17th International Semantic Web Conference_. 1–4.
  * Team ([n.d.]) Keras Team. [n.d.].  Keras Documentation.  https://keras.io/about/
  * TutorialsPoint ([n.d.]) TutorialsPoint. [n.d.].  Types of databases.  https://www.tutorialspoint.com/Types-of-databases
  * Vanhoucke et al. (2011) Vincent Vanhoucke, Andrew Senior, and Mark Z. Mao. 2011\.  Improving the speed of neural networks on CPUs. In _Deep Learning and Unsupervised Feature Learning Workshop, NIPS 2011_.
  * Zhang (2019) Jiawei Zhang. 2019\.  Gradient Descent based Optimization Algorithms for Deep Learning Models Training.  arXiv:1903.03614 [cs.LG]


# Explainable AI models for predicting liquefaction-induced lateral spreading

Cheng-Hsi Hsiao 1,∗, Krishna Kumar 1 and Ellen Rathje 1

###### Abstract

## 1

Earthquake-induced liquefaction can cause substantial lateral spreading,
posing threats to infrastructure. Machine learning (ML) can improve lateral
spreading prediction models by capturing complex soil characteristics and site
conditions. However, the “black box” nature of ML models can hinder their
adoption in critical decision-making. This study addresses this limitation by
using SHapley Additive exPlanations (SHAP) to interpret an eXtreme Gradient
Boosting (XGB) model for lateral spreading prediction, trained on data from
the 2011 Christchurch Earthquake. SHAP analysis reveals the factors driving
the model’s predictions, enhancing transparency and allowing for comparison
with established engineering knowledge. The results demonstrate that the XGB
model successfully identifies the importance of soil characteristics derived
from Cone Penetration Test (CPT) data in predicting lateral spreading,
validating its alignment with domain understanding. This work highlights the
value of explainable machine learning for reliable and informed decision-
making in geotechnical engineering and hazard assessment.

##  2 Keywords:

Explainable AI, SHAP, XGBoost, Lateral spreading, 2011 Christchurch earthquake

##  3 Introduction

Lateral spreading, characterized by the horizontal movement of liquefied soil
triggered by earthquakes, can cause severe damage to infrastructure and
present significant hazards in urban environments. The movement results from
rapid soil transition from a solid to a liquid state during earthquake
excitation, influenced by the initial state of soil and amplitude of shaking.
While a site might undergo liquefaction, it will only experience lateral
spreading if the topography facilitates it. Given the complex factors
involved, advanced, data-driven approaches are essential for predicting the
extent and impact of lateral spreading. In recent years, machine learning (ML)
applications have shown promising results in modeling lateral spreading (Demir
and Şahin, 2022; Kaya, 2016; Durante and Rathje, 2021).

Although ML models extract relationships between the various influencing
factors and failure occurrence, they pose interpretability and trustworthiness
challenges (Miller, 2019; Doshi‐Velez and Kim, 2017). The term “black box”
refers to the opaque nature of these ML models, where the internal logic of
their decision-making processes remains hidden. This opacity can have
significant consequences, primarily when using ML models for policy decisions
related to natural hazards, where inaccurate predictions could lead to
substantial damage. Furthermore, ML models generally excel in interpolation,
meaning they perform well when predicting scenarios similar to their training
data. However, they may struggle with extrapolation, failing to accurately
predict outcomes in new, unseen conditions, posing challenges in their
generalizability and application to sites outside their training regime.

The complexity of lateral spreading exacerbates these challenges. Factors such
as site geometry, soil type, and loading conditions influence lateral
spreading, making predictions more complex. In ML models that do not consider
the physics of these processes, there is a risk of including redundant or
excessive features, termed feature coherence. These features could lead to
model overfitting, where the model memorizes behavior rather than learning the
general trend, limiting its applicability to different datasets (Paris et al.,
2004; Xue, 2019). Excessive coherence among features indicates that multiple
inputs might overlap in their influence, causing redundancy within the model’s
decision-making framework and further complicating the interpretation and
reliability of the model’s predictions.

To enhance the transparency and interpretability of ML models, researchers
have developed a variety of explainable AI (XAI) methodologies (Arrieta et
al., 2020; Linardatos et al., 2020). These approaches fall into two broad
categories: model-agnostic methods and model-specific methods. Model-agnostic
methods, such as Local Interpretable Model-agnostic Explanations (LIME;
Ribeiro et al., 2016) and SHapley Additive exPlanations (SHAP; Lundberg and
Lee, 2017), provide versatile techniques applicable across different types of
ML models. These methods generate explanations by observing how the model
output changes as it perturbs the input data. This process allows model-
agnostic explainers to approximate how each feature influences the model’s
output, offering insights into its behavior without requiring access to its
internal workings, making them universally applicable to any ML model.
Conversely, model-specific methods, like Layer-wise Relevance Propagation
(LRP; Bach et al., 2015), are tailored to specific model architectures. Such
XAI techniques play a pivotal role in demystifying the decision-making process
of machine learning models, highlighting influential features, and
facilitating the incorporation of domain knowledge for model improvement.

In this study, we utilize a model-agnostic SHAP approach to delve into the
decision-making process of an XGB (Chen and Guestrin, 2016) model designed to
predict the occurrence of lateral spreading. Our analysis examines the 2011
Christchurch earthquake dataset (Durante and Rathje, 2022). By employing SHAP,
we aim to improve the interpretability of our developed model by exploring
local and global explanations for each feature. SHAP allows us to uncover
complex relationships between input variables and model predictions. Notably,
SHAP helps identify instances where the impact of specific features goes
against established physics principles or domain knowledge. This approach
contributes to a deeper understanding of lateral spreading prediction and
enhances the accessibility of complex ML models for researchers and
practitioners.

The 2011 Christchurch Earthquake was a devastating event characterized by
severe liquefaction and lateral spreading Cubrinovski et al. (2011). These
phenomena, induced by the earthquake ground shaking, led to significant
horizontal movement of saturated soil, causing widespread damage to
infrastructure and posing substantial risks to urban areas.

Figure 1() illustrates the horizontal displacement amplitudes measured using
remote sensing techniques, as documented by Rathje et al. (2017b). The map
categorizes the lateral displacements into three zones: green signifies 0 to
0.3 meters, orange indicates 0.3 to 0.5 meters, and red represents
displacements exceeding 0.5 meters. Our study leverages the dataset developed
by Durante and Rathje (2022) on the Rathje et al. (2017b) measurements to
predict lateral spreading. This dataset encompasses 6,505 locations with
displacement data and co-located CPT profiles. We employ a displacement
threshold of 0.3 meters, as depicted by the orange and red regions in Figure
1(), to distinguish between the occurrence and non-occurrence of lateral
spreading. This threshold is based on research by Rathje et al. (2017b), which
demonstrated that displacements below this value could not be statistically
distinguished from areas with zero movement. Figure 1() shows the spatial
distribution of the Yes/No lateral spreading data based on the threshold
including 2752 Yes data (red triangle; 42.31%) and 3753 No data (green circle;
57.69%).

figure1  (a) Observed horizontal displacement (Rathje et al., 2017b) (b)
Lateral spreading occurrence (Durante and Rathje, 2022) (c) Observed
liquefaction-related damage (NZGD)

() Lateral spreading occurrence map of 2011 Christchurch earthquake. (A)
Lateral spreading displacement from remote sensing data (data from Rathje et
al., 2017b).(B) Spatial distribution of lateral spreading occurrence (data
from Durante and Rathje, 2022).(C) Observed liquefaction-related damage (data
from New Zealand Geotechnical Database, 2013)

figure2  (a) GWD (b) PGA (c) Distance to the river (d) Slope (e) Elevation

() Distribution of input features and outliers analysis. (A) GWD (B) PGA (C)
Distance to the river (D) Slope (E) Elevation

Similar to Durante and Rathje (2021), we develop an ML model to predict
lateral spreading based on five features: the distance between the site and
the nearest river location (L), ground slope, ground elevation, groundwater
depth (GWD), and peak ground acceleration (PGA). Our paper aims to develop a
comprehensive explanation of the ML model for lateral spreading occurrence
using the SHAP analysis.

Figure 2() shows the distributions of these predictive features. Each
histogram is accompanied by a box diagram summarizing the feature
distributions, with outliers identified as red dots using the interquartile
range (IQR) method. Specifically, observations exceeding 1.5 times the IQR
below the first quartile (Q1) or above the third quartile (Q3) are flagged as
outliers, with the IQR representing the difference between Q1 and Q3.

We notice that groundwater depth (GWD) and slope data (Figures 2() and 2())
vary significantly in the region of interest. Specifically, we have identified
117 data points with GWD values exceeding the median value of 3.70 meters and
500 data points with slope median values exceeding 2.84%. While these values
deviate significantly from the median value, they are characteristics of the
sites. Hence, we have opted to retain them in our analysis. Our rationale is
to enable our trained model to learn from a broader range of feature values,
enhancing its ability to generalize and make accurate predictions across
various scenarios.

In this study, we conducted a comprehensive SHAP analysis to gain insights
into the behavior of the XGB models developed for lateral spreading
prediction. We introduce the XGB and SHAP approaches in the Methodology
section. The Model Development section outlines the model configuration,
training procedures, and overall performance metrics. Finally, in the Results
section, we delve into the details of SHAP’s local and global explanations for
our trained XGB models. Through this SHAP analysis, we unveiled key insights
derived from our models and assessed the significance of the CPT features in
our predictive models.

##  4 Methodology

###  4.1 eXtreme Gradient Boosting

XGB (Chen and Guestrin, 2016) is an ensemble learning method extensively
employed in various predictive modeling tasks. Fundamentally, XGB combines
multiple weak classifiers, i.e., single decision trees, to construct a robust
predictive model. This approach mirrors other tree-based ensemble methods like
random forest, where the prediction is derived from the summation score of
leaves across all trees. However, XGB distinguishes itself from the random
forest in its tree generation process. While random forest concurrently builds
all decision trees in the ensemble, XGB adopts a sequential approach through
the gradient boosting technique (Figure 1). This methodology entails
constructing new trees based on the loss (see Equation 1) computed from the
previous trees. Additionally, XGB incorporates regularization terms into its
loss function to mitigate overfitting, enhancing its predictive accuracy and
generalization capabilities.

$Loss=\sum_{i=1}^{n}l(y_{i}+y_{i}^{\prime})+\lambda_{1}\sum_{j=1}^{t}\mid
w_{j}\mid+\lambda_{2}\sum_{j=1}^{t}w_{j}^{2}\,,$ (1)

where $l$ denotes the predictive loss function, $y_{i}$ and $y_{i}^{\prime}$
denote the target value and the predicted value of the training sample $i$,
$w_{j}$ denotes the weights from tree $j$, and 1 and 2 denotes the coefficient
for L1 and L2 regularization terms, respectively. Through this iterative
refinement process, where the model continually improves based on the errors
of previous iterations, XGB excels in achieving superior predictive
performance compared to standalone weak classifiers or alternative ensemble
methods. This capability enables XGB to capture intricate patterns and
relationships inherent in the data, rendering it particularly suitable for
tasks where predictive accuracy is paramount.

Figure 2: Developing sequence of an XGB model

###  4.2 SHAP

In this study, we employ the SHapley Additive exPlanations (SHAP; Lundberg and
Lee, 2017) approach to provide insights into the decision-making process of
the XGB model for predicting lateral spreading. The Shapley value (Shapley,
1953), rooted in cooperative game theory, distributes the “contributions” of
players within a game equitably. It achieves this by calculating the marginal
contribution when a player $i$ joins a player subset $S\subseteq F$, where $F$
represents the set of all players. Equation 2 illustrates the computation of
the Shapley value $\phi_{i}$ for player $i$:

$\phi_{i}=\sum_{S\subseteq F}\frac{\mid S\mid!(\mid F\mid-\mid S\mid-1)!}{\mid
F\mid!}[f(S\cup i)-f(S)]\,,$ (2)

where $f$ denotes the model for evaluating the “payout” of a given subset of
players and $[f(S\cup i)-f(S)]$ represents the marginal contribution when a
player $i$ joins the player subset S. The coefficient $\mid S\mid!(\mid
F\mid-\mid S\mid-1)!|F|!$, is applied to the marginal contribution, where
$\mid F\mid$ denotes the total number of players, $\mid F\mid!$ represents the
total permutations of all players, $\mid S\mid!$ accounts for the permutations
of players already in the team before the player $i$ joins, and $(\mid
F\mid-\mid S\mid-1)!$ represents the permutations of players who have not yet
joined the team (excluding player $i$). This coefficient effectively weighs
the impact of the permutation of a subset $S$. The Shapley value of player $i$
is then obtained as the weighted average of marginal contributions across all
possible player subsets.

However, evaluating “payout” from different subsets of input features requires
retraining models, as the inputs vary across subsets. To address this
challenge, SHAP offers an innovative approach by estimating Shapley values,
termed SHAP values, through the conditional expectation of the original model
denoted as $E[f(z)\mid z_{S}]$. Here, $z$ represents the model input with
missing feature values for features not in the subset $S$, representing
unknown information, while $z_{S}$ denotes the feature values for features
within subset $S$. To compute the expected model output with missing features,
we create sample instances with missing features filled by feature values
sampled from the dataset. The random sampled feature values approximate the
effect of missing. The expectation is then computed by averaging the model
outputs of sampled instances. The SHAP value for each feature is calculated as
the change in the expected model prediction when conditioning on that feature.

In Figure 3, we provide an example of SHAP calculation of predicting the risk
of lateral spreading based on three factors: L (Distance to the river), PGA
(Peak Ground Acceleration), and GWD (Groundwater Depth). Consider a model that
predicts the risk of lateral spreading with a risk score ranging from 0 (no
risk) to 100 (highest risk). The model assigns a risk score of 60 for a site
with these characteristics: L = 150 m, PGA = 0.35g, and GWD = 1.5 meters.

Figure 3: An example of SHAP calculation for factors influencing lateral
spreading.

The SHAP analysis begins with a baseline risk score of 50 ($v\\{\\}=50$), here
$v$ represents a function that returns the expected risk, which is the model’s
output without specific data for the factors, based only on the training data.
When assessing the individual impact of each factor, we find that introducing
L (150 m) alone results in a risk score increase to 70, indicating a positive
SHAP value for L ($SHAP_{L}:v\\{L\\}-v\\{\\}=70-50=20$). However, when we
introduce PGA (0.35g) alone, the risk score decreases to 5, suggesting a
negative SHAP value for PGA ($SHAP_{PGA}:v\\{PGA\\}-v\\{\\}=5-50=-45$),
implying that this PGA level reduces the model’s perceived risk. Similarly,
when GWD (1.5 m) is introduced alone, the risk score decreases to 20,
indicating a negative SHAP value for GWD
($SHAP_{GWD}:v\\{GWD\\}-v\\{\\}=20-50=-30$).

Next, we evaluate the pairwise combinations. With both L and GWD, the risk
score increases to 45, showing a combined positive effect, even though GWD
alone has a negative contribution. With L and PGA, the risk score reaches 50,
demonstrating a lower negative contribution from PGA when considered alongside
L, indicating an interaction effect where PGA’s negative impact is mitigated.
With PGA and GWD, the risk score is 40, suggesting that their combined effect
is less negative than their individual contributions.

Considering all three factors, the model predicts a risk score of 60. The
combined influence of L, PGA, and GWD results in a total increase of 10 points
from the baseline, with the interplay between the factors affecting their
individual contributions.

Finally, we average the contributions across all permutations to calculate the
SHAP values. This process yields average SHAP values, with positive values
indicating factors increasing the risk and negative values indicating factors
that reduce the risk score. The SHAP values in order of highest to lowest
importance are $SHAP_{L}=25$ (strongest positive impact), $SHAP_{PGA}=-10$
(negative impact, but less when combined with L), and $SHAP_{GWD}=-5$
(negative impact).

These SHAP values explain the model’s risk prediction for lateral spreading by
illustrating how each factor contributes and how their interactions affect the
overall risk score. The presence of negative SHAP values helps to understand
the conditions under which certain factors may decrease the perceived risk,
providing a more nuanced explanation of the model’s behavior.

One of the key advantages of SHAP is its model-agnostic nature. Regardless of
the model’s underlying structure or type, whether a tree-based, neural
network, or any other functional model, SHAP can compute Shapley values by
analyzing the input permutations and model outputs. Another notable advantage
of SHAP is its ability to provide global and local explanations. While other
explainers, such as permutation feature importance (Breiman, 2001) and Partial
Dependence Plots (PDP; Friedman, 2001) offer only global explanations that
describe the average behavior of the model. SHAP and methods like LIME can
provide explanations for individual predictions, thereby offering local
insights. Furthermore, SHAP explanations exhibit a high degree of stability.
The results remain consistent when rerunning the SHAP explainer due to its
computation of Shapley values. In contrast, LIME explanations may vary as they
involve creating a surrogate model through a sampling process, which can lead
to inconsistencies.

However, one potential drawback of SHAP is its computational time. The
intensive computations required for a single instance, as demonstrated in the
example referred to in Figure 3, can make it impractical to utilize SHAP for a
large number of instances or instances with numerous features. Nevertheless,
given the manageable dataset size in our study, with only 6500 instances and
five to nine features, the computation time remains within a few minutes,
allowing us to select SHAP as the explainer for our analysis confidently.

##  5 Model Development

We develop an XGBoost classifier, Model A, trained on five features (PGA, GWD,
L, slope, and elevation) to predict lateral spreading. The dataset is split
into three subsets with a ratio of 80-5-15 for training, validation, and
testing. Table 1 presents the distribution of Yes/No data across each dataset.
In the training data, the ratio of Yes to No instances is 42-58, while in the
validation data, it is 46-54, and in the testing data, it is 43-57. Although
there is slight variation in the proportion of Yes and No instances across
datasets due to randomness, it closely mirrors the ratio observed in the total
dataset (i.e., 42-58), with no significant imbalance detected.

Model A is trained on training data using the log loss function (Equation 3)
as loss function. The equation for the log loss function is represented as:

$l=-\frac{1}{n}\sum_{i=1}^{n}y_{i}\cdot log(y_{i}^{\prime})+(1-y_{i})\cdot
log(1-y_{i}^{\prime})\,,$ (3)

where $n$ denotes the number of training data points, $y_{i}$ denotes the
binary label indicating lateral spreading for data point $i$ (0 for no lateral
spreading, 1 for lateral spreading), $y_{i}^{\prime}$ represents the predicted
probability, ranging from 0 to 1, generated by Model A.

We optimize the configuration of Model A using cross-validation preprocessing.
Table 2 shows the optimized parameters for Model A. The training also
restricts the depth of the decision trees in the ensemble to avoid over-
fitting the data because our interest is in explaining the predictions rather
than the most accurate model possible.

Table 1: Number of Yes/No data in each dataset. Dataset | Yes lateral spread | No lateral spread
---|---|---
Total | 2752 (42.31%) | 3753 (57.69%)
Train | 2182 (41.93%) | 3022 (58.07%)
Val | 150 (46.01%) | 176 (53.99%)
Test | 420 (43.08%) | 555 (56.92%)
Table 2: Model information: Training set features, model configurations, and performance. |  | Model A | Model B | Model C
---|---|---|---|---
| Distance from the river (L) in km |  |  |
| Peak Ground Acceleration (PGA) in g |  |  |
| Elevation in m |  |  |
| Groundwater Depth (GWD) in m |  |  |
Input Feature | Slope angle in % |  |  | $\bigcirc$
|  $I_{c}$ Median | $\bigcirc$ |  |
|  $q_{c1Ncs}$ Median | $\bigcirc$ |  | $\bigcirc$
|  $I_{c}$ STD | $\bigcirc$ |  |
|  $q_{c1Ncs}$ STD | $\bigcirc$ |  | $\bigcirc$
| Loss function | log-loss
| Max Tree depth | 6
| $\lambda_{1}$ | 0.0
Parameter | $\lambda_{2}$ | 1.0
| Early stop | 50
| Learning rate | 0.300
| Number of estimators | 52 | 29 | 90
| Training (%) | 94.37 | 91.99 | 98.33
Accuracy | Validation (%) | 87.12 | 86.50 | 88.04
| Testing (%) | 84.21 | 82.77 | 85.54

##  6 Results

Model A yields 87.12% accuracy on the validation data and 84.21% on the
testing data (see Table 2). Figure 5() shows the confusion matrix on the
testing data, showing the number of cases in four categories: true positive
(TP), true negative (TN), false positive (FP), and false negative (FN). The
model correctly classifies 488 no-lateral-spreading cases (50.1% TN) and 333
lateral-spreading cases (34.2% TP). The model misclassified 67 no-lateral-
spreading cases (6.9%) as lateral spreading cases (FP) and 87 lateral-
spreading cases (8.9%) as no-lateral-spreading cases (FN). Figure 5() shows
the spatial distribution of these four categories with different markers. We
highlight two regions of false prediction to discuss in the result section:
Region A for FNs and Region B for FPs. In Region A, only 25% of lateral
spreading cases are correctly predicted in the testing data, while the
predictions in the training data match 50% of the cases. Meanwhile, Model A
only predicts 40% of no-lateral-spreading cases in the Region B testing
dataset. This underscores the importance of unpacking the ’black box’ and the
model training dataset to understand the underlying physics or mechanisms
learned by the model.

Figure 6() shows the distribution of predictive probabilities within the
testing dataset. Predictive probability serves as an indicator of the
classification predicted by the model, where a value exceeding 0.5 signifies a
prediction of lateral spreading. In contrast, a value below 0.5 indicates a
classification of no-lateral spreading. Additionally, the predictive
probability is the model’s confidence in its predictions. A probability
closing 1.0 or 0.0 suggests high confidence, whereas a probability near 0.5
indicates that the model is uncertain.

In Figure 6(), the distribution of TPs is depicted in dark gray, peaking at a
probability of 0.99 with 63 instances (19% of TPs). 51.9% of TP cases fall
within the range of 0.9 to 1.0, indicating high confidence in the TP category.
A similar pattern of high confidence is observed in the distribution of TNs.
In contrast, the predictive probabilities associated with FPs and FNs are
closer to 0.5, as shown in Figure 6(). Notably, 34% of FPs fall within the
range of 0.5 to 0.6, while only 10.8% of TPs are within the same range. No FPs
exhibit high predictive probabilities exceeding 0.96. Similarly, 19.5% of FNs
fall within 0.4 to 0.5, in contrast to only 5% of TNs. This comparison
underscores that the model tends to be more confident in its correct
predictions, as evidenced by higher predictive probabilities while displaying
less confidence in incorrect predictions, as indicated by the clustering
around 0.5.

figure5  (a) Confusion matrix of testing data (b) Spatial distribution of
predictions

() Performance of the Model A: (a) Confusion matrix of testing data. (b)
Spatial distribution of predictions of the whole dataset.

figure6  (a) TP and TN (b) FP and FN

() The predictive probability distribution of Model A in the testing set is
categorized as (a) TP and TN (b) FP and FN.

###  6.1 Local explanations

####  6.1.1 True positives and true negatives

We focus on specific sites to demonstrate the predictive ability of XGBoost
Model A and explain the predictions. Consider Site 1138 (see Figure 5()),
which is situated near the river (L = 6 m) and is thus prone to lateral
spreading. Model A correctly predicts lateral spreading with a 93.1%
probability. Figure 7() illustrates the SHAP values as a waterfall plot for
various features influencing this prediction. A positive SHAP value suggests a
higher likelihood of lateral spreading, whereas a negative value indicates no
lateral spreading. Site 1138 has a high positive SHAP value of 2.6.03,
demonstrating a significant probability of lateral spreading. SHAP prioritizes
the five input features based on their predictive importance for this site,
with proximity to the river (L = 6 m) being the most crucial and PGA (0.444 g)
the least. The SHAP explanations suggest that the site’s location near the
river and a low GWD of 1.055 m are the reasons the site will likely experience
lateral spreading. This finding aligns with established domain knowledge.

On the other hand, Site 43457 (see Figure 5()), situated further from the
river, represents a different scenario. The SHAP waterfall plot (see Figure
7()) indicates that its location (L = 2619 m from the river) is a critical
factor in predicting the absence of lateral spreading. The prediction aligns
with our field observations in the training data, where fewer than 17.5% of
sites located more than 2 km from the river experienced liquefaction. The
site’s high elevation, low GWD, and low slope (having a minor positive
influence) contribute to this prediction, resulting in an overall negative
SHAP value of -3.713. Note that the training dataset is slightly unbalanced
relative to the yes/no observations of lateral spreading, with the ratio of no
sites to yes sites, termed the scale positive weight, equal to 1.38. The
imbalance in the dataset leads to an initial bias in the expected SHAP value
($E[f(x)]=v\\{\\}=-0.369$); before considering any features, the model expects
a given site not to experience lateral spreading, i.e., a baseline prediction
is no lateral spreading.

figure7  (a) Site 1138 (True positive); 0.941 m displacement; 93.10% (b) Site
43457 (True negative); 0.214 m displacement; 2.38% (c) Site 32896 (False
positive); 0.259 m displacement; 76.60% (d) Site 56297 (False negative); 0.388
m displacement; 15.86%

() Local explanations of Model A: (a) Site 1138 (True positive); (b) Site
43457 (True negative); (c) Site 32896 (False positive); (d) Site 56297 (False
negative).

####  6.1.2 False positive and false negative

Next, we look into the false positive in Region B. Site 32896 (see Figure
5()), with an observed displacement of 0.259 m and labeled as no lateral
spread (we use a threshold of 0.3 m to label as lateral spreading), is
incorrectly predicted by the XGB model with a 76.60% probability of lateral
spreading occurrence. The SHAP plot (see Figure 7()) shows that L, PGA, and
GWD are the major contributors to the positive total SHAP value of 1.186. The
site is near the river (L = 465m), providing a positive SHAP value of 0.62.
This value is relatively small compared to the counterpart in Figure 7()
because a larger distance decreases the likelihood of lateral spreading. The
SHAP value for the low GWD (1.822 m) is 0.61, while Elevation, PGA, and slope
contribute minor effects to the total SHAP values. Although the prediction is
incorrect, we argue that the label of no lateral spreading for the site is
ambiguous because the observed displacement is close to the defined 0.3 m
lateral spreading threshold. Thus, while the prediction is incorrect, the
behavior of the site and the SHAP values is aligned with our scientific
expectations.

Site 56297 (see Figure 5()), situated further from the river (L = 2104 m) in
Region A, experienced lateral spreading with a displacement of 0.388 m. Figure
7() shows the local SHAP explanation waterfall plot for this site. The XGB
model predicts no lateral spreading with a 15.86% probability of occurrence.
The SHAP value for L is shown as a negative value of -1.6 due to its large
value far from the river. However, we find that the false negatives in Region
A align with the Shirley stream (see Figure 5()), a small stream that was not
considered when computing L. Instead, L was computed based on the distance to
the larger Avon River, which overestimates L, thus reducing its SHAP value.
Recognizing additional water bodies as significant features could enhance the
model’s predictive accuracy. Nonetheless, lateral spreading near smaller
creeks like Dudley Creek and St. Alban Stream are rare, indicating the need
for further research to accurately identify influential factors and reduce
misclassifications due to data inaccuracies. Additionally, for this site (ID:
56297), the PGA has the second most prominent effect, a positive SHAP value of
+0.4, despite a low PGA (0.367 g). In contrast, the SHAP values for larger
PGAs of 0.444 g and 0.473 g (see Figures Figure 7() and Figure 7()) are 0.06
and 0.13, significantly smaller than 0.4. Considering only PGAs, we expect the
SHAP values to increase with increasing PGAs. The difference in SHAP values
with PGAs implies a non-linear and non-monotonic relation between PGA and
other factors. To comprehend the effect of various factors and their
corresponding SHAP values, we assess them globally.

###  6.2 Global explanations

Figure 8() shows a collective plot of SHAP values of all local explanations
from all datasets. Each point represents the SHAP value of a site for a
specific feature. The x-axis refers to SHAP values, and the color denotes the
magnitude of the feature value (red indicates the feature has a high value,
and blue means a low value). The plot highlights how the SHAP value evolves
for each feature, which allows us to explain if the model predictions agree
with domain knowledge. For example, the SHAP values for L are larger than +2
for L smaller than 0.175 km (i.e., blue points), indicating that sites near
the river are predicted to experience lateral spreading. The SHAP values are
smaller than -2 for L greater than 1.71 km (i.e., red points), indicating no
lateral spreading at distances far from the river. However, sites that are at
a moderate distance from the river are scattered close to the center (i.e.,
SHAP 0.0), indicating no clear distinction of lateral spreading for medium L
ranges, using only L.

Similarly, SHAP values shift from positive for shallow GWDs (see Figure 8()),
indicating higher lateral spreading likelihood, to negative for deeper GWDs,
suggesting lower risk. The SHAP trend for GWD aligns with domain knowledge.
Figure 8() shows that neither elevation nor slope angle relates clearly to
SHAP values. Figure 8() shows the global ranking of each feature based on
their mean SHAP value for each site. PGA is the most critical feature among
the five input features responsible for a mean SHAP value of 0.84 in all
predictions. The elevation feature is the third significant, contributing a
mean SHAP value of 0.66 SHAP. In contrast, the slope is the least important
feature, contributing an average of 0.19 SHAP value.

While PGA is the most important feature (see Figure 8() ranking of global
features), it does not have a monotonic relationship with SHAP like L and GWD.
Contrary to our domain knowledge, we observe high PGAs have large negative
SHAP values and a low likelihood of lateral spreading. Also, a few low PGA
values have high positive values and a high likelihood of lateral spreading.
This contribution of PGA to SHAP value is counterintuitive.

figure8  (a) Global explanations (b) Feature importance

() Global explanations and feature importance of Model A. (a) Global
explanations. (b) Feature importance.

###  6.3 Incorrect learning of PGA

Let us explore how Model A incorrectly learned feature-behavior relation from
the dataset using local explanations. We have seen that PGA plays a prominent
role in predictive factors. However, examining a TP site prediction (Figure
9()), we see that a low PGA of 0.376 g positively influences the SHAP value
(+1.29), thus predicting lateral spreading. In contrast, a high PGA of 0.532 g
(Figure 9()) results in a large negative SHAP value (-2.05), thus a low chance
of lateral spreading (22.82%). Even in these TP and TN categories, the model
has incorrectly predicted that high PGAs cause no lateral spreading and that
low PGAs are responsible for lateral spreading. Of course, this incorrect
learning arises from the observations in the training dataset. Inspecting the
training set, we see that 46% of data with PGA in the 0.37 g - 0.375 g
experience lateral spreading (48 sites), resulting in positive SHAP values for
PGA within this range. On the other hand, 83% of training data with high PGA
($>0.52g$) experience no lateral spreading (119 sites), which explains the
counterintuitively negative SHAP value for these sites.

figure9  (a) Site 3781 (True positive) (b) Site 170 (True negative)

() Local waterfall plot explanations for Model A’s incorrect learning of the
feature-behavior relation. (a) Site 3781 (b) Site 170

figure10  (a) Feature relation between Elevation and PGA and its effect on the
SHAP values. (b) Sites with PGA: 0.43 - 0.47 g and Elevation $>$ 3.0 m (c)
Sites with PGA: 0.47 - 0.50 g and Elevation $>$ 4.0 m

() Feature dependency of PGA and Elevation (a) Feature relation between
Elevation and PGA and its effect on the SHAP values. (b) Sites with PGA: 0.43
- 0.47 g and Elevation $>$ 3.0 m (c) Sites with PGA: 0.47 - 0.50 g and
Elevation $>$ 4.0 m

As multiple factors influence the prediction for each site, we plot the impact
of PGA on SHAP with its strongly correlated elevation feature. Figure 10()
shows the relationship between SHAP and PGA values. Contrary to expectations,
the SHAP values for PGA do not monotonically increase with rising PGA values.
Notably, positive SHAP values are concentrated around a PGA of 0.37.
Conversely, there is a distinct cluster of negative SHAP values at higher PGA
levels ($>0.52g$). These patterns corroborate the incorrect model predictions
described earlier.

In addition to the relationship between PGA and its SHAP value, Elevation
influences the SHAP value (as shown by the color gradient in Figure 10()). We
observe a clear separation between high (red) and low (blue) elevation values
from 0.43 g to 0.47 g of PGA. Within that range of PGA, high elevation values
have high SHAP values for PGA, and low elevation values have lower SHAP values
for PGA. Figure 10() shows the spatial distribution of lateral spreading
observations in training data with PGA ranging from 0.43 g to 0.47 g and with
Elevation larger than 3.0 m (corresponding to the purple dots and the red dots
within the yellow ellipse in Figure 10()). We observe that 77% of these sites
experience lateral spreading, while only 44% experience lateral spreading when
the elevation is lower than 3.0 m. An opposite dependency of elevation is
observed in the ranges of 0.47 g to 0.50 g of PGA, where high elevation values
have negative SHAP values and low elevation values have positive SHAP values.
Figure 10() shows the observations with PGA ranging from 0.47 g to 0.50 g and
elevation larger than 4.0 m (corresponding to the red dots within the green
ellipse in Figure 10()). 84% of data experience no lateral spreading, while
only 42% experience no lateral spreading when elevation is below 4.0 m. These
observations agree with the separation observed in Figure 10(), meaning
elevation impacts the SHAP values for PGA.

Figure 10() shows a non-linear relation between PGA and SHAP values. Low PGAs
(0.35 - 0.4 g) show no lateral spreading, while PGAs in the 0.43 - 0.47 g
range show a higher chance of lateral spreading (positive SHAP values), and as
PGA increases above 0.47 g, the sites show no signs of lateral spreading. The
chance of lateral spreading decreasing with PGA is counterintuitive. Consider
the region within the yellow circle in Figure 10(), which shows positive SHAP
values for PGAs between 0.43 - 0.47 g and elevation larger than 3 m. This
cluster of lateral spreading occurrences (sites marked within the yellow
circle in Figure 10()) is shown on a map in Figure 10(). While this cluster
shows signs of lateral spreading based on optical image correlation (Figure
1()), these sites do not exhibit traditional signs of lateral spreading, such
as observation of cracks, and disagree with the observations (Figure 1())
recorded by the New Zealand Geotechnical Database (NZGD, 2013). The
manifestations recorded by NZGD in this area are categorized as “minor
liquefaction with no apparent lateral movement.” This disparity arises from
differences in the methods used to measure lateral movement. NZGD’s
measurements are based on the displacement observed based on cracks, whereas
the lateral spreading occurrences documented in this study rely on optical
image correlation of horizontal displacement. This difference explains the
reasoning behind positive SHAP values for PGA between 0.43 - 0.47 g with
elevation larger than 3.0 m.

Let us explore the reasoning behind no lateral spreading cases for high PGAs
($>0.47g$), marked by a green circle in Figure 10(). The region of interest is
mapped in Figure 10() and experiences significant no lateral spreading.
Although this region experiences higher PGA, it is situated at a much higher
elevation and located further away from the river, thus showing no signs of
lateral spreading.

####  6.3.1 Effect of CPT on Lateral Spreading Predictions

In addition to topological features, the likelihood of lateral spreading also
depends on the site characteristics – i.e., its soil properties. Therefore, we
introduce CPT features to account for the soil properties of the site. This
study uses the same CPT features as Durante and Rathje (2021): the median and
the standard deviation of $q_{c1Ncs}$ and $I_{c}$ over 4 m below the
groundwater table. The parameters $q_{c1Ncs}$ and $I_{c}$ represent the soil
strength and type. A site with high $q_{c1Ncs}$ or $I_{c}$ means the soil is
strong or clayey to liquefy and, consequently, has a low risk of lateral
spreading. The medians and standard deviations represent the profile’s general
trend and variance over 4 m below the GWD. Figure 4 shows the pair plots of
the four CPT features—the plots on the diagonal show the histogram of each
feature. The plots in the lower triangle refer to the scatter plot of each
pair of features. We observe that the median of $q_{c1Ncs}$ and $I_{c}$ have a
strong negative correlation, while other pairs of features do not show strong
correlations. In Figure 4, blue means no lateral spreading, and orange
represents lateral spreading.

Figure 10: Pairplots of CPT features with lateral spreading observations.

We notice a higher frequency of no lateral spreading cases for sites with
fine-grained soil ($I_{c\\_med}$ value exceeds 2.5). Fine-grained soils are
less prone to generating excess pore water pressure during dynamic loading,
making it difficult to liquefy or experience lateral spreading. Similarly,
when $q_{c1Ncs}$ med surpasses 175, we also observe a greater prevalence of no
lateral spreading cases. High $q_{c1Ncs}$ med values indicate dense soil,
which offers higher shear strength against dynamic loading, reducing the
likelihood of lateral spreading. Despite the expectation that regions with low
$q_{c1Ncs}$ med ($<75$), i.e., loose soil, are more susceptible to
liquefaction, we find a clustering of cases without lateral spreading in this
region. However, most low $q_{c1Ncs}$ med sites in our study region have high
$I_{c\\_med}$ values (see red circle in Figure 4), indicating that the loose
soils are fine-grained and less susceptible to liquefaction. Hence it is
important to correlate CPT features to lateral spreading to improve our model
predictions.

We develop another XGBoost model, Model B, using the same dataset as Model A
with four additional CPT features as input parameters (see Table 2). After
optimization, all hyperparameters are the same, except Model B uses only 26
estimators, as the validation performance does not improve with more
estimators. Model B yields an accuracy of 86.5% on the validation dataset and
82.77% on the testing dataset. Table 2 compares the performance of Models A
and B. Adding CPT features did not improve the performance of Model B compared
to Model A. This result agrees with Durante and Rathje (2021). We will dive
into local and global explanations of Model B to evaluate why adding CPT
features did not improve the model’s performance.

#### Global Explanation of Model B:

To evaluate the global performance of Model B, we summarize all the local
explanations of Model B predictions in Figure 12(). The trend of SHAP values
for Model B aligns with Model A. Both GWD and L show behavior that aligns with
our expectations, and the SHAP values of PGA show a counterintuitive trend.
However, the CPT features do not show a clear trend. Even after adding site-
specific soil properties through CPT, Figure 12() shows that PGA remains the
most influential factor in Model B, with a mean SHAP value of 0.68. GWD and
Elevation exhibit similar contributions to predictions, each with a mean SHAP
value of 0.54. L becomes less critical in Model B than in Model A, ranking
fourth. On the other hand, including all four CPT features only make minor
contributions to predictions.

figure12  (a) Global explanation (b) Feature importance

() Global explanations and feature importance of Model B (a) Global
explanation (b) Feature importance

#### Local Explanation of Model B:

Let us explore the local explanations of Model B to understand the impact of
CPT features. Site 5252 (see Figure 5()), situated near the river (L = 99 m)
with a moderate PGA of 0.42 g and a deep GWD of 2.569 m, experiences lateral
spreading. Model A accurately predicts this occurrence with a 68.46%
probability. The SHAP waterfall plot for Model A (Figure 13()) highlights the
proximity to the river, contributing a decisive SHAP value of 1.414, leading
to a true positive prediction. However, including CPT features results in an
incorrect prediction of no-lateral-spreading by Model B (Figure 13()). Despite
L contributing a positive SHAP value of 0.59, the median values of $I_{c}$
(2.546) and $q_{c1Ncs}$ (69.031) counteract its effect with SHAP values of
-0.8 and -0.55, respectively. High $I_{c}$ and low $q_{c1Ncs}$ indicate that
clayey soil is less prone to liquefaction and lateral spreading. While the
negative SHAP values align with our understanding of clayey soil, they result
in a false negative prediction. In contrast, Site 17049 (see Figure 5())
benefits from considering the CPT features. Model A misclassifies the site as
prone to lateral spreading, but Model B accurately predicts lateral spreading
by including CPT data. The site is located 1 km from the Avon River with a GWD
of 2 m, and Model A incorrectly predicted a 61% chance of lateral spreading
due to positive SHAP values in both L and GWD of 0.52 and 0.31 (see Figure
13()). However, Model B (see Figure 13()), incorporating CPT features like the
low $I_{c}$ (1.319) and high $q_{c1Ncs}$ (199.992), yields negative SHAP
values of -0.56 and -0.35, indicating strong sandy soil resistant to
liquefaction, leading to a correct prediction of no lateral spreading.

figure13  (a) Site 5252 in Model A (True positive) (b) Site 5252 in Model B
(False negative) (c) Site 17049 in Model A (False positive) (d) Site 17049 in
Model B (True negative)

() Local explanation of Site 5252 and Site 17049 (a) Site 5252 in Model A
(True positive); (b) Site 5252 in Model B (False negative); (c) Site 17049 in
Model A (False positive); (d) Site 17049 in Model B (True negative).

#### What Model B has learned from CPT features:

While incorporating CPT features does not notably enhance predictive
performance, understanding Model B’s interpretation of these features is
interesting. Figure 14() presents the SHAP values of features plotted in
median $I_{c}$ and median $q_{c1Ncs}$ space, where the color indicates the
sign and magnitude of the SHAP value.

Figure 14() illustrates the SHAP value of median $I_{c}$. A distinct
separation between positive (red) and negative (blue) SHAP values is
noticeable at $I_{c}=1.5$ and $I_{c}=2.4$. When $I_{c}$ is less than 1.5,
indicating coarse-grained/sandy soil, the SHAP value ranges from 0 to -0.5,
while for $I_{c}$ greater than 2.4, indicating fine-grained soil, the SHAP
value ranges from -0.5 to -1.0. These negative SHAP values represent no
lateral spreading, which is expected for fine-grained soil but not sandy soil.
However, these small $I_{c}$ values are associated with large $q_{c1Ncs}$, so
the soil does not liquefy, and lateral spreading does not occur. In the
intermediate range ($1.5<I_{c}<2.4$), the SHAP value ranges from -0.1 to 0.3.
For these silty/sandy soils, the potential for lateral spreading is not
directly associated with $I_{c}$, but rather with the relative density of the
soil as indicated by $q_{c1Ncs}$. A similar pattern appears in the SHAP values
for median $q_{c1Ncs}$ (Figure 14()). When the median $q_{c1Ncs}$ falls
between 160 and 210, the SHAP value ranges from -0.5 to 0, indicating no
lateral spreading due to the dense nature of the soil. For $q_{c1Ncs}$ values
below 75, the SHAP value ranges from 0 to -1.0, indicating no lateral
spreading due to the large $I_{c}$ associated with these small values of
$q_{c1Ncs}$. For $q_{c1Ncs}$ between 75 and 160, the SHAP value is
approximately 0 to 0.2.

figure14  (a) SHAP values for $I_{c}$ median (b) SHAP values for $q_{c1Ncs}$
median (c) Combined SHAP values

() SHAP values for medians of Ic and qc1Ncs in Model B (a) SHAP values for
$I_{c}$ median; (b) SHAP values for $q_{c1Ncs}$ median; (c) Combined SHAP
values of $I_{c}$ median and $q_{c1Ncs}$ median.

Comparing Figures 14() and 14(), we observe a similar trend in SHAP values for
median $I_{c}$ and $q_{c1Ncs}$. Low median $I_{c}$ data often corresponds to
high median $q_{c1Ncs}$, and vice versa, resulting in a similar trend of SHAP
values. Figure 14() combines the SHAP values for median $I_{c}$ and
$q_{c1Ncs}$, revealing regions of positive (red) and negative (blue) SHAP
values consistent with our understanding of liquefaction and soil properties,
as explained earlier in this section.

#### Improving model performance with CPT data:

Although Model B effectively learns the relationship between SHAP values for
$I_{c}$ and $q_{c1Ncs}$, their incorporation does not improve predictive
performance. Additionally, Figure 12() indicates that the CPT features
($q_{c1Ncs}$ ) are least important for the model predictions. Hence, we
develop a new Model C by excluding the three least important features
identified in Figure 12(): the median and standard deviation of $q_{c1Ncs}$,
and slope. Model C exhibits improved performance, achieving accuracy rates of
88.04% and 85.54% on the validation and testing data, which are higher than
Model B (testing accuracy is 82.77%). Since $q_{c1Ncs}$ and $I_{c}$ are
negatively correlated, including only the $I_{c}$ factors improves the model’s
performance.

figure15  (a) Global explanation (b) Feature importance (c) Site 5252 in Model
C (True Positive) (d) Site 17049 in Model C (True negative)

() Global and local explanation of Model C (a) Global explanation; (b) Feature
importance; (c) Site 5252 in Model C (True positive); (d) Site 17049 in Model
C (True negative).  Figure 15: SHAP values for medians of $I_{c}$ in Model C

In Figures 15() and 15(), we present the global explanation for Model C. The
results show that the most important feature is PGA, followed by L, elevation,
and GWD. Similar to the ranking in Figure 12(), the CPT features, the median
and standard deviation of $I_{c}$ are the least important. To further examine
Model C’s behavior, we look at the local explanations for Site 5252 and Site
17049 (see Figures 15() and 15()).

For Site 5252, the vicinity to the river (L = 99m) significantly impacts the
prediction, with a SHAP value of 1.4. However, the median $I_{c}$ value of
2.546 has a negative contribution with a SHAP value of -0.69. Despite this,
the effect of proximity to the river outweighs the $I_{c}$ feature, leading to
a true positive prediction, which was incorrectly classified as a false
negative by Model C.

In the case of Site 17049, Model C has learned that a low $I_{c}$ implies a
high $q_{c1Ncs}$ and thus, low likelihood of lateral spreading. This is
reflected in the decisive SHAP value of -1.8 for the $I_{c}$ feature,
resulting in a true negative prediction. This suggests that the high
correlation between median $I_{c}$ and median $q_{c1Ncs}$ allows Model C to
learn the underlying physics, as shown in Figure 11, using $I_{c}$ only.

We aimed to avoid feature coherence and overfitting by minimizing the number
of input features, developing simpler and more explainable models. However,
this approach may make Model C less generalizable to other sites, as its
performance relies on the negative correlation between $q_{c1Ncs}$ and $I_{c}$
observed in the training dataset.

##  7 Conclusion

Our investigation into using SHapley Additive exPlanations (SHAP) for
interpretability in machine learning models targeting lateral spreading
predictions highlights four major insights. First, SHAP’s application
clarifies the XGB model’s decision-making process, enhancing model
transparency and fostering informed decision-making. Second, we identified
incorrect learning patterns within the model, especially regarding PGA’s
impact, underscoring the need to align model predictions with domain
knowledge. Third, incorporating CPT data into the XGB model, although did not
change the predictive accuracy, improved model explanations. Notably, the
model identified negative SHAP values for firm sandy and clayey soil, aligning
with our understanding of soil types and liquefaction behavior. The XGB model
accurately learned the relationship between soil characteristics inferred from
CPT data and the likelihood of lateral spreading events conforming to the
domain knowledge. Lastly, the study emphasizes the need for model optimization
to enhance efficiency and generalizability. This research underlines the
critical role of explainable AI in improving the reliability and utility of ML
models for geotechnical engineering, contributing to more accurate hazard
assessments and infrastructure protection.

## Conflict of Interest Statement

The authors declare that the research was conducted without any commercial or
financial relationships that could be construed as a potential conflict of
interest.

## Author Contributions

CH: Writing – original draft, Methodology, Formal analysis; KK: Writing –
review & editing, Conceptualization, Supervision, Funding acquisition; ER:
Writing – review & editing, Supervision.

## Acknowledgments

The authors thank the support of NHERI DesignSafe-CI for computing resources
and data storage Rathje et al. (2017a).

## Data Availability Statement

The SHAP code is available under the MIT license on GitHub
(https://github.com/geoelements/xai-lateral-spreading). The training dataset
and trained models are published under a CC-By license on DesignSafe Data
Depot (Durante and Rathje, 2022).

## References

  * Arrieta et al. (2020) Arrieta, A. B., Díaz-Rodríguez, N., Del Ser, J., Bennetot, A., Tabik, S., Barbado, A., et al. (2020).  Explainable Artificial Intelligence (XAI): Concepts, taxonomies, opportunities and challenges toward responsible AI.  _Information Fusion_ 58, 82–115.  10.1016/j.inffus.2019.12.012
  * Bach et al. (2015) Bach, S., Binder, A., Montavon, G., Klauschen, F., Müller, K. R., and Samek, W. (2015).  On Pixel-Wise explanations for Non-Linear Classifier decisions by Layer-Wise relevance propagation.  _PLOS ONE_ 10, e0130140.  10.1371/journal.pone.0130140
  * Breiman (2001) Breiman, L. (2001).  Random Forests.  _Machine Learning_ 45, 5–32.  10.1023/a:1010933404324
  * Chen and Guestrin (2016) Chen, T. and Guestrin, C. (2016).  XGBoost: A scalable tree boosting system.  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ (ACM), 785–794.  10.1145/2939672.2939785
  * Cubrinovski et al. (2011) Cubrinovski, M., Green, R. A., and Wotherspoon, L. (2011).  Geotechnical reconnaissance of the 2011 christchurch, new zealand earthquake.  _Geotechnical Extreme Events Reconnaissance Association_ 10.18118/G68G65
  * Demir and Şahin (2022) Demir, S. and Şahin, E. K. (2022).  Liquefaction prediction with robust machine learning algorithms (SVM, RF, and XGBoost) supported by genetic algorithm-based feature selection and parameter optimization from the perspective of data processing.  _Environmental Earth Sciences_ 81\.  10.1007/s12665-022-10578-4
  * Doshi‐Velez and Kim (2017) Doshi‐Velez, F. and Kim, B. (2017).  Towards a rigorous science of interpretable machine learning.  _arXiv (Cornell University)_ 10.48550/arxiv.1702.08608
  * Durante and Rathje (2022) [Dataset] Durante, M. G. and Rathje, E. (2022).  Machine learning models for the evaluation of the lateral spreading hazard in the avon river area following the 2011 christchurch earthquake.  10.17603/DS2-3ZDJ-4937
  * Durante and Rathje (2021) Durante, M. G. and Rathje, E. M. (2021).  An exploration of the use of machine learning to predict lateral spreading.  _Earthquake Spectra_ 37, 2288–2314.  10.1177/87552930211004613
  * Friedman (2001) Friedman, J. H. (2001).  Greedy function approximation: A gradient boosting machine.  _The Annals of Statistics_ 29, 1189 – 1232.  10.1214/aos/1013203451
  * Kaya (2016) Kaya, Z. (2016).  Predicting Liquefaction-Induced lateral spreading by using neural network and Neuro-Fuzzy techniques.  _International Journal of Geomechanics_ 16\.  10.1061/(asce)gm.1943-5622.0000607
  * Linardatos et al. (2020) Linardatos, P., Papastefanopoulos, V., and Kotsiantis, S. (2020).  Explainable AI: A review of Machine Learning Interpretability Methods.  _Entropy_ 23, 18.  10.3390/e23010018
  * Lundberg and Lee (2017) Lundberg, S. and Lee, S.-I. (2017).  A unified approach to interpreting model predictions.  _neural information processing systems_ 30, 4768–4777
  * Miller (2019) Miller, T. (2019).  Explanation in artificial intelligence: Insights from the social sciences.  _Artificial Intelligence_ 267, 1–38.  10.1016/j.artint.2018.07.007
  * New Zealand Geotechnical Database (2013) [Dataset] New Zealand Geotechnical Database (2013).  Liquefaction and lateral spreading observations. Map layer CGD0300 – 11 February 2013.  Accessed 5 May 2022
  * Paris et al. (2004) Paris, G., Robilliard, D., and Fonlupt, C. (2004).  _Exploring overfitting in genetic programming_ (Springer Berlin Heidelberg).  10.1007/978-3-540-24621-3_22
  * Rathje et al. (2017a) Rathje, E. M., Dawson, C., Padgett, J. E., Pinelli, J., Stanzione, D., Adair, A., et al. (2017a).  DesignSafe: New cyberinfrastructure for Natural Hazards engineering.  _Natural Hazards Review_ 18\.  10.1061/(asce)nh.1527-6996.0000246
  * Rathje et al. (2017b) Rathje, E. M., Secara, S. S., Martin, J. G., Van Ballegooy, S., and Russell, J. K. (2017b).  Liquefaction-Induced Horizontal Displacements from the Canterbury Earthquake Sequence in New Zealand Measured from Remote Sensing Techniques.  _Earthquake Spectra_ 33, 1475–1494.  10.1193/080816eqs127m
  * Ribeiro et al. (2016) Ribeiro, M. T., Singh, S., and Guestrin, C. (2016).  ”Why should I trust you?”: explaining the predictions of any classifier.  _arXiv (Cornell University)_
  * Shapley (1953) Shapley, L. S. (1953).  _17\. A value for N-Person games_ (Princeton University Press).  10.1515/9781400881970-018
  * Xue (2019) Xue, Y. (2019).  An Overview of Overfitting and its Solutions.  _Journal of Physics: Conference Series_ 1168, 022022.  10.1088/1742-6596/1168/2/022022