# Novel Building Detection and Location Intelligence in Aerial Satellite
Imagery

Sandeep Singh
<EMAIL_ADDRESS>
Christian Wiles
<EMAIL_ADDRESS>
Ahmed Bilal
<EMAIL_ADDRESS>

###### Abstract

Building structures detection and information about these buildings in aerial
images is an important solution for city planning and management, land use
analysis. It can be the center-piece to answer important questions such as
planning evacuation routes in case of an earthquake, flood management, etc.
These applications rely on being able to accurately retrieve up-to-date
information. Being able to accurately detect buildings in a bounding box
centered on a specific latitude-longitude value can help greatly. The key
challenge is to be able to detect buildings which can be commercial,
industrial, hut settlements, or skyscrapers. Once we are able to detect such
buildings, our goal will be to cluster and categorize similar types of
buildings together.

##  1 Introduction

We plan on reproducing semantic segmentation CNN models based on U-Net
Ronneberger et al. (2015) and Res-U-net Diakogiannis et al. (2020) algorithms,
trained by transfer learning using the ImageNet dataset. In addition, we will
optimize these different models with focal loss, dice loss, cross entropy
loss, hierarchical loss and differently weighted intersection-over-union (IoU)
loss to overcome issues of scale difference [3] in building detection. To
clearly delineate each individual member’s contribution, we will be organizing
our paper by each member’s contribution.

##  2 Methodology

We have decided to work on three different models in parallel, which are
inherently different on following parameters:

  1. 1.

Architecture of models: We have decided to use different variations of UNet
architectures, which might help us focus of different patterns easily.

  2. 2.

Pre-training usage of the models: We have employed training from scratch as
well using pre-trained weight for the encoder layers across out 3 of the
models.

  3. 3.

Loss functions being employed by models: We have used different loss functions
in all the different models.

  4. 4.

Image sampling and augmentations employed: We have used totally different
techniques to samples the images from the data sets available to introduce the
component of randomness in data distribution across models. While model 1 and
2 have used random resized crops of size 224x224 from big 5000x5000 images and
masks; model 3 has resorted to fixed size split of big 5000x5000 images in
512x512 tiles. Both these sampling techniques have again used totally
different set of transforms with different quantifiers.

While we have details of all the models provided in the sections 3,4 and 5,
here is a quick overview of the ensemble technique exploited by us:

  1. 1.

Get prediction from all 3 models. All predictions must of same size as input
image. We output thresholded mask as tensor of softmax probabilities.

  2. 2.

Chose the most confident pixel from each of these masks from 3 models for
every pixel location’s softmax probability. This is will be single merged mask
of same size as input.

  3. 3.

From merged mask predicted, drop anything that is not at least 0.75 confident.
Rational behind doing this is that we have already chose most confident pixel
locations from 3 models masks, So, all the softmax value at this stage must be
pretty confident ones at least in one of the models.

As seen in the fig 1. We have final mask, which is without any less confident
pixel.

Figure 1: Ensemble method being employed. Output softmax probabilities are
compared and chosen to select most confident pixels of 2 models and
thresholded to be at least 75% confident to be considered for final mask.

##  3 Standard U-Net Architecture

As all three models are derived from U-net, we will begin with a discussion of
the canonical U-net architecture. U-net is a modular model largely consisted
of encoder blocks and decoder blocks. A graphical depiction of a u-net model
can be found in Figure 5.

###  3.1 Encoder blocks

The encoder layer is primarily responsible for detecting the ’what’ elements
of the images. The goal is to be able to extract features in the image at
different scales and different levels of abstraction. As such, at every steps
of the encoder, two 2D convolutional blocks are used to extract information
from the image and double the size of the feature space. At each encoder
layer, we used a maxpool layers of 2x2 kernel size and a stride of 2 for down
sampling spatially. This allowed us to increase the number of filters at each
of our encoder layers without being extremely computationally expensive and
increase the receptive field of our filters with deeper layers allowing for
segment detection at multiple scales.

###  3.2 Decoder blocks

The decoder layers are the up sampling layers in the model. The primary
purpose of these layers is to localize the features extracted in the encoder
block. This information is essential in our semantic segmentation in order to
be able to output an image with the buildings detected localized in the right
spaces in our output mask. For up sampling, we used Transposed Convolution
layers. This allows us to ultimately assign class labels to each pixel in our
image as part of our semantic segmentation.
At each of our decoder layers, we also make use of skipped connection given to
us by the respective encoder layer for our decoder layer. The skipped
connections cross from same sized part in the encoders to the decoders. The
skipped connections allow us to overcome problem of vanishing gradient,
increasing dimensionality and help regain the initial spatial information that
we lost during the encoding path.

###  3.3 Integration

A full u-net model is composed of N encoder blocks, and N-1 decoder blocks.
The feature space of the first encoder block is a hyperparameter but seems to
be often set to 64 or 128. Save for the last encoder block, the output of the
final convolutional layer in each encoder block is cropped and concatenated
with the output of the transposed convolutional layer of the decoder block. At
the end of the model, a 1x1 convolutional layer is used to create a classifier
head with the same features space as the number of classes. This can be passed
to a softmax layer to produce class probabilities.

##  4 Model 1 approach (Christian)

###  4.1 Data Pipeline and Exploration

For this project, we used the Inria Aerial Image Labelling Dataset for
training Maggiori et al. (2017). The dataset consists of 360 (180 train and
180 test) 5000x5000 pixel full-color images with corresponding masks
indicating the presence of building or non-building pixels. Only the training
set had ground-truth segmentation masks available. Images were taken from a
variety of settings, including rural and urban cities from different
continents. A few problems had to be solved to enable training with this
dataset: data augmentation and quick random access.

####  4.1.1 Data pre-processing and augmentation

To generate more data for training, data augmentation was undertaken.
Pytorch’s standard transformation library does not make allowances for
maintaining consistent transformations between an image and a segmentation
mask, so the functional transforms library was used, which allows for the
randomness to be provided by external variables, which can be held static
between the ground truth and full-color images. For each of the 180 input
training images, 350 224x224x3 image patches were created. This was intended
to allow transfer learning for networks trained on ImageNet. A patch was taken
from the image with random width (between 100 and 500 px), height (within +/-
10% of the random width), and image origin. This image was then randomly
flipped horizontally and/or vertically and normalized by ImageNet standard
deviation and mean. The intention of this was to teach the model scale and
orientation invariant features. Once the incoming data was processed, it was
split 80/20 into train and validation sets. Data from all 5 cities in the
dataset was randomly selected for both validation and training set, as the
test set consists of different cities and would be usable for testing how well
the model can generalize. A manual seed was set such that this split was
repeatable if restarting training from a checkpoint.

####  4.1.2 Caching

It was found that performing the loading of the full-color images and
performing transformations on the fly was too computationally intensive. To
alleviate this problem, after the first time the transformations were done,
the resulting input and target tensors were saved to disk. This reduced the
time to train significantly, as a 70 MB image did not have to be loaded and
manipulated thousands of times per epoch, but instead a 600 KB tensor could be
used.

####  4.1.3 Other Considerations

It was also deemed important to add support for visualization of training-
related metrics. Tensorboard support was added to the project to track
training and validation set loss and accuracy, precision-recall curves for the
validation set, and visualization of the forward pass of the model on the
validation set. One last consideration made was the use of a seed when
splitting the training and validation set. It was noticed that when resuming
training from a checkpoint that the validation set was not the same as before
the checkpoint. By maintaining the constant seed, a barrier was maintained
between the two sets.

####  4.1.4 Dataset Statistics

As is often the case with segmentation tasks, the dataset was not balanced
between building pixels and non-building pixels. The training dataset was
analyzed to determine the prevelance of each class. The findings are below.

###  4.2 U-net from scratch in Pytorch

To test the hypothesis that building detection was a sufficiently specific
domain to merit training from scratch, a u-net was created with random Xavier-
initialized weights. 4 encoder and 3 decoder blocks were used, with the first
encoder block having a hidden dimension of 64 features. Both the original
U-net paper Ronneberger et al. (2015) and Johannes Schmidt’s blog posts
Schmidt (2021) were consulted in the creation of the model.

Two major deviations were attempted from the models mentioned above. Both
u-nets resulted in a cropped image with every convolution due to the use of
valid padding. By using same padding, on convolutions and transposed
convolutions, we can return an image that is of equal size to the input. This
may result in slightly worse accuracy in the extremities of the image due to
the extrapolation employed by same padding, but does simplify some aspects of
the analysis, as every mask pixel has a corresponding prediction.

Secondly as the output of a 2-class softmax classifier only has 1 degree of
freedom, it was attempted to perform classification as a single-class
regression, with the output of the regression put through a sigmoid function.
This one-channel output can then be interpreted as p(building). This approach
was eventually discarded, as it had a very small impact on model size due to
only affecting the final 1x1 convolutional layer, and adding a second
classifier dimension increased model performance by a few percentage points.

Finally, batch normalization was added between the convolutional and
activation layers in encoder and decoder blocks. These recenter the
distribution of the output of the convolutional layers and add to stability in
training as seen in Santurkar et al. (2019).

To discourage overfitting, dropout was added between the output of the final
decoder layer and the 1x1 convolution.

####  4.2.1 Loss function

3 different loss funcitons were attempted with this model. With the
regression-based approach, weighted mean square error was used due to the
class imbalance. On a per-batch basis, the effective number of building and
non-building pixels was calculated on a per-batch basis, similar to the
methodology in Cui et al. (2019). This weighting was weigh the loss on the
minority class (buildings) more heavily.

Once the model had progressed to a two-class method, two losses were pursued:
weighted cross-entropy loss and dice loss. Dice loss was pursued due to its
background in segmentation tasks and invariance with respect to class
imbalance (as dice loss is related to the size of the true positive region).
Binary cross-entropy weighted by the inverse of effective number was also
explored. This provides a more convex loss function that should be easier and
more stable to train.

####  4.2.2 Training

Training was performed over 20 epochs with the Adam optimizer. No learning
rate scheduler was used to govern learning rate as epochs progressed, as Adam
should manage its own learning rates on a per-parameter basis Kingma and Ba
(2017). A batch size of 20 was used to fit in GPU memory. The model state was
saved to disk whenever validation accuracy exceeded the previous maximum to
allow for training to be resumed later.

####  4.2.3 Results

Loss Fn | Accuracy | IOU Score | F1-Score
---|---|---|---
Dice | 94.9% | 0.717 | 0.836
BCE | 95.0% | 0.726 | 0.841
Table 1: Model 1 Validation Set results after 20 Epochs Figure 2: Training
accuracy over 20 epochs for model 1. Weighted CE in orange, Dice Loss in red.
Figure 3: Validation accuracy over 20 epochs for model 1. Weighted CE in
orange, Dice Loss in red.

Weighted cross-entropy resulted in marginally better training efficiency and
overall metrics, but by an almost negligible amount. For both approaches,
overfitting does not appear to be a concern, as the validation and training
accuracy are almost identical.

Though this model performed worse than Model 2, it is hard to say whether it
is due to differences in pre-training or the squeeze-and-attention layers. It
is likely, however, that due to learning features from scratch, it may provide
diversity in the ensemble that can help in overall accuracy.

##  5 Model 2 approach (Ahmed)

###  5.1 Model Specifications

####  5.1.1 Double Convolution Blocks

The Unet build consisted of a double convolution layer, where each convolution
layer consisted of a kernel size of 3, stride and padding of 1. We set the
bias to false in order to add a BatchNorm layer, which is then followed by a
ReLU activation layer. We settled on a small 3X3 kernel receptive field in our
convolution layers in order to be able to detect very small edges and shapes
in our aerial images. Doing so is especially relevant for our aerial images as
there is a lot of noise in the images and our model needs to be able to use
small edges and shapes to detect buildings as buildings appear in many
different sizes in our input images.

####  5.1.2 Encoder Layers

The authors in U-Net: Convolutional Networks for Biomedical Image Segmentation
Ronneberger et al. (2015) recommend encoding layers with output channels
64,128,256, 512 and 1024. However, we found more success with output channel
layers 16,32,64,128 and 256. We believe this is because lower output channels
of 16 and 32 in the start allow us to detect really small building segments
with a small receptive field. In addition, the 512 and 1024 channel layers
were not leading to any significant performance gains in our testing.

####  5.1.3 Decoder Layers

At each of our decoder layers, we also make use of skipped connection given to
us by the respective encoder layer for our decoder layer. The skipped
connections cross from same sized part in the encoders to the decoders. The
skipped connections allow us to overcome problem of vanishing gradient,
increasing dimensionality and help regain the initial spatial information that
we lost during the encoding path.

###  5.2 Pre-trained ResNet34 Encoder Specifications

We now add ResNet 34 Encoder layers to the model. As such, we are now
performing the Double Convolution blocks 3,4,6, and 3 times at each encoder
layer level, using skipped connections between encoder layers, and using a
higher stride to down sample instead of max pooling. These encoder layers are
also pre-trained on image-net dataset.

####  5.2.1 Encoder Modifications

After finding success in our U-Net built with 16,32 and 64 output channel
initial encoder layers, we replace the initial ResNet convolution, ReLU and
Max Pool layers with our U-Net 16,32 and 64 output channel encoder layers with
skipped connections in order to preserve a lot of the small shapes and edges
information in our images.

####  5.2.2 Attention Mechanism

For the loss function, we will be using the Dice Loss to create cleaner mask
segments to represent the buildings. In order to supplement our model in
reducing the Dice Loss, we also include an attention mechanism using spatial
and channel ’squeeze & excitation’ Blocks. This is done to aid our encoder
layers in spatial encoding for more accurate mask prediction and better
network flow. The authors in Recalibrating Fully Convolutional Networks with
Spatial and Channel ’Squeeze & Excitation’ Blocks Roy et al. (2018) found a
reduction of 4-9% in the Dice Loss. We see similar results in our testing.

####  5.2.3 Results

In our testing, we saw the pre-trained image-net backbone significantly
increase the model performance. After 15 epochs, we saw the following results.

Set | Acc | Loss | IOU Score | F1-Score
---|---|---|---|---
Train | 0.965 | 0.116 | 0.749 | 0.856
Val | 0.961 | 0.129 | 0.702 | 0.824
Table 2: Validation Set results after 15 Epochs

Detailed Model 2 results and mask outputs available in appendix Section F

##  6 Model 3 Distinctive Approach (Sandeep)

###  6.1 Data Sampling Strategy

For this model, we have used “Progressive Resizing.” This is a training
technique where we purposefully change the contents of image by resizing the
images to contain more area. Instead of randomized crops of size 224x224, we
created 512x512 non-overlapping and contiguous tiles, only then resizing them
to 224x224 input images. On average, each tile has almost 4 times more
buildings in each tile compared to model 1 and 2. Hence, the model more easily
learns smaller buildings in more crowded areas Howard (2018).

###  6.2 Architectural Considerations

For model 3, we have evaluated 3 types of encoders (Resnet18/34/50) and choose
ResNet34 as ResNet18 has shown to be struggling to encode features of smaller
buildings successfully. ResNet34 and 50 have shown very similar, but resnet 50
slow performance in detected buildings without any marginal increase in
performance. One more significant improvement in model 3 was the use of we
used Pixel Shuffle up-sampling in the decoder blocks, as provided by
shuffleblock implementation Aitken et al. (2017).

###  6.3 Training and Validation Split

For model 3, we have done the data split on the basis of geography, instead of
random ratio split. Out of training data from 5 cities as: Austin, Chicago,
Kitsap, Tyrol, Vienna. Different cities are included in each of the subsets.
e.g., images over Chicago are included in the training set, but not on the
test set. Also, images over San Francisco are included on the test set but not
on the training set. At the same time, we have tried to include the training
data all type of structures of building. e.g. low rise vs high rise vs
community living buildings apartment complexes.

###  6.4 Custom Loss Function Design

Model 3 did not used CrossEntropy loss. Instead, we have written our own
custom loss function, which has helped us predict foreground pixel with higher
softmax confidence. We implemented Combined Loss of Dice Loss and Focal Loss
with equal weights. Focal loss penalizing more confident wrong predictions
more heavily. We have used gamma value of focal loss as 2. Also, Dice score
has provided feed back to strive to keep precision and recall both highest
possible. Also, for model 3, we have used the Dice Score metric.

###  6.5 Training Convergence

For model 3, we used a LR finder scheduler before starting to fine tune the
pretrained weights of ResNet34 encoder. This has helped of find the most
appropriate maximum LR value. The second innovative technique employed by us
was ”Fit-one-cycle” Smith (2018) to achieve super convergence. In this
technique, we increase LR to maximum value in initial batches before start to
anneal the learning rates. Please refer to figure for LR finder and fit one
cycle both Sylvain Gugger (2018). This technique is taken from Leslie Smith
iconic Super Convergence paper. Also, We trained with total 40 epochs without
over-fitting and saved the model only when better score on validation was seen
without over-fitting, while Dice score was approaching 92%.

Train_Loss | Valid_Loss | Dice_Score
---|---|---
0.102414 | 0.115400 | 0.920870
Table 3:  Model 3 Losses and Metrics Values Figure 4: Model 3 Training Curves.

###  6.6 Post Processing Enhancements

As can be seen in some our validation images, the masks created by this model
are not always crisp and polygonal. Please refer to appendix figure for Model
3’s traditional computer vision techniques of segmentation and attempt to fuse
the results with Unet predicted mask. We have tried to experiment with Otsu’s
threshold, Watershed segmentation, SLIC Super pixel algorithm for
segmentation. We have selected all the coincidental mask segments from super
pixel algorithms’ output with Unet’s mask and tried to shape correct polygon
for buildings with sharp edges (see fig. 8).

##  7 Experiences and Challenges

###  7.1 Challenges

Satellite images are very noisy and affected by many factors such as weather,
zoom level, resolution, trees, and cost to obtain Wikipedia (2018). After
evaluating, we used Google Static maps API for pulling additional data because
of quality and ease of use.

###  7.2 Project Success Criterion

We were able to build model with Dice score more than 90 percent on test set
and model 2 achieving detection of more 96 percent ground truth pixels in
valid set. With help of Softmax based adaptive selection and ensemble, we have
achieved detection of more 93 percent ground truth pixels in test set.

###  7.3 Conclusion and Future Aspirations

We have explored and confirmed that deep learning based segmentation is very
effective in segmenting buildings. We could extend these models for clustering
similar buildings, classifying residential vs commercial, predicting future
constructions or detecting illegal construction or activities. Potential
technical improvements include learned polygonization, elastic transformations
during training, and model compression for cheaper inference.

## References

  * Aitken et al. (2017) Andrew Aitken, Christian Ledig, Lucas Theis, Jose Caballero, Zehan Wang, and Wenzhe Shi. 2017.  Checkerboard artifact free sub-pixel convolution: A note on sub-pixel convolution, resize convolution and convolution resize.
  * Cui et al. (2019) Yin Cui, Menglin Jia, Tsung-Yi Lin, Yang Song, and Serge Belongie. 2019.  Class-balanced loss based on effective number of samples.
  * Diakogiannis et al. (2020) Foivos I. Diakogiannis, François Waldner, Peter Caccetta, and Chen Wu. 2020.  Resunet-a: A deep learning framework for semantic segmentation of remotely sensed data.  _ISPRS Journal of Photogrammetry and Remote Sensing_ , 162:94–114.
  * Howard (2018) Jeremy Howard. 2018.  Progressive resizing for better generalization of the computer vision models.
  * Kingma and Ba (2017) Diederik P. Kingma and Jimmy Ba. 2017.  Adam: A method for stochastic optimization.
  * Maggiori et al. (2017) Emmanuel Maggiori, Yuliya Tarabalka, Guillaume Charpiat, and Pierre Alliez. 2017\.  Can semantic labeling methods generalize to any city? the inria aerial image labeling benchmark.  In _IEEE International Geoscience and Remote Sensing Symposium (IGARSS)_. IEEE.
  * Ronneberger et al. (2015) Olaf Ronneberger, Philipp Fischer, and Thomas Brox. 2015.  U-net: Convolutional networks for biomedical image segmentation.
  * Roy et al. (2018) Abhijit Guha Roy, Nassir Navab, and Christian Wachinger. 2018.  Recalibrating fully convolutional networks with spatial and channel ’squeeze & excitation’ blocks.  _CoRR_ , abs/1808.08127.
  * Santurkar et al. (2019) Shibani Santurkar, Dimitris Tsipras, Andrew Ilyas, and Aleksander Madry. 2019.  How does batch normalization help optimization?
  * Schmidt (2021) Johannes Schmidt. 2021.  Creating and training a u-net model with pytorch for 2d & 3d semantic segmentation: Dataset…
  * Smith (2018) Leslie N. Smith. 2018.  A disciplined approach to neural network hyper-parameters: Part 1 – learning rate,batch size, momentum, and weight decay.
  * Sylvain Gugger (2018) Jeremy Howard Sylvain Gugger. 2018.  Adamw and super-convergence is now the fastest way to train neural nets.
  * Wikipedia (2018) Wikipedia. 2018.  Satellite imagery and challenges associated with them.

##  Appendix A Code Repository

We have work on total of 4 repositories during our project life cycle. 3
repositories where used by each of us individually and 1 final repository was
created as a place to perform ensemble and integration of all models together
and do all the needed post-processing. Here the repositories as below:

  * •

Integrated Final repo:
https://github.com/sandeepsign/building_footprint_ensemble

  * •

Ahmed’s repo:
https://github.com/abilal19/DL_FinalProject_Draft

  * •

Christian’s repo:
https://github.com/cswksu/aerialDetection

  * •

Sandeep’s repo:
https://github.com/sandeepsign/building_footprints_cs7643

##  Appendix B Individual Contributions

Contributor | Contribution
---|---
Ahmed | Full Model 2 (best performing),
| Project Report
Christian | Full Model 1, data preprocessing
| and visualization for models
| 1 and 2, Project Report
Sandeep | Full Model 3, ensembling,
| post processing, Geo Coding, VM Setup,
| Clustering of polygons, Project Report

Table 4: Individual team member contributions.

##  Appendix C Data Source

We have used INRIA’s spacenet challenge data from:
https://project.inria.fr/aerialimagelabeling/

##  Appendix D Train Infrastructure

We have used combination of techniques to execute this project.

Each had used individual hardware for setup and eventually to run long time
training on google cloud VM instance.

GPUs Used are:
nVIDIA QUADRO RTX 5000 16GB
nVIDIA T4 16GB
nVIDIA GTX 1080ti 11GB

##  Appendix E Enlarged Figures

Figure 5: Example of overall organization of u-net model from Ronneberger et
al. (2015)

Class | Number of Pixels | Percentage | Effective Number
---|---|---|---
|  |  | ($\beta$ = 1 - 10E-9)
Building | 7.1E8 | 1.58% | 5.08E8
Not building | 4.4E10 | 98.4% | 1.00E9
Total | 4.5E10 | 100% | 1.00E9

Table 5: Class distribution of pixels in training set images.

##  Appendix F Model 2 Results

####  F.0.1 Model 2 Validation Set Mask Samples

Validation Set Target
Validation Set Output

####  F.0.2 Accuracy Results Per Epoch

Figure 6: ’Model 2 Train Accuracy’ Figure 7: ’Model 2 Valid Accuracy’

##  Appendix G Model 3 figures

Figure 8: Post processing and merging of traditional segmentation techniques
with Unet’s mask and shape correcting polygons for buildings. Figure 9: Model
3 learning rate finder Figure 10: Model 3 fit one cycle scheduler: Momentum is
moved opposite to learning rate during increase or annealing phase


# Decomposition Based Refinement for the Network Interdiction Problem

Krish Matta  School of Computer Science
Carnegie Mellon University
Pittsburgh, USA
<EMAIL_ADDRESS>Xiaoyuan Liu  Fujitsu Research of America
Sunnyvale, USA
<EMAIL_ADDRESS>Ilya Safro  Department of Computer and Information Sciences
University of Delaware
Newark, USA
<EMAIL_ADDRESS>
###### Abstract

The shortest path network interdiction (SPNI) problem poses significant
computational challenges due to its NP-hardness. Current solutions, primarily
based on integer programming methods, are inefficient for large-scale
instances. In this paper, we introduce a novel hybrid algorithm that can
utilize Ising Processing Units (IPUs) alongside classical solvers. This
approach decomposes the problem into manageable sub-problems, which are then
offloaded to the slow but high-quality classical solvers or IPU. Results are
subsequently recombined to form a global solution. Our method demonstrates
comparable quality to existing whole problem solvers while reducing
computational time for large-scale instances. Furthermore, our approach is
amenable to parallelization, allowing for simultaneous processing of
decomposed sub-problems.
Reproducibility: Our source code and experimental results are available at
https://github.com/krishxmatta/network-interdiction.

###### Index Terms:

decomposition, Ising processing hardware, network interdiction, refinement,
integer programming

##  I Introduction

Network interdiction refers to a class of challenging combinatorial
optimization problems that involve strategically disrupting flows in a network
to achieve specific objectives [30]. These problems can be characterized as a
game between attackers and defenders of a network. The defender seeks to
optimize a predefined objective across the network, such as maximizing flow
between two nodes, while the attacker aims to impede the defender’s objective
by inflicting maximum disruption upon the network. This disruption may
manifest in a variety of forms, such as targeted attacks on the network arcs
aiming to destroy or impair them.

The topic of network interdiction has gained significant attention due to its
ability to model complex situations in the defense domain. For example, Ghare
et al. [9] demonstrated how network interdiction can be used in wartime to
capacitate an enemy’s supply network to maximally disrupt the flow of enemy
troops. Network interdiction can additionally be used to identify weaknesses
in critical infrastructure to make them more resilient to terrorist attacks
and natural disasters. An example of such an application can be found in a
study by Salmerón et al. [27], which details the use of network interdiction
to model attacks on large-scale power grids. The study then utilizes this
model to determine how best to minimize economic loss resulting from these
attacks. Another study by Nandi et al. [23] details the usage of a network
interdiction model on cyber attack graphs to pinpoint vulnerabilities in
cyberinfrastructures, thereby helping prevent organizations from cyberattacks.
Such problems are even extremely useful outside of military applications, e.g.
in this line of work [19, 2, 11] the authors show how network interdiction can
be used to minimize the generally defined infection spread and attacks.

Within this paper, we study the Shortest Path Network Interdiction (SPNI)
problem. In this problem, the defender wishes to traverse the minimum length
path between two specified nodes $s$ and $t$ in a directed network. The
attacker uses their limited resources to destroy certain arcs, or increase
their effective length, to increase the defender’s shortest $s$-$t$ path
length. Each arc is either interdicted or not, and additionally, it is known
beforehand how much an arc costs to interdict and how much interdiction of an
arc causes its effective length to increase. When solving SPNI, we take the
viewpoint of the attacker, and thus our goal is to maximize the shortest
$s$-$t$ path length.

Our contribution The SPNI is an $\mathcal{NP}$-Hard problem [16]. To the best
of our knowledge, almost all current SPNI algorithms rely on integer
programming methods. On large problem instances, these methods are time
costly—for example, in one of the most important works regarding SPNI, the
largest network solved by Israeli et al. had 240 nodes and 1,042 arcs [16], so
these methods are either not scalable or produce large gaps or low quality
results. To ameliorate this issue, we propose a novel algorithm which may
leverage Ising processing units (IPUs) [5], specialized computational devices
specifically tailored to solve the Ising model, as well as classical solvers
that can exhibit high quality solutions for small instances in a reasonable
computational time. Examples of IPUs include but are not limited to quantum
annealers, gate-based quantum machines equipped with the Ising model solvers,
and both digital and analog annealers [22, 20]. Since physical limitations of
existing IPU hardware severely restrict the size of problems they can handle
and connectivity between variables, we adopt a hybrid IPU-classical approach
in which a classical computer decomposes SPNI into subproblems small enough
for an IPU, offloads them for computation, then combines the IPU’s results
into a global solution. For proof of concept, we use exact solvers instead of
a real IPU. Our decomposition based results demonstrate an ability to achieve
almost the same quality as the whole problem slow solvers that are prohibitive
for large-scale instances and IPU hardware. Additionally, our approach is
parallelization friendly, namely, sub-problems obtained as a result of SPNI
decomposition can be tackled in parallel.

##  II Related Work

In one of the earliest works studying SPNI, Malik et al. focused on a variant
termed the $k$ most vital arcs problem in which the interdiction of each arc
requires exactly one unit of the attacker’s budget and interdiction of an arc
results in its complete removal [21]. Ball et al. show that the $k$ most vital
arcs problem is $\mathcal{NP}$-Hard. Since this is a variant of SPNI, it
follows that SPNI is $\mathcal{NP}$-Hard as well [3]. Corley and Shaw [6]
investigated the single-most-vital-arc problem, which is the $k$ most vital
arcs problem where $k=1$, but this problem is a simple variant of SPNI that
lies in $\mathcal{P}$. Rather than study the full removal of arcs, Fulkerson
and Harding [8] and Golden [12] study related problems in which the length of
an arc increases as more budget is allocated to its interdiction. In one of
the most prominent works studying SNPI, Israeli and Wood [16] study the
problem in its more general form, proposing two integer-programming based
algorithms of different quality depending on whether arcs are fully removed or
if arc length is increased as a result of interdiction.

More recent work on SPNI has since shifted away from integer-programming
algorithms. Huang et al., for example, approximate SPNI using a reinforcement
learning framework [15]. However, this approach generates quite significant
gaps and the reinforcement learning is not particularly scalable. Rocco et al.
[25] study an extension of SPNI where the attacker is not solely focused on
maximizing the $s$-$t$ path length, but also wishes to minimize interdiction
cost, and propose an evolutionary algorithm to do so. This is a multiobjective
problem that is different than ours. For more details, we refer the reader to
the recent comprehensive survey on the network interdiction problem [30]. To
our knowledge, however, no attempts have been made for parallelization
friendly decomposition of SPNI for IPU hardware.

Several frameworks have been proposed to solve large instances of
combinatorial optimization problems on small IPU hardware. The quadratic
unconstrained binary optimization (QUBO) formulation is a popular formulation
for combinatorial optimization problems that is equivalent to the Ising model,
and thus is solvable on all IPU hardware [10]. QUBO formulations contain only
quadratic binary variables with no constraints. A large number of work in the
quantum computing space focuses on formulating problems into one large QUBO,
then identifying sub-QUBOs that can be solved directly on an IPU. The qbsolv
tool developed by D-Wave accomplishes exactly this, once allowing researchers
to solve a 1254 binary variable problem using quantum hardware that can only
solve problems with at most 64 binary variables [24]. An alternative approach
is, instead of first formulating the problem as a QUBO then decomposing the
QUBO itself, to decompose the problem and then formulate these subproblems
into QUBOs. For example, Shaydulin et al. [29] take this approach to solve the
community detection problem on quantum hardware. This is the approach we opt
to take in this paper, applied to SPNI.

##  III Problem Formulation

Let $G=(N,A)$ be a directed network with node set $N$ and arc set $A$. The
length of arc $k\in A$ is $c_{k}\geq 0$, and the added interdiction length is
$d_{k}\geq 0$, meaning that if arc $k$ is interdicted, then its effective
length becomes $c_{k}+d_{k}$. If interdiction destroys arc $k$, then $d_{k}$
can be set to a sufficiently large value. For any node $i\in N$, let
$A^{+}(i)\subseteq A$ denote the set of all arcs directed out of node $i$, and
similarly let $A^{-}(i)\subseteq A$ denote the set of all arcs directed into
node $i$. While other forms of SNPI allow for multiple resource types and
variant interdiction costs, for our algorithmic specializations we model a
single resource type and the scenario in which interdiction of any arc in $A$
costs only one unit of this resource. Thus, let $r_{0}$ denote the amount of
budget we have available. Let $\mathbf{x}=\\{x_{k}\\}_{k=1}^{|A|}$ where for
each arc $k\in A$, $x_{k}$ is a binary variable denoting whether arc $k$ is
interdicted or not by the attacker, and let
$\mathbf{y}=\\{y_{k}\\}_{k=1}^{|A|}$ where $y_{k}$ is a binary variable
denoting whether arc $k$ is traversed by the defender. Finally, let $s,t\in N$
denote the source and sink node in $G$, respectively. The integer programming
formulation of SNPI is

$\begin{array}[]{ll@{}ll}\displaystyle\max_{\mathbf{x}\in
X}\min_{\mathbf{y}}&\displaystyle\sum_{k\in A}(c_{k}+x_{k}d_{k})y_{k}\\\
\text{s.t.}&\displaystyle\sum_{k\in A^{+}(s)}y_{k}-\sum_{k\in
A^{-}(s)}y_{k}=1\\\ &\displaystyle\sum_{k\in A^{+}(t)}y_{k}-\sum_{k\in
A^{-}(t)}y_{k}=-1\\\ &\displaystyle\sum_{k\in A^{+}(i)}y_{k}-\sum_{k\in
A^{-}(i)}y_{k}=0,\forall i\in N\setminus\\{s,t\\},\\\ \end{array}$

where
$X=\\{\mathbf{x}\in\\{0,1\\}^{|A|}~{}\mid~{}\sum_{i=1}^{|A|}\mathbf{x}_{i}\leq
r_{0}\\}$.

Note that the above formulation is biobjective, and thus can not be directly
converted into a QUBO model. As Israeli et al. [16] show, we may fix
$\mathbf{x}$, take the dual of the inner minimization, make some
modifications, then release $\mathbf{x}$, and obtain a single objective
formulation

$\begin{array}[]{ll@{}ll}\displaystyle\max_{\mathbf{x},\mathbf{\pi}}&\pi_{t}\\\
\displaystyle\text{s.t.}&\pi_{j}-\pi_{t}-d_{k}x_{k}\leq c_{k},\forall
k=(i,j)\in A\\\ \displaystyle&\pi_{s}=0\\\ \displaystyle&\mathbf{x}\in
X.\end{array}$ (1)

We may interpret $\pi_{i}$ as the post-interdiction shortest-path length from
$s$ to $i$. As such, we may impose bounds on
$\pi_{i}\in[0,|N|\max_{k^{\prime}\in A}\\{c_{k^{\prime}}+d_{k^{\prime}}\\}]$.
Since each variable is constrained and we have a single objective, we may now
convert each constraint to a penalty function and decompose bounded variables
into several binary variables, for example using the mapping proposed by
Karimi et al. [17]. Letting $P$ be a sufficiently large positive penalty
value, our QUBO formulation is

$\begin{array}[]{ll@{}ll}\displaystyle\max_{\mathbf{x},\mathbf{\pi},\mathbf{m},n}&\displaystyle\pi_{t}-P\sum_{k=(i,j)\in
A}(c_{k}-(m_{k}+\pi_{j}-\pi_{i}-d_{k}x_{k}))^{2}\\\
&\displaystyle-P(\pi_{s})^{2}-P(r_{0}-(n+\sum_{k\in A}x_{k}))^{2}\\\
\text{s.t.}&\displaystyle\pi_{i}\in[0,|N|\max_{k^{\prime}\in
A}\\{c_{k^{\prime}}+d_{k^{\prime}}\\}],\forall i\in N\\\ &\displaystyle
m_{k}\in[0,|N|\max_{k^{\prime}\in
A}\\{c_{k^{\prime}}+d_{k^{\prime}}\\}+d_{k}],\forall k\in A\\\ &\displaystyle
n\in[0,r_{0}]\\\ &\displaystyle\mathbf{x}\in\\{0,1\\}^{|A|}.\end{array}$ (2)

Note that we have included bounded non-binary variables for simplicity.
Additionally, we have included variables $\mathbf{m}$ and $n$ to act as slack
variables for inequality constraints. In the actual QUBO formulation, all
bounded variables and slack variables should be mapped to several binary
variables. The given QUBO formulation is directly solvable on an IPU.

We now describe how to formulate subproblems of SPNI for the purposes of our
algorithm. We define a connected partition of a network $G=(N,A)$ as a subset
of nodes $N$ such that its induced subgraph is weakly connected. Say we have a
connected partition of $G$, i.e., a subset of nodes $N_{p}\subseteq N$, a sink
node $t^{\prime}\in N_{p}$, and let $A_{p}$ represent all arcs associated with
this subset of nodes, i.e. $A_{p}=\\{(u,v)\in A|v\in N_{p}\\}$. Note that we
may partition $A_{p}=A_{p,1}\cup A_{p,2}$ where $A_{p,1}=\\{(u,v)\in A|u\in
N_{p}\land v\in N_{p}\\}$ and $A_{p,2}=\\{(u,v)\in A|u\notin N_{p}\land v\in
N_{p}\\}$. We define a subproblem of SPNI over this partition $N_{p}$ and sink
node $t^{\prime}$ by restricting our free variables to only those associated
with $N_{p}$—rather than define $\pi_{i}$ for all $i\in N$, we will only
define $\pi_{i}$ for all $i\in N_{p}$, and similarly only define $x_{k}$ for
all $k\in A_{p}$. Additionally, we change the objective formula to maximize
$\pi_{t^{\prime}}$ rather than $\pi_{t}$.

For those variables not associated with $N_{p}$, we may fix their values to
constants. Assume that we have a current solution $A^{\prime}\subseteq A$—by
“solution,” we refer to a subset of arcs chosen for interdiction. We use
$\Gamma_{i}$ to denote a constant representing the shortest-path length from
node $s$ to $i$ after all arcs in $A^{\prime}$ have been interdicted. Since
$\pi_{i}$ represents the postinterdiction shortest-path length from nodes $s$
to $i$, we can fix $\pi_{i}$ to the value $\Gamma_{i}$ for all nodes $i$ not
in $N_{p}$. Additionally, we limit our budget from $r_{0}$ to
$r_{0}-|A^{\prime}\setminus A_{p}|$, ensuring we use only the budget we have
allocated within $N_{p}$. The subproblem formulation over $N_{p}$ with sink
node $t^{\prime}$ is thus

$\begin{array}[]{ll@{}ll}\displaystyle\max_{\mathbf{x},\mathbf{\pi},\mathbf{m},n}&\displaystyle\pi_{t^{\prime}}-P\sum_{k=(i,j)\in
A_{p,1}}(c_{k}-(m_{k}+\pi_{j}-\pi_{i}-d_{k}x_{k}))^{2}\\\
&\displaystyle-P\sum_{k=(i,j)\in
A_{p,2}}(c_{k}-(m_{k}+\pi_{j}-\Gamma_{i}-d_{k}x_{k}))^{2}\\\
&\displaystyle-P(\pi_{s})^{2}\\\ &\displaystyle-P((r_{0}-|A^{\prime}\setminus
A_{p}|)-(n+\sum_{k\in A_{p}}x_{k}))^{2}\\\
\text{s.t.}&\displaystyle\pi_{i}\in[0,|N|\max_{k^{\prime}\in
A}\\{c_{k^{\prime}}+d_{k^{\prime}}\\}],\forall i\in N_{p}\\\ &\displaystyle
m_{k}\in[0,|N|\max_{k^{\prime}\in
A}\\{c_{k^{\prime}}+d_{k^{\prime}}\\}+d_{k}],\forall k\in A_{p}\\\
&\displaystyle n\in[0,r_{0}].\end{array}$ (3)

Note that in the above formulation, we have restricted our free variables to
those associated with $N_{p}$, thereby decreasing computational load and
potentially allowing the problem to be solved on an IPU depending on the size
of $N_{p}$.

##  IV Algorithm

In this section we present an algorithm for approximating SPNI by leveraging
small IPU hardware. Our algorithm first proceeds by generating an initial
greedy solution, then refining this solution over several iterations.

First we define terminology and helper functions. A problem instance of SPNI
is a tuple $\mathcal{P}=(G,\mathbf{c},\mathbf{d},s,t,r_{0})$ where $G=(N,A)$
is a directed network, $\mathbf{c}$ is a vector representing arc lengths of
$A$, $\mathbf{d}$ is a vector representing added interdiction arc lengths of
$A$, $s\in N$ is the source node, $t$ is the sink node, and $r_{0}$ is the
interdiction budget.

We define PARTITION to be a function that takes in a graph $G=(N,A)$ and an
integer $n$ as inputs, then returns $P\subseteq\mathcal{P}(V)$, a connected
partitioning of $G$ where $\forall N_{p}\in P$, $|N_{p}|\approx n$. These
partitions are how we will choose to define subproblems. The purpose of the
size restriction is to ensure that subproblems are sufficiently small for IPU
hardware. Additionally, PARTITION is necessarily non-deterministic ensuring
different partitions on each call to prevent the solver from being stuck in a
local attraction basin.

We define CALC-LENGTH, a function that takes in a problem instance
$\mathcal{P}=(G=(N,A),\mathbf{c},\mathbf{d},s,t,r_{0})$, and a set
$A^{\prime}\subseteq A$ where $|A^{\prime}|\leq r_{0}$. It returns the
shortest $s$-$t$ path length in $G$ after all arcs in $A^{\prime}$ have been
interdicted.

We also define CALC-PATH that takes in a problem instance
$\mathcal{P}=(G=(N,A),\mathbf{c},\mathbf{d},s,t,r_{0})$, and a set
$A^{\prime}\subseteq A$ where $|A^{\prime}|\leq r_{0}$. It returns the set
$N^{\prime}\subseteq N$ consisting of all nodes in the shortest $s$-$t$ path
in $G$ after all arcs in $A^{\prime}$ have been interdicted.

Finally, function SOLVE takes in a problem instance
$\mathcal{P}=(G=(N,A),\mathbf{c},\mathbf{d},s,t,r_{0})$, a set $N_{p}\subseteq
N$ representing a continuous partition of $G$, a node $t^{\prime}\in N_{p}$
for the sink node of the subproblem, and a set of arcs $A^{\prime}\subseteq A$
where $|A^{\prime}|\leq r_{0}$ representing the current working solution. It
formulates a subproblem based off these parameters according to (3), solves
the subproblem, and returns a new arc set $A^{\prime\prime}\subseteq A$
representing which arcs to interdict to maximize $\pi_{t^{\prime}}$. Here is
where computational work may be offloaded to an IPU. Note that due to the
subproblem restriction, SOLVE can only modify the interdiction status of arcs
associated with $N_{p}$, and thus all arcs not associated with $N_{p}$ but in
$A^{\prime}$ will remain in $A^{\prime\prime}$.

Initial solution generation The pseudocode is given in Algorithm 1. It
receives in input a problem instance
$\mathcal{P}=(G=(N,A),\mathbf{c},\mathbf{d},s,t,r_{0})$ alongside an integer
$n$ representing the approximate nodes allowed per partition. To generate an
initial solution to SPNI, the algorithm takes a greedy approach. As each arc
costs exactly one unit of the budget to interdict, the algorithm iterates over
each unit of the budget attempting to decide which arc to allocate this unit
to. To make this decision, the algorithm first partitions the network using
PARTITION and additionally calculates the shortest $s$-$t$ path with CALC-
PATH. Note that it is sufficient to only consider edges along the $s$-$t$ path
for interdiction, since edges outside of this path will not affect the $s$-$t$
path length after interdiction. For each node in the $s$-$t$ path, the
algorithm determines which partition this node belongs to, then solves a
subpartition using this node as the sink node by calling SOLVE. It then takes
the output solution, and calls CALC-LENGTH to determine this solution’s impact
on the $s$-$t$ path length. If this path length is greater than the current
maximum path length, the algorithm notes it as a potential candidate solution
and updates the maximum. After repeating this process for each node in the
$s$-$t$ path, the algorithm chooses a candidate solution randomly, then
continues this process until the budget has ran out. Note that in each
iteration, the algorithm incrementally adds one arc to the current solution.

Solution refinement The pseudocode is given in Algorithm 2. It receives in
input a problem instance
$\mathcal{P}=(G=(N,A),\mathbf{c},\mathbf{d},s,t,r_{0})$, an integer $n$
representing the approximate number of nodes per partition, an integer
$\lambda$ for the number of iterations to run the algorithm, and a solution
$A^{\prime}\subseteq A$ to improve. Note that the given solution $A^{\prime}$
is expected to be the output of Algorithm 1. This algorithm attempts to
improve the initial solution. To do so, it iterates $\lambda$ times, and in
each iteration, partitions the graph using PARTITION and calculates nodes in
the current $s$-$t$ path given the current interdiction solution using CALC-
PATH. For each of these nodes, alongside all nodes found in the shortest path
for the previous iteration, it determines what partition the node lies in and
solves a subproblem where this node is the target node using SOLVE. CALC-
LENGTH is then called to determine whether the outputted solution’s $s$-$t$
path length is greater than the current $s$-$t$ path length, and if so, marks
this solution as a potential candidate solution. After iterating over each
node, if a better $s$-$t$ path length is found, the current solution is chosen
to be a randomly selected candidate solution. If a better $s$-$t$ path length
is not found, the algorithm attempts to explore alternate solutions by
iterating over each arc in the current solution, temporarily deleting that arc
from the solution, and observing whether this deletion can result in a better
solution through a partitioning process similar to the above. If a better
solution is found, the deletion of this arc becomes permanent and the better
solution is adopted. Note that to prevent excessive computation, the algorithm
keeps track of arcs which, upon being temporarily deleted, have not previously
increased the $s$-$t$ path length in a set named good-arcs. In future
iterations, the algorithm does not consider temporarily deleting arcs found in
good-arcs.

Require: A problem instance
$\mathcal{P}=(G=(N,A),\mathbf{c},\mathbf{d},s,t,r_{0})$

Require: Approximate number of nodes per partition $n$

1 $A^{\prime}\leftarrow\\{\\}$

2prev-nodes $\leftarrow$ CALC-PATH$(\mathcal{P},A^{\prime})$

3$r^{\prime}\longleftarrow 1$

4while  _$r^{\prime}\leq r_{0}$_ do

5 $P\leftarrow\text{PARTITION}(G,n)$

6 curr-nodes $\leftarrow$ CALC-PATH$(\mathcal{P},A^{\prime})$

7 next-sols $\leftarrow\\{A^{\prime}\\}$

8 best-length $\leftarrow$ CALC-LENGTH$(\mathcal{P},A^{\prime})$

9 foreach  _node $\in$ curr-nodes $\cup$ prev-nodes_ do

10 $N_{p}\leftarrow\emptyset$

11 foreach  _$N_{p}^{\prime}\in P$_ do

12 if  _node $\in N_{p}^{\prime}$_ then

13 $N_{p}\leftarrow N_{p}^{\prime}$

14 end if

15

16 end foreach

17

18 curr-sol $\leftarrow$ SOLVE$(\mathcal{P},N_{p},\text{node},A^{\prime})$

19 curr-length $\leftarrow$ CALC-LENGTH$(\mathcal{P},\text{curr-sol})$

20 if  _curr-length $>$ best-length_ then

21 next-sols $\leftarrow\emptyset$

22 end if

23 if  _curr-length $==$ best-length_ then

24 next-sols $\leftarrow$ next-sols $\cup\\{\text{curr-sol}\\}$

25 end if

26

27 end foreach

28

29 $A^{\prime}\stackrel{{\scriptstyle R}}{{\leftarrow}}$ next-sols
$\triangleright$ Pick a solution randomly

30 $r^{\prime}\leftarrow r^{\prime}+1$

31 end while

return $A^{\prime}$

Algorithm 1 Initial solution generation

Require: A problem instance
$\mathcal{P}=(G=(N,A),\mathbf{c},\mathbf{d},s,t,r_{0})$

Require: Approximate number of nodes per partition $n$

Require: Number of iterations $\lambda$

Require: Current solution $A^{\prime}$

1 prev-nodes $\leftarrow$ CALC-PATH$(\mathcal{P},A^{\prime})$

2next-sols $\leftarrow\\{A^{\prime}\\}$

3best-length $\leftarrow$ CALC-LENGTH$(\mathcal{P},A^{\prime})$

4good-arcs $\leftarrow\emptyset$

5$i\leftarrow 0$

6while  _$i <\lambda$_ do

7 $P\leftarrow$ PARTITION$(G,n)$

8 curr-nodes $\leftarrow$ CALC-PATH$(\mathcal{P},A^{\prime})$

9 foreach  _node $\in$ curr-nodes $\cup$ prev-nodes_ do

10 $N_{p}\leftarrow\emptyset$

11 foreach  _$N_{p}^{\prime}\in P$_ do

12 if  _node $\in N_{p}^{\prime}$_ then

13 $N_{p}\leftarrow N_{p}^{\prime}$

14 end if

15

16 end foreach

17

18 curr-sol $\leftarrow$ SOLVE$(\mathcal{P},N_{p},\text{node},A^{\prime})$

19 curr-length $\leftarrow$ CALC-LENGTH$(\mathcal{P},\text{curr-sol})$

20 if  _curr-length $>$ best-length_ then

21 next-sols $\leftarrow\emptyset$

22 end if

23 if  _curr-length $==$ best-length_ then

24 next-sols $\leftarrow$ next-sols $\cup\\{\text{curr-sol}\\}$

25 end if

26

27 end foreach

28 prev-length $\leftarrow$ CALC-LENGTH$(\mathcal{P},A^{\prime})$

29 if  _best-length $>$ prev-length_ then

30 $A^{\prime}\stackrel{{\scriptstyle R}}{{\leftarrow}}$ next-sols
$\triangleright$ Pick a solution randomly

31 else

32 foreach  _arc $\in A^{\prime}\setminus$ good-arcs_ do

33 $A^{\prime\prime}\leftarrow A^{\prime}\setminus\\{\text{arc}\\}$

34 foreach  _node $\in$ curr-nodes $\cup$ prev-nodes_ do

35 partition = $\emptyset$

36 foreach  _$p\in P$_ do

37 if  _node $\in p$_ then

38 partition $\leftarrow p$

39 end if

40

41 end foreach

42

43 curr-sol $\leftarrow$
SOLVE$(\mathcal{P},N_{p},\text{node},A^{\prime\prime})$

44 curr-length $\leftarrow$ CALC-LENGTH$(\mathcal{P},\text{curr-sol})$

45 if  _curr-length $>$ prev-length_ then

46 $A^{\prime}\leftarrow$ curr-sol

47 else

48 good-arcs $\leftarrow$ good-arcs $\cup\\{$ arc $\\}$

49 end if

50

51 end foreach

52

53 end foreach

54

55

56 end while

return A’

Algorithm 2 Solution refinement

##  V Computational Results

To evaluate the performance of our algorithm and similar to other network
interdiction literature, we have generated grid networks with randomized edge
weights and randomized interdiction weights. To each grid we added the source
node $s$ connected to all nodes in the first column, and the sink node $t$
connected to all nodes in the last column, with all of these arcs having
length $0$ and interdiction length $0$ to make the problem more challenging.
All horizontal arcs are oriented forward and all vertical arcs are oriented
down. All arcs not connected to $s$ or $t$ have a randomly assigned length and
interdiction length, each in the range $[1,10]$. Figure 1 depicts an example
of grid 3x3.

Figure 1: Example of a generated 3x3 grid

Models were built using Pyomo, and were solved using the open-source solver
CBC [7]. No significant difference was observed in experiments with Gurobi
solver [13]. Graph representation and utilities such as Dijkstra’s was done
using NetworkX [14]. Graph partitioning was done using METIS [18]. All
problems were ran on the University of Delaware’s Caviness cluster which uses
a Linux system and has exclusively Intel E5-2695 v4 18 core processors in its
pool.

Below we plot the results of our computational experiments. The “Size” axis
represents the number of nodes in the experiment’s network ($|N|$). Note that
by the structure of the grid network, the number of arcs can be calculated by
$|A|=2(|N|-2)$. The “Quality” axis displays a metric to evaluate the quality
of our solutions (defined in Eq. (4)).

We run both our refinement algorithm and CBC attempting to solve Equation (1)
across the entire network—we shall hereby refer to the latter as the “full
problem solver.” Note that the full problem solver is significantly slower
than the refinement algorithm, and thus, we timeout the full problem solver
once it has exceeded the time that the refinement algorithm takes, then
proceed to take the current best solution. _In every problem instance below,
the full problem solver always timed out_. For any given experiment, if $r$ is
the shortest $s$-$t$ path length after interdicting the arcs from the solution
of our refinement algorithm and $f$ the shortest $s$-$t$ path length after
interdicting the arcs from the full problem solver’s solution, the ratio can
be calculated by

$(r-f)/\max(r,f).$ (4)

Thus, the higher the quality, the better—positive quality indicates that the
refinement solver performed better than the full problem solver.

We run three sets of experiments. In each experiment, we set the number of
refinement iterations $\lambda$ to 50. (However, the convergence is often
observed earlier.) In the first set of experiments, represented by Figure 2,
we set the interdiction budget $r_{0}$ to $0.25\%$ of the number of arcs in
the network and limit our subproblem size to $n=20$ nodes. The average quality
for this set of experiments is $-0.008$, implying that a subproblem size of
only 20 nodes is sufficient for a good estimation of SPNI. In the second set
of experiments, represented by Figure 3, we set the interdiction budget
$r_{0}$ to $0.5\%$ of the number of arcs in the network and limit our
subproblem size to $n=20$ nodes. The average quality for this set of
experiments is $-0.013$, indicating relative consistency in estimation
regardless of budget. In the final set of experiments, we set the interdiction
budget $r_{0}$ to $0.25\%$ of the number of arcs in the network and limit our
subproblem size to $n=40$ nodes. The average quality for this set of
experiments is $-0.012$, implying that _an increase in subproblem size is not
necessary for estimation and limited capacity IPU machines can already be
integrated to solve this problem using decomposition-based refinement_.

Out of all 39 experiments ran, the difference between the $s$-$t$ path length
produced by the refinement solver and the path length produced by the full
problem solver was $\geq-10$ for 35 of the experiments. In 10 of these
experiments, the refinement solver performed better than the full problem
solver. _These results are significant as they indicate that in the cases when
full problem solver was better the refinement missed the value of a single
arc._

Figure 2: Computational results for different size networks with parameters:
$\lambda=50$, $n=20$, $r_{0}=0.0025|A|$. Each point corresponds to a graph.
The horizontal axis corresponds to the number of nodes. The vertical axis is
the quality defined in Eq. (4).

Figure 3: Computational results for different size networks with parameters:
$\lambda=50$, $n=20$, $r_{0}=0.05|A|$. Each point corresponds to a graph. The
horizontal axis corresponds to the number of nodes. The vertical axis is the
quality defined in Eq. (4).

Figure 4: Computational results for different size networks with parameters:
$\lambda=50$, $n=40$, $r_{0}=0.0025|A|$. Each point corresponds to a graph.
The horizontal axis corresponds to the number of nodes. The vertical axis is
the quality defined in Eq. (4).

##  VI Discussion

###  VI-A Solving on Ising Processing Hardware

The current results shown in this paper rely on using an exact solver to solve
subproblems. We have previously attempted to utilize classical _heuristic_
QUBO solvers, such as D-Wave’s _qbsolv_ and D-Wave Ocean, to solve these
subproblems, but such solvers provided poor quality solutions even on much
smaller problem sizes that those we experiment with. Additionally, due to lack
of accessibility to massive experiments with the IPU hardware, utilizing exact
solvers is necessary to demonstrate the efficiency of our algorithm. We hope
that as IPUs continue to develop in capacity, their superior solution quality
and speed when compared to classical heuristic solvers will be able to make
full use of the algorithm presented in this paper, generating solutions that
are able to beat classical exact solvers consistently.

###  VI-B Algorithmic Challenges and Obstacles

Refinement algorithms for solving SPNI may face several challenges, the most
notable being the problem of proper budget allocation amongst partitions. If
an algorithm decides to improperly allocate budget to a subproblem, for
example deciding to interdict an arc where this unit of budget would ideally
be spent much further away in the graph, it proves to be very difficult to
shift this budget significantly far away from its current position. Within our
refinement approach, the movement of interdicted arcs is primarily dictated by
random partitioning—for example, in one iteration, the subproblem solver may
determine that a specific arc should be interdicted under the context of the
current partititioning configuration, i.e. within the context of the current
subproblem the arc lies in. But when we repartition in the next iteration, the
subproblem solver reevaluates the placement of interdicted arcs—the cost unit
used to interdict the previous arc may be relocated elsewhere based on the new
subproblem formulation. Since the cost unit used for an interdicted arc may
only shift to arcs relatively near this arc, it is unlikely for the cost unit
to travel far away in the graph. To resolve this issue, a multilevel approach
widely known in combinatorial scientific computing and applied graph
algorithms [4] may be taken. Multilevel algorithms create a hierarchy of
coarse problem representations, find a solution to the coarsest (smallest and
most compressed) problem, then gradually derive a solution to the original
finest problem by projecting the solution back from the coarse to fine levels
and refining it using our algorithm, thereby potentially aiding in proper
budget allocation. This approach has seen great success in solving various
(hyper)graph problems [1, 26, 28] including hybrid quantum-classical multi-
and single-level refinement for partitioning and community detection [29, 31],
indicating that its application to SPNI may improve solution quality.

###  VI-C Performance

One of our main goals was to develop a trivial IPU parallelization friendly
algorithm. Our current implementation of the proposed refinement does not take
advantage of the fact that subproblems produced by a given partitioning can be
solved independently of one another, since they only depend on the current
working solution. Consequently, subproblems can be solved in parallel rather
than sequentially, offering a potentially large optimization to the current
performance of the algorithm.

###  VI-D Generalizability

The algorithm introduced in this paper does not handle SPNI problem instances
in which the binary interdiction of a given arc may cost more than one unit of
the interdictor’s budget. The algorithm can currently handle instances where
an arc between two nodes can take on multiple values of interdiction (e.g., if
interdicted once add $1$ to the length, if interdicted twice add $2$ to the
length) through the introduction of multi-edges with appropriate lengths, but
instances in which a binary decision needs to be made (i.e., interdicted or
not) and differing levels of interdiction cost are not supported.

Another scenario which the algorithm is unable to handle is in the case of
different resource types. For example of such a scenario, let $r_{0}$ and
$r_{1}$ represent the budget for two different types of resources, then arc
$k_{0}$ may cost $1$ of $r_{0}$ and $2$ of $r_{1}$, but arc $k_{1}$ may cost
$0$ of $r_{0}$ and $3$ of $r_{1}$. This restriction is not much of a concern,
however, considering that in one of the most prominent works regarding SPNI,
the authors of [16] also do not accommodate this case for their algorithmic
cases.

##  VII Conclusion

Solving SPNI holds crucial importance in the defense, engineering and other
domains, addressing challenges in securing critical infrastructure from
terrorist attacks, natural disasters and disrupting enemy supply networks
among other applications. Current solutions to SPNI, however, are often too
slow to be scalable to large networks, and thus prove to be impractical for
real world purposes. In this paper, we have introduced a novel decomposition
approach addressing SPNI harnessing the rapid solving abilities of IPU
hardware to yield an approximate solution. We have further shown that
solutions generated by our decomposition algorithm are extremely close in
quality to state-of-the-art integer programming methods. We anticipate that as
IPU hardware advancements continue to improve speed in solving optimization
problems, our algorithm will enable an even quicker approximation of SPNI
thereby facilitating a more efficient approach for real world scenarios of
SPNI. More work is under way to produce higher quality approximations in even
shorter amounts of time and extend the algorithm to more general forms of
SPNI.
Reproducibility:  Our results and code are available at
https://github.com/krishxmatta/network-interdiction.

## Acknowledgements

We would like to thank J. Cole Smith and Yongjia Song, the authors of [30],
for their input and discussion about the current status of the network
interdiction solvers.

This work was supported in part with funding from the Defense Advanced
Research Projects Agency (DARPA). The views, opinions and/or findings
expressed are those of the author and should not be interpreted as
representing the official views or policies of the Department of Defense or
the U.S.Government.

This research was supported in part through the use of DARWIN computing
system: DARWIN - A Resource for Computational and Data-intensive Research at
the University of Delaware and in the Delaware Region, which is supported by
NSF Grant #1919839.

## References

  * [1] A. Abou-Rjeili and G. Karypis.  Multilevel algorithms for partitioning power-law graphs.  In Proceedings 20th IEEE International Parallel & Distributed Processing Symposium, pages 10 pp.–, 2006.
  * [2] Mine Altunay, Sven Leyffer, Jeffrey T Linderoth, and Zhen Xie.  Optimal response to attacks on the open science grid.  Computer Networks, 55(1):61–73, 2011.
  * [3] Michael O. Ball, Bruce L. Golden, and Rakesh V. Vohra.  Finding the most vital arcs in a network.  Operations Research Letters, 8(2):73–76, 1989.
  * [4] Cédric Chevalier and Ilya Safro.  Comparison of coarsening schemes for multilevel graph partitioning.  pages 191–205, 01 2009.
  * [5] Carleton Coffrin, Harsha Nagarajan, and Russell Bent.  Evaluating ising processing units with integer programming.  In Integration of Constraint Programming, Artificial Intelligence, and Operations Research, pages 163–181. Springer International Publishing, 2019.
  * [6] H.W Corley and David Y Sha.  Most vital links and nodes in weighted networks.  Operations Research Letters, 1(4):157–160, 1982.
  * [7] John Forrest and Robin Lougee-Heimer.  Cbc user guide.  In Emerging theory, methods, and applications, pages 257–277. INFORMS, 2005.
  * [8] D. R. Fulkerson and Gary C. Harding.  Maximizing the minimum source-sink path subject to a budget constraint.  Math. Program., 13(1):116–118, dec 1977.
  * [9] P. M. Ghare, D. C. Montgomery, and W. C. Turner.  Optimal interdiction policy for a flow network.  Naval Research Logistics Quarterly, 18(1):37–45, 1971.
  * [10] Fred W. Glover and Gary A. Kochenberger.  A tutorial on formulating QUBO models.  CoRR, abs/1811.11538, 2018.
  * [11] Noam Goldberg, Sven Leyffer, and Ilya Safro.  Optimal response to epidemics and cyber attacks in networks.  Networks, 66(2):145–158, 2015.
  * [12] Bruce Golden.  A problem in network interdiction.  Naval Research Logistics Quarterly, 25(4):711–713, 1978.
  * [13] Gurobi Optimization, LLC.  Gurobi Optimizer Reference Manual, 2023.
  * [14] Aric Hagberg, Pieter Swart, and Daniel S Chult.  Exploring network structure, dynamics, and function using networkx.  Technical report, Los Alamos National Lab.(LANL), Los Alamos, NM (United States), 2008.
  * [15] Dian Huang, Zhaofang Mao, Kan Fang, and Lin Chen.  Solving the shortest path interdiction problem via reinforcement learning.  International Journal of Production Research, 61(1):31–48, 2023\.
  * [16] Eitan Israeli and R. Kevin Wood.  Shortest-path network interdiction.  Networks, 40(2):97–111, 2002.
  * [17] Sahar Karimi and Pooya Ronagh.  Practical integer-to-binary mapping for quantum annealers.  Quantum Information Processing, 18(4), feb 2019.
  * [18] George Karypis and Vipin Kumar.  METIS: A Software Package for Partitioning Unstructured Graphs, Partitioning Meshes, and Computing Fill-Reducing Orderings of Sparse Matrices, September 1998.
  * [19] Sven Leyffer and Ilya Safro.  Fast response to infection spread and cyber attacks on large-scale networks.  Journal of Complex Networks, 1(2):183–199, 2013.
  * [20] Xiaoyuan Liu, Hayato Ushijima-Mwesigwa, Avradip Mandal, Sarvagya Upadhyay, Ilya Safro, and Arnab Roy.  Leveraging special-purpose hardware for local search heuristics.  Computational Optimization and Applications, 82(1):1–29, 2022.
  * [21] K. Malik, A.K. Mittal, and S.K. Gupta.  The k most vital arcs in the shortest path problem.  Operations Research Letters, 8(4):223–227, 1989.
  * [22] Naeimeh Mohseni, Peter L McMahon, and Tim Byrnes.  Ising machines as hardware solvers of combinatorial optimization problems.  Nature Reviews Physics, 4(6):363–379, 2022.
  * [23] Apurba K. Nandi, Hugh R. Medal, and Satish Vadlamani.  Interdicting attack graphs to protect organizations from cyber attacks: A bi-level defender–attacker model.  Computers & Operations Research, 75:118–131, 2016.
  * [24] Florian Neukart, Gabriele Compostella, Christian Seidel, David von Dollen, Sheir Yarkoni, and Bob Parney.  Traffic flow optimization using a quantum annealer.  Frontiers in ICT, 4, 2017.
  * [25] Claudio M. Rocco S. and José Emmanuel Ramirez-Marquez.  A bi-objective approach for shortest-path network interdiction.  Computers & Industrial Engineering, 59(2):232–240, 2010.
  * [26] Ilya Safro and Boris Temkin.  Multiscale approach for the network compression-friendly ordering.  Journal of Discrete Algorithms, 9(2):190–202, 2011.
  * [27] Javier Salmeron, Kevin Wood, and Ross Baldick.  Worst-case interdiction analysis of large-scale electric power grids.  IEEE Transactions on Power Systems, 24(1):96–104, 2009.
  * [28] Ruslan Shaydulin, Jie Chen, and Ilya Safro.  Relaxation-based coarsening for multilevel hypergraph partitioning.  Multiscale Modeling & Simulation, 17(1):482–506, 2019.
  * [29] Ruslan Shaydulin, Hayato Ushijima-Mwesigwa, Ilya Safro, Susan Mniszewski, and Yuri Alexeev.  Network community detection on small quantum computers.  Advanced Quantum Technologies, 2(9):1900029, 2019.
  * [30] J Cole Smith and Yongjia Song.  A survey of network interdiction models and algorithms.  European Journal of Operational Research, 283(3):797–811, 2020\.
  * [31] Hayato Ushijima-Mwesigwa, Ruslan Shaydulin, Christian FA Negre, Susan M Mniszewski, Yuri Alexeev, and Ilya Safro.  Multilevel combinatorial optimization across quantum architectures.  ACM Transactions on Quantum Computing, 2(1):1–29, 2021.