The analysis of the accuracy of results (usually called “responses”) computed by models relies fundamentally on the functional derivatives (usually called “sensitivities”) of the respective model responses with respect to the parameters in the respective computational model. Such sensitivities are needed for many purposes, including: (i) understanding the model by ranking the importance of the various parameters; (ii) performing “reduced-order modeling” by eliminating unimportant parameters and/or processes; (iii) quantifying the uncertainties induced in a model response due to model parameter uncertainties; (iv) performing “model validation,” by comparing computations to experiments to address the question “does the model represent reality?” (v) prioritizing improvements in the model; (vi) performing data assimilation and model calibration as part of forward “predictive modeling” to obtain best-estimate predicted results with reduced predicted uncertainties; (vii) performing inverse “predictive modeling”; (viii) designing and optimizing the system.

Response sensitivities are computed by using either deterministic or statistical methods. The simplest deterministic method for computing response sensitivities is to use finite-difference schemes in conjunction with re-computations using the model with “judiciously chosen” altered parameter values. Evidently, such methods can at best compute approximate values of a very limited number of sensitivities. Deterministic methods that can compute more exactly the values of first-order sensitivities include the “Green’s function method” [1], the “forward sensitivity analysis methodology” [2], and the “direct method” [3], which rely on analytical or numerical differentiation of the computational model under investigation to compute local response sensitivities exactly. However, for a computational model comprising many parameters, the conventional deterministic methods become impractical for computing sensitivities higher than first-order because they are subject to the “curse of dimensionality,” a term coined by Belmann [4] to describe phenomena in which the number of computations increases exponentially in the respective phase-space. In the particular case of sensitivity analysis using conventional deterministic methods, the number of large-scale computations increases exponentially in the phase-space of model parameter as the order of sensitivities increases.

The alternatives to “deterministic methods” are the “statistical methods”, which construct an approximate response distribution (often called “response surface”) in the parameters space, and subsequently use scatter plots, regression, rank transformation, correlations, and so-called “partial correlation analysis,” in order to identify approximate expectation values, variances and covariances for the responses. These statistical quantities are subsequently used to construct quantities that play the role of (approximate) first-order response sensitivities. Thus, statistical methods commence with “uncertainty analysis” and subsequently attempt an approximate “sensitivity analysis” of the approximately computed model response (called a “response surface”) in the phase-space of the parameters under consideration. The currently popular statistical methods for uncertainty and sensitivity analysis are broadly categorized as sampling-based methods [5,6], variance-based methods [7,8], and Bayesian methods [9]. Various variants of the statistical methods for uncertainty and sensitivity analysis are reviewed in the book edited by Saltarelli et al. [10]. The main advantage of using statistical methods for uncertainty and sensitivity analysis is that they are conceptually easy to implement. On the other hand, statistical methods for uncertainty and sensitivity analysis have three major inherent practical drawbacks, as follows:

It is known that the “adjoint method of sensitivity analysis” has been the most efficient method for computing exactly first-order sensitivities, since it requires a single large-scale (adjoint) computation for computing all of the first-order sensitivities, regardless of the number of model parameters. The idea underlying the computation of response sensitivities with respect to model parameters using adjoint operators was first used by Wigner [11] to analyze first-order perturbations in nuclear reactor physics and shielding models based on the linear neutron transport (or diffusion) equation, as subsequently described in textbooks on these subjects [12,13,14,15,16]. Cacuci [2] is credited (see, e.g., [17,18]) for having conceived the rigorous “1^st-order adjoint sensitivity analysis methodology” for generic large-scale nonlinear (as opposed to linearized) systems involving generic operator responses and having introduced these principles to the earth, atmospheric and other sciences.

Cacuci [19,20] has extended his 1^st-order adjoint sensitivity analysis methodology to enable the comprehensive computation of 2^nd-order sensitivities of model responses to model parameters (including imprecisely known domain boundaries and interfaces) for large-scale linear and nonlinear systems. The 2^nd-order adjoint sensitivity analysis methodology for linear systems was applied [21] to compute exactly the 21,976 first-order sensitivities and 482,944,576 second-order sensitivities (of which 241,483,276 are distinct from each other) for an OECD/NEA reactor physics benchmark [22]. This benchmark is modeled by the neutron transport equation involving 21,976 uncertain parameters, the solving of which is representative of “large-scale computations.” The neutron transport equation was solved using the software package PARTISN [23] in conjunction with the MENDF71X cross section library [24], which comprises 618-group cross sections based on ENDF/B-VII.1 nuclear data [25]. The spontaneous fission source was computed using the code SOURCES4C [26]. This work has demonstrated that, contrary to the widely held belief that second- and higher-order sensitivities are negligeable for reactor physics systems, many 2^nd-order sensitivities of the OECD benchmark’s response to the benchmark’s uncertain parameters were much larger than the largest 1^st-order ones. This finding has motivated the investigation of the largest 3^rd-order sensitivities, many of which were found to be even larger than the 2^nd-order ones. Subsequently, the mathematical framework for determining and computing the 4^th-order sensitivities was developed, and many of these were found to be larger than the 3^rd-order ones. This sequence of findings has motivated the development by Cacuci [27] of the “n^th-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Response-Coupled Forward/Adjoint Linear Systems” (abbreviated as “n^th-CASAM-L”), which was developed specifically for linear systems because important model responses produced by such systems are various Lagrangian functionals which depend simultaneously on both the forward and adjoint state functions governing the respective linear system. Among the most important such responses are the Raleigh quotient for computing eigenvalues and/or separation constants when solving partial differential equations, and the Schwinger functional for first-order “normalization-free” solutions [28,29]. These functionals play a fundamental role in optimization and control procedures, derivation of numerical methods for solving equations (differential, integral, integro-differential), etc.

In parallel with developing the n^th-CASAM-L, Cacuci [30] has also developed the n^th-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (n^th-CASAM-N). Just like the n^th-CASAM-L, the n^th-CASAM-N is also formulated in linearly increasing higher-dimensional Hilbert spaces (as opposed to exponentially increasing parameter-dimensional spaces), thus overcoming the curse of dimensionality in sensitivity analysis of nonlinear systems, enabling the most efficient computation of exactly-determined expressions of arbitrarily high-order sensitivities of generic nonlinear system responses with respect to model parameters, uncertain boundaries and internal interfaces in the model’s phase-space.

Recently, Cacuci [31] has introduced the “Second-Order Function/Feature Adjoint Sensitivity Analysis Methodology for Nonlinear Systems” (2^nd-FASAM-N), which enables a considerable reduction (by comparison to the 2^nd-CASAM-N) of the number of large-scale computations needed to compute the second-order sensitivities of a model response with respect to the model parameters, thereby becoming the most efficient methodology known for computing second-order sensitivities exactly. Paralleling the construction of the 2^nd-FASAM-N, this work introduces the “First- and Second-Order Function/Feature Adjoint Sensitivity Analysis Methodology for Response-Coupled Adjoint/Forward Linear Systems” (1^st and 2^nd-FASAM-L). The mathematical methodology of the 1^st-FASAM-L is presented in Section 3, while the mathematical methodology of the 2^nd-FASAM-L is presented in Section 4. The application of the 1^st-FASAM-L and the 2^nd-FASAM-L is illustrated in Section 5 by means of a simplified yet representative energy-dependent neutron-slowing down model which is fundamental importance to reactor physics and design [32,33,34]. The concluding discussion presented in Section 6 prepares the ground for the subsequent generalization of the present work to enable the most efficient possible computation of exact sensitivities of any (arbitrarily-high) order with respect to “feature functions” of model parameters and, hence, to the model’s parameters.

The generic mathematical model considered in this work is fundamentally the same as was considered in [27], but with the major difference that functions (‘features”) of the primary model parameters will be generically identified within the model. The primary model parameters will be denoted as ${\alpha }_1$,…,${\alpha }_{TP}$, where the subscript “TP” indicates “Total number of Primary Parameters;” the qualifier “primary” indicates that these parameters do not depend on any other parameters within the model. These model parameters are considered to include imprecisely known geometrical parameters that characterize the physical system’s boundaries in the phase-space of the model’s independent variables. These boundaries depend on the physical system’s geometrical dimensions, which may be imprecisely known because of manufacturing tolerances. In practice, these primary model parameters are subject to uncertainties. It will be convenient to consider that these parameters are components of a “vector of primary parameters” denoted as $\alpha \triangleq {\left({\alpha }_1,\dots ,{\alpha }_{TP}\right)}^{†}\in {ℝ}^{TP}$, where ${ℝ}^{TP}$ denotes the TP-dimensional subset of the set of real scalars. For subsequent developments, matrices and vectors will be denoted using capital and lower-case bold letters, respectively. The symbol “$\triangleq$” will be used to denote “is defined as” or “is by definition equal to.” Transposition will be indicated by a dagger $(†)$ superscript. The nominal parameter values will be denoted as ${\alpha }^0\triangleq {\left[{\alpha }_1^0,\mathrm{...},{\alpha }_i^0,\mathrm{..},{\alpha }_{TP}^0\right]}^{†}$; the superscript “0” will be used throughout this work to denote “nominal” or “mean” values.

The model is considered to comprise $TI$ independent variables which will be denoted as $x_i,\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TI$, and are considered to be the components of a $TI$-dimensional column vector denoted as $x\triangleq {\left(x_1,\dots ,x_{TI}\right)}^{†}\in {ℝ}^{TI}$, where the sub/superscript “$TI$” denotes the “Total number of Independent variables.” The vector $x\in {ℝ}^{TI}$ of independent variables is considered to be defined on a phase-space domain, denoted as $\Omega \left(\alpha \right)$, $\Omega \left(\alpha \right)\triangleq \left\{-\infty \le {\lambda }_i\left(\alpha \right)\le x_i\le {\omega }_i\left(\alpha \right)\le \infty ;\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TI\right\}$, the boundaries of which may depend on some of the model parameters $\alpha$. The lower boundary-point of an independent variable is denoted as ${\lambda }_i\left(\alpha \right)$ (e.g., the inner radius of a sphere or cylinder, the lower range of an energy-variable, the initial time-value, etc.), while the corresponding upper boundary-point is denoted as ${\omega }_i\left(\alpha \right)$ (e.g., the outer radius of a sphere or cylinder, the upper range of an energy-variable, the final time-value, etc.). A typical example of boundary conditions that depend on imprecisely-known parameters that pertain to the geometry of the model and also to parameters that pertain to the material properties of the respective model occur when modeling particle diffusion within a medium, the boundaries of which are facing vacuum. For such models, the boundary conditions for the respective state (dependent) variables (i.e., particle flux and/or current) are imposed not on the physical boundary but on the “extrapolated boundary” of the respective spatial domain. The “extrapolated boundary” depends both on the imprecisely known physical dimensions of the medium’s domain/extent and also on the medium’s properties, i.e., atomic number densities and microscopic transport cross sections. The boundary of $\Omega \left(\alpha \right)$, which will be denoted as $\partial \Omega \left[\lambda \left(\alpha \right);\omega \left(\alpha \right)\right]$, comprises the set of all of the endpoints ${\lambda }_i\left(\alpha \right),\hspace{0.17em}\hspace{0.17em}{\omega }_i\left(\alpha \right),\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TI,$of the respective intervals on which the components of $x$ are defined, i.e., $\partial \Omega \left[\lambda \left(\alpha \right);\omega \left(\alpha \right)\right]\triangleq \left\{{\lambda }_i\left(\alpha \right)\cup \hspace{0.17em}{\omega }_i\left(\alpha \right),\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TI\right\}$.

The mathematical model that underlies the numerical evaluation of a process and/or state of a physical system comprises equations that relate the system's independent variables and parameters to the system's state/dependent variables. A linear physical system can generally be modeled by a system of coupled operator-equations as follows:

In Equation (1), the vector $\phi \left(x\right)\triangleq {\left[{\phi }_1\left(x\right),\dots ,{\phi }_{TD}\left(x\right)\right]}^{†}$ is a $TD$-dimensional column vector of dependent variables and where the sub/superscript “$TD$” denotes the “Total (number of) Dependent variables.” The functions ${\phi }_i\left(x\right),\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TD$, denote the system’s “dependent variables” (also called “state functions”). The matrix $L\left(x;\alpha \right)\hspace{0.17em}\triangleq \left[L_{ij}\left(x;\alpha \right)\hspace{0.17em}\right],$ $\hspace{0.17em}i,j=1,\mathrm{...},TD$, has dimensions $TD\times TD$. The components $L_{ij}\left(x;\alpha \right)\hspace{0.17em}$ are operators that act linearly on the dependent variables ${\phi }_j\left(x\right)$ and also depend (in general, nonlinearly) on the uncertain model parameters $\alpha$. Furthermore, the vector $g\left(\alpha \right)\triangleq \left[g_1\left(\alpha \right),\dots ,g_{TG}\left(\alpha \right)\right]$ is a $TG$-dimensional vector having components $g_i\left(\alpha \right),\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TG$, which are real-valued functions of (some of) the primary model parameters $\alpha \in {ℝ}^{TP}$. The quantity $TG$ denotes the total number of such functions which appear exclusively in the definition of the model’s underlying equations. Such functions customarily appear in models in the form of correlations that describe “features” of the system under consideration, such as material properties, flow regimes. etc. Usually, the number of functions $g_i\left(\alpha \right)$ is considerably smaller than the total number of model parameters, i.e., $TG\ll TP$. For example, the numerical model (Cacuci and Fang, 2023) of the OECD/NEA reactor physics benchmark (Valentine 2006) comprises 21,976 uncertain primary model parameters (including microscopic cross sections and isotopic number densities) but the neutron transport equation, which is solved to determine the neutron flux distribution within the benchmark, does not use these primary parameters directly but instead uses just several hundreds of “group-averaged macroscopic cross sections” which are functions/features of the microscopic cross sections and isotopic number densities (which in turn are uncertain quantities that would be components of the vector of primary model parameters). In particular, a component $g_j\left(\alpha \right)$ may simply be one of the primary model parameters ${\alpha }_j$, i.e., $g_j\left(\alpha \right)\equiv {\alpha }_j$.

The $TD$-dimensional column vector $Q\left[x;g\left(\alpha \right)\right]\triangleq {\left(q_1,\mathrm{....},q_{TD}\right)}^{†}$, having components $q_i\left[x;g\left(\alpha \right)\right],\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TD$, denotes inhomogeneous source terms, which usually depend nonlinearly on the uncertain parameters $\alpha$. Since the right-side of Equation (1) may contain distributions, the equality in this equation is considered to hold in the weak (i.e., “distributional”) sense. Similarly, all of the equalities that involve differential equations in this work will be considered to hold in the distributional sense.

When $L\left[x;g\left(\alpha \right)\right]$ contains differential operators, a set of boundary and initial conditions which define the domain of $L\left[x;g\left(\alpha \right)\right]$ must also be given. Since the complete mathematical model is considered to be linear in $\phi \left(x\right)$, the boundary and/or initial conditions needed to define the domain of $L\left[x;g\left(\alpha \right)\right]$ must also be linear in $\phi \left(x\right)$. Such linear boundary and initial conditions are represented in the following operator form:

In Equation (2), the quantity $B\left[x;g\left(\alpha \right);\hspace{0.17em}\hspace{0.17em}\lambda \left(\alpha \right);\omega \left(\alpha \right)\right]$ denotes a matrix of dimensions $N_B\times TD$ having components denoted as $B_{ij}\left(x;\alpha \right);\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},N_B;\hspace{0.17em}\hspace{0.17em}j=1,\mathrm{...},TD$, which are operators that act linearly on $\phi \left(x\right)$ and nonlinearly on the components of $g\left(\alpha \right)$; the quantity $N_B$ denotes the total number of boundary and initial conditions. The $N_B$-dimensional column vector $C\left[x;g\left(\alpha \right);\hspace{0.17em}\hspace{0.17em}\lambda \left(\alpha \right);\omega \left(\alpha \right)\right]$ comprises components that are operators which, in general, act nonlinearly on the components of $g\left(\alpha \right)$.

Physical problems modeled by linear systems and/or operators are naturally defined in Hilbert spaces. The dependent variables ${\phi }_i\left(x\right),\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TD$, for the physical system represented by Eqs. (1) and (2) are considered to be square-integrable functions of the independent variables and are considered to belong to a Hilbert space which will be denoted as $H_0\left(\Omega \right)$, where the subscript “zero” denotes “zeroth-level“ or “original.” Higher-level Hilbert spaces, which will be denoted as $H_1\left(\Omega \right)\hspace{0.17em}$ and $H_2\left(\Omega \right)$, will also be used in this work. The Hilbert space $H_0\left(\Omega \right)$ is considered to be endowed with the following inner product, denoted as ${\langle \phi \left(x\right),\psi \left(x\right)\rangle }_0$, between two elements $\phi \left(x\right)\in H_0\left(\Omega \right)$ and $\psi \left(x\right)\in H_0\left(\Omega \right)$:

The “dot” in Equation (3) indicates the “scalar product of two vectors,” which is defined as follows:

The product-notation ${\displaystyle \prod _{i=1}^{TI}{\displaystyle \underset{{\lambda }_i\left(\alpha \right)}{\overset{{\omega }_i\left(\alpha \right)}{\int }}\left[\right]\hspace{0.17em}d{x}_i}}$ in Equation (3) denotes the respective multiple integrals.

The linear operator $L\left[x;g\left(\alpha \right)\right]$ is considered to admit an adjoint operator, which will be denoted as $L^{*}\left[x;g\left(\alpha \right)\right]$ and which is defined through the following relation for a vector $\psi \left(x\right)\in H_0$

In Equation (5), the formal adjoint operator $L^{*}\left[x;g\left(\alpha \right)\right]$ is the $TD\times TD$ matrix comprising elements $L_{ji}^{*}\left[x;g\left(\alpha \right)\right]$ which are obtained by transposing the formal adjoints of the forward operators $L_{ij}\left[x;g\left(\alpha \right)\right]$. Hence, the system adjoint to the linear system represented by (1) and (2) can generally be represented as follows:

When the forward operator $L\left[x;g\left(\alpha \right)\right]$ comprises differential operators, the operations (e.g., integration by parts) that implement the transition from the left-side to the right side of Equation (5) give rise to boundary terms which are collectively called the “bilinear concomitant.” The domain of $L^{*}\left[x;g\left(\alpha \right)\right]$ is determined by selecting adjoint boundary and/or initial conditions so as to ensure that the adjoint system is well-posed mathematically. It is also desirable that the selected adjoint boundary conditions should cause the bilinear concomitant to vanish when implemented in Equation (5) together with the forward boundary conditions given in Equation (2). The adjoint boundary conditions thus selected are represented in operator form by Equation (7).

The relationship shown in Equation (5), which is the basis for defining the adjoint operator, also provides the following fundamental “reciprocity-like” relation between the sources of the forward and the adjoint equations, i.e. Eqs. (1) and (6), respectively:

The functional on the right-side of Equation (8) represents a “detector response”, i.e., a reaction-rate between the particles and the medium represented by $Q^{*}\left[x;g\left(\alpha \right)\right]$ which is equivalent to the “number of counts” of particles incident on a detector of particles that “measures” the particle flux $\phi \left(x\right)$. In view of the relation provided in (8), the vector-valued source term $Q^{*}\left[x;g\left(\alpha \right)\right]\hspace{0.17em}\triangleq {\left\{q_1^{*}\left[x;g\left(\alpha \right)\right],\mathrm{....},q_{TD}^{*}\left[x;g\left(\alpha \right)\right]\right\}}^{†}$ in the adjoint equation Equation (6) is usually associated with the “result of interest” to be measured and/or computed, which is customarily called the system’s “response.” In particular, if $q_i^{*}\left[x;g\left(\alpha \right)\right]=\delta \left(x-x_d\right)$ and $q_{j\ne i}^{*}\left[x;g\left(\alpha \right)\right]=0$, then ${\langle Q^{*}\left[x;g\left(\alpha \right)\right],\hspace{0.17em}\phi \left(x\right)\rangle }_0={\phi }_i\left(x_d\right)$, which means that, in such a case, the right-side of Equation (8) provides the value of the i^th-dependent variable (particle flux, temperature, velocity, etc.) at the point in phase-space where the respective measurement is performed.

The results computed using a mathematical model are customarily called “model responses” (or “system responses” or “objective functions” or “indices of performance”). For linear physical systems, the system’s response may depend not only on the model’s state-functions and on the system parameters but may simultaneously also depend on the adjoint state function. As has been discussed by Cacuci [27,30], any response of a linear system can be formally represented (using expansions or interpolation, if necessary) and fundamentally analyzed in terms of the following generic integral representation: where $S\left[\phi \left(x\right),\psi \left(x\right);g\left(\alpha \right);h\left(\alpha \right);x\right]$ is a suitably differentiable nonlinear function of $\phi \left(x\right),\psi \left(x\right)$, and $\alpha$. The integral representation of the response provided in Equation (9) can represent “averaged” and/or “point-valued” quantities in the phase-space of independent variables. For example, if $R\left[\phi \left(x\right),\psi \left(x\right);f\left(\alpha \right)\right]$ represents the computation or the measurement (which would be a “detector-response”) of a quantity of interest at a point $x_d$ in the phase-space of independent variables, then $S\left[\phi \left(x\right),\psi \left(x\right);g\left(\alpha \right);h\left(\alpha \right);x\right]$ would contain a Dirac-delta functional of the form $\delta \left(x-x_d\right)$. Responses that represent “differentials/derivatives of quantities” would contain derivatives of Dirac-delta functionals in the definition of $S\left[\phi \left(x\right),\psi \left(x\right);g\left(\alpha \right);h\left(\alpha \right);x\right]$. The vector $h\left(\alpha \right)\triangleq \left[h_1\left(\alpha \right),\dots ,h_{TH}\left(\alpha \right)\right]$, having components $h_i\left(\alpha \right),\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TH$, which appears among the arguments of the function $S\left[\phi \left(x\right),\psi \left(x\right);g\left(\alpha \right);h\left(\alpha \right);x\right]$, represents functions of primary parameters that often appear solely in the definition of the response but do not appear in the mathematical definition of the model, i.e., in Eqs. (1), (2), (6) and (7). The quantity $TH$ denotes the total number of such functions which appear exclusively in the definition of the model’s response. Evidently, the response will depend directly and/or indirectly (through the “feature”-functions) on all of the primary model parameters. This fact has been indicated in Equation (9) by using the vector-valued function $f\left(\alpha \right)$ as an argument in the definition of the response $R\left[\phi \left(x\right),\psi \left(x\right);f\left(\alpha \right)\right]$ to represent the concatenation of all of the “features” of the model and response under consideration. The vector $f\left(\alpha \right)$ of “model features” is thus defined as follows:

As defined in Equation (10), the quantity $TF$ denotes the total number of “feature functions of the model’s parameters” which appear in the definition of the nonlinear model’s underlying equations and response.

Solving Eqs. (1) and (2), at the nominal (or mean) values, denoted as ${\alpha }^0\triangleq {\left[{\alpha }_1^0,\mathrm{...},{\alpha }_i^0,\mathrm{..},{\alpha }_{TP}^0\right]}^{†}$, of the model parameters, yields the nominal forward solution, which will be denoted as ${\phi }^0\left(x\right)$. Solving Eqs. (6) and (7) at the nominal values, ${\alpha }^0$, of the model parameters yields the nominal adjoint solution, which will be denoted as ${\psi }^0\left(x\right)$. The nominal value of the response, $R\left[{\phi }^0\left(x\right),{\psi }^0\left(x\right);f\left({\alpha }^0\right)\right]$, is determined by using the nominal parameter values ${\alpha }^0$, the nominal value ${\phi }^0\left(x\right)$ of the forward state function, and the nominal value ${\psi }^0\left(x\right)$ of the adjoint state function.

The definition provided by Equation (9) implies that the model response $R\left[\phi \left(x\right),\psi \left(x\right);f\left(\alpha \right)\right]$ depends on the components of the feature function $f\left(\alpha \right)$, and would therefore admit a Taylor-series expansion around the nominal value $f^0\triangleq f\left({\alpha }^0\right)$, having the following form: where $\delta f_j\triangleq \left[f_j\left(\alpha \right)-f_j^0\right];\hspace{0.17em}\hspace{0.17em}{f}_j^0\triangleq f_j\left({\alpha }^0\right)\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=1,\mathrm{...},TF$. The “sensitivities of the model response with respect to the (feature) functions” are naturally defined as being the functional derivatives of $R\left[f\left(\alpha \right)\right]$ with respect to the components (“features”) $f_j\left(\alpha \right)$ of $f\left(\alpha \right)$. The notation ${\{·\}}_{f^0}$ indicates that the quantity enclosed within the braces is to be evaluated at the nominal values $f^0\triangleq f\left({\alpha }^0\right)$. Since $TF\ll TP$, the computations of the functional derivatives of $R_k\left[f\left(\alpha \right)\right]$ with respect to the functions $f_j\left(\alpha \right)$, which appear in Equation (11), will be considerably less expensive computationally than the computation of the functional derivatives involved in the Taylor-series of the response with respect to the model parameters. The functional derivatives of the response with respect to the parameters can be obtained from the functional derivatives of the response with respect to the “feature” functions $f_j\left(\alpha \right)$ by simply using the chain rule, i.e.: and so on. The evaluation/computation of the functional derivatives $\partial f_{i_1}\left(\alpha \right)/\partial {\alpha }_{j_1}$, ${\partial }^2{f}_{i_1}\left(\alpha \right)/\partial {\alpha }_{j_1}\partial {\alpha }_{j_2}$, etc., does not require computations involving the model, and is therefore computationally trivial by comparison to the evaluation of the functional derivatives (“sensitivities”) of the response with respect to either the functions (“features”) $f_j\left(\alpha \right)$ or the model parameters ${\alpha }_i,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\mathrm{...},TP$.

The range of validity of the Taylor-series shown in Equation (11) is defined by its radius of convergence. The accuracy −as opposed to the “validity”− of the Taylor-series in predicting the value of the response at an arbitrary point in the phase-space of model parameters depends on the order of sensitivities retained in the Taylor-expansion: the higher the respective order, the more accurate the respective response value predicted by the Taylor-series. In the particular cases when the response happens to be a polynomial function of the “feature” functions $f_j\left(\alpha \right)$, the Taylor series represented by Equation (11) is finite and exactly represents the respective model response.

In turn, the functions $f_i\left(\alpha \right)$ can also be formally expanded in a multivariate Taylor-series around the nominal (mean) parameter values ${\alpha }^0$, namely:

The choice of feature functions $f_i\left(\alpha \right)$ is not unique but can be tailored by the user to the problem at hand. The two most important guiding principles for constructing the feature functions $f_i\left(\alpha \right)$ based on the primary parameters are as follows:

The domain of validity of the Taylor-series in Equation (13) is defined by its own radius of convergence. Of course, in the extreme case when no feature function can be constructed, the feature functions will be the primary parameters themselves, in which case the n^th-FASAM-L methodology becomes identical to the previously established n^th-CASAM-L methodology [27].

The “First-Order Function/Feature Adjoint Sensitivity Analysis Methodology for Response-Coupled Forward/Adjoint Linear Systems” (1^st-FASAM-L) aims at enabling the most efficient computation of the first-order sensitivities of a generic model response of the form $R\left[\phi \left(x\right),\psi \left(x\right);\alpha \right]$ with respect to the components of the “features” function $f\left(\alpha \right)$. In preparation for subsequent generalizations towards establishing the generic pattern for computing sensitivities of arbitrarily high-order, the function $u^{\left(1\right)}\left(2;x\right)\triangleq {\left[\phi \left(x\right),\psi \left(x\right)\right]}^{†}$ will be called the “1^st-level forward/adjoint function” and the system of equations satisfied by this function (which is obtained by concatenating the original forward and adjoint equations together with their respective boundary/initial conditions) will be called “the 1^st-Level Forward/Adjoint System (1^st-LFAS)” and will be re-written in the following concatenated matrix-form: where the following definitions were used:

In the list of arguments of the matrix $F^{\left(1\right)}\left[2\times 2;x;f\right]$, the argument “$2\times 2$” indicates that this square matrix comprises four component sub-matrices, as indicated in Equation (16). Similarly, the argument “2” that appears in the block-vectors $\hspace{0.17em}{u}^{\left(1\right)}\left(2;x\right)$, $q_F^{\left(1\right)}\left(2;x;f\right)$, and $\hspace{0.17em}{b}_F^{\left(1\right)}\left[2;u^{\left(1\right)}\left(2;x\right);\hspace{0.17em}f\right]$ indicates that each of these column block-vectors comprises two sub-vectors as components. Also, throughout this work, the quantity “$0$” will be used to denote either a vector or a matrix with zero-valued components, depending on the context. For example, the vector “$0$” in Equation (15) is considered to have as many components as the vector $\hspace{0.17em}{b}_F^{\left(1\right)}\left[u^{\left(1\right)}\left(2;x\right);\hspace{0.17em}f\right]$. On the other hand, the quantity “$0$” which appears in Equation (16) may represent either a (sub) matrix or a vector of the requisite dimensions.

The nominal (or mean) parameter values, ${\alpha }^0$, are considered to be known, but these values will differ from the true values $\alpha$, which are unknown, by variations $\delta \alpha \triangleq {\left(\delta {\alpha }_1,\dots ,\delta {\alpha }_{TP}\right)}^{†}$, where $\delta {\alpha }_i\triangleq {\alpha }_i-{\alpha }_i^0$. The parameter variations $\delta \alpha$ will induce variations $\delta f\left(\alpha \right)\triangleq {\left[\delta f_1\left(\alpha \right),\dots ,\delta f_{TF}\left(\alpha \right)\right]}^{†}$ in the vector-valued function $f\left(\alpha \right)$, around the nominal value $f^0\triangleq f\left({\alpha }^0\right)$, and will also induce variations $\delta \phi \left(x\right)$ and $\delta \psi \left(x\right)$, respectively, around the nominal solution $\left({\phi }^0,{\psi }^0\right)$ through the equations underlying the model. All of these variations will induce variations in the model response $R\left[u^{\left(1\right)}\left(2;x\right);f\right]\equiv R\left[\phi \left(x\right),\psi \left(x\right);f\left(\alpha \right)\right]$, in a neighborhood $\left[{\phi }^0\left(x\right)+\epsilon \delta \phi \left(x\right),{\psi }^0\left(x\right)+\epsilon \delta \psi \left(x\right);f^0+\epsilon \delta f\right]$ around $\left({\phi }^0,{\psi }^0;f^0\right)$, where $\epsilon$ is a real-valued scalar.

Formally, the first-order sensitivities of the response $R\left[u^{\left(1\right)}\left(2;x\right);f\right]$ with respect to the components of the feature function $f\left(\alpha \right)$ are provided by the first-order Gateaux (G-)variation of $R\left(\phi ,\psi ,f\right)$ at the phase-space point $\left({\phi }^0,{\psi }^0,f^0\right)$, which is defined as follows: where the following definitions were used:

In general, the G-variation $\delta R\left({\phi }^0,{\psi }^0,f^0;\delta \phi ,\delta \psi ,\delta f\right)$ is nonlinear in the variations $\delta f\left(\alpha \right)$, $\delta \phi \left(x\right)$ and/or $\delta \psi \left(x\right)$. In such cases, the partial functional Gateaux (G-)derivatives of the response $R\left(\phi ,\psi ,f\right)$ with respect to the functions $\phi ,\psi ,f$ do not exist, which implies that the response sensitivities to the model parameters do not exist, either. Therefore, it will be henceforth assumed in this work that $\delta R\left({\phi }^0,{\psi }^0,f^0;\delta \phi ,\delta \psi ,\delta f\right)$ is linear in the respective variations, so the corresponding partial G-derivatives exist and $\delta R\left({\phi }^0,{\psi }^0,f^0;\delta \phi ,\delta \psi ,\delta f\right)$ is actually the first-order G-differential of the response. The usual numerical methods (e.g., Newton’s method and variants thereof) for solving the equations underlying the model also require the existence of the first-order G-derivatives of the original model equations; these will also be assumed to exist. When the 1^st-order G-derivatives exists, the G-differential $\delta R\left[u^{\left(1,0\right)}\left(2;x\right);f^0;v^{\left(1\right)}\left(2;x\right),\delta f\right]$ can be written as follows:

In Equation (20), the “direct-effect” term ${\left\{\delta R\left[u^{\left(1\right)}\left(2;x\right);f;\delta f\right]\right\}}_{dir}$ comprises only dependencies on $\delta f\left(\alpha \right)$ and is defined as follows:

The following convention/definition was used in Equation (21):

The above convention implies that:

(a) For $j=1,\mathrm{...},TG$:

(b) For $j=TG+1,\mathrm{...},TG+TH$:

(c) For $j=TG+TH+1,\mathrm{...},TG+TH+TI$:

(d) For $j=TG+TH+TI+1,\mathrm{...},TG+TH+2TI$:

The notation on the left-side of Equation (22) represents the inner product between two vectors, but the symbol “($†$)” which indicates “transposition” has been omitted in order to keep the notation as simple as possible. “Daggers” indicating transposition will also be omitted in other inner products, whenever possible, while avoiding ambiguities.

In Equation (20), the “indirect-effect” term ${\left\{\delta R\left[u^{\left(1\right)}\left(2;x\right);f;v^{\left(1\right)}\left(2;x\right)\right]\right\}}_{ind}$ depends only on the variations $v^{\left(1\right)}\left(2;x\right)\triangleq {\left[\delta \phi \left(x\right),\delta \psi \right]}^{†}$ in the state functions, and is defined as follows:

In Eqs. (21) and (27), the notation ${\left\{\right\}}_{{\alpha }^0}$ has been used to indicate that the quantity within the brackets is to be evaluated at the nominal values of the parameters and state functions. This simplified notation is justified by the fact that when the parameters take on their nominal values, it implicitly means that the corresponding state functions also take on their corresponding nominal values. This simplified notation will be used throughout this work.

The direct-effect term can be computed after having solved the forward system modeled by Eqs. (1) and (2), as well as the adjoint system modeled by Eqs. (6) and (7), to obtain the nominal values ${\phi }^0,{\psi }^0$ of the forward and adjoint dependent variables.

On the other hand, the indirect-effect term ${\left\{\delta R\left[u^{\left(1\right)}\left(2;x\right);f;v^{\left(1\right)}\left(2;x\right)\right]\right\}}_{ind}$ defined in Equation (27) can be quantified only after having determined the variations $v^{\left(1\right)}\left(2;x\right)\triangleq {\left[\delta \phi \left(x\right),\delta \psi \right]}^{†}$ in the state functions of the 1^st-Level Forward/Adjoint System (1^st-LFAS). The variations $v^{\left(1\right)}\left(2;x\right)$ are obtained as the solutions of the system of equations obtained by taking the first-order G-differentials of the 1^st-LFAS defined by Eqs. (14) and (15), which are obtained by definition as follows: