Given multi-fidelity data $D^{i}=\{{\bf x}^{i}_{n},y^{i}_{n}\}_{n=1}^{N^{i}}$, for $i=1,\dots,\tau$, where ${\bf x}\in\mathbb{R}^{l}$ denotes the system inputs (a vector of parameters that appear in the system of equations and/or in the initial-boundary conditions for a simulation); ${\bf y}^{i}\in\mathbb{R}^{d^{i}}$ denotes the corresponding outputs, where $d^{i}$ is the dimension for $i$ fidelity; $\tau$ is the total number of fidelities. Generally speaking, higher-fidelity are closer to the ground truth and are expensive to obtain, and thus we have fewer samples for high fidelity, i.e., $N^{1}>N^{2}>\dots>N^{\tau}$. The dimensionality $d^{i}$ is not necessary the same or aligned across different fidelities. In most work e.g., [15, 16, 11], the system inputs of higher-fidelity are chosen to be the subset of the lower-fidelity, i.e., ${\bf X}^{\tau}\subset\dots\subset{\bf X}^{2}\subset{\bf X}^{1}$. We call this the subset structure for a multi-fidelity dataset, as opposed to arbitrary data structures, which we will resolve in Section section 3.3 with a closed-form solution and extend it to the classic AR. Our goal is to estimate the function ${\bf f}^{\tau}:\mathbb{R}^{l}\rightarrow\mathbb{R}^{d^{\tau}}$ given the observation data across different fidelities $\{D^{i}\}_{i=1}^{\tau}$.

For the sake of clarity, we consider a two-fidelity case with superscript $h$ indicating high-fidelity and $l$ low-fidelity. Nevertheless, the formulation can be generalized to problems with multiple fidelities recursively. Considering a simple scalar output for all fidelities, the AR [3] assumes

where $\rho$ is a factor transferring knowledge from the low fidelity in a linear fashion whereas $f^{r}({\bf x})$ tries to capture the residual information. If we assume a zero mean Gaussian process (GP) prior [17] (see Appendix A for a brief description) for $f^{l}({\bf x})$ and $f^{r}({\bf x})$, i.e., $f^{l}({\bf x})\sim{\mathcal{N}}(0,k^{l}({\bf x},{\bf x}^{\prime}))$ and $f^{r}({\bf x})\sim{\mathcal{N}}(0,k^{r}({\bf x},{\bf x}^{\prime}))$, the high-fidelity function also follows a GP. This gives an elegant joint GP for the joint observations ${\bf y}=[{{\bf y}^{l}};{{\bf y}^{h}}]^{T}$,

where ${\bf y}^{l}\in\mathbb{R}^{N_{l}\times 1}$ is the low-fidelity observations corresponding to input ${\bf X}^{l}\in\mathbb{R}^{N_{l}\times L}$ and ${\bf y}_{h}\in\mathbb{R}^{N_{h}\times 1}$ is the high-fidelity observations; $[{\bf K}^{l}({\bf X}^{l},{\bf X}^{l})]_{ij}=k^{l}({\bf x}_{i},{\bf x}_{j})$ is the covariance matrix of the low-fidelity inputs ${\bf x}_{i},{\bf x}_{j}\in{\bf X}^{l}$; $[{\bf K}^{r}({\bf X}^{l},{\bf X}^{l})]_{ij}=k^{r}({\bf x}_{i},{\bf x}_{j})$ is for the high-fidelity inputs ${\bf x}_{i},{\bf x}_{j}\in{\bf X}^{h}$; $[{\bf K}^{l}({\bf X}^{l},{\bf X}^{h})]_{ij}=k^{l}({\bf x}_{i},{\bf x}_{j})$ is the cross-fidelity covariance matrix of the low-fidelity inputs ${\bf x}_{i}\in{\bf X}^{l}$ and high-fidelity inputs ${\bf x}_{j}\in{\bf X}^{h}$, and ${\bf K}^{l}({\bf X}^{h},{\bf X}^{l})=({\bf K}^{l}({\bf X}^{l},{\bf X}^{h}))^{T}$. One immediate advantage of the AR is that the joint Gaussian form allows not only joint training for all low- and high-fidelity data but also predictions for any given new inputs by conditional Gaussian (as the posterior is derived in a standard GP [17]). Furthermore, Le Gratiet [15] derive Lemma 1 to reduce the complexity from $O((N^{l}+N^{h})^{3})$ to $O((N^{l})^{3}+(N^{h})^{3})$ with a subset data structure.

To resolve the scalability issue, we rearrange all the output into a multidimensional space (i.e., a tensor space) and introduce latent coordinate features to index the outputs to capture their correlations as in HOGP [18]. More specifically, we organize the low-fidelity outputs as a $M$-mode tensor, ${\bm{\mathsfit{Z}}}^{l}\in\mathbb{R}^{d^{l}_{1}\times\dots\times d^{l}_{M}}$, where the output dimension $d^{l}=\prod_{m=1}^{M}d^{l}_{m}$. The element ${\mathsfit{Z}}^{l}$ is indexed based on its coordinates ${\bf c}=(c_{1},\dots,c_{M})$($1\leqslant c_{k}\leqslant d_{k}$ for $k=1,\dots,M$). If the original data indeed admit a multi-array structure, we can use their original index with actual meaning to index the coordinates. For instance, a 2D spatial dataset can use its original spatial coordinate to index a single location (pixel). To improve our model flexibility, we do not have to limit ourselves from using the original index, particularly for the cases where the original data does not admit a multi-array structure or the multi-array structure is of too large size. In such case, we can use arbitrary tensorization and a latent feature vector ${\bf v}^{l}_{c_{m}}$ (whose values are inferred in model training) for each coordinate $c_{m}$ in mode $m$. This way, the element ${\mathsfit{Z}}^{l}$ is indexed by the vector $({\bf v}^{l}_{c_{1}},\dots,{\bf v}^{l}_{c_{M}})$. Following the linear transformation of Eq. (1), we first introduce the tensor-matrix product [19] at mode $m$,

where ${\bm{\mathsfit{F}}}^{h}({\bf x})$ denotes target function ${\bf f}^{h}({\bf x})$ with its output organized into multi-array ${\bm{\mathsfit{Z}}}^{h}$, and the same concept applies to ${\bm{\mathsfit{F}}}^{l}({\bf x})$ and ${\bm{\mathsfit{F}}}^{r}({\bf x})$; $\times_{m}$ denotes the tensor-matrix product at mode $m$. To give a concrete example, considering an arbitrary tensor ${\bm{\mathsfit{Z}}}^{l}\in\mathbb{R}^{d^{l}_{1}\times\dots\times d_{M}^{l}}$ and a matrix ${\bf W}_{m}\in\mathbb{R}^{s\times d_{m}}$, the $\times_{m}$ product is calculated as $[{\bm{\mathsfit{Z}}}^{l}\times_{m}{\bf W}_{m}]_{i_{1},\dots,i_{m-1},j,i_{m+1}\dots,i_{M}}=\sum_{k=1}^{c_{m}}w_{jk}{\mathsfit{Z}}_{i_{1},\dots,i_{m-1},k,i_{m+1}\dots,i_{M}}$, which becomes an $d_{1}^{l}\times\dots\times d_{m-1}^{l}\times s\times d_{m+1}^{l}\times\dots\times d^{l}_{M}$ tensor. We can further denote the group of $M$ linear transformation matrixes as a Tucker tensor ${\bm{\mathsfit{W}}}=[{\bf W}_{1},\dots,{\bf W}_{M}]$ and represent Eq. (3) compactly using a Tucker operator [19], ${\bm{\mathsfit{F}}}^{l}({\bf x})\times{\bm{\mathsfit{W}}}$, which has an important property:

Inspired by the AR of Eq. (1), we place a tensor-variate GP (TGP) prior [20] for the low-fidelity tensor function ${\bm{\mathsfit{F}}}^{l}({\bf x})$ and the residual tensor function ${\bm{\mathsfit{F}}}^{r}({\bf x})$:

where ${\bf S}_{m}^{i}\in\mathbb{R}^{d_{m}\times d_{m}}$ are the output correlation matrix with $[{\bf S}_{m}^{i}]_{jk}=\tilde{k}^{i}_{m}({\bf v}^{i}_{c_{i}},{\bf v}^{i}_{c_{k}})$ and $\tilde{k}^{i}_{m}(\cdot,\cdot)$ being the kernel function (with unknown hyperparameters). A TGP is a generalization of a multivariate GP that essentially represents a joint GP prior $\mathrm{vec}({{\bm{\mathsfit{Y}}}^{l}})\sim{\mathcal{N}}\left(0,{\bf K}^{l}({\bf X}^{l},{\bf X}^{l})\bigotimes_{m=1}^{M}{\bf S}_{m}\right)$. Similar to the joint probability of (2), we can derive the joint probability for ${\bf y}=[\mathrm{vec}({{\bm{\mathsfit{Y}}}^{l}})^{T},\mathrm{vec}({\bm{\mathsfit{Y}}}^{h})^{T}]^{T}$ based on Tucker transformation of (3); we keep the the proof in the Appendix for clarity.

Lemma 2 admits any arbitrary outputs (living in different spaces, having different dimensions and/or modes, and being unaligned) at different fidelity. Also, it does not require a subset dataset to hold.

Lemma 2 defines our GAR model, a generalized AR with special tensor structures. The covariance for low-fidelity is $\mathrm{cov}({\mathsfit{Z}}^{l}_{\bf c}({\bf x}),{\mathsfit{Z}}^{l}_{{\bf c}^{\prime}}({\bf x}^{\prime}))=k^{l}({\bf x},{\bf x}^{\prime})\prod_{m=1}^{M}\tilde{k}^{l}_{m}({\bf v}^{m}_{c_{m}},{\bf v}^{m}_{c^{\prime}_{m}})$, cross-fidelity $\mathrm{cov}({\mathsfit{Z}}^{l}_{\bf c}({\bf x}),{\mathsfit{Z}}^{h}_{{\bf c}^{\prime}}({\bf x}^{\prime}))=k^{l}({\bf x},{\bf x}^{\prime})\prod_{m=1}^{M}\tilde{k}^{l}_{m}({\bf v}^{l}_{c_{m}},{\bf v}^{l}_{c^{\prime}_{m}})w^{m}_{c,c^{\prime}}$ (where $w^{m}_{c,c^{\prime}}$ is the ${c,c^{\prime}}$-th element of ${\bf W}^{m}$), and high-fidelity $\mathrm{cov}({\mathsfit{Z}}^{h}_{\bf c}({\bf x}),{\mathsfit{Z}}^{h}_{{\bf c}^{\prime}}({\bf x}^{\prime}))=k^{l}({\bf x},{\bf x}^{\prime})\prod_{m=1}^{M}\tilde{k}^{l}_{m}({\bf v}^{l}_{c_{m}},{\bf v}^{l}_{c^{\prime}_{m}})(w^{m}_{c,c^{\prime}})^{2}+k^{h}({\bf x},{\bf x}^{\prime})\prod_{m=1}^{M}\tilde{k}^{r}_{m}({\bf u}^{m}_{c_{m}},{\bf u}^{m}_{c^{\prime}_{m}})(w^{m}_{c,c^{\prime}})^{2}$. The complex between-fidelity output correlations are captured using latent features $\{{\bf V}^{m},{\bf U}^{m}\}_{m=1}^{M}$ with arbitrary kernel function $\tilde{k}_{m}^{i}$ whereas the cross-fidelity output correlations are captured in a composite manner. This combination overcomes the simple linear correlations assumed in previous work that simply decomposes the output as a dimension reduction preprocess [12]. When the dimensionality aligns for ${\bm{\mathsfit{Z}}}^{l}$ and ${\bm{\mathsfit{Z}}}^{h}$ and thus $d^{l}_{m}=d^{h}_{m}$, we can share the same latent features across the two fidelity data by letting ${\bf v}^{m}_{j}={\bf u}^{m}_{j}$ while keeping the kernel functions different. This way, the latent features are more resistant to overfitting. For non-aligned data with explicit indexing, we can use kernel interpolation [21] for the same purpose. To further encourage sparsity in the latent feature, we impose a Laplace prior, i.e., ${\bf v}^{m}_{j}\sim\mathrm{Laplace}(\lambda)\propto\exp(-\lambda||{\bf v}^{m}_{j}||_{1})$.

With the model fully defined, we can now train the model to obtain all unknown model parameters. For compactness, we use the following compact notation ${\bf S}^{l}=\bigotimes_{m=1}^{M}{\bf S}^{l}_{m}$, ${\bf S}^{h}=\bigotimes_{m=1}^{M}{\bf S}^{h}_{m}$, ${\bf W}=\bigotimes_{m=1}^{M}{\bf W}_{m}$, ${\bf K}^{l}={\bf K}^{l}({\bf X}^{l},{\bf X}^{l})$, ${\bf K}^{lh}={\bf K}^{l}({\bf X}^{l},{\bf X}^{h})$, ${\bf K}^{hl}={\bf K}^{l}({\bf X}^{h},{\bf X}^{l})$, ${\bf K}^{lr}={\bf K}^{l}({\bf X}^{h},{\bf X}^{h})$, and, ${\bf K}^{r}={\bf K}^{r}({\bf X}^{h},{\bf X}^{h})$ (with a slight abuse of notation).

We preserve the proof in the Appendix for clarity. Note that $\mathcal{L}^{l}$ and $\mathcal{L}^{r}$ are HOGP likelihoods for ${\bm{\mathsfit{Y}}}^{l}$ and residual ${\bm{\mathsfit{Y}}}^{h}-{\bm{\mathsfit{Y}}}^{l}\times\hat{{\bm{\mathsfit{W}}}}$, respectively. Since the computational of $\mathcal{L}^{l}$ and $\mathcal{L}^{r}$ are independent, the model training can be conducted efficiently in parallel.

Predictive posterior. Similarly, we can derive the concrete predictive posterior for high-fidelity outputs by integrating out the latent functions after some tedious linear algebra (see Appendix), which is also Gaussian, $\mathrm{vec}({\bm{\mathsfit{Z}}}^{h}_{*})\sim\mathcal{N}(\mathrm{vec}(\bar{{\bm{\mathsfit{Z}}}}^{h}_{*}),{{\bf S}^{h}_{*}})$, where

where ${\bf k}^{l}_{*}={\bf k}^{l}({\bf x}_{*},{\bf X}^{l})$ is the vector of covariance between the give input ${\bf x}_{*}$ and low-fidelity observation inputs ${\bf X}^{l}$ and similarly, ${\bf k}^{l}_{**}={\bf k}^{l}({\bf x}_{*},{\bf x}_{*})$, ${\bf k}^{r}_{*}={\bf k}^{r}({\bf x}_{*},{\bf X}^{h})$, ${\bf k}^{r}_{**}={\bf k}^{r}({\bf x}_{*},{\bf x}_{*})$.

In practice, it is sometimes difficult to ask the multi-fidelity data to preserve a subset structure, particularly in the case of multi-fidelity Bayesian optimization [22, 23]. This presents the challenge for most SOTA multi-fidelity models e.g., NAR [16], ResGP [9], stochastic collocation [24]. In contrast, the advantage of AR is that even if the multi-fidelity data does not admit a subset data structure, the model can still be trained using all available data based on the joint likelihood of (2). However, this method lacks scalability due to the inversion of the large joint covariance matrix $\bm{\Sigma}$. The situation will get worse if we deal with multi-fidelity data with more than two fidelities. To resolve this issue, we propose a fast inference method based on imaginary subset. More specifically, considering the missing low-fidelity data as latent variables $\hat{{\bm{\mathsfit{Y}}}}^{l}$, the joint likelihood function is

where $p({\bm{\mathsfit{Y}}}^{h}|\hat{{\bm{\mathsfit{Y}}}}^{l},{\bm{\mathsfit{Y}}}^{l})$ is the likelihood in Lemma 3 given the complementary imaginary subset, and $p(\hat{{\bm{\mathsfit{Y}}}}^{l}|{\bm{\mathsfit{Y}}}^{l})\sim{\mathcal{N}}(\bar{{\bm{\mathsfit{Y}}}}^{l},\hat{{\bf S}}^{l}\otimes{\bf S}^{l})$ is the imaginary posterior with the given low-fidelity observations ${\bm{\mathsfit{Y}}}^{l}$. The integral can be calculated using Gaussian quadrature or other sampling methods as in [8, 25], which is slow and inaccurate.

We preserve the proof in the appendix. Notice that $\hat{{\bf E}}\hat{{\bf S}}^{l}\hat{{\bf E}}^{T}\otimes{\bf W}^{T}{\bf S}^{l}{{\bf W}}=\left(\begin{array}[]{cc}{\bf 0}&{\bf 0}\\ {\bf 0}&\hat{{\bf S}}^{l}\end{array}\right)\otimes{\bf W}^{T}{\bf S}^{l}{{\bf W}}$ is the low-right block of the predictive variance for the missing low-fidelity observations $\mathrm{vec}(\hat{{\bm{\mathsfit{Y}}}})$; We can easily understand the last part of the likelihood as a GP with accumulated uncertainty(variance) added to the corresponding missing points. Lemma 4 naturally applied to AR when the output is a scalar, where ${\bf W}=\rho$, ${\bf S}^{l}=1$, and ${\bf S}^{r}=1$

Predictive posterior. Surprisingly, the posterior also turns out to be a Gaussian distribution,

where and $\mathbf{\Gamma}$ and the mean of the predictive posterior $\bar{{\bm{\mathsfit{Z}}}}_{*}$ are given as follows,

where ${\bf E}_{m}$ and ${\bf E}_{n}$ are the selection matrices that $\hat{{\bf X}}^{h}={\bf E}_{m}^{T}[{\bf X}^{l},\hat{{\bf X}}^{h}]$, ${\bf X}^{h}={\bf E}_{n}^{T}[{\bf X}^{l},\hat{{\bf X}}^{h}]$, $\hat{{\bf W}}={\bf E}_{n}^{T}\otimes{\bf W}$ and $\hat{{\bf K}}^{l}$ is the covariance matrix that $\hat{{\bf K}}^{l}={\bf K}^{l}([{\bf X}^{l},\hat{{\bf X}}^{h}],[{\bf X}^{l},\hat{{\bf X}}^{h}])$.

For subset structured data, the computational complexity of GAR is decomposed into two independent TGPs for likelihood and predictive posterior. Thanks to the tensor algebra (mainly $({\bf K}\otimes{\bf S})^{-1}={\bf K}^{-1}\otimes{\bf S}^{-1}$), the complexity of the $i$-fidelity kernel matrix inversion is reduced to $\mathcal{O}(\sum_{m=1}^{M}(d^{i}_{m})^{3}+(N^{i})^{3})$ instead of $\mathcal{O}((N^{i}d^{i})^{3})$. For the non-subset case, the computational complexity in Eq. (8) is unfortunately $\mathcal{O}((N_{m}^{i}d^{i})^{3})$ where $N_{m}$ is the number of the imaginary low-fidelity points. Nevertheless, due to the tensor structure, we can still use conjugate gradient [26] to solve the linear system efficiently.

Notice that the mean prediction $\bar{{\bm{\mathsfit{Z}}}}^{h}_{*}$ in Eq. (A.60) and Eq. (9) does not depend on any output covariance matrixes $\{{\bf S}^{h}_{m},{\bf S}^{l}_{m}\}_{m=1}^{M}$, which reassemble the autokrigeability (no knowledge transfer in noiseless cases for mean predictions) [14, 9] based on the GAR framework. For applications where the predictive variation is not of interest, we can introduce a conditional independent output-correlation, i.e., ${\bf S}^{h}_{m}={\bf I},{\bf S}^{l}_{m}={\bf I}$ and orthogonal weight matrices, i.e., ${\bf W}_{m}^{T}{\bf W}_{m}={\bf I}$, to reduce the computationally complexity further down to $\mathcal{O}((N^{i})^{3})$ (see Appendix for detailed proof). We call this CIGAR, an abbreviation for conditional independent GAR. In practice, CIGAR is slightly worse than GAR due to the difficulty of ensuring ${\bf W}_{m}^{T}{\bf W}_{m}={\bf I}$.

We first access GAR in canonical PDE simulation benchmarks, which produce high-dimensional spatial/spatial-temporal fields as model outputs. Specifically, we test on Burgers’ , Poisson’s and Heat equations as in [12, 51, 52, 53]. The high fidelity results are obtained by solving these equations using finite difference on a $32\times 32$ mesh, whereas low fidelity by an $8\times 8$ coarse mesh. The solutions on these grid points are recorded and vectorized as outputs. Because the mesh differs, the dimensionality varies. To compare with the standard multi-fidelity method that can only deal with aligned outputs, we use interpolation to upscale the low-fidelity and record them at the high-fidelity grid nodes. The corresponding inputs are PDE parameters and parameterized initial or boundary conditions. Detailed experimental setups can be found in Appendix E.1

We uniformly generate $128$ samples for testing and 32 for training. We increase the high-fidelity training samples to the number of low-fidelity training samples of 32. The comparisons are conducted five times with shuffled samples. The statistical results (mean and std) of the RMSE are reported in Fig. 1. GAR and CIGAR outperform competitors with a significant margin, up to 6x reductions in RMSE and also reaching the optimal with a maximum of 8 high-fidelity samples, indicating a successful transfer from low- to high-fidelity. CIGAR is slightly worse than GAR possibility due to a lack of hard constraint on the orthogonality of its weight matrixes during implementation. As we have discovered in the literature, AR consistently performs well. With a flexible linear transformation, GAR outperforms the AR while inheriting its robustness, leading to the best performance. For the unaligned output, MF-BNN showed slightly worse performance than in the aligned cases, highlighting the challenges of the unaligned outputs. In contrast, GAR and CIGAR show almost identical performance for both cases. Nevertheless, MF-BNN also shows good performance compared to the rest of the other methods, which is consistent with the finding in [13]. It is interesting to see that for the non-subset data, the capable methods show better performances than in the subset cases. GAR and CIGAR still outperform the competitors with a clear margin.

To approximately assess the performance under an active learning process. We instead generate training samples in a Sobol sequence [54]. The results are shown in Appendix E.2 where GAR and CIGAR also outperform the other methods by a large margin.

Optimal topology structure is the optimized layout of materials, e.g., alloy and concrete, given some design specifications, e.g., external force and angle. This topology optimization is a key technique in mechanical designs, but it is also known for its high computational cost, which renders the need for multi-fidelity fusion. We consider the topology optimization of a cantilever beam with the location of point load, the angle of point load, and the filter radius [55] as system inputs. The low-fidelity use a $16\times 16$ regular mesh for the finite element solver, whereas the high-fidelity $64\times 64$. Please see Appendix E.3 for detailed setup.

As in the previous experiment, the RMSE statistics against an increasing number of high-fidelity training samples are shown in Fig. 2. It is clear that GAR outperforms other competitors with a large margin consistently. CIGAR can be as good as GAR when the number of training samples is large.

Steady-state 3D solid oxide fuel cell which solves complex coupled PDEs including Ohm’s law, Navier-Stokes equations, Brinkman equation, and Maxwell-Stefan diffusion and convection simultaneously, is a key model for modern fuel cell optimization. The model was solved using finite element in COMSOL. The inputs were taken to be the electrode porosities, cell voltage, temperature, and pressure in the channels. The low-fidelity experiment is conducted using 3164 elements and relative tolerance of 0.1, whereas the high-fidelity uses 37064 elements and relative tolerance of 0.001. The outputs are the coupled fields of electrolyte current density (ECD) and ionic potential (IP) in the $x-z$ plane located at the center of the channel. Please see Appendix E.4 for detailed experiment setup.

The RMSE statistics are shown in Fig. 3(a), which, again, highlights the superiority of the proposed method with only four high-fidelity data points. To further assess the model capacity for non-structured outputs, we keep only the ECD (Fig. 3(b)) and IP (Fig. 3(c)) in the low-fidelity training data to rise the challenges of predicting high fidelity ECD+IP fields. We can see that removing some low-fidelity information indeed increases the difficulties, especially when removing ECD, where MF-BNN outperforms GAR and CIGAR with a small number of training data. As soon as the training number increases, GAR and CIGAR become superior again.

Plasmonic nanoparticle arrays is a complex physical simulation that calculate the extinction and scattering efficiencies $Q_{ext}$ and $Q_{sc}$ for plasmonic systems with varying numbers of scatterers using the Coupled Dipole Approximation (CDA) approach. CDA is a method for mimicking the optical response of an array of similar, non-magnetic metallic nanoparticles with dimensions far smaller than the wavelength of light (here 25 nm). $Q_{ext}$ and $Q_{sc}$ are defined as the QoIs in this work. Please see Appendix E.5 for detailed experiment setup.

We conducted the experiments 5 times with shuffled samples, and we fixed the number of low-fidelity training samples to 32, 64, and 128 and gradually increase the high-fidelity training data from 4 to 32, 64, and 128. We can see in Fig.4, GAR outperforms others by a clear margin, especially when the high-fidelity data is only 4 samples. When there is a large training sample dataset, CIGAR can be as excellent as GAR.