# CatFLW: Cat Facial Landmarks in the Wild Dataset

George Martvel Nareed Farhat Ilan Shimshoni Anna Zamansky
University of Haifa, Israel

###### Abstract

Animal affective computing is a quickly growing field of research, where only
recently first efforts to go beyond animal tracking into recognizing their
internal states, such as pain and emotions, have emerged. In most mammals,
facial expressions are an important channel for communicating information
about these states. However, unlike the human domain, there is an acute lack
of datasets that make automation of facial analysis of animals feasible.

This paper aims to fill this gap by presenting a dataset called Cat Facial
Landmarks in the Wild (CatFLW) which contains 2016 images of cat faces in
different environments and conditions, annotated with 48 facial landmarks
specifically chosen for their relationship with underlying musculature, and
relevance to cat-specific facial Action Units (CatFACS). To the best of our
knowledge, this dataset has the largest amount of cat facial landmarks
available.

In addition, we describe a semi-supervised (human-in-the-loop) method of
annotating images with landmarks, used for creating this dataset, which
significantly reduces the annotation time and could be used for creating
similar datasets for other animals.

The dataset is available on request.

##  1 Introduction

Emotion recognition based on facial expression and body language analysis are
important and challenging topics addressed within the field of human affective
computing. Comprehensive surveys cover analysis of facial expressions [27],
and body behaviors [33], with the recent trend being multi-modal emotion
recognition approaches [35]. Automated facial expression analysis is also
studied in pain research [2], and is most crucial for addressing pain in
infants[40].

In animals, facial expressions are known to be produced by most mammalian
species[31, 14]. Analogously to the human domain, they are assumed to convey
information about emotional states, and are increasingly studied as potential
indicators of subjective states in animal emotion and welfare research [5, 29,
3, 19].

Although research in animal affective computing has so far lagged behind the
human domain, this is rapidly changing in recent years [22]. Broome et al [7]
provide a comprehensive review of more than twenty studies addressing
automated recognition of animals’ internal states such as emotions and pain,
with the majority of these studies focusing on the analysis of animal faces in
different species, such as rodents, dogs, horses and cats.

Cats are of specific interest in the context of pain, as they are one of the
most challenging species in terms of pain assessment and management due to a
reduced physiological tolerance and adverse effects to common veterinary
analgesics[25], a lack of strong consensus over key behavioural pain
indicators[30] and human limitations in accurately interpreting feline facial
expressions[12]. Three different manual pain assessment scales have been
developed and validated in English for domestic cats: the UNESP-Botucatu
multidimensional composite pain scale[6], the Glasgow composite measure pain
scale (CMPS, [34]) and the Feline Grimace Scale (FGS) [16]. The latter was
further used for a comparative study in which human’s assignment of FGS to
cats during real time observations and then subsequent FGS scoring of the same
cats from still images were compared. It was shown that there was no
significant difference between the scoring methods [15], indicating that
facial images can be a reliable medium from which to assess pain.

Finka and colleagues[18] were the first to apply 48 geometric landmarks to
identify and quantify facial shape change associated with pain in cats. These
landmarks were based on both the anatomy of cat facial musculature, and the
range of facial expressions generated as a result of facial action units
(CatFACS)[8]. The authors manually annotated the landmarks, and used
statistical methods (PCA analysis) to establish a relationship between PC
scores and a well-validated measure of pain in cats. Feighelstein et al [17]
used the dataset and landmark annotations from [18] to automate pain
recognition in cats using machine learning models based on geometric landmarks
grouped into multivectors, reaching accuracy of above 72%, which was
comparable to an alternative approach using deep learning. This indicates that
the use of the 48 landmarks from [18] can be a reliable method for pain
recognition from cat facial images. The main challenge in this method remains
the automation of facial landmarks detection. In general, one of the most
crucial gaps in animal affective computing today, as emphasized in [7], is the
lack of publically available datasets. To the best of our knowledge, all
datasets containing animal faces available today have very few facial
landmarks, as discussed in the next section.

Figure 1: Annotated Cat’s Face. Image of a cat with face bounding box and 48
facial landmarks.

To address this gap for the feline species, where facial analysis can help
address acute problems in the context of pain assessment, we present the Cat
Facial Landmarks in the Wild (CatFLW) dataset, containing images of cat faces
in various conditions and environments with face bounding boxes and 48 facial
landmarks of [18, 17] for each image (see Figure 1). We hope that this dataset
can serve as a starting point for creating similar datasets for various
animals, as well as push forward automatic detection of pain and emotions in
cats and other feline species.

Our main contributions are, therefore, the following:

  * •

We present a dataset containing 2016 annotated images of cats with 48
landmarks, which are based on the cat’s facial anatomy and allows the study of
complex morphological features of cats’ faces.

  * •

We describe an AI-based ‘human-in-the-loop’ (also known as cooperative machine
learning[21]) method of facial landmark annotation, and show that it
significantly reduces annotation time per image.

##  2 Related Datasets

Our dataset is based on the original dataset collected by Zhang et al. [41],
that contains 10,000 images of cats annotated with 9 facial landmarks. It
includes a wide variety of different cat breeds in different conditions, which
can provide good generalization when training computer vision models, however
some images from it depict several animals, have visual interference in front
of animal faces or cropped (according to our estimates 10-15%). In [41] Zhang
et al. 9 landmarks are labeled for each image, which is sufficient for
detecting and analyzing general information about animal faces (for example,
the tilt of the head or the direction of movement), but is not enough for
analyzing complex movements of facial muscles.

Dataset | Animal | Size | Facial Landmarks
---|---|---|---
Khan et al. [23] | Various | 21900 | 9
Zhang et al. [41] | Cat | 10000 | 9
Liu et al. [28] | Dog | 8351 | 8
Cao et al. [9] | Various | 5517 | 5
Mougeot et al. [32] | Dog | 3148 | 3
Sun et al. [37] | Cat | 1706 | 15
Yang et al. [39] | Sheep | 600 | 8
CatFLW | Cat | 2016 | 48
Table 1: Comparison of animal facial landmarks datasets

Sun et al. [37] used the same dataset in their work, expanding the annotation
to 15 facial landmarks. To the best of our knowledge, only 1,706 out of the
declared 10,000 images are publicly available.

Khan et al. [23] collected the AnimalWeb dataset, consisting of an impressive
number of 21,900 images annotated with 9 landmarks. Despite the wide range of
represented species, only approximately 450 images can be attributed to feline
species.

Other relevant datasets, containing less than 9 facial landmarks, are
mentioned in Table 1, including those in [28, 9, 32, 39]. For comparison,
popular datasets for human facial landmark detection [4, 26] have several
dozens of landmarks. There are also other datasets suitable for face detection
and recognition of various animal species [11, 13, 20, 24, 10], but they have
no facial landmark annotations.

##  3 The CatFLW Dataset

The CatFLW dataset consists of 2016 images selected from the dataset in [41]
using the following inclusion criterion which optimize the training of
landmark detection models: image should contain a single fully visible cat
face, where the cat is in non-laboratory conditions (‘in the wild’). Thus many
images from the original dataset were discarded (Figure 2 shows the examples).

Figure 2: Examples of not suitable images from [41]. The image is considered
inapplicable when it is impossible to visually determine the location of the
facial landmark with full confidence or uncertainty is created due to the
presence of several animals. Figure 3: IOD distributions for CatFLW and
AnimalWeb. Both datasets have approximately the same mean inter-ocular
distance, however, the latter has both a high number of small-sized faces and
a wide variation in large faces.

The images in this dataset have a large variation of colors, breeds and
environments, as well as of scale, position and angle of head tilts. Figure 3
shows the distribution of inter-ocular distances (IOD — distance between the
two medial canthi of each eye) on CatFLW and AnimalWeb [23]. CatFLW dataset
has similar characteristics to AnimalWeb, except that the latter has about
twice the standard deviation. This can be explained by the large species
diversity, as well as a greater variation in the size of the images.

###  3.1 The Landmark Annotation Process

After selection and filtering of images, a bounding box and 48 facial
landmarks were placed on each image using the Labelbox platform[1].

Bounding box. In the current dataset, the bounding box is not tied to facial
landmarks (this is usually done by selecting the outermost landmarks and
calculating the difference in their coordinates), but is visually placed so
that the entire face of the animal fits into the bounding box as well as about
10% of the space around the face. This margin is made because when studying
detection models on the initial version of the dataset, it was noticed that
some of them tend to crop the faces and choose smaller bounding boxes. In the
case when face detection is performed for the further localization of facial
landmarks, the clipping of faces can lead to the disappearance from the image
of important parts of the cat’s face, such as the tips of the ears or the
mouth. Such an annotation still makes it possible to construct bounding boxes
by outermost landmarks if needed.

Landmarks. 48 facial landmarks introduced in Finka et al. [18] and used in
Feighelstein et al. [17] were manually placed on each image. These landmarks
were specifically chosen for their relationship with underlying musculature
and anatomical features, and relevance to cat-specific Facial Action Units
(CatFACS) [8].

Annotation Process. The process of landmark annotation had several stages: at
first, 10% images were annotated by an annotator who has an extensive
experience labeling facial landmarks for various animals, then they were
annotated by second expert with the same annotation instructions. The
landmarks were then compared to verify the internal validity and reliability
via the Inter Class Correlation Coefficient ICC2 [36], and reached a strong
agreement between the annotators with a score of 0.998.

Finally, the remaining part of images were annotated by the first annotator
using the ‘human-in-the-loop’ method described below.

After the annotation, a review and correction of landmarks and bounding boxes
were performed, all images were re-filtered and verified.

Figure 4: Left: The distribution of annotation time per image for different
batches. The predictions of the first model significantly reduced the
annotation time, while the subsequent refinement by the second model reduced
the variance, without changing the median value too much. Right: The
percentage of images from the batch for which a predicted landmark was not
shifted. The more accurate the model’s predictions are, the less time is
needed for annotation: correctly predicted landmarks do not need to be
shifted.

###  3.2 AI-assisted Annotation

Khan et al.[23] spent approximately 5,408 man-hours annotating all the images
in the AnimalWeb dataset (each annotation was obtained by taking the median
value of the annotations of five or more error-prone volunteers). Roughly, it
is 3 minutes for the annotation of one image (9 landmarks) or 20 seconds for
one landmark (taken per person). Our annotation process took $\sim$140 hours,
which is about 4.16 minutes per image (48 landmarks) and only 5.2 seconds per
landmark. Such performance is achieved due to a semi-supervised ’human-in-the-
loop’ method that uses the predictions of a gradually-trained model as a basis
for annotation.

To assess the impact of AI on the annotation process, 3 batches were created
from the data. The first one consisted of two hundred images and was annotated
without AI assistance. Next, we trained the EfficientNet-based [38] model on
the dataset from [18] (about 600 images) combined with the first annotated
batch ($\sim$800 images total). After that, we used this model (v1) to predict
the landmarks on the second batch ($\sim$900 images). Based on the
predictions, landmarks were annotated for the second batch. This approach made
it possible to reduce the average time for annotation of one image between the
first and second batches by 35%, since for annotation it became necessary only
to adjust the position of the landmark, and not to place it from scratch.
However, the predictions obtained turned out to be inaccurate, so the first
and second batches were combined to train the next version of the model (v2)
with the same architecture ($\sim$1700 images total). After training, the v2
model was used to predict the landmarks on the third batch, and then an
annotation was performed based on these predictions. The results obtained for
the time spent on annotation of one image are shown in Figure 4(left).
Increasing the accuracy of predicting the position of a landmark reduces the
time required to adjust it, sometimes to zero if the annotator considers the
position to be completely correct. Figure 4(right) shows the distribution of
the percentage of images in which a specific facial landmark was not shifted
by the annotator: between the two versions of the model, the increase is
approximately doubled.

##  4 Conclusion

In this paper we present a CatFLW dataset of annotated cat faces. It includes
2016 images of cats in various conditions with face bounding boxes and 48
facial landmarks for each cat. Each of the landmarks corresponds to a certain
anatomical feature, which makes it possible to use these landmarks for various
kinds of tasks on evaluating the internal state of cats. The proposed AI-
assisted method of image annotation has significantly reduced the time for
annotation of landmarks and reduced the amount of manual work during the
preparation of the dataset. The created dataset opens up opportunities for the
development of various computer vision models for detecting facial landmarks
of cats and can be the starting point for research in the field of behaviour
observation and well-being of cats and animals in general.

## Acknowledgements

The research was supported by the Data Science Research Center at the
University of Haifa. We thank Ephantus Kanyugi for his contribution with data
annotation and management. We thank Yaron Yossef for his technical support.

## References

  * [1] Labelbox, ”labelbox,” online, 2023. [online]. available:.  https://labelbox.com.
  * [2] Rasha M Al-Eidan, Hend Suliman Al-Khalifa, and AbdulMalik S Al-Salman.  Deep-learning-based models for pain recognition: A systematic review.  Applied Sciences, 10:5984, 2020.
  * [3] Niek Andresen, Manuel Wöllhaf, Katharina Hohlbaum, Lars Lewejohann, Olaf Hellwich, Christa Thöne-Reineke, and Vitaly Belik.  Towards a fully automated surveillance of well-being status in laboratory mice using deep learning: Starting with facial expression analysis.  PLoS One, 15(4):e0228059, 2020.
  * [4] Peter N Belhumeur, David W Jacobs, David J Kriegman, and Neeraj Kumar.  Localizing parts of faces using a consensus of exemplars.  IEEE transactions on pattern analysis and machine intelligence, 2013\.
  * [5] Tali Boneh-Shitrit, Marcelo Feighelstein, Annika Bremhorst, Shir Amir, Tomer Distelfeld, Yaniv Dassa, Sharon Yaroshetsky, Stefanie Riemer, Ilan Shimshoni, Daniel S Mills, et al.  Explainable automated recognition of emotional states from canine facial expressions: the case of positive anticipation and frustration.  Scientific reports, 12(1):22611, 2022.
  * [6] Juliana T Brondani, Khursheed R Mama, Stelio PL Luna, Bonnie D Wright, Sirirat Niyom, Jennifer Ambrosio, Pamela R Vogel, and Carlos R Padovani.  Validation of the english version of the unesp-botucatu multidimensional composite pain scale for assessing postoperative pain in cats.  BMC Veterinary Research, 9(1):1–15, 2013.
  * [7] Sofia Broome, Marcelo Feighelstein, Anna Zamansky, Carreira G Lencioni, Haubro P Andersen, Francisca Pessanha, Marwa Mahmoud, Hedvig Kjellström, and Albert Ali Salah.  Going deeper than tracking: A survey of computer-vision based recognition of animal pain and emotions.  International Journal of Computer Vision, 131(2):572–590, 2023\.
  * [8] Cátia C Caeiro, Anne M Burrows, and Bridget M Waller.  Development and application of catfacs: Are human cat adopters influenced by cat facial expressions?  Applied Animal Behaviour Science, 2017.
  * [9] Jinkun Cao, Hongyang Tang, Haoshu Fang, Xiaoyong Shen, Cewu Lu, and Yu-Wing Tai.  Cross-domain adaptation for animal pose estimation.  CoRR, abs/1908.05806, 2019.
  * [10] Peng Chen, Pranjal Swarup, Wojciech M Matkowski, Kong Adams WK, Su Han, Zhihe Zhang, and Hou Rong.  A study on giant panda recognition based on images of a large proportion of captive pandas.  Ecol Evol., 2020.
  * [11] Melanie Clapham, Ed Miller, Mary Nguyen, and Russel C Van Horn.  Multispecies facial detection for individual identification of wildlife: a case study across ursids.  Mamm Biol 102, 943–955, 2022.
  * [12] Lauren Dawson, Joanna Cheal, Lauren Niel, and Georgia Mason.  Humans can identify cats’ affective states from subtle facial expressions.  Animal Welfare, 28(4):519–531, 2019.
  * [13] Debayan Deb, Susan Wiper, Sixue Gong, Yichun Shi, Cori Tymoszek, Alison Fletcher, and Anil K. Jain.  Face recognition: Primates in the wild.  In 2018 IEEE 9th International Conference on Biometrics Theory, Applications and Systems (BTAS), pages 1–10, 2018.
  * [14] Rui Diogo, Virginia Abdala, N Lonergan, and BA Wood.  From fish to modern humans–comparative anatomy, homologies and evolution of the head and neck musculature.  Journal of Anatomy, 213(4):391–424, 2008.
  * [15] Marina C Evangelista, Javier Benito, Beatriz P Monteiro, Ryota Watanabe, Graeme M Doodnaught, Daniel SJ Pang, and Paulo V Steagall.  Clinical applicability of the feline grimace scale: real-time versus image scoring and the influence of sedation and surgery.  PeerJ, 8:e8967, 2020.
  * [16] Marina C Evangelista, Ryota Watanabe, Vivian SY Leung, Beatriz P Monteiro, Elizabeth O’Toole, Daniel SJ Pang, and Paulo V Steagall.  Facial expressions of pain in cats: the development and validation of a feline grimace scale.  Scientific reports, 9(1):1–11, 2019.
  * [17] Marcelo Feighelstein, Ilan Shimshoni, Lauren Finka, Stelio P. Luna, Daniel Mills, and Anna Zamansky.  Automated recognition of pain in cats.  Scientific Reports, 12, 2022.
  * [18] Lauren R Finka, Stelio P Luna, Juliana T Brondani, Yorgos Tzimiropoulos, John McDonagh, Mark J Farnworth, Marcello Ruta, and Daniel S Mills.  Geometric morphometrics for the study of facial expressions in non-human animals, using the domestic cat as an exemplar.  Scientific reports, 9(1):1–12, 2019.
  * [19] Karina B Gleerup, Björn Forkman, Casper Lindegaard, and Pia H Andersen.  An equine pain face.  Veterinary anaesthesia and analgesia, 42(1):103–114, 2015.
  * [20] Songtao Guo, Pengfei Xu, Qiguang Miao, Guofan Shao, Colin A Chapman, Xiaojiang Chen, Gang He, Dingyi Fang, He Zhang, Yewen Sun, Zhihui Shi, and Baoguo Li.  Automatic identification of individual primates with deep learning techniques.  iScience, 2020.
  * [21] Alexander Heimerl, Katharina Weitz, Tobias Baur, and Elisabeth André.  Unraveling ml models of emotion with nova: Multi-level explainable ai for non-experts.  IEEE Transactions on Affective Computing, 13(3):1155–1167, 2020\.
  * [22] Hilde I Hummel, Francisca Pessanha, Albert Ali Salah, Thijs JPAM van Loon, and Remco C Veltkamp.  Automatic pain detection on horse and donkey faces.  In FG, 2020.
  * [23] Muhammad Haris Khan, John McDonagh, Salman H Khan, Muhammad Shahabuddin, Aditya Arora, Fahad Shahbaz Khan, Ling Shao, and Georgios Tzimiropoulos.  Animalweb: A large-scale hierarchical dataset of annotated animal faces.  CVRR, abs/1909.04951, 2019.
  * [24] Matthias Körschens, Björn Barz, and Joachim Denzler.  Towards automatic identification of elephants in the wild.  CoRR, abs/1812.04418, 2018.
  * [25] B Duncan X Lascelles and Sheilah A Robertson.  Djd-associated pain in cats: what can we do to promote patient comfort?  Journal of Feline Medicine & Surgery, 12(3):200–212, 2010.
  * [26] Vuong Le, Jonathan Brandt, Zhe Lin, Lubomir Bourdev, and Thomas S Huang.  Interactive facial feature localization.  Proceedings of European Conference on Computer Vision (ECCV), Springer, 2012.
  * [27] Shan Li and Weihong Deng.  Deep facial expression recognition: A survey.  IEEE transactions on affective computing, 13(3):1195–1215, 2022\.
  * [28] Jiongxin Liu, Angjoo Kanazawa, David Jacobs, and Peter Belhumeur.  Dog breed classification using part localization.  In European conference on computer vision, pages 172–185. Springer, 2012.
  * [29] Katrina Merkies, Chloe Ready, Leanne Farkas, and Abigail Hodder.  Eye blink rates and eyelid twitches as a non-invasive measure of stress in the domestic horse.  Animals (Basel), 2019.
  * [30] Isabella Merola and Daniel S Mills.  Behavioural signs of pain in cats: an expert consensus.  PloS one, 11(2):e0150040, 2016.
  * [31] Alexander Mielke, Bridget M Waller, Claire Pérez, Alan V Rincon, Julie Duboscq, and Jérôme Micheletta.  Netfacs: Using network science to understand facial communication systems.  Behavior Research Methods 54, 1912–1927, 2022.
  * [32] Guillaume Mougeot, Dewei Li, and Shuai Jia.  A deep learning approach for dog face verification and recognition.  Lecture Notes in Computer Science, 2019.
  * [33] Fatemeh Noroozi, Ciprian Adrian Corneanu, Dorota Kamińska, Tomasz Sapiński, Sergio Escalera, and Gholamreza Anbarjafari.  Survey on emotional body gesture recognition.  IEEE transactions on affective computing, 12(2):505–523, 2018.
  * [34] Jacky Reid, Mariann E Scott, Gillian Calvo, and Andrea M Nolan.  Definitive glasgow acute pain scale for cats: validation and intervention level.  Veterinary Record, 180(18):449–449, 2017.
  * [35] Garima Sharma and Abhinav Dhall.  A survey on automatic multimodal emotion recognition in the wild.  In Gloria Phillips-Wren, Anna Esposito, and Lakhmi C. Jain, editors, Advances in Data Science: Methodologies and Applications, pages 35–64. Springer International Publishing, Cham, 2021.
  * [36] Patrick E Shrout and Joseph L Fleiss.  Intraclass correlations: Uses in assessing rater reliability.  Psychological Bulletin, 86(2):420–428, 1979.
  * [37] Yifan Sun and Noboru Murata.  Cafm: A 3d morphable model for animals.  IEEE Winter Applications of Computer Vision Workshops (WACVW), 2020\.
  * [38] Mingxing Tan and Quoc V Le.  Efficientnet: Rethinking model scaling for convolutional neural networks.  CoRR, abs/1905.11946, 2019.
  * [39] Heng Yang, Renqiao Zhang, and Peter Robinson.  Human and sheep facial landmarks localisation by triplet interpolated features.  CVRR, abs/1509.04954, 2015.
  * [40] Ghada Zamzmi, Rangachar Kasturi, Dmitry Goldgof, Ruicong Zhi, Terri Ashmeade, and Yu Sun.  A review of automated pain assessment in infants: features, classification tasks, and databases.  IEEE reviews in biomedical engineering, 11:77–96, 2017.
  * [41] Weiwei Zhang, Jian Sun, and Xiaoou Tang.  Cat head detection - how to effectively exploit shape and texture features.  ECCV, 2008.


# Quantum computer based Feature Selection in Machine Learning

Gerhard Hellstern DHBW Ravensburg, Marienplatz 2, D-88212 Ravensburg, Germany
Vanessa Dehn Fraunhofer-Institut für Angewandte Festkörperphysik IAF,
Tullastr. 72, D-79108 Freiburg, Germany  Martin Zaefferer DHBW Ravensburg,
Marienplatz 2, D-88212 Ravensburg, Germany

###### Abstract

The problem of selecting an appropriate number of features in supervised
learning problems is investigated in this paper. Starting with common methods
in machine learning, we treat the feature selection task as a quadratic
unconstrained optimization problem (QUBO), which can be tackled with classical
numerical methods as well as within a quantum computing framework. We compare
the different results in small-sized problem setups. According to the results
of our study, whether the QUBO method outperforms other feature selection
methods depends on the data set. In an extension to a larger data set with 27
features, we compare the convergence behavior of the QUBO methods via quantum
computing with classical stochastic optimization methods. Due to persisting
error rates, the classical stochastic optimization methods are still superior.

## 1\. Introduction and literature review

When using supervised machine learning to tackle real-world problems one is
faced with the task of selecting the "right" features for the learning
procedure. When too few features are chosen, this may lead to poor results for
the algorithms. When too many features are selected, the training time is
negatively affected, which may lead to unstable results. Moreover, the more
features are used, the more effort is required to maintain and control these
features for a deployed machine learning model. Just testing all possible
combinations of features with a brute-force approach is usually not feasible,
since the number of combinations grows exponentially with the number of
features.

In the past, different methods have been proposed to tackle this problem.
Brown et al. provide a compilation of several approaches where mutual
information is used as a dependency measure, introduced within an information-
theoretic framework [6]. An extensive survey of feature selection methods that
are applied to different data sets to demonstrate their behavior is given in
[7] and [27].

The idea of using quantum computers for solving this problem is introduced in
[20]: There, the feature selection problem has been reformulated as a
quadratic unconstrained optimization problem (QUBO). In addition, this problem
was solved via classical optimization methods as well as with a quantum
annealing approach. The different methods were applied to a well-known data
set, the "German Credit Data" of UC Irvine [10].

The annealing approach of solving the QUBO problem for feature selection was
also employed in [22], where the authors investigate feature selection for
recommender systems. With the same approach, other data sets were explored in
[9]. Results that combine results obtained with annealing methods as well as
with a gate-based approach to quantum computing are reported in [21]. An
approach that examines an unconstrained black box binary optimization and
applies it to feature selection is presented in [30]. Finally, [26] compares
different algorithms (annealing, gate-based and classical) to solve the QUBO
problem of feature selection in a small-size setup.

In this paper, we report on the following approach to investigate the feature
selection problem: Using small-sized data sets we first compare well-known
machine-learning methods for feature selection with the brute-force method.

After reformulating the feature selection problem as a QUBO, hereby allowing
different dependency measures between the features and also between the
features and the target, we investigate the stability of the solutions when
using different values of the undetermined weighting parameter with classical
optimization methods. Since for this problem size the exact solution for the
features selection problem can be determined, we investigate to what extent
this solution can be achieved on the one hand with the commonly used feature
selection methods and on the other hand with the way via QUBO optimization.

In order to solve the QUBO with the gate-based approach of quantum computing,
we use the QAOA- and the VQE-ansatz. In addition to using the out-of-the-box
QAOA algorithm of IBM, we compare the results with a customized solution
scheme for QAOA. Hereby it can be shown whether and to what extent the QAOA
algorithm is a feasible way to solve the QUBO problem and how strong the
goodness of the solutions depends on the QAOA depth as well as on the
dependence measures used.

After discussing our results on small-sized data sets, we scale our methods up
to 27 features (i.e., 27 qubits) and report optimization results on different
physical quantum computers of IBM. For this setup, we also compare the results
to several classical stochastic optimization methods. At this point, we
investigate, on the one hand, whether 27 qubits on a real quantum computer can
yield any meaningful results beyond random noise, given its error rates. On
the other hand, we focus on comparing the convergence of the QAOA algorithm
with a selection of classical stochastic optimization methods.

## 2\. Feature Selection in Supervised Learning

### 2.1. Machine Learning approaches

In the following, we restrict ourselves to binary classification as one of the
most important topics in supervised learning:

Given are a real $m\times n$ matrix $X$ consisting of $n$ features and $m$
training examples and a target vector $y$; $y$ is of dimension $m$ and each
element $y_{i}\in\left\\{0,1\right\\}$. A classification model is a mapping
$f:X\in A\rightarrow B$, i.e. from an input space $A$ to an output space $B$.
In order to find the best model $f$, we use a loss function $g:A\times
B\rightarrow\mathbb{R}$ that has to be minimized. In the case of a binary
classification problem, the binary cross-entropy is mostly used as a loss
function. For a thorough discussion of machine learning, see e.g. [2].

Given a training set $(X,y)$, the following problem arises in practice: How
many and which features should be included in the training of a particular
model in order to do the best job? Naively, one would expect that considering
all features is naturally the best choice. However, there are several reasons
why this is in general not the case:

  * •

Features with a high dependence on each other can lead to an unstable training
process.

  * •

Each feature used by an algorithm must be prepared and pre-processed and
controlled in a productive environment. So, from an economical point of view,
one should not consider more features than necessary.

  * •

Models for data with more features require more computing time.

  * •

It may happen, that some of the features have no predictive power at all.
Identifying and removing them can save costs.

Selecting the optimal number of features by brute force, i.e. testing all
possible combinations, is in general not feasible since the combinatorial
possibilities grow exponentially with the number of features. Therefore,
several methods for feature selection have been proposed and applied in the
past:

  * •

Using a penalty term $\alpha\sum_{i=1}^{n}\mid w_{i}\mid$ with an adjustable
hyper-parameter $\alpha$ in the cost function of the training algorithm
effectively reduces the number of features. This is called LASSO [25] in the
machine learning context.

  * •

The recursive feature elimination procedure (RFE) starts with all possible
features and then reduces the number of features according to their importance
one by one in a greedy manner. Instead of comparing all combinations of
features, the number of models trained and compared corresponds to the number
of features.

  * •

One may also use unsupervised methods, e.g., principal component analysis
(PCA), to reduce the number of features in the model. However, this is done
independently by the target and can therefore not take into account which
features may have high predictive power. Also, PCA and similar methods
generate a (lesser number of) new features from the original features, rather
than merely selecting. Hence, the original features would still have to be
available when deploying this approach in practice.

### 2.2. QUBO-formulation of the problem

Here, following the suggestion of [20], we propose and explore a different
approach for feature selection: Let $\rho_{i,j}$ be a dependency measure
between column $i$ and column $j$ of the matrix $X$ and let $\rho_{i,Y}$ be
the dependence between column $i$ and the target vector. Later, we will
specify which dependency measures are used.

In addition, let $z_{i}\in\\{0,1\\}$ be a binary variable, indicating if
column $i$ should be selected ($z_{i}=1$) in the model or not ($z_{i}=0$).
Intuitively, features should be selected so that the dependence

$\sum_{i,j=1,i\neq j}^{n}z_{i}z_{j}\mid\rho_{i,j}\mid,$ (1)

between them is minimized, whereas the dependence between the features and the
target vector

$\sum_{i=1}^{n}z_{i}\mid\rho_{i,Y}\mid.$ (2)

is maximized. Note, that we take the absolute values of the dependency
measures since for the dependence measures considered in the following (e.g.
correlation) a high negative dependence is as relevant as a high positive
dependence.

Putting both criteria together, we get the following optimization (here:
minimization) problem:

$h(\boldmath
z)=-\left[\phi\sum_{i=1}^{n}z_{i}\mid\rho_{i,Y}\mid-(1-\phi)\sum_{i,j=1,i\neq
j}^{n}z_{i}z_{j}\mid\rho_{i,j}\mid\right]$ (3)

The parameter $\phi\in[0,1]$ has to be considered as a hyper-parameter and
governs the weighting between the two conditions. By tuning $\phi$ one can
decide which of the two conditions is more important. Using the fact that
$z_{i}$ is a binary variable, i.e. $z_{i}z_{i}=z_{i}$, the above equation can
be transformed to

$h(\boldmath z)=-\left[{\boldmath z^{T}}Q{\boldmath z}\right]$ (4)

and the solution to our problem is given by

${\boldmath z}^{*}=\mathop{\mathrm{arg\,min}}_{z}\left[{h(\boldmath
z)}\right]$ (5)

Such a problem, known as quadratic unconstrained binary optimization (QUBO)
problem can be treated with classical numerical methods but is also suited for
a solution with a quantum computer.

Up to now the dependency measure is not yet determined. The most natural
choice would be to use correlation as a linear measure, ${\rho}^{Correl}$, as
shown in Fig. 1. Alternatively, one could choose the rank correlation.

Figure 1: Illustration of the objective function for lending club data.

Other possible dependency measures are the mutual information $\rho^{MI}$
which has been discussed in [6], the univariate ROC-value $\rho^{ROC}$ and the
Anova F-Statistic $\rho^{Anova}$. All these dependency measures could be used
either between the features or between the features and the target variable.

In the following, for the inter-feature dependence we restrict ourselves to
correlation and the mutual information measure. For the dependence on the
target, we tested all measures. In the following, we will denote the selected
dependence by the tuple $(\rho_{Feature}^{X},\rho_{Target}^{Y})$ where
$X\in\\{Correl,MI\\}$ and $Y\in\\{Correl,MI,ROC,Anova\\}$.

## 3\. Classical feature selection for a small-sized data set

### 3.1. Machine Learning method

In this chapter, a comparison of the described feature selection methods is
performed. To that end, we consider two data sets: a) a subset of Kaggle’s
lending club data [18] and b) the Wisconsin breast cancer data [10]. In the
case of a) we use 8 features and 1000 observations; in the case of b) we have
10 features and 569 observations in the data set. For both data sets we seek a
binary classification model with the "optimal" number of features. To keep the
analysis transparent we use logistic regression as the machine learning model
and use the whole data set for training. We are aware, that this may lead to
overfitting, however, logistic regression, in comparison to more sophisticated
models, proves often to be robust against overfitting (as long as the model
term is restricted to low-order polynomials).

For comparing models with different input features, we use the overall
accuracy (ACC) of the model as well as the area under curve (AUROC) as
performance measures. Accuracy is the proportion of correctly classified cases
by the total number of all data points. AUROC is the area under the ROC curve
which is constructed by plotting the true positive rate (TPR) against 1-FPR,
i.e. false positive rate (FPR). Intuitively AUROC measures the ability of the
model to separate the two classes which is desirable in many use cases.

To find the optimal number of features, we use these two measures as relevant
criteria. Since the accuracy of a model may provide misleading results in case
of imbalanced data sets, we chose the AUROC as the measure to optimize.

For a small number of features $n$, it is possible to determine the optimal
feature selection by brute force, i.e. simply trying all possible combinations
of features and then training the corresponding model. The determined optimal
choice is compared to the results obtained with the RFE- and the LASSO-method.
For LASSO we used cross-validation to find the optimal value of $\alpha$.

In the case of the lending club data, we observed the results shown in Table
1.

Table 1: Lending club data: performance of the logistic regression model for different feature selection approaches. | AUROC | Accuracy | Number of features
---|---|---|---
All features | 0.6651 | 0.83 | 8
Best choice of features | 0.6682 | 0.83 | 6
Selection with RFE | 0.6593 | 0.83 | 1
Selection with LASSO | 0.6681 | 0.83 | 3

While the accuracy is robust against different choices of features, this is
not the case for the AUROC. The best choice of features (obtained with brute
force) can not be obtained via RFE or LASSO: The optimal choice is six out of
eight features; RFE recommends keeping only one feature while LASSO recommends
keeping three features out of eight.

In the case of the breast cancer data, the results are shown in Table 2.

Table 2: Breast cancer data: performance of the logistic regression model for different feature selection approaches. | AUROC | Accuracy | Number of features
---|---|---|---
All features | 0.9844 | 0.9297 | 10
Best choice of features | 0.9863 | 0.9297 | 5
Selection with RFE | 0.9839 | 0.9297 | 6
Selection with LASSO | 0.9845 | 0.9262 | 3

The optimal choice is five out of ten features; RFE recommends to keep six
features and LASSO recommends to keep three features out of eight. Again, the
best choice of features is not obtained with RFE or LASSO.

The results obtained so far should be considered as a benchmark and it has to
be shown, whether feature selection via optimization as described above, is
able to compete with these benchmarks.

### 3.2. Optimization methods

To answer the previous question, we treat the feature selection task as an
optimization problem. Again, due to the low dimensionality of the two
problems, the exact solution of the optimization, c.f. equation (5), can be
obtained by brute force.

However, it is not obvious, which value of $\phi$ should be chosen to obtain
the best solution. To tackle this problem and to explore the dependence of the
results on $\phi$, we solve the optimization problem for different values of
$\phi$ and compare the results to each other. Here, we also consider different
combinations of dependency measures - between the features, but also between
each feature and the target vector.

In Fig. 3 and Fig. 3 the results for the lending club data are shown. Several
findings can be observed in these figures:

  * •

Using mutual information as the dependence measure between features,
irrespective of the measure to the target, leads to the worst AUROC values.

  * •

When using ANOVA as the dependency measure, we observe that the AUROC is
largely independent of $\phi$.

  * •

In most cases, we see that the AUROC depends on $\phi$; the AUROC tends to
increase with increasing values of $\phi$.

(a) Lending club: AUROC as a function of $\phi$ for different dependency
measures (1) (b) Lending club: AUROC as a function of $\phi$ for different
dependency measures (2)

(c) Lending club: AUROC as a function of $\phi$; best obtainable result

Figure 2: Lending club data

(d) Breast cancer: AUROC as a function of $\phi$ for different dependency
measures (1) (e) Breast cancer: AUROC as a function of $\phi$ for different
dependency measures (2)

(f) Breast cancer: AUROC as a function of $\phi$; best obtainable result

Figure 3: Breast cancer data

In all of our experiments, we find the highest value of the AUROC when using
the ROC as the dependency measure between the features and the correlation or
mutual information between each feature and the target. In addition, the $\
phi$ value should be chosen larger than 0.85. Note, $\phi=1$ is the case where
only the dependence between features and target is considered. The
$\phi-$dependence is displayed in Fig. 3. Here, we show for comparison - as a
straight line - the best AUROC of the best model (obtained by brute force),
the AUROC of the RFE- and LASSO models. It can be seen, that RFE is
outperformed by our proposed optimization approach with a suitable choice of
dependency measure. The LASSO result is only slightly worse than the best
possible. The best solution of our optimization approach coincides with the
result, obtained by selecting all features, here denoted by "ALL".

From this analysis, we conclude, that a high value for the AUROC can be
obtained when using the ROC-dependency measure between the features and the
correlation dependency measure between each feature and the target vector,
i.e., the combination (ROC, Correl). Furthermore, we find that a rather high
value of $\phi\approx 0.95$ should be used.

The breast cancer data set leads to somewhat different results, which shows
the well-known dependency of machine learning results on individual data sets.
In Fig. 3, the overall $\phi$-dependence of the AUROC for different dependency
measures is shown. Here, we find the highest value of the AUROC for the
combination (MI, MI), but only for a narrow band of $\phi$-values around 0.75.
The next best combination is (ROC, Correl), as in the case of lending club
data.

Fig. 3 shows the best and the second best constellation - compared to the
brute-force, RFE, and LASSO results. Again, with all selection methods
investigated the best value of AUROC can not be achieved. However, with the
optimization approach (and a suitable choice of parameters) we can beat RFE as
well as LASSO.

To summarize our results so far:

  * •

For both data sets it was not possible to determine the best subset of
features (as obtained by brute-force) with the classical selection methods
(RFE and LASSO) considered.

  * •

Using feature selection via optimization may be superior over these classical
methods but the results are highly dependent on the dependency measures used
and the hyper-parameter $\phi$.

  * •

For both data sets we found better results for higher values of $\phi$ which
indicates that the dependence of the features with the target seems to be more
important than the intra-feature dependence.

## 4\. Feature Selection within the gate-based quantum computer framework

### 4.1. Quantum Algorithms for Optimization

The optimization problem in equation (5) is the starting point for
transferring the problem to a quantum computer in the framework of the gate-
based approach. The binary variables $z_{i}$ are converted to operators with
eigenvalues $\\{1;-1\\}$ and eigenstates $|0\rangle$ and $|1\rangle$. In this
way, by applying the transformation $z_{i}=(1+s_{i})/2$ our problem is mapped
to finding the ground state of the Ising-like model

${\cal H}=\sum_{i,j}J_{i,j}s_{i}s_{j}\leavevmode\nobreak\ ,$ (6)

where the coupling matrix $J_{i,j}$ is derived from the above introduced
matrix $Q$.

For solving these types of problems on a gate-based computer several
algorithms have been proposed. Here, we mainly use the QAOA algorithm, which
has been suggested by Farhi et al. [12].

The key idea of QAOA is to generate the following quantum state
$|\psi_{\vec{\gamma},\vec{\beta}}\rangle$ depending on the parameters
$\vec{\gamma}=(\gamma_{1},\dots,\gamma_{p})$ and
$\vec{\beta}=(\beta_{1},\dots,\beta_{p})$ (with a number of iterations $p$):

$|\psi_{\vec{\gamma},\vec{\beta}}\rangle_{M}=\hat{U}_{M}(\beta_{p})e^{-i\gamma_{p}\hat{F}}\dots\hat{U}_{M}(\beta_{2})e^{-i\gamma_{2}\hat{F}}\hat{U}_{M}(\beta_{1})e^{-i\gamma_{1}\hat{F}}|\psi_{0}\rangle_{M}\leavevmode\nobreak\
,$ (7)

where the initial state $|\psi_{0}\rangle_{M}$ and the operator
$\hat{U}_{M}(\beta)$ depend on the choice of the mixer $M$. Here, we use the
so-called standard mixer with the Pauli X-matrix

$\hat{U}_{\rm
standard}(\beta)=e^{i\beta\sum_{i=1}^{n}\hat{X}_{i}}\leavevmode\nobreak\ .$
(8)

After applying the parameterized gates, all qubits are measured with respect
to the standard basis which leads to a classical bit-string of zeros and ones.
Each bit-string corresponds to a different state and appears with a certain
probability. There are $2^{n}$ different states.

These bit-strings of several thousand "shots" (we use 8192 shots) are plugged
into the objective function eq. (3) to calculate the expectation value. Then,
a classical optimizer is used to update the parameters $\vec{\gamma}$,
$\vec{\beta}$, and in order to minimize the expectation value and thus the
objective function. In other words, to find the ground state of its converted
problem Hamiltonian. A thorough benchmarking study of the QAOA approach,
including a description of the solution method for QAOA, is described in great
detail by Brandhofer et al. [5].

As a measure to gauge the performance of this optimization procedure we use
the approximation ratio:

$r(z_{1},\dots,z_{n})=\frac{h(z_{1},\dots,z_{n})-h_{\rm max}}{h_{\rm
min}-h_{\rm max}}\leavevmode\nobreak\ .$ (9)

Here $h_{\rm max}$ and $h_{\rm min}$ denote the worst and the best bit-string
solution of the problem. If the trained quantum algorithm would always produce
the optimal solution with probability 1, the approximation ratio would be $1$.

To compare the quality of the results, in addition, we use the VQE algorithm
[17] with the ansatz of real-amplitudes [16]. Basically, the ansatz of real-
amplitudes consists of one parametrized R(y)-rotation per qubit and entangling
neighboring qubits with a CNOT-gate. Repeating this procedure several times (=
number of layers) increases the depth of the circuit and thereby the
complexity.

All quantum-based calculations in this paper have been performed with IBM’s
Qiskit framework [23]. Especially, we use out-of-the-box versions of the QAOA-
and VQE-algorithm [15], the state vector- and QASM-simulator, as well as
different physical quantum backends. The formulation of the optimization
problem with IBM’s Qiskit [23] is straightforward. Within this framework,
different solution methods are available for comparison:

  * •

Solution with IBM’s classical optimizer CPLEX [8].

  * •

Small problems can be solved by a direct diagonalization of the problem’s
Hamiltonian.

  * •

Using the state vector simulator allows to perform the quantum mechanical
calculations directly without destroying the quantum state. Again, this is
only feasible for small problem sizes.

  * •

Using a quantum simulator that mimics real hardware, where qubits are
prepared, transformed by several gates, and then measured. From the
measurements, the solution to the problem can be derived. Here, we use the
noise-free QASM simulator within IBM’s framework which is available in the IBM
cloud for sizes up to 32 qubits.

  * •

Using a physical quantum computer in the cloud to perform the calculation on a
real device.

For small problem sizes, as the one considered up to now, all solution
strategies described can be used to benchmark the results. However, the bigger
the problem, the fewer are the range of usable strategies. From a problem size
of approx. 50 qubits, the exact diagonalization, and the state vector and
quantum simulation methods are not available anymore and the only remaining
possibility is using CPLEX (running classical optimization routines under the
hood) or using real quantum hardware.

### 4.2. Results on quantum simulators for different dependency measures and
$\phi$-values

#### 4.2.1. Lending Club Data

As a first step, we compare the results shown in Table 3 for the case of the
lending club data set with different tuples
$(\rho_{Feature}^{X},\rho_{Target}^{Y})$. Interestingly, for all the applied
solution methods and algorithms (CPLEX, exact diagonalization, state vector-
simulator, QASM-simulator, QAOA- and VQE-algorithm) the same optimization
value (given $(\rho_{Feature}^{X}$ and $\rho_{Target}^{Y})$) is obtained.
Furthermore, the calculated accuracy is for all dependency measures equal to
$0.830$. However, as shown in the table above, the number of selected features
and the area under curve (AUROC) differ.

Table 3: Lending club data: performance of the optimization approaches when selecting features for the logistic regression model with different dependency measures. $\rho_{Feature}^{X}$ | $\rho_{Target}^{Y}$ | Optimization value | AUROC | Number of features
---|---|---|---|---
$X=Correl$ | $Y=Correl$ | -0.3542 | 0.6635 | 5
$X=Correl$ | $Y=MI$ | -0,0617 | 0.6405 | 3
$X=Correl$ | $Y=ROC$ | -3.7803 | 0.6651 | 8
$X=Correl$ | $Y=ANOVA$ | -87.8750 | 0.6634 | 7
$X=MI$ | $Y=Correl$ | -0.4117 | 0.6639 | 6
$X=MI$ | $Y=MI$ | -0.0445 | 0.5790 | 1
$X=MI$ | $Y=ROC$ | -4.0131 | 0.6651 | 8
$X=MI$ | $Y=ANOVA$ | -88.0025 | 0.6634 | 7

In Fig. 4, we show the dependency of the approximation ratio on the number of
layers - for QAOA as well as for VQE. For comparison, we also show the
approximation ratio of random search. This value is obtained by drawing 1000
uniform random bit-strings and calculating the mean approximation ratio.

Figure 4: Lending Club data: Approximation ratio obtained with VQE, QAOA and
tuned QAOA for different values of p, $\phi=0.95$ and for the following
dependency measures: (a) $(\rho_{Feature}^{Correl},\rho_{Target}^{Correl})$,
(b) $(\rho_{Feature}^{Correl},\rho_{Target}^{MI})$, (c)
$(\rho_{Feature}^{Correl},\rho_{Target}^{ROC})$, (d)
$(\rho_{Feature}^{Correl},\rho_{Target}^{ANOVA})$, (e)
$(\rho_{Feature}^{MI},\rho_{Target}^{Correl})$, (f)
$(\rho_{Feature}^{MI},\rho_{Target}^{MI})$, (g)
$(\rho_{Feature}^{MI},\rho_{Target}^{ANOVA})$, (h)
$(\rho_{Feature}^{MI},\rho_{Target}^{ROC})$

The approximation ratios obtained with the VQE algorithm are in almost all
cases close to 1 and are in general higher than the results obtained with
QAOA. Further, the QAOA results are distinct: In some cases, the approximation
ratios are - independent of $p$ \- near 1, in other cases, they increase with
increasing $p$. This supports the conjecture raised in [5] that there are easy
and not-so-easy parameter combinations. From a practical point of view, it
would be desirable to have on the one hand a QUBO matrix that leads to a high
value of AUROC (see table 3) and on the other hand, a QUBO matrix that is easy
to handle for the calculation, i.e. it leads to high approximation ratios with
only a few iterations $p$. Taking this into account,
$(\rho_{Feature}^{Correl},\rho_{Target}^{ROC})$ seems to be a good choice for
this problem.

In order to improve the performance of QAOA, we make use of four heuristic
methods ((i) extrapolation, (ii) linear ansatz, (iii) quadratic ansatz, (iv)
adding zero angles) to choose suitable candidates for initial values of the
circuit parameters $\gamma$ and $\beta$ for increasing QAOA layers $p$ as
proposed in [5]. By repeating the optimization runs with different initial
values one can avoid getting stuck in a local minimum. Further, it has been
shown that the required time scales exponentially faster to achieve a similar
performance when using the observed patterns of the optimal parameters instead
of using arbitrary starting values [29]. The QAOA routine as proposed in [5]
and which is used here to solve the feature selection problem is depicted in
Fig. 5.

Figure 5: Flowchart of the QAOA routine for solving the Feature Selection
problem, whereby the particularity is to find suitable candidates for initial
values of the circuit parameters using four heuristic methods described in [5]

In Fig. 6, the approximation ratios for different values of $\phi$ for the
lending club data are given as a function of the number of QAOA layers $p$.
For this data set, it is noticeable that the majority of the dependency
measures (Fig. 6 (a), (b), (d), (e), (g)) a value of $\phi=0$ lead to the
fastest convergence towards the optimal solution, whereas higher values of
$\phi$ (e.g. $0.75\leq\phi<0.95$) lead to a slower convergence. In other
words, we find that the dependence between the features is more significant
than the dependence between features and the target. For the other dependency
measures (Fig. 6 (c), (f), (h)), we find that large values of $\phi$ (e.g.
$0.75\leq\phi<1$) perform better than values of $\phi\leq 0.5$. Consequently,
the weight for the two conditions (intra-feature dependence and dependence
between features and target) depends on the choice of the dependency measure
and it is not useful to be set globally in order to achieve the best possible
performance - it is rather suggested to set a local value for the weighting
parameter depending on the individual dependency measure.

Figure 6: Landing Club data: Approximation ratio for different values of
$\phi$ shown for different values of QAOA layers p evaluated for the following
dependency measures: (a) $(\rho_{Feature}^{Correl},\rho_{Target}^{Correl})$,
(b) $(\rho_{Feature}^{Correl},\rho_{Target}^{MI})$, (c)
$(\rho_{Feature}^{Correl},\rho_{Target}^{ROC})$, (d)
$(\rho_{Feature}^{Correl},\rho_{Target}^{ANOVA})$, (e)
$(\rho_{Feature}^{MI},\rho_{Target}^{Correl})$, (f)
$(\rho_{Feature}^{MI},\rho_{Target}^{MI})$, (g)
$(\rho_{Feature}^{MI},\rho_{Target}^{ANOVA})$, (h)
$(\rho_{Feature}^{MI},\rho_{Target}^{ROC})$

Lastly in this section, we present the results obtained using a physical
quantum backend of IBM for a fixed parameter setup and with the ROC-dependence
between the features and the target, and the linear correlation between all
features. Furthermore, we set $\phi=0.95$.

Using the lending club data, we obtain the same optimal value for the
objective function, $obj=-3.7803$, with the classical optimizer, the direct
diagonalization of the Hamiltonian, the state vector simulator, and with the
quantum simulator. Further, we also obtain the same value for the AUROC, i.e.
$AUROC=0.6651$, c.f. (3). This corresponds to the selection of all features.

When using a physical quantum backend (here: the IBM system Montreal) with
$p=1$, we obtain the same objective function value of $-3.7803$. However, the
approximation ratio drops to 0.59. With $p=5$, the approximation ratio
increases to 0.66. These values are to be compared with 0.52, the
approximation ratio obtained by randomly selecting bit-strings. While these
results reflect the errors of existing real-world QC backends, they show
nevertheless the expected dependence on $p$.

#### 4.2.2. Breast Cancer Data

We now repeat the calculation from the last subsection for the breast cancer
data. The results for the different tuples
$(\rho_{Feature}^{X},\rho_{Target}^{Y})$ are reported in Table 4.

Table 4: Breast cancer data: performance of the optimization approach when selecting features for the logistic regression model with different dependency measures. $\rho_{Feature}^{X}$ | $\rho_{Target}^{Y}$ | Optimization value | AUC | ACC | Number of features
---|---|---|---|---|---
$X=Correl$ | $Y=Correl$ | -0.91611 | 0.9806 | 0.9209 | 3
$X=Correl$ | $Y=MI$ | -0.3300 | 0.9644 | 0.8910 | 1
$X=Correl$ | $Y=ROC$ | -1.7031 | 0.9822 | 0.9192 | 4
$X=Correl$ | $Y=ANOVA$ | -2913.6696 | 0.9839 | 0.9315 | 9
$X=MI$ | $Y=Correl$ | -1.438 | 0.9796 | 0.9262 | 5
$X=MI$ | $Y=MI$ | -0.4139 | 0.9859 | 0.9297 | 3
$X=MI$ | $Y=ROC$ | -2.6366 | 0.9828 | 0.9279 | 6
$X=MI$ | $Y=ANOVA$ | -2914.5083 | 0.9839 | 0.9315 | 9

While the optimization values, given $(\rho_{Feature}^{X},\rho_{Target}^{Y})$)
coincides in general for CPLEX, exact diagonalization, and QAOA (state vector-
simulator and QASM-simulator), we detect that VQE is not always able to find
the same optimization value. Interestingly, this behavior depends on the
number of layers. I.e. in the case
$(\rho_{Feature}^{Correl},\rho_{Target}^{ROC})$) it occurs that only for
$p\leq 3$ is the VQE algorithm able to find the optimal value.

Here, in contrast to the lending club data, we see that the accuracy depends
on the dependency measure and the highest value for the accuracy can not be
obtained together with the highest value of the AUROC.

Figure 7: Breast cancer data: Approximation ratio obtained with VQE, QAOA, and
tuned QAOA for different values of p, $\phi=0.75$ and for the following
dependency measures: (a) $(\rho_{Feature}^{Correl},\rho_{Target}^{Correl})$,
(b) $(\rho_{Feature}^{Correl},\rho_{Target}^{MI})$, (c)
$(\rho_{Feature}^{Correl},\rho_{Target}^{ROC})$, (d)
$(\rho_{Feature}^{Correl},\rho_{Target}^{ANOVA})$, (e)
$(\rho_{Feature}^{MI},\rho_{Target}^{Correl})$, (f)
$(\rho_{Feature}^{MI},\rho_{Target}^{MI})$, (g)
$(\rho_{Feature}^{MI},\rho_{Target}^{ANOVA})$, (h)
$(\rho_{Feature}^{MI},\rho_{Target}^{ROC})$

For comparison, the approximation ratios for the breast cancer data are shown
in Fig. 8. Again, different values of $\phi$ are used and we show the
dependence on the number of QAOA layers. For the dependency measures
$(\rho_{\text{Feature}}^{\text{Correl}},\rho_{\text{Target}}^{\text{Correl}})$
and $(\rho_{\text{Feature}}^{\text{MI}},\rho_{\text{Target}}^{\text{Correl}})$
we find that setting the parameter $\phi$ to $0.75$ leads to the slowest
convergence towards the optimal solution, whereas for the other measures a
value of $\phi=0.95$ is showing this behavior. The best convergence is
achieved in all cases for an intra-feature setting ($\phi=0$) that opens up
the question of what influence the condition between feature and target should
have, while the results show that this condition does not play a role in order
to achieve the best-optimized solution. Compared to the lending club data set,
we find that the approximation ratios for the different values of $\phi$ show
a more continuous behavior in terms of the choice of the dependency measure.

Figure 8: Breast cancer data: Approximation ratio for different values of
$\phi$ shown for different values of QAOA layers p evaluated for the following
dependency measures: (a) $(\rho_{Feature}^{Correl},\rho_{Target}^{Correl})$,
(b) $(\rho_{Feature}^{Correl},\rho_{Target}^{MI})$, (c)
$(\rho_{Feature}^{Correl},\rho_{Target}^{ROC})$, (d)
$(\rho_{Feature}^{Correl},\rho_{Target}^{ANOVA})$, (e)
$(\rho_{Feature}^{MI},\rho_{Target}^{Correl})$, (f)
$(\rho_{Feature}^{MI},\rho_{Target}^{MI})$, (g)
$(\rho_{Feature}^{MI},\rho_{Target}^{ANOVA})$, (h)
$(\rho_{Feature}^{MI},\rho_{Target}^{ROC})$

As a final result, we present the solution for a fixed parameter combination
on a physical quantum backend of IBM. The results are obtained with the mutual
information dependence between the features and also to the target.
Furthermore, we fix $\phi=0.75$. Due to limitations in the use of physical
quantum backends, we restrict ourselves to two cases and solved the problem
with the QAOA algorithms with different values of $p$.

1) $(\rho_{Feature}^{Correl},\rho_{Target}^{Correl})$: As backends, we use the
IBM systems Toronto ($p=1,5$), Auckland ($p=2,3$), and Hanoi ($p=4$). We are
aware of the fact that results from different machines are not directly
comparable. However, due to the availability of system resources, it was not
possible to perform all calculations on the same hardware. As approximation
ratios we get: $p=1:0.8981;p=2:0.8734;p=3:0.8917;p=4:0.8562;p=5:0.8706$. These
values have to be compared with the random approximation ratio of $0.7984$,
which is significantly lower for all values of $p$. The approximation ratio
obtained on the physical backends is in the order of magnitude as obtained
with the QASM simulator.

2) $(\rho_{Feature}^{Correl},\rho_{Target}^{MI})$: Again, we obtained the same
value for the objective function and the same features selected as reported in
Table 4. As backends, we use the IBMQ systems Toronto ($p=1$) and Geneva
($p=2,3,4,5$). As approximation ratios we get:
$p=1:0.7897;p=2:0.8348;p=3:0.7311;p=4:0.7620;p=5:0.7491$. These values have to
be compared with the random approximation ratio of $0.7743$. The approximation
ratio obtained on the physical backends is in the order of magnitude obtained
with the QASM simulator.

In the next section, we will explore the results obtained in the context of a
larger problem set.

## 5\. Extension to 27 features aka qubits

The results obtained so far are encouraging and we extend the calculation to a
larger problem. Here, we use a subset of 27 features of the data considered in
[20].

For this problem size, using a brute-force approach to determine the best
selection of features is not possible anymore. In addition, the state vector-
simulator and the exact diagonalization of the Hamiltonian are infeasible.
Therefore, we compare the following approaches before applying feature
selection via optimization:

  * •

Calculate the AUROC with all features,

  * •

Select the most relevant features via RFE and calculate the AUROC,

  * •

Select the most relevant features via LASSO and calculate the AUROC.

Table 5: Model performance for the data with 27 features. | AUROC | Accuracy | Number of features
---|---|---|---
All features | 0.7913 | 0.764 | 27
Selection with RFE | 0.7904 | 0.757 | 17
Selection with LASSO | 0.7768 | 0.732 | 17

The results are shown in Table 5. Again, these values serve as benchmarks and
have to be compared with the results obtained according to the QUBO scheme. We
translate the feature selection problem to an optimization problem, using
correlation as the relevant dependence measure between the features and
between features and target and we fix $\phi=0.9$.

Then, the following strategies are used to solve the optimization problem:

  * •

Select the features with a greedy algorithm, starting from a random
configuration and then flipping the qubits to improve the objective function

  * •

Solve the optimization problem with a classical optimizer (here: IBM’s CPLEX).

  * •

Solve the problem with a quantum simulator, available in the IBM cloud with
using the QAOA algorithm.

  * •

Solve the optimization problem using the QAOA algorithm with a 27-qubit real
quantum device.

In contrast to the last chapter, where we used a tuned QAOA algorithm, we use
here IBM’s standard implementation: First, the advantages of the tuned
algorithm are only observable at higher values of p. Since we are now looking
at 27 qubits, a useful computation for higher values of p is no longer
possible, so we can no longer take advantage of this variant and choose to use
IBM’s standard implementation instead. Second, when using the IBM algorithms,
the cloud infrastructure and the interplay between classical and quantum
algorithms via the so-called Runtime can be optimally exploited.

Table 6: Results for the data with 27 features, with classical optimization approaches. | Objective | AUROC | Accuracy | Number of features
---|---|---|---|---
Greedy optimization | -0.7232 | 0.7469 | 0.709 | 8
IBM’s CPLEX | -0.7232 | 0.7469 | 0.709 | 8

The results using the greedy algorithm and IBM’s CPLEX optimizer are presented
in Table 6. The calculation with the real quantum device shown below should be
considered a proof of concept. For the calculation of the approximation ratio,
we use the optimal CPLEX result as ground truth.

We use different US devices available in the IBM cloud - Auckland and Kolkata
as well as the German system Ehningen. Results are shown in Table 7.

Table 7: Results for the data with 27 features, with QC optimization approaches. | Objective | AUROC | Accuracy | Number of features
---|---|---|---|---
QUASM-Simulator, p=1 | -0.6821 | 0.7412 | 0.731 | 7
QUASM-Simulator, p=2 | -0.6974 | 0.718 | 0.707 | 7
QUASM-Simulator, p=3 | -0.70421 | 0.7435 | 0.707 | 7
Real Device (Auckland), p=1 | -0.6097 | 0.7323 | 0.722 | 5
Real Device (Auckland), p=2 | -0.6218 | 0.7444 | 0.709 | 8
Real Device (Auckland), p=3 | -0.6348 | 0.7469 | 0.7225 | 6
Real Device (Kolkata), p=4 | -0.6199 | 0.7482 | 0.72 | 9
Real Device (Ehningen), p=1 | -0.6073 | 0.7374 | 0.702 | 7
Real Device (Ehningen), p=2 | -0.5982 | 0.7304 | 0.695 | 7
Real Device (Ehningen), p=3 | -0.3543 | 0.7452 | 0.712 | 7

Again, we also calculate the random approximation ratio, calculated with a
sample of 1000 bit-strings.

Due to the fact that now we do not necessarily get the same optimal objective
value compared with CPLEX, we use two different prescriptions:

  * •

Approximation ratio rel. CPLEX: We take the best and the worst objective value
out of CPLEX as ground truth and calculate the approximation ratios relative
to these values.

  * •

Approximation ratio rel. Sim.: From the outcome of the simulations (either
random sampling as a baseline or the quantum simulations) we take the best and
the worst value of the objective function and use those in the formula for the
approximation ratio. This number describes the spread of the different
simulation results.

The respective values are shown in Table 8.

Table 8: Approximation ratios for the data with 27 features, with QC optimization approaches. | Approx. ratio rel. CPLEX | Approx. ratio rel. Sim.
---|---|---
Random Sampling | 0.7314 | 0.6899
Real Device (Auckland), p=1 | 0.8351 | 0.7086
Real Device (Auckland), p=2 | 0.8189 | 0.7466
Real Device (Auckland), p=3 | 0.8274 | 0.7159
Real Device (Kolkata), p=4 | 0.8271 | 0.7272
Real Device (Ehningen), p=1 | 0.8356 | 0.7185
Real Device (Ehningen), p=2 | 0.8734 | 0.6941
Real Device (Ehningen), p=3 | 0.8541 | 0.7299

The results obtained with the physical quantum devices are surprisingly good,
considering the existing error rates. The minimum objective function as
determined with CPLEX is almost achieved in most cases. The approximation
ratios are well above the random sampling result and we see a slight
improvement when increasing the value of $p$.

### 5.1. Comparison to classical/metaheuristic optimization algorithms

In order to gain more insight into the details of this optimization problem,
we perform another experiment. Here, we are interested in a comparison between
classical, metaheuristic algorithms, such as evolutionary algorithms (EA) or
estimation of distribution algorithms (EDA) [11, 3, 14, 4].

One core element of such algorithms are stochastic components, e.g., the
random mutation of solutions in an EA or the sampling from a distribution in
an EDA. Similarly to non-quantum metaheuristics, QAOA also has random
components (potentially as part of the classical update step, as well as
stochastic elements in the physical quantum circuit). As such, a comparison to
(mostly) deterministic solvers such as CPLEX is less intuitive. Moreover,
globally optimal solutions might not be necessary (and in fact not guaranteed
to be found within a limited time frame), and hence it may be of interest to
investigate how intermediate (non-optimal) results from a classical stochastic
algorithm compare to those of QAOA.

Another interesting aspect of this comparison is that the metaheuristics we
consider are not limited to quadratic (or polynomial) objective functions:
they are also applicable to arbitrary non-linear optimization problems (of
course, with the respective consequences for performance).

### 5.2. Classical, stochastic optimization algorithms

We tested four metaheuristics: A simple evolutionary algorithm (EA), a self-
adaptive evolutionary algorithm (SAEA), a univariate estimation-of-
distribution algorithm (UEDA), and a discrete version of the covariance matrix
adaption evolution strategy (DCMA).

#### 5.2.1. Evolutionary algorithm

Evolutionary algorithms (EA) iteratively generate new candidate solutions
("offspring") from existing candidate solutions ("parents") via random
variation operators ("mutation" and "recombination"). Well-performing
candidates influence candidates in the next iteration. For a more in-depth
introduction see, e.g., [11].

The simple EA employed in this experiment uses a fixed mutation and
recombination rate, considering a population of $\mu$ parents. Initially,
candidate solutions are generated by first sampling the number $n_{bit}$
uniformly from ${1,2,...,27}$, then uniformly and randomly generating a bit-
string with exactly $n_{bit}$ ones.

In each iteration, new candidate solutions are generated via a mutation
operator and a recombination operator. The following operators are used:

  * •

Mutation: bit-flip; inverting a random number of bits of an existing candidate
solution.

  * •

Mutation: block-inversion; inverting a random block (sequence)
${i,i+1,...j-1,j}$ of bits.

  * •

Mutation: cycle; cyclically shifting bit-values to the right or left.

  * •

Recombination: one-point crossover; all bits up to bit $i$ are taken from one
parent, the remaining bits from the second parent.

  * •

Recombination: two-point crossover; selecting a random sequence of bits
${i,i+1,...j-1,j}$ from one parent, the remaining bits from the second parent.

  * •

Recombination: uniform crossover; selecting bits uniform randomly from both
parents.

The mutation rate $r_{m}$ determines how many bits are changed with each
mutation operator (i.e., number of bits flipped, size of the inverted block,
step size of the cyclical shift). We will investigate four parameters for
tuning: population size $\mu$, mutation operator choice for recombination
$o_{r}$ and mutation $o_{m}$, and the mutation rate $r_{m}$. The remaining
parameters of the EA are left at default values. For performance reasons, an
archive of all candidate solutions is not retained. The employed
implementation is available via the R-package CEGO [28].

#### 5.2.2. Self-adaptive evolutionary algorithm

One issue of the ’simple’ EA described above is that it requires setting fixed
algorithm parameters. The performance of the EA is rather sensitive to some of
these parameters (e.g., population size, operator choices, and mutation rate).
Hence, these parameters may require tuning or else, may lead to unsatisfactory
performance.

The self-adaptive EA (SAEA) tries to partially alleviate this issue by
attaching several parameters $\theta$ to each solution candidate (i.e., to
each individual in a population). Our SAEA variant follows a similar approach
as described for the Mixed Integer Evolution Strategy [19].

Subsequently, the parameters $\theta$ are evolved alongside the actual binary
decision variables $z$: each candidate solution is then composed as
$z^{\prime}=\\{z,\theta\\}$. Thus, if some solution candidate is successful
(i.e., has a good objective function value), the associated parameters
$\theta$ (e.g., the mutation rate) that were used to generate it will be more
likely to be used in subsequent iterations as well.

Note, that some algorithm parameters will necessarily be fixed (at least
during a single algorithm run), and not attached to each candidate solution.
For example, we have to specify a set of operators, ranges of parameters, and
learning rates to configure the self-adaptive procedure itself. The underlying
motivation is that the algorithm should be less sensitive to the parameters
that control the self-adaptive procedure (e.g., learning rate) when compared
to the sensitivity to changes in the parameters controlled by that procedure
(e.g., mutation rate).

In detail, the SAEA will self-adapt the following three parameters: the
mutation rate $\theta_{1}=r_{m}\in[1/n,1]$, the mutation operator choice
$\theta_{2}=o_{m}$ (bit-flip, block-inversion, cycle), and the recombination
operator choice $\theta_{3}=o_{r}$ (1-point, 2-point, uniform).

The choices of operators for mutation and recombination are themselves also
mutated, choosing randomly with probability $p_{r}$ from the parents.
Similarly, the mutation rate is mutated in each iteration, via
$\theta_{1}^{*}=\theta_{1}e^{\tau\epsilon}$. Here, $\tau$ is a learning rate
(hyperparameter) and $\epsilon$ is a normal-distributed random sample with
zero mean and unit variance. Intermediate crossover is used to recombine
$\theta_{1}^{*}=f(\theta_{1}^{*(1)},\theta_{1}^{*(2)})$ values from parent
solutions $\\{z^{\prime(1)},z^{\prime(2)}\\}$.

Like the simple EA, the SAEA implementation is available via the R-package
CEGO [28], and uses the same initialization strategy.

#### 5.2.3. Univariate estimation-of-distribution algorithm

The third algorithm we employ is an estimation of distribution algorithm
(EDA). EDAs essentially model the set of candidate solutions with a
distribution, iteratively sampling from that distribution to generate new
candidates, then updating the distribution parameters based on the evaluated
candidates.

Here, we use a separate univariate binary distribution to model each solution
bit $z_{i}$. Specifically, our univariate EDA (UEDA) is a variant of the
population-based incremental learning algorithm (PBIL) [1]:

  * •

Each bit $z_{i}$ is represented by an independent probability $p_{z_{i}}$.

  * •

New candidate solutions are sampled based on these probabilities. We refer to
the number of generated candidate solutions (bit-strings) as the population
size $\mu$.

  * •

Each candidate solution is evaluated with the objective function.

  * •

The better 50% of the candidate solutions are selected.

  * •

Based on that selection, new probabilities $p^{*}_{z_{i}}$ are estimated (via
the mean value for each bit, over all selected candidates).

  * •

The probabilities for the next iteration are set via a ’learning’ update
$p_{z_{i}}=p_{z_{i}}(1-\tau)+p^{*}_{z_{i}}\tau$.

  * •

To avoid collapsing to $p_{z_{i}}=0$ or $p_{z_{i}}=1$, the probabilities are
limited to $p_{z_{i}}\in[\frac{1}{n},1-\frac{1}{n}]$.

An implementation of this simple UEDA was written by the authors in R. It uses
the same initialization strategy as the EA.

#### 5.2.4. Discrete covariance matrix adaption evolution strategy (DCMA)

Finally, we also employ a variant of the covariance matrix adaption evolution
strategy (CMAE-ES), which can handle discrete (including binary) variables by
introducing a lower bound on the marginal probabilities [13]. In essence, this
algorithm can also be viewed as an EDA. But in contrast to the UEDA,
correlations between variables $z_{i}$ are taken into account, as the internal
distribution is based on an evolving covariance matrix.

We use the implementation of this margin-based CMA-ES provided by the cmaes
python library [24] and refer to this algorithm as discrete covariance matrix
adaption (DCMA). Besides the population size $\mu$, this algorithm is also
impacted by an initial standard deviation $\sigma_{\text{init}}$, as well as
the margin parameter $\alpha_{m}$.

#### 5.2.5. Random and Greedy Search

Finally, we employ two simple algorithms as baselines to compare against: a
simple greedy local search with restarts (GRS) and random search (RS). The
former is at its core a purely exploitative, local search, while the latter is
purely explorative, global search.

The local search procedure of the GRS starts at a random initial solution
(determined in the same way as the initial solutions in the EA). Neighboring
candidate solutions are produced by bit-flips, in a random order. Once a
neighbor improves on the current solution, it immediately becomes the current
solution. If no neighbor improves on the current solution, the procedure is
restarted with a new initial solution. We do not tune this procedure for the
sake of simplicity. Note that we could in principle make this procedure more
complex, e.g., by redefining neighboring solutions to mean something else than
single bit-flips, by introducing a more clever way to generate new initial
solutions for restarts, or by hybridizing it with the EA or similar
algorithms.

The random search simply tries out new candidate solutions at random, where
each solution is generated in the same way as random initial solutions in the
EAs.

### 5.3. Experiment

The whole data set is randomly split into two subsets. Each subset is used to
generate a corresponding instance of the QUBO objective function, using
correlation as the respective dependency measure. One instance will be used
for tuning, another for validating / comparing the (tuned) algorithms.

For the quantum algorithm QAOA, we run the algorithm on real QC devices
(mainly Geneva and Montreal) nine times for $p=3$, solving only the validation
instance. Similarly, CPLEX is run on the validation instance of the objective
function (and will be used to provide a baseline, as it solves this problem
exactly).

#### 5.3.1. Algorithm Tuning

Tuning is performed on the first of the two instances of our objective
function. The goal is to provide just a very rough pre-configuration for the
meta-heuristics, to avoid detrimental performance due to a poor algorithm
configuration and allow for a reasonably fair comparison between the
metaheuristics.

Therefore, the employed tuning method is uniform random search: For each
algorithm, 100 algorithm configurations are sampled uniformly at random. For
each algorithm configuration, 10 independent algorithm runs are performed.
Each run uses a maximum of 2000 objective function evaluations. Afterward, the
configuration with the best median performance over the 10 runs is returned as
the ’tuned’ configuration. Here, ’performance’ means the median sum (overall
objective function evaluations (OFE)) of the cumulative minimum objective
value recorded by the tuning procedure. We tune two to four ‘main’ parameters
of each algorithm, as specified in Table 9 (omitting the more simple
baselines, and omitting QAOA due to computing constraints).

Table 9: Tuning: bounds and results for each parameter. $b_{t}$ is the budget, i.e., the maximum number of evaluations per run during tuning. Mutation operator codes are 1: bit-flip, 2: block-inversion, 3: cycle. Recombination operator codes are 1: one-point, 2: two-point, 3: uniform. Where values are specified with exponents, random search applies uniformly in the exponents (instead of uniform in terms of actual values). The column ’performance’ lists the median sum of the cumulative minimum objective value recorded by the tuning procedure. In brackets, the standard deviation over all random configurations is given (not the standard deviation of the best performance only). algorithm | parameter | range | tuned | performance
---|---|---|---|---
|  | of values | value | (sd-all)
EA | population size $\mu$ | $[4,200]$ | 52 | -565
| mutation rate $r_{m}$ | $[0,1]$ | 0.042 | (73.2)
| mutation operator $o_{m}$ | $\\{1,2,3\\}$ | 1 |
| recombination operator $o_{r}$ | $\\{1,2,3\\}$ | 3 |
SAEA | population size $\mu$ | $[4,200]$ | 14 | -593
| learning rate $\tau$ | $[10^{-4},10^{0}]$ | $10^{-1.32}$ | (25.6)
| probability $p_{r}$ | $[0,1]$ | 0.30 |
UEDA | population size $\mu$ | $[4,200]$ | 15 | -586
| learning rate $\tau$ | $[0,1]$ | 0.95 | (66.5)
DCMA | population size $\mu$ | $[4,200]$ | 19 | -576
| initial stand. dev. $\sigma_{\text{init}}$ | $[10^{-4},10^{4}]$ | $10^{-2.03}$ | (59.0)
| margin $\alpha_{m}$ | $[(\mu n)^{-1.5},(\mu n)^{-0.5}]$ | $(\mu n)^{-1.23}$ |

The tuning results in Table 9 indicate that SAEA performs best on the
respective tuning instance of the objective function, followed by UEDA. The
simple EA performs worst during tuning. In addition, the low standard
deviation implies that SAEA is less sensitive to the tuned parameters, which
may reasonably be attributed to the self-adaptive mechanism.

In terms of tuned parameters, the mutation rate for the EA is close to $1/n$,
and the respective operator is bit-flip. This suggests that disturbing
solutions as little as possible are preferable. The chosen recombination
operator is a uniform crossover. This makes sense, as there is no obvious
structure in the order of bits $z_{i}$ (i.e., the underlying data features,
whose order is arbitrary) which would have to be preserved by using 1- or
2-point crossover. The population size is the largest in the tested set of
algorithms, which may be a reasonable measure to account for the relatively
static behavior of the EA.

The learning rate $\tau=0.048$ of the SAEA is set rather low, resulting in an
average change of the mutation rate at only 3.8 % per iteration. Considering
that the initial mutation rate is at $1/n$, this means that it will probably
remain fairly low throughout the complete run. Comparatively, the probability
for changes of the operators is set rather high, leading to more frequent
changes of the operator (which may either indicate that these operators have
to change dynamically during the optimization run, or else, that the algorithm
is not that sensitive to frequent changes of the operators).

The UEDA receives a rather high learning rate. This means, new probabilities
$p^{*}_{z_{i}}$ estimated in an iteration contribute more to the respective
probabilities in the next iteration than the probabilities employed in the
current iteration of the UEDA. Seemingly, the UEDA requires changing the
probabilities rather dynamically throughout each algorithm run.

Finally, the DCMA receives a population size that is slightly larger than the
default at 13. The initial standard deviation is set rather low, and the
margin parameter is set close to the suggested default $(\mu n)^{-1}$.

Only the more complex non-QC metaheuristics are tuned. The QAOA algorithm is
not tuned before running it on the physical QC devices, since the
accessibility of the system does not allow for the respective larger number of
runs required for tuning. Also, the metaheuristics, unlike QAOA, are not
specifically developed for quadratic (or polynomial) objective functions and
are hence, by default, less specialized to this problem class.

#### 5.3.2. Validation experiment

Validation is performed on the second of the two instances of the objective
function. The goal of this is an unbiased comparison, to avoid a potential
overfit of algorithm parameters to the tuning instance, although this is
relatively unlikely given the rather rough tuning procedure.

For each algorithm, a maximum of 5000 objective function evaluations are
allowed per independent run. 20 independent runs are performed for each
metaheuristic algorithm (each run with different initial seeds).

Due to limitations in terms of computing access, QAOA does only 9 independent
runs on the validation instance, with 58 to 79 iterations. Note, that each run
takes three to four hours to complete (time-in-queue not included). In
comparison, a run of any of the metaheuristics takes 2 to 6 seconds, depending
on hardware details (using no parallel computing resources). CPLEX usually
requires less than one second but makes use of multiple CPU cores, if
available.

#### 5.3.3. Validation result

Figure 9: Top left: best found objective value after 100 objective function
evaluations (OFE) per run. Top right: same after 5000 OFE. Note: QAOA is best
overall, not counting OFE, hence the same in top left and top right. Bottom
left: approximation ratio (of the last 50 of 5000 OFE in case of the non-
quantum algorithms). Bottom right: raw objective function values of the last
50 OFE.

Figure 9 shows some first results from the validation experiment. The upper
two boxplots show the best solution founds (aggregated over the independent
runs of each algorithm), early after just 100 classical objective function
evaluations (OFE) and finally after 5000 OFE. Note that the numbers shown for
QAOA are in each case recorded at the end of a terminated run. Here, we see
that while QAOA results are better than early results from each classical
algorithm, the classical algorithms eventually find better results, given that
enough OFE is invested.

Between the classical algorithms, GRS performs best in terms of early
performance. SAEA and GRS are the only algorithms to attain the global optimum
in all runs within 5000 OFE, with GRS converging considerably faster. This
clear advantage of GRS is also visible in Fig. 10, and may indicate that the
problem is ‘easy’ to solve at least from the perspective of local search: only
a few repeated local searches are sufficient to consistently determine the
optimum.

The lower left of Fig. 9 shows the approximation ratio (AE) after 5000 OFE.
Note, that unlike the QAOA, AR values of the classical algorithms simply
consider the last 50 evaluations of each run, weighting each observed
candidate solution equally. The EA receives AR values comparable to QAOA, all
other algorithms receive higher values. But note, that AR is not necessarily a
comprehensive quality measure for the classical algorithms, since good
performance (in terms of the best solution found) can be achieved by all four
tested algorithms, despite significant differences in terms of AR. The reason
here is that while final performances are equal, the variation of results
differs a lot between the four algorithms, as the bottom right plot shows.
Here, the EA and SAEA variants show much larger variations than the UEDA and
DCMA. This is not (necessarily) a negative observation for EA and SAEA. In
fact, larger variation within the sampled candidate solutions may imply a
better ability to escape local optima.

Figure 10: Optimality gap of the metaheuristics over the number of OFE (with
log-scaled axes). Colored thick line is the median of each metaheuristic. The
colored ribbon shows 5 and 95 percentiles. Black dashed lines show 5, 50, and
95 percentiles for QAOA. The thin red line at the bottom shows the global
optimum value of the objective function (determined via CPLEX). The optimality
gap is the difference of the best found objective function value to the
optimum. The optimum is rounded down on the third digit to avoid infinite
values. Hence, the global optimum is just below the 0.001 mark. Figure 11:
Optimality gap of the metaheuristics and QAOA over the number of OFE (log-
scaled axes), similar to Figure 10, for the first 1000 OFE. Since QAOA does
not consume OFE, it is instead reported against the number of QAOA iterations,
where each QAOA iteration is (somewhat arbitrarily) valued at 86 OFE, to
enable a comparison to the metaheuristics (i.e., scaled up to 5000 OFE). The
optimality gap is the difference of the best found objective function value to
the optimum. The optimum is rounded down on the third digit to avoid infinite
values. Hence, the global optimum is just below the 0.001 mark.

As seen in Fig. 10, the metaheuristics require a few hundred of OFE to meet
the performance of the final QAOA results. Since QAOA does not consume OFE, it
is instead reported against the number of QAOA iterations in Fig. 11, where
each QAOA iteration is (somewhat arbitrarily) valued at 86 OFE, to enable a
visual comparison to the metaheuristic solvers. The corresponding plot shows
that QAOA itself starts with better objective function values than the
competing algorithms, but after some initial gains, progresses only a little.
Only the simple RS is consistently outperformed by QAOA.

While QAOA (when run on actual QC hardware) achieves decent results fairly
early, it is not yet able to compete with non-QC optimization algorithms. To
determine whether this may change eventually will require a more in-depth
investigation, that would have to make strong assumptions about, among others,
how costs/effort of QC and non-QC algorithms are to be compared reasonably,
and how they scale with aspects such as problem dimension or available
parallel computing power.

## 6\. Discussion

In this paper, the problem of selecting appropriate features for a supervised
learning algorithm has been transformed to an optimization problem. This
optimization problem has been solved with classical numerical methods as well
as with quantum computing.

During the transformation of the feature selection problem, several
ambiguities arise: One may choose between different dependency measures and
choose whether the dependence between the features or the dependence between
each feature and the target is more important. While the choice of the
dependency measure depends on the data set considered, we found that in
general a high value of $\phi$ is favored, which means, that the dependence
between the target and the features is more important.

When comparing the performance of the optimization method (with a proper
choice of the dependency measure and $\phi$), this method can compete with
established methods like RFE or LASSO. However, all feature selection methods
seem to be not able to find the global optimum, which could be obtained by
brute force, for our small data sets.

After confirming that the optimization method for feature selection works in
principle, we explored different methods to solve the resulting quadratic
binary optimization problem: A simple greedy search algorithm, a well-
established commercial optimizer (IBM CPLEX), and a cutting-edge method of
quantum computing.

For the quantum computing solution, we used the QAOA algorithm within the
gate-based approach and used quantum simulators, as well as real quantum
hardware available on the IBM quantum cloud. Whereas the solution of the
small-size data sets merely served to demonstrate the feasibility, the
solution of the 27-feature problem is undoubtedly more interesting. Although
we were limited to $k<=3$ within QAOA, on real hardware we obtained
approximation ratios above 0.80; well above the random sampling result. With
the number of qubits available at the moment, the classical and meta-
heuristical methods are not outperformed by quantum computing. An open
question that will be addressed in future work is the scaling of both
approaches, i.e. if and for which number of qubits there may be a quantum
advantage.

Acknowledgements This work is funded by the Ministry of Economic Affairs,
Labour and Tourism Baden-Württemberg in the frame of the Competence Center
Quantum Computing Baden-Württemberg (project ‘QORA II’).

Useful discussions and the careful reading of the manuscript by PD Dr. Thomas
Wellens, IAF, are greatly acknowledged.

Data availability The data that support the findings of this study are
available upon reasonable request.

Declarations The authors have no relevant financial or non-financial interests
to disclose.

## References

  * [1] Shummet Baluja.  Population-based incremental learning: A method for integrating genetic search based function optimization and competitive learning.  Technical report, Carnegie Mellon University, Pittsburgh, PA, USA, 1994\.
  * [2] Christopher M. Bishop.  Pattern Recognition and Machine Learning.  Springer, New York, NY, 2006.
  * [3] Christian Blum and Andrea Roli.  Metaheuristics in combinatorial optimization: Overview and conceptual comparison.  ACM Computing Surveys, 35(3):268–308, September 2003.
  * [4] Ilhem Boussaid, Julien Lepagnot, and Patrick Siarry.  A survey on optimization metaheuristics.  Information Sciences, 237:82–117, July 2013.
  * [5] Sebastian Brandhofer, Daniel Braun, Vanessa Dehn, Gerhard Hellstern, Matthias Hüls, Yanjun Ji, Ilia Polian, Amandeep Singh Bhatia, and Thomas Wellens.  Benchmarking the performance of portfolio optimization with QAOA.  Quantum Information Processing, 22(1), December 2022.
  * [6] Gavin Brown, Adam Pocock, Ming-Jie Zhao, and Mikel Luj.  Conditional likelihood maximisation: A unifying framework for information theoretic feature selection.  Journal ofMachine Learning Research, 2012.
  * [7] Girish Chandrashekar and Ferat Sahin.  A survey on feature selection methods.  Computers & Electrical Engineering, 40(1):16–28, 2014.  40th-year commemorative issue.
  * [8] IBM ILOG Cplex.  V12. 1: User’s manual for cplex.  International Business Machines Corporation, 46(53):157, 2009.
  * [9] Maurizio Ferrari Dacrema, Fabio Moroni, Riccardo Nembrini, Nicola Ferro, Guglielmo Faggioli, and Paolo Cremonesi.  Towards feature selection for ranking and classification exploiting quantum annealers.  Proceedings of the 45th International ACM SIGIR Conference on Research and Development in Information Retrieval, 2022.
  * [10] Dheeru Dua and Casey Graff.  UCI machine learning repository, 2017.
  * [11] A. E. Eiben and J. E. Smith.  Introduction to Evolutionary Computing.  Natural Computing Series. Springer, Berlin, 2nd edition, 2015.
  * [12] Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.  A quantum approximate optimization algorithm, 2014.  Available at https://arxiv.org/abs/2207.10555.
  * [13] Ryoki Hamano, Shota Saito, Masahiro Nomura, and Shinichi Shirakawa.  Cma-es with margin: Lower-bounding marginal probability for mixed-integer black-box optimization.  In Proceedings of the Genetic and Evolutionary Computation Conference GECCO’22, page 639–647, New York, NY, USA, 2022. Association for Computing Machinery.
  * [14] Mark Hauschild and Martin Pelikan.  An introduction and survey of estimation of distribution algorithms.  Swarm and Evolutionary Computation, 1(3):111–128, September 2011\.
  * [15] IBM.  Qaoa - qiskit 0.39 documentation, 2022.  Available at https://qiskit.org/documentation/stubs/qiskit.algorithms.QAOA.html.
  * [16] IBM.  Real amplitudes - qiskit 0.39 documentation, 2022.  Available at https://qiskit.org/documentation/stubs/qiskit.circuit.library.RealAmplitudes.html.
  * [17] IBM.  ’vqe - qiskit 0.39 documentation’, 2022.  Available at https://qiskit.org/documentation/stubs/qiskit.algorithms.VQE.html.
  * [18] Kaggle.  All lending club loan data, 2022.  Available at https://www.kaggle.com/datasets/wordsforthewise/lending-club.
  * [19] Rui Li, Michael T.M. Emmerich, Jeroen Eggermont, Thomas Bäck, M. Schütz, J. Dijkstra, and J.H.C. Reiber.  Mixed integer evolution strategies for parameter optimization.  Evolutionary Computation, 21(1):29–64, March 2013.
  * [20] Andrew Milne, Maxwell Rounds, and Phil Goddard.  Optimal feature selection in credit scoring and classification using a quantum annealer.  Technical report, 1QBit, 2017.
  * [21] Sascha Mücke, Raoul Heese, Sabine Müller, Moritz Wolter, and Nico Piatkowski.  Feature selection on quantum computers.  Quantum Machine Intelligence, 5(1), February 2023.
  * [22] Riccardo Nembrini, Maurizio Ferrari Dacrema, and Paolo Cremonesi.  Feature selection for recommender systems with quantum computing.  Entropy (Basel, Switzerland), 23(8), 2021.
  * [23] Qiskit contributors.  Qiskit: An open-source framework for quantum computing, 2023.
  * [24] Masashi Shibata.  cmaes: Lightweight Covariance Matrix Adaptation Evolution Strategy (CMA-ES) implementation., 2022.  Python library version 0.9.1. Available at https://github.com/CyberAgentAILab/cmaes.
  * [25] Robert Tibshirani.  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society. Series B (Methodological), 58(1):267–288, 1996.
  * [26] Gloria Turati, Maurizio Ferrari Dacrema, and Paolo Cremonesi.  Feature selection for classification with QAOA.  In 2022 IEEE International Conference on Quantum Computing and Engineering (QCE). IEEE, sep 2022.
  * [27] Hemanth Venkateswara, Prasanth Lade, Binbin Lin, Jieping Ye, and Sethuraman Panchanathan.  Efficient approximate solutions to mutual information based global feature selection.  In ICDM’15: Proceedings of the 2015 IEEE International Conference on Data Mining (ICDM), pages 1009–1014, Atlantic City, NJ, USA, 2015\. IEEE.
  * [28] Martin Zaefferer.  CEGO: Combinatorial Efficient Global Optimization in R, 2021.  R package version 2.4.2. Available at https://CRAN.R-project.org/package=CEGO.
  * [29] Leo Zhou, Sheng-Tao Wang, Soonwon Choi, Hannes Pichler, and Mikhail D. Lukin.  Quantum approximate optimization algorithm: Performance, mechanism, and implementation on near-term devices.  Phys. Rev. X, 10:021067, 2020.
  * [30] Christa Zoufal, Ryan V. Mishmash, Nitin Sharma, Niraj Kumar, Aashish Sheshadri, Amol Deshmukh, Noelle Ibrahim, Julien Gacon, and Stefan Woerner.  Variational quantum algorithm for unconstrained black box binary optimization: Application to feature selection.  Quantum, 7:909, January 2023.