# 4D Millimeter-Wave Radar in Autonomous Driving: A Survey

Zeyu Han1, Jiahao Wang1, Zikun Xu1, Shuocheng Yang2, Lei He1,∗, Shaobing
Xu1,∗, Jianqiang Wang1,∗ Zeyu Han and Jiahao Wang contribute equally to this
work.This work was supported by the National Natural Science Foundation of
China (NSFC) under grant number 52221005.1School of Vehicle and Mobility,
Tsinghua University, Beijing, China2Xingjian College, Tsinghua University,
Beijing, China∗Correspondence<EMAIL_ADDRESS>(L.H.);
<EMAIL_ADDRESS>(S.X<EMAIL_ADDRESS>(J.W.)

###### Abstract

The 4D millimeter-wave (mmWave) radar, capable of measuring the range,
azimuth, elevation, and velocity of targets, has attracted considerable
interest in the autonomous driving community. This is attributed to its
robustness in extreme environments and outstanding velocity and elevation
measurement capabilities. However, despite the rapid development of research
related to its sensing theory and application, there is a notable lack of
surveys on the topic of 4D mmWave radar. To address this gap and foster future
research in this area, this paper presents a comprehensive survey on the use
of 4D mmWave radar in autonomous driving. Reviews on the theoretical
background and progress of 4D mmWave radars are presented first, including the
signal processing flow, resolution improvement ways, extrinsic calibration
process, and point cloud generation methods. Then it introduces related
datasets and application algorithms in autonomous driving perception and
localization and mapping tasks. Finally, this paper concludes by predicting
future trends in the field of 4D mmWave radar. To the best of our knowledge,
this is the first survey specifically for the 4D mmWave radar.

##  I Introduction

Autonomous driving technology, which aims to provide safe, convenient and
comfortable transportation experiences, is going through rapid development. To
realize high-level autonomous driving, the capabilities of environment
perception, localization, and mapping are crucial. Therefore, the sensors on
autonomous vehicles, such as cameras, LiDARs, and radars, as well as their
algorithms, are attracting increasing research interest.

Among the various sensors, mmWave radars, benefiting from their recognized
advantages of small-size, low-cost, all-weather operation, velocity-measuring
ability, and high range resolution, etc., [1] have always been widely used for
autonomous driving. However, traditional mmWave radars, also known as 3D
mmWave radars, exhibit weak performance in measuring the elevation of targets,
and their data typically only includes range, azimuth, and velocity
information. Additionally, 3D mmWave radars suffer from clutter, noise, and
low resolution, particularly in the angular dimension, which further limits
their applicability to complex perception tasks.

The recent advancement of multiple-input multiple-output (MIMO) antenna
technology has improved elevation resolution, leading to the emergence of 4D
mmWave radar. As the name suggests, 4D mmWave radar can measure four types of
target information: range, azimuth, elevation, and velocity.

The 4D mmWave radar not only serves as an improved version of mmWave radar,
but also introduces numerous significant research topics. The raw data size of
4D mmWave radars is much larger than that of traditional ones, which poses
challenges in signal processing and data generation, not to mention the
clutter and noise. The sparsity and noise of 4D mmWave radar point clouds
generated in the existing signal processing flow are more severe than those of
LiDAR point clouds, necessitating the careful design of perception,
localization and mapping algorithms that account for the 4D mmWave radar’s
inherent characteristic.

Researchers have conducted a number of surveys on the theory and application
of mmWave radars. In recent years, Bilik et al. [2] review the challenges of
mmWave radar in autonomous vehicles and its future trends. Venon et al. [3]
comprehensively summarize the theory and existing perception algorithms of
mmWave radar in autonomous driving, while Harlow et al. [4] focus on mmWave
radar applications in robotics for their survey. It is evident that most
reviews are centered on 3D mmWave radars.

Despite the revolutionary rise of the 4D mmWave radar and correlated
algorithms, there have been few specialized surveys. To address this gap, this
paper presents a thorough review of 4D mmWave radar in autonomous driving. The
main contributions of this work can be summarized as follows:

  * •

To the best of our knowledge, this is the first survey concentrating on 4D
mmWave radar in autonomous driving.

  * •

Given the uniqueness of the 4D mmWave radar, this survey specifically
introduces its theoretical background and signal processing pipeline.

  * •

This paper provides an exhaustive survey of 4D mmWave radar application
algorithms in autonomous driving, covering research on perception,
localization and mapping.

The remainder of this paper is organized as follows: Section II introduces the
basic theory of 4D mmWave radars, including the signal processing flow, data
formats, and the methods for improving resolution. Section III outlines some
extrinsic calibration algorithms. Section IV summarizes some learning-based
methods for generating point clouds. Section V lists available 4D mmWave radar
datasets for researchers’ convenience. Section VI reviews 4D mmWave radar
perception applications, categorized into 4D-radar-only methods and multi-
modal methods. 4D mmWave radar applications in localization and mapping are
presented in Section VII by the sort of odometry, relocalization, and
simultaneous localization and mapping (SLAM). Section VIII discusses future
trends of 4D mmWave radar in autonomous driving, and Section IX draws the
conclusion.

##  II Background of 4D mmWave Radars

For researchers focusing on autonomous driving, fundamental knowledge about 4D
mmWave radars may often be somewhat overlooked. This section briefly reviews
the basic theory and resolution-improving approaches of 4D mmWave radars as
the foundation for the following sections.

###  II-A Signal Processing Flow

Figure 1: The traditional signal processing flow and corresponding data
formats of 4D mmWave radars [5] [6] [7]

The traditional signal processing flow and corresponding data formats of 4D
mmWave radars are shown in Fig.1. In step 1, millimeter waves are transmitted
from transmission (TX) antennas. After reaching surrounding targets, waves are
reflected to reception (RX) antennas. The waveform of most existing 4D mmWave
radars is the Frequency Modulated Continuous Wave (FMCW), which offers
superior resolution compared to other waveforms. In every working cycle (i.e.
chirp) of the transmission antennas of FMCW radars, the frequency of the
signal increases linearly with a starting frequency $f_{c}$, a bandwidth $B$,
a frequency slope $S$, and a time duration $T_{c}$. By measuring the frequency
of the received signal, the range $r$ of the target can be calculated as
follows:

$r=\frac{ct}{2},\quad t=\frac{\Delta f}{S},$ (1)

where $t$ is the time interval between transmission and reception, $c$ is the
light speed, and $\Delta f$ is the frequency difference between the
transmitted and received signals. Meanwhile, one frame of an FMCW radar
contains $N_{c}$ chirps and has a time duration $T_{f}$. To avoid interference
between adjacent chirps, the transmitted and received signals are considered
within one chirp, so the maximum unambiguous detection range of 4D mmWave
radars is restricted by $T_{c}$. Taking an AWR1843 from Texas Instrument as an
instance, its $T_{c}=0.33\mu s$, accordingly its maximum unambiguous range is
50m. Under the assumption that the range of a target in one frame is constant,
the radial relative velocity $v$ of the target can be calculated by the
Doppler effect as follows:

$v=\frac{\Delta fc}{2f_{c}},\quad\Delta f=\frac{\Delta\varphi}{2\pi T_{c}},$
(2)

where the first equation is the Doppler effect formula, $\Delta f$ and
$\Delta\varphi$ are the frequency and phase drift between the received signals
of adjacent two chirps, respectively. It is evident that the range and Doppler
resolutions depend on $f_{c},T_{c},N_{c}$, etc., for which readers can refer
[3] for details.

To estimate the direction-of-arrival (DOA) of the target, a multiple-input
multiple-output (MIMO) antenna design is typically applied in mmWave FMCW
radars. The $n$ TX antennas and $m$ RX antennas consist of $n\times m$ virtual
TX-RX pairs. To separate transmit signals, different TX antennas should
transmit orthogonal signals. By comparing the phase drift between different
TX-RX pairs, distance differences between different pairs to the same target
can be calculated. Further considering the positional relationship along TX
and RX antennas, the DOA of the target can be obtained.

The signals of each pair are mixed by a mixer at step 2 and converted by an
Analog-to-Digital Converter (ADC) at step 3 to obtain raw ADC data. It should
be noted that the coordinates of the pair matrix in Fig. 1 represent sample
timestamps within a chirp and a frame, respectively, and the value in each
cell corresponds to reflection intensity. Then, at step 4, a 2D Fast Fourier
Transformation (FFT) is conducted along the range and Doppler dimensions to
generate the Range-Doppler map, the coordinates of which are range and
velocity. Subsequently, there are two mainstream signal processing flows. The
former one is to first conduct FFT along different TX-RX pairs to induce
azimuth and elevation information (step 5a). acquiring a 4D range-azimuth-
elevation-Doppler tensor, while for 3D mmWave radars, the result is a 3D
range-azimuth-Doppler tensor. At step 6a, the constant false alarm rate (CFAR)
algorithm [8] is usually applied in the four dimensions to filter the tensor
by the intensity of every cell and obtain real targets in the format of point
cloud for downstream tasks [9]. In contrast, the latter signal processing flow
first filters RD maps to generate target cells using also CFAR algorithm (step
5b), then digital beamforming (DBF) method is employed at step 6b to restore
angular information and generate a radar point cloud [6].

###  II-B Methods to Improve Resolution

As mentioned above, the most crucial ability of 4D mmWave radars is to measure
the elevation of targets, which is equal to improving the elevation
resolution. Specific methods can be divided into hardware and software levels
as follows:

####  II-B1 Hardware

At the hardware level, increasing the number of TX-RX pairs or the aperture of
antennas are two primary ways to improve resolution, including:

  * •

Cascading: simply cascading several standard mmWave radar chips [10] can
increase TX-RX pairs, thus improving angular resolution. For instance, a
12TX-16RX (192 pairs) 4D mmWave radar can be formed by cascading four standard
3TX-4RX (12 pairs) radar chips. It is the most straightforward approach, but
the size and power dissipation are also increased.

  * •

Chip integration: integrating more antennas on a chip is another promising
technique [11]. It has the potential to replace cascading, but the disturbance
between antennas remains an unsolved problem.

  * •

Meta-material: apertures constructed by meta-material can significantly
increase the angular resolution while controlling the size [12], but such
methods have not been mature enough to be widely applied.

####  II-B2 Software

By virtually realizing hardware improvement or optimizing signal processing
algorithms along the processing flow, radar resolution can be improved at the
software level.

  * •

Virtual aperture imaging: inspired by the traditional synthetic aperture radar
(SAR), some researchers try to virtually expand the aperture of antennas
through software design, thus enhancing the angular resolution [13]. This
method has a remarkable effect on angular resolution improvement but usually
needs the help of cascading to reduce noise.

  * •

Super-resolution: super-resolution can be achieved by replacing traditional
methods like FFT in the signal processing flow with innovative algorithms
[14], even learning-based algorithms [15]. However, it also requires deeper
research for practical application.

##  III Extrinsic Calibration

Radar point clouds are relatively sparse, and spectrum data is not
sufficiently intuitive. Due to multi-path effects and clutter interference,
noise is also considerable, posing challenges for calibration. For 4D mmWave
radar, the higher resolution alleviates this issue, but there is still a lack
of sufficiently robust online calibration methods.

Following traditional calibration methods of 3D mmWave radars, retro-
reflectors are commonly used to improve calibration accuracy. By carefully
placing several retro-reflectors in specific locations, analyzing sensing
results of the 4D mmWave radar, and comparing them with LiDAR and camera data,
the extrinsic parameters can be calibrated [16]. Instead of successively
calibrating each sensor pair, Domhof et al. calibrate all sensors together
directly with respect to the mobile robot’s body and achieve a median rotation
error of only 0.02 ∘ [17].

However, the practicability of retro-reflectors in real scenarios is limited.
In recent years, some researchers have designed calibration methods for 4D
mmWave radars that do not require specially placed retro-reflectors. Instead,
radar motion measurement is utilized to conduct calibration for radars [18] or
radar-camera pairs [19]. The convenience of these methods is guaranteed, but
verification in extreme weather conditions still needs to be implemented.

##  IV Learning-based Radar Point Cloud Generation

One primary reason for the sparsity of 4D radar point clouds is the
substantial information loss caused by CFAR. To address this problem, an
increasing number of learning-based methods have been proposed to replace
CFAR, and work directly with RD maps or the 4D tensor to improve the quality
of 4D radar point clouds and the performance of downstream autonomous driving
tasks, such as perception and localization. Related works, as well as
datasets, perception, localization and mapping algorithms of 4D mmWave radars
that will be illustrated in the following, are uniformly shown on a timeline
in Fig. 2.

Figure 2: Timeline of 4D mmWave radar point cloud generation, datasets,
perception, localization and mapping algorithms

Generally, CFAR is an optimal detection algorithm if cells are independent and
identically distributed. However, as the targets in the real world usually
have different shapes and occupy multiple cells, the CFAR-type methods will
lead to masking effects, reducing the resolution of point clouds and suffering
from information loss.

Brodeski et al. [20] first apply CNN to RD maps for detection and localization
of multiple objects, which is called DRD(deep radar detection) net. They
formulate the target detection in RD maps as a segmentation task and adopt a
model structure similar to 2D-U-Net[21]. Facing the lack of well-annotated
datasets, specifically for radar RD maps, they refer to the radar calibration
process, arrange retro-reflectors in an anechoic chamber to collect
corresponding data, and map them back to RD maps as labels. Experiments show
that DRD net can operate in real-time ($\sim$20ms for inference) and
outperform classic methods in detection accuracy.

Figure 3: Overview of the radar signal processing chain with [22]

However, the labeling challenge for RD maps still exists since the data
collected in the anechoic chamber differ from those collected in the real-
world driving scene, which is more challenging with the multi-path
reflections, interference, attenuation, etc. To address this challenge, Cheng
et al. [6][22] use LiDAR point clouds as the supervision and successively
design network structures based on U-Net[23] and GAN[24]. In complex roadway
scenes, the generated 4D radar point clouds by [22] not only contain fewer
clutter points but also provide denser point clouds of real targets compared
to the classical CFAR detectors.

##  V Datasets

Currently, available datasets with 4D mmWave radars are introduced in this
section, which are summarized in Table I.

TABLE I: 4D mmWave radar datasets

Dataset |  Num of Labeled data | Data Formats ∗ | 3D Bbox |  Tracking ID | Odometry
---|---|---|---|---|---
Astyx [25] | 0.5K | RPC | ✓ | ✓ | $\times$
Coloradar [9] | 0K∗∗ | ADC, RT, RPC | $\times$ | $\times$ | ✓
VoD [26] | 8.7K | RPC | ✓ | ✓ | ✓
TJ4DRadSet [16] | 7.8K | RPC | ✓ | ✓ | ✓
K-Radar [27] | 35K | RT, RPC | ✓ | ✓ | ✓

  * •

∗ADC: raw ADC data; RT: radar tensor; RPC: radar point cloud;

  * •

∗∗Coloradar is a dataset for radar odometry, containing 2 hours of data
without bounding box labels.

Public datasets are crucial for 4D mmWave radar-based algorithms as they
provide a platform for developing and comparing different algorithms and
stimulate related research. Therefore, we curated published datasets up to the
time of writing hoping these datasets will facilitate the research of 4D
radar-based algorithms.

Public 4D Radar data are fairly rare, Astyx HiRes 2019 is the first existing
dataset [25]. The data provided for free consist of 500 synchronized frames
(radar, LiDAR, camera) containing about 3,000 precisely annotated 3D object
annotations. It can be seen that the amount of data in this dataset is
relatively small.

ColoRadar is a dataset dedicated to research on localization and mapping,
containing approximately 2 hours of data from radar, LiDAR, and the pose
groundtruth[9]. It provides three levels of processing for radar data: raw ADC
data, 3D range-azimuth-elevation tensors derived by compressing the Doppler
dimension of 4D radar tensors, and radar point clouds. This dataset collects
data in several unique environments, both indoors and outdoors, providing a
diverse range of sensor data.

The VoD dataset is a new multi-sensor automotive dataset for multi-class 3D
object detection, which consists of calibrated and synchronized LiDAR, camera,
and radar data[26]. It contains 8693 frames of data acquired in complex urban
traffic which consists of 123106 3D bounding box annotations of both moving
and static objects and tracking IDs for each annotated object which is useful
for tracking. Similarly, the Tj4DRadSet dataset contains 44 consecutive
sequences with a total of 7757 synchronized frames, well-labeled using 3D
bounding boxes and trajectory IDs [16]. Furthermore, it covers much richer and
more challenging driving scenario clips (e.g. Urban roads, highways,
industrial parks).

To the best of our knowledge, K-Radar is currently the largest large-scale
dataset based on 4D mmWave radar and collects 35k frame conditions (e.g.
sunny, foggy, rainy, snowy) [27]. In addition to 4D radar data, K-Radar
provides high-resolution LiDAR point clouds, surround RGB images from four
stereo cameras, and RTK-GPS and IMU data from the ego vehicle. It is worth
mentioning that K-radar is currently the only dataset that provides 4D radar
tensors. In order to facilitate experiments on various neural network
structures, K-radar also provides a visual program to modularize the neural
network training code.

Although 4D mmWave radar is receiving more and more attention from the
academic community, and more and more datasets have been released, compared
with vision or LiDAR, the amount of data is still not large enough.

##  VI Perception Applications

Currently, the point cloud density of 4D mmWave radar has already reached a
level comparable to that of low-beam LiDAR, and 4D mmWave radar exhibits
superior robustness under low visibility and adverse weather conditions.
Therefore, researchers have been attempting to transfer LiDAR point cloud
processing models to 4D mmWave radar for target detection, scene flow
prediction, and other tasks. Furthermore, as described in Section IV, pre-CFAR
radar data contain richer information, promoting some researchers to work
directly with RD maps or 4D tensors, bypassing point cloud generation tasks.

Related work can be further divided into those relying solely on 4D radar or
on multi-modal sensor fusion.

###  VI-A 4D-Radar-only Methods

Naturally, most related 4D-radar-only methods are derived from LiDAR-based
methods. However, due to the sparsity and noise characteristics of mmWave
radar, specific network designs are still required.

####  VI-A1 3D detection

As for the 3D object detection task, according to the difference between model
architectures, perception methods can be divided into CNN-based and
Transformer-based.

Palffy et al. [26] first apply PointPillars to 4D radar point clouds for 3D
detection of multi-class road users. The performance is improved by temporal
integration and by introducing additional features, such as elevation and
Doppler velocity. However, the result of the proposed method (mAP 47.0) is
still far inferior to the LiDAR detector on 64-beam LiDAR (mAP 62.1).

RPFA-Net [28] achieves progress by introducing a radar pillar features
attention (PFA) module. It utilizes self-attention instead of the simplified
PointNet [29] to extract the global context feature from pillars, which aims
to effectively capture the information over a long distance and improve
heading angle estimation accuracy.

As set operators, attention-based Transformer has inherent superiority in
processing these point sets, which have permutation invariance. Therefore, to
cope with sparse and noisy data of 4D mmWave radar, Tan et al. [30] propose a
3D object-detection framework based on multi-frame point clouds. They obtain
the ego vehicle velocity and compensated velocity information from point
clouds first, then accumulate nearby frames to the last frame.

In addition to directly processing perception tasks on radar point clouds,
some studies have turned their attention to RD maps or 4D tensors, aiming to
utilize more underlying hidden information. K-Radar dataset[27] proposes a 3D
object detection baseline that directly consumes 4D tensors as input and
verifies that the elevation information of 4D tensor is essential for 3D
object detection. The proposed model also demonstrates the robustness of 4D
tensor-based perception, especially under adverse weather conditions.

####  VI-A2 Scene flow estimation

Scene flow estimation aims to calculate a 3D motion vector field that
represents the motion of both static and dynamic elements within the
environment concerning the ego agent. While several research has traditionally
relied on different sensing modalities such as cameras or LiDAR for scene flow
estimation, there are also approaches that utilize 4D mmWave radar data to
accomplish this task.

Representatively, Ding et al. [31] propose a novel approach to 4D radar-based
scene flow estimation via cross-modal learning, motivated by the co-located
sensing redundancy in modern autonomous vehicles. Such redundancy implicitly
provides various forms of supervision cues to the radar scene flow estimation,
which can effectively solve the difficulty in labeling radar point clouds.
Specifically, this work introduces a multitask model architecture for the
cross-modal learning problem. Extensive experiments show the state-of-the-art
performance of this method and demonstrate the effectiveness of cross-modal
supervised learning to infer more accurate 4D mmWave radar scene flow.

###  VI-B Fusion Methods

Considering that 4D mmWave radar can already provide point cloud information,
some scholars have fused it with the camera or LiDAR for target detection,
hoping to improve the accuracy and robustness of the perception model.
Generally, there are three fusion levels for different modalities: data level,
feature level, and decision level. Existing research primarily focuses on
feature-level fusion.

As for 4DRV (4D Radar and Vision) fusion, 4D mmWave radar can provide high-
precision depth and velocity information in a low-cost way, compensating for
the shortcomings of cameras and thereby improving the accuracy of 3D
detection. In recent studies, 4D mmWave radar signals are usually transformed
into 2D image-like features so that they can be practically deployed together
with camera images.

Representatively, Meyer et al. [32] apply a network based on [33] to the
fusion of 4D mmWave radar and camera, which is originally developed for
camera-LiDAR fusion. In order to make up for the data format difference, they
discard the Doppler information and only retain the position information and
the reflection intensity of 4D mmWave radar point clouds. The point clouds of
each frame are used to generate a BEV image and then 3D proposals.
Surprisingly, the precision of the fusion network is higher when radar data is
used instead of LiDAR, and it reaches an average precision of 61% AP on Astyx
dataset[25]. The authors argue that the reason may be the LiDAR sensor has
only 16 beams, but further studies are still required.

A subsequent study is performed by Cui et al.[34] with a newly added self-
supervised model adaptation block[35], which dynamically adapts the fusion of
different modalities according to the object properties. Besides, a front view
map is generated from the radar point clouds together with the BEV image. The
presented method outperforms the former study[32] by up to 9.5% in 3D AP. The
front view map can make better use of the elevation information provided by
the 4D mmWave radar, and it is easier to fuse with the monocular camera
feature, balancing between detection accuracy and computational efficiency.

Despite the above advantages of 4DRV fusion, the visual-based branch may still
struggle to work facing aggressive lighting changes or adverse weather
conditions, which in turn affects the overall performance of the model. Thus,
Wang et al.[36] first explore the advantages of 4DRL(4D Radar and LiDAR)
fusion with InterFusion, an interaction-based fusion framework. They design an
InterRAL module(Interaction of Radar and LiDAR) and update pillars from two
modalities to enhance feature expression. Ablation experiments are carried out
to prove its effectiveness.

Figure 4: Overview of the $M^{2}$ 4DRL fusion model [37]

In a following study, Wang et al. [37] propose $M^{2}$-Fusion network that
integrates an interaction-based multi-modal fusion method termed IMMF and a
center-based multi-scale fusion method termed CMSF. Evaluated using the Astyx
dataset[25], it outperforms mainstream LiDAR-based object detection methods
significantly. As LiDAR can accurately detect objects at close range, 4D
mmWave radar has a far detection range due to its penetrability, so 4DRL
fusion has the potential to be a reliable technical solution with low cost and
high quality.

##  VII Localization and Mapping Applications

In severe environments where satellite positioning information is inaccurate,
or high-definition maps are unavailable, localization and mapping by
perception sensors are necessary. Some relevant research is carried out using
the emerging 4D mmWave radar. As the radar point clouds is much lighter than
the tensor, and studies for LiDAR can be transferred to it with minor
adjustments, there is hardly any research on radar tensors about localization
and mapping.

###  VII-A Odometry

Radar odometry estimation is the core of localization, and is also a key
component of SLAM. There is quite a few related research on 4D mmWave radars.
However, due to the sparsity and noise of radars, odometry is mostly generated
with the help of the Inertial Measurement Unit (IMU).

Doer and Trommer make plenty of contributions to this topic using an Unmanned
Aerial Vehicle (UAV). In [38], they estimate the height of the UAV by the
barometer, then utilize a Random Sample and Consensus (RANSAC) based Least
Squares to estimate ego velocity leveraging the Doppler information of the
radar point clouds. The IMU data are fused finally to construct Extended
Kalman Filter (EKF) based radar inertial odometry. On this basis, they
consider Manhattan world assumptions, which assume planes in the environment
are orthogonal to each other, and achieve comparable accuracy as the state-of-
the-art visual-inertial odometry [39]. This contribution is then extended to
multiple radars and shows satisfactory performance under degraded visual
conditions with very little computational resource requirement[40] [41].
Besides, they also research the fusion of 4D mmWave radars with visual and
thermal information to further enhance the result [42]. The only downside of
these researches is that the EKF-based algorithm may face difficulty coping
with large-scale environments since the odometry drift will increase. And in
most works, Manhattan world assumptions may restrict the applicability in
severe outdoor environments.

The EKF framework is also applied by Michalczyk et al. [43]. Instead of
directly fusing IMU with ego velocity estimated by radar point clouds, they
realize 3D point matching across sparse, noisy radar scans to measure the
displacement between radar scans. The estimation of the 3D trajectory of the
UAV reaches a 3.32% drift at the end of the total traveled distance.

Learning-based radar odometry estimation is also explored. Lu et al. [44]
extract the features of radar point clouds and IMU by CNN and RNN encoders,
respectively. Then a two-stage cross-modal attention layer is designed to fuse
these features, and an RNN is used to model the long-term dynamics. The output
of the whole network is the 6-DOF odometry, which achieves a 0.8m Absolute
Trajectory Error (ATE). The performance demonstrates further upgrades with
assistance from the RGB camera or depth camera.

###  VII-B Relocalization

Relocalization depends on high-accuracy online mapping, and is significant
when using high-definition maps or detecting loop closure in SLAM. Considering
the sparsity and noisiness of 4D mmWave radar point clouds, Cheng et al. [45]
make use of Doppler velocity to remove moving objects, and then enhance the
point clouds by merging multiple consecutive scans. Inspired by the famous
PointNet, A multi‐layer perception (MLP) based network is employed to increase
the dimension of each point from 4 to 1024 for pointwise feature extraction.
By comparing extracted features of the current scan and the global map,
relocalization can be realized.

###  VII-C SLAM

Since the above odometry and relocalization are indispensable for SLAM,
research on 4D mmWave radar SLAM has emerged quite recently. Zhuang et al.
[46] develop a whole-process SLAM for 4D mmWave radar point clouds based on
iterative EKF. To avoid the sparsity caused by RANSAC-like methods, they
conduct ego-velocity estimation by an iterative reweighted least squares. The
weight of each radar point also reflects its dynamics, thus can help remove
moving objects. The unconventional distribution-to-distribution matching
between scan and submap further decreases the influence of sparsity. The
effects are shown in Fig. 5.

Figure 5: the mapping performance of 4D iRIOM [46]

Zhuge[47] applies the pose graph to construct a SLAM system adapted from a
LiDAR SLAM method named hdl_graph_slam [48]. Some impressive experiments are
conducted under the influence of smoke and rain to prove the robustness of 4D
mmWave radar in extreme environments.

##  VIII Future Trends

The 4D mmWave radar has the potential to bring about profound changes to
autonomous vehicles. Nonetheless, it is far from mature at the moment. The
future trends of 4D mmWave radar may rely mainly on the following areas.

####  VIII-1 Point cloud enhancement

As the most commonly used data format, 4D mmWave radar point clouds undergo
evident low quality compared with other sensors. The data quality of point
cloud is severely impacted by the characteristic of radars, such as the multi-
path effect. Furthermore, there is an urgent need to refine the information
loss during the signal processing flow, particularly by replacing CFAR with
precisely designed learning-based methods. Learning-based methods for DOA
estimation, instead of DBF methods, can also be explored for super-resolution
angle estimation.

####  VIII-2 Application algorithms redesign

In addition to improving 4D mmWave radar point clouds, the application
algorithms after signal processing is another non-neglectable issue. Up till
now, a number of application algorithms for 4D mmWave radars are modified from
corresponding LiDAR algorithms. The specialty of 4D mmWave radars, such as
velocity measuring ability and adaptive capability in extreme environments
should be further explored by future research. For perception tasks, multi-
modal fusion is undoubtedly the future development direction. However, it
remains to be explored whether the robustness of 4D radar in extreme weather
conditions will be weakened by the integration of other sensor branches. For
4D mmWave radar localization and mapping, the sensor fusion with LiDARs and
cameras has a great deal of room for discovery.

####  VIII-3 Pre-point cloud data utilizing

As for the unique data formats along the 4D mmWave radar signal processing
flow, such as raw ADC data, RD maps, and 4D tensors, utilizing these data to
perform perception, localization and mapping is an interesting but almost
untouched topic. Learning-based models that take advantage of ample
information from these data while keeping acceptable real-time performance may
be a research hot spot.

####  VIII-4 Dataset enriching

Same as all other data-driven research areas, the datasets of 4D mmWave radars
play a significant role in correlating studies. Existing datasets containing
4D mmWave radars are relatively rare. Data formats and scenario richness are
two main fields waiting for expansion.

##  IX Conclusion

This paper offers a comprehensive overview of 4D mmWave radar in autonomous
driving. It successively covers the signal processing theory, datasets,
extrinsic calibration methods, learning-based radar point cloud generation
algorithms, applications in perception, localization and mapping, as well as
future trends. Research on 4D mmWave radars in the field of autonomous driving
is still in progress and holds great potential for future advancements.

## References

  * [1] M. Jiang, G. Xu, H. Pei, Z. Feng, S. Ma, H. Zhang, and W. Hong, “4d high-resolution imagery of point clouds for automotive mmwave radar,” _IEEE Transactions on Intelligent Transportation Systems_ , pp. 1–15, 2023\.
  * [2] I. Bilik, O. Longman, S. Villeval, and J. Tabrikian, “The rise of radar for autonomous vehicles: Signal processing solutions and future research directions,” _IEEE Signal Processing Magazine_ , vol. 36, no. 5, pp. 20–31, Sep. 2019.
  * [3] A. Venon, Y. Dupuis, P. Vasseur, and P. Merriaux, “Millimeter wave fmcw radars for perception, recognition and localization in automotive applications: A survey,” _IEEE Transactions on Intelligent Vehicles_ , vol. 7, no. 3, pp. 533–555, Sep. 2022.
  * [4] K. Harlow, H. Jang, T. D. Barfoot, A. Kim, and C. Heckman, “A new wave in robotics: Survey on recent mmwave radar applications in robotics,” May 2023.
  * [5] S. Abdulatif, Q. Wei, F. Aziz, B. Kleiner, and U. Schneider, “Micro-doppler based human-robot classification using ensemble and deep learning approaches,” in _2018 IEEE Radar Conference (RadarConf18)_ , Apr. 2018, pp. 1043–1048.
  * [6] Y. Cheng, J. Su, H. Chen, and Y. Liu, “A new automotive radar 4d point clouds detector by using deep learning,” in _ICASSP 2021 - 2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , Jun. 2021, pp. 8398–8402.
  * [7] G. Hakobyan and B. Yang, “High-performance automotive radar: A review of signal processing algorithms and modulation schemes,” _IEEE Signal Processing Magazine_ , vol. 36, no. 5, pp. 32–44, Sep. 2019.
  * [8] H. Rohling, “Radar cfar thresholding in clutter and multiple target situations,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. AES-19, no. 4, pp. 608–621, Jul. 1983.
  * [9] A. Kramer, K. Harlow, C. Williams, and C. Heckman, “Coloradar: The direct 3d millimeter wave radar dataset,” _The International Journal of Robotics Research_ , vol. 41, no. 4, pp. 351–360, Apr. 2022.
  * [10] A. Och, C. Pfeffer, J. Schrattenecker, S. Schuster, and R. Weigel, “A scalable 77 ghz massive mimo fmcw radar by cascading fully-integrated transceivers,” in _2018 Asia-Pacific Microwave Conference (APMC)_ , Nov. 2018, pp. 1235–1237.
  * [11] P. Ritter, M. Geyer, T. Gloekler, X. Gai, T. Schwarzenberger, G. Tretter, Y. Yu, and G. Vogel, “A fully integrated 78 ghz automotive radar system-an-chip in 22nm fd-soi cmos,” in _2020 17th European Radar Conference (EuRAD)_ , Jan. 2021, pp. 57–60.
  * [12] J. Jiang, Y. Li, L. Zhao, and X. Liu, “Wideband mimo directional antenna array with a simple meta-material decoupling structure for x-band applications,” _The Applied Computational Electromagnetics Society Journal (ACES)_ , pp. 556–566, May 2020.
  * [13] Z. Wu, L. Zhang, and H. Liu, “Generalized three-dimensional imaging algorithms for synthetic aperture radar with metamaterial apertures-based antenna,” _IEEE Access_ , vol. 7, pp. 59 716–59 727, 2019.
  * [14] F. Zhang, C. Wu, B. Wang, and K. J. R. Liu, “mmeye: Super-resolution millimeter wave imaging,” _IEEE Internet of Things Journal_ , vol. 8, no. 8, pp. 6995–7008, Apr. 2021.
  * [15] H.-W. Cho, W. Kim, S. Choi, M. Eo, S. Khang, and J. Kim, “Guided generative adversarial network for super resolution of imaging radar,” in _2020 17th European Radar Conference (EuRAD)_. New York: Ieee, Jan. 2021, pp. 144–147.
  * [16] L. Zheng, Z. Ma, X. Zhu, B. Tan, S. Li, K. Long, W. Sun, S. Chen, L. Zhang, M. Wan, L. Huang, and J. Bai, “Tj4dradset: A 4d radar dataset for autonomous driving,” in _2022 IEEE 25th International Conference on Intelligent Transportation Systems (ITSC)_ , Oct. 2022, pp. 493–498.
  * [17] J. Domhof, J. F. P. Kooij, and D. M. Gavrila, “A joint extrinsic calibration tool for radar, camera and lidar,” _IEEE Transactions on Intelligent Vehicles_ , vol. 6, no. 3, pp. 571–582, Sep. 2021.
  * [18] Y. Bao, T. Mahler, A. Pieper, A. Schreiber, and M. Schulze, “Motion based online calibration for 4d imaging radar in autonomous driving applications,” in _2020 German Microwave Conference (GeMiC)_ , Mar. 2020, pp. 108–111.
  * [19] E. Wise, J. Persic, C. Grebe, I. Petrovic, and J. Kelly, “A continuous-time approach for 3d radar-to-camera extrinsic calibration,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , May 2021, pp. 13 164–13 170.
  * [20] D. Brodeski, I. Bilik, and R. Giryes, “Deep radar detector,” in _2019 IEEE Radar Conference (RadarConf)_ , Apr. 2019, pp. 1–6.
  * [21] C. Zhang, W. Luo, and R. Urtasun, “Efficient convolutions for real-time semantic segmentation of 3d point clouds,” in _2018 International Conference on 3D Vision (3DV)_ , Sep. 2018, pp. 399–408.
  * [22] Y. Cheng, J. Su, M. Jiang, and Y. Liu, “A novel radar point cloud generation method for robot environment perception,” _IEEE Transactions on Robotics_ , vol. 38, no. 6, pp. 3754–3773, Dec. 2022.
  * [23] O. Ronneberger, P. Fischer, and T. Brox, “U-net: Convolutional networks for biomedical image segmentation,” in _Lecture Notes in Computer Science_ , ser. Lecture Notes in Computer Science, N. Navab, J. Hornegger, W. M. Wells, and A. F. Frangi, Eds. Cham: Springer International Publishing, 2015, pp. 234–241.
  * [24] P. Luc, C. Couprie, S. Chintala, and J. Verbeek, “Semantic segmentation using adversarial networks,” in _NIPS Workshop on Adversarial Training_ , 2016\.
  * [25] M. Meyer and G. Kuschk, “Automotive radar dataset for deep learning based 3d object detection,” in _2019 16th European Radar Conference (EuRAD)_ , Oct. 2019, pp. 129–132.
  * [26] A. Palffy, E. Pool, S. Baratam, J. F. P. Kooij, and D. M. Gavrila, “Multi-class road user detection with 3+1d radar in the view-of-delft dataset,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 2, pp. 4961–4968, Apr. 2022.
  * [27] D.-H. Paek, S.-H. Kong, and K. T. Wijaya, “K-radar: 4d radar object detection for autonomous driving in various weather conditions,” in _Thirty-Sixth Conference on Neural Information Processing Systems Datasets and Benchmarks Track_ , 2022.
  * [28] B. Xu, X. Zhang, L. Wang, X. Hu, Z. Li, S. Pan, J. Li, and Y. Deng, “Rpfa-net: A 4d radar pillar feature attention network for 3d object detection,” in _2021 IEEE International Intelligent Transportation Systems Conference (ITSC)_ , Sep. 2021, pp. 3061–3066.
  * [29] R. Q. Charles, H. Su, M. Kaichun, and L. J. Guibas, “Pointnet: Deep learning on point sets for 3d classification and segmentation,” in _2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_. Honolulu, HI: IEEE, Jul. 2017, pp. 77–85.
  * [30] B. Tan, Z. Ma, X. Zhu, S. Li, L. Zheng, S. Chen, L. Huang, and J. Bai, “3d object detection for multi-frame 4d automotive millimeter-wave radar point cloud,” _IEEE Sensors Journal_ , pp. 1–1, 2022.
  * [31] F. Ding, A. Palffy, D. M. Gavrila, and C. X. Lu, “Hidden gems: 4d radar scene flow learning using cross-modal supervision,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2023, pp. 9340–9349.
  * [32] M. Meyer and G. Kuschk, “Deep learning based 3d object detection for automotive radar and camera,” in _2019 16th European Radar Conference (EuRAD)_ , Oct. 2019, pp. 133–136.
  * [33] J. Ku, M. Mozifian, J. Lee, A. Harakeh, and S. L. Waslander, “Joint 3d proposal generation and object detection from view aggregation,” in _2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , Oct. 2018, pp. 1–8.
  * [34] H. Cui, J. Wu, J. Zhang, G. Chowdhary, and W. R. Norris, “3d detection and tracking for on-road vehicles with a monovision camera and dual low-cost 4d mmwave radars,” in _2021 IEEE International Intelligent Transportation Systems Conference (ITSC)_ , Sep. 2021, pp. 2931–2937.
  * [35] A. Valada, R. Mohan, and W. Burgard, “Self-supervised model adaptation for multimodal semantic segmentation,” _International Journal of Computer Vision_ , vol. 128, no. 5, pp. 1239–1285, May 2020.
  * [36] L. Wang, X. Zhang, B. Xv, J. Zhang, R. Fu, X. Wang, L. Zhu, H. Ren, P. Lu, J. Li, and H. Liu, “Interfusion: Interaction-based 4d radar and lidar fusion for 3d object detection,” in _2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , Oct. 2022, pp. 12 247–12 253.
  * [37] L. Wang, X. Zhang, J. Li, B. Xv, R. Fu, H. Chen, L. Yang, D. Jin, and L. Zhao, “Multi-modal and multi-scale fusion 3d object detection of 4d radar and lidar for autonomous driving,” _IEEE Transactions on Vehicular Technology_ , vol. 72, no. 5, pp. 5628–5641, May 2023.
  * [38] C. Doer and G. F. Trommer, “An ekf based approach to radar inertial odometry,” in _2020 IEEE International Conference on Multisensor Fusion and Integration for Intelligent Systems (MFI)_ , Sep. 2020, pp. 152–159.
  * [39] ——, “Yaw aided radar inertial odometry using manhattan world assumptions,” in _2021 28th Saint Petersburg International Conference on Integrated Navigation Systems (ICINS)_ , May 2021, pp. 1–9.
  * [40] ——, “X-rio: Radar inertial odometry with multiple radar sensors and yaw aiding,” _Gyroscopy and Navigation_ , vol. 12, no. 4, pp. 329–339, Dec. 2021\.
  * [41] C. Doer, J. Atman, and G. F. Trnmmer, “Gnss aided radar inertial odometry for uas flights in challenging conditions,” in _2022 IEEE Aerospace Conference (AERO)_ , Mar. 2022, pp. 1–10.
  * [42] C. Doer and G. F. Trommer, “Radar visual inertial odometry and radar thermal inertial odometry: Robust navigation even in challenging visual conditions,” in _Gyroscopy and Navigation_ , Sep. 2021, pp. 331–338.
  * [43] J. Michalczyk, R. Jung, and S. Weiss, “Tightly-coupled ekf-based radar-inertial odometry,” in _2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. Kyoto, Japan: IEEE, Oct. 2022, pp. 12 336–12 343.
  * [44] C. X. Lu, M. R. U. Saputra, P. Zhao, Y. Almalioglu, P. P. B. de Gusmao, C. Chen, K. Sun, N. Trigoni, and A. Markham, “milliego: Single-chip mmwave radar aided egomotion estimation via deep sensor fusion,” Oct. 2020.
  * [45] Y. Cheng, C. Pang, M. Jiang, and Y. Liu, “Relocalization based on millimeter wave radar point cloud for visually degraded environments,” _Journal of Field Robotics_ , vol. 40, no. 4, pp. 901–918, Jun. 2023.
  * [46] Y. Zhuang, B. Wang, J. Huai, and M. Li, “4d iriom: 4d imaging radar inertial odometry and mapping,” Mar. 2023.
  * [47] H. Zhuge, “A 4d imaging radar slam system for large-scale environments based on pose graph optimization,” Master’s thesis, Nanyang Technological University, Singapore, 2023.
  * [48] K. Koide, J. Miura, and E. Menegatti, “A portable three-dimensional lidar-based system for long-term and wide-area people behavior measurement,” _International Journal of Advanced Robotic Systems_ , vol. 16, Feb. 2019.


# Optimal design for kernel interpolation: applications to uncertainty
quantification

Akil Narayan  ,  Liang Yan  and  Tao Zhou

###### Abstract.

The paper is concerned with classic kernel interpolation methods, in addition
to approximation methods that are augmented by gradient measurements. To apply
kernel interpolation using radial basis functions (RBFs) in a stable way, we
propose a type of quasi-optimal interpolation points, searching from a large
set of candidate points, using a procedure similar to designing Fekete points
or power function maximizing points that use pivot from a Cholesky
decomposition. The proposed quasi-optimal points results in smaller condition
number, and thus mitigates the instability of the interpolation procedure when
the number of points becomes large. Applications to parametric uncertainty
quantification are presented, and it is shown that the proposed interpolation
method can outperform sparse grid methods in many interesting cases. We also
demonstrate the new procedure can be applied to constructing gradient-enhanced
Gaussian process emulators.

Department of Mathematics and Scientific Computing and Imaging Institute,
University of Utah, Salt Lake City, UT 84112. Email<EMAIL_ADDRESS>A.
Narayan is partially supported by NSF DMS-1848508 and AFOSR FA9550-20-1-0338.

School of Mathematics, Southeast University, Nanjing, China. Email:
<EMAIL_ADDRESS>L. Yan is supported by NSF of China (No.11771081), the
science challenge project (No. TZ2018001) and Zhishan Young Scholar Program of
SEU

LSEC, Institute of Computational Mathematics, Academy of Mathematics and
Systems Science, Chinese Academy of Sciences, Beijing 100190, China. Email:
<EMAIL_ADDRESS>T. Zhou is partially supported by the National Key R$\&$D
Program of China(No. 2020YFA0712000), the NSF of China (under grant numbers
11822111, 11688101and 11731006), the science challenge project (No.
TZ2018001), the Strategic Priority Research Program of Chinese Academy of
Sciences (No. XDA25000404) and youth innovation promotion association (CAS)

##  1\. Introduction

The field of uncertainty quantification (UQ) has received intensive attention
in recent years as theoreticians and practitioners have tackled problems in
the diverse areas of stochastic analysis, high-dimensional approximation, and
Bayesian learning. A fundamental problems in UQ is to approximation of a
multivariate function $u(Z)$, where the parameters
$Z=(Z^{(1)},\cdots,Z^{(d)})$ is a $d$-dimensional random vector. The function
$u$ might be a solution resulting from a stochastic PDE problem or a derived
quantities of interest (QoI) from such a system. The most popular approach is
to expand the function $u$ in a generalized polynomial chaos (gPC) basis [10,
36] and approximate its coefficients with a Galerkin projection. In this case,
a serious drawback is that existing deterministic solvers must be rewritten.

In recent years, more attention has been devoted to the construction of a
surrogate of $u$ based on data
$D_{N}\coloneqq\\{(z_{j},u(z_{j}))\\}_{j=1}^{N}.$ This in general refers to
stochastic collocation methods, as the sample points $\\{z_{j}\\}$ are chosen
from the parametric domain. The procedure is also known as non-intrusive
response construction.“Non-intrusive” effectively means that existing black-
box tools can be used in their current form. Most existing work concerns the
collocation based polynomial approximations, this include the interpolation
technic based on sparse grid [9, 21, 23, 35], the least-squares projection
onto polynomial spaces [3, 11, 22, 13], the compressive sampling method with
$\ell_{1}$ minimization [5, 16, 12, 37]. Another approach to construct the
surrogate using $D_{N}$ is the so called Gaussian process (GP) regression [2,
1, 26, 31]. The GP provides an analytically tractable Bayesian framework where
prior information about the $f$ can be encoded in the covariance function, and
the uncertainty about the prediction is easily quantified given measurements.
One drawback of GP is the number of samples required for an accurate surrogate
increased exponentially as the number of input parameters grows. One method of
mitigating this limitation is the incorporation of gradient information into
the training of the surrogate [4, 18, 20, 32]. By incorporating derivative
values, the cost associated with training an accurate surrogate can be greatly
reduced. Unfortunately, gradient-enhanced GP emulators raise more
computational challenges compared to ordinary GP emulators [14]. Especially,
the condition number of the gradient-enhanced GP emulator will increase much
faster than the original GP emulator.

In this work, we approximate $u(Z)$ using a kernel interpolation through a
limited number of support points. Over the last four decades, kernel methods
using radial basis functions (RBFs) have been successfully applied to
scattered data interpolation/approximation in hight dimensions [7, 34]. It is
known that the kernel interpolation process is equivalent to finding the
unbiased estimator for a Gaussian process with covariance from measurements
[30, 28]. The theoretical connections between kernel methods and Gaussian
process regressions for the classical interpolations has recently been
highlighted [6, 30]. Similar to GP progression, the numerical performance of
kernel interpolation is depended on the number of sample points $N$ and the
so-called shape parameter $\epsilon$. If $\epsilon$ is too large, the
interpolant or solution will usually be inaccurate. On the other hand, if
$\epsilon$ is too small, the condition number of the resultant matrix system
will become so bad that linear solvers may perform poorly. Searching for the
best shape parameter is still an open problem. On the other hand, for some
fixed $\epsilon$, the problem of ill-conditioning can still arise either if
the total number of points (commonly called centers) is too large or if some
points are closely clustered. How to select the optimal centers in the kernel
methods is still a challenging work. In this paper, we consider a more
practical and solvable problem:

Given a large set of candidate points, how can we choose an “optimal” subset
of samples for interpolative kernel approximation?

To this end, we borrow the notion of Fekete points in the polynomial
interpolation community. The selection strategy we propose only involves
linear algebra, and essentially amounts to Cholesky-type decompositions for
the design matrix. More precisely, the strategy is data-independent, so it can
proceed offline before data is collected. We shown that the selected quasi-
optimal points results in a smaller condition number, and thus postpone
dramatically the instability of the interpolation procedure when the number of
points goes to large. Comparisons with other sampling strategies are
presented, and it is noticed that our method demonstrates much improved
stability property. Applications to parametric uncertainty quantification are
also presented, and it is shown that the proposed kernel interpolation method
over perform the sparse grid methods in many interesting cases.

The remained of this paper is organized as follows. Section 2 summaries the
problem formulation, where we start with an abstract problem and provide
detains on the kernel interpolation. Our approach for choosing the optimal
sub-set of samples is introduced in Section 3. Several numerical examples are
considered in Section 4 to study the empirical performance of the new
approach. The conclusions are drawn in the final section.

##  2\. Problem setup and kernel interpolation

Here we consider a general setting for PDEs with random input parameters. Let
$Z=(Z^{(1)},\cdots,Z^{(d)})\in I_{Z}\subseteq\mathbb{R}^{d},d\geq 1$ be a
parameter domain representing the uncertain inputs of the system. Consider

(1)
$\displaystyle\begin{cases}u_{t}(x,t,Z)=\mathcal{L}(u),&\mathcal{D}\times(0,T]\times
I_{Z}\\\ \mathcal{B}(u)=0,&\partial{\mathcal{D}}\times(0,T]\times I_{Z}\\\
u=u_{0}&\mathcal{D}\times\\{t=0\\}\times I_{Z}\end{cases}$

where $\mathcal{D}\in\mathbb{R}^{l},l=1,2,3$, is the physical domain, and
$T>0$ is the terminal time. Here $\mathcal{L}$ is a differential operator and
$\mathcal{B}$ a boundary operator that may involve differential operators with
respect to the spacial variable $x.$ In a probabilistic framework, $Z$ is
modeled as a random variable on a complete probability space
$(\Omega,\mathcal{F},\mathcal{P})$ equipped with the probability distribution
function $F_{Z}(z)=P(Z\leq z)$, with $z\in\mathbb{R}^{d}$. Note that $Z$
reflects the complete parameterization of uncertain inputs and can include
both physical parameters of the system and hyperparameters that characterize
certain input processes.

Throughout the paper, we assume that conditioned on the $i$th independent
sample of $Z$, denoted by $z_{i}$, a numerical solution to this problem may be
identified by a fixed deterministic solver, such as finite element solves,
finite difference solvers, ect. Let $\Xi=\\{z_{1},\cdots,z_{N}\\}\subset
I_{Z},N\geq 1$, be a set of sample points in $I_{Z}$ . For any fixed
$x\in\mathcal{D}$ and $t>0$, we shall denote $u_{j}=u(x,t,z_{j})$ for
simplicity. Hereafter we will suppress the notions of $x$ and $t$ whenever
possible, with the understanding that our statements are made for all fixed
$x$ and $t$. Once the pairings $(z_{j},u_{j}),j=1,\cdots,N$, are obtained, the
objective is to construct a function $u_{N}(Z)$ such that $u_{N}(Z)\approx
u(Z)$ in a proper sense.

###  2.1. Kernel interpolation in parameter space

In the current work, we shall consider the classic kernel interpolations using
RBF to construct $u_{N}$ based on unstructured meshes. The pairing information
will be enforced exactly by requiring $u_{N}(z_{j})=u_{j}$ for all
$j=1,\cdots,N$. In the kernel interpolation framework, we first pick a
translation-invariant radial kernel
$K=\Phi(\epsilon\|\cdot-\cdot\|):\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$.
Here, $\Phi:\mathbb{R}\rightarrow\mathbb{R}$ known as the radial basis
functions (RBF), $\epsilon$ is a shape parameter and $\|\cdot\|$ is, for
example, Euclidean distance. Using the set $\Xi=\\{z_{j}\\}^{N}_{j=1}\subset
I_{Z}$, which is usually referred to as the trial centers in kernel
interpolation, we can define a finite-dimensional trial space in the form of

(2)
$\displaystyle\mathcal{U}_{\Xi}=\mathcal{U}_{\Xi,K}:=\textmd{span}\\{\Phi(\epsilon\|\cdot-
z_{j}\|)|z_{j}\in\Xi\\}.$

In applications, it is trivial that having a suitable trial space is essential
for the performance of kernel methods. A good trial space should contain some
good approximation to the solution. In the present work we shall confine
ourselves to the case of positive definite kernels/RBFs. The important role of
positive definiteness for kernel interpolation was pointed out in [19] and,
for example, guarantees existence and uniqueness of the approximation. Widely
used positive definite kernels/RBFs are the Gaussians, with
$\Phi(r)=\exp(-r^{2})$, the inverse multiquadrics (IMQ) with
$\Phi(r)=1/\sqrt{1+r^{2}}$, and the family of compactly supported (CS) RBFs
[33, 34].

Now, let us consider kernel interpolation problems. Note that any trial
function $u_{N}$ is a linear combination of the basis used in defining (2) and
is in the form of

(3) $\displaystyle u_{N}(Z)=\sum^{N}_{j=1}c_{j}\,\Phi(\epsilon\|Z-z_{j}\|),$

for some coefficients $\mathbf{c}=[c_{1},\cdots,c_{N}]^{T}\in\mathbb{R}^{N}$.
Using the interpolation conditions $u_{N}(z_{i})=u_{i},i=1,\cdots,N$, we
obtain the following linear system

(4)
$\displaystyle\begin{bmatrix}\Phi(\epsilon\|z_{1}-z_{1}\|)&\cdots&\Phi(\epsilon\|z_{1}-z_{N}\|)\\\
\vdots&\ddots&\vdots\\\
\Phi(\epsilon\|z_{N}-z_{1}\|)&\cdots&\Phi(\epsilon\|z_{N}-z_{N}\|)\end{bmatrix}\begin{bmatrix}c_{1}\\\
\vdots\\\ c_{N}\end{bmatrix}=\begin{bmatrix}u_{1}\\\ \vdots\\\
u_{N}\end{bmatrix}\quad\mbox{or}\quad\mathbf{Ac=u},$

where the matrix $\mathbf{A}=K(\Xi,\Xi)$ is symmetric with entries
$\mathbf{A}_{ij}=\Phi(\epsilon\|z_{i}-z_{j}\|)$ for $z_{i},z_{j}\in\Xi$. The
coefficients $\mathbf{c}$ are unique, if the interpolation matrix $\mathbf{A}$
is invertible.

Remark 1: It is known that the kernel interpolation process in trial space (2)
is equivalent to finding the unbiased estimator for a Gaussian process with
covariance $K$ from realization $\mathbf{u}$ at $\Xi$, see A.

Remark 2: If two trial centers $z_{i},z_{j}\in\Xi$ are too closed to each
other, identically shaped trial basis functions centered at these centers have
nearly equal values at $\Xi$. This leads to two columns of nearly identical
values and the problem of ill-conditioning. Yet, we need to decide what the
experimental configuration $\Xi$ is in order to effectively use the kernel
interpolation.

###  2.2. Leave-One-Out Cross Validation (LOOCV)

Noted that the accuracy of the solution of Eq. (4) and the well-conditioning
of the matrix $\mathbf{A}$ depends on the shape parameter $\epsilon$. Much
work addressese this issue, see, e.g. [7, 8, 25, 29, 38], and the references
therein. It is out of our scope of this paper to thoroughly analyze choosing
optimal shape parameters. In this work, we use the leave-one out cross
validation (LOOCV) algorithm [8, 25] to select the value of the shape
parameter. In the RBF interpolation setting the LOOCV algorithm targets
solution of the underlying linear system $\mathbf{Ac}=\mathbf{u}$, where the
entries of the matrix $\mathbf{A}$ depend on the shape parameter $\epsilon$
which we seek to optimize. The cost function defined as the norm of the error
vector $e(\epsilon)$ has entries [8, 25],

$e_{i}(\epsilon)=\frac{c_{i}}{\mathbf{A}^{-1}_{ii}},$

where $c_{i}$ is the $i$th component of the vector $\mathbf{c}$ and
$\mathbf{A}^{-1}_{ii}\coloneqq\left(\mathbf{A}^{-1}\right)_{ii}$ is the $i$th
diagonal element of the inverse of the coefficient matrix. Thus, the optimal
value of the shape parameter is considered as the one which minimizes the cost
function $e(\epsilon)$. We define the optimal shape parameter as the value
$\epsilon^{*}$ such that

(5) $\|e(\epsilon^{*})\|=\min_{\epsilon}\|e(\epsilon)\|.$

###  2.3. A Gradient enhanced approach

We consider inclusion of gradient measurements in the kernel interpolation.
Consider the availability of the following data:
$\\{z_{i},u(z_{i})\\}^{N}_{i=1}$ and
$\\{z_{i},u^{\prime}_{m}(z_{i})\\}_{i=1}^{N}$ for all $m=1,2,\ldots,d)$. Here
$u^{\prime}_{m}(z)=\frac{\partial f(z)}{\partial z^{(m)}}$ stands for the
partial derivative with respect to the m-th variable $z^{(m)}$. We try to find
an interpolant of the form

$u_{N}(Z)=\sum^{N}_{j=1}c_{j}\Phi(\epsilon\|Z-z_{j}\|)-\sum^{d}_{m=1}\sum^{N}_{j=1}\beta_{m,j}\Phi^{\prime}_{m}(\epsilon\|Z-z_{j}\|)$

with appropriate radial basis functions $\Phi$ so that $u_{N}$ satisfies the
generalized interpolation conditions

$\lambda_{i}u=\lambda_{i}u_{N},\,i=1,\ldots,N(d+1).$

Here, $\lambda_{i}$ could denote point evaluation at the point $z_{i}$, or it
could denote evaluation of some derivative at the point $z_{i}$.

After enforcing the interpolation conditions the system matrix is given by

$\displaystyle\mathbf{B}=\begin{bmatrix}\mathbf{A}_{0,0}&\mathbf{A}_{0,1}&\cdots&\mathbf{A}_{0,d}\\\
\mathbf{A}_{1,0}&\mathbf{A}_{1,1}&\cdots&\mathbf{A}_{1,d}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathbf{A}_{d,0}&\mathbf{A}_{d,1}&\cdots&\mathbf{A}_{d,d}\\\ \end{bmatrix}$

where $\mathbf{A}_{m,n}$ is the $N\times N$ matrix with the $(i,j)$th element

(6)
$\displaystyle\left(\boldsymbol{A}_{m,n}\right)_{i,j}=\begin{cases}\Phi(\epsilon\|z_{i}-z_{j}\|),&m=n=0,\\\
-\Phi^{\prime}_{n}(\epsilon\|z_{i}-z_{j}\|),&m=0,n\neq 0,\\\
\Phi^{\prime}_{m}(\epsilon\|z_{i}-z_{j}\|),&m\neq 0,n=0,\\\
-\Phi^{\prime\prime}_{m,n}(\epsilon\|z_{i}-z_{j}\|),&m\neq 0,n\neq
0,\end{cases}$

and $\mathbf{B}$ is a symmetric $(d+1)N\times(d+1)N$ matrix.

Note that the matrix $\mathbf{B}$ coincides with the joint covariance matrix
of the gradient-enhanced GP emulator using the covariance kernel $K$, see B.
Similar to the gradient-enhanced GP emulator, the condition numbers of the
matrix $\mathbf{B}$ are much larger than those of the ordinary Lagrange
interpolation matrix $\mathbf{A}$ [14]. Figs. 1 and 2 present the condition
numbers of the system matrix for the design matrixes $\mathbf{A}$ and
$\mathbf{B}$. Clearly, the design matrix becomes more singular as the shape
parameters are closer to zero or the total number of points becomes large.

Figure 1. Curves of the condition numbers for the Lagrange interpolation
(solid curves) and the Hermite interpolation (dashed curves) w.r.t. the shape
paramter using Guassian, IMQ and CS, respectively. Left: $N=50$; Middle:
$N=100$; Right: $N=300$; $d=2$.

Figure 2. Curves of the condition numbers for the Lagrange interpolation
(solid curves) and the Hermite interpolation (dashed curves) w.r.t. the number
of sampling points $N$ using Guassian, IMQ and CS, respectively. Left:
$\epsilon=0.1$; Middle: $\epsilon=1$; Right: $\epsilon=3$; $d=2$.

##  3\. Methodology

It is obvious that design of the collocation points plays a important role in
kernel interpolation. However, selecting optimal locations for trial centers
in $\Xi$ is highly non-trivial. These are the motivation of the development of
adaptive solvers [15, 17, 27]. Nevertheless, it is also believed that
uniformly distributed interpolation points can lead to better approximation
results [29]. Popular choices of such points include the Sobol’ points, the
low discrepancy points or even the random distributed points. However, it is
observed that all these choices have the instability issue when the number of
collocation points becomes large, say $N\sim\mathcal{O}(10^{2})$, see Fig. 2.

In this work, we shall propose a type of quasi-optimal collocation set for
kernel interpolations, by searching from a large set of candidate points. The
idea is inspired by the computation of “Fekete points” in polynomial
approximation.

###  3.1. Fekete-type interpolation points

We recall the kernel interpolation (3)

$\displaystyle
u_{N}(Z)=\sum^{N}_{j=1}c_{j}\,\Phi(\epsilon\|Z-z_{j}\|)=\sum_{j=1}^{N}d_{j}\,\ell_{j}\,(Z),$

where the latter equality writes the solution in terms of cardinal Lagrange
interpolants, given by

(7) $\displaystyle\ell_{j}\,(Z)$
$\displaystyle=\frac{\det\mathbf{A}(z_{1},...,z_{j-1},Z,z_{j+1},...,z_{N})}{\det\mathbf{A}(z_{1},...,z_{j-1},z_{j},z_{j+1},...,z_{N}))},$
$\displaystyle j$ $\displaystyle=1,...,N.$

The Lebesgue constant, corresponding to the operator norm of interpolation, is
defined as

$\displaystyle\Lambda_{N}:=\max_{Z\in
I_{Z}}\sum_{j=1}^{N}\left|\ell_{j}(Z)\right|.$

Small values of $\Lambda_{N}$ indicate that interpolation is a stable
operation. The so called Fekete points are a configuration of
$(z_{1},\ldots,z_{N})$ that maximizes the Vandermonde-like matrix determinant,
namely,

(8) $\displaystyle\Xi^{*}=\max_{\Xi=\\{z_{1},...,z_{N}\\}\subset
I_{Z}}\Big{\\{}\left|\det\mathbf{A}(z_{1},...,z_{N})\right|\Big{\\}}$

With the relation (7), one observes that the configuration $\Xi^{\ast}$
ensures that $|\ell_{j}(Z)|\leq 1$ for all $Z\in I_{Z}$, and thus ensures
$\Lambda_{N}\leq N$. In practice, the behavior of Lebesgue constant for Fekete
points is frequently sublinear. Thus, Fekete points provide a strategy for
stable interpolation. Unfortunately, there is no known way to explicitly
characterize or compute these points outside of special polynomial
approximation cases. In addition, computationally solving the optimization
problem (8) is a daunting task. A typical approach is to relax the
optimization problem by seeking a greedy (iterative) solution, i.e.,

$\displaystyle z_{N+1}=\operatorname*{arg\,max}_{z\in
I_{Z}}\left|\det\mathbf{A}(z_{1},\ldots,z_{N})\right|.$

In the polynomial approximation literature such an appropriate is titled the
“approximate Fekete point” approach, and for kernel interpolation it is a
power function maximization approach. We will utilize this procedure in our
gradient-enhanced setting to choose centers for approximation.

###  3.2. Generating quasi-optimal points via greedy selection

In this section, we present a data-independent algorithm to select optimal
distribution of centers. The data-independent feature is useful for practical
problems, as it allows one to determine, prior to conducting expensive
simulations or experimentations, where to collect data samples.

The general goal is to achieve a set of determinant-maximizing (i.e., Fekete)
points associated to the gradient-enhanced design matrix $\boldsymbol{B}$.
However, direct optimization is likely to be computationally challenging since
the optimization problem is non-convex, and the dimension of the optimization
variables can be very large when a large interpolation grid is used.
Therefore, we opt for a greedy strategy, performing the sequential
optimization,

$\displaystyle z_{N+1}$ $\displaystyle\coloneqq\operatorname*{arg\,max}_{z\in
I_{Z}}\det\boldsymbol{B}(z_{1},\ldots,z_{N},z),$

where $(z_{1},\ldots,z_{N})\subset\mathbb{R}^{d}$ is a given set of $N$
centers. The procedure above defines an iterative scheme. Although in
principle computing this determinant requires construction of the matrix
$\boldsymbol{B}(z_{1},\ldots,z_{N},z)$, such an expensive construction can be
avoided by exercising low-rank updates of $LU$ factorizations. First we note
that given $z_{1},\ldots,z_{N}$, we have

$\displaystyle\boldsymbol{B}(z_{1},\ldots,z_{N},z)$
$\displaystyle\coloneqq\left(\begin{array}[]{cccc}\boldsymbol{A}_{0,0}&\boldsymbol{A}_{0,1}&\cdots&\boldsymbol{A}_{0,d}\\\
\boldsymbol{A}_{1,0}&\boldsymbol{A}_{1,1}&\cdots&\boldsymbol{A}_{1,d}\\\
\vdots&\vdots&\ddots&\vdots\\\
\boldsymbol{A}_{d,0}&\boldsymbol{A}_{d,1}&\cdots&\boldsymbol{A}_{d,d}\end{array}\right),$

where above $\boldsymbol{A}_{m,n}=\boldsymbol{A}_{m,n}(z_{1},\ldots,z_{N},z)$.
Then if $\boldsymbol{P}_{N+1}\in\mathbb{R}^{(N+1)(d+1)\times(N+1)(d+1)}$ is a
permutation matrix that maps element $(N+1)i$ to element $N(d+1)+i$ for
$i=1,\ldots,d+1$ (and retains the order for the remaining $N(d+1)$ elements in
the first $N(d+1)$ entries), then

$\displaystyle\boldsymbol{B}(z_{1},\ldots,z_{N},z)=\boldsymbol{P}_{N+1}\boldsymbol{B}\boldsymbol{P}_{N+1}^{-T}$
$\displaystyle=\left(\begin{array}[]{cc}\boldsymbol{B}_{N}&\boldsymbol{W}(z)\\\
\boldsymbol{W}^{T}(z)&\boldsymbol{B}(z)\end{array}\right)$

Then the determinant can be written in terms of the Schur complement of this
latter expression,

$\displaystyle\det\boldsymbol{B}(z_{1},\ldots,z_{N},z)$
$\displaystyle=\det\left(\boldsymbol{P}\boldsymbol{B}\boldsymbol{P}^{T}\right)=\det\left(\begin{array}[]{cc}\boldsymbol{B}_{N}&\boldsymbol{W}\\\
\boldsymbol{W}^{T}&\boldsymbol{B}(z)\end{array}\right)$
$\displaystyle=\det\boldsymbol{B}_{N}\det\left(\boldsymbol{B}(z)-\boldsymbol{W}^{T}\boldsymbol{B}_{N}^{-1}\boldsymbol{W}\right),$

where we have introduced the notation
$\boldsymbol{B}_{N}\coloneqq\boldsymbol{B}(z_{1},\ldots,z_{N})$. Thus,

$\displaystyle\operatorname*{arg\,max}_{z\in
I_{Z}}\det\boldsymbol{B}(z_{1},\ldots,z_{N},z)$
$\displaystyle=\operatorname*{arg\,max}_{z\in I_{Z}}F(z),$ $\displaystyle
F(z)$
$\displaystyle\coloneqq\det\left(\boldsymbol{B}(z)-\boldsymbol{W}^{T}\boldsymbol{B}_{N}^{-1}\boldsymbol{W}\right).$

For a fixed $z$, the major cost of evaluting $F$ involves computing the
solution $\boldsymbol{X}\in\mathbb{R}^{N(d+1)\times(d+1)}$ to the linear
system $\boldsymbol{B}_{N}\boldsymbol{X}=\boldsymbol{W}$. There is a
secondary, comparatively negligible, cost of computing the determinant of a
$(d+1)\times(d+1)$ matrix. For the linear system, we can write the solution in
terms of the Cholesky factor of the permuted version of $\boldsymbol{B}_{N}$,

$\displaystyle\boldsymbol{P}_{N}\boldsymbol{B}_{N}\boldsymbol{P}_{N}^{T}=\boldsymbol{L}_{N}\boldsymbol{L}_{N}^{T}$

so that, if $\boldsymbol{L}_{N}^{-1}$ is available and stored, then evaluation
of $F$ requires application of $\boldsymbol{L}_{N}^{-1}$, with a
$\mathcal{O}(d^{3}N^{2})$ cost, and we can rewrite the optimization as,

(9) $\displaystyle z_{N+1}=\operatorname*{arg\,max}_{z\in I_{Z}}F(z)$
$\displaystyle=\operatorname*{arg\,max}_{z\in
I_{Z}}\det\left(\boldsymbol{B}(z)-\boldsymbol{V}^{T}\boldsymbol{V}\right),$
$\displaystyle\boldsymbol{V}$
$\displaystyle\coloneqq\boldsymbol{L}_{N}^{-1}\boldsymbol{W}(z).$

To aid in continuing the iteration, an efficient update that transforms
$\boldsymbol{L}^{-1}_{N}$ into $\boldsymbol{L}^{-1}_{N+1}$ (the Cholesky
factor for the pivoted $\boldsymbol{B}_{N+1}$) is helpful. Again exercising
Schur complements, we find that

$\displaystyle\boldsymbol{P}_{N+1}\boldsymbol{B}_{N+1}\boldsymbol{P}_{N+1}^{T}=\boldsymbol{L}_{N+1}\boldsymbol{L}_{N+1}^{T}\hskip
10.0pt\Longrightarrow\hskip
10.0pt\boldsymbol{L}_{N+1}=\left(\begin{array}[]{cc}\boldsymbol{L}_{N}&\boldsymbol{0}_{N(d+1)\times(d+1)}\\\
\boldsymbol{V}^{T}&\widetilde{\boldsymbol{L}}\end{array}\right),$

where $\widetilde{\boldsymbol{L}}=\widetilde{\boldsymbol{L}}(z)$ is the
Cholesky factor of the Schur complement of the block $\boldsymbol{B}(z_{N+1})$
of the matrix $\boldsymbol{B}_{N+1}$,

(10)
$\displaystyle\widetilde{\boldsymbol{L}}\widetilde{\boldsymbol{L}}^{T}=\boldsymbol{B}(z_{N+1})-\boldsymbol{V}^{T}\boldsymbol{V}.$

Then a computation shows that

(13)
$\displaystyle\boldsymbol{L}^{-1}_{N+1}=\left(\begin{array}[]{cc}\boldsymbol{L}^{-1}_{N}&\boldsymbol{0}_{N(d+1)\times(d+1)}\\\
-\widetilde{\boldsymbol{L}}^{-1}\boldsymbol{V}^{T}\boldsymbol{L}_{N}^{-1}&\widetilde{\boldsymbol{L}}^{-1}\end{array}\right).$

In summary, given $\boldsymbol{L}_{N}^{-1}$ that identifies $z_{N+1}$ from the
optimization (9), then we compute $\boldsymbol{L}_{N+1}$ by (i) computing
$\boldsymbol{V}(z_{N+1})$ from (9), (ii) computing the inverse of the
$(d+1)\times(d+1)$ Cholesky factor $\widetilde{\boldsymbol{L}}(z_{N+1})$ of
the Schur complement from (10), (iii) assembling the matrix in (13), which
requires an additional application of $\boldsymbol{L}_{N}^{-1}$.

Finally, we emphasize again that in practice we do not compute the solution to
(9), but rather we solve the easier optimization problem that replaces the
continuum $I_{Z}$ by a discrete candidate set.

##  4\. Numerical tests

In this section we use several numerical tests to demonstrate the benefit of
the quasi-optimal interpolation points. In order to implement our proposed
method, we first pick $M$ random centers on $I_{Z}$ as a candidate set, then
we select $N$ optimal points from Section 3.2 as RBF centers. In the following
examples, we choose $M=10^{4}$.

For comparisons against other methods, we employ the following sampling
methods.

  * •

Random sampling: $N$ samples are taken randomly in $I_{Z}$;

  * •

Sobol sampling: choose $N$ values of the Sobol sequence;

  * •

Halton sampling: choose $N$ values of the Halton sequence.

To measure the performance of an approximation, we will use the $l_{2}$ error
(RMSE). Specifically given a set of $Q=1000$ random samples
$\Xi_{test}=\\{z_{i}\\}^{Q}_{i=1}\in I_{Z}$ and samples of the true function
$u(z_{i})$ and the kernel approximation $s(z_{i})$ we compute

(14)
$E_{l_{2}}=\Big{(}\frac{1}{Q}\sum^{Q}_{i=1}|s(z_{i})-u(z_{i})|^{2}\Big{)}^{1/2}.$

In what follows we will use the term “Cholesky” for our proposed algorithm,
“random” to denote random sampling points, “sobol” to denote Sobol points, and
“halton” to denote Halton points. Furthermore, in all examples that follow we
perform 50 trials of each procedure and report the mean results along with
$20\%$ and $80\%$ quantiles. For CS RBFs, we use the following Wendlands
compactly supported function $\Phi_{d}(r)$, which has the form

$\Phi_{d}(r)=(1-r)_{+}^{l+3}\Big{(}(l^{3}+9l^{2}+23l+15)r^{3}+(6l^{2}+36l+45)r^{2}+(15l+45)r+15\Big{)}/15$

with $l=\left\lfloor\frac{d}{2}\right\rfloor+4.$

###  4.1. Kernel interpolation for UQ

In this section, some numerical examples are shown to demonstrate the
stability and convergent behaviors of our proposed algorithms for UQ.

####  4.1.1. Matrix stability

We first investigate the condition number
$Cond(\mathbf{A})=\frac{\sigma_{max}(\mathbf{A})}{\sigma_{min}(\mathbf{A})}$
of the design matrix $\mathbf{A}$. In Fig. 3, we show the condition number
with respect to the different shape parameters. In Fig. 4, we show the
condition number with respect to the number of the selected RBF centers $N$.
For both Gaussian and IMQ, our algorithm is much more stable compared to the
other sampling methods. In fact, the other choice of the candidates points do
not dramatically affect the performance of the proposed method. Fig.5 show the
corresponding results. In the figure, “uniform” utilizes the evenly spaced
points in $I_{Z}$.

Figure 3. Condition numbers with respect to shape parameters using Gaussian,
IMQ and CS. Left: $N=100$; Right: $N=300$; $d=2.$

Figure 4. Condition numbers with respect to shape parameters with respect to
the number of sample points $N$ using Gaussian, IMQ and CS. Left:
$\epsilon=3$; Right: $\epsilon=5$; $d=2$.

Figure 5. Condition number with respect to the number of sample points $N$ for
different choice of the candidates. Left: $\epsilon=3$; Right: $\epsilon=5$.

####  4.1.2. Numerical accuracy

In this section, we will compare Cholseky, random, sobol and halton in terms
of their ability to approximate test function.

We first consider the benchmark Franke’s function
$u(z),z=(z_{1},z_{2})\in[0,1]^{2}$, given by

(15)
$\displaystyle\begin{split}u(z)=&\frac{3}{4}e^{-((9z^{(1)}-2)^{2}+(9z^{(2)}-2)^{2})/4}+\frac{3}{4}e^{-(9z^{(1)}+1)^{2}/49-(9z^{(2)}+1)^{2})/10}\\\
&+\frac{1}{2}e^{-((9z^{(1)}-7)^{2}+(9z^{(2)}-3)^{2})/4}-\frac{1}{5}e^{-(9z^{(1)}-4)^{2}-(9z^{(2)}-7)^{2}}.\end{split}$

The corresponding numerical errors for the test function in (15) are shown in
Fig. 6. We observe that our proposed algorithm produces superior results than
the other sampling methods. Moreover, the new algorithm shows clear
convergence patterns as $N$ increases. This is considered to be an improvement
over the other sampling methods: it is clear in the error profile of the other
sampling methods that providing more sampling points does not always result in
better accuracy. Again, the choice of the candidates points do not
dramatically affect the performance of the proposed method. Fig.7 show the
corresponding results.

Figure 6.  Approximation error with respect to the number of sample points $N$
using Gaussian, IMQ and CS. Left: $\epsilon=3$; Right: $\epsilon=5$.

Figure 7. Numerical results with respect to the number of sample points $N$
for different choice of the candidates: Left: $\epsilon=3$; Right:
$\epsilon=5$.

Figure 8. Numerical results with respect to the number of sample points $N$
for different shape parameters. Left: condition number; Right: errors.

In Fig.8 we show the numerical results as a function of the sample points,
$N$, obtained using the Cholesky algorithm with various values of the shape
parameter $\epsilon$. The left plot shows the condition number while the right
plot shows the corresponding errors. It can be seen from this figure that as
$\epsilon$ gets smaller, the resulting matrix becomes ill-conditioned. The
numerical results when the LOOCV algorithm is used are also presented in Fig.
8. From these figures, we observe that the Cholesky algorithm using LOOCV to
choose the good value of shape parameter can provide a very good numerical
results.

We now consider the stochastic elliptic equation, one of the most used
benchmark problems in UQ, in one spatial dimension,

(16)
$-\frac{d}{dx}[\kappa(x,z)\frac{du}{dx}(x,z)]=f,\quad(x,z)\in(0,1)\times\mathbb{R}^{d},$

with boundary conditions

$\displaystyle u(0,z)=0,\quad u(1,z)=0,$

and $f=2$. The random diffusivity takes the following form

(17) $\kappa(x,z)=1+\sigma\sum^{d}_{k=1}\frac{1}{k^{2}\pi^{2}}\cos(2\pi
kx)z^{(k)},$

where $z=(z^{(1)},\dots,z^{(d)})$ is a random vector with independent and
identically distributed components.

Here we approximate the solution $u(z)=u(0.5,z)$, while the uncertain inputs
$z^{(i)}\sim U[-1,1],k=1,\cdots,d$. For further comparison, we also employ the
sparse grids stochastic collocation method using Legendre polynomial of total
order $k=8(d=3)$ and $k=4(d=6)$ to solve the problem. The numerical results
are presented in Figs. 9 and 10. The RBFs approximation methods, are notably
superior to the sparse grids method.

Figure 9.  Approximation error against the number of sample points $N$ (d=3).
Left: Gaussian; Middle: IMQ; Right: CS.

Figure 10.  Approximation error against the number of sample points $N$ (d=6).
Left: Gaussian; Middle: IMQ; Right: CS.

###  4.2. Gradient-enhanced approach for UQ

We now illustrate the symmetric approach to Hermite interpolation with a set
of numerical experiments.

Example 1:  We use the corner peak

$u(z)=(1+\sum^{d}_{i=1}\omega_{i}z^{(i)})^{-(d+1)},\,z\in\Omega=[0,1]^{d}$

to generate the data. In this example, we use
$d=2,\omega_{i}=\frac{1}{i^{2}}$.

Example 2: We consider the following 2-dimensional Rastrigin function:

$u(z)=20+\sum^{2}_{i=1}\Big{(}(z^{(i)})^{2}-10\cos(2\pi z^{(i)})\Big{)},\quad
z^{(i)}\in[-4,4].$

The 2-D Rastrigin function is plotted in Fig. 11.

Figure 11. 2-dimensional Rastrigin function.

Example 3: Consider the following 5-dimensional Friedman function

$u(z)=10\sin(\pi z^{(1)}z^{(2)})+20(z^{(3)}-0.5)^{2}+10z^{(4)}+5z^{(5)}.$

The corresponding numerical results are shown in Figs. 12 \- 14. We can
conclude that the design matrix $\mathbf{B}$ can be well conditioned under the
proposed algorithm. Again, our proposed algorithm produces superior results
than the other sampling methods.

Figure 12. Numerical results with respect to the number of sample points $N$
for 2-dimensional corner peak.

Figure 13. Numerical results with respect to the number of sample points $N$
for 2-dimensional Rastrigin function.

Figure 14. Numerical results with respect to the number of sample points $N$
for 5-dimensional Friedman function.

##  5\. Summary

In this work, we presented a strategy for selecting a quasi-optimal sample
points for kernel interpolation. We first demonstrate that the traditional
sampling methods suffer from the numerical instability in the sense that the
condition number will increase fast as the total number of points becomes
large. To improve this, we propose to use the efficient Cholesky decomposition
with pivoting for the Vandermonde-like interpolation matrix to choose the
sample points. Then the quasi-optimal sample points are used to control the
condition number of the design matrices. It is demonstrated that with the new
approach the stability can be much improved. On the other hand, for problems
inclusion of gradient measurements in the kernel interpolation are also
considered. Applications to parametric UQ problems are illustrated.

##  Appendix A Gaussian process and kernel interpolation

In order to help readers get deeper insights into the connections between
kernel methods and Gaussian process regressions, we give a simple example in
this section. More detains can be found in [30].

Suppose that we have a set of observed data values
$\\{(z_{j},u_{j})_{j=1}^{N}|z_{j}\in\Xi\\}$. We aim to estimate the unknown
value at some other locations $Z$ based on the observed data at $\Xi$. In the
Gaussian process regression [24, 30], a.k.a the simple kriging method, $u$ is
assumed as a realization of a random field $\mathbf{Y}$, which is a collection
$\\{\mathbf{Y}(Z):Z\in I_{Z}\\}$ of random variables over the probability
space $(\Omega,\mathcal{F},\mathcal{P})$. The observations $\\{u(z_{j})\\}$
are then realizations of the random variables $\mathbf{Y}(z_{j})$. To predict
$\mathbf{Y}$ at some location $Z\in I_{Z}$, one considers all linear
predictors of the form

$\displaystyle\mathbf{Y}(Z)=\sum^{N}_{i=1}\omega_{i}(Z)\mathbf{Y}(z_{i})$

which are themselves random variables.

To determine optimal weights
$\mathbf{\omega(Z)}:=(\omega_{1}(Z),\cdots,\omega_{N}(Z))^{T}$, additional
structural assumptions on $\mathbf{Y}$ are needed. In this work, we only
discuss the simple kriging, due to their close connection to kernel
interpolation. Suppose $\mathbf{Y}$ is a Gaussian process with a mean 0 and
symmetric covariance kernel $K$, i.e. $\mathbb{E}(\mathbf{Y}(Z))=0$,
$Cov(\mathbf{Y}(z_{i}),\mathbf{Y}(z_{j}))=K(z_{i},z_{j}):=\phi(z_{i}-z_{j})$.
We usually assume that such covariance kernel $K$ is positive definite. The
simple kriging method provides the best linear unbiased estimator $\hat{f}$ of
$\mathbf{Y}$ in the form of

$\hat{f}:=\omega_{1}(Z)u_{1}+\cdots+\omega_{N}(Z)u_{N}$

where the weighting $\mathbf{\omega(Z)}$ is uniquely determined by the linear
system ([30])

$\displaystyle\begin{bmatrix}K(z_{1},z_{1})&\cdots&K(z_{1},z_{N})\\\
\vdots&\ddots&\vdots\\\
K(z_{N},z_{1})&\cdots&K(z_{N},z_{N})\end{bmatrix}\begin{bmatrix}\omega_{1}(Z)\\\
\vdots\\\ \omega_{N}(Z)\end{bmatrix}=\begin{bmatrix}K(Z,z_{1})\\\ \vdots\\\
K(Z,z_{N})\end{bmatrix}$

or in matrix form $K(\Xi,\Xi)\mathbf{\omega(Z)}=K(Z,\Xi)^{T}$. Thus, the
estimator $\hat{f}(Z)$ can be written as

(18) $\displaystyle\hat{f}(Z):=\mathbf{\omega(Z)u}$

Note that if the RBF kernel and the covariance kernel coincide, then (3) and
(18) suggest that
$u_{N}(Z)=K(Z,\Xi)^{T}K(\Xi,\Xi)^{-1}\mathbf{u}=(K(\Xi,\Xi)^{-1}K(Z,\Xi))^{T}\mathbf{u}=\hat{f}(Z)$.
It means that the RBF estimator and the unbiased estimator are indeed
identical.

##  Appendix B Gradient-enhanced Gaussian process

If the gradient information
$\\{z_{i},u^{\prime}_{m}(z_{i})\\}_{i=1}^{N}\,(m=1,2,\ldots,d)$ are included
in the ordinary GP model, the covariance matrix becomes a block matrix that
includes the covariance between derivative observations in addition to the
covariance between function observations. The block matrix can be represented
as [18]:

$\displaystyle\mathbf{K}=\begin{bmatrix}Cov(\mathbf{Y},\mathbf{Y})&Cov(\mathbf{Y},\nabla\mathbf{Y})\\\
Cov(\nabla\mathbf{Y},\mathbf{Y})&Cov(\nabla\mathbf{Y},\nabla\mathbf{Y})\end{bmatrix},$

where $Cov(\mathbf{Y},\mathbf{Y})$ represents an $N\times N$ covariance matrix
of the function values at the sample points,
$Cov(\mathbf{Y},\nabla\mathbf{Y})$ is an $Nd\times N$ cross covariance matrix
between gradient components and the function values at the sample points, and
$Cov(\nabla\mathbf{Y},\nabla\mathbf{Y})$ is an $Nd\times Nd$ covariance matrix
of the gradients. If $\mathbf{Y}$ is a Gaussian process with a mean 0 and
symmetric covariance kernel $K$ (or $\phi$), then the joint covariance matrix
of the outputs is given by

$\displaystyle\mathbf{K}=\begin{bmatrix}\mathbf{P}_{0,0}&\mathbf{P}_{0,1}&\cdots&\mathbf{P}_{0,d}\\\
\mathbf{P}_{1,0}&\mathbf{P}_{1,1}&\cdots&\mathbf{P}_{1,d}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathbf{P}_{d,0}&\mathbf{P}_{d,1}&\cdots&\mathbf{P}_{d,d}\\\ \end{bmatrix}$

where $\mathbf{P}_{k,l}$ is the $N\times N$ matrix with the $(i,j)$th element

(19) $\displaystyle\begin{cases}\phi(z_{i}-z_{j}),&k=l=0,\\\
-\phi^{\prime}_{l}(z_{i}-z_{j}),&k=0,l\neq 0,\\\
\phi^{\prime}_{k}(z_{i}-z_{j}),&k\neq 0,l=0,\\\
-\phi^{\prime\prime}_{k,l}(z_{i}-z_{j}),&k\neq 0,l\neq 0,\end{cases}$

and $\mathbf{K}$ is a symmetric $(d+1)N\times(d+1)N$ matrix.

Note that if the RBF kernel $\Phi$ and the covariance kernel $\phi$ coincide,
then the joint covariance matrix $\mathbf{K}$ and the design matrix
$\mathbf{B}$ are indeed identical.

## References

  * [1] Ilias Bilionis and Nicholas Zabaras.  Multi-output local gaussian process regression: Applications to uncertainty quantification.  Journal of Computational Physics, 231(17):5718–5746, 2012.
  * [2] Ilias Bilionis and Nicholas Zabaras.  Bayesian uncertainty propagation using gaussian processes.  Handbook of Uncertainty Quantification, pages 1–45, 2016.
  * [3] A. Chkifa, A. Cohen, and C. Schwab.  High-dimensional adaptive sparse polynomial interpolation and applications to parametric PDEs.  Found. Comput. Math., 14(4):601–633, 2014.
  * [4] Jouke HS de Baar, Richard P Dwight, and Hester Bijl.  Improvements to gradient-enhanced kriging using a bayesian interpretation.  International Journal for Uncertainty Quantification, 4(3), 2014\.
  * [5] A. Doostan and H. Owhadi.  A non-adaptive sparse approximation for pdes with stochastic inputs.  J. Comput. Phys., 230(8):3015–3034, 2011.
  * [6] Gregory Fasshauer and Michael McCourt.  Kernel-based approximation methods using MATLAB, volume 19.  World Scientific Publishing Company, 2015.
  * [7] Gregory E. Fasshauer.  Meshfree approximation methods with MATLAB, volume 6 of Interdisciplinary Mathematical Sciences.  World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2007.
  * [8] Gregory E. Fasshauer and Jack G. Zhang.  On choosing “optimal” shape parameters for RBF approximation.  Numer. Algorithms, 45(1-4):345–368, 2007.
  * [9] B. Ganapathysubramanian and N. Zabaras.  Sparse grid collocation methods for stochastic natural convection problems.  J. Comput. Phys., 225(1):652–685, 2007.
  * [10] R. G. Ghanem and P. D. Spanos.  Stochastic finite elements: a spectral approach.  Springer-Verlag New York, Inc., 1991.
  * [11] L. Guo, A. Narayan, L. Yan, and T. Zhou.  Weighted approximate fekete points: sampling for least-squares polynomial approximation.  SIAM Journal on Scientific Computing, 40(1):A366–A387, 2018.
  * [12] Ling Guo, Yongle Liu, Akil Narayan, and Tao Zhou.  Sparse approximation of data-drien PCEs: an induced sampling approach.  Commun. Math. Res.,, 36:128–153, 2020.
  * [13] Ling Guo, Akil Narayan, and Tao Zhou.  Constructing least-squares polynomial approximations.  SIAM Review, 62(2):483–508, 2020.
  * [14] Xu He and Peter Chien.  On the instability issue of gradient-enhanced gaussian process emulators for computer experiments.  SIAM/ASA Journal on Uncertainty Quantification, 6(2):627–644, 2018\.
  * [15] Y. C. Hon, R. Schaback, and X. Zhou.  An adaptive greedy algorithm for solving large RBF collocation problems.  Numer. Algorithms, 32(1):13–25, 2003.
  * [16] J. Jackman, A. Narayan, and T. Zhou.  A generalized sampling and preconditioner scheme for sparse approximation of polynomial chaos expansions.  SIAM J. Sci. Comput, 39(3):A 1114 –1144, 2017.
  * [17] Leevan Ling and Robert Schaback.  An improved subspace selection algorithm for meshless collocation methods.  Internat. J. Numer. Methods Engrg., 80(13):1623–1639, 2009.
  * [18] Brian A Lockwood and Mihai Anitescu.  Gradient-enhanced universal kriging for uncertainty propagation.  Nuclear Science and Engineering, 170(2):168–195, 2012.
  * [19] Charles A Micchelli.  Interpolation of scattered data: Distance matrices and conditionally positive definite functions.  Constructive Approximation, 2(1):11–22, 1986.
  * [20] Max D Morris, Toby J Mitchell, and Donald Ylvisaker.  Bayesian design and analysis of computer experiments: use of derivatives in surface prediction.  Technometrics, 35(3):243–255, 1993.
  * [21] A. Narayan and J. Jakeman.  Adaptive Leja sparse grid constructions for stochastic collocation and high-dimensional approximation.  SIAM J. Sci. Comput., 36(6):A2952–A2983, January 2014.
  * [22] A. Narayan, J. D. Jakeman, and T. Zhou.  A Christoffel function weighted least squares algorithm for collocation approximations.  Math. Comput., 86:1913–1947, 2017.
  * [23] F. Nobile, R. Tempone, and C. G. Webster.  An Anisotropic Sparse Grid Stochastic Collocation Method for Partial Differential Equations with Random Input Data.  SIAM J. Numer. Anal., 46(5):2411–2442, January 2008.
  * [24] Carl Edward Rasmussen and Christopher KI Williams.  Gaussian process for machine learning.  MIT press, 2006.
  * [25] Shmuel Rippa.  An algorithm for selecting a good value for the parameter $c$ in radial basis function interpolation.  Adv. Comput. Math., 11(2-3):193–210, 1999.
  * [26] Jerome Sacks, Susannah B Schiller, and William J Welch.  Designs for computer experiments.  Technometrics, 31(1):41–47, 1989.
  * [27] Robert Schaback and Holger Wendland.  Adaptive greedy techniques for approximate solution of large RBF systems.  Numer. Algorithms, 24(3):239–254, 2000.
  * [28] Robert Schaback and Holger Wendland.  Kernel techniques: from machine learning to meshless methods.  Acta numerica, 15:543–639, 2006.
  * [29] Michael Scheuerer.  An alternative procedure for selecting a good value for the parameter $c$ in RBF-interpolation.  Adv. Comput. Math., 34(1):105–126, 2011.
  * [30] Michael Scheuerer, Robert Schaback, and Martin Schlather.  Interpolation of spatial data–a stochastic or a deterministic problem?  European Journal of Applied Mathematics, 24(4):601–629, 2013.
  * [31] Rohit Tripathy, Ilias Bilionis, and Marcial Gonzalez.  Gaussian processes with built-in dimensionality reduction: Applications to high-dimensional uncertainty propagation.  Journal of Computational Physics, 321:191–223, 2016.
  * [32] Selvakumar Ulaganathan, Ivo Couckuyt, Tom Dhaene, Joris Degroote, and Eric Laermans.  Performance study of gradient-enhanced kriging.  Engineering with computers, 32(1):15–34, 2016.
  * [33] Holger Wendland.  Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree.  Adv. Comput. Math., 4(4):389–396, 1995.
  * [34] Holger Wendland.  Scattered data approximation, volume 17 of Cambridge Monographs on Applied and Computational Mathematics.  Cambridge University Press, Cambridge, 2005.
  * [35] D. Xiu and J. S. Hesthaven.  High-order collocation methods for differential equations with random inputs.  SIAM J. Sci. Comput., 27(3):1118–1139, 2005.
  * [36] D. Xiu and G.E. Karniadakis.  The Wiener-Askey polynomial chaos for stochastic differential equations.  SIAM J. Sci. Comput., 24(2):619–644, 2002.
  * [37] L. Yan, L. Guo, and D. Xiu.  Stochastic collocation algorithms using $\ell_{1}$-minimization.  International Journal for Uncertainty Quantification, 2(3):279–293, 2012.
  * [38] F. L. Yang, L. Yan, and L. Ling.  Doubly stochastic radial basis function methods.  Journal of Computational Physics, 363:87–97, 2018.