# Convolutional Sparse Support Estimator Based Covid-19 Recognition from X-ray
Images

Mehmet Yamaç  Tampere University, Faculty of Information Technology and
Communication Sciences, Tampere, Finland  Huawei Technologies Oy (Finland) Co.
Ltd, Tampere, Finland  Mete Ahishali  Tampere University, Faculty of
Information Technology and Communication Sciences, Tampere, Finland  Aysen
Degerli  Tampere University, Faculty of Information Technology and
Communication Sciences, Tampere, Finland  Serkan Kiranyaz  Department of
Electrical Engineering, Qatar University, Qatar  Muhammad E. H. Chowdhury
Department of Electrical Engineering, Qatar University, Qatar  Moncef Gabbouj
Tampere University, Faculty of Information Technology and Communication
Sciences, Tampere, Finland

###### Abstract

Coronavirus disease (Covid-19) has been the main agenda of the whole world
since it came in sight in December, 2019. It has already caused thousands of
causalities and infected several millions worldwide. Any technological tool
that can be provided to healthcare practitioners to save time, effort, and
possibly lives has crucial importance. The main tools practitioners currently
use to diagnose Covid-19 are Reverse transcription-polymerase chain reaction
(RT-PCR) and Computed Tomography (CT), which require significant time,
resources and acknowledged experts. X-ray imaging is a common and easily
accessible tool that has great potential for Covid-19 diagnosis. In this
study, we propose a novel approach for Covid-19 recognition from chest X-ray
images. Despite the importance of the problem recent studies in this domain
produced not so satisfactory results due to the limited datasets available for
training. Recall that Deep Learning techniques can generally provide state-of-
the-art performance in many classification tasks when trained properly over
large datasets, such data scarcity can be a crucial obstacle when using them
for Covid-19 detection. Alternative approaches such as representation-based
classification (collaborative or sparse representation) might provide
satisfactory performance with limited size datasets, but they generally fall
short in performance or speed compared to Machine Learning methods. To address
this deficiency, Convolution Support Estimation Network (CSEN) has recently
been proposed as a bridge between model-based and Deep Learning approaches by
providing a non-iterative real-time mapping from query sample to ideally
sparse representation coefficient’ support, which is critical information for
class decision in representation based techniques.

Main premises of this study can be summarized as follows: (i) a benchmark
X-ray dataset, namely QaTa-Cov19, containing over 6200 X-ray images is
created. Up to date, this is the largest dataset covering 462 X-ray images
from Covid-19 patients along with three other classes; bacterial pneumonia,
viral pneumonia, and normal. (ii) In such a scarce and imbalanced dataset the
proposed CSEN based classification scheme equipped with feature extraction
from a state-of-the-art deep neural network solution for X-ray images,
CheXNet, achieves over $98\%$ sensitivity and over $95\%$ specificity for
Covid-19 recognition directly from raw X-ray images without any pre- or post-
processing. (iii) Having such an elegant Covid-19 assistive diagnosis
performance, this study further provides solid evidence that Covid-19 induces
a unique pattern in X-rays that can be discriminated with a high accuracy.

###### Index Terms:

Covid-19 Recognition, SARS-CoV-2 virus, Transfer Learning, Representation
based Classification.

##  I Introduction

Coronavirus disease 2019 (Covid-19) has been declared as a pandemic by the
World Health Organization (WHO) two months after its first appearance in
December, 2019 in Wuhan, China. It has infected more than 3 million people,
caused thousands of causalities and has so far paralyzed the mobility all
around the World. The spreading rate of Covid-19 is so high that the number of
cases is expected to be doubled every three days if the social distancing is
not strictly observed to slow this accretion [1]. Roughly around half of
Covid-19 positive patients exhibit also a comorbidity [2], making difficult to
differentiate Covid-19 from other lung diseases. Automated and accurate
Covid-19 diagnosis is critical for both saving lives and preventing its rapid
spread in the community. Currently, RT-PCR (Reverse transcription polymerase
chain reaction) and CT (computed tomography) are the common diagnosis
techniques used today. RT-PCR results are ready at the earliest 24 hours for
critical cases and generally take several days to conclude a decision [3]. CT
may be an alternative at initial presentation; however, it is expensive and
not easily accessible [4]. The most common tool that medical experts use for
both diagnostic and monitoring the course of the disease is X-ray imaging.
Compared to RT-PCR or CT test, having an X-ray image is an extremely low cost
and a fast process, usually taking only few seconds. Recently, WHO reported
that even RT-PCR may give false results in Covid-19 cases due to several
reasons such as poor quality specimen from the patient, inappropriate
processing of the specimen, taking the specimen at an early or late stage of
the disease [5]. For this reason, X-ray imaging has a great potential to be an
alternative technological tool to be used along with the other tests for an
accurate diagnosis.

Accordingly, there are several recent works [6, 7, 8, 9] that have been
proposed for Covid-19 detection/ classification from X-ray images. However,
they use a rather small dataset (the largest containing only a few hundreds of
X-ray images), with only a few Covid-19 samples. This makes it difficult to
generalize their results in practice. To address this deficiency and provide
reliable results, in this study the researchers of Qatar University and
Tampere University have compiled the largest Covid-19 dataset, called QaTa-
Cov19. Compared to the earlier benchmark dataset created in this domain, such
as COVID Chestxray Dataset [10] or Covid-19 DATASET [11], QaTa-Cov19 has the
followıng unique benchmarking properties. First, it is the largest dataset,
not only in terms of the number of images (more than 6200 images) but its
versatility i.e., QaTa-Cov-19 contains additional major pneumonia categories,
such as Viral and Bacterial, along with the control (normal) class. Moreover,
this is the most diverse dataset encapsulating X-ray images from several
countries (e.g. Italy, Spain, China, etc.) produced by different X-ray
machines. Finally, the images are in different quality, resolution and SNR
levels as shown in Fig. 1.

Figure 1: Sample Covid-19 X-ray images from QaTa-Cov19.

QaTa-Cov19 contains many X-ray images from the Covid-19 patients who are in
the early stages; therefore, their X-ray images show mild or no-sign of
Covid-19 infestation by the naked eye. Some sample images are shown in Fig.
2-(b). Another fact which makes the diagnosis far more challenging is that
inter-class similarity can be very high for many X-ray images as some samples
shown in Fig. 2-(a). Against such high inter-class similarities and intra-
class variations, in this study we aim for a high robustness level. Our
primary objective is to achieve the highest sensitivity possible in the
diagnosis of Covid-19 induced pneumonia with an acceptable false-alarm rate
(e.g. specificity $>95\%$). In particular, the misdiagnosis of a Covid-19
X-ray image as a normal case should be minimized whilst a small number of
false negatives is tolerable.

(a)

(b)

Figure 2: Sample QaTa-Cov19 X-ray images: (a) X-ray images from different
classes. (b) X-ray images from the Covid-19 patients who are in the different
stages.

In numerous classification tasks, Deep Learning techniques have been shown to
achieve state-of-the-art performance in term of both recognition accuracy and
their parallelizable computing structures which play an important role
especially in real-time applications. Despite their advantages, in order to
achieve a desired performance level in a deep model, a proper training over a
massive training dataset is usually needed. Nevertheless, this is
unfortunately not an option yet for this problem since the available data is
still rather limited.

An alternative supervised approach, which requires a limited number of
training samples to achieve satisfactory classification accuracy is
representation-based classification [12, 13, 14]. In representation-based
classification systems, a dictionary, whose columns consist of the training
samples that are stacked in such a way that a subset of them corresponding to
a class, is pre-defined. A test sample is expected to be a linear combination
of all points from the same class as the test sample. Therefore, given a
predefined dictionary matrix, $\mathbf{D}$ and a test sample $\mathbf{y}$, we
expect the solution $\mathbf{\hat{x}}$ from $\mathbf{y}=\mathbf{Dx}$, carry
enough information about the class of $\mathbf{y}$. The two well-known
representation based classification methodologies are sparse representation-
based classification (SRC) [13] and collaborative representation based
classification (CRC) [12]. Out of these two, SRC provides slightly improved
accuracy by solving a sparse representation problem, i.e., producing a sparse
solution $\mathbf{\hat{x}}$ from $\mathbf{y}=\mathbf{Dx}$. Then, the location
of the non-zero elements of $\mathbf{\hat{x}}$, which is also known as support
set, provides us with the class of the query $\mathbf{y}$. Despite improved
recognition accuracy, SRC solutions are iterative solutions and can be
computational demanding compared to CRC. In a recent work [15], a compact
neural network design that can be considered as a bridge between learning-
based and representation-based methodologies was proposed. The so-called
Convolutional Support Estimation Network (CSEN) uses a pre-defined dictionary
and learns a direct mapping using moderate/low size training set, which maps
query samples, $\mathbf{y}$, directly to the support set of representation
coefficients, $\mathbf{x}$ (as it should be purely sparse in the ideal case).

In this study, to address the aforementioned limitations in Covid-19 diagnosis
from X-ray images we propose a CSEN-based approach. Since the largest set of
Covid-19 X-ray images ever compiled is used in this study, the proposed
approach can be evaluated rigorously against a high-level of diversity to
obtain a reliable analysis. The general pipeline of the proposed CSEN based
recognition scheme is illustrated in Fig. 3. In order to obtain highly
discriminative features, we use the recently proposed CheXNet [16], which is
the fine-tuned version of $121$ layer Dense Convolutional Network
(DenseNet-121) [17] by using over $100000$ frontal view X-ray images form $14$
classes. Having the pre-trained CheXNet for feature extraction, we develop two
different strategies to obtain the classes of query X-ray images: (i) using
collaborative representation-based classification with a proper pre-
processing; (ii) a slightly modified version of our recently proposed
convolution support estimator (CSEN) models. The proposed CSEN scheme
outperforms the competing methods and achieves over $98\%$ of sensitivity and
over $95\%$ for specificity in this challenging dataset.

Figure 3:  The proposed approach for Covid recognition from X-ray images. The
proposed convolution support estimator network (CSEN) which can be trained
from a moderate size training set. The pipeline employs the pre-trained deep
neural network for feature extraction. $\mathbf{A}$ is the dimensional
reduction (PCA) matrix, the coarse estimation of representation coefficient
(sparse in ideal case), $\hat{\mathbf{x}}$ is obtained via the denoiser
matrix,
$\mathbf{B}=\left(\mathbf{D}^{T}\mathbf{D}+\lambda\mathbf{I}\right)^{-1}\mathbf{D}^{T}$,
where $\mathbf{D}=\mathbf{A}\mathbf{\Phi}$ and $\mathbf{\Phi}$ is the pre-
defined dictionary matrix of training samples (before dimensional reduction).

The rest of the paper is organized as follows. In Section II, notations and
mathematical preliminaries are given with emphasis on sparse representation
and sparse support estimation. Then in Section III, a literature review on
deep learning models over X-ray images and representation based classification
is presented. The proposed CSEN-based Covid-19 recognition system is
introduced in Section IV along with two recent alternative approaches that are
used as the competing methods. The data collection is also explained in this
section. Experimental setup and the main results are provided in Section V.
Finally, Section VI concludes the paper and suggests topics for future
research.

##  II Preliminaries and Mathematical Notations

###  II-A Notations

In this study, the $\ell_{p}$-norm of a vector $\mathbf{x}\in\mathbb{R}^{n}$
is defined as
$\left\|\mathbf{x}\right\|_{\ell_{p}^{n}}=\left(\sum_{i=1}^{n}\left|x_{i}\right|^{p}\right)^{1/p}$
for $p\geq 1$. On the other hand, the $\ell_{0}$-norm of the vector
$\mathbf{x}\in\mathbb{R}^{n}$ is defined as
$\left\|\mathbf{x}\right\|_{\ell_{0}^{n}}=\lim_{p\to
0}\sum_{i=1}^{n}\left|x_{i}\right|^{p}=\\#\\{j:x_{j}\neq 0\\}$ and the
$\ell_{\infty}$-norm is defined as
$\left\|\mathbf{x}\right\|_{\ell_{\infty}^{n}}=\max_{i=1,...,n}\left(\left|x_{i}\right|\right)$.
A signal $\mathbf{s}$ is called strictly $k$-sparse if
$\left\|\mathbf{x}\right\|_{0}\leq k$. Sparse support set or simply support
set, $\Lambda\subset\\{1,2,3,...,n\\}$ of sparse signal $\mathbf{x}$ can be
defined as the set of non-zero coefficients’ location, i.e.,
$\Lambda:=\left\\{i:x_{i}=0\right\\}$.

###  II-B Sparse Signal Representation

Sparse representation (SR) of a signal $\mathbf{s}\in\mathbb{R}^{d}$ in a pre-
defined set of waveforms, $\mathbf{\Phi}\in\mathbb{R}^{d\times n}$, can be
defined as representing $\mathbf{s}$ as a linear combination of only a small
subset of atoms of in the dictionary $\mathbf{\Phi}$, i.e,
$\mathbf{s}=\mathbf{\Phi}\mathbf{x}$. Defining these sets, which dates back to
Fourier’s pioneering work [18], has been excessively studied in the
literature. In the early approaches, these sets of waveforms have been
selected as a collection of linearly independent and generally orthogonal
waveforms (which are called a complete dictionary or basis i.e, $d=n$) such as
Fourier Transform, DCT and Wavelet Transform, until the pioneering work of
Mallat [19] on overcomplete dictionaries ($n>>d$). In the last decade,
interest in SR research increased tremendously and their wide range of
applications includes denoising [20], classification [21], anomaly detection
[22, 23], Deep Learning [24] and Compressive Sensing (CS) [25, 26].

With a possible dimensional reduction that can be satisfied via a compression
matrix $\mathbf{A}\in\mathbb{R}^{m\times d}$ ($m<<d$), sample can be obtained
from $\mathbf{s}$,

$\mathbf{y}=\mathbf{A}\mathbf{s}=\mathbf{A}\mathbf{\Phi}\mathbf{x}=\mathbf{D}\mathbf{x},$
(1)

where $\mathbf{D}\in\mathbb{R}^{m\times n}$ can be called the equivalent
dictionary. Because Eq. (1) describes an under-determined system of linear
equations, finding the representation coefficient vector $\mathbf{x}$ requires
at least one more constraint to have a unique solution. Using the prior
information about sparsity, the following representation

$\min_{\mathbf{x}}~{}\left\|\mathbf{x}\right\|_{0}~{}\text{subject
to}~{}\mathbf{D}\mathbf{x}=\mathbf{y}$ (2)

which is also a sparse representation of $\mathbf{x}$ has a unique solution
provided that $\mathbf{D}$ satisfies some required properties [27]. However,
the optimization problem in Eq. (2) is a NP-hard. Fortunately, the following
relaxation

$\min_{\mathbf{x}}~{}\left\|\mathbf{x}\right\|_{1}~{}\text{subject
to}~{}\mathbf{D}\mathbf{x}=\mathbf{y}$ (3)

produces exactly the same solution as that of Eq. (2) provided that
$\mathbf{D}$ obeys some criteria [28] and $m>k(\log(n/k))$. In addition, real
world applications generally exhibit not exact sparsity but approximate
sparsity. Furthermore, the query sample $\mathbf{y}$ can be corrupted with an
additive noise pattern. In this case, the equality constraint in Eq. (3) can
be further relaxed such as in the Basis Pursuit Denoising (BPDN) [29]:
$\min_{\mathbf{x}}\left\|\mathbf{x}\right\|~{}\text{s.t.}~{}\left\|\mathbf{y}-\mathbf{Dx}\right\|\leq\epsilon$,
where $\epsilon$ is a small constant that depends on the noise level.

We may refer to the Sparse Support Estimation (SE) problem as finding the
indices a set, $\Lambda$, of non-zero elements of $\mathbf{x}$ [30, 31].
Indeed, in many applications, SE can be more important than finding the
magnitude and sign of $\mathbf{x}$ as well as $\Lambda$, which refers to the
sparse Signal Recovery (SSR) via a recovery technique, such as Eq. (3). For
example, in a sparse representation based classification system, a query
sample $\mathbf{y}$ can be represented with sparse coefficient vector,
$\mathbf{x}$, in the dictionary, $\mathbf{D}$ in such a way that when we
recover this representation coefficient from
$\mathbf{y}=\mathbf{D}\mathbf{x}$, the solution vector $\mathbf{\hat{x}}$ is
expected to have a significant number of non-zero coefficients coming from the
particular locations corresponding to the class of $\mathbf{y}$.

Readers are referred to [15] for more detailed literature review on SE and its
applications. In the sequel, we briefly summarize the building blocks of the
proposed approach.

##  III Background and Prior Art

###  III-A CheXNet

In the proposed approach, we first use the pre-trained deep network, CheXNet,
to extract discriminative features from raw X-ray images. CheXNet was
developed for pneumonia detection from the chest X-ray images [16]. In [16],
it was claimed that their CheXNet can perform even better than expert
radiologist in the pneumonia detection problem. This deep neural network
design is based on previously proposed DenseNet [17] that consists of 121
layers. It is first pre-trained over ImageNet dataset [32] and performed
transfer learning over $112120$ frontal-view chest X-ray images in the ChestX-
ray14 dataset [33].

###  III-B Representation Based Classification

Given a test sample $\mathbf{y}$, which represents either the extracted
features, $\mathbf{s}$, or their dimensionally reduced version, i.e.,
$\mathbf{y}=\mathbf{A}\mathbf{s}$. In developing the dictionary, training
samples are stacked in $\mathbf{D}$ with particular locations in such a way
that the optimal support for a given query $\mathbf{y}$ should be the set of
all points coming from the same class as $\mathbf{y}$. Therefore, a solution
vector, $\mathbf{\hat{x}}$ of $\mathbf{y}=\mathbf{D}\mathbf{x}$ is supposed to
have enough information, i.e., the sparse support should be the set of
location indices of the training sample from the same class as $\mathbf{y}$.
This strategy is generally known as representation-based classification.
However, a typical solution $\mathbf{\hat{x}}$ of
$\mathbf{y}=\mathbf{D}\mathbf{x}$ is not necessarily a sparse one especially
when its size grows with more training samples, which results in a highly
under-determined system of linear equations. Fortunately, if one estimates the
representation coefficient vector with a sparse recovery design such as
$\ell_{1}$-minimization as in Eq. (3), we can expect that the important non-
zero entries of the solution, $\mathbf{\hat{x}}$, are grouped in the
particular locations that correspond to the locations of the training samples
from the same class as $\mathbf{y}$. This can be a typical example of
scenarios where support estimation can be more valuable than the magnitudes
and sign recovery as explained in Section II-B.

For instance, [14] proposed a systematic way of determining the identity of
face images using $\ell_{1}$-minimization. The authors develop a three-step
classification technique that includes: (i) normalization of all the atoms in
$\mathbf{D}$ and $\mathbf{y}$ to have unit $\ell_{2}$-norm; (ii) estimating
the representation coefficient vector via sparse recovery, i.e.,
$\hat{\mathbf{x}}=\arg\min_{\mathbf{x}}\left\|\mathbf{x}\right\|_{1}\text{s.t}\left\|\mathbf{y}-\mathbf{D}\mathbf{x}\right\|_{2}$;
and (iii) finding the residuals corresponding to each class via
$\mathbf{e_{i}}=\left\|\mathbf{y}-\mathbf{D_{i}}\mathbf{\hat{x}_{i}}\right\|_{2}$,
where $\mathbf{\hat{x}_{i}}$ is the group of the estimated coefficients,
$\mathbf{\hat{x}}$, that correspond to class $i$.

This technique, which is known as Sparse Representation based Classification
(SRC), and its variants have been applied to a wide range of applications in
literature [34, 35], e.g., human action recognition [36], and hyperspecral
image classification [37], to name a few. Despite the good recognition
accuracy performance of SRC systems, their main drawbacks is the fact that
their sparse recovery algorithms (e.g., $\ell_{1}$-minimization) is iterative
methods and computationally costly, rendering them infeasible in real time
applications. Later, the authors of [12] introduced Collaborative
Representation based Classification (CRC), which is similar to SRC except for
the use of traditional $\ell_{2}$-minimization in the second step;
$\mathbf{\hat{x}}=\arg\min_{\mathbf{x}}\left\\{\left\|\mathbf{y}-\mathbf{D}\mathbf{x}\right\|_{2}^{2}
+\lambda\left\|\mathbf{x}\right\|_{2}^{2} \right\\}$. Thus, CRC does not
require an iterative solution to obtain representation coefficient thanks to
that $\ell_{2}$-minimization has a closed form solution,
$\hat{\mathbf{x}}=\left(\mathbf{D}^{T} \mathbf{D}+\lambda\mathbf{I}_{n\times
n} \right)^{-1}\mathbf{D}^{T} \mathbf{y}$. Although, the sparsity in
$\mathbf{\hat{x}}$ cannot be guaranteed, it has often been reported to achieve
a comparable classification performance, especially in small-size training
datasets.

##  IV Proposed Approach

Figure 4: Baseline Approach I: collaborative representation based
classification is fed by deep learning based extracted features that are pre-
processed.  Figure 5: Baseline Approach II: A 5-layer MLP layer is used over
the features of CheXNet.

###  IV-A The Benchmark Dataset: QaTa-Cov19

Covid-19 chest X-ray images were gathered from different publicly available
but scattered image sources. However, the major sources of Covid-19 images are
Italian Society of Medical and Interventional Radiology (SIRM) COVID-19
Database [11], Radiopaedia [38], Chest Imaging (Spain) at thread reader [39]
and online articles and news-portals. The authors have carried out the task of
collecting and indexing the X-ray images for Covid-19 positive cases reported
in the published and preprint articles from China, South Korea, USA, Taiwan,
Spain, and Italy, as well as online news-portals (up to 20th April 2020).
Therefore, these X-ray images represent different age groups, gender,
ethnicity and country. Negative Covid19 cases were normal, viral and bacterial
pneumonia chest X-ray images and collected from the Kaggle chest X-ray
database. Kaggle chest X-ray database contains 5863 chest X-ray images of
normal, viral and bacterial pneumonia with varying resolutions [40]. Out of
these $5863$ chest X-ray images, $1583$ images are normal images and the
remaining are bacterial and viral pneumonia images. Sample X-ray images from
QaTa-Cov19 dataset are shown in Fig. 6.

Figure 6: Samples from the benchmark QU-Chest dataset.

###  IV-B Feature Extraction

With their outstanding performance in image classification along with other
inference tasks, deep neural networks became a dominant paradigm. However,
these techniques usually necessitates a large number of training samples
(e.g., several hundred-thousand to millions depending on the network size) to
achieve an adequate generalization capability. That is to say, the
aforementioned problem of the data scarcity with the Covid-19 case prevents us
from training a deep learning technique from scratch. Albeit, we can still
leverage their power by finding properly pre-trained models for similar
problems. To this end, we use a-state-of-the art pneumonia detection network,
CheXNet, whose details are summarized in Section III-A. With the pre-trained
model, we extract $1024$-long vectors, right after the last average pooling
layer. After data normalization (zero mean and unit variance), we obtain a
feature vector $\mathbf{s}\in\mathbb{R}^{d=1024}$.

A dimensionality reduction PCA is applied to $\mathbf{s}$ in order to get the
query sample, $\mathbf{y}=\mathbf{A}\mathbf{s}\in\mathbb{R}^{m}$, where
$\mathbf{A}\in\mathbb{R}^{m\times d}$ is PCA matrix ($m<d$).

###  IV-C The proposed CSEN-based Classification

Considering the limited number of training data in our Covid-19 dataset, a
representation-based classification can be applied hereafter to obtain the
class of $\mathbf{y}$ using the dictionary $\mathbf{\Phi}$ (in the form of
$\mathbf{D}=\mathbf{A\Phi}$), whose columns are stacked training samples with
class-specific locations.

As discussed earlier, sparse representation-based classification is a support
estimation problem which is expected to be an easier task than a sparse signal
recovery problem. On the other hand, even if the exact signal recovery is not
possible in noisy cases or in cases where $\mathbf{\hat{x}}$ is not exactly
but approximately sparse (which is the case in almost all the time in
dictionary-based classification problems), it is still possible to recover the
support set exactly [9, 30, 41, 42] or partially [42, 43, 44]. However, many
works in the literature dealing with SE problems tend to first apply a sparse
recovery technique on $\mathbf{y}$ to first get $\mathbf{\hat{x}}$, then use
simple thresholding over $\mathbf{\hat{x}}$ to obtain a sparse support
estimation, $\hat{\Lambda}$. Nevertheless, SSR techniques such as
$\ell_{1}$-minimization are rather slow and their performance varies from one
SRR tool to another [15]. In our previous work [15], we proposed an
alternative solution for this handcrafted sparse recovery approach which aims
to learn a direct map from test sample $\mathbf{y}$ to the support set
$\hat{\Lambda}$. Along with the speed and stability compared to conventional
SSR based techniques, and recent deep learning based solutions to SRR problem,
CSEN has a crucial advantage of having a compact design that can achieve a
good performance level even over scarce training data.

Mathematically speaking, an ideal CSEN is supposed to yield a binary mask
$\mathbf{v}\in\left\\{0,1\right\\}^{n}$:

${v_{i}=} 1 ~{}\text{ if $i\in\Lambda$ }$ (4)

which indicates the true support i.e.,
${\Lambda}=\left\\{i\in\left\\{1,2,..,n\right\\}:{v}_{i}=1\right\\}$. In order
to approximate this ideal case, a CSEN network,
$\mathcal{P}\left(\mathbf{y},\mathbf{D}\right)$ produces a probability vector
$\mathbf{p}$ which returns a measure about the probability of each index being
in ${\Lambda}$ such that $p_{i}\in\left[0,1\right]$. Having the estimated
probability map, estimating the support can easily be done via
$\hat{\Lambda}=\left\\{ i\in \left\\{ 1,2,..,n\right\\}:p_{i}>\tau \right\\}$,
by thresholding $\mathbf{p}$ with $\tau$ where $\tau$ is a fixed threshold.

A CSEN is composed of fully convolutional layers, and as input it takes a
proxy, $\mathbf{\tilde{x}}$, of sparse coefficient vector, which is a coarse
estimation of $\mathbf{x}$ i.e.,
$\left(\mathbf{D}^{T}\mathbf{D}+\lambda\mathbf{I}\right)^{-1}\mathbf{D}^{T}\mathbf{y}$
or simply $\mathbf{\tilde{x}=D^{T}y}$. Using such a proxy of $\mathbf{x}$,
instead of making inference directly on $\mathbf{y}$ has also studied in a few
more recent studies. For instance, In [45, 46], the authors proposed
reconstruction-free image classification from compressively sensed images.

The input vector $\mathbf{\tilde{x}}$ is reshaped to a 2-D plane in order to
use is with 2-D convolutional layers. This transformation is performed via re-
ordering the indices of the atoms in such a way that the non-zero elements of
the representation vector $\mathbf{x}$ for a specific class come together in
the 2-D plane. A representative illustration of the proposed dictionary design
with compared to the traditional one is shown in Fig. 7.

Hereafter the proxy $\mathbf{\tilde{x}}$ is convolved with the weight kernels,
connecting the input with the next layer with $N$ filters to yield the inputs
of the next layer, with the biases $\mathbf{b_{1}}$ as follows:

$\mathbf{f_{1}}=\\{S(ReLu(b_{1}^{i}+\mathbf{w}_{1}^{i}*\tilde{\mathbf{x}}))\\}_{i=1}^{N},$
(5)

where $\mathbf{b_{1}}$ is the weight bias, $S(.)$ is the down- or up-sampling
operation and $ReLu(x)=max(0,x)$. In more general form, the $k^{th}$ feature
map of layer $l$ is defined as,

$\mathbf{f_{l}^{k}}=\textsc{S}(\textsc{ReLu}(b_{l}^{k}+\sum_{i=1}^{N_{l-1}}\textsc{conv2D}(\mathbf{w}_{l}^{ik},\mathbf{f}_{l-1}^{i},^{\prime}\textsc{ZeroPad}^{\prime}))).$
(6)

Therefore, the trainable parameters of CSEN will be:
$\mathbf{\Theta_{CSEN}}=\big{\\{}\\{\mathbf{w}_{1}^{i},b_{1}^{i}\\}_{i=1}^{N_{1}},\\{\mathbf{w}_{2}^{i},b_{2}^{i}\\}_{i=1}^{N_{2}},...\\{\mathbf{w}_{L}^{i},b_{L}^{i}\\}_{i=1}^{N_{L}}\big{\\}}$
for a L layer CSEN design.

In developing the dictionary that is to be used in the sparse representation
based classification, the training samples are stacked-in by grouping of them
according to their classes. Thus, instead of using traditional
$\ell_{1}$-minimization formulation as in Eq. (3), the following group
$\ell_{1}$-minimization formulation may result in increased classification
accuracy,

$\min_{\mathbf{x}}\left\\{\left\|\mathbf{D}\mathbf{x}-\mathbf{y}\right\|_{2}^{2}+\lambda\sum_{i=1}^{c}\left\|\mathbf{x_{Gi}}\right\|_{2}\right\\}$
(7)

where $\mathbf{x_{Gi}}$ is the group of coefficients from the $i^{t}h$ class.
In this manner, one possible cost function for a SE network would be,

$E(\mathbf{x})=\sum_{p}(\mathcal{P}_{\Theta}\left(\mathbf{\tilde{x}}\right)_{p}-v_{p})^{2}+\lambda\sum_{i=1}^{c}\left\|\mathcal{P}_{\Theta}\left(\mathbf{\tilde{x}}\right)_{Gi}\right\|_{2}.$
(8)

where $\mathcal{P}_{\Theta}\left(\mathbf{\tilde{x}}\right)_{p}$ is network
output at location $p$ and $v_{p}$ is the ground truth binary mask of the
sparse code $\mathbf{x}$. Due to its high computational complexity, we
approximate the cost function in (8) with a simpler average pooling layer
after convolutional layer, which can produce directly the estimated class in
our CSEN design. An illustration of proposed CSEN-based Covid-19 recognition
is shown in Fig. 3.

Figure 7: The illustration of proposed dictionary design vs. conventional
design in representation based classifiers.

###  IV-D Competing Methods

This section summarizes the competing methods that are selected among numerous
alternatives due to their superior performance levels obtained in similar
problems. For a fair comparative evaluations, all classification methods have
the same input feature vectors fed to the proposed CSENs.

####  IV-D1 Collaborative representation-based classification

As a possible competing technique to the proposed CSEN based technique which
is a hybrid method, CRC [12] is a direct and representation-based
classification method. It is a non-iterative support estimation technique,
that satisfies faster and comparable classification performance with SRC while
it is more stable compared to existing iterative sparse recovery tools as it
is shown in [15]. In the first step of CRC, the trade-off parameter of
regularized least square solution is set as $\lambda=2*10^{-12}$.

####  IV-D2 Multi-layer Perceptron (MLP) classification

As one of the most-common classifiers, a 4-hidden layer MLP is used for this
problem. For training we used Back-Propagation (BP) with Adam optimization
technique [47]. The network and training hyper-parameters are as follows:
learning rate, $\alpha=10^{-4}$, and moment updates $\beta_{1}=0.9$,
$\beta_{2}=0.999$, and $50$ as the number of epochs. Fig. 8 illustrates the
network configuration in detail. This network configuration has achieved the
best performance among others (deeper and shallower) where deep configurations
have suffered from over-fitting while the shallow ones exhibit an inferior
learning performance.

Figure 8: The MLP configuration.

####  IV-D3 Support Vector Machines (SVMs)

For a multi-class problem, the first objective is to select the SVM topology
for ensemble learning: one-vs-one or one-vs-all. In order to find the optimal
topology and the hyper-parameters (e.g. kernel type and its parameters) we
first performed a grid-search with the following variations and setting:
kernel function {linear, radial basis function (RBF)}, box constraint ($C$
parameter) in the range $[1,10^{3}]$ with log scale, and kernel scale
($\gamma$ for the RBF kernel) in the range $[10^{-4},10^{-2}]$ with log scale.

####  IV-D4 k-Nearest-Neighbor (k-NN)

Finally, we use a traditional approach, k-Nearest Neighbor (k-NN) is used with
PCA dimensionality reduction. In a similar fashion, the distance metric and
the k-value are optimized by a prior grid-search. The following distance
metrics are evaluated: City-block, Chebyshev, correlation, cosine, Euclidean,
Hamming, Jaccard, Mahalanobis, Minkowski, standardized Euclidean, and Spearman
metrics. The k-value is varied within the range of $[1,4416]$ with log scale.

TABLE I: Classification Performances of the proposed CSEN and competing methods. The best Covid-19 recognition (sensitivity) rates are highlighted. | Bacterial | Viral | Normal | Covid-19 | Bacterial | Viral | Normal | Covid-19 | Bacterial | Viral | Normal | Covid-19
---|---|---|---|---|---|---|---|---|---|---|---|---
| Accuracy | Sensitivity | Specificity
NN | 0.777 | 0.801 | 0.903 | 0.950 | 0.623 | 0.612 | 0.899 | 0.965 | 0.898 | 0.859 | 0.904 | 0.949
SVM | 0.771 | 0.788 | 0.928 | 0.928 | 0.586 | 0.632 | 0.911 | 0.981 | 0.916 | 0.837 | 0.933 | 0.924
MLP | 0.761 | 0.765 | 0.923 | 0.947 | 0.620 | 0.561 | 0.885 | 0.965 | 0.872 | 0.828 | 0.936 | 0.946
CRC | 0.820 | 0.827 | 0.928 | 0.955 | 0.758 | 0.550 | 0.922 | 0.968 | 0.869 | 0.913 | 0.930 | 0.954
ReconNET | 0.765 | 0.785 | 0.918 | 0.936 | 0.590 | 0.625 | 0.891 | 0.970 | 0.902 | 0.834 | 0.927 | 0.933
CSEN1 | 0.793 | 0.805 | 0.926 | 0.955 | 0.656 | 0.642 | 0.906 | 0.985 | 0.901 | 0.856 | 0.932 | 0.953
CSEN2 | 0.794 | 0.803 | 0.927 | 0.959 | 0.659 | 0.646 | 0.904 | 0.985 | 0.900 | 0.852 | 0.934 | 0.957

##  V Experimental Results

###  V-A Experimental Setup

We have performed our experiments over the QaTa-Cov19 dataset, which consists
of normal and three pneumonia classes: bacterial, viral, and Covid-19. The
proposed approach is evaluated using a stratified 5-fold cross-validation (CV)
scheme with a ratio of 80% for training and 20% for the test (unseen folds)
splits, respectively.

TABLE II: Number of images per class and per-fold before and after data augmentation. Class | # of Samples |  | Training
---
Samples
| Augmented
---
Training Samples
| Test
---
Samples
| Bacterial
---
Pneumonia
2760 | 2208 | 2208 | 552
| Viral
---
Pneumonia
1485 | 1188 | 2208 | 297
Normal | 1579 | 1263 | 2208 | 316
Covid-19 | 462 | 370 | 2208 | 92
Total | 6286 | 5029 | 8832 | 1257

Table II shows the number of X-ray images per class in the QaTa-Cov19 dataset.
Since the dataset is unbalanced, we have applied data augmentation to the
training set in order to balance the size of each class in the train set.
Therefore, the X-ray images in viral and Covid-19 pneumonia, and normal
classes are augmented up to the same number as the bacterial pneumonia class
in the train set. We use Image Data Generator by Keras to perform data
augmentation by applying ZCA whitening with epsilon of $10^{-6}$, randomly
rotating the X-ray images in a range of 10 degrees, randomly shifting images
both horizontally and vertically within the interval of $[-0.1,+0.1]$. In each
CV fold, we use a total of 8832 and 1257 images in the train and test (unseen
in the fold) sets, respectively.

The experimental evaluations of SVM, k-NN and CRC are performed using MATLAB
version 2019a, running on PC with Intel ® i7-8650U CPU and 32 GB system
memory. On the other hand, MLP and CSEN methods are implemented using
Tensorflow library [48] with Python on NVidia ® TITAN-X GPU card. For the CSEN
training, ADAM optimizer [47] is used with the proposed default learning
parameters: learning rate, $\alpha=10^{-3}$, and moment updates
$\beta_{1}=0.9$, $\beta_{2}=0.999$ with only 15 Back-Propagation epochs.
Neither grid-search nor any other parameter or configuration optimization was
performed for CSEN.

###  V-B Experimental Results

The same network configurations are used for CSEN as in [15]. Accordingly, we
use two compact CSEN designs: CSEN1 and CSEN2, respectively. The first CSEN
network consists of only two hidden convolutional layers, the first layer has
48 neurons and the second has 24. ReLu activation function is used in the
hidden layers and the filter size was $3\times 3$. On the other hand CSEN2
uses max-pooling and has one additional hidden layer with 24 neurons to
perform transposed-convolution. CSEN1 and CSEN2 are compared against the 6
competing methods under the same experimental setup.

For the dictionary construction in $\mathbf{\Phi}$ each CSEN design, $625$
images for each class (from the augmented training samples per fold) are
stacked in a such way that the representation coefficient in the 2-D plane,
$\mathbf{X}$ has $50\times 50$ size as shown in Fig. 7. The rest of the images
in the training set are used to train each CSEN i.e., $1583$ samples from each
class. We use PCA dimensional reduction matrix, $\mathbf{A}$ with the
compression ratio, $CR=\frac{m}{d}=0.5$. Therefore, we have $512\times 2500$
equivalent dictionary, $\mathbf{D}$, and $2500\times 512$ denoiser
$\mathbf{B}=\left(\mathbf{D}^{T}\mathbf{D}+\lambda\mathbf{I}\right)^{-1}\mathbf{D}^{T}$
to obtain a coarse estimation of the representation (sparse in ideal case)
coefficients, $\mathbf{\tilde{x}}\in\mathbb{R}^{n=2500}$. Hereafter, the CSEN
networks are trained to have class of $\mathbf{y}$ from input
$\mathbf{\tilde{x}}$ as illustrated in Fig., 3.

Due to the lack of other learning-based SE studies in the literature, we chose
a deeper network compared to CSEN designs to investigate the role of network
depth in this problem. ReconNet [49] was proposed as a non-iterative deep
learning solution to compressive sensing problem i.e.,
$\mathbf{\hat{s}}\leftarrow\mathcal{P}\left(\mathbf{y}\right)$ and it is one
of the state-of-the-art in compressively sensed image recognition task. It
consists of 6 fully convolutional layers and one dense layer in front of the
convolutional ones, which act as the learned denoiser for the mapping from
$\mathbf{y}\in\mathbb{R}^{m}$ to $\mathbf{\tilde{s}}\in\mathbb{R}^{d}$. Then,
the convolutional layers are responsible for producing the reconstructed
signal, $\mathbf{\hat{s}}$ from $\mathbf{\tilde{s}}$. Therefore, by replacing
this dense layer with the denoiser matrix $\mathbf{B}$, this network can be
used as a competing method.

Both CSEN and the modified ReconNet use $\mathbf{\tilde{x}}$ as a input, which
is produced using an equivalent dictionary $\mathbf{D}$ and its pseudo-inverse
matrix $\mathbf{B}$.

Figure 9: False Negatives of proposed Covid-19 recognition scheme.

In designing the dictionary of CRC system, all training samples are stacked in
the dictionary, $\mathbf{\Phi}$, i.e., 2208 samples from each class. The same
PCA matrix used in CSEN based recognition, $\mathbf{A}$ is applied to
features, $\mathbf{s}\in\mathbb{R}^{d=1024}$. Therefore, a dictionary
$\mathbf{D}$ of size $512\times 8832$ and the corresponding denoiser matrix
$\mathbf{B}$ of size $8832\times 512$ are used in the CRC framework.

TABLE III: The number of network parameters of each method.  | MLP | CSEN1 | CSEN2 | ReconNet
---|---|---|---|---
| # of trainable
---
parameters
672,836 | 11,089 | 16,297 | 22,914
TABLE IV: Computation times (sec) of each method over 1257 test images. |  | CRC
---
(light)
CRC | CSEN1 | CSEN2 | ReconNet | MLP
| Computation
---
Time (in sec.)
13.4176 | 40.7878 | 0.2196 | 0.2272 | 0.2993 | 0.2935

The classification performance of the proposed CSEN-based approach and the
competing methods is presented in Table I. As can be easily observed from the
Table I, the proposed approaches surpass all competing methods in Covid-19
recognition performance by achieving $98.5\%$ sensitivity, and over $95\%$
specificity. As shown in Table III, compared to MLP and ReconNet, the proposed
CSEN designs are very compact, and computationally efficient. This is evident
in Table IV where the computational complexity (measured as total computation,
time over the 1257 test images) is reported.

When compared against CRC in particular, CSEN-based classification has two
advantages; computational efficiency and, a superior Covid-19 recognition
performance. The computational efficiency comes from the fact that a larger
size dictionary matrix (of size of $512\times 8832$) is used in CRC and hence,
this requires more computations in terms of matrix-vector multiplications.
Furthermore, saving the trainable parameters ($\sim 16k$) and a light
dictionary matrix coefficients ($\sim 1280k$) in the test device is more
memory efficient compared to saving coefficients ($\sim 4521k$) of larger size
dictionary used in CRC.

TABLE V: Performance of CRC algorithm when the dictionary (size of 625 per class) that is used in CSEN is used.  | CRC (Light)
---|---
| Accuracy | Sensitivity | Specificity
Bacterial | 0.8129 | 0.7464 | 0.8650
Viral | 0.8163 | 0.5461 | 0.8998
Normal | 0.9267 | 0.9170 | 0.9299
Covid-19 | 0.9564 | 0.9394 | 0.9578

For further analysis, we also tested the CRC framework by using the light
dictionary (of size $512\times 2500$) used in CSEN based recognition. We
called it CRC (light), and as it can be seen in Table V, the performance of
CRC further reduced, and there was no significant improvement concerning the
computational cost. When it comes to creating deeper convolutional layers
instead of using CSEN designs, such as the modified ReconNet, the results
presented in Table I shows us that compact CSEN structures are indeed
preferable to achieve superior classification performances compared to deeper
networks.

TABLE VI: The overall (cumulative) confusion matrix of the proposed recognition scheme. CSEN2 | Predicted
---|---
|  | Bacterial | Viral | Normal | Covid-19
| Bacterial | 1818 | 636 | 180 | 126
| Viral | 338 | 959 | 127 | 61
| Normal | 15 | 71 | 1428 | 65
Real | Covid-19 | 0 | 3 | 4 | 455

Finally, Table VI presents the overall (cumulative) confusion matrix of the
proposed CSEN-based Covid-19 recognition approach over the new QaTa-Cov19
Dataset. The most critical mis-classifications are the false-positives, that
is, the mis-classified Covid-19 X-ray images. The confusion matrix shows that
the proposed approach has mis-classified 7 Covid-19 images (out of 462). The 3
out of 7 misclassifications are still in “Viral Pneumonia” category, which can
be an expected confusion due to the viral nature of Covid-19. However, the
other four cases are mis-classified as “Normal” which is indeed a severe
clinical misdiagnosis. A close look to these false-negatives in Fig. 9 reveals
the fact that they are indeed very similar to normal images where typical
Covid-19 patterns are hardly visible even by an expert’s naked eye. It is
possible that these images come from the patients who were in the very early
stages of Covid-19.

##  VI Conclusions

The commonly used methods in Covid-19 diagnosis, namely Reverse Transcription-
Polymerase Chain Reaction and Computed Tomogrophy have certain limitations and
drawbacks such as long processing times and unacceptably high mis-diagnosis
rates. These drawbacks are also shared by most of the recent works in the
literature based on deep learning due to the data scarcity from the Covid-19
cases. Although Deep Learning based recognition techniques are dominant in
Computer Vision where they achieved state-of-the-art performance, their
performance degrades fast due to data scarcity, which is the reality in this
problem at hand. This study aims to address such limitations by proposing a
robust and highly accurate Covid-19 recognition approach directly from raw
X-ray images without any pre- or post-processing. The proposed approach is
based on the CSEN that can be seen as a bridge between Deep Learning models
and representation-based methods. CSEN uses both a dictionary and a set of
training samples to train direct map from the query samples to the sparse
support set of representation coefficients. With this unique ability and
having the advantage of a compact network, the proposed CSEN-based Covid-19
recognition systems surpass the competing methods and achieve over $98\%$
sensitivity and over $95\%$ specificity. Furthermore, they yield the most
computationally efficient scheme in terms of speed and memory. Finally, the
largest dataset of X-ray images, QaTa-Cov19 will be released along with this
study as a benchmark dataset in this domain. This will, henceforth, accelerate
the research efforts globally and support the fight against Covid-19
worldwide.

## References

  * [1] L. Pellis, F. Scarabel, H. B. Stage, C. E. Overton, L. H. Chappell, K. A. Lythgoe, E. Fearon, E. Bennett, J. Curran-Sebastian, R. Das _et al._ , “Challenges in control of covid-19: short doubling time and long delay to effect of interventions,” _arXiv preprint arXiv:2004.00117_ , 2020.
  * [2] F. Zhou, T. Yu, R. Du, G. Fan, Y. Liu, Z. Liu, J. Xiang, Y. Wang, B. Song, X. Gu _et al._ , “Clinical course and risk factors for mortality of adult inpatients with covid-19 in wuhan, china: a retrospective cohort study,” _The Lancet_ , 2020.
  * [3] Y. Fang, H. Zhang, J. Xie, M. Lin, L. Ying, P. Pang, and W. Ji, “Sensitivity of chest ct for covid-19: comparison to rt-pcr,” _Radiology_ , p. 200432, 2020.
  * [4] K. A. Erickson, K. Mackenzie, and A. Marshall, “Advanced but expensive technology. balancing affordability with access in rural areas.” _Canadian family physician Medecin de famille canadien_ , vol. 39, pp. 28–30, 1993.
  * [5] W. H. Organization _et al._ , “Laboratory testing for coronavirus disease 2019 (covid-19) in suspected human cases: interim guidance, 2 march 2020,” World Health Organization, Tech. Rep., 2020.
  * [6] M. E. Chowdhury, T. Rahman, A. Khandakar, R. Mazhar, M. A. Kadir, Z. B. Mahbub, K. R. Islam, M. S. Khan, A. Iqbal, N. Al-Emadi _et al._ , “Can ai help in screening viral and covid-19 pneumonia?” _arXiv preprint arXiv:2003.13145_ , 2020.
  * [7] I. D. Apostolopoulos and T. A. Mpesiana, “Covid-19: automatic detection from x-ray images utilizing transfer learning with convolutional neural networks,” _Physical and Engineering Sciences in Medicine_ , p. 1, 2020.
  * [8] L. O. Hall, R. Paul, D. B. Goldgof, and G. M. Goldgof, “Finding covid-19 from chest x-rays using deep learning on a small dataset,” _arXiv preprint arXiv:2004.02060_ , 2020.
  * [9] M. Wainwright, “Information-theoretic bounds on sparsity recovery in the high-dimensional and noisy setting,” in _2007 IEEE International Symposium on Information Theory_. IEEE, 2007, pp. 961–965.
  * [10] J. P. Cohen, P. Morrison, and L. Dao, “Covid-19 image data collection,” _arXiv preprint arXiv:2003.11597_ , 2020.
  * [11] “Covid-19 database,” 2020. [Online]. Available: https://www.sirm.org/category/senza-categoria/covid-19/
  * [12] L. Zhang, M. Yang, and X. Feng, “Sparse representation or collaborative representation: Which helps face recognition?” in _2011 International conference on computer vision_. IEEE, 2011, pp. 471–478.
  * [13] J. Wright, A. Y. Yang, A. Ganesh, S. S. Sastry, and Y. Ma, “Robust face recognition via sparse representation,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 31, no. 2, pp. 210–227, 2008.
  * [14] J. Wright, Y. Ma, J. Mairal, G. Sapiro, T. S. Huang, and S. Yan, “Sparse representation for computer vision and pattern recognition,” _Proceedings of the IEEE_ , vol. 98, no. 6, pp. 1031–1044, 2010.
  * [15] M. Yamac, M. Ahishali, S. Kiranyaz, and M. Gabbouj, “Convolutional sparse support estimator network (csen) from energy efficient support estimation to learning-aided compressive sensing,” _arXiv preprint arXiv:2003.00768_ , 2020\.
  * [16] P. Rajpurkar, J. Irvin, K. Zhu, B. Yang, H. Mehta, T. Duan, D. Ding, A. Bagul, C. Langlotz, K. Shpanskaya _et al._ , “Chexnet: Radiologist-level pneumonia detection on chest x-rays with deep learning,” _arXiv preprint arXiv:1711.05225_ , 2017.
  * [17] G. Huang, Z. Liu, L. Van Der Maaten, and K. Q. Weinberger, “Densely connected convolutional networks,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2017, pp. 4700–4708.
  * [18] O. D. Escoda, L. Granai, and P. Vandergheynst, “Mémoire sur la propagation de la chaleur dans les corps solides,” _Nouveau Bulletin des Sciences de la Société Philomathique de Paris 6_ , pp. 112–116, 1808.
  * [19] S. G. Mallat and Z. Zhang, “Matching pursuits with time-frequency dictionaries,” _IEEE Transactions on signal processing_ , vol. 41, no. 12, pp. 3397–3415, 1993.
  * [20] J.-L. Starck, E. J. Candès, and D. L. Donoho, “The curvelet transform for image denoising,” _IEEE Transactions on image processing_ , vol. 11, no. 6, pp. 670–684, 2002.
  * [21] J. Yang, K. Yu, Y. Gong, T. S. Huang _et al._ , “Linear spatial pyramid matching using sparse coding for image classification.” in _CVPR_ , vol. 1, no. 2, 2009, p. 6.
  * [22] A. Adler, M. Elad, Y. Hel-Or, and E. Rivlin, “Sparse coding with anomaly detection,” _Journal of Signal Processing Systems_ , vol. 79, no. 2, pp. 179–188, 2015.
  * [23] D. Carrera, G. Boracchi, A. Foi, and B. Wohlberg, “Detecting anomalous structures by convolutional sparse models,” in _2015 International Joint Conference on Neural Networks (IJCNN)_. IEEE, 2015, pp. 1–8.
  * [24] W. Wen, C. Wu, Y. Wang, Y. Chen, and H. Li, “Learning structured sparsity in deep neural networks,” in _Advances in neural information processing systems_ , 2016, pp. 2074–2082.
  * [25] D. L. Donoho _et al._ , “Compressed sensing,” _IEEE Transactions on information theory_ , vol. 52, no. 4, pp. 1289–1306, 2006.
  * [26] E. J. Candès _et al._ , “Compressive sampling,” in _Proceedings of the International Congress of Mathematicians_ , vol. 3, 2006, pp. 1433–1452.
  * [27] D. L. Donoho and M. Elad, “Optimally sparse representation in general (nonorthogonal) dictionaries via $ell_{1}$ minimization,” _Proceedings of the National Academy of Sciences_ , vol. 100, no. 5, pp. 2197–2202, 2003.
  * [28] E. J. Candes, “The restricted isometry property and its implications for compressed sensing,” _Comptes rendus mathematique_ , vol. 346, no. 9-10, pp. 589–592, 2008.
  * [29] S. S. Chen, D. L. Donoho, and M. A. Saunders, “Atomic decomposition by basis pursuit,” _SIAM review_ , vol. 43, no. 1, pp. 129–159, 2001.
  * [30] W. Wang, M. J. Wainwright, and K. Ramchandran, “Information-theoretic limits on sparse support recovery: Dense versus sparse measurements,” in _2008 IEEE International Symposium on Information Theory_. IEEE, 2008, pp. 2197–2201.
  * [31] J. Haupt and R. Baraniuk, “Robust support recovery using sparse compressive sensing matrices,” in _2011 45th Annual Conference on Information Sciences and Systems_. IEEE, 2011, pp. 1–6.
  * [32] J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and L. Fei-Fei, “Imagenet: A large-scale hierarchical image database,” in _2009 IEEE conference on computer vision and pattern recognition_. Ieee, 2009, pp. 248–255.
  * [33] X. Wang, Y. Peng, L. Lu, Z. Lu, M. Bagheri, and R. M. Summers, “Chestx-ray8: Hospital-scale chest x-ray database and benchmarks on weakly-supervised classification and localization of common thorax diseases,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2017, pp. 2097–2106.
  * [34] S. Shekhar, V. M. Patel, N. M. Nasrabadi, and R. Chellappa, “Joint sparse representation for robust multimodal biometrics recognition,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 36, no. 1, pp. 113–126, 2013.
  * [35] X. Mei and H. Ling, “Robust visual tracking and vehicle classification via sparse representation,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 33, no. 11, pp. 2259–2272, 2011.
  * [36] T. Guha and R. K. Ward, “Learning sparse representations for human action recognition,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 34, no. 8, pp. 1576–1588, 2011.
  * [37] W. Li and Q. Du, “A survey on representation-based classification and detection in hyperspectral remote sensing imagery,” _Pattern Recognition Letters_ , vol. 83, pp. 115–123, 2016.
  * [38] 2020\. [Online]. Available: https://radiopaedia.org/playlists/25975?lang=us
  * [39] 2020\. [Online]. Available: https://threadreaderapp.com/thread/1243928581983670272.html
  * [40] P. Mooney, “Chest x-ray images (pneumonia),” _kaggle, Marzo_ , 2018. [Online]. Available: https://www.kaggle.com/paultimothymooney/chest-xray-pneumonia
  * [41] K. R. Rad, “Nearly sharp sufficient conditions on exact sparsity pattern recovery,” _IEEE Transactions on Information Theory_ , vol. 57, no. 7, pp. 4672–4679, 2011.
  * [42] J. Scarlett and V. Cevher, “Limits on support recovery with probabilistic models: An information-theoretic framework,” _IEEE Transactions on Information Theory_ , vol. 63, no. 1, pp. 593–620, 2016.
  * [43] G. Reeves and M. Gastpar, “Sampling bounds for sparse support recovery in the presence of noise,” in _2008 IEEE International Symposium on Information Theory_. IEEE, 2008, pp. 2187–2191.
  * [44] G. Reeves and M. C. Gastpar, “Approximate sparsity pattern recovery: Information-theoretic lower bounds,” _IEEE Transactions on Information Theory_ , vol. 59, no. 6, pp. 3451–3465, 2013.
  * [45] A. Değerli, S. Aslan, M. Yamac, B. Sankur, and M. Gabbouj, “Compressively sensed image recognition,” in _2018 7th European Workshop on Visual Information Processing (EUVIP)_. IEEE, 2018, pp. 1–6.
  * [46] S. Lohit, K. Kulkarni, and P. Turaga, “Direct inference on compressive measurements using convolutional neural networks,” in _2016 IEEE International Conference on Image Processing (ICIP)_ , Sep. 2016, pp. 1913–1917.
  * [47] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * [48] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin _et al._ , “Tensorflow: Large-scale machine learning on heterogeneous distributed systems,” _arXiv preprint arXiv:1603.04467_ , 2016.
  * [49] K. Kulkarni, S. Lohit, P. Turaga, R. Kerviche, and A. Ashok, “Reconnet: Non-iterative reconstruction of images from compressively sensed measurements,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2016, pp. 449–458.


# Expert-guided Regularization
via Distance Metric Learning

Shouvik Mani Mehdi Maasoumy Sina Pakazad Henrik Ohlsson

C3.ai
Redwood City, CA
<EMAIL_ADDRESS>

###### Abstract

High-dimensional prediction is a challenging problem setting for traditional
statistical models. Although regularization improves model performance in high
dimensions, it does not sufficiently leverage knowledge on feature importances
held by domain experts. As an alternative to standard regularization
techniques, we propose Distance Metric Learning Regularization (DMLreg), an
approach for eliciting prior knowledge from domain experts and integrating
that knowledge into a regularized linear model. First, we learn a Mahalanobis
distance metric between observations from pairwise similarity comparisons
provided by an expert. Then, we use the learned distance metric to place prior
distributions on coefficients in a linear model. Through experimental results
on a simulated high-dimensional prediction problem, we show that DMLreg leads
to improvements in model performance when the domain expert is knowledgeable.

##  1 Introduction

Machine learning models exhibit poor predictive performance in high-
dimensional problems, where the number of features is relatively large
compared to the number of observations. In the high-dimensional setting,
models suffer from high variance and overfit to the training data, a
phenomenon known as the curse of dimensionality. Most methods used to mitigate
this issue rely solely on data and fail to leverage knowledge from domain
experts. In this paper, we propose a method that elicits an expert’s knowledge
on feature importances and integrates it into a regularized linear model to
improve performance in high-dimensional settings.

High-dimensional datasets are ubiquitous in practice. Certain types of data
such as images, text, and time series are inherently high-dimensional. In many
domains, features are often cheap and numerous – for instance, one may collect
measurements from many sensors without knowing which ones are relevant to the
prediction task at hand. In contrast, observations are often expensive to
obtain and label. As a result, many problems have a large number of features,
but relatively few observations. Under such circumstances, classical
approaches are no longer valid (Donoho, 2000).

Given the wide presence of high-dimensional data, techniques such as
regularization and dimensionality reduction are often used to fit simpler,
less flexible models. These methods reduce the variance of the model at the
cost of increased bias. For example, ridge regression and lasso are
regularized variants of least squares linear regression, and they reduce the
model’s variance by shrinking coefficient estimates toward zero (Hoerl &
Kennard, 1970; Tibshirani, 1996).

Figure 1: Regularization is effective for high-dimensional problems: ridge
regression and lasso (with optimal $\lambda$ chosen by cross-validation)
outperform k-nearest neighbors and linear regression. Still, even ridge
regression and lasso suffer from overfitting as the number of noise features
increases.

Before introducing our approach, it is worth discussing the strengths and
weaknesses of regularization. In Figure 1, we compare the performance of
k-nearest neighbors, linear regression, ridge regression, and lasso on a
simulated regression problem where there are $n\in\\{10,25,50,75,100\\}$
features, out of which only the first $10$ features are truly associated with
the response. The remaining $n-10$ features are noise features with true
coefficient values equal to zero. The training set consists of $m=100$
observations and the validation MSE is evaluated on a validation set of $900$
observations111For complete details on the data generation procedure, see
Section 4.. Cross-validation is used for the ridge regression and lasso models
to select the optimal regularization parameter from
$\lambda\in\\{10^{-4},10^{-3},...,10^{2}\\}$ at each $n$. As $n$ approaches
$m$, the performance of k-nearest neighbors and linear regression rapidly
deteriorate, but that of ridge regression and lasso remain relatively stable.

Still, even ridge regression and lasso suffer from overfitting, indicated by a
steady increase in their validation MSEs as $n$ approaches $m$. With enough
noise features, even regularization does not prevent overfitting, since chance
associations between the features and the response on the training set result
in some noise features being assigned nonzero coefficient estimates, even
though those features are not truly associated with the response (James et
al., 2014). Moreover, truly associated features with high true coefficient
values have their coefficient estimates driven toward zero due to the
regularization penalty.

To better understand these limitations and motivate our approach, it is
helpful to look at regularization from a Bayesian perspective. Ridge
regression and lasso have simple Bayesian interpretations (Murphy, 2012; Keng,
2016):

  * •

Ridge regression is the maximum a posteriori (MAP) solution from assuming a
standard linear model with a Gaussian prior on the coefficients, with mean
zero and variance parameterized by $\lambda$.

  * •

Lasso is the MAP solution from assuming a standard linear model with a Laplace
prior on the coefficients, with mean zero and scale parameterized by
$\lambda$.

In both cases, all model coefficients share the same prior distribution, since
they are all centered at zero and parameterized by the same regularization
parameter $\lambda$. If we expect coefficients to be small or sparse in
general, then ridge regression and lasso provide effective ways to encode that
prior knowledge into the model. But if we have more informative prior
knowledge about individual feature importances, then ridge regression and
lasso do not help to encode that knowledge.

We propose a new regularization technique called Distance Metric Learning
Regularization (DMLreg) to elicit prior knowledge on feature importances and
incorporate that knowledge into a linear (classification or regression) model.
In situations where training data are limited yet domain experts possess a
wealth of prior knowledge about the problem, DMLreg combines human-driven
priors with data-driven parameter estimation to fit a better regularized
model.

##  2 Related work

This paper relies on distance metric learning, a concept pioneered by Xing et
al. (2002) to automatically learn a distance metric from data and user
guidance. Distance metric learning algorithms optimize a distance function,
most commonly the Mahalanobis distance, to respect the user’s notion of
similarity between instances. As input to the algorithm, the user provides
weak supervision through pairs of similar and dissimilar instances (e.g.
“$x_{i}$ is similar/dissimilar to $x_{j}$”). Alternatively, the user can
provide relative comparisons between instances (e.g. “$x_{i}$ is more similar
to $x_{j}$ than to $x_{k}$”) (Schultz & Joachims, 2003). Surveys by Kulis
(2013) and Bellet et al. (2013) review the vast literature on distance metric
learning.

Typically, the learned distance metric is used as an alternative to Euclidean
distance to improve the performance of nearest-neighbor or kernel regression
methods (Weinberger et al., 2006; Weinberger & Tesauro, 2007); however, in
this paper, we integrate the learned metric into a regularized linear model.
This offers advantages over nearest-neighbor methods such as greater
interpretability, robustness to an incorrect learned distance metric, and a
better model fit when a linear model is more appropriate for the underlying
data.

Our work has parallels to research in knowledge elicitation. Daee et al.
(2017) have a similar objective to elicit and incorporate expert knowledge,
but they do so by using expected information gain to identify the most
important features on which to query expert feedback. Unlike their approach
which elicits knowledge at the feature level, our approach elicits knowledge
at the sample level. This gives us advantages in high-dimensional settings
where individual features may not be interpretable (e.g. pixel intensities in
an image), but comparing samples for similarity may be easier.

More generally, a thorough discussion on the challenges and approaches in
high-dimensional problems is provided by Hastie et al. (2009) and James et al.
(2014). These resources describe the curse of dimensionality in both kernel
methods and linear models. Additionally, Murphy (2012) and Wasserman (2010)
offer good references for Bayesian inference, a key framework used in this
paper.

##  3 Methodology

Our DMLreg regularization approach has two key steps: (i) elicit domain
knowledge on feature importances and (ii) incorporate that knowledge into a
linear model. We describe this process in detail and derive MAP estimates for
two models (linear regression and logistic regression) combined with DMLreg.

###  3.1 Eliciting domain knowledge

The first objective of DMLreg is to elicit knowledge on feature importances
held by the domain expert. Here, one may ask: why not simply ask the expert
which features are relevant to the task at hand? While ideal, this is
impractical in the high-dimensional setting – it may be overwhelming or even
impossible for an expert to rank or select good features in a feature space
with hundreds or thousands of features.

In contrast, in many problems, it is relatively easy for an expert to look at
pairs of observations and determine whether they are similar or dissimilar. In
DMLreg, we exploit this ability of experts to do pairwise similarity
comparisons in order to learn their tacit knowledge about feature importances.
In particular, we aim to learn a weighting of the features that reflect their
relative importances, and we achieve this task through distance metric
learning.

We use the original distance metric learning formulation (Xing et al., 2002),
through which we learn a Mahalanobis distance metric $d_{A}(x_{i},x_{j})$
between points $x_{i}$ and $x_{j}$ that respects the expert’s notion of
similarity between observations. In other words, the learned distance metric
assigns small distances between observations that the expert considers similar
and large distances between observations that the expert considers dissimilar.

$\displaystyle d_{A}(x_{i},x_{j})=\sqrt{(x_{i}-x_{j})^{T}A(x_{i}-x_{j})}$ (1)

To learn this distance metric, the domain expert must provide weak supervision
on pairwise similarities in the form of sets $S$ and $D$, each containing
pairs of observations which the expert considers similar and dissimilar,
respectively.

$\displaystyle S$ $\displaystyle=\\{(x_{i},x_{j}):x_{i}\text{ and }x_{j}\text{
should be similar}\\}$ (2) $\displaystyle D$
$\displaystyle=\\{(x_{i},x_{j}):x_{i}\text{ and }x_{j}\text{ should be
dissimilar}\\}$ (3)

Learning the distance metric can then be posed as a convex optimization
problem: find an optimal matrix $A$ to minimize the sum of squared distances
between the pairs of points in set $S$ that the expert has identified as
similar.

$\displaystyle\underset{A}{\text{minimize}}\quad$
$\displaystyle\textstyle\sum_{(x_{i},x_{j})\in S}\|x_{i}-x_{j}\|^{2}_{A}$ (4)
subject to $\displaystyle\textstyle\sum_{(x_{i},x_{j})\in
D}\|x_{i}-x_{j}\|_{A}\geq 1,$ (5) $\displaystyle A\succeq 0,$ (6)
$\displaystyle A_{ij}=0,\quad i\neq j.$ (7)

The constraint in (5) maintains a lower bound on the total of distances
between the dissimilar points in $D$ and prevents the trivial solution to (4)
of $A=0$. Constraint (6) requires $A$ to be a positive semi-definite matrix so
that it is a valid distance metric that satisfies non-negativity and triangle
inequality. Finally, constraint (7) restricts $A$ to be a diagonal matrix,
which will make it easier to integrate $A$ with the model.

The diagonal elements of the learned matrix $A$ represent a weighting of the
features based on their relative importances as indicated by the expert (at
least with regard to computing distances between observations). For
simplicity, we will refer to $A$ as the “distance metric” in the following
discussion.

###  3.2 Incorporating domain knowledge into the model

Once we have captured the expert’s knowledge on feature importances in the
form of a learned distance metric $A$, DMLreg integrates this knowledge into a
regularized linear model. A natural way to combine prior knowledge with data
is through Bayesian inference. In Bayesian inference, we specify a prior
probability distribution $P(\theta)$ and a likelihood function $P(y|\theta)$.
Using Bayes’ theorem, we can write the posterior distribution $P(\theta|y)$ as
proportional to the likelihood times the prior.

$\displaystyle P(\theta|y)=\frac{P(y|\theta)P(\theta)}{P(y)}\propto
P(y|\theta)P(\theta)$ (8)

Specifically, let us first consider the case of Bayesian linear regression,
which assumes a linear model $y=\theta^{T}x+\epsilon$ with Gaussian errors
$\epsilon$. This implies that $y$ is also a Gaussian centered at the
regression line with a corresponding likelihood function $P(y|\theta,x)$.

$\displaystyle y$
$\displaystyle=\theta^{T}x+\epsilon,\quad\epsilon\sim\mathcal{N}(0,\sigma^{2})$
(9) $\displaystyle\iff y|\theta,x\sim\mathcal{N}(\theta^{T}x,\sigma^{2})$ (10)
$\displaystyle P(y|\theta,x)$
$\displaystyle=\frac{1}{\sqrt{2\pi\sigma^{2}}}exp\Big{(}\frac{-(y-\theta^{T}x)^{2}}{2\sigma^{2}}\Big{)}$
(11)

The prior distribution $P(\theta)$ in Bayesian linear regression lets us
express our knowledge about the coefficients $\theta$. DMLreg encodes the
knowledge held in the learned distance metric $A$ through a Gaussian prior on
each coefficient. Like ridge regression, the Gaussian priors are all centered
at zero, since $A$ does not provide any information on the direction of the
relationship between each feature and the response.

However, the diagonal elements of $A$ represent the relative importances of
the features. For $i\in 1,…,n$, if $A_{ii}$ is large, then we can expect a
priori that feature $i$ has a large positive or negative coefficient
$\theta_{i}$; in other words, we can expect $\theta_{i}$ to have a large
variance around zero. So, we encode this information through the variance of
each coefficient’s prior. This leads to placing a Gaussian prior distribution
on each coefficient $\theta_{i}$ with mean zero and variance $A_{ii}$, and a
corresponding density function $P(\theta_{i})$.

$\displaystyle\theta_{i\in 1,...,n}$ $\displaystyle\sim\mathcal{N}(0,A_{ii})$
(12) $\displaystyle P(\theta_{i})$ $\displaystyle=\frac{1}{\sqrt{2\pi
A_{ii}^{2}}}exp\Big{(}\frac{-\theta_{i}^{2}}{2A_{ii}}\Big{)}$ (13)

With our likelihood function $P(y|\theta,x)$ and prior distributions
$P(\theta_{i})$ defined, we can use Bayes’ theorem to calculate the posterior
distribution of $\theta$. Since we are only interested in the coefficient
point estimates and not the full distribution of $\theta$, we can use maximum
a posteriori (MAP) estimation to find the posterior mode $\hat{\theta}_{MAP}$:
the value of $\theta$ that is most likely given the data and the priors.

$\displaystyle\hat{\theta}_{MAP}$
$\displaystyle=\operatorname*{argmax}_{\theta}P(\theta|y_{1},...,y_{m},x_{1},...,x_{m})$
(14)
$\displaystyle=\operatorname*{argmax}_{\theta}\frac{P(y_{1},...,y_{m}|\theta,x_{1},...,x_{m})P(\theta)}{P(y_{1},...,y_{m})}$
(15)
$\displaystyle=\operatorname*{argmax}_{\theta}\prod_{i=1}^{m}P(y_{i}|\theta,x_{i})\prod_{j=1}^{n}P(\theta_{j})$
(16)
$\displaystyle=\operatorname*{argmax}_{\theta}\sum_{i=1}^{m}logP(y_{i}|\theta,x_{i})+\sum_{j=1}^{n}logP(\theta_{j})$
(17)
$\displaystyle=\operatorname*{argmax}_{\theta}\sum_{i=1}^{m}\frac{-(y_{i}-\theta^{T}x_{i})^{2}}{2\sigma^{2}}+\sum_{j=1}^{n}\frac{-\theta_{j}^{2}}{2A_{jj}}$
(18)
$\displaystyle=\operatorname*{argmax}_{\theta}-\frac{1}{2\sigma^{2}}\sum_{i=1}^{m}(y_{i}-\theta^{T}x_{i})^{2}-\frac{1}{2}\sum_{j=1}^{n}\frac{\theta_{j}^{2}}{A_{jj}}$
(19)
$\displaystyle=\operatorname*{argmin}_{\theta}\sum_{i=1}^{m}(y_{i}-\theta^{T}x_{i})^{2}+\sum_{j=1}^{n}\frac{\theta_{j}^{2}}{A_{jj}}$
(20)
$\displaystyle=\operatorname*{argmin}_{\theta}\|y-X\theta\|_{2}^{2}+\theta^{T}A^{-1}\theta$
(21)

Line (21) is the loss function for linear regression with DMLreg (Gaussian
prior). We can solve for $\hat{\theta}$ to obtain a closed-form solution for
the coefficient estimates.

$\displaystyle\nabla_{\hat{\theta}}\Big{(}\|y-X\hat{\theta}\|_{2}^{2}+\hat{\theta}^{T}A^{-1}\hat{\theta}\Big{)}=0$
(22) $\displaystyle-2X^{T}(y-X\hat{\theta})+2A^{-1}\hat{\theta}=0$ (23)
$\displaystyle-X^{T}y+X^{T}X\hat{\theta}+A^{-1}\hat{\theta}=0$ (24)
$\displaystyle(X^{T}X+A^{-1})\hat{\theta}=X^{T}y$ (25)
$\displaystyle\hat{\theta}=(X^{T}X+A^{-1})^{-1}X^{T}y$ (26)

Thus, for linear regression with DMLreg (Gaussian prior), we can compute
coefficient estimates $\hat{\theta}$ in terms of $X,y$, and the learned
distance metric $A$. This is equivalent to the coefficient estimates for ridge
regression when $A=\lambda I$, for regularization parameter $\lambda$. An
interpretation of this result is that linear regression with DMLreg (Gaussian
prior) is a generalization of ridge regression with a separate regularization
parameter $A_{ii}$ for each feature $i$, instead of a single regularization
parameter $\lambda$ across all features.

Apart from linear regression, DMLreg can be easily applied to any generalized
linear model. In Table 1, we show how to fit linear regression and logistic
regression models with DMLreg for both Gaussian and Laplace priors on the
coefficients. 222For the Logistic Regression + DMLreg models in Table 1,
$\sigma(\cdot)$ is the Sigmoid function $\sigma(x)=\frac{1}{1+exp(-x)}$.

Table 1: Fitting linear models with DMLreg using Gaussian and Laplace priors on coefficients. The diagonal matrix $A$ represents the learned distance metric from the first step. Model Name | Assumptions | Optimization Method
---|---|---
Linear Regression \+ DMLreg (Gaussian prior) |  $\theta_{i}\sim\mathcal{N}(0,A_{ii})$ $y\sim\mathcal{N}(\theta^{T}x,\sigma^{2})$ |  Analytic solution: $\theta=(X^{T}X+A^{-1})^{-1}X^{T}y$
Linear Regression \+ DMLreg (Laplace prior) |  $\theta_{i}\sim Laplace(0,A_{ii})$ $y\sim\mathcal{N}(\theta^{T}x,\sigma^{2})$ |  Subgradient method: $\operatorname*{argmin}_{\theta}\|y-X\theta\|_{2}^{2}+\|A^{-1}\theta\|_{1}$
Logistic Regression \+ DMLreg (Gaussian prior) |  $\theta_{i}\sim\mathcal{N}(0,A_{ii})$ $y\sim Bern(p=\sigma(\theta^{T}x))$ |  Gradient descent: $\theta:=\theta+\alpha[y-\sigma(X\theta)]X-\alpha A^{-1}\theta$
Logistic Regression \+ DMLreg (Laplace prior) |  $\theta_{i}\sim Laplace(0,A_{ii})$ $y\sim Bern(p=\sigma(\theta^{T}x))$ |  Subgradient method: $\operatorname*{argmin}_{\theta}-ylog\sigma(X\theta)-$ $(1-y)log[1-\sigma(X\theta)]+\|A^{-1}\theta\|_{1}$

##  4 Results

We evaluate DMLreg through an experiment on an artificial dataset with
simulated domain knowledge. The experiment is outlined in Figure 2.

Figure 2: Experiment used to evaluate DMLreg on artificially-generated data
$X,y$ using simulated domain knowledge in $A$.

First, we generate artificial high-dimensional datasets with $m=100$
observations and $n\in\\{10,25,50,75,100\\}$ features. We sample an
$n$-dimensional coefficient vector $\theta$ with the first 10 coefficients
drawn from a Uniform distribution between -10 and 10 and the remaining $n-10$
coefficients set exactly to zero (line 27). Keeping $\theta$ fixed, we
simulate 100 datasets $X,y$ through the procedure in lines 28-30. To be
precise, our training set has $m=100$ observations; we hold out a separate
validation set with 900 observations.

$\displaystyle\theta$
$\displaystyle\in\mathbb{R}^{n},\quad\theta_{1,...,10}\sim\textit{Unif}(-10,10),\quad\theta_{11,...n}=0$
(27) $\displaystyle X$ $\displaystyle\in\mathbb{R}^{m\times n},\quad
X\sim\textit{Unif}(0,1)$ (28) $\displaystyle\epsilon$
$\displaystyle\in\mathbb{R}^{m},\quad\epsilon\sim\mathcal{N}(0,1)$ (29)
$\displaystyle y$ $\displaystyle\in\mathbb{R}^{m},\quad y=X\theta+\epsilon$
(30)

To simulate an expert’s domain knowledge on feature importances, we create
three true distance metrics $A$, listed in Table 2, which represent varying
levels of expertise: perfect, noisy, and incorrect. These true metrics $A$ are
used for computing pairwise distances between observations in order to
generate sets $S$ and $D$ for distance metric learning (to simulate an expert
manually comparing observations and creating these sets).

Table 2: True metrics used in the experiment to simulate domain knowledge on
feature importances with varying levels of expertise.
Knowledge type | True metric value
---|---
Perfect | $A=\textit{diag}(|\theta_{1}|,...,|\theta_{n}|)$
Noisy |  |  $A=\textit{diag}(|\theta_{1}+\epsilon_{1}|,...,|\theta_{n}+\epsilon_{n}|)$, $\epsilon_{i}\sim\mathcal{N}(0,0.5)$
---
Incorrect |  |  $A=\textit{diag}(|\alpha_{1}|,...,|\alpha_{n}|)$, $\alpha_{i}\sim\textit{Unif}(-5,5)$
---

In practice, we expect real-world expert knowledge to be similar to a noisy or
incorrect knowledge metric. The noisy knowledge metric is a diagonal matrix
which weighs each feature by the magnitude of its true coefficient value
$\theta_{i}$ plus some error $\epsilon_{i}$. It represents an expert who knows
the correct feature importances up to some error. In contrast, the incorrect
knowledge metric weighs the features with completely random weights
$\alpha_{i}$.

Now, we explain how to generate sets $S$ and $D$ in the experiment. For a
feature matrix $X$ and a true metric $A$, we compute pairwise distances
between all observation pairs $(x_{i},x_{j})$ in $X$ using the Mahalanobis
distance $d_{A}(x_{i},x_{j})=\sqrt{(x_{i}-x_{j})^{T}A(x_{i}-x_{j})}$. The $p$
pairs with the highest distances are added to set $D$, and the $p$ pairs with
the lowest distances are added to set $S$.

The sets $S$ and $D$ are then passed on to the distance metric learning
algorithm to learn an approximation $\hat{A}$ for $A$. The learned metric
$\hat{A}$, along with $X$ and $y$ are used as inputs to DMLreg, which fits a
regularized linear model.

We evaluate the performance of the Linear Regression + DMLreg (Laplace prior)
model on the validation set for each learned metric $\hat{A}$ and for each
$n\in\\{10,25,50,75,100\\}$ features. Beside the metrics $\hat{A}$ learned to
approximate the three true metrics $A$, we also consider the Euclidean
distance metric $A=I$, which is identical to the lasso ($\lambda=1$) model
when used in DMLreg with a Laplace prior. For this evaluation, we fix $p=300$
pairs. The validation set MSEs are aggregated across the 100 $X,y$ datasets
and displayed in Figure 3.

Figure 3: Performance of Linear Regression + DMLreg (Laplace prior) model
using various learned distance metrics as the number of features increases,
for a fixed $m=100$ observations.

For high values of $n$, DMLreg using metrics learned from perfect or noisy
knowledge performs better than DMLreg using the Euclidean distance metric
(lasso with $\lambda=1$). DMLreg substantially outperforms k-nearest neighbors
and linear regression evaluated on the same data (refer to Figure 1). It also
performs slightly better than the lasso model which had $\lambda$ tuned
through grid search and cross-validation (refer to Figure 1), even though no
hyperparameter search was required for its DMLreg counterpart. Nevertheless,
when provided a metric learned from incorrect knowledge, DMLreg performs worse
than lasso. These results suggest that DMLreg is effective for high-
dimensional problems as long as the domain knowledge provided through the
learned metric is reasonably correct.

Figure 4: DMLreg performance with $n=100$ features against the number of
observation pairs $p$.

An important practical concern when using DMLreg is: how many pairs of
observations $p$ does one need to provide in sets $S$ and $D$ in order to get
good performance? Figure 4 shows the validation MSEs for DMLreg on a single
dataset $X,y$ with a fixed $n=100$ features as $p$ is varied from 25 to 700
observations pairs. For DMLreg using the metric learned from noisy knowledge,
providing just $p=25$ observation pairs yields nearly the same performance as
$p=300$ pairs, as in the experiment above.

To understand why the DMLreg model using the metric learned from noisy
knowledge outperforms lasso, it is helpful to examine the coefficient
estimates in Figure 5. For a single dataset $X,y$, Figure 5 shows the
coefficients for the first 10 (relevant) features on the left and the
remaining (noise) features on the right. It is clear that the coefficient
estimates for the DMLreg model (with noisy knowledge and Laplace prior) are
closer to the true coefficients than the coefficient estimates for the lasso
($\lambda=1$) model. In fact, the lasso model slightly underestimates the
coefficients of the truly associated features and overestimates the
coefficients of the noise features when compared to the DMLreg model.

Figure 5: True coefficients vs. coefficient estimates for the Linear
Regression + DMLreg (Laplace prior) model and the Lasso ($\lambda=1$) model.

##  5 Conclusion

When training data is limited but domain knowledge on feature importances is
available, DMLreg can be used to fit a better regularized model. Using
distance metric learning and regularization, DMLreg elicits and integrates
expert knowledge into a linear model. Through an experiment on artificial data
using simulated domain knowledge, we demonstrated that DMLreg outperforms
ridge regression and lasso when the knowledge elicited is approximately
correct.

## References

  * Bellet et al. (2013) Bellet, A., Habrard, A., and Sebban, M.  A survey on metric learning for feature vectors and structured data.  _CoRR_ , abs/1306.6709, 2013.
  * Daee et al. (2017) Daee, P., Peltola, T., Soare, M., and Kaski, S.  Knowledge elicitation via sequential probabilistic inference for high-dimensional prediction.  _Machine Learning_ , 106(9):1599–1620, Oct 2017\.
  * Donoho (2000) Donoho, D. L.  High-dimensional data analysis: The curses and blessings of dimensionality.  In _AMS Conference on Math Challenges of the 21st Century_ , 2000\.
  * Hastie et al. (2009) Hastie, T., Tibshirani, R., and Friedman, J.  _The Elements of Statistical Learning: Data Mining, Inference, and Prediction_.  Springer, 2009.
  * Hoerl & Kennard (1970) Hoerl, A. E. and Kennard, R. W.  Ridge regression: Biased estimation for nonorthogonal problems.  _Technometrics_ , 12(1):55–67, 1970.
  * James et al. (2014) James, G., Witten, D., Hastie, T., and Tibshirani, R.  _An Introduction to Statistical Learning: With Applications in R_ , pp. 238–244.  Springer, 2014.
  * Keng (2016) Keng, B.  A probabilistic interpretation of regularization, 2016.
  * Kulis (2013) Kulis, B.  _Metric Learning: A Survey_.  Foundations and trends in machine learning. Now Publishers, 2013.
  * Murphy (2012) Murphy, K. P.  _Machine Learning: A Probabilistic Perspective_.  The MIT Press, 2012.
  * Schultz & Joachims (2003) Schultz, M. and Joachims, T.  Learning a distance metric from relative comparisons.  In _Proceedings of the 16th International Conference on Neural Information Processing Systems_ , 2003.
  * Tibshirani (1996) Tibshirani, R.  Regression shrinkage and selection via the lasso.  _Journal of the Royal Statistical Society. Series B (Methodological)_ , 58(1):267–288, 1996.
  * Wasserman (2010) Wasserman, L.  _All of Statistics: A Concise Course in Statistical Inference_.  Springer, 2010.
  * Weinberger & Tesauro (2007) Weinberger, K. Q. and Tesauro, G.  Metric learning for kernel regression.  In _Artificial Intelligence and Statistics_ , pp. 612–619, 2007\.
  * Weinberger et al. (2006) Weinberger, K. Q., Blitzer, J., and Saul, L. K.  Distance metric learning for large margin nearest neighbor classification.  In _Advances in neural information processing systems_ , 2006.
  * Xing et al. (2002) Xing, E. P., Ng, A. Y., Jordan, M. I., and Russell, S.  Distance metric learning, with application to clustering with side-information.  In _Proceedings of the 15th International Conference on Neural Information Processing Systems_ , 2002.