{'fraction_non_alphanumeric': 0.03242924528301887, 'fraction_numerical': 0.011658447684391081, 'mean_word_length': 4.798445998445999, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Subspace-based dimension reduction for chemical kinetics applications with epistemic uncertainty

We focus on an efficient approach for quantification of uncertainty in complex chemical reaction networks with a large number of uncertain parameters. Parameter dimension reduction is accomplished by computing an active subspace that predominantly captures the variability in the quantity of interest (QoI). In the present work, we compute the active subspace for a H2/O2 mechanism that involves 19 chemical reactions, using an efficient iterative strategy. The active subspace is first computed for a 19-parameter problem wherein only the uncertainty in the pre-exponents of the individual reaction rates is considered. This is followed by the analysis of a 33-dimensional case wherein the activation energies are also considered uncertain. In both cases, a 1-dimensional active subspace is identified, which indicates enormous potential for efficient statistical analysis of complex chemical systems. In addition, we explore links between active subspaces and global sensitivity analysis, and exploit these links for identification of key contributors to the variability in the model response.


Introduction
Time evolution of a chemically reacting system is largely dependent upon rate constants associated with individual reactions. The rate constants are typically assumed to exhibit a certain correlation with temperature (e.g., Arrhenius-type). Hence, accurate specification of the rate-controlling parameters is critical to the fidelity of simulations. However, in practical applications, these parameters are either specified using expert knowledge or estimated based on a regression fit to a set of sparse and noisy data [1][2][3][4]. Intensive research efforts in recent years within the field of uncertainty quantification (UQ) address the quantification and propagation of uncertainty in system models due to inadequate data, parametric uncertainty, and model errors [5][6][7][8][9][10][11].
In complex mechanisms involving a large number of reactions, characterizing the propagation of uncertainty from a large set of inputs to the model output is challenging due to the associated computational effort. A major focus of this article is the implementation of a robust framework that aims to identify important directions in the input space that predominantly capture the variability in the model output. These directions, which constitute the so called active subspace [12], are the dominant eigenvectors of a matrix derived from the gradient information of the model output with respect to an input. The active subspace methodology thus focuses on reducing the dimensionality of the problem, and hence the computational effort associated with uncertainty propagation.
The application problem considered in this work is the H 2 /O 2 reaction mechanism from [13]. This mechanism has gained a lot of attention as a potential source of clean energy for locomotive applications [14], and more recently in fuel cells [15,16]. The mechanism involves 19 reactions including chain reactions, dissociation/recombination reactions, and formation and consumption of intermediate species; see Table 1. For each reaction, the reaction rate is assumed to follow an Arrhenius correlation with temperature: where A i is the pre-exponent, n i is the index of T , E a,i is the activation energy corresponding to the i th reaction, and R is the universal gas constant. mechanism can be considered as follows: The equivalence ratio (φ) is given as follows: where the numerator on the right-hand-side denotes the ratio of the fuel (H 2 ) and oxidizer (O 2 ) at a given condition to the same quantity under stoichiometric conditions. In this study, computations were performed at fuel-rich conditions, φ = 2.0. Homogeneous ignition at constant pressure is simulated using the TChem software package [17], and the time required for the rate of temperature increase to exceed a given threshold, regarded as ignition delay is recorded. We seek to understand the impact of uncertainty in the rate-controlling parameters, preexponent (A i ) and the activation energy (E a,i ) on the ignition delay. The A i 's associated with all reactions and the E a,i 's with non-zero nominal estimates are considered to be uniformly distributed about their nominal estimates provided in [13]. The total number of uncertain inputs is 33 which makes the present problem computationally challenging due to the large number of uncertain parameters. To address this challenge, we focus on reducing the dimensionality of the problem by computing the aforementioned active subspace. This involves repeated evaluations of the gradient of a model output with respect to the input parameters. Several numerical techniques are available for computing the gradient. Perturbation-based techniques include finite differences and more advanced methods involving adjoints [18][19][20]. The adjoint-based method requires a solution of the state equation (forward solve) and the corresponding adjoint equation. Hence, it is limited by the availability of an adjoint solver. Additional model evaluations at neighboring points are required if finite difference is used which increases the computational effort. Regression-based techniques aim to estimate the gradient by approximating the model output using a regression fit. These are computationally less intensive than the former. However, as expected, there is a trade-off between computational effort and accuracy in the two approaches for estimating the gradient.
In this work, we adopt an iterative strategy to reduce the computational effort associated with active subspace computation. Moreover, we explore two approaches for estimating the gradient of the ignition delay with respect to the uncertain rate-controlling parameters: pre-exponent, A i and the activation energy, E a,i , for i th reaction in the H 2 /O 2 mechanism. As mentioned above, the first approach uses finite difference to estimate the gradient and is regarded as the perturbation-based approach, whereas, the second approach involves a linear regression-fit to the available set of model evaluations and is regarded as the regression-based approach.
An alternate strategy to dimension reduction involves computing the global sensitivity measures associated with the uncertain inputs of a model. Depending upon the estimates of the sensitivity measures, only the important inputs are varied for the purpose of uncertainty quantification (UQ). Sobol' indices are commonly used as global sensitivity measures [21]. They are used to quantify the relative contributions of the uncertain inputs to the variance in the output, either individually, or in combination with other inputs. Multiple efforts have focused on efficient computation of the Sobol' indices [22][23][24][25] including the derivative-based global sensitivity measures (DGSMs), developed to compute approximate upper bounds for the Sobol' indices with much fewer computations [26,27]. It was noted in [28,29] that DGSMs can be approximated by exploiting their links with active subspaces. This led to the definition of the so-called activity scores. In Section 3, we build on these ideas to provide a complete analysis of links between Sobol indices, DGSMs, and activity scores for functions of independent random inputs whose distribution law belongs to the class of Boltzmann probability measures. It is worth mentioning that computing global sensitivity measures provides important information about a model that go beyond dimension reduction. By identifying parameters with significant impact on the model output, we can assess regimes of validity of the model formulation, and gain critical insight into the underlying physics in many cases.
The main contributions of this paper are as follows: • Active subspace discovery in a high-dimensional H 2 /O 2 kinetics problem involving 33 uncertain inputs: The methodology presented in this work resulted in a 1-dimensional active subspace, indicating immense potential for computational savings. The presented analysis can also guide practitioners in other problems of chemical kinetics on using the method of active subspaces to achieve efficiency in uncertainty propagation.
• Comprehensive numerical investigation of the perturbation-based and the regression-based approaches: We investigate the suitability of both approaches for estimating the gradient of ignition delay in the H 2 /O 2 mechanism. Specifically, we compare resulting active subspaces, surrogate models, and the ability to approximate global sensitivity measures through a comprehensive set of numerical experiments. Our results reveal insight into the merits of the methods as well as their shortcomings.
• Analysis of the links between global sensitivity measures: By connecting the recent theoretical advances in variance-based and derivative-based global sensitivity analysis to active subspaces, we provide a complete analysis of the links between total Sobol' indices, DGSMs, and activity scores for a broad class of probability distributions. Our analysis is concluded by a result quantifying approximation errors incurred due to fixing unimportant parameters, deemed so by computing their activity scores.
This article is organized as follows. In section 2, a brief theoretical background on the active subspace methodology is provided. In section 3, it is shown that the activity scores provide a reasonable approximation to the DGSMs especially in a high-dimensional setting. Additionally, a relationship between the three global sensitivity measures, namely, the activity scores, DGSMs, and the total Sobol' indices is established. In section 4, a systematic framework for computing the active subspace is provided. Numerical results based on the perturbation-based approach are compared with those obtained using the regression-based approach. The active subspace is initially computed for a 19-dimensional H 2 /O 2 reaction kinetics problem wherein only the A i 's are considered as uncertain. In this case, the two approaches are observed to yield consistent results. We further compute the active for a 33-dimensional H 2 /O 2 reaction kinetics problem in section 5. Once again, active subspace is found to be 1-dimensional. However, in this case, it is observed that although the regression-based estimation of the gradient results in a reasonable approximation of the mean and mode, the variability in the ignition delay due to uncertain rate-controlling parameters is under-estimated. Additionally, the sensitivity towards one of the most important parameter is not captured. Hence, a trade-off between efficiency and accuracy between the two approaches is apparent in the higher-dimensional problem. Finally, a summary and discussion based on our findings is included in section 6.

Active subspaces
Herein, we use a random vector ξ ∈ Ω ∈ R Np to parameterize model uncertainties, where N p is the number of uncertain inputs. In practical computations, the canonical variables ξ i , i = 1, . . . , N p , are mapped to physical ranges meaningful in a given mathematical model. As mentioned in the introduction, an active subspace is a low-dimensional subspace that consists of important directions in a model's input parameter space [12]. The effective variability in a model output f due to uncertain inputs is predominantly captured along these directions. The directions constituting the active subspace are the dominant eigenvectors of the positive semidefinite matrix with µ(dξ) = π(ξ)dξ, where π(ξ) is the joint probability distribution of ξ. Herein, f is assumed to be a square integrable function with continuous partial derivatives with respect to the input parameters; moreover, we assume the partial derivatives are square integrable. Since C is symmetric and positive semidefinite, it admits a spectral decomposition: Here Λ = diag(λ 1 , . . . , λ Np ) with the eigenvalues λ i 's sorted in descending order and W has the (orthonormal) eigenvectors w 1 , . . . , w Np as its columns. The eigenpairs are partitioned about the rth eigenvalue such that λ r /λ r+1 1, The columns of W 1 = w r · · · w r span the dominant eigenspace of C and define the active subspace, and Λ 1 is a diagonal matrix with the corresponding set of eigenvalues, λ 1 , . . . , λ r , on its diagonal. Once the active subspace is computed, dimension reduction is accomplished by transforming the parameter vector ξ into y = W 1 ξ ∈ R r . The elements of y are referred to as the set of active variables. Consider the function Following [12], we use the approximation That is, the model output f (ξ), in the original parameter space, is approximated by G(W 1 ξ) in the active subspace. We could confine uncertainty analysis to the inputs in the active subspace whose dimension is typically much smaller (in applications that admit such a subspace) than the dimension of the original input parameter. To further expedite uncertainty analysis, one could fit a regression surface to G using the following sequence of steps, as outlined in [12, chapter 4].
1. Consider a given set of N data points, In practice, the matrix C defined in (4) is approximated using pseudo-random sampling techniques such as Monte Carlo or Latin hypercube sampling: Clearly the computational effort associated with constructing the matrixĈ scales with the number of samples, N . Hence, an iterative computational approach is adopted in this work to gradually increase N until the dominant eigenpairs are approximated with sufficient accuracy; see Section 4.

GSA measures and their links with active subspaces
Consider a function f = f (ξ 1 , ξ 2 , . . . , ξ Np ). While the active subspace framework described above does not make any assumptions about statistical independence of the inputs ξ i , i = 1, . . . , N p , the classical framework of variance based sensitivity analysis [21,30] assumes that the inputs are statistically independent. While extensions to the cases of correlated inputs exist [31][32][33][34], we limit the discussion in this section to the case of random inputs that are statistically independent and are either uniformly distributed or distributed according to the Boltzmann probability distribution.
Note that a measure µ on R is referred to as a Boltzmann measure if it is absolutely continuous with respect to the Lebesgue measure and admits a density of the form π( where V is a continuous function and C a normalization constant [27]. The total Sobol' index (T i (f )) of a model output, f (ξ) quantifies the total contribution of the input, ξ i to the variance of the output [21]. Mathematically, this can be expressed as follows: where ξ ∼i is the input parameter vector with the i th entry removed, and V denotes the variance. The total Sobol' index accounts for the contribution of a given input to the variability in the output by itself as well as due to its interaction or coupling with other inputs. Determining accurate estimates of T i (f ) typically involves a large number of model runs and is therefore can be prohibitive in the case of compute-intensive applications. Derivative based global sensitivity measures (DGSMs) [26] provide a means for approximating informative upper bounds on T i (f ) at a lower cost; see also [35]. For f : Ω → R, we consider the DGSMs, Here π is the joint PDF of ξ. Let C be as in (4), and consider the spectral decomposition written as C = Np k=1 λ k w k w k . Herein, we use the notation ·, · for the Euclidean inner product. The following result provides a representation of DGSMs in terms of the spectral representation of C: In the case where the eigenvalues decay rapidly to zero, we can obtain accurate approximations of ν i (f ) by truncating the summation: The quantities ν i,r (f ) are called activity scores in [28,29], where links between GSA measures and active subspaces is explored. The following result, which can also be found in [28,29], quantifies the error in this approximation. We provide a short proof for completeness.
The last inequality holds because 1 = e i The utility of this result is realized in problems with high-dimensional parameters in which the eigenvalues λ i , i = 1, . . . , N p , decay rapidly to zero; in such cases, this result implies that where r is the numerical rank of C. This will be especially effective if there is a large gap in the eigenvalues.
The relations recorded in the following lemma will be useful in the discussion that follows.
It was shown in [27] that the total Sobol' index T i (f ) can be bounded in terms of ν i (f ): where for each i, C i is an appropriate Poincaré constant that depends on the distribution of ξ i . For instance, if ξ i is uniformly distributed on [−1, 1], then C i = 4/π 2 ; and in the case ξ i is normally distributed with variance σ 2 i , then C i = σ 2 i . Note that (9) for the special cases of uniformly distributed or normally distributed inputs was established first in [26]. The bound (9) provides a strong theoretical basis for using DGSMs to identify unimportant inputs.
Combining Proposition 3.1 and (9), shows an interesting link between the activity scores and total Sobol' indices. Specifically, by computing the activity scores, we can identify the unimportant inputs. Subsequently, one can attempt to reduce parameter dimension by fixing unimportant inputs at nominal values.
Suppose activity scores are used to approximate DGSMs, and suppose ξ i is deemed unimportant as a result, due to a small activity score. We want to estimate the approximation error that occurs once ξ i is fixed at a nominal value. To formalize this process, we proceed as follows. Let ξ be given and let z be a nominal value for ξ i . Consider the reduced model, obtained by fixing ξ i at the nominal value: and consider the following relative error indicator: .
This error indicator is a function of z with z distributed according to the distribution of ξ i .
where the equality can be shown using arguments similar to the proof of the main result in [36]. Using this, along with (9) and Proposition 3.1, we have In [35] the screening metricν was shown to be useful for detecting unimportant inputs. We can also bound the normalized DGSMs using activity scores as follows. It is straightforward to see that In the case where where λ r+1 ≈ 0, this motivates definition of normalized activity scoresν Remark 3.1. If the random inputs ξ i , i = 1, . . . , N p , are iid, then the C i 's in the definition of the normalized screening metric will cancel and The expression for the denominator follows from Lemma 3.2(b). Also, in the iid case, using Lemma 3.2(a) we can simplify the normalized activity scores as follows.
The significance of the developments in this section are as follows. Theorem 3.1 provides a theoretical basis for parameter dimension reduction using activity scores. This is done by providing an estimate of the error between the reduced model and the original model. If a precise ranking of parameter importance based on total Sobol' indices is desired, one can first identify unimportant inputs by computing activity scores, and then perform a detailed variance based GSA of the remaining model parameters. This approach will provide great computational savings as variance based GSA will now be performed only for a small number of inputs deemed important based on their activity scores. Moreover, the presented result covers a broad class of input distributions coming from the Boltzmann family of distributions. Additionally, the normalized activity scores discussed above provide practical screening metrics that require only computing the activity scores. This is in contrast to the bound in Theorem 3.1 that requires the variance V(f ) of the model output.

Methodology
In this section, we present an efficient framework for computing the active subspace. The framework is employed to the H 2 /O 2 reaction kinetics problem whereby the pre-exponent (A i ) in the rate law associated with individual reactions provided in Table 1 is considered to be uniformly distributed in the interval, [0.9A * i , 1.1A * i ]; A * i is the nominal estimate provided in [13]. Two approaches are explored for estimating the gradient of ignition delay with respect to uncertain A i 's: a perturbationbased approach that involves computation of model gradients using finite difference in order to construct the matrixĈ in (7), and a regression-based approach that involves a linear regression fit to the available set of model evaluations in order to approximate the gradient. The active subspace is computed in an iterative manner to avoid unnecessary model evaluations once converged is established.
As discussed earlier, gradient estimation using finite differences requires additional model evaluations at the neighboring points in the input domain. Hence, for N samples in a d-dimensional parameter space, N (d + 1) model evaluations are needed. On the other hand, gradient estimation using the regression-based approach involves a series of linear regression fits to subsets of available evaluations as discussed in [12] (see Algorithm 1.2). Hence, the computational effort is reduced by a factor (d + 1) when using the regression-based approach. The specific sequence of steps for computing the active subspace is discussed as follows.
We begin by evaluating the gradient of the model output, ∇ ξ f , at an initial set of n 0 samples (generated using Monte Carlo sampling) denoted by ξ i , i = 1, . . . , n 0 . Using the gradient evaluations, the matrix,Ĉ is computed. Eigenvalue decomposition ofĈ yields an initial estimate of the dominant eigenspace,Ŵ 1 and the set of corresponding eigenvalues,Λ 1 . Note thatŴ 1 is obtained by partitioning the eigenspace around λ j such that the ratio of subsequent eigenvalues, λ j λ j+1 > O(10 2 ). At each subsequent iteration, model evaluations are generated at a new set of n k samples. The new set of gradient evaluations are augmented with the available set to re-construct C followed by its eigenvalue decomposition. The relative change in the norm of the difference in squared value of individual components of the dominant eigenvectors between subsequent iterations is evaluated. The process is terminated and the resulting eigenspace is considered to have converged once the maximum relative change at iteration k, δŴ (k) 1,j (j is used as an index for the eigenvectors), is smaller than a given tolerance, τ . A regression fit to G(Ŵ 1 ξ) is used as a surrogate to characterize and quantify the uncertainty in the model output. Moreover, the components of the dominant eigenvectors are used to compute the activity scores, ν r (f ), which provide an insight into the relative importance of the uncertain inputs. Note that the index, r, corresponds to the number of dominant eigenvectors inŴ 1 . The sequence of steps as discussed are outlined in Algorithm 1.

4:
Set N total = n k

11:
Set N total = N total + n k

17:
if max δŴ To assess its feasibility and suitability, we implement Algorithm 1 to compute the active subspace for the 19-dimensional H 2 /O 2 reaction kinetics problem with uncertain A i 's. For the purpose of verification,Ĉ was initially constructed using a large set of samples (N = 1000) in the input domain. The gradient was estimated using finite difference, and hence, a total of 20,000 model runs were performed. In Figure 1, we illustrate the comparison of the resulting normalized eigenvalue spectrum (λ i /λ 0 , i = 1, . . . , 19) with the same using a much smaller set of samples, n = {20, 40, 80, 120}. We observe that the dominant eigenvalues, λ 1 , . . . , λ 4 , are approximated reasonably well with just 20 samples. As expected, the accuracy of higher-order eigenvalues is observed to improve with the sample size. Since λ 1 is roughly three orders of magnitude larger than λ 2 , we expect a 1-dimensional active subspace. Hence, the first eigenvector sufficiently captures the uncertainty in the model output. To further confirm this, we evaluate a relative L 2 norm of the difference (ε N −n L 2 ) between the squared value of corresponding components of the dominant eigenvector, computed using N = 1000 (w 1,N ) and n = {20, 40, 80, 120} (w 1,n ) as follows: The quantity, ε N −n L 2 , was found to be O(10 −2 ) for n = 20, 40, and 80; and O(10 −3 ) for n = 120. Thus, even a small sample size, n = 20, seems to approximate the dominant eigenspace with reasonable accuracy in this case. The iterative strategy therefore offers a significant potential for computational gains.
The active subspace for the 19-dimensional problem was also computed using regression-based estimates of the gradient that does not require model evaluations at neighboring points as discussed inant eigenvector from the two approaches are observed to be in excellent agreement with each other. Moreover, as expected, max(δŴ (r) 1,j ) is observed to decrease with iterations, and requires model evaluations at only 30 samples to converge using τ = 0.1.
As mentioned earlier, the model output f (ξ) i.e. the ignition delay in the H 2 /O 2 reaction in this case, varies predominantly in the active subspace. Hence, f (ξ) can be approximated as G(Ŵ 1 ξ) in the subspace. The plot of G versusŴ 1 ξ, regarded as the sufficient summary plot (SSP), obtained using the perturbation-based and regression-based gradient estimates are compared in Figure 3.
The two SSPs are observed to be in excellent agreement with each other. Moreover, it is interesting Figure 3: An illustrative comparison of the SSPs generated using the perturbation-based and the regression-based strategies for computing the active subspace.
to note that the response in ignition delay based on the considered probability distributions for A i 's is observed to be approximately linear in the active subspace. This observation further indicates that the active subspace is 1-dimensional in this case.
We further estimate the normalized activity scores for individual uncertain inputs (ν i,r ; r=1 since the active subspace is 1-dimensional) using the components of the dominant eigenvector as shown in Algorithm 1. The activity scores for the 19 uncertain pre-exponents (A i 's), estimated using the perturbation-based and regression-based approaches are plotted in Figure 4. The activity scores based on the two approaches for gradient estimation agree favorably with each other as well as those based on the screening metric involving the DGSMs in [35]. It is observed that the uncertainty associated with the ignition delay is largely due to the uncertainty in A 9 and A 1 . Sensitivity towards A 15 and A 17 is also found to be significant.
The above comparisons indicate that the gradient of the ignition delay with respect to the uncertain A i 's is reasonably approximated using both, perturbation-based and regression-based approaches in this case. Since the number of model evaluations required in the case of regression are much fewer, this approach seems favorable for the 19-dimensional problem. In the following section, we shift our focus to the higher-dimensional H 2 /O 2 reaction kinetics application involving 33 uncertain inputs.

H 2 /O 2 reaction kinetics: higher-dimensional case
For the high-dimensional case, we aim to investigate the impact of the uncertain pre-exponents (A i 's) as well as the activation energies (E a,i 's) on ignition delay during the H 2 /O 2 reaction. Prior distributions for A i 's were considered to be the same as provided earlier in section 4. The E a,i 's for all reactions except R 6 -R 9 and R 13 (due to zero nominal values for E a ) were considered to be uncertain and uniformly distributed in the interval: [0.99E * a,i , 1.01E * a,i ]. Hence, the total number of uncertain rate parameters is 33. The nominal values, E * a,i corresponding to the different reaction rates are provided in [13].

Computing the active subspace
The active subspace was computed using the iterative procedure outlined in Algorithm 1. The convergence of the eigenvectors was examined by tracking max(δŴ  In Figure 6 (left), we illustrate converged estimates of the eigenvalues. The second eigenvalue is found to be roughly two orders of magnitude smaller than the first, indicating that the active subspace is 1-dimensional. This is further confirmed by the SSP plots in Figure 6 (right), generated using both, perturbation-based and regression-based approaches for estimating the gradient. Specifically in the case of perturbation-based approach, the SSP is observed to exhibit a linear variation that is reasonably captured by the corresponding straight-line fit,G, constructed using the sequence of steps outlined in Section 2. On the other hand, the SSP obtained by estimating the gradient using regression, is relatively more scattered. Consequently, the two surrogates (straightline fits) are observed to exhibit a discrepancy unlike our earlier observation for the 19-dimensional problem in Figure 3.

Surrogate Assessment
The 1-dimensional surrogate (G) shown in Figure 6 (right) is investigated for its accuracy as well as the ability to capture the uncertainty in the model output in two ways. Firstly, we estimate the relative L 2 norm of the discrepancy (ε d ) between estimates of ignition delay in the case of H 2 /O 2 reaction, obtained using the model output and the surrogate in the following equation: Model evaluations (G) at 10 4 samples in the input domain were used to evaluate ε d . It was estimated to be 1.35×10 −2 and 1.02×10 −1 usingG based on the perturbation-based and the regression-based approaches respectively. Hence, ε d was found to increase by an order of magnitude when the gradient is estimated using the regression-based approach. Therefore, in the norm-sense, it can be said that the 1-dimensional surrogate in the active subspace, obtained using the perturbation-based approach is remarkably accurate with a relative error of less than 1.5%, whereas, the regressionbased approach resulted in a relative error of about 10%. Secondly, we verified the accuracy of the two surrogates in a probabilistic setting. In particular, we compared probability density functions (PDFs) obtained using the true set of model evaluations, and 1-dimensional surrogates (G's) based on the two approaches for gradient estimation as shown in Figure 7. Note that the three PDFs were evaluated using the same set of 10 4 samples in the crossvalidation set. The PDFs based on the model evaluations and the perturbation-based approach are nearly identical, confirming the existence of a 1-dimensional active subspace as well as the accuracy of gradient estimates in this case. On the other hand, although the PDF obtained using the regression-based approach is observed to show consistency with the true PDF in the modal estimate, it clearly underestimates the spread in the ignition delay. Specific values of the mean and the standard deviation in each case are provided in Table 2. The three PDFs are found to be  Table 2: The mean (µ), and the standard deviation (σ), computed using the model (G), and the surrogate (G) based on the two strategies at 10 4 samples in the cross-validation set.
consistent in their prediction of mean value of the ignition delay. The standard deviation estimates using model evaluations and the perturbation-based approach are found to be in close agreement, whereas, it is accurately estimated in the case of regression-based approach only upto the first significant digit.

GSA consistency check
The normalized activity scores (ν i,r ) based on the dominant eigenspace, obtained using the two approaches for gradient estimation (perturbation-based and regression-based), are compared with the total Sobol' indices in Figure 8. Note that the Sobol' indices were computed using the verified 1dimensional surrogate (G) in the active subspace, obtained using the perturbation-based approach. Several useful inferences can be drawn. Firstly, the total Sobol indices are found to be consistent with the normalized activity scores as obtained when the gradient is estimated using a perturbation- based approach. The regression-based approach, however, yields consistent results for the A i 's but fails capture the sensitivity associated with one of the most important parameters, E a, 15 . This is because numerical errors are incurred when approximating the model gradients using a linear regression fit in this case. Secondly, our findings indicate that the ignition delay is predominantly sensitive towards A 1 , A 9 , E a,1 , and E a,15 and moderately sensitive towards A 15 and A 17 . Sensitivity towards the remaining uncertain inputs is found to be low or negligible.

Summary and Discussion
In this work, we focused on the uncertainty associated with the rate-controlling parameters in the H 2 /O 2 reaction mechanism and its impact on ignition delay predictions. The mechanism involves 19 different reactions and in each case, the reaction rate depends upon the choice of a pre-exponent and an activation energy. Hence, in theory, the evolution of the chemical system depends upon 38 inputs. However, we considered epistemic uncertainty in all pre-exponents and activation energies with non-zero nominal values i.e. a total of 33 parameters instead of 38. To facilitate efficient uncertainty analysis, we focused our efforts on reducing the dimensionality of the problem by identifying important directions in the parameter space such that the model output predominantly varies along these directions. These important directions constitute the active subspace. Additionally, we demonstrated that the activity scores, computed using the components of the dominant eigenvectors provide an efficient means for approximating derivative based global sensitivity measures (DGSMs). Furthermore, we established generalized mathematical linkages between the different global sensitivity measures: activity scores, DGSMs, and total Sobol' index which could be exploited to reduce computational effort associated with global sensitivity analysis. Active subspace computation requires repeated evaluations of the gradient of the QoI i.e. the ignition delay. For this purpose, we explored two approaches, namely, perturbation-based (e.g. finite difference -used in this work, adjoint-based methods, etc.) and regression-based. Both approaches were shown to yield consistent results for the 19-dimensional problem wherein only the pre-exponents were considered to be uncertain. Additionally, the activity scores were also shown to be consistent with the screening metric estimates based on DGSMs in [35]. An iterative procedure was adopted to enhance the computational efficiency.
The active subspace was further computed for a 33-dimensional problem wherein all preexponents and activation energies with non-zero nominal estimates were considered uncertain. Once again, the active subspace was found to be 1-dimensional. However, unlike the 19-dimensional case, gradient estimation using the regression-based approach resulted in an under-estimation of the variability in the ignition delay and also failed to capture the sensitivity towards one of the most important parameters, E a, 15 . We attribute this discrepancy to two potential sources of numerical errors: (1) Regression-based approximation of the gradient of the model output, and (2) Scatter in the corresponding SSP (see Figure 6) leading to numerical error incurred by the linear regression fit in this case. Global sensitivity estimates using the perturbation-based approximation of the gradient indicated that the variability in the ignition delay is largely dependent on A 1 , A 9 , E a,1 and E a,15 , although contributions from A 15 and A 17 were also found to be significant. Hence, while the regression-based approach reduces computational effort, the present investigations clearly reveal potential pitfalls associated with this approach.
Therefore, in general, perturbation-based approaches should be preferred for active subspace computation. The regression-based approach can be explored in situations involving intensive simulations where the perturbation-based methods are intractable. We also mention that alternate regression based approaches such as ones based on computing a global quadratic model have been proposed and used in the literature; see e.g., [37]. The applicability of such an approach in the context of high-dimensional chemical reaction networks is subject to future work.
The computational framework presented in this work is agnostic to the choice of the chemical system and can be easily adapted for other systems as long as the quantity of interest is continuously differentiable in the considered domain of the inputs. We have demonstrated that the active subspace could be exploited for efficient forward propagation of the uncertainty from inputs to the output. The resulting activity scores and the low-dimensional surrogate could further guide optimal allocation of computational resources for calibration of the important rate-controlling parameters in a Bayesian setting. Additionally, dimension reduction using active subspaces could assist in developing robust formulations for predicting discrepancy between simulations and measurements due to epistemic uncertainty in the model inputs.