# Partitioned Active Learning for Heterogeneous Systems

Cheolhei Lee a, Kaiwen Wang b, Jianguo Wu c, Wenjun Cai b, and Xiaowei Yue a
a Grado Department of Industrial and Systems Engineering,
Virginia Tech, Blacksburg, VA, USA
b Department of Materials Science and Engineering,
Virginia Tech, Blacksburg, VA, USA
c Department of Industrial Engineering and Management,
Peking University, Beijing, China

###### Abstract

Active learning is a subfield of machine learning that focuses on improving
the data collection efficiency of expensive-to-evaluate systems. Especially,
active learning integrated surrogate modeling has shown remarkable performance
in computationally demanding engineering systems. However, the existence of
heterogeneity in underlying systems may adversely affect the performance of
active learning. In order to improve the learning efficiency under this
regime, we propose the partitioned active learning that seeks the most
informative design points for partitioned Gaussian process modeling of
heterogeneous systems. The proposed active learning consists of two systematic
subsequent steps: the global searching scheme accelerates the exploration of
active learning by investigating the most uncertain design space, and the
local searching exploits the circumscribed information induced by the local
GP. We also propose Cholesky update driven numerical remedies for our active
learning to address the computational complexity challenge. The proposed
method is applied to numerical simulations and two real-world case studies
about (i) the cost-efficient automatic fuselage shape control in aerospace
manufacturing; and (ii) the optimal design of tribocorrosion-resistant alloys
in materials science. The results show that our approach outperforms benchmark
methods with respect to prediction accuracy and computational efficiency.

Keywords: Active Learning; Sequential Design; Gaussian Process; Heterogeneous
System

##  1 Introduction

Active learning is a subfield of machine learning and artificial intelligence
that maximizes information acquisition to train models data-efficiently
(Settles, 2009). Contrary to nonsequential passive learning such as Latin
Hypercube design (LHD) and factorial design (Santner et al., 2018), active
learning sequentially selects design points in the modeling phase after
observing intermediate models and outputs. It is also called query learning,
sequential design, adaptive sampling, or optimal design in other literature,
while they pursue the same objective: finding the best subset of inputs from
the design space according to information criteria that evaluates the
informativeness of input referring to uncertainty, disagreement, etc. Active
learning has got increasing attention from various applications where sampling
is timely and costly demanding, such as quality engineering, response surface
investigation, and image recognition (Alaeddini et al., 2019; Cao et al.,
2020). Especially, active learning for Gaussian processes (GPs) has been
vastly addressed in both methodological perspective and applications, and has
shown remarkable performance in practice (Seo et al., 2000; Gramacy and Lee,
2009; Yue et al., 2020).

Gaussian processes have been extensively utilized for modeling of various
systems spanning robotics, aerospace, and manufacturing process (Deisenroth et
al., 2013; Yue et al., 2018; Cho et al., 2019) due to the capability of
uncertainty quantification (UQ) and simplicity (Rasmussen and Williams, 2006).
Also, the UQ and diversity of kernel choice in GPs have enabled development of
various active learning algorithms, while they have been mostly studied with
single GPs (SGPs) that may be inappropriate for systems with heterogeneity
(e.g., discontinuity, and abrupt variations in gradient norms or frequencies).
Such heterogeneity is ubiquitous in industrial and engineering systems. For
example, composite materials, one of the most versatile materials in various
contemporary products, are anisotropic and highly nonlinear to external
treatments (Hyer and White, 2009), so they exhibit different behaviors in the
design space (Lee et al., 2020). Another example is the corrosive rates of
aluminum alloys. Fig. 1 illustrates the corrosive rates of alloys emulated
with the finite element method (FEM) over two pairs of control variables. We
can observe that the response surface shows spatial heterogeneity, so the
design space can be partitioned into three subregions according to the levels
of variation. In both cases, the efficiency of active learning driven by SGPs
can be significantly deteriorated, since SGPs cannot learn their heterogeneity
data-efficiently, thereby misleading the information criteria in active
learning.

Partitioned GPs (PGPs) can overcome the modeling limitations of SGPs in
heterogeneous systems by allocating multiple independent local GPs on disjoint
subregions. Subregions are defined or estimated according to distinguishable
characteristics of target systems, so PGPs can accommodate the heterogeneity
with local GPs. Moreover, the partitioning improves the scalability of GPs,
one of the main drawbacks, by introducing sparsity in their covariance
matrices. Due to these advantages, a few approaches for PGP modeling have been
proposed (Kim et al., 2005; Heaton et al., 2017; Pope et al., 2021). However,
most of them have been utilized without active learning, and the existing
algorithms for SGPs are suboptimal for PGPs in terms of learning and
computational efficiency. More explicitly, the SGP active learning strategy
does not take into account the heterogeneity in local GPs, so they are
computationally inefficient and cannot conduct macroscopic exploration.
Considering the predominance of heterogeneity in industrial and engineering
systems, PGPs are urgently in need of a bespoke active learning.

Figure 1: Corrosive rates of alloys with different control variables, Left:
response surface with Tafel slope and anodic exchange current, Right: response
surface with Young’s modulus and yield strength.

Motivated from the aforementioned limitations of existing active learning
strategies, we propose an active learning method for data acquisition and
modeling of heterogeneous systems. The proposed method seeks the most
uncertain subregion in the design space with the global searching scheme, and
the localized integrated mean squared error (IMSE) criterion is subsequently
utilized to select the most informative design point within the subregion. The
contributions of this paper are summarized as follows. First, we develop an
information criterion for capturing the heterogeneity of response surface.
Second, we propose novel active learning with a two-step searching scheme that
exploits the partitioned structure of models so as to improve the sampling
efficiency for the heterogeneous systems. Third, we propose Cholesky update
driven numerical remedies to reduce the computational cost of information
criterion in our active learning algorithm.

The remainder of this paper is organized as follows. In Section 2, we review
existing active learning methods for GPs, and discuss several PGP modeling
approaches with some active learning strategies therein. In Section 3, we
elucidate the new partitioned active learning algorithm, and provide
applicable techniques for improving the learning efficiency and numerical
costs. In Section 4, we implement our method on the approximation of two
simulation functions in different dimensions, predictive modeling of residual
stress in composite structures and corrosive rates of alloys in Section 5. To
evaluate our method in both simulation and case studies, we compare the
performance to other benchmark methods. Finally, a brief summary of this paper
is provided in Section 6.

##  2 Literature Review

In this section, we review related literature of existing active learning
strategies for SGPs and recent applications in real-world cases. Afterward, we
discuss literature of the PGP method and some active learning schemes therein.

###  2.1 _Active Learning for Gaussian Processes_

Due to the capability of UQ in GPs, the predictive uncertainty has been
frequently involved in the construction of active learning criteria for GPs.
Seo et al. (2000) compared two criteria for GP regression. Both are suggested
by MacKay (1992) and Cohn et al. (1996), and called active learning Mackay
(ALM) and active learning Cohn (ALC), respectively. The ALM refers to the
variance-based criterion that selects unobserved points with the largest
predictive uncertainty quantified with variance or standard deviation. A
number of applications exist that employ the ALM due to its simplicity and
straightforwardness (Uteva et al., 2018). Meanwhile, the ALC refers to the
integrated mean squared prediction error (IMSE) that seeks a point expected to
reduce the most the integrated variance of the GP model over the design space,
and it has also been widely applied to different GP frameworks. Chen and Zhou
(2017) utilized the IMSE criterion for stochastic kriging (Ankenman et al.,
2010), so as to balance between exploration over the design space and
exploitation with additional replications. Yue et al. (2020) also proposed two
active learning algorithms for GP surrogate models of multimodal systems. The
variance-based weighted active learning uses the weighted sum of variance-
based criteria of modes, and the D-optimal weighted active learning uses the
fisher information matrix with the same engineering-driven weighting
procedure. The Bayesian optimization (BO) (Snoek et al., 2012) is another
interesting trend of sequential sampling that utilizes the predictive
uncertainty of GPs for global optimization problems. Although the BO and
active learning show subtle nuance in their objectives, they share the common
principle, reduction in sampling cost. There are several criteria (also called
acquisition functions) for the BO, such as the expected improvement (EI) and
the probability of improvement (PI). They refer to the GP surrogate model of
blackbox objective function, thereby choosing a point to be queried most
likely to be optimal.

For strategies that do not take the predictive uncertainty into account or
less consider, Pasolli and Melgani (2011) suggested using space-filling
designs in the Kernel Hilbert space induced by intermediate GP models. Freytag
et al. (2013) and Käding et al. (2018) used expected model output change as
their criterion for active learning. The strategy chooses the location where
is expected to induce the largest change of model outputs over the design
space. Erickson et al. (2018) and Marmin et al. (2018) involve the gradient of
GP model’s response in active learning so as to draw more samples from where
the response abruptly changes. Consequently, the gradient-based criterion
focuses more on local variation than the entire design space, thereby reducing
the prediction error more efficiently. Uteva et al. (2018) and Kim et al.
(2009) referred to the discrepancy of multiple GP models trained on different
subsets of data. However, the aforementioned active learning strategies mostly
devised for SGPs that are inappropriate for modeling heterogeneous systems.
Moreover, it is inefficient to directly apply them to PGPs, since they cannot
consider the partitioned structure of PGP models. Liu et al. (2018) provided a
comprehensive active learning strategy for GP surrogate modeling of
engineering systems.

###  2.2 _Partitioned Gaussian Processes_

Although there are techniques other than partitioning design spaces to
overcome the limitation of stationary GP models such as input-dependent length
scale, warping and convolution kernels (Heinonen et al., 2016; Marmin et al.,
2018; Higdon, 1998), we mainly focus on PGPs that explicitly partition the
design space for multiple GP models. Kim et al. (2005) proposed the piecewise
GP using Voronoi tessellation with training dataset for partitioning the
design space. Advantages of Voronoi tessellation are simplicity, consistency,
and distance-based algorithm that coincides with stationary GPs. They
estimated the number of regions and centers with the Monte-Carlo approach, and
fitted independent GP models for subregions. Subsequently, Pope et al. (2021)
generalized the partitioning by merging convex Voronoi cells in order to
generate nonconvex subregions, and relaxed centers of cells. Gramacy and Lee
(2008) proposed the Treed GP (TGP) using decision trees for partitioning. They
fitted independent local GPs on each leave corresponds to partitioned
subregions. Heaton et al. (2017) proposed to partition the design space prior
to GP modeling by hierarchical clustering referring to finite differences of
samples. They insist that the approach allows avoiding the expensive Markov
Chain Monte-Carlo (MCMC) algorithm in the aforementioned approaches. The
mixture of GP experts proposed by Rasmussen and Ghahramani (2002) is another
approach using multiple stationary GPs. Although components in their modeling
process have own terminologies, the underlying idea is very close to the
aforementioned methodologies.

For active learning in PGPs, Pope et al. (2021) proposed to an active learning
algorithm with their PGP modeling that takes a point within regional
boundaries and maximizes the space-filling property of design points. However,
it focuses more on detecting discontinuity in design space rather than
reducing the prediction error over the design space, although such
discontinuity is rare in industrial and engineering applications. Gramacy and
Lee (2009) and Konomi et al. (2014) provided active learning algorithms that
are modified versions of Seo et al. (2000) for the TGP. They provided more
elaborate techniques considering the posterior structure of partitioned design
space. However, their approaches are highly dependent on the tree classifier,
while the tree partitioning mostly induces boundaries parallel to axes that
may not be realistic in practice (Heaton et al., 2017; Pope et al., 2021).
Moreover, the choice of design point candidates are dependent to the areas of
partitioned regions, so it can be irrelevant when a plausible partition is not
realizable with a few tree-partitioned regions. If we expand our field of
interest to classification with generic machine learning models, Cortes et al.
(2019) has a similar concept of ours in the macroscopic aspect, while there
are main differences: (i) they consider classification problems, thus
theoretical bounds they provided may be intractable in practice with
regression problems (e.g., model capacity, bounded loss); and (ii) they
consider the disagreement-based criterion within a model class, while we
advocate the IMSE criterion, which belongs to variance reduction criterion.

##  3 Partitioned Active Learning

In this section, we propose our active learning algorithm for the PGP
modeling. First, we briefly discuss a generic framework of PGPs as the
initialization step of active learning, and elucidate our strategy
consecutively. Lastly, applicable remedies by Cholesky update are proposed for
improving the computational complexity of the algorithm.

As mathematical notations, we use lower letters for scalars, and distinguish
vectors with boldface. Upper letters indicate sets or matrices, and a set of
subregion indices are denoted as $[M]=\\{1,\ldots,M\\}$. In subscriptions, we
use parentheses numbers for indices of subregions, and the subsequent normal
letter as the indices of letter in the set or matrix. If the index of
subregion is apparent, we may omit the region index and only use the index in
the set.

###  3.1 _Surrogate Modeling for Heterogeneous Systems_

We aim to build a PGP model of a heterogeneous system which is expensive to
operate. Suppose that the system can be expressed as a target function $f$
defined over a design space $\Omega\subset\mathbb{R}^{d}$ mapping to
$\mathbb{R}$, and the function is endowed with the heterogeneity in $\Omega$.
Let $h:\Omega\to\mathbb{R}$ be our estimated PGP model of $f$, and
$g:\Omega\to[M]$ be the region classifier that partitions $\Omega$ into $m$
disjoint regions such that $\Omega=\bigcup_{m=1}^{M}\Omega_{(m)}$ accordance
with the heterogeneity. Then, the PGP composed of $M$ local GPs is

$\displaystyle
h(\mathbf{x})=\sum_{m=1}^{M}1_{\\{g(\mathbf{x})=m\\}}h_{(m)}(\mathbf{x}),\qquad\mathbf{x}\in\Omega,$
(1)

where $1_{\\{C\\}}$ is an indicator function which has value 1 when $C$ is
true and 0 otherwise, and $h_{(m)}$ is a local GP assigned to $\Omega_{(m)}$.
Although PGPs can take any valid kernel, we mainly consider stationary kernel
family such as radial basis function (RBF) and Matérn. Suppose any
$\mathbf{x},\,\mathbf{x}^{\prime}\in\Omega$, and $\mathbf{x}=\left[\>x_{1}\
\ldots\ x_{d}\>\right]^{\top}$. The local GP defined over $\Omega_{(m)}$ with
the RBF kernel is defined as

$\displaystyle
h_{(m)}(\mathbf{x})\sim\mathcal{GP}\left(\mu_{(m)}(\mathbf{x}),\
k_{(m)}(\mathbf{x},\,\mathbf{x}^{\prime})\right),\quad\mathbf{x},\,\mathbf{x}^{\prime}\in\Omega_{(m)},$
(2) $\displaystyle
k_{(m)}(\mathbf{x},\,\mathbf{x}^{\prime})=\tau_{(m)}^{2}\prod_{j=1}^{d}\exp\left(-\frac{(x_{j}-x_{j}^{\prime})^{2}}{l_{(m),j}^{2}}\right)+\sigma^{2}_{(m)}1_{\\{\mathbf{x}=\mathbf{x}^{\prime}\\}},$
(3)

where $\mu_{(m)}(\mathbf{x})$ is the mean function assumed to be 0 without
loss of generality, and $k_{(m)}$ is the kernel in which nonnegative
$\tau_{(m)},\ l_{(m),j}$ and $\sigma^{2}$ are referred as scale, length and
noise hyperparameters. Note that $\mathbf{l}_{(m)}=\left[\>l_{(m),1}\ \ldots\
l_{(m),d}\>\right]^{\top}$ implies that the local GP implements automatic
relevance determination (ARD) (Rasmussen and Williams, 2006). We use $\Theta$
as the set of hyperparameters of PGP, that is
$\Theta=\\{\tau,\,\mathbf{l},\,\sigma\\}$, and
$\Theta_{(m)}=\\{\tau_{(m)},\,\mathbf{l}_{(m)},\,\sigma_{(m)}\\}$.

Suppose we have finite $n$ observations on
$X=\\{\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\\}\subset\Omega$ such that
$D=\\{(\mathbf{x}_{1},y_{1}),\ldots,(\mathbf{x}_{n},y_{n})\\}$, where $y_{i}$
is the observation at $\mathbf{x}_{i}$ with noise,
$=f(\mathbf{x}_{i})+\epsilon$ and $\epsilon\sim\mathcal{N}(0,\ \sigma^{2})$.
According to the region classifier, $X$ and $D$ can be partitioned as
$X_{(m)}=\\{\mathbf{x}_{(m),1},\ldots,\mathbf{x}_{(m),n}\\}$, and
$D_{(m)}=\\{(\mathbf{x}_{(m),1},y_{(m),1}),\ldots,(\mathbf{x}_{(m),1},y_{(m),1})\\}$.
Let the covariance matrix associated with $X_{(m)}$ be $K_{(m)}$ such that
$(K_{(m)})_{ij}=k_{(m)}(\mathbf{x}_{(m),i},\,\mathbf{x}_{(m),j})$ for
$i,\,j\in[n]$. The PGP model can be trained with $D$ by maximizing the log
marginal likelihoods of local GPs, which is proportional to

$\displaystyle\ell(\Theta)\propto\sum_{m=1}^{M}\left(\mathbf{y}_{(m)}^{\top}K_{(m)}^{-1}\mathbf{y}_{(m)}-\log\det{K_{(m)}}-n_{(m)}\log{2\pi}\right),$
(4)

where $\mathbf{y}_{(m)}=\left[\>y_{(m),1}\ \ldots\ y_{(m),n}\>\right]^{\top}$,
and $n_{(m)}$ is the number of samples in $D_{(m)}$. Note that the PGP’s
likelihood is the multiplication of that of local GPs due to their
independence. Consequently, possibly with some ordering process, the PGP
produces a block diagonal covariance matrix which implies that the numerical
advantage of PGPs comes from the sparsity. Although the construction of entire
covariance matrix is generally unnecessary in practice, it informs us that the
model can be manipulated more efficiently by treating each local GP
independently.

The performance of PGP is determined by the credibility of the region
classifier $g$ and estimated hyperparameters of local GPs
$\boldsymbol{\Theta}$. When the underlying partitioning rule is unknown, $g$
must be estimated, which is the major focus of most PGP approaches. The MCMC
approach is generally employed, while the expensive MCMC inference can be
avoided by clustering finite differences (Heaton et al., 2017) and input-
output pairs (Nguyen and Smeulders, 2004). We assume that the region
classifier $g$ is already available or can be estimated via one of the
aforementioned approaches, thereby focusing more on the estimation of
$\Theta$. Since $\Theta$ would be inferred by maximizing Equation (4) given
$D_{n}$, which is generally realized with nonsequential designs, it is unwise
to exhaust the sampling budget at the first state. Hence, we take only a small
portion in the initial step.

Unfortunately, there is no universal concrete theorem for the optimal portion
of initial sampling, while some empirical suggestions can be found in (Liu et
al., 2018; Yue et al., 2020). However, we can conjecture that the number of
initial samples has a trade-off property. If the initial sampling is weighed
too much, the advantage of active learning will be diluted. Otherwise, active
learning can be hindered by low-quality information come from unreliable
intermediate models. Also, when the region classifier $g$ is not explicitly
known so that $g$ must be estimated with the initial samples, the number
should be enough to obtain an acceptable $g$.

###  3.2 _Partitioned Active Learning Strategy_

The essence of active learning is an information criterion function
$J:\Omega\to\mathbb{R}$ that quantifies the potential importance of a design
point. By optimizing $J(x)$ with respect to $x\in\Omega$, the next design
point is determined and queried to the target function, which is also called
the oracle in active learning. The variance-based and the IMSE criteria are
mostly considered in active learning for GPs due to their versatility and
simplicity, so we shortly summarize two criteria and tailor them to establish
a new criterion for PGPs.

Suppose we have $D\subset\Omega$, and intend to determine the next location
$\mathbf{x}_{n+1}\in\Omega$ with GPs without partitioning. The variance-based
criterion is defined as

$\displaystyle
J_{\text{VAR}}(\mathbf{x})=\text{Var}_{n}(\mathbf{x}):=k(\mathbf{x})-\mathbf{k}(\mathbf{x},\,X)^{\top}K^{-1}\mathbf{k}(\mathbf{x},\,X),$
(5)

where $\text{Var}_{n}(\mathbf{x})$ is the predictive variance at $\mathbf{x}$
given $n$ observations $D$. The ALM maximizes the variance-based criterion so
as to select the location with the greatest predictive uncertainty. Meanwhile,
the IMSE criterion is calculated by

$\displaystyle
J_{\text{IMSE}}(\mathbf{x})=\int_{\Omega}\text{Var}_{n+1}(\mathbf{s}|\mathbf{x})\,p(\mathbf{s})\,d\mathbf{s},$
(6)
$\displaystyle\text{Var}_{n+1}(\mathbf{s}|\mathbf{x})=k(\mathbf{s})-\mathbf{k}(\mathbf{s},\,X_{n+1})^{\top}K_{n+1}^{-1}\,\mathbf{k}(\mathbf{s},\,X_{n+1}),\quad
X_{n+1}=\left[\>X_{n}\ \ \mathbf{x}\>\right]^{\top},$

where $\mathbf{s}\in\Omega$ with density (or importance) $p(\mathbf{s})$, and
$K_{n+1}$ is the covariance matrix associated with $X_{n+1}$. Minimizing the
IMSE criterion selects the location which is expected to reduce the predictive
uncertainty the most over $\Omega$, and we refer the active learning with the
IMSE criterion as the ALC.

There are more behind derivations of both criteria, while we mention shortly
herein. It turns out that the ALM is equivalent to the maximum entropy design,
since the choice leads to maximizing the determinant of covariance matrix
given $\Theta_{n}$. Meanwhile, the ALC can be explained by minimizing the
generalization error of statistical learning, which can be decomposed into the
famous form,

$\displaystyle\mathbb{E}_{\Omega}\left[\|f-h\|^{2}\right]=\mathbb{E}_{\Omega}\left[\|f-h^{\ast}\|^{2}\right]+\mathbb{E}_{\Omega}\left[\|h^{\ast}-h\|^{2}\right],$

where $h^{\ast}$ is the best estimation of $f$ in the given PGP class.
Assuming that the first term representing the bias of the best model is
acceptable, the ALC focuses on minimizing the second term, which is the
variance. Regarding the variance of $h^{\ast}$ as a constant, it becomes
minimizing the overall variance of $h$.

Although the variance-based criterion is more straightforward and numerically
inexpensive than the IMSE criterion, the ALC empirically has shown better
performance than the ALM (Seo et al., 2000; Gramacy and Lee, 2009). Moreover,
one of the promising characteristic of the ALC is that it avoids sampling from
the boundary of design space which may provoke loss of information, while the
ALM frequently does. It makes the IMSE criterion more appealing to PGPs, since
local GPs share common boundaries, so the variance-based criterion may lead to
oversampling near the boundary shared with two adjacent local GPs. Another
advantage of the ALC is that it is possible to consider the importance of
$\mathbf{s}\in\Omega$ in the information criterion, while the ALM cannot. It
allows us to incorporate prior knowledge and to give more weight on a specific
region, thereby making the algorithm more distinguishable from the space-
filling design.

However, the IMSE criterion of (6) does not consider multiple local GPs
simultaneously, thus the block diagonal structure of covariance matrix is not
involved. Therefore, we dedicate to construction of our criterion for the PGPs
by aggregating every information criterion attained by each local GPs. When
the IMSE criterion is considered for a candidate location
$\mathbf{x}\in\Omega_{(m)}$ with PGPs, it can be written as

$\displaystyle J(\mathbf{x}|\mathbf{x}\in\Omega_{(m)})=$
$\displaystyle\sum_{i\neq
m}\int_{\Omega_{(i)}}\text{Var}_{(i)}(\mathbf{s}_{(i)})\,p(\mathbf{s}_{(i)})\,d\mathbf{s}_{(i)}$
$\displaystyle+\int_{\Omega_{(m)}}\text{Var}_{(m),n_{(m)}+1}(\mathbf{s}_{(m)}|\mathbf{x})\,p(\mathbf{s}_{(m)})\,d\mathbf{s}_{(m)},$
(7)

where $\mathbf{s}_{(m)}\in\Omega_{(m)}$ with density $p(\mathbf{s}_{(m)})$.
The interpretation of each term in (7) is quite worthwhile. The first term is
the sum of IMSEs except for $h_{(m)}$, and it is invariant to the choice of
$\mathbf{x}\in\Omega_{(m)}$. The second term is equivalent to (6), in which
$h_{(m)}$ is only considered, so that (7) will only take account of the local
region associated with the candidate location. Consequently, there are two
main differences between (6) and (7): (i) consideration of IMSEs over other
local regions; and (ii) the localized IMSE criterion. We focus on each term
subsequently considering their meanings, thereby efficiently minimizing (7).

Heuristically, (7) is more likely to be minimized when the most uncertain
local GP is taken account into the second term, since the local GP has more
potential to be reduced with additional observations. Each IMSE in the first
term indicates the regional uncertainty of PGP, so it can be used for
investigating the most uncertain region. We involve every local GP in the
first term of (7), and conduct the global searching prior to calculating the
second term by using the following criterion

$\displaystyle J_{\text{G}}(m)$
$\displaystyle=\int_{\Omega_{(m)}}\text{Var}_{(m),n_{(m)}}(\mathbf{s}_{(m)})\,p(\mathbf{s}_{(m)})\,d\mathbf{s}_{(m)},$
$\displaystyle m^{\ast}$ $\displaystyle=\arg\max_{m\in[M]}J_{G}(m).$ (8)

A possible problem can be aroused in the global searching is that if multiple
regions have similar uncertainties, then this scheme will only care about the
most uncertain one. However, it can be easily resolved by taking a percentile
or providing a threshold to consider multiple regions in the global searching.

Once the most uncertain region is determined by (8), we focus on the second
term within $\Omega_{(m^{\ast})}$ as the local searching, which is

$\displaystyle J_{\text{L}}(\mathbf{x})$
$\displaystyle=\int_{\Omega_{(m^{\ast})}}\text{Var}_{(m^{\ast}),n_{(m^{\ast})}+1}(\mathbf{s}_{(m^{\ast})}|\mathbf{x})\,p(\mathbf{s}_{(m^{\ast})})\,d\mathbf{s}_{(m^{\ast})},$
$\displaystyle\mathbf{x}^{\ast}$
$\displaystyle=\arg\min_{\mathbf{x}\in\Omega_{(m^{\ast})}}J_{L}(\mathbf{x}).$
(9)

Since the local GP in (9) reflects the local behavior of underlying function
excluding heterogeneity from other regions, it can lead to improvement in
exploitation of active learning by avoiding implausible predictive
uncertainty. Moreover, it turns out that the localized IMSE criterion reduces
the computational cost of IMSE criterion due to the reduced size of local GP.
We call the sequential criteria (8) and (9) by Partitioned IMSE (PIMSE), and
the active learning with PIMSE as Partitioned ALC (PALC).

Two termination criteria can be considered in the algorithm. One is based on
the pre-defined budget, and the other is the achievement of desired prediction
accuracy with the model. The budget criterion is straightforward, while most
active learning algorithms may not guarantee that the last model is the best
one unless the trained dataset is noise-free. To avoid the unwanted
overfitting, early-stopping with a desired prediction accuracy is one possible
solution. Note that, to utilize the prediction accuracy criterion for the
termination, a separated testing dataset or cross validation scheme is
required for the model evaluation in the learning process.

Algorithm 1 PALC

1: Prerequisite:
$M,N_{0},N_{\text{Max}},N_{\text{Ref}},N_{\text{Cand}},\Omega,(\text{Optional:
}L,D_{\text{Test}},e^{\ast})$

2: Initialize $D_{0}$ with $N_{0}$ space-filling design over $\Omega$

3: Partition $\Omega$ into $M$ subregions

4: Train PGP ($h$) on $D_{0}$

5: while card($D$) $<N_{\text{Max}}$ and
$L_{D_{\text{Test}}}(h)>\epsilon^{\ast}$ do

6: Generate $X_{\text{Ref}}\subset\Omega$ with $N_{\text{Ref}}$ space-filling
design

7: Global Searching: Solve $m^{\ast}=\arg\max_{m}J_{\text{G}}(m)$

8: Generate $X_{\text{Cand}}\subset\Omega_{(m^{\ast})}$ with $N_{\text{Cand}}$
space-filling design

9: Local Searching: Solve $\mathbf{x}^{\ast}=\arg\min_{\mathbf{x}\in
X_{\text{Cand}}}J_{\text{L}}(\mathbf{x})$

10: Observe $y^{\ast}$ at $\mathbf{x}^{\ast}$

11: $D=D\cup\\{(\mathbf{x}^{\ast},y^{\ast})\\}$

12: Train $h$ on $D$

13: if $D_{\text{Test}}$ exists then

14: Check $L_{D_{\text{Cand}}}(h)<e^{\ast}$

15: end if

16: end while

The pseudocode of PALC is provided in Algorithm 1. In prerequisites,
$(L,D_{\text{Test}},e^{\ast})$ is optionally required for the use of
predictive accuracy criterion which of $L$ is a loss function,
$D_{\text{Test}}$ is a separated dataset, and $e^{\ast}$ is a desired
prediction accuracy. $N_{0}$ and $N_{\text{Max}}$ stand for the numbers of
initial samples and the maximum attainable samples (i.e., budget),
respectively. $X_{\text{Ref}}$ is a reference set composed of $N_{\text{Ref}}$
space-filling points
($\\{\mathbf{s}_{i}\\}_{i=1}^{N_{\text{Ref}}}\subset\Omega$), which is
required to implement the global searching approximately as

$\displaystyle
J_{\text{G}}(m)=\frac{1}{N_{Ref}}\sum_{i=1}^{N_{Ref}}\text{Var}_{(m),n_{(m)}}(\mathbf{s}_{i})\,\tilde{p}(\mathbf{s}_{(m),i}),$

where $\tilde{p}(\mathbf{s}_{(m),i})$ is the approximated probability mass
function at $\mathbf{s}_{(m),i}$. The subset of $X_{\text{Ref}}$ composed of
points in $\Omega_{(m^{\ast})}$, which is
$\\{\mathbf{s}_{i}\\}_{i=1}^{N_{\text{Ref}}^{\ast}}$ will be referred to the
approximated local searching subsequently as

$\displaystyle
J_{\text{L}}(\mathbf{x})=\frac{1}{N_{\text{Ref}}^{\ast}}\sum_{i=1}^{N_{\text{Ref}}^{\ast}}\text{Var}_{(m^{\ast}),n_{(m^{\ast})}+1}(\mathbf{s}_{(m^{\ast}),i}|\mathbf{x})\,\tilde{p}(\mathbf{s}_{(m^{\ast}),i}),$
(10)

where $\mathbf{x}\in
X_{\text{Cand}}=\\{\mathbf{x}_{i}\\}_{i=1}^{N_{\text{Cand}}}$. In this paper,
we have generated new $X_{\text{Ref}}$ and $X_{\text{Cand}}$ with LHD in every
step in order to encourage exploration.

###  3.3 _Cholesky Update based Numerical Remedies to Tackle the Computational
Complexity Challenge_

The IMSE criterion is numerically more demanding than the variance-based
criterion, since it involves the inversion of $K_{n+1}$ which should be
updated for every candidate. That is, calculation of the IMSE criterion
requires $\mathcal{O}(n^{3})$ for each candidate. Moreover, when $N$
candidates are provided to the active learning module, the computational cost
is multiplied by the number. Although the significance of their effects may
vary with situations, the effect of candidate number can be more considerable
than the inversion cost. Therefore, in order to improve the numerical aspect
of PALC, we should provide some remedies for both matrix inversion and the
number of candidates.

The global searching can be used for reducing the number of candidates by
narrowing down the region of interest. Generally, candidates for active
learning are given with space-filling or dense-grid over the design space.
Therefore, given that candidates are provided uniformly over the overall
design space, taking the subset of candidates in the most uncertain region
with the global searching leads to reduction in the number of matrix
inversions proportional to

$\displaystyle\text{Pr}(\mathbf{s}\in\Omega_{(m^{\ast})})=\int_{\Omega_{(m^{\ast})}}p(\mathbf{s}_{(m^{\ast})})\,d\mathbf{s}_{(m^{\ast})},$

which is the ratio of the most uncertain region considering the importance of
input.

The matrix inversion cost of PGP is automatically alleviated by partitioning
the design space with the block diagonal covariance matrix. That is, the
inversion cost reduces from $\mathcal{O}(n^{3})$ to at most
$\mathcal{O}({n_{(m)}}^{3})$, where ${n_{(m)}}<n$ usually. Another applicable
remedy is updating the inverse of $K_{n_{(m)}+1}$ in (9) exploiting
$K_{n_{(m)}}^{-1}$ iteratively. Although it is possible to apply the Sherman-
Morrison formula to get the updated inverse matrix directly as shown in
Gramacy and Lee (2009), we propose to use the Cholesky decomposition
(Rasmussen and Williams, 2006) for solving the linear system
$K_{n_{(m)}}^{-1}\mathbf{k}$ considering the numerical stability and cost.
Given that the Cholesky factor (a lower triangular matrix) of $K_{n_{(m)}}$ is
known, the update of Cholesky factor of $K_{n_{(m)}+1}$ with that of
$K_{n_{(m)}}$ only requires the forward substitution step of the size
$n_{(m)}$ triangular system, thus it needs only $\mathcal{O}(n_{(m)}^{2})$
instead of $\mathcal{O}({n_{(m)}}^{3})$. A more detailed procedure of the
Cholesky update is provided as follows.

#### _Updating Cholesky Factor for PIMSE Criterion_

Suppose we have the Cholesky factor $L_{n}$ of $K_{n}$, which is the
covariance matrix of $X_{n}$, such that $K_{n}=L_{n}{L_{n}}^{\top}$. We aim to
get the Cholesky factor $L_{n+1}$ of

$\displaystyle K_{n+1}=\left[\begin{array}[]{c
c}K_{n}&\mathbf{k_{n}^{\ast}}\\\
\mathbf{k_{n}^{\ast}}^{\top}&k(\mathbf{x}^{\ast})\end{array}\right],$ (13)

where $\mathbf{x}^{\ast}$ is a candidate input, and
$\mathbf{k}_{n}^{\ast}=k(X_{n},\,\mathbf{x}^{\ast})$. Since $K_{n+1}$ shares
the same part of $K_{n}$, it turns out that $L_{n+1}[:n,:n]$ is equivalent to
$L_{n}$. For $L_{n+1}$, we can apply the Cholesky-Banachiewicz algorithm for
$i=1,\ldots,n$ as

$\displaystyle
L_{n+1,i}=L_{i,i}^{-1}\left(\mathbf{k_{n}^{\ast}}_{i}-\sum_{j=1}^{i-1}L_{n+1,j}L_{i,j}\right),$
$\displaystyle
L_{n+1,n+1}=\sqrt{k(\mathbf{x}^{\ast})-\sum_{j=1}^{n}L_{n+1,j}}.$

Rather than calculating the PIMSE directly, the PALC can be faster by skipping
redundant computation applying the Cholesky updating approach. The Cholesky
factor $L_{n}$ can be used for the predictive variance of GP in (5) for
$\mathbf{s}\in X_{\text{Ref}}$, and the global searching criterion of (8) as

$\displaystyle\text{Var}_{n}^{2}(\mathbf{s})=k(\mathbf{s})-{\mathbf{v}_{n}}^{\top}{\mathbf{v}_{n}},\qquad\mathbf{v}_{n}=L_{n}^{-1}{\mathbf{k}_{n}},$

where $\mathbf{k}_{n}=k(X_{n},\,\mathbf{s})$. In the same manner, the optimal
solution of the PIMSE criterion (9) can be obtained by minimizing

$\displaystyle\text{Var}_{(m^{\ast}),n_{(m^{\ast})}+1}^{2}(\mathbf{s})=k(\mathbf{s})-\mathbf{v^{\ast}}^{\top}\mathbf{v}^{\ast},$
(14) $\displaystyle\mathbf{v}^{\ast}=L_{n+1}^{-1}k(X_{n+1},\,\mathbf{s}),$
(15)

Since we already have the solution of (15) partially with $\mathbf{v}$ (i.e.,
$\mathbf{v}^{\ast}[1:n]\equiv\mathbf{v}_{n}$), we need only
$v^{\ast}_{n+1}:=\mathbf{v}^{\ast}[n+1]$, which can be calculated with a
forward substitution as

$\displaystyle
v^{\ast}_{n+1}=k(\mathbf{x}^{\ast},\mathbf{s})-\sum_{j=1}^{n}L_{n+1,k}\mathbf{v}_{j}.$

If the importance $p(\mathbf{s})$ over $\Omega$ is trivial (e.g., uniform), we
can simplify (14). Since $k(\mathbf{s})$ in (14) is invariant to $\mathbf{s}$
and $\mathbf{x}^{\ast}$, we can only consider the second term
$\mathbf{v^{\ast}}^{\top}\mathbf{v}^{\ast}$ as the simplified PIMSE criterion,
which must be maximized in the PALC.

##  4 Simulation Study

In this section, we evaluate our active learning algorithm with simulation
data. Two functions are considered that can be visualized straightforwardly.
Both functions contain heterogeneous response surfaces, and observation noise
is involved. In order to reduce the variability from the different
initialization settings, each simulation is replicated ten times. As our
benchmark methods, the ALC and the ALM are considered for SGPs, and the
following partitioned active learning methods are considered: (i) the
variance-based criterion for the local searching referred to as PALM; (ii)
PALC without global searching (PALC-NoG). Additionally, uniform random
sampling (Rand) and LHD are also considered as passive learning. For
evaluation of models, one thousand space-filling design points for each design
space are used with root mean squared error (RMSE) as a metric. Computational
times spent on considered methods are also compared excluding the labelling
and model training times.

###  4.1 _One-dimensional Data_

We apply our proposed active learning algorithm to a one-dimensional
simulation function

$\displaystyle f(x)=2x\sin{(8\pi x^{3})},$

which is defined on [0, 1] as the dotted line in Fig. 2. The observation noise
is imposed with zero-mean Gaussian noise whose variance is
$\sigma^{2}=10^{-4}$. We allocate ten samples for initial training using
Maximin LHD, and 20 samples are sequentially obtained via active learning. The
function is differentiable, with heterogeneous frequency and amplitude over
the domain. Although there is no specific boundary for this function, we
rigorously partition the design space with the following translated heaviside
function, $g(x)=\mathbf{1}_{\\{x\geq 0.5\\}}$, which is assumed as the ground
truth.

Figure 2: Results of different active learning algorithms in 1-D simulation,
Left: SGP, Right: PGP with two local GPs

Fig. 2 shows each GP model fitted with initial samples. First, the right
figure in Fig. 2 illustrates that the PGP prevents misled active learning by
providing appropriate predictive uncertainty. Also, we can observe from the
left of Fig 2 that the next design point to be queried in the SGP model is
chosen from the low frequency region with both variance-based and IMSE
criteria. Meanwhile, in the PGP model, the variance-based criterion takes the
point from the boundary, while the IMSE does not. It implies that the IMSE
criterion can be more promising when the adjacent local GPs show comparable
predictive uncertainties.

After full running of each active learning strategy, the learning curves of
considered models are presented in the left plot of Fig. 4. From the result,
we can observe that both partitioned active learning strategies outperform the
SGPs, and the PALC is slightly better than the PALM. Considering the first
query location of each method in Fig. 2, it is not so surprising that the PALC
performs the best. Table 1 summarizes the results after the full data
acquisition, and we grouped considered methods into three: (i) passive
learning (Rand; and LHD); (ii) variance-based methods (ALM; and PALM); and
(iii) IMSE-based methods (ALC; PALC-NoG; and PALC). We can observe that the
predictive accuracy of PALC outperforms the others. In the computational time,
variance-based methods (ALM and PALM) are definitely faster than the IMSE-
based methods, while they are deficient in predictive accuracy; the ALM even
worse than the random sampling and the LHD. Among the IMSE-based methods, the
PALC is faster than the others. Moreover, if we focus on the PALC and PALC-
NoG, we can observe that the global searching does not only reduce the
computational time, but also improves the learning efficiency.

Table 1: One-dimensional simulation study results Methods |  RMSE |  Time (s)
---|---|---
Rand |  0.395 (0.152) |  -
LHD |  0.332 (0.133) |  -
ALM |  0.352 (0.159) |  0.341 (0.140)
PALM |  0.072 (0.270) |  0.144 (0.049)
ALC |  0.272 (0.205) |  11.468 (0.424)
PALC-NoG |  0.051 (0.277) |  7.049 (0.024)
PALC |  0.048 (0.273) |  3.553 (0.137)

###  4.2 _Two-dimensional Data_

We expand our active learning experiment on a two-dimensional function (shown
in the left of Fig. 3), which is also used in (Gramacy and Lee, 2009; Konomi
et al., 2014). Since the purpose of the function is to evaluate the
performance of the proposed active learning approach, it is comprised of two
regions: even and uneven, and zero-mean Gaussian noise with variance
$\sigma^{2}=10^{-6}$ is imposed as the observation noise. In a similar manner,
we begin with 15 samples with LHD, and obtain 15 additional samples via active
learning.

Figure 3: IMSE criterion contour plots in 2-D simulation, Left: ground truth
and initial design points with partitioned regions, Center: IMSE of ALC,
Right: PIMSE of PALC

For partitioning the design space, we borrow some ideas from Heaton et al.
(2017), exploiting the finite differences between initial samples. Unlikely to
their original clustering approach, we use a support vector classifier (SVC)
after labeling initial samples based on their magnitudes of finite differences
to generate a more flexible boundary than the Voronoi tessellation.
Consequently, we obtain a partitioning shown in the left figure of Fig. 3.
Since the region 1 is less interesting than the region 2, we can observe that
the PIMSE criterion induced by two independent local GPs provides more
relevant information of design points as shown in the right plot of Fig. 3.
Consequently, the IMSE criterion with the SGP fails to pick from the more
interesting region due to the misled information criterion.

In Table 2, we can see that the PALC surpasses the other methods again in both
predictive accuracy and computational time among the IMSE-based methods. Also,
the learning curve of PALC in the right plot of Fig. 4 shows that the PALC is
more promising than the active learning with SGPs.

Table 2: Two-dimensional simulation study results Methods |  RMSE |  Time (s)
---|---|---
Rand |  0.044 (0.017) |  -
LHD |  0.046 (0.021) |  -
ALM |  0.060 (0.023) |  0.344 (0.081)
PALM |  0.019 (0.025) |  0.151 (0.026)
ALC |  0.039 (0.013) |  8.368 (0.161)
PALC-NoG |  0.026 (0.004) |  5.570 (0.593)
PALC |  0.016 (0.005) |  1.924 (0.016)

Figure 4: Learning curves in simulation study, Left: one-dimensional, Right:
two-dimensional

##  5 Case Study

In this section, we apply our approach to construction of surrogate models of
two real-world cases. The original purpose of surrogate models herein is to
embed them into the automatic shape control system of fuselage assembly
process, and the optimal design of aluminum alloy. Both cases include higher
input dimensions than the previous simulations, and heterogeneity in their
design spaces due to complex material properties.

###  5.1 _Residual Stress of Composite Fuselages in Aerospace Manufacturing_

We apply our proposed active learning strategy to the construction of
predictive model of residual stress in the composite fuselage assembly
process. In the aircraft manufacturing process, composite fuselages are built
in several subsections independently, so they are subject to the discrepancy
in the junction part. In consequence, the reshaping of composite fuselage is
conducted with multiple fixed actuators with the automatic shape control. In
order to achieve the optimal manufacturing process, the shape control needs to
consider not only the deformation, but also the residual stress of the
structure due to their fatal affects on the final product. The development of
highly accurate predictive model for the shape control is very challenging,
since the problem is endowed with both the heterogeneity and the demanding
cost of real experiments and physics-based model. Especially, the stress of
composite fuselage is more difficult than deformation to predict (Lee et al.,
2020), we apply our method and other benchmarks to predictive modeling of the
stress.

In order to implement our case study cost-efficiently, we utilized the FEM
model, which is well-calibrated based on the real experiment (Wang et al.,
2020). The simulation mimics the real shape adjustment process that has ten
actuators under the fuselage section as shown in the left figure of Fig. 5
(Wen et al., 2018), and the maximum magnitude of actuator’s force is 450 lbf.
The maximum residual stress on the fuselage section is our interest, which is
measured in the psi scale. The ten-dimensional design space is partitioned
into three regions using SVC based on the clustering of input and output
pairs, and the ARD Matérn kernel is used for each local GPs and SGPs. As the
initial dataset, 20 samples are given with the LHD, and additional 20 samples
are sequentially obtained from different one thousand of different LHD points
with different active learning strategies. The model evaluation is conducted
with a separated testing dataset composed of 100 LHD samples, and mean
absolute error (MAE) is used as a metric.

Figure 5: Shape adjustment of composite fuselage, Left: composite fuselage
installed upon the fixture with actuators, Right: simulated residual stress of
composite fuselage in ANSYS

Table 3 summarizes the results of each learning strategy in the case study.
The PALC surpasses the other methods in both predictive error and the
computational time among the IMSE criterion-based methods. Interestingly,
except for the PALC-NoG, the other active learning methods are worse than two
passive learning strategies. Moreover, the PALM performs the worst in this
case, even though they chose design points to spread out the design space
well. In computational times, we can observe that the PALC-NoG has taken more
time than the ALC. This is due to the domination of one subregion in this
case. More explicitly, the PALC and PALC-NoG have assigned 15 and 16 design
points of 20 additional samples from one region, thereby decreasing the
efficiency of our matrix inversion scheme. However, it also shows that the
global searching of PALC scheme has reduced not only the computational time,
but also the predictive error.

Table 3: Residual stress case study results Methods |  MAE |  Time (s)
---|---|---
Rand |  3.858 |  -
LHD |  3.319 |  -
ALM |  4.545 |  $2.6\times 10^{-4}$
PALM |  11.555 |  $2.5\times 10^{-4}$
ALC |  4.693 |  11.596
PALC-NoG |  3.691 |  12.750
PALC |  3.207 |  9.729

###  5.2 _Tribocorrosion in Aluminium Alloys_

As our second case study, the material loss rate during stress corrosion (i.e.
tribocorrosion) in aluminum alloys with six control variables are considered.
To test the tribocorrosion resistance of metals, experimental tests and FEM
simulations were carried out by scratching the surface of the samples in
corrosive environment (Wang et al., 2021; Wang and Cai, 2021) as shown in Fig.
6. During the tribocorrosion process, the mechanical deformation and the
electrochemical processes including active corrosion and passivation work
synergistically to cause material degradation. The FEM model calculates the
contact mechanics between the indenter and the sample and simulate the wear
process as well as the wear-accelerated material dissolution of the corrosion
process, and generate the volume loss results. The 6 control variables for the
FEM model are material property descriptors: young’s modulus, yield strength,
anodic Tafel slope, anodic exchange current density, cathodic Tafel slope, and
cathodic exchange current density. The former two govern the mechanical
properties while the latter four determine the corrosion behavior of the
alloy. The output of the FEM model is the tribocorrosion rate of the alloy,
expressed as volume loss per time.

Figure 6: Schematic tribocorrosion simulation setup.

The surrogate model of the FEM model is constructed to assist the optimal
design of alloys with uncertainty quantification and alleviating the high-
computational cost of the FEM model. To establish the relationship between
material property and tribocorrosion rate, a total of 106 FEM simulations were
performed by systematically varying the 6 control variables. Since scales of
variables in the dataset are inconsistent, each variable is normalized to be
within a unit interval. For evaluation, mean relative error (MRE) is used as a
metric due to the infinitesimal scale of the output. The MRE is calculated as

$\displaystyle\text{MRE}=\frac{1}{n}\sum_{i=1}^{n}\frac{|y_{i}-h(x_{i})|}{|y_{i}|}.$

The PGP in this case is composed of 3 local GPs with the ARD RBF kernel, and
the SVC model is used for partitioning based on the normalized input-output
pairs. Considering the relatively small size of samples, 5-fold cross-
validation is used. That is, about 84 samples are passed to each model as the
candidate set, and the rest of samples are used for the model evaluation.
Compared methods are trained up to 50 samples from 20 common initial samples.

Table 4 shows the results of each learning strategy in which standard
deviations are obtained from the cross-validated evaluation. We can observe
that the PALC achieves the minimum averaged predictive error and the
computational time among the IMSE-based methods. Passive learning methods are
worse than others, and the variance-based methods are also worse than the
IMSE-based methods. In this case, overall computational times are much lower
than the previous case due to the small number of samples in the candidate
pool, while the numerical remedies in PALC have significantly reduced the time
of ALC.

Table 4: Tribocorrosion case study results Methods |  MRE |  Time (s)
---|---|---
Rand |  0.028 (0.017) |  -
LHD |  0.028 (0.017) |  -
ALM |  0.026 (0.017) |  0.014 (0.003)
PALM |  0.023 (0.012) |  0.104 (0.000)
ALC |  0.022 (0.013) |  1.283 (0.127)
PALC-NoG |  0.022 (0.012) |  0.764 (0.018)
PALC |  0.020 (0.013) |  0.757 (0.017)

##  6 Conclusion

Active learning is a special machine learning that seeks to improve sampling
efficiency and lower data collection cost. Existing active learning strategies
mainly focus on investigating one response surface, so they are insufficient
for reliable and cost-efficient surrogate modeling of heterogeneous systems.
This paper dedicated establishing an efficient partitioned active leaning
strategy that adopts two-step searching schemes based on the partitioned IMSE
criterion structure. By partitioning the design space into multiple subregions
according to heterogeneity in the target system, the global searching scheme
refers to the integrated predictive uncertainties of local GPs to determine
the most uncertain subregion. The global searching scheme allows us to reduce
the region of interest, thereby not only accelerating the searching speed, but
also improving the overall learning efficiency as shown in the simulation and
case studies. The local searching scheme exploited the chosen local GPs in the
global searching phase, so the localized IMSE criterion may provide more
relevant information minimizing the interruption of heterogeneous
characteristics in the other regions.

For the numerical perspective of active learning, the following applicable
remedies are provided: reducing the number of candidates with the global
searching scheme, and the Cholesky factor update, which can be embedded into
the PALC. From two of each simulation and case study, the proposed method
outperformed the benchmark methods including passive learning, the variance-
based methods, and the IMSE-based methods. Moreover, the global searching
scheme dramatically improved the performance of PALC by comparing our method
to the PIMSE without the global searching step. The results of simulation and
case studies imply that it is also applicable to other domains where
heterogeneity exists.

## Acknowledgements

The authors acknowledge the financial support by the US National Science
Foundation CMMI-2035038 and CMMI-1855651.

## References

  * Alaeddini et al. (2019) Alaeddini, A., E. Craft, R. Meka, and S. Martinez (2019).  Sequential laplacian regularized v-optimal design of experiments for response surface modeling of expensive tests: An application in wind tunnel testing.  IISE Transactions 51(5), 559–576.
  * Ankenman et al. (2010) Ankenman, B. E., B. L. Nelson, and J. Staum (2010).  Stochastic kriging for simulation metamodeling.  Operations Research 58(2), 371–382.
  * Cao et al. (2020) Cao, X., J. Yao, Z. Xu, and D. Meng (2020).  Hyperspectral image classification with convolutional neural network and active learning.  IEEE Transactions on Geoscience and Remote Sensing 58(7), 4604–4616.
  * Chen and Zhou (2017) Chen, X. and Q. Zhou (2017).  Sequential design strategies for mean response surface metamodeling via stochastic kriging with adaptive exploration and exploitation.  European Journal of Operational Research 262(2), 575–585.
  * Cho et al. (2019) Cho, W., Y. Kim, and J. Park (2019).  Hierarchical anomaly detection using a multioutput gaussian process.  IEEE Transactions on Automation Science and Engineering 17(1), 261–272.
  * Cohn et al. (1996) Cohn, D. A., Z. Ghahramani, and M. I. Jordan (1996).  Active learning with statistical models.  Journal of Artificial Intelligence Research 4, 129–145.
  * Cortes et al. (2019) Cortes, C., G. DeSalvo, C. Gentile, M. Mohri, and N. Zhang (2019).  Region-based active learning.  In The 22nd International Conference on Artificial Intelligence and Statistics, pp. 2801–2809. PMLR.
  * Deisenroth et al. (2013) Deisenroth, M. P., D. Fox, and C. E. Rasmussen (2013).  Gaussian processes for data-efficient learning in robotics and control.  IEEE Transactions on Pattern Analysis and Machine Intelligence 37(2), 408–423.
  * Erickson et al. (2018) Erickson, C. B., B. E. Ankenman, M. Plumlee, and S. M. Sanchez (2018).  Gradient based criteria for sequential design.  In 2018 Winter Simulation Conference (WSC), pp. 467–478. IEEE.
  * Freytag et al. (2013) Freytag, A., E. Rodner, P. Bodesheim, and J. Denzler (2013).  Labeling examples that matter: Relevance-based active learning with gaussian processes.  In German Conference on Pattern Recognition, pp. 282–291. Springer.
  * Gramacy and Lee (2009) Gramacy, R. B. and H. K. Lee (2009).  Adaptive design and analysis of supercomputer experiments.  Technometrics 51(2), 130–145.
  * Gramacy and Lee (2008) Gramacy, R. B. and H. K. H. Lee (2008).  Bayesian treed gaussian process models with an application to computer modeling.  Journal of the American Statistical Association 103(483), 1119–1130.
  * Heaton et al. (2017) Heaton, M. J., W. F. Christensen, and M. A. Terres (2017).  Nonstationary gaussian process models using spatial hierarchical clustering from finite differences.  Technometrics 59(1), 93–101.
  * Heinonen et al. (2016) Heinonen, M., H. Mannerström, J. Rousu, S. Kaski, and H. Lähdesmäki (2016).  Non-stationary gaussian process regression with hamiltonian monte carlo.  In Artificial Intelligence and Statistics, pp. 732–740. PMLR.
  * Higdon (1998) Higdon, D. (1998).  A process-convolution approach to modelling temperatures in the north atlantic ocean.  Environmental and Ecological Statistics 5(2), 173–190.
  * Hyer and White (2009) Hyer, M. W. and S. R. White (2009).  Stress analysis of fiber-reinforced composite materials.  DEStech Publications, Inc.
  * Käding et al. (2018) Käding, C., E. Rodner, A. Freytag, O. Mothes, B. Barz, J. Denzler, and C. Z. AG (2018).  Active learning for regression tasks with expected model output changes.  In BMVC, pp. 103.
  * Kim et al. (2009) Kim, B., Y. Lee, and D.-H. Choi (2009).  Construction of the radial basis function based on a sequential sampling approach using cross-validation.  Journal of Mechanical Science and Technology 23(12), 3357–3365.
  * Kim et al. (2005) Kim, H.-M., B. K. Mallick, and C. Holmes (2005).  Analyzing nonstationary spatial data using piecewise gaussian processes.  Journal of the American Statistical Association 100(470), 653–668.
  * Konomi et al. (2014) Konomi, B., G. Karagiannis, A. Sarkar, X. Sun, and G. Lin (2014).  Bayesian treed multivariate gaussian process with adaptive design: Application to a carbon capture unit.  Technometrics 56(2), 145–158.
  * Lee et al. (2020) Lee, C., J. Wu, W. Wang, and X. Yue (2020).  Neural network gaussian process considering input uncertainty for composite structures assembly.  IEEE/ASME Transactions on Mechatronics early access(26 Nov 2020), 1–1.
  * Liu et al. (2018) Liu, H., Y.-S. Ong, and J. Cai (2018).  A survey of adaptive sampling for global metamodeling in support of simulation-based complex engineering design.  Structural and Multidisciplinary Optimization 57(1), 393–416.
  * MacKay (1992) MacKay, D. J. (1992).  Information-based objective functions for active data selection.  Neural Computation 4(4), 590–604.
  * Marmin et al. (2018) Marmin, S., D. Ginsbourger, J. Baccou, and J. Liandrat (2018).  Warped gaussian processes and derivative-based sequential designs for functions with heterogeneous variations.  SIAM/ASA Journal on Uncertainty Quantification 6(3), 991–1018.
  * Nguyen and Smeulders (2004) Nguyen, H. T. and A. Smeulders (2004).  Active learning using pre-clustering.  In Proceedings of the Twenty-first International Conference on Machine Learning, pp. 79.
  * Pasolli and Melgani (2011) Pasolli, E. and F. Melgani (2011).  Gaussian process regression within an active learning scheme.  In 2011 IEEE International Geoscience and Remote Sensing Symposium, pp. 3574–3577. IEEE.
  * Pope et al. (2021) Pope, C. A., J. P. Gosling, S. Barber, J. S. Johnson, T. Yamaguchi, G. Feingold, and P. G. Blackwell (2021).  Gaussian process modeling of heterogeneity and discontinuities using voronoi tessellations.  Technometrics 63(1), 53–63.
  * Rasmussen and Williams (2006) Rasmussen, C. and C. Williams (2006).  Gaussian Processes for Machine Learning.  Cambridge, MA, USA: MIT Press.
  * Rasmussen and Ghahramani (2002) Rasmussen, C. E. and Z. Ghahramani (2002).  Infinite mixtures of gaussian process experts.  Advances in Neural Information Processing Systems 2, 881–888.
  * Santner et al. (2018) Santner, T. J., B. J. Williams, and W. I. Notz (2018).  The design and analysis of computer experiments, Volume 2.  Springer.
  * Seo et al. (2000) Seo, S., M. Wallat, T. Graepel, and K. Obermayer (2000).  Gaussian process regression: Active data selection and test point rejection.  In Mustererkennung 2000, pp. 27–34. Springer.
  * Settles (2009) Settles, B. (2009).  Active learning literature survey.  Computer Sciences Technical Report 1648, University of Wisconsin–Madison.
  * Snoek et al. (2012) Snoek, J., H. Larochelle, and R. P. Adams (2012).  Practical bayesian optimization of machine learning algorithms.  In F. Pereira, C. J. C. Burges, L. Bottou, and K. Q. Weinberger (Eds.), Advances in Neural Information Processing Systems, Volume 25. Curran Associates, Inc.
  * Uteva et al. (2018) Uteva, E., R. S. Graham, R. D. Wilkinson, and R. J. Wheatley (2018).  Active learning in gaussian process interpolation of potential energy surfaces.  The Journal of Chemical Physics 149(17), 174114.
  * Wang and Cai (2021) Wang, K. and W. Cai (2021).  Modeling the effects of individual layer thickness and orientation on the tribocorrosion behavior of al/cu nanostructured metallic multilayers.  Wear 477, 203849.
  * Wang et al. (2021) Wang, K., Y. Wang, X. Yue, and W. Cai (2021).  Multiphysics modeling and uncertainty quantification of tribocorrosion in aluminum alloys.  Corrosion Science 178, 109095.
  * Wang et al. (2020) Wang, Y., X. Yue, R. Tuo, J. H. Hunt, J. Shi, et al. (2020).  Effective model calibration via sensible variable identification and adjustment with application to composite fuselage simulation.  Annals of Applied Statistics 14(4), 1759–1776.
  * Wen et al. (2018) Wen, Y., X. Yue, J. H. Hunt, and J. Shi (2018).  Feasibility analysis of composite fuselage shape control via finite element analysis.  Journal of Manufacturing Systems 46, 272–281.
  * Yue et al. (2018) Yue, X., Y. Wen, J. H. Hunt, and J. Shi (2018, 02).  Surrogate Model-Based Control Considering Uncertainties for Composite Fuselage Assembly.  Journal of Manufacturing Science and Engineering 140(4), MANU–17–1157.  041017\.
  * Yue et al. (2020) Yue, X., Y. Wen, J. H. Hunt, and J. Shi (2020).  Active learning for gaussian process considering uncertainties with application to shape control of composite fuselage.  IEEE Transactions on Automation Science and Engineering 18(1), 36–46.