# One-shot Representational Learning for Joint Biometric and Device
Authentication

Sudipta Banerjee and Arun Ross  Michigan State University
Email: {banerj24<EMAIL_ADDRESS>
###### Abstract

In this work, we propose a method to simultaneously perform (i) biometric
recognition (i.e., identify the individual), and (ii) device recognition,
(i.e., identify the device) from a single biometric image, say, a face image,
using a one-shot schema. Such a joint recognition scheme can be useful in
devices such as smartphones for enhancing security as well as privacy. We
propose to automatically learn a joint representation that encapsulates both
biometric-specific and sensor-specific features. We evaluate the proposed
approach using iris, face and periocular images acquired using near-infrared
iris sensors and smartphone cameras. Experiments conducted using 14,451 images
from 15 sensors resulted in a rank-1 identification accuracy of upto 99.81%
and a verification accuracy of upto 100% at a false match rate of 1%.

##  I Introduction

Biometric data such as face, fingerprint or iris images reveal information
about the identity of the individual as well as the identity of the device
used to acquire the data [1, 2]. In some applications such as smartphone
banking, it is necessary to authenticate both the user as well as the device
in order to enhance security [3, 4]. This can be done by invoking two separate
modules: one for biometric recognition and the other for device or sensor
recognition.111The terms “device” and “sensor” are interchangeably used in
this paper. Thus, determining the identity of a smartphone camera (i.e.,
sensor) is akin to determining the identity of the smartphone (i.e., device)
itself. In such cases, the system has to store two distinct templates: a
biometric template denoting the identity of the user and sensor template
denoting the identity of the device.

In this paper, we approach this problem by designing a joint template that can
be used to authenticate both the user and the device simultaneously. Our
objective is as follows: Given a biometric image we would like to
simultaneously recognize the individual and the acquisition device. In the
process of accomplishing this objective, we address the following questions:

(a)

(b)

Figure 1: Difference between (a) methods that use separate modules for
computing biometric and sensor representations, and (b) the proposed method
that uses an embedding network to generate a joint biometric-sensor
representation.

  1. 1.

Why do we need to combine biometric and device recognition?
Smartphones are increasingly using biometrics for access control and monetary
transactions. Examples include fingerprint and face recognition on iPhones and
iris recognition on Samsung Galaxy S9. Device verification222Typically used in
two factor authentication (2FA) protocol that combines any two of the three
factors: ‘something you are’ (biometrics), ‘something you have’ (a code on the
authorized device) and ‘something you know’ (a password) for additional
security. can provide assurance that the biometric sample is being acquired by
an authorized device. A combined biometric and device recognition system can
therefore guarantee that the right person is accessing the remote service
(e.g., banking) using an authorized device.

  2. 2.

Can existing device verification techniques be used in the smartphone
application scenario?
Device identification can be performed using the MAC (media access control)
address, a unique networking address assigned to each device. However, in case
of smartphones that have multiple network interfaces, such as Wi-Fi, 4G,
bluetooth, etc., there can be multiple MAC addresses which may be broadcasted,
making them vulnerable. Alternatively, SRAM cells can be used to deduce
physically unclonable cues for device identification [4]; this is a hardware-
based solution and requires access to the physical device. In a mobile banking
scenario, where the verification is conducted remotely, the customer provides
a biometric sample in the form of an image, and some device information, but
not necessarily the physical device itself. In this scenario, hardware-based
solutions will be ineffective.

  3. 3.

Why do we need a joint representation?
Existing literature uses separate modules to tease out the biometric-specific
and sensor-specific details from an image and perform feature-level or score-
level fusion [3, 4]. However, they suffer from the following limitations: (i)
the overall performance is limited by the weakest recognition module, and (ii)
the process may not generalize well across different biometric modalities and
multi-spectral sensors. Therefore, a joint representation that combines both
biometric and sensor-specific features present in a biometric image can offer
the following advantages: (i) the joint representation is not constrained by
the performance of the individual recognition module, and the same method can
be employed across different biometric modalities, and (ii) the joint
representation integrates the biometric and sensor representations into a
compact template, such that, the individual templates cannot be easily de-
coupled; this implicitly imparts privacy to the biometric component.

##  II Related Work

Biometric recognition systems comprise a feature extraction module that
elicits a salient feature representation from the acquired biometric data, and
a comparator module that compares two sets of feature representations to
compute a match score [5]. On the other hand, sensor recognition systems
extract sensor pattern noise [6] from a set of training images obtained from
different sensors to generate sensor reference patterns. To deduce the sensor
identity of an unknown test image, first its sensor pattern noise is
extracted, and then it is correlated with the reference patterns. The test
image is assigned to the device whose reference pattern yields the highest
correlation value.

In [7], the authors used partial face images acquired using smartphones and
employed a weighted sum fusion rule at the score level to combine sensor and
biometric recognition. Later, they extended their work to include feature
level fusion in [3] and concluded that score level fusion performed
comparatively better. In [8], the authors performed HOG-based face recognition
and combined it with Photo Response Non-Uniformity-based sensor recognition at
the score level. In [4], the authors combined fingerprint recognition with
device recognition by performing feature level fusion of minutiae-cylinder-
codes with SRAM start-up values. Fusion at the score or feature level is often
dependent on the specific biometric modality and the device sensor used. A
specific fusion rule producing the best results on a particular biometric and
sensor modality (e.g, iris and near-infrared sensors) may not yield optimal
results on a different modality (e.g, face and RGB sensors), and therefore,
needs to be tuned separately for each pair of biometric and sensor modalities.
Furthermore, feature-level fusion retains the individual biometric and sensor-
specific components that can be recovered from the fused representation using
appropriate measures. Obtaining the biometric template may compromise the
privacy aspect of biometrics. In contrast, the proposed joint representation
non-trivially unifies the biometric and sensor-specific features. As a result,
typical countermeasures will be ineffective in disentangling the biometric
component from the joint representation. This will implicitly preserve the
privacy of the biometric component.

The remainder of the paper is organized as follows. Section III describes the
proposed method. Section IV describes the dataset and experimental protocols
used in this work. Section V reports the results. Section VI summarizes the
findings and concludes the paper.

##  III Motivation and Proposed Method

An image contains both low frequency and high frequency components. For
example, in a face image, the low frequency components capture the
illumination details while the high frequency components capture the
structural details present in the face that are useful for biometric
recognition. Recently, sensor recognition has been successully accomplished
using Photo Response Non-Uniformity (PRNU) [6] for different types of sensors,
such as DSLR sensors [9], smartphone sensors [10], and also near-infrared iris
sensors [11]. PRNU is a form of sensor pattern noise in an image that
manifests due to anomalies during the fabrication process and is, therefore,
unique to each sensor. Typically, PRNU resides in the high frequencies that
can be useful for sensor recognition [6]. Since the high frequencies dominate
in both biometric and sensor representations, we hypothesize that there is a
joint representation, that, if effectively extracted, can be utilized for both
tasks of biometric and sensor recognition. Our objective is to learn this
joint representation that lies at the intersection of the sensor and biometric
space. Mathematically, it can be represented as $\displaystyle
J(\bm{X})=B(\bm{X})\cap S(\bm{X})$, where $\bm{X}$ is an input biometric
image, $B(\cdot)$ is the biometric representation extracted from $\bm{X}$,
$S(\cdot)$ is the sensor representation computed from the same input $\bm{X}$,
and $J(\cdot)$ is the joint representation. Existing methods process $\bm{X}$
using two independent routines to extract the two representations, and can
optionally perform fusion, either at feature level or at score level, to make
a decision. However, we propose to leverage the two representations to derive
a joint representation (see Figure 1). The joint space can be best
approximated using an embedding network that can convert images to compact
representations [12]. The embedding network $\mathscr{E}$, takes two inputs,
$\bm{X}$ and the dimensionality ($k$) of the embedding to be generated, such
that $\displaystyle J(\bm{X})=\mathscr{E}(\bm{X},k)\approx B(\bm{X})\cap
S(\bm{X})$. The second argument $k$, allows us to regulate the dimensionality
of the joint representation, which will be much lesser than the original
dimensionality of the image, as well as the combined dimensionality of the two
representations computed separately, i.e., if $\bm{X}\in\mathbb{R}^{d}$,
$B(\bm{X})\in\mathbb{R}^{m}$ and $S(\bm{X})\in\mathbb{R}^{n}$, then the joint
representation $J(\bm{X})\in\mathbb{R}^{k}$, where, $k<<d$ and $k<(m+n)$.

Figure 2: Outline of the proposed method used for computing the joint
biometric and sensor representation. Input: A single image, or a pair of
images, or 3-tuple images to the embedding network. Output: Joint biometric-
sensor representation. The embedding network is trained in three mutually
exclusive modes, viz., classical mode (top row), siamese mode (middle row) and
triplet mode (bottom row). The switching circuit selects only one training
mode at a time.

In this work, we used a deep convolutional neural network that serves the role
of the embedding network (see Figure 2). The embedding network consists of two
2-D convolutional layers and three linear layers. We used max-pooling for
down-sampling the feature map and a parametric-rectified linear activation
unit (PReLU) as the activation function. The embedder accepts an image,
resized to $48\times 48$, as the input and produces a 8-dimensional output,
which is the joint representation. The choice of the dimensionality of the
representation along with the experimental setup is described later (see
Section IV-C).

The main contributions of this work are as follows:

  1. 1.

We propose a method to learn a joint biometric and sensor representation using
a one-shot approach that can be used in joint identification and joint
verification scenarios. A correct joint identification/verification occurs
only if both subject identity and device identity yield correct matches.

  2. 2.

We employ an embedding network that can learn the joint representation
irrespective of the biometric modality and the sensor used for acquisition. In
this context, we evaluate the proposed method using three different biometric
modalities (face, iris and periocular), and different types of sensors (iris
sensors operating in the near-infrared spectrum and smartphone camera sensors
operating in the visible spectrum).

  3. 3.

We perform extensive experiments using different training paradigms and loss
functions, and compare the proposed method with existing state-of-the-art
algorithms for biometric and sensor recognition.

##  IV Datasets and Experiments

###  IV-A Datasets

In this work, we focused on three different biometric modalities, viz., iris,
periocular and face. To this end, we used three different datasets - (i)
CASIA-Iris Image Database Version 2 [13] which contains near-infrared iris
images acquired using two sensors, (ii) Mobile Iris Challenge Evaluation
(MICHE-I) dataset [14] which contains partial face images acquired using two
smartphones (front and rear sensors separately) and front camera of a tablet,
and (iii) OULU-NPU dataset [15] which contains face images acquired using the
front sensors of six smartphones. We used only bonafide images from all three
sessions in the OULU-NPU dataset. Table I describes the datasets used in this
work. Note that the smartphone datasets (MICHE-I and OULU-NPU) contain images
acquired in the visible spectrum.

TABLE I: Dataset specifications used in this work. We used three datasets corresponding to 3 biometric modalities viz., iris, periocular and face. Here, we perform joint biometric and sensor recognition, so total $\\#Classes$ is computed as the product of $\\#Subjects$ and $\\#Sensors$. (∗MICHE-I dataset has a total 75 subjects, out of which the first 48 subjects were imaged using iPhone 5S UNIT I and the remaining 27 subjects were imaged using iPhone 5S UNIT II, as observed in [3]. Here, ‘UNIT’ refers to two different units of the same brand and model iPhone 5S, and therefore, should be treated as two different smartphones. In this case, $\\#Classes$ = 375 since only a subset of the total 75 subjects were imaged using either of the two units of iPhone 5S smartphone at a time. Therefore, $75\times 3$ [Samsung] $+48\times 2$ [UNIT I] $+27\times 2$ [UNIT II] $=375$.)  Modality | Dataset |  | Name of sensors
---
| (# Subjects, # Sensors,
---
# Classes)
Split | # Images
Iris | CASIA-Iris V2 | CASIA IrisCAM-V2, OKI IrisPass-h | (60, 2, 120) | Train | 1,680
Test | 720
Periocular | MICHE-I |  Apple iPhone 5S (Front and Rear) UNIT I and UNIT II, Samsung Galaxy S4 (Front and Rear), Samsung Galxy Tab GT2 (Front) | (75, 7, 375∗) | Train | 2,278
Test | 863
Face | OULU-NPU |  HTC Desire EYE, Sony XPERIA C5 Ultra Dual, MEIZU X5, Oppo N3, Samsung Galaxy S6 Edge, ASUS Zenfne Selfie | (55, 6, 330) | Train | 5,940
Test | 2,970
TOTAL | (190, 15, 825) |  | 14,451

###  IV-B Evaluation Protocol

Before we describe the experiments, we present the protocol that is used to
evaluate the proposed approach. We evaluate the method in two scenarios, viz.,
(i) joint identification and (ii) joint verification. The terms joint
identification and joint verification are different from the terms used
conventionally in the biometric literature. In the case of joint
identification, a correct identification occurs only when both sensor and
subject labels of the test sample match with the ground truth labels. To
perform evaluation in the joint identification scenario, we select one
embedding from each class (combines both sensor and subject label) to form the
gallery, and the remaining embeddings are used as probes. We use two metrics
to compute the distance or similarity between the probe and gallery embeddings
and select the top three matches: (i) standardized Euclidean distance
(computes the pairwise euclidean distance divided by the standard deviation)
and (ii) cosine similarity. We plot the cumulative match characteristics (CMC)
curves corresponding to the top three ranks. In the case of joint
verification, two joint representations will yield a match if both the
embeddings belong to the same sensor and same subject, otherwise a mismatch
occurs. Incorrect match can occur in three cases as follows: (i) if the two
joint representations belong to the same subject, but different sensors, (ii)
if the two joint representations belong to the same sensor, but different
subjects, and (iii) if the two joint representations belong to different
subjects and different sensors. To perform evaluation in the joint
verification scenario, we compute the distance or similarity between all the
test embeddings and present receiver operating characteristics (ROC) curves to
indicate the joint verification performance. We also report the true match
rate (TMR) values @1% and 5% false match rates (FMR).

###  IV-C Experimental Settings

In this work, we designed the experimental settings using three different
modes of training. See Figure 2. Say, $\bm{O}$ denotes the output of an
embedding network for input $\bm{X}$, i.e., $\bm{O}=\mathscr{E}(\bm{X},k)$. In
the first mode, referred to as the classical mode, the embedding $\bm{O}$ is
fed to a classification network which minimizes the cross-entropy loss
computed between the ground truth label and the predicted label. The
classification network in our case is a shallow network which applies PReLU
activation on the embedding, followed by a fully-connected layer and then
applies softmax to compute a probability value. We assigned the ground truth
label for the $i^{th}$ image, such that $l_{i}\in Sub_{i}\otimes\ Sen_{i}$,
where $Sub_{i}$ denotes the subject identifier of image $i$, $Sen_{i}$ denotes
the sensor identifier for the same image and $\otimes$ denotes the tensor
product. The cardinality of the set of labels $|L|=|Sub\times Sen|$. In the
second mode, referred to as the siamese mode, a siamese network [16] is used
which feeds a pair of images to the embedding network. The embedding network
then computes a pair of embeddings ($\bm{O}_{i},\bm{O}_{j}$) and the siamese
network is trained by minimizing the contrastive loss [17] computed between
the pair of embeddings. We used single margin (SMCL) and double margin (DMCL)
contrastive losses. Finally, in the third mode, referred to as the triplet
mode, a triplet network [18] is trained using embeddings generated from an
anchor ($\bm{O}_{a}$), a positive ($\bm{O}_{p}$) and a negative ($\bm{O}_{n}$)
sample by minimizing the triplet loss [19]. We performed offline triplet
mining as well as online triplet mining [20] with different triplet selection
strategies (random negative triplet selection, semi hard negative triplet
selection and hardest negative triplet selection). The triplet loss considers
only one negative example at a time. Alternatively, multi-class N-pair loss
function [21] considers multiple negative instances from several classes. In
this work, we consider a positive example as one which belongs to the same
class as the anchor (same subject and same sensor), whereas there can be three
types of negative examples, viz., same subject but different sensor, same
sensor but different subject and different sensor with different subject.
Therefore, the number of negative classes in this work is significantly high,
so we used multi-class N-pair loss using two mining techniques: (i) all
positive pairs and (ii) hard negative pairs. Table II summarizes the different
loss functions used in the three training modes in this work. Note that each
input to the embedding network as shown in Figure 2 is mutually exclusive,
i.e., the embedding network can operate independently in any of the three
training modes. We modified the design of an existing embedding network for
implementing the different training paradigms [22]. We used learning rate =
$1\times\exp{(-4)}$, batch size = 4, Adam optimizer, and a step decay to
reduce the learning rate by a factor $\gamma=0.1$ every 8 epochs. The proposed
network is shallow so we trained only for 50 epochs. The margin values in
single margin contrastive loss and triplet losses are set to 1, while in
double margin contrastive loss, both margins are set to 0.5.

TABLE II: Description of the training modes and the loss functions used in
this work.

Training mode | Loss function
---|---
Classical | Cross entropy
Siamese | Single margin contrastive loss (SMCL)
Double margin contrastive loss (DMCL)
Triplet | Offline triplet mining
Online triplet mining | Random negative
Semi-hard negative
Hardest negative
Multi-class N-pair | All positive pair
Hard negative pair

For each dataset, we used a training set and a test set (see Table I). The
number of classes is computed as the product of the number of sensors and
number of subjects in that dataset. For example, CASIA-Iris V2 dataset has 60
subjects and 2 sensors, so total number of classes is $60\times 2=120$. Each
class has 20 images, therefore, the total number of images (samples) is 2,400
($20\times 120$). The training and test partitions follows a 70:30 split. So,
for a single class, out of 20 samples, 14 samples are randomly selected as the
training set and the remaining 6 samples form the test set. Similar protocol
is followed for the remaining datasets.

Figure 3: Variation in the joint verification performance as a function of the
dimensionality of the joint representation. Experiment is conducted on the
validation set using 50 images from the MICHE-I dataset and four
dimensionality values viz., $\\{4,8,16,32\\}$. 8-dimensional embedding
resulted in the highest joint verification accuracy, and is therefore selected
in this work.

Next, in the training phase, the embedding network accepts an image (resized
to $48\times 48$) as input. Different image resolutions were used $\\{28\times
28,48\times 48,96\times 96\\}$, but $48\times 48$ provided optimal trade-off
between accuracy and training time. The embeddings are trained in (i)
classical, (ii) siamese and (iii) triplet modes. Then, in the testing phase,
we computed the embeddings from the test set. We evaluate the test embeddings
in joint identification and joint verification scenarios.

Although deep learning-based sensor identification methods exist in the
literature [23, 24, 25], we used Enhanced PRNU [26] (with enhancement Model
III) as the sensor identification baseline for all three modalities due to its
low computational burden and effectiveness against multi-spectral images [27].
Enhanced PRNU requires creation of sensor reference patterns, that serve as
gallery and test (probe) noise residuals, that are correlated with the
reference patterns. We used training images to compute the sensor reference
patterns and test images for correlation. A test image is assigned to the
sensor class resulting in the highest correlation value. See [10] for more
details. Test noise residuals computed from JPEG images can be matched
successfully against sensor reference patterns computed from RAW images [28],
thereby, justifying the use of PRNU as a state-of-the-art sensor
identification baseline. We used COTS matcher as the biometric recognition
baseline for iris and face modalities. For the periocular modality, we used a
pretrained ResNet-101 architecture [29] and used the features from layer 170
as the biometric representation for the test samples. This particular
architecture is used because it has demonstrated good performance in biometric
verification on the MICHE-I dataset [10]. The gallery comprises the training
images and the probes are the test images. Since, PRNU can only be used for
the task of sensor identification, we selected to implement both the baselines
only in identification scenario.

We further conducted an experiment using a validation set comprising 50 images
from the MICHE-I dataset (excluded from the test set) to analyze the effect of
the dimensionality of the embedding on the verification performance. To this
end, we used four values of $k=\\{4,8,16,32\\}$, and then selected that value
which results in the highest performance for the remaining experiments.

Figure 4: 2-D projection of the embeddings using t-SNE used for sensor
identification in the OULU-NPU dataset. Each sensor class is sufficiently
discriminated from the rest of the sensors.

##  V Results and Analysis

###  V-A Selection of the metric and dimensionality of embedding

In terms of the choice of the distance/similarity metric, we observed that
standardized euclidean distance metric resulted in better performance compared
to the cosine similarity metric. This can be attributed to the standardization
process which takes into account the intra-class and inter-class variations in
the embeddings. In terms of the choice of the dimensionality of the embedding,
we observed that 8 was the optimal value, since, it resulted in the best
performance (64% on the MICHE-I validation set) as indicated in Figure 3.
Therefore, we used 8-dimensional embedding and standardized Euclidean distance
metric for all the experiments. Furthermore, we presented the t-SNE [30]
visualization of the performance of the embedding network in terms of sensor
identification for the OULU-NPU dataset in Figure 4. The well-separable
clusters corresponding to the six sensors demonstrate the capability of the
embedding network used in this work.

###  V-B Performance of each of the three training modes

In terms of training algorithms, the overall results in both joint
identification and joint verification scenarios indicate that the embedding
network trained in siamese mode outperformed the remaining training paradigms
(see Figures 5 and 6). The reason for the superior performance of siamese
network can be attributed to the use of contrastive loss. Out of the two
contrastive losses, single margin contrastive loss outperformed double margin
contrastive loss. The contrastive loss considers a pair of embeddings at a
time, and tries to either minimize the distance between them if they belong to
the same class, or increases the distance between them by some margin if they
belong to different classes. On the other hand, triplet loss tries to
simultaneously minimize the distance between the anchor and positive sample,
whereas, maximize the distance between the anchor and negative samples. In
this work, the number of negative classes is very high (in a 330 class
dataset, 1 class is positive and the remaining 329 classes are negative). This
makes the task of triplet loss much more complex as compared to contrastive
loss. Given the huge variation in the possible combination of negative
triplets (see Figure 2), we suspect that the triplet loss struggled to
determine the accurate decision boundary between the positive and negative
classes, resulting in an overall reduction in performance. We investigated
different types of triplet mining strategies, and observed that online triplet
mining outperformed offline triplet mining and multi-class N-pair in a
majority of the cases.

TABLE III: Results in the joint identification scenario. Results are reported
in terms of Rank 1 identification accuracies (%). A correct joint
identification implies that both sensor and subject resulted in a match.
Mismatch of either subject or sensor or both will result in an incorrect joint
identification.

Dataset |  Method for baseline | Baseline performance |  Proposed method (%)
---|---|---|---
| Sensor
---
identification
(%)
| Biometric
---
identification
(%)
| Joint
---
identification
(%)
| CASIA-
---
Iris V2
| PRNU
---
COTS
100.00 | 56.52 | 56.52 | 89.67
MICHE-I |  | PRNU
---
ResNet-101
99.86 | 18.05 | 18.05 | 47.53
OULU-NPU |  | PRNU
---
COTS
98.48 | 84.24 | 83.13 | 99.81

###  V-C Results of the joint identification experiment

In terms of the performance in joint identification scenario, Table III
compares the results with the baseline performance for all the datasets. We
reported the baselines for sensor identification (PRNU), biometric
identification (COTS or ResNet), followed by joint identification, separately.
We reiterate that joint identification involves a correct match only if both
sensor and subject labels are correct to allow fair comparison with the
proposed method. Results indicate that the proposed method outperformed the
baseline (joint identification) by 26.41% averaged across all three datasets
computed at Rank 1. The poor performance for the MICHE-I dataset can be
attributed to two factors - firstly, the large number of classes (= 375)
compared to rest of the datasets (see Table I), and secondly, the diverse
acquisition settings (indoor vs. outdoor) resulting in degraded biometric
recognition, and subsequently leading to overall poor performance.
Surprisingly, the proposed method can still outperform the baseline by
$\sim$30%. We have further analyzed this performance in Section V-E. CMC
curves indicate the superior performance of the siamese network in contrast to
classical and triplet networks.

(a)  (b)  (c)

Figure 5: Cumulative Matching Characteristics (CMC) curves for the proposed
method in the joint identification scenario for the following datasets: (a)
CASIA-Iris V2 (b) MICHE-I and (c) OULU-NPU. Refer to Table II for the
different training networks and loss functions indicated in the legend in an
identical order.

(a)  (b)  (c)

Figure 6: Receiver Operating Characteristics (ROC) curves for the proposed
method in the joint verification scenario for the following datasets: (a)
CASIA-Iris V2 (b) MICHE-I and (c) OULU-NPU. Refer to Table II for the
different training networks and loss functions indicated in the legend in an
identical order.  TABLE IV: Results in the joint verification scenario.
Results are reported in terms of true match rate (TMR) at false match rates
(FMRs) of 1% and 5%.

Dataset | TMR@FMR=1% | TMR@FMR=5%
---|---|---
CASIA-Iris V2 | 90.00 | 98.00
MICHE-I | 62.00 | 90.00
OULU-NPU | 100.00 | 100.00

###  V-D Results of the joint verification experiment

In terms of the performance in joint verification scenario, Table IV reports
the results. Results indicate that the proposed method achieved an average
joint TMR of 84% @1% FMR, and an average TMR of 96% @5% FMR, indicating the
strong representative capability of the joint representation. ROC curves in
Figure 6 indicate that the joint representation learnt using siamese network
trained with single margin contrastive loss (see the curve marked Siamese-
SMCL-Emb[Joint]) outperformed the remaining joint representations. We would
like to point out that in [3], the authors achieved 23% (by using feature
level fusion) and 86% (by using score level fusion) at 5% FMR on the MICHE-I
dataset (the authors excluded the Samsung Galaxy Tab 2 subset of the MICHE-I
dataset, which we included in our evaluations). Although their objectives were
different compared to the proposed work (they adopted a fusion rule for
integrating their proposed biometric and sensor recognition performances), we
would like to indicate that the task of joint recognition is difficult. In
spite of that, the proposed method performed reasonably well.

###  V-E Analysis of the performance of the proposed method on MICHE-I dataset

In both cases of joint identification and joint verification experiments, we
observed that the performance of the proposed method evaluated on the MICHE-I
dataset was relatively worse compared to the remaining two datasets. We
hypothesize that the poor performance can be attributed to two reasons: (i)
the image characteristics, and (ii) the variation in the performance across
different lateralities, i.e., left vs. right periocular images. MICHE-I
dataset was assembled as a part of an iris challenge evaluation and contains
images acquired in unconstrained settings (indoor and outdoor settings) having
occlusions (specular reflection and downward gaze). See some challenging
images from the MICHE-I dataset images in Figure 7. In contrast, CASIA and
OULU datasets contain images acquired in controlled settings.

(a)

(b)

(c)

(d)

Figure 7: Example images from the challenging MICHE-I dataset. (a) Occlusion,
(b) Downward gaze and specular reflection, (c) Prominent background in the
outdoor setting and (d) A single image containing both eyes but labeled as
right eye image (061_GT2_OU_F_RI_01_3 where, RI indicates right eye).

We presented the CMC curves corresponding to joint identification results for
two lateralities separately in Figure 8. Results indicate that the proposed
method performed better on left periocular images compared to right periocular
images. This variation in the performance across the two lateralities resulted
in the overall poor performance on the entire MICHE-I dataset. MICHE-I dataset
has an imbalanced distribution of lateralities. Only 30% of the total number
of images are of left periocular images. We hypothesize that the imbalanced
distribution coupled with some mislabeled test case (see Figure 7(d)) may have
further compounded the challenges, resulting in overall poor performance.

(a)  (b)

Figure 8: Cumulative Matching Characteristics (CMC) curves for the proposed
method in the joint identification scenario for the MICHE-I dataset evaluated
separately on the two lateralities. (a) Left periocular images and (b) Right
periocular images. Results indicate that the proposed method performs better
on left periocular images compared to right periocular images.

The main findings from the experiments are as follows:

  1. 1.

The joint biometric and sensor representation performed well in both joint
identification scenario, with an average identification accuracy of $\sim
80\%$ computed at Rank 1, and an average joint verification accuracy of $96\%$
at a false match rate of 5%, averaged across the three biometric modalities.

  2. 2.

The representation is robust across three modalities (iris, face and
periocular), and different sensors (near-infrared iris sensors and visible
smartphone sensors).

  3. 3.

The joint embedding outperformed baselines that used state-of-the-art
commercial biometric matchers and sensor identification schemes across three
datasets corresponding to three biometric modalities and multi-spectral
sensors.

##  VI Summary and Future Work

In this paper, we proposed a one-shot method to simultaneously authenticate
the user and the device from a single image, say a face or an iris image. To
accomplish this task, we developed a method to learn a joint representation
that can be used for combined biometric and sensor (device) recognition. The
joint representation will be useful in remote application scenarios that
employ multiple factor authentication. Examples include mobile banking on
smartphones and accessing secure services. Additionally, the joint
representation will be implicitly privacy-preserving as the biometric and
sensor representations cannot be trivially separated. We used an embedding
network to generate the joint representation. We evaluated the proposed
approach by (i) exploring different training schemes, and (ii) testing the
method on multiple datasets belonging to three different biometric modalities
(iris, face and periocular) in both the identification and verification
scenarios. We observed best performing results of identification accuracy of
99.81% at Rank1 and a verification accuracy of TMR=100% at 1% FMR using the
proposed method.

Future work will involve improving the current framework by incorporating
domain knowledge (sensor characteristics) to further increase the recognition
performances. We will also evaluate the proposed method on datasets with
larger number of subjects and sensors.

## References

  * [1] N. Bartlow, N. Kalka, B. Cukic, and A. Ross, “Identifying sensors from fingerprint images,” in _IEEE Computer Society Conference on Computer Vision and Pattern Recognition Workshops_ , 2009, pp. 78–84.
  * [2] S. Prabhakar, A. Ivanisov, and A. Jain, “Biometric recognition: Sensor characteristics and image quality,” _Instrumentation & Measurement Magazine, IEEE_, vol. 14, pp. 10 – 16, 07 2011.
  * [3] C. Galdi, M. Nappi, and J. L. Dugelay, “Multimodal authentication on smartphones: Combining iris and sensor recognition for a double check of user identity,” _Pattern Recognition Letters_ , vol. 3, pp. 34–40, 2015\.
  * [4] R. Arjona, M. A. Prada-Delgado, I. Baturone, and A. Ross, “Securing minutia cylinder codes for fingerprints through physically unclonable functions: An exploratory study,” in _Proc. of 11th IAPR International Conference on Biometrics (ICB)_ , Gold Coast, Australia, June 2018.
  * [5] A. K. Jain, A. Ross, and K. Nandakumar, “Introduction to biometrics,” _Springer_ , 2011.
  * [6] J. Lukas, J. Fridrich, and M. Goljan, “Digital camera identification from sensor pattern noise,” _IEEE Transactions on Information Forensics and Security_ , vol. 1, no. 2, pp. 205–214, June 2006.
  * [7] C. Galdi, M. Nappi, and J.-L. Dugelay, “Combining hardwaremetry and biometry for human authentication via smartphones,” in _Image Analysis and Processing (ICIAP)_ , V. Murino and E. Puppo, Eds. Springer International Publishing, 2015, pp. 406–416.
  * [8] ——, “Secure user authentication on smartphones via sensor and face recognition on short video clips,” in _Green, Pervasive, and Cloud Computing_ , M. H. A. Au, A. Castiglione, K.-K. R. Choo, F. Palmieri, and K.-C. Li, Eds. Springer International Publishing, 2017, pp. 15–22.
  * [9] M. Chen, J. Fridrich, M. Goljan, and J. Lukas, “Determining image origin and integrity using sensor noise,” _IEEE Transactions on Information Forensics and Security_ , vol. 3, no. 1, pp. 74–90, March 2008.
  * [10] S. Banerjee and A. Ross, “Smartphone camera de-identification while preserving biometric utility,” in _Proc. of 10th IEEE International Conference on Biometrics: Theory, Applications and Systems (BTAS)_ , Tampa, USA, September 2019.
  * [11] A. Uhl and Y. Höller, “Iris-sensor authentication using camera PRNU fingerprints,” in _5th IAPR International Conference on Biometrics (ICB)_ , March 2012, pp. 230–237.
  * [12] L. Wang, Y. Li, and S. Lazebnik, “Learning deep structure-preserving image-text embeddings,” _IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pp. 5005–5013, 2016.
  * [13] “CASIA Iris Database Version 2,” http://biometrics.idealtest.org/dbDetailForUser.do?id=2, [Online accessed: 18th December 2019].
  * [14] M. D. Marsico, M. Nappi, F. Narducci, and H. Proença, “Insights into the results of MICHE I - mobile iris challenge evaluation,” _Pattern Recognition_ , vol. 74, pp. 286 – 304, 2018.
  * [15] Z. Boulkenafet, J. Komulainen, L. Li, X. Feng, and A. Hadid, “OULU-NPU: A mobile face presentation attack database with real-world variations,” _IEEE International Conference on Automatic Face and Gesture Recognition_ , 2017.
  * [16] J. Bromley, I. Guyon, Y. LeCun, E. Säckinger, and R. Shah, “Signature verification using a ”siamese” time delay neural network,” in _Proc. of 6th International Conference on Neural Information Processing Systems_ , 1993, pp. 737–744.
  * [17] S. Chopra, R. Hadsell, and Y. Lecun, “Learning a similarity metric discriminatively, with application to face verification,” _IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , vol. 1, pp. 539– 546, 07 2005.
  * [18] E. Hoffar and N. Ailon, “Deep metric learning using triplet network,” in _International Workshop on Similarity-Based Pattern Recognition_. Springer, 2015, pp. 84–92.
  * [19] M. Schultz and T. Joachims, “Learning a distance metric from relative comparisons,” in _Proceedings of the 16th International Conference on Neural Information Processing Systems_ , S. Thrun, L. K. Saul, and B. Schölkopf, Eds., 2003, pp. 41–48.
  * [20] F. Schroff, D. Kalenichenko, and J. Philbin, “FaceNet: A unified embedding for face recognition and clustering,” in _IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2015, pp. 815–823.
  * [21] K. Sohn, “Improved Deep Metric Learning with Multi-class N-pair Loss Objective,” in _Proceedings of the 30th International Conference on Neural Information Processing Systems_ , D. D. Lee, M. Sugiyama, U. V. Luxburg, I. Guyon, and R. Garnett, Eds., 2016, pp. 1857–1865.
  * [22] “Embedding network tutorial,” https://github.com/adambielski/siamese-triplet/blob/master/Experiments_MNIST.ipynb, [Online accessed: 12th January, 2020].
  * [23] F. Marra, G. Poggi, C. Sansone, and L. Verdoliva, “A deep learning approach for iris sensor model identification,” _Pattern Recognition Letters_ , vol. 113, pp. 46 – 53, 2018.
  * [24] D. Freire-Obregón, F. Narducci, S. Barra, and M. Castrillón-Santana, “Deep learning for source camera identification on mobile devices,” _Pattern Recognition Letters_ , vol. 126, pp. 86 – 91, 2019.
  * [25] A. Agarwal, R. Keshari, M. Wadhwa, M. Vijh, C. Parmar, R. Singh, and M. Vatsa, “Iris sensor identification in multi-camera environment,” _Information Fusion_ , vol. 45, pp. 333 – 345, 2019.
  * [26] C. T. Li, “Source camera identification using enhanced sensor pattern noise,” _IEEE Transactions on Information Forensics and Security_ , vol. 5, no. 2, pp. 280–287, June 2010.
  * [27] S. Banerjee and A. Ross, “Impact of photometric transformations on PRNU estimation schemes: A case study using near infrared ocular images,” in _International Workshop on Biometrics and Forensics (IWBF)_ , June 2018, pp. 1–8.
  * [28] D. Valsesia, G. Coluccia, T. Bianchi, and E. Magli, “User authentication via PRNU-based physical unclonable functions,” _IEEE Transactions on Information Forensics and Security_ , vol. 12, no. 8, pp. 1941–1956, 2017.
  * [29] K. Hernandez-Diaz, F. Alonso-Fernandez, and J. Bigun, “Periocular recognition using CNN features off-the-shelf,” in _Proc. of 17th International Conference of the Biometrics Special Interest Group (BIOSIG)_ , September 2018, pp. 1–5.
  * [30] L. van der Maaten and G. E. Hinton, “Visualizing data using t-SNE,” _Journal of Machine Learning Research_ , vol. 9, pp. 2431–2556, November 2008\.


# Confluent Vandermonde with Arnoldi

Qiang Niu Hui Zhang∗ Youzhou Zhou
School of Mathematics and Physics, Xi’an Jiaotong-Liverpool University, Ren’ai
Rd. 111, Suzhou 215123, China
∗ Corresponding author. Emails: {qiang.niu, hui.zhang,
<EMAIL_ADDRESS>
###### Abstract

In this note, we extend the Vandermonde with Arnoldi method recently advocated
by P. D. Brubeck, Y. Nakatsukasa and L. N. Trefethen [SIAM Review, 63 (2021)
405-415] to dealing with the confluent Vandermonde matrix. To apply the
Arnoldi process, it is critical to find a Krylov subspace which generates the
column space of the confluent Vandermonde matrix. A theorem is established for
such Krylov subspaces for any order derivatives. This enables us to compute
the derivatives of high degree polynomials to high precision. It also makes
many applications involving derivatives possible, as illustrated by numerical
examples. We note that one of the approaches orthogonalizes only the function
values and is equivalent to the formula given by P. D. Brubeck and L. N.
Trefethen [SIAM J. Sci. Comput., 44 (2022) A1205–A1226]. The other approach
orthogonalizes the Hermite data. About which approach is preferable to
another, we made the comparison, and the result is problem dependent.

Keywords Hermite interpolation, Vandermonde, confluent, Krylov subspace,
Arnoldi

MSC 2020 65D15, 65F15, 65N25, 65N35

##  1 Introduction.

The recent work [1] presents a linear algebra method for the polynomial
interpolation and fitting on arbitrary grids. It does not need the prior
knowledge of a well-conditioned polynomial basis, and works easily on two
intervals and subsets of complex plane, for example. The method starts with
the ill-conditioned Vandermonde system

$\sum_{k=0}^{n}x_{j}^{k}c_{k}=f(x_{j}),\quad j=1,\ldots,m$

such that $f(x)\approx\sum_{k=0}^{n}x^{k}c_{k}$, and transforms it into a much
more well-conditioned system

$\sum_{k=0}^{n}q_{jk}d_{k}=f(x_{j}),\quad j=1,\ldots,m$

such that $f(s_{j})\approx\sum_{k=0}^{n}w_{jk}d_{k}$. In matrix form, the
Vandermonde system $V\mathbf{c}=\mathbf{f}$, i.e.,

$\begin{pmatrix}1&x_{1}&x_{1}^{2}&\cdots&x_{1}^{n}\\\
1&x_{2}&x_{2}^{2}&\cdots&x_{2}^{n}\\\ \vdots&\vdots&\vdots&\cdots&\vdots\\\
1&x_{m}&x_{m}^{2}&\cdots&x_{m}^{n}\end{pmatrix}\begin{pmatrix}c_{0}\\\
c_{1}\\\ \vdots\\\ c_{n}\end{pmatrix}=\begin{pmatrix}f(x_{1})\\\ f(x_{2})\\\
\vdots\\\ f(x_{m})\end{pmatrix}$

is transformed into $Q\mathbf{d}=\mathbf{f}$ where $Q=(q_{jk})$ is a matrix
and $\mathbf{d}=(d_{k})$ is a column vector, and the approximation at $s_{j}$,
$j=1,\ldots,M$ is evaluated as $(f(s_{j}))\approx W\mathbf{d}$ where
$W=(w_{jk})$. The columns of $Q$ form an orthogonal basis of the column space
of $V$. To construct $Q$, it is noted that $\mathrm{Col}\,V$ is just the
Krylov subspace
$\mathrm{span}\\{\mathbf{e},X\mathbf{e},X^{2}\mathbf{e},\ldots,X^{n}\mathbf{e}\\}$
with $X=\mathrm{diag}(x_{1},x_{2},\ldots,x_{m})$,
$\mathbf{e}=(1,\ldots,1)^{T}$. So, the Arnoldi process for the Krylov subspace
is exploited to generate the columns of $Q$ as well as a Hessenberg matrix $H$
such that $XQ_{-}=QH$ (here $Q_{-}$ is obtained by deleting the last column of
$Q$). The evaluation matrix $W$ has the first column all one’s and is
constructed column by column such that $SW_{-}=WH$, where
$S=\mathrm{diag}(s_{1},\ldots,s_{M})$ and $W_{-}$ is $W$ without the
$(n+1)$-th column.

Now, we would like to employ the same idea for the Hermite interpolation and
fitting. We start with the confluent Vandermonde system

$\sum_{k=0}^{n}x_{j}^{k}c_{k}=f(x_{j}),\quad\sum_{k=0}^{n}kx_{j}^{k-1}c_{k}=f^{\prime}(x_{j}),\quad
j=1,\ldots,m$

or in matrix form (with $V$ as before)

$\begin{pmatrix}V\\\
V_{{}^{\prime}}\end{pmatrix}\mathbf{c}=\begin{pmatrix}\mathbf{f_{\phantom{{}^{\prime}}}}\\\
\mathbf{f}_{{}^{\prime}}\end{pmatrix},\quad
V_{{}^{\prime}}=\begin{pmatrix}0&1&2x_{1}&3x_{1}^{2}&\ldots&nx_{1}^{n-1}\\\
0&1&2x_{2}&3x_{2}^{2}&\ldots&nx_{2}^{n-1}\\\
\vdots&\vdots&\vdots&\cdots&\vdots\\\
0&1&2x_{m}&3x_{m}^{2}&\ldots&nx_{m}^{n-1}\end{pmatrix},\quad\mathbf{f}_{{}^{\prime}}=\begin{pmatrix}f^{\prime}(x_{1})\\\
f^{\prime}(x_{2})\\\ \vdots\\\ f^{\prime}(x_{m})\end{pmatrix}.$

We notice that the column space of $[V;V_{{}^{\prime}}]$ (semicolon means as
in MATLAB) is a Krylov subspace:

$\begin{pmatrix}V\\\
V_{{}^{\prime}}\end{pmatrix}=\left(\begin{pmatrix}\mathbf{e}\\\
\mathbf{0}\end{pmatrix},\;\;\begin{pmatrix}X&O\\\
I&X\end{pmatrix}\begin{pmatrix}\mathbf{e}\\\
\mathbf{0}\end{pmatrix},\;\;\begin{pmatrix}X&O\\\
I&X\end{pmatrix}^{2}\begin{pmatrix}\mathbf{e}\\\
\mathbf{0}\end{pmatrix},\;\;\ldots,\begin{pmatrix}X&O\\\
I&X\end{pmatrix}^{n}\begin{pmatrix}\mathbf{e}\\\
\mathbf{0}\end{pmatrix}\right).$

So, the Arnoldi process applies and can be used to generate a matrix $Q$ with
orthogonal columns and a Hessenberg matrix $H$. Then, the evaluation matrix
$W$ can be constructed with

$\text{the first column }\begin{pmatrix}\mathbf{e}\\\
\mathbf{0}\end{pmatrix},\quad\begin{pmatrix}S&O\\\ I&S\end{pmatrix}W_{-}=WH.$

Following [1], we give the MATLAB codes for the method of confluent
Vandermonde with Arnoldi.

         function [d,H] = polyfitAh(x,f,fp,n)
         m = size(x,1); Q = zeros(2*m,n+1); Q(1:m,1) = 1; H = zeros(n+1,n);
         for k = 1:n
             q = [x.*Q(1:m,k); Q(1:m,k)+x.*Q(m+1:2*m,k)];
             for j = 1:k
                 H(j,k) = Q(:,j)’*q/m; q = q - H(j,k)*Q(:,j);
             end
             H(k+1,k) = norm(q)/sqrt(m); Q(:,k+1) = q/H(k+1,k);
         end
         d = Q\[f;fp];

         function [y,yp] = polyvalAh(d, H, s)
         M = size(s,1);  n = size(d,1)-1; W = zeros(2*M,n+1); W(1:M,1) = 1;
         for k = 1:n
             w = [s.*W(1:M,k); W(1:M,k)+s.*W(M+1:2*M,k)];
             for j = 1:k
                 w = w - H(j,k)*W(:,j);
             end
             W(:,k+1) = w/H(k+1,k);
         end
         y = W*d;  yp = y(M+1:2*M,:);  y = y(1:M,:);

To compare, we give also the naive method of confluent Vandermonde without
Arnoldi.

         function c = polyfith(x,f,fp,n)
         A = x.^(0:n); o = zeros(size(x)); D = diag(1:n); A = [A; o A(:,1:n)*D]; c = A\[f;fp];

         function [y,yp] = polyvalh(c,s)
         n = length(c)-1;  B = s.^(0:n);  y = B*c;
         o = zeros(size(s));  D = diag(1:n);  yp = [o B(:,1:n)*D]*c;

Other than the Hermite scenario for both fitting and evaluation, we can also
fit only the function values using `polyfitA` (or `polyfit`) from [1] but
evaluates the function values and derivatives using `polyvalAh` (or
`polyvalh`) here. This scenario is more interesting if the data of derivatives
are not available for fitting. In other words, we construct first a discrete
orthogonal basis for the function values by `polyfitA` and then the derivative
matrix `W(M+1:2*M,:)` under the same basis by `polyvalAh`. We notice that an
equivalent formula for doing this has been given in [2], which we should
detail in the Discussion section. There is a third possibility: one can still
construct a discrete orthogonal basis for the Hermite data, but fit only the
function values, and then construct the derivative evaluation matrix under
that basis; however we found no advantage of it. The fourth possibility is to
use the function values orthogonalized basis for Hermite fitting and
evaluation, which we found performs equally well as using the Hermite data
orthogonalized basis for Examples 1–2 in the following sections. We will see
also for the other examples the comparison results of Hermite and function
values orthogonalized bases.

Both the fitting and evaluation can be generalized to high-order derivatives.
For example, let $V_{{}^{\prime\prime}}$ be the second-order derivative matrix
under the standard basis. We have

$\begin{pmatrix}V\\\ V_{{}^{\prime}}\\\
V_{{}^{\prime\prime}}\end{pmatrix}=\left(\begin{pmatrix}\mathbf{e}\\\
\mathbf{0}\\\ \mathbf{0}\end{pmatrix},\;\;\begin{pmatrix}X&O&O\\\ I&X&O\\\
O&2I&X\end{pmatrix}\begin{pmatrix}\mathbf{e}\\\ \mathbf{0}\\\
\mathbf{0}\end{pmatrix},\;\;\begin{pmatrix}X&O&O\\\ I&X&O\\\
O&2I&X\end{pmatrix}^{2}\begin{pmatrix}\mathbf{e}\\\ \mathbf{0}\\\
\mathbf{0}\end{pmatrix},\;\;\ldots,\begin{pmatrix}X&O&O\\\ I&X&O\\\
O&2I&X\end{pmatrix}^{n}\begin{pmatrix}\mathbf{e}\\\ \mathbf{0}\\\
\mathbf{0}\end{pmatrix}\right).$

It may be worthwhile to state the general case as follows.

######  Theorem 1.

Let $V$ be the Vandermonde matrix and $V_{(\ell)}$ be the $\ell$th order
derivative matrix corresponding to $V$. Then, the matrix
$[V;V_{(1)};V_{(2)};\ldots;V_{(\ell)}]$ has the $(k+1)$th column exactly equal
to

$\begin{pmatrix}X&O&\cdots&\cdots&O\\\ I&X&\ddots&\vdots&\vdots\\\
O&2I&\ddots&\ddots&\vdots\\\ \vdots&\ddots&\ddots&\ddots&O\\\ O&\cdots&O&\ell
I&X\end{pmatrix}^{k}\begin{pmatrix}\mathbf{e}\vphantom{X}\\\
\mathbf{0}\vphantom{\ddots}\\\ \vdots\vphantom{X}\\\
\mathbf{0}\vphantom{\ddots}\\\ \mathbf{0}\vphantom{X}\end{pmatrix}.$

##  2 Example 1: interpolation in Chebyshev points.

This example is used in [1] for the Lagrange interpolation. Now we interpolate
the function $f(x)=1/(1+25x^{2})$ and its derivative simultaneously
(`polyfitAh`+`polyvalAh`) by a degree $n$ (odd) polynomial in $m=(n+1)/2$
Chebyshev points $x_{j}=\cos((m-j)\pi/(m-1))$, $1\leq j\leq m$. We also do the
Lagrange interpolation in $n+1$ points as in [1] and evaluate the derivative
as proposed here (`polyfitA`+`polyvalAh`). The results are shown in Figure 1.
We see that with Arnoldi both $f$ and $f^{\prime}$ can be approximated to high
accuracy with the exponential convergence rate, while without Arnoldi the
method is prone to round-off errors as $n>50$. The derivative evaluation
enabled by the confluent Vandermonde with Arnoldi (`polyvalAh`) works equally
well for the Hermite and Lagrange interpolation. Throughout the paper, the
function norm is the maximum norm computed on a very fine grid.

Figure 1: Hermite (solid lines/points) and Lagrange (dashed lines/circles)
interpolation of $f(x)=1/(1+25x^{2})$ on Chebyshev grids of $[-1,1]$ by the
confluent Vandermonde without (left column) and with (right column) Arnoldi.

##  3 Example 2: least-squares fitting on two intervals.

For the function $f(x)=\sqrt{|x|}$ on $[-1,-1/3]\cup[1/5,1]$, we sample the
function values (and its derivatives for the Hermite version) on every
interval with equally spaced points, and fit by least-squares these data with
a degree $n$ polynomials. The number of sampling points $m$ on each interval
is chosen as $m=5(n+1)$ for the Hermite version and $m=10(n+1)$ for the
Lagrange version. The results are shown in Figure 2.

Figure 2: Least-squares fitting of $f(x)=\sqrt{|x|}$ on $[-1,-1/3]\cup[1/5,1]$
(dashed lines/circles) and simultaneously its derivatives (solid lines/points)
in equally spaced points by the confluent Vandermonde without (left column)
and with (right column) Arnoldi.

##  4 Example 3: Hermite Fourier extension.

This example is used in [1] for approximation of $f(x)=1/(10-9x)$ on $[-1,1]$
by the Fourier series on the larger interval $[-2,2]$. Now we consider both
the function values and its derivatives:

$f(x)\approx\mathrm{Re}\left(\sum_{k=0}^{n}\mathrm{e}^{\mathrm{i}k\pi
x/2}c_{k}\right),\quad-\frac{2}{\pi}f^{\prime}(x)\approx\mathrm{Im}\left(\sum_{k=0}^{n}k\mathrm{e}^{\mathrm{i}k\pi
x/2}c_{k}\right).$

Following [1], this is realized with $z:=\mathrm{e}^{\mathrm{i}\pi x/2}$,
$c_{k}=a_{k}+\mathrm{i}b_{k}$ ($a_{k},b_{k}\in\mathbb{R}$) and

$f(x)\approx\sum_{k=0}^{n}(\mathrm{Re}z^{k}a_{k}-\mathrm{Im}z^{k}b_{k}),\quad-\frac{2}{\pi}f^{\prime}(x)\approx\sum_{k=0}^{n}\left(\mathrm{Im}(kz^{k})a_{k}+\mathrm{Re}(kz^{k})b_{k}\right).$

The last line of `polyfith` is changed to

         m = length(x);  c = [real(A(1:m,:))        imag(A(1:m,2:n+1)); ...
                              imag(A(m+1:2*m,:))   -real(A(m+1:2*m,2:n+1))]\[f;fp];
         c = c(1:n+1) - 1i*[0; c(n+2:2*n+1)];

where `fp` stores the scaled derivatives $-\frac{2}{\pi}f^{\prime}(x_{j})$.
Similar changes should be made also in `polyfitAh` but with $Q$ instead of
$A$. In `polyvalh` and `polyvalAh`, now we take the real part for `y` and
imaginary part for `yp`. Note again that the output `yp` stores the scaled
derivatives $-\frac{2}{\pi}p^{\prime}(s_{j})$. Following [1], 500 Chebyshev
points on $[-1,1]$ are used for the Hermite fitting. The results are shown in
Figure 3. The Lagrange fitting is not present in the figure since the
derivative evaluation for it gave totally wrong results.

Figure 3: Fourier extension of $f(x)=1/(10-9x)$ and simultaneously its
derivatives on Chebyshev grids of $[-1,1]$ by the truncated Fourier series on
$[-2,2]$ using the confluent Vandermonde without (left column) and with (right
column) Arnoldi. Bottom right: Hermite (dots) vs function values (circles)
orthogonalized bases.

##  5 Example 4: Laplace equation and Dirichlet-to-Neumann map.

Given a domain $\Omega\subset\mathbb{R}^{2}$ with boundary $\partial\Omega$
and a function $f:\partial\Omega\to\mathbb{R}$, the Dirichlet-to-Neumann map
$T$ involves first solving the Laplace equation $-\Delta u=0$ in $\Omega$,
$u=f$ on $\partial\Omega$ and then evaluating the outer normal derivative
$\partial_{\boldsymbol{\nu}}u=:Tf$ on $\partial\Omega$. By identifying
$(x,y)\in\mathbb{R}^{2}$ and $z=x+\mathrm{i}y\in\mathbb{C}$, $u$ can be viewed
as the real part of an analytic function $h$ in the complex plane, so that
$-\Delta u=0$ is automatically satisfied. Following [1], we approximate $h(z)$
by a polynomial $\sum_{k=0}^{n}c_{k}z^{k}$. Let $c_{k}=a_{k}+\mathrm{i}b_{k}$
($a_{k}$, $b_{k}$ $\in\mathbb{R}$). Then,
$u(z)\approx\sum_{k=0}^{n}(\mathrm{Re}z^{k}a_{k}-\mathrm{Im}z^{k}b_{k})$. To
compute $Tf$, we need only to fit the boundary condition $u=f$ on
$\partial\Omega$ and then evaluate
$\partial_{\boldsymbol{\nu}}u=\mathrm{Re}(\nu(z)h^{\prime}(z))$ (since
$h^{\prime}(z)=\partial_{x}u-\mathrm{i}\partial_{y}u$ and
$\nu(z)=\nu_{1}(z)+\mathrm{i}\nu_{2}(z)$ corresponds to the unit outer normal
$\boldsymbol{\nu}=(\nu_{1},\nu_{2})$ on $\partial\Omega$), or more precisely

$\partial_{\boldsymbol{\nu}}u=\sum_{k=1}^{n}(\nu_{1}(z)\mathrm{Re}z^{k-1}k-\nu_{2}(z)\mathrm{Im}z^{k-1}k)a_{k}-\sum_{k=1}^{n}(\nu_{2}(z)\mathrm{Re}z^{k-1}k+\nu_{1}(z)\mathrm{Im}z^{k-1}k)b_{k}.$

A concrete example is demonstrated in Figure 4, where $m=10n$ points with the
argument of $z\in\partial\Omega$ equally spaced are used for the fitting. We
see that the Hermite orthogonalized basis leads to a bit more accurate results
than the function values orthogonalized basis does.

Figure 4: Dirichlet-to-Neumann map of the Dirichlet data
$f(z)=\mathrm{Re}\,(\log(0.8+z))^{2}$ at
$z\in\partial\Omega=\\{\mathrm{e}^{\mathrm{i}t}(0.7+0.25\cos(4t-2)+0.05\cos(8t-4)):0\leq
t=\arg z\leq 2\pi\\}$ by the confluent Vandermonde without (left column) and
with (right column) Arnoldi. Bottom-right: Hermite (dots) vs function values
(circles) orthogonalized bases.

##  6 Example 5: Steklov eigenvalue problem in an ellipse.

This is a continuation of the previous section. We look for the eigenvalue
$\lambda$ of $T$ and the assoicated eigenfunction
$f:\partial\Omega\to\mathbb{R}$ such that $Tf=\lambda f$. In other words, we
want the real part $u$ of an analytic function in $\Omega$ to satisfy
$\partial_{\boldsymbol{\nu}}u=\lambda u$ on $\partial\Omega$. Let $Q$ be the
matrix from the Arnoldi process on the confluent Vandermonde matrix; see the
Introduction section. We partition the $2m$ rows of $Q$ in half:
$Q=(Q_{0};Q_{1})$ with $Q_{0}$ containing the first $m$ rows that represents
the sampling of basis function values and $Q_{1}$ representing the sampling of
derivatives of basis functions. Let the matrices be separated into real parts
and imaginary parts: $Q_{0}=Q_{0,r}+\mathrm{i}Q_{0,i}$,
$Q_{1}=Q_{1,r}+\mathrm{i}Q_{1,i}$, and $\nu_{1}$, $\nu_{2}$ be diagonal
matrices containing the sampled values of the unit outer normal of
$\partial\Omega$. Then, the discrete Steklov eigenvalue problem is to find
$\mathbf{a}\in\mathbb{R}^{n+1}$, $\mathbf{b}=(0,b_{1},\ldots,b_{n})^{T}$ and
$\lambda\in\mathbb{R}$ such that

$(\nu_{1}Q_{1,r}-\nu_{2}Q_{1,i})\mathbf{a}-(\nu_{1}Q_{1,i}+\nu_{2}Q_{1,r})\mathbf{b}=\lambda(Q_{0,r}\mathbf{a}-Q_{0,i}\mathbf{b}),$

denoted by $Q_{\nu}\boldsymbol{\beta}=\lambda Q_{00}\boldsymbol{\beta}$ with
$\boldsymbol{\beta}=(\mathbf{a};-\hat{\mathbf{b}})$
($\hat{\mathbf{b}}=(b_{1},\ldots,b_{n})^{T}$), which is a rectangular
eigenvalue problem when $m>2n+1$. Following [6], we first do the economic QR
decomposition $Q_{00}=\tilde{Q}_{00}\tilde{R}_{00}$, then multiply the
eigenvalue equation with $\tilde{Q}_{00}^{*}$ to get
$\tilde{Q}_{00}^{*}Q_{\nu}\boldsymbol{\beta}=\lambda\tilde{R}_{00}\boldsymbol{\beta}$,
which is a generalized eigenvalue problem of square matrices and can be solved
by the standard method (using the MATLAB command `eig`). For the solved
example, see Figure 5. The $m=10n+1$ points are sampled at equally spaced
elliptic angle coordinates. The benchmark solution uses $n=400$.

Figure 5: Laplace Steklov eigenvalue problem in the ellipse $z=\cos
t+\mathrm{i}\frac{1}{5}\sin t$ by the confluent Vandermonde without (left
column) and with (right column) Arnoldi. Top: the $k$th eigenfunction $p_{k}$
($\|p_{k}\|=1$), $k=20,40$, for the eigenvalues
$0=\lambda_{1}\leq\ldots\leq\lambda_{k}\leq\ldots$. Bottom: errors of
$p_{20}$, $\lambda_{20}$ ($-$, $\cdot$) and $p_{40}$, $\lambda_{40}$ ($--$,
$\circ$). Bottom right: Hermite vs function values (in black) orthogonalized
bases.

##  7 Example 6: sloshing modes on top of a cup.

This example is taken from [8]. The domain is $\Omega=(0,1)^{2}$. On the left,
right and bottom sides, the homogeneous Neumann condition
$\partial_{\boldsymbol{\nu}}u=0$ is imposed, while on the top side is the
eigenvalue equation $\partial_{\boldsymbol{\nu}}u=\lambda u$. The analytic
solution is $\lambda_{n}=n\pi\tanh(n\pi)$, $u_{n}=\cos(n\pi x)\sinh(n\pi y)$;
see [8]. Let $I:u\mapsto u$ be the identity operator on the top side and
$O:u\mapsto 0u$ be the zero operator on the other sides. Define $B$ on
$\partial\Omega$ piecewise equal to $I$ and $O$. Then, we can write the
boundary equation uniformly as $\partial_{\boldsymbol{\nu}}u=\lambda Bu$ on
$\partial\Omega$. Using the notation from the previous section, the discrete
rectangular eigenvalue problem becomes $Q_{\nu}\boldsymbol{\beta}=\lambda
BQ_{00}\boldsymbol{\beta}$. Now, testing the equation with the column space
basis of $BQ_{00}$ as before is not a good idea because $BQ_{00}$ is rank
deficient. Note that the two sides of the equation live in the sum of the
column spaces of $Q_{\nu}$ and $BQ_{00}$. So, we follow [5] to do the economic
SVD: $(Q_{\nu},BQ_{00})=USV^{*}$, but we keep the first $2n+1$ columns of $U$
as $\hat{U}$ to match the length of $\boldsymbol{\beta}$, then multiply the
eigenvalue equation with $\hat{U}^{*}$ to get
$\hat{U}^{*}Q_{\nu}\boldsymbol{\beta}=\hat{U}^{*}BQ_{00}\boldsymbol{\beta}$.
The results are shown in Figure 6. On each side of the square, we used
$20(n+1)$ Chebyshev points of the first kind for discretization of the
boundary equation. From the bottom right, we can see that the function values
orthogonalized basis leads to instability. Convergence of the eigenfunctions
from the method with Hermite Arnoldi basis slowers down after some $n$, which
could be fixed with more sampling points. For example, for $n=60$ and
$100(n+1)$ points per side, we got the errors $2.4\times 10^{-14}$, $1.2\times
10^{-13}$ for $p_{5}$, $p_{10}$.

Figure 6: Sloshing modes on the top side of a square cup by the confluent
Vandermonde without (left column) and with (right column) Arnoldi. Top: the
$k$th eigenfunction $p_{k}$ ($\|p_{k}\|=1$), $k=5,10$, for the eigenvalues
$0=\lambda_{1}\leq\ldots\leq\lambda_{k}\leq\ldots$. Bottom: errors of the 5th
($-$, $\cdot$) and 10th ($--$, $\circ$) eigenfunctions/eigenvalues. Bottom
right: Hermite vs function values (in black) orthogonalized bases.

##  8 Example 7: indefinite integral on two intervals.

We want to find a polynomial approximation of $\int f(x)\mathrm{d}x\approx
p(x)+C$, $C$ is a general constant. For example, if
$g(x)=\int_{a}^{x}u(t)~{}\mathrm{d}t$, we have $g^{\prime}(x)=f(x)$, $g(a)=0$,
the simplest differential equation with a point value condition. The goal is
to find $p(x)\approx g(x)$. The implementation is mostly the same as
`polyfitAh` but the function header and the last line: `function [d, H] =
polyfitAih(x,f,n)` and `d = Q(m+1:2*m,2:n+1)\f; d = [0; d];` alternatively,
one can orthogonalize only the function values and from the resulting
Hessenberg matrix generate the derivative matrix. The solved example is shown
in Figure 7. From the bottom right, we can see that, for both the Hermite
Arnoldi basis and the function values Arnoldi basis, the error of $p$ can grow
with $n$ in some range but still decreases to $10^{-15}$ afterwards. Such
transient instability does not appear for $f(x)=1/(1+25x^{2})$ on the
Chebyshev grids of $[-1,1]$.

Figure 7: Polynomial indefinite integral of $f(x)=\mathrm{sign}\ x$ on
equidistant grids of $[-1,-1/10]\cup[1/10,1]$ by the confluent Vandermonde
without (left column) and with (right column) Arnoldi. Bottom right: Hermite
(dots) vs function values (circles) orthogonalized bases.

##  9 Discussion.

This work is greatly inspired by the fascinating paper Vandermonde with
Arnoldi [1] by Brubeck, Nakatsukasa and Trefethen. We saw a door open to the
well-conditioned bases of polynomial spaces on arbitrary grids, and to many
interesting applications, e.g., the lightning solvers [2, 3, 9, 10]. In the
spirit of Chebfun [4], it is natural then to consider differential, integral
and other operations of polynomials under the Arnoldi basis. Under the
standard basis, the differential operation or Hermite interpolation leads to
the confluent Vandermonde matrix, which is ill-conditioned [7]. The Arnoldi
process can come to the rescue. However, it is not so straightforward to apply
the method from [1]. The key idea behind the new method is the discovery of a
connection of the confluent Vandermonde matrix to a Krylov subspace. This
enables us to represent and compute the derivatives in a stable way. Two
approaches are proposed to accomplish this task. One is constructing the basis
for the column space of the confluent Vandermonde matrix directly by the
Confluent Vandermonde with Arnoldi method based on the confluent Krylov
subspace that we discovered. Another approach is using the basis generated by
Vandermonde with Arnoldi and constructing the derivative matrix of the basis
by exploiting the confluent Krylov subspace. With the new functionality of
computing derivatives, we are able to do the Hermite interpolation/evaluation
to high precision, to accomodate the Neumann/Robin boundary conditions, to
compute the Dirichlet-to-Neumann map and its eigenmodes, and to do the
indefinite integral.

We realized just upon finalizing this note that a formula of the derivative
matrix has been given in [9, 2]. It is based on the recursive relation of the
polynomial basis from Vandermonde with Arnoldi. Recall that the Arnoldi
process produces the basis matrix $Q$ and the Hessenberg matrix $H$ that
satisfy $XQ_{-}=QH$. If we imagine that the number of rows of $X$ and $Q$ goes
to infinity, then we would get the polynomial recursion
$xp_{k}(x)=\sum_{j=1}^{k}p_{j}(x)h_{jk}+p_{k+1}(x)h_{k+1,k}$ or
$p_{k+1}(x)=(xp_{k}(x)-\sum_{j=1}^{k}p_{j}(x)h_{jk})/h_{k+1,k}$. Taking
derivative gives
$p_{k+1}^{\prime}(x)=(p_{k}(x)+xp_{k}^{\prime}(x)-\sum_{j=1}^{k}p_{j}^{\prime}(x)h_{jk})/h_{k+1,k}$.
It can be seen that our second approach (e.g. `polyfitA` \+ `polyvalAh`) is
equivalent to this formula.

Comparison of the Hermite orthogonalized basis and the function values
orthogonalized basis has been made on each example. It is found that they
perform equally well for the Hermite interpolation/evaluation task. But the
function values orthogonalized basis is unstable for the example of sloshing
modes, while the Hermite orthogonalized basis works stably. We also observed
instability of both the approaches for the Fourier extension example and the
indefinite integral example. A theoretical understanding of the comparison
results and the instability phenomena could be interesting.

## References

  * [1] Pablo D Brubeck, Yuji Nakatsukasa, and Lloyd N Trefethen.  Vandermonde with Arnoldi.  SIAM Review, 63(2):405–415, 2021.
  * [2] Pablo D Brubeck and Lloyd N Trefethen.  Lightning Stokes solver.  SIAM Journal on Scientific Computing, 44(3):A1205–A1226, 2022.
  * [3] Stefano Costa and Lloyd N Trefethen.  AAA-least squares rational approximation and solution of Laplace problems.  arXiv preprint arXiv:2107.01574, 2021.
  * [4] Tobin A Driscoll, Nicholas Hale, and Lloyd N Trefethen.  Chebfun Guide.  Pafnuty Publications, Oxford, 2014.
  * [5] Behnam Hashemi and Yuji Nakatsukasa.  Least-squares spectral methods for ODE eigenvalue problems.  arXiv preprint arXiv:2109.05384, 2021.
  * [6] Behnam Hashemi, Yuji Nakatsukasa, and Lloyd N Trefethen.  Rectangular eigenvalue problems.  arXiv preprint arXiv:2112.13698, 2021.
  * [7] Nicholas J Higham.  Stability analysis of algorithms for solving confluent Vandermonde-like systems.  SIAM Journal on Matrix Analysis and Applications, 11(1):23–41, 1990\.
  * [8] David Mora, Gonzalo Rivera, and Rodolfo Rodríguez.  A virtual element method for the Steklov eigenvalue problem.  Mathematical Models and Methods in Applied Sciences, 25(08):1421–1445, 2015.
  * [9] Yuji Nakatsukasa and Lloyd N Trefethen.  Reciprocal-log approximation and planar PDE solvers.  SIAM Journal on Numerical Analysis, 59(6):2801–2822, 2021.
  * [10] Lloyd N Trefethen.  Numerical conformal mapping with rational functions.  Computational Methods and Function Theory, 20(3):369–387, 2020\.