1 Introduction↩︎

Uncertainty is pervasive in engineering practice, arising from manufacturing imperfections, model inadequacy, and material variability. Accurate quantification of such uncertainties is fundamental to reliability assessment and robust design [1]–[3]. In general, uncertainty is commonly categorized into two fundamentally different types: aleatory and epistemic uncertainty [4], [5]. Aleatory uncertainty refers to inherent randomness in physical systems that cannot be eliminated and is typically represented using well-defined probability distributions. In contrast, epistemic uncertainty arises from incomplete knowledge or limited data and is, in principle, reducible as information increases [6]. Maintaining a clear distinction between these two sources of uncertainty is essential for rigorous uncertainty quantification [7], since conflating them may lead to overconfident probabilistic representations and misleading reliability estimates [8], [9].

In practical engineering systems, aleatory and epistemic uncertainties coexist, motivating the development of hybrid uncertainty frameworks [10]. To represent epistemic uncertainty without imposing overly restrictive probabilistic assumptions, several mathematical frameworks have been proposed, including interval analysis, evidence theory, and probability boxes (p-boxes) [11]–[16]. These approaches replace a single probability distribution with a family of admissible distributions, thereby enabling bounded probabilistic descriptions under incomplete information while preserving stochastic modeling of aleatory variability.

Given such hybrid uncertainty representations, Global Sensitivity Analysis (GSA) provides a systematic framework to quantify how input uncertainties propagate to model outputs [17]. Variance-based methods, particularly Sobol’ indices, are widely used due to their rigorous decomposition of output variance into orthogonal contributions of inputs and their interactions [18]–[20]. However, classical GSA is fundamentally designed for scalar-valued outputs under fully probabilistic settings.

When hybrid uncertainty is present, two fundamental difficulties arise. First, the output is no longer characterized by a single probability distribution or scalar variance, but by a family of distributions or bounded response sets [21]. Second, and more critically, in many modern engineering applications the output is vector-valued or spatially distributed, and existing sensitivity measures often collapse such responses into a single scalar index. This aggregation leads to a loss of structural information and prevents the identification of component-wise or mode-wise sensitivity patterns.

To address these limitations, several extensions have been proposed, including imprecise-probability-based indices, entropy-based measures, and distribution-distance approaches such as the Bhattacharyya metric [17], [22], [23]. Nevertheless, most of these methods rely on nested Monte Carlo simulation, where epistemic variables are sampled in an outer loop and conditional aleatory propagation is performed in an inner loop. Such double-loop procedures are computationally expensive and, more importantly, entangle epistemic and aleatory effects, making their analytical separation difficult [24]–[27]. This indicates that a principled decoupling of uncertainty sources is a prerequisite for scalable sensitivity analysis.

To overcome this entanglement, an auxiliary-variable representation based on the inverse probability integral transform has been proposed to map hybrid uncertainty into a latent space where aleatory and epistemic components become independent [28], [29]. This transformation restores the independence structure required by kernel-based and ANOVA-type decompositions and provides a foundation for unified sensitivity analysis.

Building on this decoupled representation, we next consider global sensitivity analysis for high-dimensional and multivariate outputs. Classical variance-based methods rely on functional ANOVA decompositions, leading to Sobol’ indices that quantify variance contributions of individual inputs and their interactions [18], [20]. However, variance-based measures inherently capture only second-order dispersion and are insufficient to characterize full distributional effects, especially in multivariate or spatially distributed systems.

To overcome this limitation, moment-independent sensitivity measures have been proposed to quantify input importance based on distributional dissimilarities rather than variance decomposition [30]. More recently, kernel-based dependence measures, particularly the Hilbert–Schmidt Independence Criterion (HSIC), have emerged as a unified framework for global sensitivity analysis. By embedding probability distributions into reproducing kernel Hilbert spaces (RKHS) [31], HSIC admits an ANOVA-style decomposition under the assumption of mutually independent inputs [32], thereby connecting dependence measures with classical variance-based sensitivity indices.

However, existing kernel-based approaches typically reduce vector-valued outputs to scalar dependence measures via global kernels or aggregated distances. This scalarization removes internal output structure and prevents the resolution of how different output components respond to specific combinations of aleatory and epistemic uncertainty. This limitation is particularly restrictive in high-dimensional engineering systems, where different output dimensions often correspond to distinct physical mechanisms and exhibit heterogeneous sensitivity patterns.

To address these challenges, this paper proposes a double-space tensor-product RKHS framework for global sensitivity analysis under hybrid uncertainty. The method leverages the decoupled latent representation to construct separable zero-mean kernels over both input and output spaces, enabling a unified functional analytic treatment of uncertainty propagation and output structure. A dual-domain Möbius inversion scheme is introduced to decompose the total dependence measure across both input and output dimensions, yielding dimension-resolved sensitivity indices. This allows explicit characterization of cross-domain interactions between aleatory and epistemic uncertainty at the level of individual output components. From a computational perspective, the resulting formulation admits a fully vectorized single-loop implementation, eliminating the need for nested Monte Carlo simulation. This significantly improves computational efficiency while maintaining accuracy, and enables scalable application to high-dimensional and multi-output engineering problems.

The remainder of this paper is organized as follows. Section 2 introduces the theoretical background, including hybrid uncertainty representation via probability integral transforms, HSIC-based dependence measures, and properties of zero-mean kernels. Section 3 develops the proposed tensor-product RKHS framework and dual Möbius inversion decomposition. Section 4 presents two case studies, including a modified Ishigami function and a NACA 2412 aerodynamic pressure field. Finally, the appendix provides the proof of the double ANOVA decomposition theorem.

2 Backgroud↩︎

2.1 Hybrid Uncertainty Modeling Based on the Probability Integral Transform↩︎

In practical engineering systems, uncertainties are described within a fixed model universe consisting of basic random variables, probabilistic sub-models, and physical response models. Within this framework, it is convenient to distinguish between aleatory and epistemic uncertainty, where the distinction depends on how the model is constructed rather than on intrinsic properties of the physical world.

Let \(\theta \in \Theta \subset \mathbb{R}^d\) denote uncertain parameters associated with the probabilistic description of the system, reflecting incomplete knowledge in the modeling process. Conditional on \(\theta\), the system input is represented by a probabilistic model of the form \[X \sim p(X \mid \theta),\] where the distribution is interpreted as part of the chosen model structure within the model universe.

Within the same framework, the system response is defined through a deterministic physical model \[Y = g(X),\] where all uncertainty is embedded in the probabilistic characterization of \(X\) and the uncertainty in \(\theta\). In this representation, aleatory uncertainty refers to the variability of \(X\) under a fixed model specification, while epistemic uncertainty corresponds to the uncertainty in the model parameters \(\theta\), which may in principle be reduced through additional information or model refinement [5].

In the context of epistemic uncertainty, a key idea is that when information is incomplete or data are limited, probabilities do not need to be represented as single precise values. Instead, they can be modeled as intervals or sets, where uncertainty is expressed through lower and upper bounds or a family of admissible probability distributions, reflecting incomplete knowledge about the true probability measure [13].

At this stage, a useful interpretation emerges: the separation between aleatory variability and epistemic uncertainty can be further formalized through the introduction of an auxiliary variable, which enables a fully deterministic reformulation of the stochastic system[29]. In this context, the auxiliary variable is introduced to explicitly represent the intrinsic variability of the input random variable while decoupling it from distributional parameter uncertainty. For each uncertain input variable \(X\), the probability integral transform establishes a one-to-one mapping between realizations of \(X\) and a standard uniform variable \(\xi \sim \mathcal{U}(0,1)\), such that \[\xi = F_X(X \mid \theta).\]

Consequently, the input can be equivalently expressed as \[X = F_X^{-1}(\xi \mid \theta),\] where \(\xi\) serves as an auxiliary variable that is statistically independent of the epistemic parameters \(\theta\). This transformation provides two important advantages: (i) it yields an explicit representation of the aleatory variability via a distribution-free latent variable \(\xi\), and (ii) it ensures a deterministic mapping from the joint space \((\xi, \theta)\) to the system response.

In this formulation, the overall uncertainty in \(X\) is naturally decomposed into two components: the intrinsic variability captured by \(\xi\), and the epistemic uncertainty associated with the parameter vector \(\theta\). This separation aligns with the classical double-loop interpretation in uncertainty propagation, while simultaneously enabling a unified deterministic representation of the model response. As a result, the system output can be written as \[Y = g\!\left(F_X^{-1}(\xi \mid \theta)\right) = h(\xi,\theta),\] which forms the basis for global sensitivity analysis and uncertainty decomposition.

2.2 Hilbert-Schmid Independence Criterion based sensitivity analysis technique↩︎

Let \(\mathcal{X}\) and \(\mathcal{Y}\) be topological spaces endowed with Borel probability measures \(P_X\) and \(P_Y\). By introducing positive definite kernels \(k_{\mathcal{X}}\) and \(k_{\mathcal{Y}}\) with their associated Reproducing Kernel Hilbert Spaces (RKHS) \(\mathcal{H}_X\) and \(\mathcal{H}_Y\), the corresponding feature maps are given by \(\phi(X) \in \mathcal{H}_X\) and \(\psi(Y) \in \mathcal{H}_Y\). The mean embeddings of the marginal distributions are defined as[33]: \[\mu_{X} := \mathbb{E}_{X \sim P_{X}}[\phi(X)], \qquad \mu_{Y} := \mathbb{E}_{Y \sim P_{Y}}[\psi(Y)].\]

To capture the statistical dependence between \(X\) and \(Y\), the centered cross-covariance operator \(C_{XY}: \mathcal{H}_Y \rightarrow \mathcal{H}_X\) is introduced[34]: \[C_{XY} := \mathbb{E}_{XY}\big[(\phi(X) - \mu_X) \otimes (\psi(Y) - \mu_Y)\big] = \mathbb{E}_{XY}[\phi(X) \otimes \psi(Y)] - \mu_X \otimes \mu_Y,\] where \(\otimes\) denotes the tensor product.

The Hilbert-Schmidt Independence Criterion (HSIC) is defined as the squared Hilbert-Schmidt norm of this cross-covariance operator. From a distributional perspective, this operator-theoretic definition is mathematically equivalent to the Maximum Mean Discrepancy (MMD) between the joint distribution \(P_{XY}\) and the product of its marginals \(P_X P_Y\)[35], [36]: \[\mathrm{HSIC}(X,Y) := \|C_{XY}\|_{\mathrm{HS}}^{2} = \|\mu_{XY} - \mu_X \otimes \mu_Y\|_{\mathcal{H}_X \otimes \mathcal{H}_Y}^{2},\] where \(\mu_{XY} := \mathbb{E}_{XY}[\phi(X) \otimes \psi(Y)]\) is the mean embedding of the joint distribution. For characteristic kernels, independence (\(P_{XY} = P_X P_Y\)) is achieved if and only if \(\mathrm{HSIC}(X,Y) = 0\) .

Although HSIC provides a robust dependence measure, obtaining normalized sensitivity indices analogous to the classical variance-based ANOVA decomposition requires an orthogonal construction of the kernels. Specifically, the multivariate input kernel \(k_\mathcal{X}\) must be structured as a tensor product of univariate kernels [32]: \[k_\mathcal{X}(\mathbf{x},\mathbf{x}') = \prod_{l=1}^d \left(1 + k_l(x_l,x_l')\right),\] where \(d\) is the dimensionality of the inputs. Furthermore, to ensure orthogonality, each univariate kernel \(k_l(\cdot,\cdot)\) must satisfy the zero-mean condition with respect to its marginal distribution: \[\int_{\mathcal{X}_l} k_l(x_l,x_l') d\mathrm{P}_{X_l}(x_l') = 0, \quad \forall x_l \in \mathcal{X}_l.\]

Once the zero-mean condition is satisfied, the total HSIC dependence measure admits a unique orthogonal decomposition over all possible subsets of the input variables. Let \(\mathcal{D} = \{1, \ldots, d\}\) denote the full index set of the inputs. The total HSIC can be expanded as: \[\mathrm{HSIC}\left(\mathbf{X},Y\right) = \sum_{A\subseteq \mathcal{D}} \mathrm{HSIC}_A,\] where \(\mathrm{HSIC}_A\) isolates the partial dependence contribution exclusively attributed to the subset of variables \(\mathbf{X}_A\). Using the inclusion-exclusion principle, this partial contribution is computed as: \[\mathrm{HSIC}_A = \sum_{B \subseteq A}(-1)^{|A|-|B|} \mathrm{HSIC}\left(\mathbf{X}_B,Y\right).\] Here, the operator \(\mathrm{HSIC}\left(\mathbf{X}_B,Y\right)\) measures the dependence using only the subset of inputs \(B\). It is explicitly computed by replacing the full input kernel \(k_\mathcal{X}\) with the restricted sub-kernel \(k_B\): \[k_B(\mathbf{x}_B,\mathbf{x}_B') = \prod_{l\in B}\left(1+k_l(x_l,x_l')\right).\]

Analogous to the classical ANOVA decomposition, this formulation expresses the total dependence as a sum of partial dependencies. Consequently, the normalized first-order and total-effect HSIC-based sensitivity indices for any input subset \(A\) are explicitly defined as: \[\begin{align} S_A^{\mathrm{HSIC}} &= \frac{\mathrm{HSIC}_A}{\mathrm{HSIC}\left(\mathbf{X},Y\right)}, \\ S_A^{\text{Total},\mathrm{HSIC}} &= \sum_{B\subseteq \mathcal{D},\, B\cap A \neq \emptyset}S_B^{\mathrm{HSIC}} = 1 - \frac{\mathrm{HSIC}(\mathbf{X}_{\sim A},Y)}{\mathrm{HSIC}\left(\mathbf{X},Y\right)}, \end{align}\] where \(\mathbf{X}_{\sim A}\) denotes the vector of all input variables except those in subset \(A\). These indices are strictly bounded within \([0, 1]\) and satisfy the fundamental sum-to-one property: \[\sum_{A\subseteq \mathcal{D}}S_A^{\mathrm{HSIC}}=1.\]

An unbiased estimate of HSIC based on U-statistics is formulated as [37]: \[\mathrm{HSIC}_{u}(X, Y) = \frac{1}{n(n-3)} \left[ \mathrm{tr}(\tilde{K}_X \tilde{K}_Y) + \frac{\mathbf{1}^{\mathbf{T}} \tilde{K}_X \mathbf{1} \mathbf{1}^{\mathbf{T}} \tilde{K}_Y \mathbf{1}}{(n-1)(n-2)} - \frac{2}{n-2} \mathbf{1}^{\mathbf{T}} \tilde{K}_X \tilde{K}_Y \mathbf{1} \right],\] where \(\mathbf{1}\) is a vector of ones, \(n\) is the sample size, and \(\tilde{K}_X\) and \(\tilde{K}_Y\) correspond to the Gram matrices where all the diagonal elements have been set equal to zero: \[\tilde{K}_{ij} = (1 - \delta_{ij})K_{ij}.\]

2.3 Zero-mean Characteristic Kernels↩︎

In kernel-based learning, a valid kernel function must be positive definite to guarantee the existence of a Reproducing Kernel Hilbert Space (RKHS). Furthermore, for independence testing and dependency quantification (such as HSIC), it is crucial that the feature map has a zero expectation. This requires the application of zero-mean (or centered) kernels.

Analytically, this zero-mean condition can be fulfilled starting from an arbitrary positive definite kernel \(k(\cdot,\cdot)\). A theoretically centered kernel \(k_0^D(\cdot,\cdot)\) can be derived via integration over the marginal distribution \(\mathrm{P}_{X_l}\) [38]: \[k_0^D(x_l,x_l') = k(x_l,x_l') - \frac{\int k(x_l,s)d\mathrm{P}_{X_l}(s)\int k(x_l',s)d\mathrm{P}_{X_l}(s)}{ k(s,t)d\mathrm{P}_{X_l}(s)d\mathrm{P}_{X_l}(t)}. \label{eq:zero95mean95kernel}\tag{1}\]

Alternatively, instead of explicitly computing the integral, a zero-mean RKHS can be directly defined by choosing specific kernel families [39]. More recently, several works have made use of the Stein operator to define the Stein discrepancy in an RKHS [40], which shows great potential for Monte-Carlo integration [41] or goodness-of-fit tests [42], [43] when the target distribution is either impossible to sample or is known only up to a normalization constant.

However, computing these integrals explicitly is generally intractable. In practice, to capture complex dependencies, we rely on powerful characteristic kernels. These can be categorized into kernels that are naturally zero-mean and those that require explicit empirical centering [44].

A class of distance-induced characteristic kernels inherently satisfies the zero-mean property without requiring explicit centering operations [45]. They are defined as: \[k_d^{\alpha}(r,r') := \frac{1}{2} \left( \|r\|_2^{2\alpha} + \|r'\|_2^{2\alpha} - \|r-r'\|_2^{2\alpha} \right), \quad \forall \alpha \in [0,2].\] This elegant formulation avoids Equation (1 ) entirely.

3 Methodology↩︎

3.1 Extension to High-Dimensional Responses via Tensor-Product Kernels↩︎

To isolate dependence contributions associated with different subsets of variables, we adopt the factorized ANOVA kernel over the independent space: \[k_Z(z,z') = \prod_{i=1}^{m} \left( 1+k_{\xi_i}(\xi_i,\xi_i') \right) \prod_{j=1}^{d} \left( 1+k_{\Theta_j}(\theta_j,\theta_j') \right).\] Crucially, because \(\xi\) and \(\Theta\) are mutually independent by construction, each univariate kernel can be strictly centered with respect to its marginal distribution: \[\int k_{\xi_i}(\xi,\xi')\, dP_{\xi_i}(\xi') = 0, \qquad \int k_{\Theta_j}(\theta,\theta')\, dP_{\Theta_j}(\theta') = 0,\] for all \(i=1,\dots,m\) and \(j=1,\dots,d\). This guarantees that the cross-covariance operators associated with different subsets are strictly orthogonal, thereby salvaging the mathematical foundation of the HSIC-ANOVA sensitivity analysis.

When the system response is high-dimensional, i.e., \(Y = (Y_1, \dots, Y_p) \in \mathbb{R}^p\), utilizing a single global radial basis function kernel with a lumped Euclidean distance metric becomes suboptimal due to distance concentration in high-dimensional spaces. To maintain strict mathematical symmetry and consistency with the input space decomposition, we construct the output kernel using the same orthogonal ANOVA tensor-product architecture.

Following the identical logic used for the independent input space \(Z\), the global reproducing kernel for the multidimensional output space \(\mathcal{Y}\) is constructed as a tensor product of the component-wise kernels: \[k_Y(y, y') = \prod_{k=1}^p \left( 1 + k_{Y_k}(y_k, y_k') \right).\] This factorized formulation explicitly embeds the high-dimensional response into a tensor-product Reproducing Kernel Hilbert Space (RKHS).

By virtue of the factorized tensor-product architecture applied to both the input space \(Z\) and the output space \(Y\), the HSIC dependence measure can be structurally expanded across both domains. Let \(V \subseteq \{1,\dots,p\}\) denote a subset of the output dimensions, and let \(Y_V\) be the corresponding subvector of the system response.

Similar to the input space, we define the subset-specific cross-dependence measure: \[\mathrm{HSIC}(Z_A, Y_V) := \left\| C_{Z_A Y_V} \right\|_{\mathrm{HS}}^2.\]

The classical ANOVA decomposition extracts orthogonal interaction effects for multivariate functions using Möbius inversion over a single subset lattice [46]. In this paper, we extend this fundamental principle into a dual-space reproducing kernel framework. By applying the Möbius inversion concurrently over the subset lattices of both the input dimensions \(A\) and the output dimensions \(V\), the pure cross-domain interaction effect between a specific subset of inputs and a specific subset of outputs can be isolated as follows.

To obtain a scale-free and comparable measure of dependence contribution mapping the propagation of specific hybrid uncertainties onto specific structural modes of the high-dimensional response, we define the normalized index: \[S_{A, V}^{\mathrm{HSIC}} := \frac{\mathrm{HSIC}_{A, V}}{\mathrm{HSIC}(Z,Y)},\] and \[\sum_{A \neq \emptyset}\sum_{V \neq \emptyset} S_{A,V}^{\mathrm{HSIC}} = 1.\]

3.2 Multi-Output HSIC Sensitivity Indices and Statistical Testing↩︎

Building upon the double ANOVA decomposition established in Theorem 1, the empirical estimation of the sensitivity indices is executed through the V-statistics of the distance-induced kernels. Let \(Z = (\xi_1, \dots, \xi_m, \Theta_1, \dots, \Theta_d)\) denote the \(m+d\) independent inputs (indexed by \(\mathcal{N}\)), and \(Y = (Y_1, \dots, Y_p)\) denote the \(p\)-dimensional system response (indexed by \(\mathcal{P}\)). The total variance of the joint multi-output system is quantified by the global dependence measure \(V_{\mathrm{global}} = \mathrm{HSIC}(Z, Y)\).

To comprehensively disentangle the uncertainty propagation across different structural levels, the algorithm evaluates three categories of sensitivity indices derived directly from the Möbius inversion framework:

1. Global Input Decomposition Indices: This category evaluates the contribution of individual input factors to the holistic system state, treating the multidimensional output \(Y\) as an indivisible joint entity. By setting the output subset \(V = \mathcal{P}\), the system-level main effect for a specific input variable \(Z_i\) (\(i \in \mathcal{N}\)) is given by: \[S_{Z_i}^{\mathrm{sys}} = \frac{\mathrm{HSIC}(\{Z_i\}, Y)}{V_{\mathrm{global}}}.\] The aggregated pure higher-order interaction effect encompassing all cross-domain (aleatory-epistemic) and intra-domain couplings is seamlessly isolated by subtracting all first-order main effects from the total normalized variance: \[S_{\mathrm{Interaction}}^{\mathrm{sys}} = 1 - \sum_{i=1}^{m+d} S_{Z_i}^{\mathrm{sys}}.\]

2. Global Output Decomposition Indices: To reveal the intrinsic topological coupling within the high-dimensional response, the total system variance is further decomposed from the perspective of the output space. By setting the input subset \(A = \mathcal{N}\), the independent main effect of a specific output dimension \(Y_k\) (\(k \in \mathcal{P}\)), and the pure synergistic interaction effect between outputs \(Y_k\) and \(Y_l\), are quantified respectively as: \[S_{Y_k}^{\mathrm{main}} = \frac{\mathrm{HSIC}(Z, \{Y_k\})}{V_{\mathrm{global}}}, \qquad S_{Y_k, Y_l}^{\mathrm{int}} = \frac{\mathrm{HSIC}_{\mathcal{N}, \{k, l\}}}{V_{\mathrm{global}}},\] where \(\mathrm{HSIC}_{\mathcal{N}, \{k, l\}}\) is explicitly extracted via the Möbius inversion. These indices isolate the fraction of total uncertainty absorbed by the independent fluctuation of individual outputs versus the coupled resonance between them.

3. Dimension-wise Decomposition Indices with Universal Scale: For granular cross-dimensional comparison, it is mathematically imperative to establish a universal denominator that avoids the distortion caused by inter-output dependencies (i.e., the synergistic terms \(S_{Y_k, Y_l}^{\mathrm{int}}\)). We define the universal scale as the sum of all independent output variances: \[V_{\mathrm{uni}} = \sum_{k=1}^p \mathrm{HSIC}(Z, \{Y_k\}).\] For any target output dimension \(Y_k\), the dimension-wise main effect of input \(Z_i\) and the explicit second-order interaction between inputs \(Z_i\) and \(Z_j\) are formulated as: \[S_{Z_i \to Y_k} = \frac{\mathrm{HSIC}(\{Z_i\}, \{Y_k\})}{V_{\mathrm{uni}}}, \qquad S_{Z_i, Z_j \to Y_k} = \frac{\mathrm{HSIC}_{\{i, j\}, \{k\}}}{V_{\mathrm{uni}}}.\] By scaling against the identical baseline \(V_{\mathrm{uni}}\), the indices \(S_{Z_i \to Y_k}\) inherently preserve the absolute magnitude of uncertainty transferred to each specific dimension. This guarantees that the derived dimension-wise interaction matrix \(\mathbf{I}_{i,j}^{(k)} = S_{Z_i, Z_j \to Y_k}\) can be directly and fairly compared across all \(p\) output dimensions, ensuring an unbiased assessment of critical structural regimes.

The proposed Double-Space Tensor-Product RKHS sensitivity analysis procedure is summarized in Algorithm 1.

For statistical validation, we employ a non-parametric permutation test. Under the null hypothesis \(H_0: S_i = 0\), we construct an empirical null distribution by performing \(N_p\) random permutations of the response kernel matrix \(\widetilde{K}_Y\), yielding the p-value: \[\text{p-value} = \frac{1}{N_p} \sum_{p=1}^{N_p} \mathbb{I}\left( \widehat{S}_i^{(p)} \geq \widehat{S}_i \right).\]

To quantify sampling variability, we further construct \((1-\alpha)\) confidence intervals via bootstrap resampling. By drawing \(N_b\) samples with replacement, the confidence interval is given by empirical quantiles: \[\text{CI}_{1-\alpha} = \left[ q_{\alpha/2}\left( \{\widehat{S}_i^{(b)}\}_{b=1}^{N_b} \right), \, q_{1-\alpha/2}\left( \{\widehat{S}_i^{(b)}\}_{b=1}^{N_b} \right) \right].\]

3.3 Isoprobabilistic Transformation-based Algorithmic Implementation↩︎

To circumvent the intractable computational bottleneck imposed by traditional double-loop nested Monte Carlo simulations, we propose a fully vectorized, single-loop algorithmic architecture based on the isoprobabilistic transformation (Inverse Probability Integral Transform, PIT) [47]. This framework evaluates the physics model in the physical parameter space while conducting the rigorous RKHS variance decomposition strictly within the decoupled, orthogonal latent space. The implementation consists of four sequential stages:

1. Joint Sampling in the Latent Space. To ensure optimal space-filling properties and avoid spurious correlations, a joint sampling strategy is executed in a standard uniform latent space [48]. Let \(N\) be the total computational budget. A Latin Hypercube Sampling (LHS) design is utilized to generate an \(N \times (d_\xi + d_\theta)\) sample matrix \(\mathbf{U}\) drawn from an independent standard uniform distribution \(\mathcal{U}(0,1)^{d_\xi + d_\theta}\). This matrix is structurally partitioned into two decoupled blocks: the latent aleatory subset \(\mathbf{U}_\xi \in \mathbb{R}^{N \times d_\xi}\) and the latent epistemic subset \(\mathbf{U}_\Theta \in \mathbb{R}^{N \times d_\theta}\).

2. Inverse Isoprobabilistic Mapping. The latent uniform samples are deterministically mapped to the physical domain, formulating the physical aleatory matrix \(\mathbf{X}\) and epistemic matrix \(\mathbf{\Theta}\). Depending on the complexity of the uncertainty models, three mapping paradigms are incorporated:

Independent Marginals: For variables with arbitrary independent distributions (e.g., Beta, Weibull, or Uniform), the physical column vector \(\mathbf{x}_i \in \mathbb{R}^{N \times 1}\) is obtained directly via the element-wise inverse Cumulative Distribution Function (CDF) applied to the corresponding latent column \(\mathbf{u}_{\xi, i}\): \[\mathbf{x}_i = F_{X_i}^{-1}(\mathbf{u}_{\xi, i}).\]

Hierarchical/Conditioned Variables: For aleatory variables whose distribution bounds or moments are dictated by epistemic parameters (forming a probability box or hierarchical uncertainty), the mapping is conditioned on the mapped epistemic realization \(\mathbf{\Theta}\): \[\mathbf{x}_i = F_{X_i|\Theta}^{-1}(\mathbf{u}_{\xi, i} \mid \mathbf{\Theta}).\]

Correlated Multivariate Distributions: Let \(\mathbf{U}_c \in \mathbb{R}^{N \times d_c}\) denote a \(d_c\)-dimensional subset of the latent uniform samples intended for correlated aleatory variables. To embed the dependence structure via a Gaussian copula [49], a component-wise inverse probability integral transform is first applied to obtain an independent standard normal representation: \[\mathbf{Z}_{\mathrm{ind}} = \Phi^{-1}(\mathbf{U}_c),\] where \(\Phi^{-1}(\cdot)\) is the inverse standard normal CDF, yielding \(\mathbf{Z}_{\mathrm{ind}} \sim \mathcal{N}(\mathbf{0}, \mathbf{I})\).

Given the target correlation matrix \(\mathbf{R} \in \mathbb{R}^{d_c \times d_c}\), the covariance structure is factorized via Cholesky decomposition [50] as \(\mathbf{R} = \mathbf{L}\mathbf{L}^T\). The correlated Gaussian matrix is then generated via the forward affine transformation: \[\mathbf{Z}_{\mathrm{corr}} = \mathbf{Z}_{\mathrm{ind}} \mathbf{L}^T,\] ensuring the column vectors follow \(\mathcal{N}(\mathbf{0}, \mathbf{R})\). Finally, the physical samples \(\mathbf{X}_c\) are recovered through the target inverse marginal distributions applied to the correlated uniform space: \[\mathbf{X}_c = F_c^{-1}(\Phi(\mathbf{Z}_{\mathrm{corr}})).\]

3. Vectorized System Evaluation. The fully formulated physical sample matrix \(\mathbf{X}\) and the epistemic parameters \(\mathbf{\Theta}\) are fed into the deterministic computational model \(\mathcal{M}\) to yield the high-dimensional system response matrix \(\mathbf{Y} \in \mathbb{R}^{N \times d_y}\): \[\mathbf{Y} = \mathcal{M}(\mathbf{X}, \mathbf{\Theta}).\] Because the input matrix is generated via joint LHS and mapped via explicit matrix operations, this stage is entirely free of nested for-loops, enabling extreme parallelization and vectorized matrix operations.

4. Latent-Space Kernel Estimation and Statistical Inference. To strictly preserve the mutual independence required by the double-space ANOVA decomposition (as established in Theorem 1), the sensitivity indices must be computed using the latent independent matrices (i.e., \(\mathbf{U}_\xi\) and \(\mathbf{\Theta}\)) alongside the evaluated system output matrix \(\mathbf{Y}\). It is crucial to emphasize that the physical aleatory matrix \(\mathbf{X}\) is intentionally bypassed during the RKHS kernel construction. This strategic decoupling precludes any spurious collinearity or structural dependencies induced by the hierarchical epistemic bounds or Copula correlations. The standardized latent inputs and physical outputs are subsequently fed into Algorithm 1 to compute the unbiased V-statistics of the distance-induced Gram matrices. Finally, the empirical estimators are subjected to the aforementioned non-parametric permutation tests for significance screening, and the \(N_b\)-iteration bootstrap resampling is executed to establish the \((1-\alpha)\) confidence intervals, thereby completing the robust uncertainty quantification cycle.

4 Case↩︎

4.1 Case Study 1: Multi-Output Modified Ishigami Function↩︎

To verify the proposed Double-Space Tensor-Product RKHS sensitivity analysis framework, we first investigate a modified multi-output Ishigami function [51]. The classical Ishigami function is widely used as a benchmark in global sensitivity analysis due to its strong nonlinearity and non-monotonicity. Here, it is extended to a multi-dimensional response space subjected to complex hybrid uncertainty. In this hybrid formulation, the epistemic parameters act as hyper-parameters that directly govern the distribution boundaries of the intrinsic aleatory variables.

The system involves two epistemic parameters \(\mathbf{\Theta} = [\theta_1, \theta_2]\) (denoted as \(a\) and \(b\)) and three intrinsic aleatory variables \(\mathbf{\Xi} = [\xi_1, \xi_2, \xi_3]\). The epistemic parameters are assigned independent uniform distributions: \[a \sim \mathcal{U}(6, 8), \qquad b \sim \mathcal{U}(0.05, 0.15).\] The physical aleatory variables \(\mathbf{X} = [X_1, X_2, X_3]\) are subsequently mapped via conditional inverse probability integral transforms (Inverse PIT). To construct the hybrid uncertainty structure, their physical domains are strictly dictated by the realizations of the epistemic parameters \(a\) and \(b\): \[X_1 \sim \mathcal{U}(-a, a), \qquad X_2 \sim \mathcal{U}(-\pi+b, \pi+b), \qquad X_3 \sim \mathcal{U}(-ab, ab).\]

For the algorithmic implementation, a joint Latin Hypercube Sampling (LHS) design with \(N = 1000\) is executed in the 5-dimensional standard uniform latent space to generate independent \(\mathbf{\Xi}\) and \(\mathbf{\Theta}\).

The corresponding multi-output nonlinear system is defined as a black-box forward model with three response components, capturing heterogeneous nonlinear interactions between epistemic and aleatory uncertainties: \[\begin{align} Y_1 &= \sin(X_1) + a \, \sin^2(X_2) + b \, X_3^4 \sin(X_1), \nonumber \\ Y_2 &= \cos(a X_1) - b X_2, \nonumber \\ Y_3 &= X_1^2 + X_2^2 + X_3^2. \end{align}\]

This construction extends the classical Ishigami function by introducing (i) epistemic modulation of nonlinear amplitude through \((a,b)\), (ii) cross-scale coupling between epistemic and aleatory variables, and (iii) a mixture of periodic, polynomial, and quadratic response structures, thereby providing a challenging benchmark for tensor-product RKHS-based sensitivity decomposition under hybrid uncertainty.

To rigorously quantify the statistical sampling variability, a bootstrap resampling procedure with \(N_b = 300\) iterations is performed to establish the 95% confidence intervals for all sensitivity estimators.

Figure 2 illustrates the holistic system-level sensitivity indices evaluated via the global tensor-product kernels (\(K_{\mathcal{X}}\) and \(K_{\mathcal{Y}}\)). As shown in Figure 2 (a), the global input decomposition robustly identifies the overall influence of the coupled inputs. Most remarkably, the interaction term (“Inter”) single-handedly dominates the system variance (capturing approximately \(70\%\)). This validates the profound structural impact of the hybrid uncertainty framework, where the hierarchical epistemic bounds heavily intertwine with the intrinsic aleatory fluctuations. Figure 2 (b) provides a novel topological perspective by decomposing the output space. It delineates how the total system uncertainty is partitioned. It is evident that the independent fluctuation of \(Y_3\) governs the largest share of variance (\(>0.3\)). Furthermore, the synergistic cross-couplings involving \(Y_3\) (i.e., \(Y_{13}\) and \(Y_{23}\)) exhibit higher variance contributions than the independent main effects of \(Y_1\) and \(Y_2\) themselves. This macro-level output decomposition justifies the necessity of granular dimension-wise analyses, as it proves the output dimensions resonate on vastly disparate scales.

To ensure a mathematically rigorous cross-dimensional comparison, the dimension-wise input sensitivities are normalized against the Universal Scale (\(V_{\mathrm{uni}}\)) and presented in Figure 3. Unlike traditional local normalization that distorts absolute magnitudes, this scaling preserves the true share of uncertainty absorbed by each output. The grouped bar charts intuitively reveal that the inputs dynamically shift their dominance across different structural modes. For instance, the aleatory variable \(X_1\) dictates a massive variance proportion exclusively in \(Y_3\) (approaching 0.17) due to its quadratic formulation, while remaining negligible in \(Y_1\) and \(Y_2\). Conversely, \(X_2\) exhibits moderate, comparable main effects on both \(Y_1\) and \(Y_2\). Across all dimensions, the interaction effects unequivocally govern the uncertainty propagation, peaking significantly in \(Y_3\) (\(\sim 0.4\)). The exceptionally narrow 95% confidence intervals (black error bars) confirm the high numerical stability of the proposed vectorized V-statistics estimator at \(N=1000\).

The explicit second-order interaction effects extracted via the dual-space Möbius inversion are visualized as heatmaps in Figure 4. The strictly lower-triangular mask ensures clarity, with the absolute values of \(S_{i,j}^{(k)}\) annotated. The algorithm successfully isolates the theoretical physical couplings induced by the hybrid uncertainties. For output \(Y_1\) (Figure 4a), the strongest inter-domain interaction is accurately detected between the latent variables corresponding to \(X_2\) and \(\theta_1\) (\(0.020\)), perfectly mirroring the physical multiplication mechanism in the \(a \sin^2(X_2)\) term.

More importantly, the heatmap for \(Y_3\) (Figure 4c) highlights the paramount advantage of the proposed latent-space RKHS framework. A naive observation of the physical equation \(Y_3 = X_1^2 + X_2^2 + X_3^2\) might falsely presume zero cross-interactions due to its additive structure. However, because the integration bounds are epistemically conditioned (e.g., \(X_1 \sim \mathcal{U}(-a, a)\)), the isoprobabilistic mapping enforces a multiplicative coupling in the independent latent space (e.g., \(X_1^2 = a^2(2\xi_1 - 1)^2\)). The proposed algorithm flawlessly unveils this hidden hybrid dependency, capturing massive cross-domain interactions, prominently the \(0.058\) interaction magnitude between \(X_1\) and \(\theta_1\). This confirms the strict orthogonality and unprecedented diagnostic power of the proposed tensor-product HSIC-ANOVA methodology.

4.2 Case Study 2: NACA 2412 Aerodynamic Pressure Field Analysis↩︎

To validate the capability of the proposed method in handling high-dimensional, continuous spatial field outputs driven by computationally intensive black-box solvers, an aerodynamic uncertainty quantification (UQ) benchmark problem is investigated. In modern aerodynamic design, the accurate prediction of the surface pressure coefficient (\(C_p\)) distribution is critical for evaluating lift, drag, and pitching moments. However, the aerodynamic system is inherently subjected to profound hybrid uncertainties: the objective operational flight conditions exhibit irreducible aleatory randomness, while the semi-empirical physical sub-models (e.g., boundary layer transition criteria) introduce inevitable epistemic imprecision due to a lack of perfect knowledge.

Prior to conducting the high-dimensional mixed sensitivity analysis, a baseline aerodynamic verification is performed to ensure the fidelity of the underlying numerical solver. Figure 5 illustrates the geometric profile of the selected airfoil and its corresponding surface pressure distribution under nominal flight conditions, where the deterministic computational model \(\mathcal{M}\) utilized in this study is XFOIL [52], an extensively validated high-fidelity subsonic viscous solver based on the panel method coupled with an integral boundary layer formulation.

As depicted in Fig. 5 (a), the continuous profile of the standard NACA 2412 airfoil is reconstructed, where the solid gray-shaded region represents the internal material domain of the wing section. The black dash-dot line delineates the mean camber line, capturing the maximum camber of \(2\%\) located at \(40\%\) of the chord length from the leading edge. To prevent any artificial visual distortion of the geometric thickness-to-chord ratio (\(12\%\)), the aspect ratio of the ordinate to the abscissa is strictly constrained to be equal.

Figure 5 (b) presents the baseline steady-state surface pressure coefficient (\(C_p\)) distribution computed by the viscous-coupled panel solver at the nominal operational condition (\(\alpha = 3.0^\circ\), \(Mach = 0.3\), and \(Re = 3\times 10^6\)). Adhering to standard aerodynamic conventions, the vertical axis representing \(C_p\) is inverted, mapping negative pressure (suction) toward the upper part of the grid.

The lower (pressure) surface and upper (suction) surface distributions are explicitly distinguished by the deep blue and orange-red solid lines, respectively. Due to the positive angle of attack, the upper surface curve exhibits a prominent, localized suction peak (\(C_p \approx -1.8\)) near the leading edge (\(x/c \approx 0.08\)), which serves as the primary mechanism for lift generation.

Furthermore, \(30\) discrete evaluation nodes uniformly distributed within the chordwise range of \(x/c \in [0.05, 0.95]\) are superimposed on the upper surface curve as black solid circles. These nodes explicitly define the extraction locations of the spatial continuous field output. Rather than reducing the flow diagnostics to a single lumped scalar property (such as the lift coefficient), the subsequent global sensitivity analysis will be executed across this 30-dimensional spatial vector to unravel the localized propagation of hybrid uncertainties.

The specific uncertainty characteristics and their corresponding transformations from the orthogonal latent space \(\mathbf{U} = [\mathbf{\Xi}, \mathbf{\Theta}]\) are summarized in Table 1.

The output of the system is the continuous chordwise pressure coefficient distribution over the upper surface of the airfoil. To embed this spatial field into the proposed multi-output framework, the continuous function is discretized at \(d_y = 30\) uniformly spaced evaluation nodes along the chord \(x/c \in [0.05, 0.95]\). Consequently, the output response is mathematically formulated as a high-dimensional matrix \(\mathbf{Y} = [C_{p}^{(1)}, \dots, C_{p}^{(30)}] \in \mathbb{R}^{N \times 30}\).

Based on the defined hybrid uncertainty framework, the proposed double-space tensor-product RKHS algorithm is executed. To ensure statistical reliability and avoid spurious evaluations from non-convergent XFOIL runs, an oversampled Latin Hypercube Sampling (LHS) strategy is utilized in the 11-dimensional standard uniform latent space, securing \(N=1000\) fully converged aerodynamic evaluations. The statistical confidence intervals are established via \(N_b=100\) Bootstrap resamplings at a \(95\%\) confidence level.

Figure 6 presents the holistic system-level total sensitivity indices, where the 30-dimensional continuous pressure field is treated as an indivisible joint entity. To provide a comprehensive measure of each variable’s overall influence, these total indices aggregate the main effects with the spatial average of all their associated higher-order interactions. The results reveal that the global aerodynamic uncertainty is profoundly dominated by the latent aleatory variable associated with the angle of attack (\(\xi_\alpha\)). In aerodynamic physics, the angle of attack exponentially alters the stagnation point and the leading-edge suction peak, thereby governing the entire pressure integration. Furthermore, the epistemic hyperparameters controlling the angle of attack (\(\mu_\alpha\) and \(\sigma_\alpha\)) exhibit substantial total effects, rigorously validating that subjective bounds on operational conditions critically shape the global uncertainty boundary. The narrow Bootstrap confidence intervals (black error bars) confirm the high numerical robustness of the proposed vectorized V-statistics estimator under complex black-box solvers.

To unravel the localized propagation mechanisms of these hybrid uncertainties, Figure 7 visualizes the spatial evolution of the dimension-wise main effect indices along the upper surface chord (\(x/c \in [0.05, 0.95]\)). By explicitly filtering out the interaction components, this stacked area plot provides a highly granular topological diagnosis of the independent contributions that traditional scalar-aggregated GSA fails to capture. Near the leading edge (\(x/c < 0.2\)), the main effect variance is predominantly dictated by \(\xi_\alpha\) and its epistemic parameters (\(\mu_\alpha\), \(\sigma_\alpha\)), corresponding to the high-gradient suction peak region where flow acceleration is most severe. As the evaluation node moves downstream towards the mid-chord and trailing edge, the aerodynamic influence undergoes a distinct regime shift. The independent sensitivities of the Mach number (\(\xi_{Ma}\)) and Reynolds number (\(\xi_{Re}\)), along with the boundary layer transition epistemic parameters (\(N_{crit}\) and \(x_{trip}\)), gradually expand. This explicitly captures the physical transition of uncertainty propagation from potential flow dominance (inviscid) at the leading edge to viscous boundary layer development and pressure recovery at the aft section.

The paramount advantage of the proposed double-space Möbius inversion is its ability to isolate pure cross-domain interaction effects. Figure 8 extracts the explicit \(11 \times 11\) second-order interaction matrices (\(S_{i,j}^{(k)}\)) at three representative stations: the leading edge (\(x/c \approx 0.1\)), mid-chord (\(x/c = 0.5\)), and trailing edge (\(x/c = 0.9\)). The heatmaps demonstrate that the coupling topologies are highly heterogeneous across the spatial field. At the leading edge (Fig. 8a), strong interactions are exclusively clustered between the latent aleatory variable \(\xi_\alpha\) and its governing epistemic parameters (\(\mu_\alpha\) and \(\sigma_\alpha\)), perfectly reflecting the hierarchical mapping mechanism of the inverse probability integral transform. In contrast, at the mid-chord and trailing edge (Fig. 8b and 8c), denser off-diagonal synergistic couplings emerge between the fluid-dynamic variables (e.g., \(\xi_{Re}\), \(\xi_{Ma}\)) and the semi-empirical epistemic bounds (\(N_{crit}\), \(x_{trip}\)). This mathematical decomposition flawlessly diagnoses the complex, hidden non-linearities embedded within the XFOIL viscous-inviscid interaction solver, proving the rigorous analytical decoupling capability of the proposed RKHS framework under high-dimensional physical constraints.

5 Conclusion↩︎

In this paper, we propose a double-space tensor-product Reproducing Kernel Hilbert Space (RKHS) framework for hybrid aleatory and epistemic uncertainty quantification in high-dimensional systems. The main contributions of this work can be summarized as follows:

1. An auxiliary-variable representation based on the inverse probability integral transform is introduced to map coupled hybrid uncertainties into a decoupled latent space, thereby satisfying the independence requirements for variance-based decomposition.

2. Classical functional ANOVA is extended simultaneously to both input and output domains via factorized zero-mean kernels, resulting in a double Möbius inversion decomposition. This enables orthogonal decomposition of dependence measures without scalarizing multivariate outputs, unlike classical HSIC-based approaches.

3. The proposed framework allows explicit decomposition of global dependence into pure aleatory effects, pure epistemic effects, and their cross-domain interactions at the level of individual output components and structural modes.

4. The framework eliminates nested Monte Carlo simulations by adopting a fully vectorized single-loop implementation. The physical model is evaluated in the original space, while kernel-based decomposition is performed in the latent space, significantly improving computational efficiency.

5. The proposed method is validated on a modified multi-output Ishigami function and a NACA 2412 aerodynamic pressure field, demonstrating its capability to capture complex hybrid interactions and scale to high-dimensional spatial outputs.

6 Proof of Theorem 1↩︎

References↩︎