Introducing the Second-Order Features Adjoint Sensitivity Analysis Methodology for Neural Integral Equations of Volterra-Type: Mathematical Methodology and Illustrative Application to Nuclear Engineering

1. Introduction

It is well known that Neural Ordinary Differential Equations (NODEs) have enabled the use of deep learning for modeling discretely sampled dynamical systems. NODEs [1,2,3] provide a flexible trade-off between efficiency, memory costs and accuracy while bridging traditional numerical modeling with modern deep learning, as demonstrated by various applications, including time-series, dynamics and control [1,2,3,4,5,6,7,8,9]. However, since each time-step is determined locally in time, NODEs are limited to describing systems that are instantaneous. On the other hand, integral equations (IE) model global “long-distance” spatiotemporal relations, and IE solvers often possess stability properties that are superior to solvers for ordinary and/or partial differential equations. Therefore, differential equations are occasionally recast in integral-equation forms that can be solved more efficiently using IE solvers, as exemplified by the applications described in [10,11,12].

Due to their non-local behavior, IE solvers are suitable for modeling complex dynamics, learning the operator underlying the system under consideration by using data sampled from the respective system. As discussed in [13], the operator learning problem is formulated on finite grids, using finite-difference methods that approximate the domain of the functions under investigation; the learning is performed by using an IE solver which samples the domain of integration continuously. As shown in [14], Neural Integral Equations (NIEs) and the Attentional Neural Integral Equations (ANIEs) can be used to generate dynamics and infer the spatiotemporal relations that initially generated the data, thus enabling the continuous learning of non-local dynamics with arbitrary time resolution. The ANIE interprets the self-attention mechanism as the Nystrom method for approximating integrals [15], which enables efficient integration over higher dimensions, as discussed in [10,11,12,13,14,15] and references therein.

Neural nets are trained by minimizing a “loss functional” chosen by the user to represent the discrepancy between the output produced by the neural net’s decoder and some user-chosen “reference solution.” However, the physical system modeled by a neural net inevitably comprises imperfectly known parameters that stem from measurements and/or computations and are therefore afflicted by uncertainties that stem from the respective experiments and/or computations. Hence, even if the neural net reproduces perfectly a given state of a physical system, the neural net’s “optimized weights” are subject to the uncertainties inherent in the parameters that characterize the underlying physical system, and these uncertainties inevitably propagate to the decoder’s output response. It is hence important to quantify the impact of parameters/weights uncertainties on the uncertainties induced in the decoder’s output response. This impact is quantified by the sensitivities of the decoder’s response with respect to the optimized weights/parameters comprised within the neural net.

Neural nets comprise not only scalar-valued weights/parameters but also functions (e.g., correlations) of such scalar model parameters, which can be conveniently called “features of primary model parameters”. Cacuci [16] has developed the “n^th-Order Features Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (n^th-FASAM-N)”, which enables the most efficient computation of the exact expressions of arbitrarily high-order sensitivities of model responses with respect to the model’s “features”. In turn, the sensitivities of the responses with respect to the primary model parameters are determined, analytically and trivially, by applying the “chain rule” to the expressions obtained for the response sensitivities with respect to the model’s “features.” The n^th-FASAM-N [16] has been applied to develop general first- and second-order sensitivity analysis methodologies for NODEs [17] and for Neural Integral Equations of Fredholm-type [18], which enable the computation, with unsurpassed efficacy, of the exact expressions of first and second-order sensitivities of decoder responses with respect to the underlying neural net’s optimized weights.

This work continues the application of the n^th-FASAM-N [16] methodology to develop the “First- and Second-Order Methodologies for Neural Integral Equations of Volterra Type” (acronyms “1^st-FASAM-NIE-V” and, respectively, “2^nd-FASAM-NIE-V”). The 1^st-FASAM-NIE-V methodology, which is presented in Section 2, enables the most efficient computation of exact expressions of all of the first-order sensitivities of NIE decoder responses with respect to all of the optimal values of the NIE-net’s parameters/weights, after the respective NIE-Volterra-net was optimized to represent the underlying physical system. The efficiency of the 1^st-FASAM-NIE-V is illustrated in Section 3 by applying it to perform a comprehensive first-order sensitivity analysis of the well-known model [19,20,21] of neutron slowing down in a homogeneous medium containing fissionable material.

The general mathematical framework of the 2^nd-FASAM-NIE-Volterra methodology, which is presented in Section 4, enables the most efficient computation of the exact expressions of the second-order sensitivities of NIE decoder responses with respect to all of the optimal values of the NIE-net’s parameters/weights. The efficiency of the 2^nd-FASAM-NIE-V is illustrated in Section 5 by applying it to perform a comprehensive second-order sensitivity analysis of the neutron slowing down model [19,20,21] considered in Section 3. Section 6 concludes this work by presenting a discussion that highlights the unparalleled efficiency of the 2^nd-FASAM-NIE-V methodology for performing sensitivity analysis of Volterra-type neural integral equations.

2. First-Order Features Adjoint Sensitivity Analysis Methodology for Neural Integral Equations of Volterra-Type (1^st-FASAM-NIE-V)

Following [14], a network of nonlinear “Neural Integral Equations of Volterra-type (NIE-Volterra)” can be represented by the system of coupled equations shown below:

$$h_i\left(t\right)=g_i\left[F\left(\theta \right);t\right]+{\phi }_i\left[F\left(\theta \right);t\right]{\displaystyle \underset{t_0}{\overset{t}{\int }}{\psi }_i\left[h\left(\tau \right);F\left(\theta \right);\tau \right]d\tau };t_0\le t\le t_f;1=1,\mathrm{...}TH.$$

(1)

The quantities appearing in Eq. (1) are defined as follows:

(i)

The real-valued scalar quantities $t$, $t_0\le t\le t_f$, and $\tau$, $t_0\le \tau \le t_f$, are time-like independent variables which parameterize the dynamics of the hidden/latent neuron units. Customarily, the variable $t$ is called the “global time” while the variable $\tau$ is called the “local time.” The initial time-value is denoted as $t_0$ while the stopping time-value is denoted as $t_f$.

(ii)

The components of the vector $\theta \triangleq {\left[{\theta }_1,\mathrm{...},{\theta }_{TW}\right]}^{†}$ represent scalar learnable adjustable weights, where $TW$ denotes the total number of adjustable weights in all of the latent neural nets. The components of the column-vector $\theta \triangleq {\left[{\theta }_1,\mathrm{...},{\theta }_{TW}\right]}^{†}$ are considered to be “primary parameters” while the components of the vector-valued function $F\left(\theta \right)\triangleq {\left[F_1\left(\theta \right),\mathrm{...},F_{TF}\left(\theta \right)\right]}^{†}$ represent the ”feature” functions of the respective weights. The quantity $TF$ denotes the “total number of feature/functions of the primary model parameters” comprised in the NIE-Volterra. In general, $F\left(\theta \right)$ is a nonlinear function of $\theta$. The total number of feature functions must necessarily be smaller than the total number of primary parameters (weights), i.e., $TF<TW$. In the extreme case when there are no feature functions, it follows that $F_i\left(\theta \right)\equiv {\theta }_i$, for all $i=1,\mathrm{...},TW\equiv TF$. In this work, all vectors are considered to be column vectors, and the dagger “$†$” symbol will be used to denote “transposition.”. The symbol “$\triangleq$” will be used to denote “is defined as” or, equivalently, “is by definition equal to.”

(iii)

The $TH$-dimensional vector-valued function $h\left(t\right)\triangleq {\left[h_1\left(t\right),\mathrm{...},h_{TH}\left(t\right)\right]}^{†}$ represents the hidden/latent neural networks. The quantity $TH$ denotes the total number components of $h\left(t\right)$. At the initial time-value $t_0$, the functions $h_i\left(t_0\right)$ take on the known values $h_i\left(t_0\right)=g_i\left[F\left(\theta \right);t_0\right]$.

(iv)

The functions $g_i\left[F\left(\theta \right);t\right]=h_i\left(t_0\right)$, $i=1,\mathrm{...},TH$, model the initial state (“encoder”) of the network. The functions ${\phi }_i\left[F\left(\theta \right);t\right]$ and ${\psi }_i\left[h\left(\tau \right);F\left(\theta \right);\tau \right]$, $i=1,\mathrm{...},TH$, depend nonlinearly on $h\left(t\right)$ and $F\left(\theta \right)$, respectively, and model the dynamics of the latent neurons.

The “training” of the NIE-Volterra net is accomplished by using the “adjoint” or other methods to minimize the user-chosen “loss functional” intended to represent the discrepancy between the output produced by the NIE-decoder and a “reference solution” chosen by the user. After the training is completed, the primary parameters (“weights”) $\theta \triangleq {\left[{\theta }_1,\mathrm{...},{\theta }_{TW}\right]}^{†}$ will have been assigned “optimal” values which are obtained as a result of having minimized the chosen loss functional. These optimal values for the primary parameters (“weights”) will be denoted using a superscript “zero,” as follows: ${\theta }^0\triangleq {\left[{\theta }_1^0,\mathrm{...},{\theta }_{TW}^0\right]}^{†}$. Using these optimal/nominal parameter values to solve the NIE-system will yield the optimal/nominal solution $h^0\left(t\right)\triangleq {\left[h_1^0\left(t\right),\mathrm{...},h_{TH}^0\left(t\right)\right]}^{†}$, which will satisfy the following form of Eq. (1):

$$h_i^0\left(t\right)=g_i^0\left[F\left({\theta }^0\right);t\right]+{\phi }_i^0\left[F\left({\theta }^0\right);t\right]{\displaystyle \underset{t_0}{\overset{t}{\int }}{\psi }_i^0\left[h^0\left(\tau \right);F\left({\theta }^0\right);\tau \right]d\tau };t_0\le t\le t_f;1=1,\mathrm{...}TH.$$

(2)

After the NIE-net is optimized to reproduce the underlying physical system as closely as possible, the subsequent responses of interest are no longer “loss functions” but become specific functionals of NIE’s “decoder” response/output. Such a decoder-response, which will be denoted as $R\left[h;F\left(\theta \right)\right]$, can be generally represented a scalar-valued functional of $h\left(t\right)$ and $F\left(\theta \right)$, defined as follows:

$$R\left[h;F\left(\theta \right)\right]={\displaystyle \underset{t_0}{\overset{t_f}{\int }}D\left[h\left(t\right);F\left(\theta \right);t\right]dt}$$

(3)

The function $D\left[h\left(t\right);F\left(\theta \right);t\right]$ models the decoder and may contain distributions (e.g., Dirac-delta and/or Heaviside functionals, etc.), if the decoder-response is to be evaluated at some particular point in time or over a subinterval within the interval $\left[t_0,t_f\right]$.

The optimal value of the decoder-response, denoted as $R\left[h^0;F\left({\theta }^0\right)\right]$, is represented by evaluating Eq. (3) at the optimal/nominal parameter values ${\theta }^0\triangleq {\left[{\theta }_1^0,\mathrm{...},{\theta }_{TW}^0\right]}^{†}$ and optimal/nominal solution $h^0\left(t\right)$, as follows:

$$R\left[h^0;F\left({\theta }^0\right)\right]={\displaystyle \underset{t_0}{\overset{t_f}{\int }}D\left[h^0\left(t\right);F\left({\theta }^0\right);t\right]dt}$$

(4)

The true values $\theta$ of the primary parameters (“weights”) that characterize the physical system modeled by the NIE-V net are afflicted by uncertainties inherent to the experimental and/or computational methodologies employed to model the original physical system. Therefore, the true values $\theta$ of the primary parameters (“weights”) will differ from the known nominal values ${\theta }^0$ (which are obtained after training the NIE-net to represent the model of the physical system) by variations denoted as $\delta \theta \triangleq \theta -{\theta }^0$. The variations $\delta \theta \triangleq \theta -{\theta }^0$ will induce corresponding variations $\delta F\triangleq F\left(\theta \right)-F^0$, $F^0\triangleq F\left({\theta }^0\right)$, in the feature functions, which in turn will induce variations $v^{\left(1\right)}\left(t\right)\triangleq {\left[v_1^{\left(1\right)}\left(t\right),\dots ,v_{TH}^{\left(1\right)}\left(t\right)\right]}^{†}$, $v_i^{\left(1\right)}\left(t\right)\triangleq h_i\left(t\right)-h_1^0\left(t\right)$, $i=1,\mathrm{...},TH$, around the nominal/optimal functions $h^0\left(t\right)$. Subsequently, the variations $\delta F$ and $v^{\left(1\right)}\left(t;x\right)$ will induce variations in the NIE decoder’s response.

The 1^st-FASAM-IDE-V methodology for computing the first-order sensitivities of the decoder’s response with respect to the NIE’s weights will be established by applying the same principles as those underlying the 1^st-FASAM-N [16] methodology. These first-order sensitivities are embodied in the first-order G-variation $\delta R\left(h^0;F^0;v^{\left(1\right)};\delta F\right)$ of the response $R\left[h;F\left(\theta \right)\right]$, for variations $v^{\left(1\right)}\left(t\right)$ and $\delta F$ around the nominal values $h^0\left(t\right)$ and $F^0$, which is by definition obtained as follows:

$$\begin{array}{l}\delta R\left(h^0;F^0;v^{\left(1\right)};\delta F\right)={\left\{{\displaystyle \frac{d}{d\epsilon }}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}D\left[h^0\left(t\right)+\epsilon v^{\left(1\right)}\left(t\right);F^0+\epsilon \delta F;t\right]dt}\right\}}_{\epsilon =0}\\ ={\left\{\delta R\left(h^0;F^0;\delta F\right)\right\}}_{dir}+{\left\{\delta R\left(h^0;F^0;v^{\left(1\right)}\right)\right\}}_{ind}.\end{array}$$

(5)

In Eq. (5), the “direct-effect term” ${\left\{\delta R\left(h^0;F^0;\delta F\right)\right\}}_{dir}$ arises directly from variations $\delta F$ (which in turn stem from parameter variations $\delta \theta$) and is defined as follows:

$$\begin{array}{l}{\left\{\delta R\left(h^0;F^0;\delta F\right)\right\}}_{dir}\triangleq {\displaystyle \underset{t_0}{\overset{t_f}{\int }}{\left\{\frac{\partial D\left[h\left(t\right);F\left(\theta \right);t\right]}{\partial F}\delta F\right\}}_{\left(h^0;F^0\right)}dt}\\ \equiv {\displaystyle \sum _{j=1}^{TF}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{\left\{\frac{\partial D\left[h\left(t\right);F\left(\theta \right);t\right]}{\partial F_j}\delta F_j\right\}}_{\left(h^0;F^0\right)}dt}},\end{array}$$

(6)

while the “indirect-effect term” ${\left\{\delta R\left(h^0;F^0;v^{\left(1\right)}\right)\right\}}_{ind}$ arises through the variations $v^{\left(1\right)}\left(t\right)$ in the hidden state functions $h\left(t\right)$, and is defined as follows:

$$\begin{array}{l}{\left\{\delta R\left(h^0;F^0;v^{\left(1\right)}\right)\right\}}_{ind}\triangleq {\displaystyle \underset{t_0}{\overset{t_f}{\int }}{\left\{\frac{\partial D\left[h\left(t\right);F\left(\theta \right);t\right]}{\partial h}{v}^{\left(1\right)}\left(t\right)\right\}}_{\left(h^0;F^0\right)}dt}\\ \equiv {\displaystyle \sum _{j=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{\left\{\frac{\partial D\left[h\left(t\right);F\left(\theta \right);t\right]}{\partial h_j}\delta h_j\right\}}_{\left(h^0;F^0\right)}dt}}.\end{array}$$

(7)

The first-order relationship between the variations $v^{\left(1\right)}\left(t\right)$ and $\delta F$ is obtained from the first-order G-variation of Eq. (1) for $i=1,\mathrm{...}TH$, as follows:

$$\begin{array}{l}{\left\{{\displaystyle \frac{d}{d\epsilon }}\left[h_i^0\left(t\right)+\epsilon v_i^{\left(1\right)}\left(t\right)\right]\right\}}_{\epsilon =0}={\left\{{\displaystyle \frac{d}{d\epsilon }}{g}_i\left(F^0+\epsilon \delta F;t\right)\right\}}_{\epsilon =0}+{\left\{{\displaystyle \frac{d}{d\epsilon }}{\phi }_i\left(F^0+\epsilon \delta F;t\right)\right\}}_{\epsilon =0}\\ \times {\left\{{\displaystyle \frac{d}{d\epsilon }}{\displaystyle \underset{t_0}{\overset{t}{\int }}{\psi }_i\left[h_1^0\left(\tau \right)+\epsilon v_1^{\left(1\right)}\left(\tau \right),\mathrm{...},h_{TH}^0\left(\tau \right)+\epsilon v_{TH}^{\left(1\right)}\left(\tau \right);F^0+\epsilon \delta F;\tau \right]d\tau }\right\}}_{\epsilon =0};t_0\le t\le t_f.\end{array}$$

(8)

Performing the operations indicated in Eq. (8) yields the following NIE-V net, which will be called the “1^st-Level Variational Sensitivity System” (1^st-LVSS), for the components $v_i^{\left(1\right)}\left(t\right)$, $i=1,\mathrm{...}TH$, of the “1^st-level variational function” $v^{\left(1\right)}\left(t\right)$:

$$v_i^{\left(1\right)}\left(t\right)={\left\{{\phi }_i\left(F;t\right){\displaystyle \sum _{j=1}^{TH}{\displaystyle \underset{t_0}{\overset{t}{\int }}\frac{\partial {\psi }_i\left[h\left(\tau \right);F;\tau \right]}{\partial h_j\left(\tau \right)}{v}_j^{\left(1\right)}\left(\tau \right)d\tau }}\right\}}_{\left(h^0,F^0\right)}+{\displaystyle \sum _{k=1}^{TF}{\left\{q_{ik}\left(F;t\right)\delta F_k\right\}}_{\left(h^0;F^0\right)}}$$

(9)

where:

$$q_{ik}\left(F;t\right)\triangleq {\displaystyle \frac{\partial g_i\left(F;t\right)}{\partial F_k}}+{\displaystyle \frac{\partial {\phi }_i\left(F;t\right)}{\partial F_k}}{\displaystyle \underset{t_0}{\overset{t}{\int }}{\psi }_i\left[h\left(\tau \right);F;\tau \right]d\tau }+{\phi }_i\left(F;t\right){\displaystyle \underset{t_0}{\overset{t}{\int }}\frac{\partial {\psi }_i\left[h\left(\tau \right);F;\tau \right]}{\partial F_k}d\tau }.$$

(10)

As indicated in Eq. (9), the 1^st-LVSS is to be computed at the nominal/optimal values for the respective model parameters. It is important to note that the 1^st-LVSS is linear in the variational function $v^{\left(1\right)}\left(t\right)$, although it generally remains nonlinear in $h\left(t\right)$.

The 1^st-LVSS would need to be solved anew to obtain the function $v^{\left(1\right)}\left(t\right)$ that would correspond to each variation $\delta F_j$, $j=1,\mathrm{...},TF$; this procedure would become prohibitively expensive computationally if $TF$ is a large number. The need for repeatedly solving the 1^st-LVSS can be avoided by recasting the indirect-effect term ${\left\{\delta R\left(h^0;F^0;v^{\left(1\right)}\right)\right\}}_{ind}$ in terms of an expression that does not involve the function $v^{\left(1\right)}\left(t\right)$. This goal can be achieved by expressing ${\left\{\delta R\left(h^0;F^0;v^{\left(1\right)}\right)\right\}}_{ind}$ in terms of another function, which will be called the “1^st-level adjoint function,” and which will be the solution of the “1^st-Level Adjoint Sensitivity System (1^st-LASS)” to be constructed next.

The 1^st-LASS will be constructed in a Hilbert space, denoted as $H_1\left({\Omega }_t\right)$, where ${\Omega }_t\triangleq \left[t_0,t_f\right]$, comprising elements of the same form as $v^{\left(1\right)}\left(t\right)\in H_1\left({\Omega }_t\right)$. The inner product of two elements ${\chi }^{(1)}\left(t\right)\triangleq {\left[{\chi }_1^{\left(1\right)}\left(t\right),\dots ,{\chi }_{TH}^{\left(1\right)}\left(t\right)\right]}^{†}\in H_1\left({\Omega }_t\right)$ and ${\eta }^{(1)}\left(t\right)\triangleq {\left[{\eta }_1^{\left(1\right)}\left(t\right),\dots ,{\eta }_{TH}^{\left(1\right)}\left(t\right)\right]}^{†}\in H_1\left({\Omega }_t\right)$ will be denoted as ${〈{\chi }^{(1)}\left(t\right),{\eta }^{(1)}\left(t\right)〉}_1$ and is defined as follows:

$${〈{\chi }^{(1)}\left(t\right),{\eta }^{(1)}\left(t\right)〉}_1\triangleq {\displaystyle \underset{t_0}{\overset{t_f}{\int }}{\left[{\chi }^{(1)}\left(t\right)\right]}^{†}{\eta }^{(1)}\left(t\right)dt}={\displaystyle \sum _{j=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{\chi }_i^{\left(1\right)}\left(t\right){\eta }_i^{\left(1\right)}\left(t\right)dt}}$$

(11)

The inner product ${〈{\chi }^{(1)}\left(t\right),{\eta }^{(1)}\left(t\right)〉}_1$ is required to hold in a neighborhood of the nominal values $\left(h^0;F^0\right)$.

The next step is to form the inner product of Eq. (9) with a vector $a^{\left(1\right)}\left(t\right)\triangleq {\left[a_1^{\left(1\right)}\left(t\right),\dots ,a_{TH}^{\left(1\right)}\left(t\right)\right]}^{†}\in H_1\left({\Omega }_t\right)$, where the superscript “(1)” indicates “1^st-level”, to obtain the following relationship:

$$\begin{array}{l}{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(t\right)v_i^{\left(1\right)}\left(t\right)dt}}-{\left\{{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(t\right){\phi }_i\left(F;t\right)dt}}{\displaystyle \sum _{j=1}^{TH}{\displaystyle \underset{t_0}{\overset{t}{\int }}\frac{\partial {\psi }_i\left[h\left(\tau \right);F;\tau \right]}{\partial h_j\left(\tau \right)}{v}_j^{\left(1\right)}\left(\tau \right)d\tau }}\right\}}_{\left(h^0,F^0\right)}\\ ={\displaystyle \sum _{i=1}^{TH}{\displaystyle \sum _{k=1}^{TF}\left[{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(t\right)q_{ik}\left(F;t\right)dt}\right]\delta F_k}}.\end{array}$$

(12)

The second term on the left-side of Eq. (12) is transformed using “integration by parts” as follows:

$$\begin{array}{l}{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(t\right){\phi }_i\left(F;t\right)dt}}{\displaystyle \sum _{j=1}^{TH}{\displaystyle \underset{t_0}{\overset{t}{\int }}\frac{\partial {\psi }_i\left[h\left(\tau \right);F;\tau \right]}{\partial h_j\left(\tau \right)}{v}_j^{\left(1\right)}\left(\tau \right)d\tau }}\\ ={{\displaystyle \sum _{i=1}^{TH}{\displaystyle \sum _{j=1}^{TH}\left\{\left[{\displaystyle \underset{t_0}{\overset{t}{\int }}\frac{\partial {\psi }_i\left[h\left(\tau \right);F;\tau \right]}{\partial h_j\left(\tau \right)}{v}_j^{\left(1\right)}\left(\tau \right)d\tau }\right]\left[{\displaystyle \underset{t_0}{\overset{t}{\int }}{a}_i^{\left(1\right)}\left(\tau \right){\phi }_i\left(F;\tau \right)d\tau }\right]\right\}}}}_{t=t_f}\\ -{{\displaystyle \sum _{i=1}^{TH}{\displaystyle \sum _{j=1}^{TH}\left\{\left[{\displaystyle \underset{t_0}{\overset{t}{\int }}\frac{\partial {\psi }_i\left[h\left(\tau \right);F;\tau \right]}{\partial h_j\left(\tau \right)}{v}_j^{\left(1\right)}\left(\tau \right)d\tau }\right]\left[{\displaystyle \underset{t_0}{\overset{t}{\int }}{a}_i^{\left(1\right)}\left(\tau \right){\phi }_i\left(F;\tau \right)d\tau }\right]\right\}}}}_{t=t_0}\\ -{\displaystyle \sum _{i=1}^{TH}{\displaystyle \sum _{j=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}\left[\frac{\partial {\psi }_i\left[h\left(t\right);F;t\right]}{\partial h_j\left(t\right)}{v}_j^{\left(1\right)}\left(t\right)dt\right]\left[{\displaystyle \underset{t_0}{\overset{t}{\int }}{a}_i^{\left(1\right)}\left(\tau \right){\phi }_i\left(F;\tau \right)d\tau }\right]}}}\\ ={\displaystyle \sum _{i=1}^{TH}{\displaystyle \sum _{j=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{v}_j^{\left(1\right)}\left(t\right)\frac{\partial {\psi }_i\left[h\left(t\right);F;t\right]}{\partial h_j\left(t\right)}dt{\displaystyle \underset{t}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(\tau \right){\phi }_i\left(F;\tau \right)d\tau }}}}.\\ \end{array}$$

(13)

Replacing the result obtained in Eq. (13) into the left-side of Eq. (12) yields the following relation for the left-side of Eq. (12):

$$\begin{array}{l}{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(t\right)v_i^{\left(1\right)}\left(t\right)dt}}-{\left\{{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(t\right){\phi }_i\left(F;t\right)dt}}{\displaystyle \sum _{j=1}^{TH}{\displaystyle \underset{t_0}{\overset{t}{\int }}\frac{\partial {\psi }_i\left[h\left(\tau \right);F;\tau \right]}{\partial h_j\left(\tau \right)}{v}_j^{\left(1\right)}\left(\tau \right)d\tau }}\right\}}_{\left(h^0,F^0\right)}\\ ={\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{v}_i^{\left(1\right)}\left(t\right){\left\{a_i^{\left(1\right)}\left(t\right)-{\displaystyle \sum _{k=1}^{TH}\frac{\partial {\psi }_k\left[h\left(t\right);F;t\right]}{\partial h_i\left(t\right)}{\displaystyle \underset{t}{\overset{t_f}{\int }}{a}_k^{\left(1\right)}\left(\tau \right){\phi }_k\left(F;\tau \right)d\tau }}\right\}}_{\left(h^0,F^0\right)}dt}}.\end{array}$$

(14)

The term on the right-side of Eq. (14) is now required to represent the “indirect-effect” term defined in Eq. (7), which is achieved by requiring that the components of the function $a^{(1)}\left(t\right)$ satisfy the following system of equations for $i=1,\mathrm{...},TH$:

$$a_i^{\left(1\right)}\left(t\right)-{\displaystyle \sum _{k=1}^{TH}\frac{\partial {\psi }_k\left[h\left(t\right);F;t\right]}{\partial h_i\left(t\right)}{\displaystyle \underset{t}{\overset{t_f}{\int }}{a}_k^{\left(1\right)}\left(\tau \right){\phi }_k\left(F;\tau \right)d\tau }}={\displaystyle \frac{\partial D\left[h\left(t\right);F\left(\theta \right);t\right]}{\partial h_i}}.$$

(15)

The Volterra-like neural system obtained in Eq. (15) will be called the “1^st-Level Adjoint Sensitivity System” and its solution, $a^{(1)}\left(t\right)$, will be called the “1^st-level adjoint sensitivity function.” The 1^st-LASS is to be solved using the nominal/optimal values for the parameters and for the function $h\left(t\right)$ but this fact has not been explicitly indicated in order to simplify the notation. The 1^st-LASS is linear in $a^{(1)}\left(t\right)$ but is, in general, nonlinear in $h\left(t;x\right)$. Notably, the 1^st-LASS is independent of any parameter variations and needs to be solved once only to determine the 1^st-level adjoint sensitivity function $a^{(1)}\left(t\right)$. The 1^st-LASS is a “final-value problem” since the computation of the adjoint function $a^{(1)}\left(t\right)$ will commence at $t=t_f$, with the known values $a_i^{\left(1\right)}\left(t_f\right)={\left\{\partial D\left[h\left(t\right);F\left(\theta \right);t\right]/\partial h_i\left(t\right)\right\}}_{t=t_f}.$

It follows from Eqs. (12)‒(15) that the indirect-effect term defined in Eq. (7) can be expressed in terms of the 1^st-level adjoint sensitivity function $a^{(1)}\left(t\right)$ as follows:

$${\left\{\delta R\left(h^0;F^0;a^{(1)}\right)\right\}}_{ind}={\displaystyle \sum _{k=1}^{TF}{\left\{{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(t\right)q_{ik}\left(F;t\right)dt}}\right\}}_{\left(h^0;F^0\right)}}\delta F_k.$$

(16)

Using the results obtained in Eqs. (16) and (6) in Eq. (5) yields the following expression for the G-variation $\delta R\left(h^0;F^0;v^{\left(1\right)};\delta F\right)$, which is seen to be linear in $\delta F$:

$$\begin{array}{l}\delta R\left(h^0;F^0;v^{\left(1\right)};\delta F\right)={\displaystyle \sum _{j=1}^{TF}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{\left\{\frac{\partial D\left[h\left(t\right);F\left(\theta \right);t\right]}{\partial F_j}\delta F_j\right\}}_{\left(h^0;F^0\right)}dt}}\\ +{\left\{{\displaystyle \sum _{i=1}^{TH}{\displaystyle \sum _{j=1}^{TF}\left[{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(t\right)q_{ij}\left(F;t\right)dt}\right]\delta F_j}}\right\}}_{\left(h^0;F^0\right)}\triangleq {\displaystyle \sum _{j=1}^{TF}{\left\{\frac{\partial R}{\partial F_j}\right\}}_{\left(h^0;F^0\right)}}\delta F_j.\end{array}$$

(17)

Identifying in Eq. (17) the expressions that multiply the variations $\delta F_j$, yields the following expressions for the sensitivities $\partial R/\partial F_j$ of the response $R\left[h;F\left(\theta \right)\right]$ with respect to the components $F_j\left(\theta \right)$ of the feature function $F\left(\theta \right)$, for $j=1,\mathrm{...},TF$:

$$\begin{array}{l}{\displaystyle \frac{\partial R\left[h;F\left(\theta \right)\right]}{\partial F_j}}={\displaystyle \underset{t_0}{\overset{t_f}{\int }}\frac{\partial D\left[h\left(t\right);F\left(\theta \right);t\right]}{\partial F_j}dt}+{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(t\right)\frac{\partial g_i\left(F;t\right)}{\partial F_j}dt}}\\ +{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(t\right)\left\{\frac{\partial {\phi }_i\left(F;t\right)}{\partial F_j}{\displaystyle \underset{t_0}{\overset{t}{\int }}{\psi }_i\left[h\left(\tau \right);F;\tau \right]d\tau +{\phi }_i\left(F;t\right){\displaystyle \underset{t_0}{\overset{t}{\int }}\frac{\partial {\psi }_i\left[h\left(\tau \right);F;\tau \right]}{\partial F_j}d\tau }}\right\}dt}}.\end{array}$$

(18)

The expression on the right-side of Eq. (18) is to be evaluated at the nominal/optimal values for the respective model parameters, but this fact has not been indicated explicitly in order to simplify the notation.

The sensitivities with respect to the primary model parameters can be obtained by using the result obtained in Eq. (18) together with the “chain rule” of differentiating compound functions, as follows:

$${\displaystyle \frac{\partial R}{\partial {\theta }_j}}={\displaystyle \sum _{i=1}^{TF}\frac{\partial R}{\partial F_i}\frac{\partial F_i}{\partial {\theta }_j}},j=1,\mathrm{...},TW.$$

(19)

The sensitivities $\partial R/\partial F_j$ are obtained from Eq. (18) while the derivatives $\partial F_i/\partial {\theta }_j$ are obtained analytically, exactly, from the known expressions of the feature functions $F_i\left(\theta \right)$.

Particular Case: The First-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Neural Integral Equations of Volterra-Type (1^st-CASAM-NIE-V)

When no feature functions can be constructed from the model parameters/weights, the feature functions become identical to the parameters, i.e., $F_i\left(\theta \right)\equiv {\theta }_i$ for all $i=1,\mathrm{...},TF\triangleq TW$. In this case, the expression obtained in Eq. (18) yields directly the first-order sensitivities $\partial R/\partial {\theta }_j$ of the decoder response with respect to the model weights/parameters, for all $j=1,\mathrm{...},TW$, taking on the following specific form:

$$\begin{array}{l}{\displaystyle \frac{\partial R\left[h;F\left(\theta \right)\right]}{\partial {\theta }_j}}={\displaystyle \underset{t_0}{\overset{t_f}{\int }}\frac{\partial D\left[h\left(t\right);F\left(\theta \right);t\right]}{\partial {\theta }_j}dt}+{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(t\right)\frac{\partial g_i\left(F;t\right)}{\partial {\theta }_j}dt}}\\ +{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(t\right)\left\{\frac{\partial {\phi }_i\left(F;t\right)}{\partial {\theta }_j}{\displaystyle \underset{t_0}{\overset{t}{\int }}{\psi }_i\left[h\left(\tau \right);F;\tau \right]d\tau +{\phi }_i\left(F;t\right){\displaystyle \underset{t_0}{\overset{t}{\int }}\frac{\partial {\psi }_i\left[h\left(\tau \right);F;\tau \right]}{\partial {\theta }_j}d\tau }}\right\}dt}}.\end{array}$$

(20)

Since the 1^st-LASS is independent of any parameter variations, the 1^st-level adjoint sensitivity function $a^{\left(1\right)}\left(t\right)\triangleq {\left[a_1^{\left(1\right)}\left(t\right),\dots ,a_{TH}^{\left(1\right)}\left(t\right)\right]}^{†}$ which appears in Eq. (20) remains the solution of the 1^st-LASS defined by Eq. (15). In this case, however, all of the sensitivities $\partial R/\partial {\theta }_j$, for all $j=1,\mathrm{...},TW$ would be obtained by computing integrals using quadrature formulas. Thus, when there are no feature functions of parameters, the 1^st-FASAM-NIE-V reduces to the “First-Order Comprehensive Adjoint Sensitivity Analysis Methodology [16] applied to Neural Integral Equations of Volterra-Type” (1^st-CASAM-NIE-V). On the other hand, when features of parameters can be constructed, only $TF$ $\left(TF<TW\right)$ numerical computations of integrals using quadrature formulas are required, using Eq. (18) to obtain the sensitivities $\partial R/\partial F_j$, $j=1,\mathrm{...},TF$. Subsequently, the sensitivities with respect to the model’s weights/parameters are obtained analytically using the chain-rule provided in Eq. (19).

3. Illustrative Application of the 1^st-CASAM-NIE-V and 1^st-FASAM-NIE-V Methodologies to Neutron Slowing Down in an Infinite Homogeneous Hydrogenous Medium

The illustrative model considered in this Section is a Volterra-type integral equation that describes the energy distribution of neutrons in a homogeneous hydrogenous medium (such as a water-moderated/cooled reactor system) containing ²³⁸U (among other materials), which is a heavy element that strongly absorbs neutrons. The distribution of collided neutrons in such a medium is described [19,20,21] by the following linear integral equation of Volterra-type, customarily called the “neutron slowing down equation” for the neutron collision density denoted as $C\left(E\right)$:

$$C\left(E\right)={\displaystyle \frac{{\Sigma }_s\left(E_s\right)S}{{\Sigma }_t\left(E_s\right)E_s}}+{\displaystyle \underset{E}{\overset{E_s}{\int }}C\left(e\right)\frac{{\Sigma }_s\left(e\right)}{{\Sigma }_t\left(e\right)}\frac{de}{e}}.$$

(21)

The various quantities that appear in Eq. (21) are defined as follows:

(i)

The quantity $S$ denotes the rate at which the source neutrons, considered to be monoenergetic, are emitted at the “source energy” $E_s$. Neutron upscattering is considered to be negligeable; therefore, $E_s$ is the highest energy in the medium.

(ii)

The quantity $E$, $0<E_{\mathcal{l}}\le E\le E_s$, denotes the instantaneous energy of the collided neutrons; $E_{\mathcal{l}}$ denotes the lowest neutron energy in the model.

(iii)

The quantity ${\Sigma }_s\left(E\right)$ denotes the medium’s macroscopic scattering cross section, which is defined as follows:

$${\Sigma }_s\left(E\right)\triangleq {\displaystyle \sum _{i=1}^M{w}_i{N}_i{\sigma }_s^{\left(i\right)}\left(E\right)},$$

(22)

where M denotes the number of materials in the medium, $w_i$ denotes the relative weighting of the i^th-material in the medium, $N_i$ denotes the number density of the i^th-material, while ${\sigma }_s^{\left(i\right)}\left(E\right)$ denotes the energy-dependent scattering microscopic cross section of the i^th-material.

(iv)

The quantity ${\Sigma }_t\left(E\right)$ denotes the medium’s macroscopic scattering cross section, which is defined as follows:

$${\Sigma }_t\left(E\right)\triangleq {\displaystyle \sum _{i=1}^M{w}_i{N}_i{\sigma }_t^{\left(i\right)}\left(E\right)},$$

(23)

where ${\sigma }_t^{\left(i\right)}\left(E\right)\ge {\sigma }_s^{\left(i\right)}\left(E\right)$ denotes the energy-dependent total microscopic cross section of the i^th-material. The quantities $w_i$, $N_i$, ${\sigma }_s^{\left(i\right)}\left(E\right)$, ${\sigma }_t^{\left(i\right)}\left(E\right)$ are subject to uncertainties since they are determined from experimentally obtained data.

Notably, the Volterra-type Eq. (21) is a “final-value problem” since the computation is started at the highest-energy value, $C\left(E_s\right)={\Sigma }_s\left(E_s\right)S/{\Sigma }_t\left(E_s\right)E_s$, and progresses towards the lowest energy value $E_{\mathcal{l}}$. Customarily, the solution of Eq. (21) is written in the following form:

$$C\left(E\right)={\displaystyle \frac{{\Sigma }_s\left(E_s\right)S}{{\Sigma }_t\left(E_s\right)E}}exp\left\{-{\displaystyle \underset{E}{\overset{E_s}{\int }}\frac{{\Sigma }_a\left(e\right)}{{\Sigma }_t\left(e\right)}\frac{de}{e}}\right\},$$

(24)

where ${\Sigma }_a\left(E\right)\triangleq {\Sigma }_t\left(E\right)-{\Sigma }_s\left(E\right)$ denotes the medium’s macroscopic absorption cross section. The expression provided in Eq. (24) is amenable to computations of the loss of neutrons due to absorbing materials, particularly in the so-called “resonance” energy region.

A typical “decoder response” for the NIE-Volterra network modeled by Eq. (21) is the energy-averaged collision density, denoted below as $R\left[C\left(E\right)\right]$, which would be measured by a detector having an interaction cross-section ${\Sigma }_d$. Mathematically, this detector-response can be expressed as follows:

$$R\left[C\left(E\right)\right]\triangleq {\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{\Sigma }_d\left(E\right)C\left(E\right)dE};{\Sigma }_d\left(E\right)\triangleq N_d{\sigma }_d\left(E\right),$$

(25)

where $N_d$ and ${\sigma }_d\left(E\right)$ denote, respectively, the detector material’s atomic number density and the microscopic cross section describing the interaction (e.g., absorption) of neutrons with the detector’s material; $N_d$ and ${\sigma }_d\left(E\right)$ can be considered as the “weights” that characterizes the neural net’s “decoder.”

Since the energy-dependence of the cross sections does not play a significant role in the sensitivity analysis of the NIE-Volterra modeled by Eq. (21), the respective microscopic cross-sections will henceforth be considered to be energy-independent for the purpose of illustrating the application of the 1^st-FASAM-NIE-V, in order to simplify the ensuing derivations. For energy-independent cross sections, Eqs. (21) and (25) take on the following forms, respectively:

$$C\left(E\right)=F\left(\theta \right)\left[{\displaystyle \frac{S}{E_s}}+{\displaystyle \underset{E}{\overset{E_s}{\int }}C\left(e\right)\frac{de}{e}}\right],$$

(26)

$$R\left[C\left(E\right)\right]={\Sigma }_d\left(\theta \right){\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}C\left(E\right)dE}.$$

(27)

In Eqs. (26) and (27), the source strength $S$ is an imprecisely-known “weight” that characterizes the neural net’s “encoder.” Furthermore, the (column) vector of parameters denoted as $\theta \triangleq {\left({\theta }_1,\mathrm{...},{\theta }_{TW}\right)}^{†}$ comprises as components the “imprecisely known primary model parameters” (or “weights”, as customarily called when referring to neural nets) and is defined as follows:

$$\theta \triangleq {\left({\theta }_1,\mathrm{...},{\theta }_{TW}\right)}^{†}\triangleq {\left(w_1,\mathrm{...},w_M;N_1,\mathrm{...},N_M;{\sigma }_t^{\left(1\right)},\mathrm{...},{\sigma }_t^{\left(M\right)};{\sigma }_s^{\left(1\right)},\mathrm{...},{\sigma }_s^{\left(M\right)};N_d,{\sigma }_d\right)}^{†};$$

(28)

where $\text{\hspace{0.17em}}TW\triangleq 4\times M+2$ denotes the “total number of imprecisely-known weights/parameters.” These primary model parameters/weights are not known exactly but are affected by uncertainties since they stem from experimental procedures, which determine the nominal/mean/optimal values and the second-order moments of their unknown joint distributions; their third- and higher-order moments are rarely known. It is convenient to denote the nominal values of these primary model parameters/weights by using the superscript “zero” as follows:

$${\theta }^0\triangleq {\left({\theta }_1^0,\mathrm{...},{\theta }_{TW}^0\right)}^{†}\triangleq {\left(w_1^0,\mathrm{...},w_M^0;N_1^0,\mathrm{...},N_M^0;{\sigma }_t^{\left(1,0\right)},\mathrm{...},{\sigma }_t^{\left(M,0\right)};{\sigma }_s^{\left(1,0\right)},\mathrm{...},{\sigma }_s^{\left(M,0\right)};N_d^0,{\sigma }_d^0\right)}^{†};$$

(29)

The “feature function of primary parameters,” $F\left(\theta \right)$, is defined as follows:

$$F\left(\theta \right)\triangleq {\Sigma }_s\left(\theta \right)/{\Sigma }_t\left(\theta \right)\le 1.$$

(30)

The closed-form solution of Eq. (26) has the following expression in terms of the feature function $F\left(\theta \right)$:

$$C\left(E\right)={\displaystyle \frac{S}{E_s}}F\left(\theta \right){\left({\displaystyle \frac{E_s}{E}}\right)}^{F\left(\theta \right)}.$$

(31)

The closed-form expression of the decoder response can be readily obtained by replacing the result obtained in Eq. (31) into Eq. (27) and performing the integration over the energy-variable to obtain:

$$R\left[C\left(E\right)\right]=S{\Sigma }_d\left(\theta \right){\displaystyle \frac{F\left(\theta \right)}{1-F\left(\theta \right)}}\left[1-{\left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)}^{F\left(\theta \right)-1}\right].$$

(32)

The expression obtained in Eq. (32) reveals that the imprecisely known quantities that affect the decoder-response $R\left[C\left(E\right)\right]$ are as follows:

(i)

the source strength $S$;

(ii)

the detector interaction macroscopic cross section ${\Sigma }_d\left(\theta \right)$, which can be considered to be a “feature function” of the model parameters $\theta$;

(iii)

the feature function $F\left(\theta \right)\triangleq {\Sigma }_s\left(\theta \right)/{\Sigma }_t\left(\theta \right)$.

3.1. Application of 1^st-CASAM-NIE-V to Directly Compute the First-Order Sensitivities of the Decoder Response with Respect to the Primary Model Parameters

The first-order sensitivities of the decoder response $R\left[C\left(E\right)\right]$ with respect to the model parameters is obtained by applying the definition of the G-differential to Eq. (26), for arbitrary parameter variations $\delta \theta \triangleq {\left(\delta {\theta }_1,\mathrm{...},\delta {\theta }_{TW}\right)}^{†}\triangleq \theta -{\theta }^0\triangleq {\left({\theta }_1-{\theta }_1^0,\mathrm{...},{\theta }_{TW}-{\theta }_{TW}^0\right)}^{†}$ around the parameters’ nominal values. These parameter variations will induce variations $\delta C\left(E\right)\triangleq C\left(E\right)-C^0\left(E\right)$ in the neutron collision density, around the nominal value $C^0\left(E\right)$ of the neutron collision density. The variations $\delta \theta$ and $\delta C\left(E\right)$ will induce variations $\delta R\left[C^0\left(E\right);{\theta }^0;\delta C\left(E\right);\delta \theta \right]$ in the decoder’s response.

The first-order Gateaux (G-)variation $\delta R\left[C^0\left(E\right);{\theta }^0;\delta C\left(E\right);\delta \theta \right]$ is obtained, by definition, from Eq. (27) as follows:

$$\begin{array}{l}\delta R\left(C^0;{\theta }^0;\delta C;\delta \theta \right)\triangleq {\displaystyle \frac{d}{d\epsilon }}{\left\{\left[{\Sigma }_d\left({\theta }^0+\epsilon \delta \theta \right)\right]{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}\left[C^0\left(E\right)+\epsilon \delta C\left(E\right)\right]dE}\right\}}_{\epsilon =0}\\ ={\left\{\delta R\left(C^0;{\theta }^0;\delta \theta \right)\right\}}_{dir}+{\left\{\delta R\left(C^0;{\theta }^0;\delta C\right)\right\}}_{ind},\end{array}$$

(33)

where the “direct effect” term ${\left\{\delta R\left(C^0;{\theta }^0;\delta \theta \right)\right\}}_{dir}$ arises directly from parameter variations $\delta \theta$ and is defined as follows:

$$\begin{array}{l}{\left\{\delta R\left(C^0;{\theta }^0;\delta \theta \right)\right\}}_{dir}\triangleq \left[{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{C}^0\left(E\right)dE}\right]{{\displaystyle \sum _{i=1}^{TW}\left\{\frac{\partial {\Sigma }_d\left(\theta \right)}{\partial {\theta }_i}\right\}}}_{\theta ={\theta }^0}\delta {\theta }_i\\ ={\left\{\left(\delta N_d\right){\sigma }_d+N_d\left(\delta {\sigma }_d\right)\right\}}_{\theta ={\theta }^0}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{C}^0\left(E\right)dE},\end{array}$$

(34)

while the indirect effect term arises from the variations $\delta C\left(E\right)$ and is defined as follows:

$${\left\{\delta R\left(C^0;{\theta }^0;\delta C\right)\right\}}_{ind}\triangleq {\Sigma }_d\left({\theta }^0\right){\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}\delta C\left(E\right)dE}.$$

(35)

As indicated in Eqs. (34) and (35), both the direct-effect and the indirect-effect term are to be evaluated at the nominal parameter values.

The first-order relation between the variation $\delta C\left(E\right)$ and the parameter variations $\delta {\theta }_i$ is obtained by evaluating the G-variation of Eq. (26) for variations $\delta \theta$ around the nominal parameter values ${\theta }^0$ which yields, by definition, the following NIE-Volterra equation for $\delta C\left(E\right)$:

$$\begin{array}{l}\delta C\left(E\right)\triangleq {\displaystyle \frac{d}{d\epsilon }}{\left\{F\left({\theta }^0+\epsilon \delta \theta \right){\displaystyle \frac{S^0+\epsilon \delta S}{E_s}}+F\left({\theta }^0+\epsilon \delta \theta \right){\displaystyle \underset{E}{\overset{E_s}{\int }}\left[C^0\left(e\right)+\epsilon \delta C\left(e\right)\right]\frac{de}{e}}\right\}}_{\epsilon =0}\\ =F\left({\theta }^0\right){\displaystyle \underset{E}{\overset{E_s}{\int }}\delta C\left(e\right)\frac{de}{e}}+Q\left(E\right),\end{array}$$

(36)

where:

$$\begin{array}{l}Q\left(E\right)\triangleq {\displaystyle \frac{F\left({\theta }^0\right)\delta S}{E_s}}+{\left\{\left[{\displaystyle \frac{S^0}{E_s}}+{\displaystyle \underset{E}{\overset{E_s}{\int }}{C}^0\left(e\right)\frac{de}{e}}\right]\left[{\displaystyle \sum _{i=1}^{TW}\frac{\partial F\left(\theta \right)}{\partial {\theta }_i}}\delta {\theta }_i\right]\right\}}_{\theta ={\theta }^0}\\ ={\displaystyle \frac{F\left({\theta }^0\right)\delta S}{E_s}}+{\left\{{\displaystyle \frac{S^0}{E_s}}{\left({\displaystyle \frac{E_s}{E}}\right)}^{F\left(\theta \right)}\left[{\displaystyle \sum _{i=1}^{TW}\frac{\partial F\left(\theta \right)}{\partial {\theta }_i}}\delta {\theta }_i\right]\right\}}_{\theta ={\theta }^0}.\end{array}$$

(37)

The second equality in Eq. (37) has been obtained by using Eqs. (26) and (31) to eliminate the integral term involving $C^0\left(e\right)$.

The particular form of the first-order derivative $\partial F\left(\theta \right)/\partial {\theta }_i$, which appears in Eq. (37), is obtained by using the definition of $F\left(\theta \right)$ provided in Eq. (30), which yields the following expression:

$${\displaystyle \frac{\partial F\left(\theta \right)}{\partial {\theta }_i}}={\displaystyle \frac{1}{{\Sigma }_t\left(\theta \right)}}{\displaystyle \frac{\partial {\Sigma }_s\left(\theta \right)}{\partial {\theta }_i}}-{\displaystyle \frac{{\Sigma }_s\left(\theta \right)}{{\Sigma }_t^2\left(\theta \right)}}{\displaystyle \frac{\partial {\Sigma }_t\left(\theta \right)}{\partial {\theta }_i}},$$

(38)

In view of the definition provided in Eq. (22), the derivatives $\partial {\Sigma }_s\left(\theta \right)/\partial {\theta }_i$ have the following particular expressions:

$$fori=1,\mathrm{...},M:{\theta }_i\triangleq w_i;\Rightarrow {\displaystyle \frac{\partial {\Sigma }_s\left(\theta \right)}{\partial {\theta }_i}}=N_i{\sigma }_s^{\left(i\right)};$$

(39)

$$fori=M+1,\mathrm{...},2M:\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{\theta }_i\triangleq N_i;\Rightarrow {\displaystyle \frac{\partial {\Sigma }_s\left(\theta \right)}{\partial {\theta }_i}}=w_i{\sigma }_s^{\left(i\right)};$$

(40)

$$fori=3M+1,\mathrm{...},4M:\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{\theta }_i\triangleq {\sigma }_s^{\left(i\right)};\Rightarrow {\displaystyle \frac{\partial {\Sigma }_s\left(\theta \right)}{\partial {\theta }_i}}=w_i{N}_i.$$

(41)

In view of the definition provided in Eq. (23), the derivatives $\partial {\Sigma }_t\left(\theta \right)/\partial {\theta }_i$ have the following particular expressions:

$$fori=1,\mathrm{...},M:{\theta }_i\triangleq w_i;\Rightarrow {\displaystyle \frac{\partial {\Sigma }_t\left(\theta \right)}{\partial {\theta }_i}}=N_i{\sigma }_t^{\left(i\right)};$$

(42)

$$fori=M+1,\mathrm{...},2M:\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{\theta }_i\triangleq N_i;\Rightarrow {\displaystyle \frac{\partial {\Sigma }_t\left(\theta \right)}{\partial {\theta }_i}}=w_i{\sigma }_t^{\left(i\right)};$$

(43)

$$fori=2M+1,\mathrm{...},3M:\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{\theta }_i\triangleq {\sigma }_t^{\left(i\right)};\Rightarrow {\displaystyle \frac{\partial {\Sigma }_t\left(\theta \right)}{\partial {\theta }_i}}=w_i{N}_i.$$

(44)

The NIE-Volterra net represented by Eq. (36) will be called the “1^st-Level Variational Sensitivity System (1^st-LVSS)” and its solution, $\delta C\left(E\right)$, will be called the “1^st-level variational sensitivity function.” It is evident that Eq. (36) would need to be solved $\left(TW+1\right)$-times in order to obtain the variation $\delta C\left(E\right)$ for the source variation $\delta S$ and for every parameter variation $\delta {\theta }_i$ This need for repeatedly solving Eq. (36) can be circumvented by applying the principles of the 1^st-CASAM-NIE-V, generally outlined in Section 2, to eliminate the appearance of the variation $\delta C\left(E\right)$ in the indirect-effect term defined in Eq. (35) while expressing this indirect-effect term as a functional of a first-level adjoint function that does not depend on any parameter variation, as follows:

• Consider that the function $\delta C\left(E\right)\in H_1\left({\Omega }_E\right)$ belongs to a Hilbert space denoted as $H_1\left({\Omega }_E\right)$, which is defined on the domain ${\Omega }_E\triangleq \left[E_{\mathcal{l}},E_s\right]$. The inner product in $H_1\left({\Omega }_E\right)$ of two functions $u\left(E\right)\in H_1\left({\Omega }_E\right)$ and $v\left(E\right)\in H_1\left({\Omega }_E\right)$ will be denoted as ${〈u\left(E\right),v\left(E\right)〉}_1$ and is defined as follows:

  $${〈u\left(E\right),v\left(E\right)〉}_1\triangleq {\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_0}{\int }}u\left(E\right)v\left(E\right)dE}.$$

  (45)
• Form the inner product of Eq. (36) with a vector $a^{\left(1\right)}\left(E\right)\in H_1\left(E\right)$, where the superscript “(1)” indicates “1^st-Level”, to obtain the following relationship:

  $${〈a^{\left(1\right)}\left(E\right),\delta C\left(E\right)-F\left({\theta }^0\right){\displaystyle \underset{E}{\overset{E_s}{\int }}\delta C\left(e\right)\frac{de}{e}}〉}_1={〈a^{\left(1\right)}\left(E\right),Q\left(E\right)〉}_1.$$

  (46)
• Transform the left-side of Eq. (46) as follows:

  $$\begin{array}{l}{〈a^{\left(1\right)}\left(E\right),\delta C\left(E\right)-F\left({\theta }^0\right){\displaystyle \underset{E}{\overset{E_s}{\int }}\delta C\left(e\right)\frac{de}{e}}〉}_1={\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right)\delta C\left(E\right)dE}\\ -F\left({\theta }^0\right){\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right)dE{\displaystyle \underset{E}{\overset{E_s}{\int }}\delta C\left(e\right)\frac{de}{e}}}={\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}\delta C\left(E\right)a^{\left(1\right)}\left(E\right)dE}\\ -F\left({\theta }^0\right){\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}\delta C\left(E\right)\frac{dE}{E}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E}{\int }}{a}^{\left(1\right)}\left(e\right)de}={\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}\delta C\left(E\right)\left[a^{\left(1\right)}\left(E\right)-\frac{F\left({\theta }^0\right)}{E}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E}{\int }}{a}^{\left(1\right)}\left(e\right)de}\right]dE}.\end{array}$$

  (47)
• Require the last term in Eq. (47) to represent the indirect-effect term defined in Eq. (35) which yields the following “1^st-Level Adjoint Sensitivity System (1^st-LASS)” for the first-level adjoint sensitivity function $a^{\left(1\right)}\left(E\right)$:

  $$a^{\left(1\right)}\left(E\right)-{\displaystyle \frac{F\left({\theta }^0\right)}{E}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E}{\int }}{a}^{\left(1\right)}\left(e\right)de}={\Sigma }_d\left({\theta }^0\right).$$

  (50)
  The 1^st-LASS represented by Eq. (50) is a linear NIE-Volterra net, which is independent of any parameter variation and needs to be solved just once to obtain the first-level adjoint sensitivity function $a^{\left(1\right)}\left(E\right)$. Notably, the 1^st-LASS is an “initial-value problem,” in that the computation of $a^{\left(1\right)}\left(E\right)$ commences at the lowest-energy value, where $a^{\left(1\right)}\left(E_{\mathcal{l}}\right)={\Sigma }_d\left({\theta }^0\right)$, and progresses towards the highest-energy value, $E_s$. For further reference, the closed-form solution of Eq. (50) can be obtained by differentiating this equation with respect to $E$ and subsequently integrating the resulting first-order linear differential equation, to obtain the following exact expression:

  $$a^{\left(1\right)}\left(E\right)={\displaystyle \frac{{\Sigma }_d\left(\theta \right)}{1-F\left(\theta \right)}}\left[1-F\left(\theta \right){\left({\displaystyle \frac{E}{E_{\mathcal{l}}}}\right)}^{F\left(\theta \right)-1}\right].$$

  (51)
  The expression on the right-side of Eq. (51) is to be evaluated at the nominal parameter values ${\theta }^0$, but the superscript “zero” has been omitted for notational simplicity.
• Using Eqs. (46), (47) and (50) yields the following expression for the indirect-effect term defined in Eq. (35):

  $$\begin{array}{l}{\left\{\delta R\left(C;\theta ;\delta C\right)\right\}}_{ind}={\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right)Q\left(E\right)dE}={\displaystyle \frac{F\left(\theta \right)\delta S}{E_s}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right)dE}\\ +{\displaystyle \frac{S}{E_s}}\left[{\displaystyle \sum _{i=1}^{TW}\frac{\partial F\left(\theta \right)}{\partial {\theta }_i}}\delta {\theta }_i\right]{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right)dE}{\left({\displaystyle \frac{E_s}{E}}\right)}^{F\left(\theta \right)}.\end{array}$$

  (52)
  The expression on the right-side of Eq. (52) is to be evaluated at the nominal parameter values ${\theta }^0$, but the superscript “zero” has been omitted for notational simplicity.
• Adding the expression obtained in Eq. (52) to the expression of the direct-effect term represented by Eq. (34) yields the following expression for the first-order G-variation $\delta R\left(C^0;{\theta }^0;\delta C;\delta \theta \right)$:

  $$\begin{array}{l}\delta R\left(C^0;{\theta }^0;\delta C;\delta \theta \right)=\left[\left(\delta N_d\right){\sigma }_d+N_d\left(\delta {\sigma }_d\right)\right]{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}C\left(E\right)dE}\\ +{\displaystyle \frac{F\left(\theta \right)\delta S}{E_s}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right)dE}+{\displaystyle \frac{S}{E_s}}\left[{\displaystyle \sum _{i=1}^{TW}\frac{\partial F\left(\theta \right)}{\partial {\theta }_i}}\delta {\theta }_i\right]{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right)dE{\left(\frac{E_s}{E}\right)}^{F\left(\theta \right)}}.\end{array}$$

  (53)
• It follows from Eq. (53) that the first-order sensitivities of the decoder response with respect to the (encoder’s) source strength and the optimal weights/parameters have the following expressions:

  $${\displaystyle \frac{\partial R}{\partial S}}={\displaystyle \frac{F\left(\theta \right)}{E_s}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right)dE};$$

  (54)

  $${\displaystyle \frac{\partial R}{\partial N_d}}={\sigma }_d{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}C\left(E\right)dE};$$

  (55)

  $${\displaystyle \frac{\partial R}{\partial {\sigma }_d}}=N_d{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}C\left(E\right)dE};$$

  (56)

  $${\displaystyle \frac{\partial R}{\partial {\theta }_i}}={\displaystyle \frac{S}{E_s}}{\displaystyle \frac{\partial F\left(\theta \right)}{\partial {\theta }_i}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right){\left(\frac{E_s}{E}\right)}^{F\left(\theta \right)}dE};i=1,\mathrm{...},TW-2.$$

  (57)

Inserting into Eqs. (54)‒(57) the closed-form expression for the neutron collision density obtained in Eq. (31) yields the following closed-form explicit expressions for the first-order sensitivities of the decoder response with respect to the (encoder’s) source strength and the optimal weights/parameters:

$${\displaystyle \frac{\partial R}{\partial S}}={\Sigma }_d\left(\theta \right){\displaystyle \frac{F\left(\theta \right)}{1-F\left(\theta \right)}}\left[1-{\left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)}^{F\left(\theta \right)-1}\right];$$

(58)

$${\displaystyle \frac{\partial R}{\partial N_d}}={\sigma }_{d}S{\displaystyle \frac{F\left(\theta \right)}{1-F\left(\theta \right)}}\left[1-{\left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)}^{F\left(\theta \right)-1}\right];$$

(59)

$${\displaystyle \frac{\partial R}{\partial {\sigma }_d}}=N_{d}S{\displaystyle \frac{F\left(\theta \right)}{1-F\left(\theta \right)}}\left[1-{\left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)}^{F\left(\theta \right)-1}\right]$$

(60)

$$\begin{array}{l}{\displaystyle \frac{\partial R}{\partial {\theta }_i}}={\displaystyle \frac{\partial F\left(\theta \right)}{\partial {\theta }_i}}{\displaystyle \frac{S{\Sigma }_d\left(\theta \right)}{1-F\left(\theta \right)}}\left\{{\displaystyle \frac{1}{1-F\left(\theta \right)}}\left[1-{\left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)}^{F\left(\theta \right)-1}\right]-F\left(\theta \right){\left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)}^{F\left(\theta \right)-1}\ln \left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)\right\};\\ i=1,\mathrm{...},TW.\end{array}$$

(61)

The correctness of the expressions obtained in Eqs. (58)‒(61) can be readily verified by differentiating the expressions of the decoder’s response obtained in Eq. (32).

In practice, only the exact mathematical expression of the 1^st-LASS, namely Eq. (50), and the exact mathematical expression of the first-order sensitivities obtained in Eqs. (54)‒(57) are available. The solution of the 1^st-LASS, which is a linear NIE-Volterra net for the first-level adjoint sensitivity function $a^{\left(1\right)}\left(E\right)$, would need to be obtained numerically, in practice. The numerical solution for $a^{\left(1\right)}\left(E\right)$ would be used to determine the first-order sensitivities stemming from the “indirect-effect” term by using quadrature formulas to evaluate the integrals obtained in Eqs. (54) and (57). It is very important to note that a single “large-scale” computation, for determining numerically the adjoint function $a^{\left(1\right)}\left(E\right)$ by solving the 1^st-LASS (a NIE-Volterra type equation), would be needed for evaluating all of the first-order sensitivities. The numerical computations using quadrature formulas for evaluating the integrals in Eqs. (54) and (57) are considered to be “small-scale” computations.

As has been already observed in the brief remarks following Eq. (37), the computation of the first-order sensitivities of the decoder response with respect to the encoder source strength S and model weights/parameters could also have been computed by numerically solving repeatedly the NIE-Volterra net (1^st-LVSS) represented by Eq. (36). This procedure would be very expensive computationally, since it would require $\left(TW+1\right)$ large scale computations to solve the 1^st-LVSS defined by Eq. (36) in order to obtain the variation $\delta C\left(E\right)$ for every parameter variation $\delta {\theta }_i$ and the source variation $\delta S$. In addition, the same amount of “quadrature” computations would need to be performed using Eq. (35) as would be needed for evaluating the first-order sensitivities using Eqs. (54) and (57).

3.2. Efficient Indirect Computation Using the 1^st-FASAM-NIE-V of the First-Order Sensitivities of the Decoder Response with Respect to Pimary Model Parameters

When feature functions of model parameters such as ${\Sigma }_d\left(\theta \right)$ and $F\left(\theta \right)$ can be identified, as is the case with the NIE-Volterra net and decoder response represented by Eqs. (26) and (27), respectively, it is considerably more efficient to determine the first-order sensitivities of the decoder response with respect to the feature functions and subsequently derive analytically the sensitivities with respect to the primary model parameters by using the “chain rule of differentiation,” as will be shown in this Section. Thus, considering arbitrary variations $\delta {\Sigma }_d\left(\theta \right)\triangleq {\Sigma }_d\left(\theta \right)-{\Sigma }_d\left({\theta }^0\right)$ and $\delta F\left(\theta \right)\triangleq F\left(\theta \right)-F\left({\theta }^0\right)$ around the nominal values ${\Sigma }_d\left({\theta }^0\right)$ and, respectively, $F\left({\theta }^0\right)$, the first-order G-variation of the decoder response has the following expression:

$$\begin{array}{l}\delta R\left(C^0;{\theta }^0;\delta C;\delta {\Sigma }_d\right)\triangleq {\displaystyle \frac{d}{d\epsilon }}{\left\{\left[{\Sigma }_d\left({\theta }^0\right)+\epsilon \delta {\Sigma }_d\left({\theta }^0\right)\right]{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}\left[C^0\left(E\right)+\epsilon \delta C\left(E\right)\right]dE}\right\}}_{\epsilon =0}\\ =\left(\delta {\Sigma }_d\right){\left\{{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}C\left(E\right)dE}\right\}}_{\theta ={\theta }^0}+{\left\{\delta R\left(C^0;{\theta }^0;\delta C\right)\right\}}_{ind},\end{array}$$

(62)

where the expression of the indirect effect term is defined in Eq. (35). The first-order relation between the variation $\delta C\left(E\right)$ and the variations $\delta {\Sigma }_d\left(\theta \right)$ and $\delta F\left(\theta \right)$ is obtained, by definition, from Eq. (26) as follows:

$$\begin{array}{l}\delta C\left(E\right)\triangleq {\displaystyle \frac{d}{d\epsilon }}\left\{\left[F\left({\theta }^0\right)+\epsilon \delta F\left(\theta \right)\right]{\displaystyle \frac{S^0+\epsilon \delta S}{E_s}}+\left[F\left({\theta }^0\right)+\epsilon \delta F\left(\theta \right)\right]\right.\\ {\left.\times {\displaystyle \underset{E}{\overset{E_s}{\int }}\left[C^0\left(e\right)+\epsilon \delta C\left(e\right)\right]\frac{de}{e}}\right\}}_{\epsilon =0}=F\left({\theta }^0\right){\displaystyle \underset{E}{\overset{E_s}{\int }}\delta C\left(e\right)\frac{de}{e}}+Q\left(E\right),\end{array}$$

(63)

where:

$$Q\left(E\right)\triangleq {\displaystyle \frac{F\left({\theta }^0\right)\delta S}{E_s}}+{\left\{{\displaystyle \frac{S^0}{E_s}}{\left({\displaystyle \frac{E_s}{E}}\right)}^{F\left(\theta \right)}\left[\delta F\left(\theta \right)\right]\right\}}_{\theta ={\theta }^0}.$$

(64)

Comparing Eq. (63) to Eq. (36) indicates that the only difference between these equations is the expression of the term $Q\left(E\right)$, which is expressed in terms of $\delta F\left(\theta \right)$ in Eq. (64). Consequently, the first-level adjoint sensitivity function that corresponds to the variational function $\delta C\left(E\right)$ is determined by following the same procedure as outlined in Eqs. (46)‒(50), ultimately obtaining the same 1^st-LASS as was obtained in Eq. (50), having as solution the same expression for $a^{\left(1\right)}\left(E\right)$ as was obtained in Eq. (51). It further follows that the expression of the indirect-effect term will have the following expression:

$${\left\{\delta R\left(C;\theta ;\delta C\right)\right\}}_{ind}=\left(\delta S\right){\displaystyle \frac{F\left(\theta \right)}{E_s}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right)dE}+\left[\delta F\left(\theta \right)\right]{\displaystyle \frac{S}{E_s}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right)dE}{\left({\displaystyle \frac{E_s}{E}}\right)}^{F\left(\theta \right)}.$$

(65)

It follows from Eqs. (62) and (65) that the first-order G-variation $\delta R\left(C^0;{\theta }^0;\delta C;\delta {\Sigma }_d\right)$ has the following expression:

$$\begin{array}{l}\delta R\left(C^0;{\theta }^0;\delta C;\delta {\Sigma }_d\right)=\left(\delta {\Sigma }_d\right){\left\{{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}C\left(E\right)dE}\right\}}_{\theta ={\theta }^0}\\ +\left(\delta S\right){\displaystyle \frac{F\left(\theta \right)}{E_s}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right)dE}+\left[\delta F\left(\theta \right)\right]{\displaystyle \frac{S}{E_s}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right)dE}{\left({\displaystyle \frac{E_s}{E}}\right)}^{F\left(\theta \right)}.\end{array}$$

(66)

As indicated by the expression obtained in Eq. (66), the first-order sensitivities of the decoder response with respect to the feature functions and the encoder’s source strength are as follows:

$${\displaystyle \frac{\partial R}{\partial F\left(\theta \right)}}={\displaystyle \frac{S}{E_s}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right){\left(\frac{E_s}{E}\right)}^{F\left(\theta \right)}dE};$$

(67)

$${\displaystyle \frac{\partial R}{\partial {\Sigma }_d\left(\theta \right)}}={\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}C\left(E\right)dE};$$

(68)

$${\displaystyle \frac{\partial R}{\partial S}}={\displaystyle \frac{F\left(\theta \right)}{E_s}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right)dE}.$$

(69)

The closed-form expressions of the above sensitivities are readily determined by using in Eqs. (67)‒(69) the expressions obtained in Eqs. (51) and (24), and by performing the respective integrations obtain:

$${\displaystyle \frac{\partial R}{\partial F\left(\theta \right)}}={\displaystyle \frac{S{\Sigma }_d\left(\theta \right)}{1-F\left(\theta \right)}}\left\{{\displaystyle \frac{1}{1-F\left(\theta \right)}}\left[1-{\left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)}^{F\left(\theta \right)-1}\right]-F\left(\theta \right){\left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)}^{F\left(\theta \right)-1}\ln \left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)\right\};$$

(70)

$${\displaystyle \frac{\partial R}{\partial {\Sigma }_d\left(\theta \right)}}=S{\displaystyle \frac{F\left(\theta \right)}{1-F\left(\theta \right)}}\left[1-{\left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)}^{F\left(\theta \right)-1}\right];$$

(71)

$${\displaystyle \frac{\partial R}{\partial S}}={\Sigma }_d\left(\theta \right){\displaystyle \frac{F\left(\theta \right)}{1-F\left(\theta \right)}}\left[1-{\left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)}^{F\left(\theta \right)-1}\right].$$

(72)

The first-order sensitivities with respect to the primary parameters are obtained analytically from Eqs. (67) and (68), respectively, by using the following “chain rule” of differentiation:

$${\displaystyle \frac{\partial R}{\partial N_d}}={\displaystyle \frac{\partial R}{\partial {\Sigma }_d\left(\theta \right)}}{\displaystyle \frac{\partial {\Sigma }_d\left(\theta \right)}{\partial N_d}}={\sigma }_d{\displaystyle \frac{\partial R}{\partial {\Sigma }_d\left(\theta \right)}};$$

(73)

$${\displaystyle \frac{\partial R}{\partial {\sigma }_d}}={\displaystyle \frac{\partial R}{\partial {\Sigma }_d\left(\theta \right)}}{\displaystyle \frac{\partial {\Sigma }_d\left(\theta \right)}{\partial {\sigma }_d}}=N_d{\displaystyle \frac{\partial R}{\partial {\Sigma }_d\left(\theta \right)}};$$

(74)

$${\displaystyle \frac{\partial R}{\partial {\theta }_i}}={\displaystyle \frac{\partial R}{\partial F\left(\theta \right)}}{\displaystyle \frac{\partial F\left(\theta \right)}{\partial {\theta }_i}};i=1,\mathrm{...},4M.$$

(75)

The specific expressions of the first-order sensitivities $\partial R/\partial {\theta }_i$, $i=1,\mathrm{...},4M$, are obtained by using Eq. (75) in conjunction with Eq. (69) and Eqs. (38)‒(44).

3.3. Discussion: Direct Versus Indirect Computation of the First-Order Sensitivities of Decoder Response with Respect to the Primary Model Parameters:

The principles of the 1^st-CASAM-NIE-V were applied in Section 3.1 to determine the first-order sensitivities of the decoder response directly with respect to the model’s primary parameters/weights. It has been shown that this procedure requires a single “large-scale” computation to solve a NIE-Volterra equation in order to determine the (single) 1^st-level adjoint sensitivity function $a^{\left(1\right)}\left(E\right)$, which is subsequently used in $4M+1$ integrals that are computed using quadrature formulas. The two additional first-order sensitivities with respect to the components of ${\Sigma }_d\left(\theta \right)$ require a single quadrature involving the forward function $C\left(E\right)$.

The principles of the 1^st-FASAM-NIE-V were applied in Subsection 3.2 to determine the first-order sensitivities of the decoder response with respect to the feature functions. This path required just two (as opposed to $4M+1$) numerical evaluations of (two) integrals using quadrature formulas involving the 1^st-level adjoint sensitivity function $a^{\left(1\right)}\left(E\right)$. The sensitivities of the decoder response with respect to the primary parameters/weights were subsequently determined analytically, using the “chain rule of differentiation” of the explicitly-known expression of the feature function $F\left(\theta \right)$. Evaluating the two additional first-order sensitivities with respect to the components of ${\Sigma }_d\left(\theta \right)$ require a single quadrature involving the forward function $C\left(E\right)$, as in Section 3.1. Evidently, the indirect path presented in Section 3.2 is computationally more efficient, since it requires substantially fewer numerical quadratures than the path presented in Section 3.1. The superiority of the indirect path, via “feature functions,” over the direct computation of sensitivities with respect to the model parameters will be considerably more evident for the computation of second-order sensitivities, as will be shown in the forthcoming Section 4 and Section 5, below.

Of course, when no feature functions can be identified, the 1^st-FASAM-NIE-V methodology becomes identical to the 1^st-CASAM-NIE-V methodology.

4. The Second-Order Features Adjoint Sensitivity Analysis Methodology for Neural Integral Equations of Volterra-Type (2^nd-FASAM-NIE-V)

The second-order sensitivities of the response $R\left[h;F\left(\theta \right)\right]$ defined in Eq. (3) will be computed by conceptually using their basic definitions as being the “first-order sensitivities of the first-order sensitivities.” Thus, the second-order sensitivities stemming from the first-order sensitivities $\partial R\left[h;F\left(\theta \right)\right]/\partial F_j$ are obtained from the first-order G-differential of Eq. (18), for $j=1,\mathrm{...},TF$, as follows:

$$\begin{array}{l}\delta \left({\displaystyle \frac{\partial R}{\partial F_j}}\right)\triangleq {\left\{{\displaystyle \frac{d}{d\epsilon }}\left[{\displaystyle \underset{t_0}{\overset{t_f}{\int }}\frac{\partial D\left[h^0\left(t\right)+\epsilon v^{\left(1\right)}\left(t\right);F\left({\theta }^0\right)+\epsilon \delta F;t\right]}{\partial F_j}dt}\right]\right\}}_{\epsilon =0}\\ +{\left\{{\displaystyle \sum _{i=1}^{TH}\frac{d}{d\epsilon }{\displaystyle \underset{t_0}{\overset{t_f}{\int }}\left[a_i^{\left(1,0\right)}\left(t\right)+\epsilon \delta a_i^{\left(1\right)}\left(t\right)\right]\frac{\partial g_i\left(F^0+\epsilon \delta F;t\right)}{\partial F_j}dt}}\right\}}_{\epsilon =0}\\ +\left\{{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}\left[a_i^{\left(1,0\right)}\left(t\right)+\epsilon \delta a_i^{\left(1\right)}\left(t\right)\right]\frac{\partial {\phi }_i\left(F^0+\epsilon \delta F;t\right)}{\partial F_j}}}\right.\\ {\left.\times {\displaystyle \underset{t_0}{\overset{t}{\int }}{\psi }_i\left[h^0\left(\tau \right)+\epsilon v^{\left(1\right)}\left(\tau \right);F^0+\epsilon \delta F;\tau \right]d\tau }\right\}}_{\epsilon =0}\\ +\left\{{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}\left[a_i^{\left(1,0\right)}\left(t\right)+\epsilon \delta a_i^{\left(1\right)}\left(t\right)\right]{\phi }_i\left(F^0+\epsilon \delta F;t\right)dt}}\right.\\ \times {\left.{\displaystyle \underset{t_0}{\overset{t}{\int }}\frac{\partial {\psi }_i\left[h^0\left(\tau \right)+\epsilon v^{\left(1\right)}\left(\tau \right);F^0+\epsilon \delta F;\tau \right]}{\partial F_j}d\tau }\right\}}_{\epsilon =0}\\ \triangleq \delta {\left(\partial R/\partial F_j\right)}_{dir}+\delta {\left(\partial R/\partial F_j\right)}_{ind};j=1,\mathrm{...},TF.\end{array}$$

(76)

In Eq. (76), the expression of the direct-effect term $\delta {\left(\partial R/\partial F_j\right)}_{dir}$ is obtained after performing the operations with respect to the scalar $\epsilon$ and comprises the variations $\delta F$ (stemming from variations in the model parameters), being defined as follows:

$$\begin{array}{l}\delta {\left({\displaystyle \frac{\partial R}{\partial F_j}}\right)}_{dir}\triangleq {\displaystyle \sum _{k=1}^{TF}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}\left[\frac{{\partial }^{2}D\left[h\left(t\right);F\left(\theta \right);t\right]}{\partial F_k\partial F_j}\delta F_k\right]dt}}+{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}\left[a_i^{\left(1\right)}\left(t\right){\displaystyle \sum _{k=1}^{TF}\frac{{\partial }^2{g}_i\left(h;F\right)}{\partial F_k\partial F_j}\delta F_k}\right]dt}}\\ +{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(t\right)\left[{\displaystyle \sum _{k=1}^{TF}\frac{{\partial }^2{\phi }_i\left(F;t\right)}{\partial F_k\partial F_j}\delta F_k}\right]{\displaystyle \underset{t_0}{\overset{t}{\int }}{\psi }_i\left[h\left(\tau \right);F;\tau \right]d\tau }}}\\ +{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(t\right)\frac{\partial {\phi }_i\left(F;t\right)}{\partial F_j}{\displaystyle \underset{t_0}{\overset{t}{\int }}\left[{\displaystyle \sum _{k=1}^{TF}\frac{\partial {\psi }_i\left[h\left(\tau \right);F;\tau \right]}{\partial F_k}\delta F_k}\right]d\tau }}}\\ +{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(t\right)\left[{\displaystyle \sum _{k=1}^{TF}\frac{\partial {\phi }_i\left(F;t\right)}{\partial F_k}\delta F_k}\right]dt}}{\displaystyle \underset{t_0}{\overset{t}{\int }}\frac{\partial {\psi }_i\left[h\left(\tau \right);F;\tau \right]}{\partial F_j}d\tau }\\ +{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(t\right){\phi }_i\left(F;t\right)dt}}\left[{\displaystyle \sum _{k=1}^{TF}\frac{{\partial }^2{\psi }_i\left[h\left(\tau \right);F;\tau \right]}{\partial F_k\partial F_j}\delta F_k}\right];j=1,\mathrm{...},TF.\end{array}$$

(77)

The expression on the right-side of Eq. (77) is to be evaluated at the nominal/optimal values for the respective model parameters, but this fact has not been indicated explicitly in order to simplify the notation.

The expression of the indirect-effect term $\delta {\left(\partial R/\partial F_j\right)}_{ind}$ defined in Eq. (76) is obtained after performing the operations with respect to the scalar $\epsilon$ and comprises the variations $v^{\left(1\right)}\left(t\right)$ and $\delta a^{\left(1\right)}\left(t\right)\triangleq {\left[\delta a_1^{\left(1\right)}\left(t\right),\dots ,\delta a_{TH}^{\left(1\right)}\left(t\right)\right]}^{†}$, as follows:

$$\begin{array}{l}\delta {\left({\displaystyle \frac{\partial R}{\partial F_j}}\right)}_{ind}\triangleq {\displaystyle \sum _{k=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}\frac{{\partial }^{2}D\left[h\left(t\right);F\left(\theta \right);t\right]}{\partial h_k\left(t\right)\partial F_j}{v}_k\left(t\right)dt}}+{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}\delta a_i^{\left(1\right)}\left(t\right)\frac{\partial g_i\left(F;t\right)}{\partial F_j}dt}}\\ +{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}\delta a_i^{\left(1\right)}\left(t\right)\frac{\partial {\phi }_i\left(F;t\right)}{\partial F_j}{\displaystyle \underset{t_0}{\overset{t}{\int }}{\psi }_i\left[h\left(\tau \right);F;\tau \right]d\tau }}}\\ +{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(t\right)\frac{\partial {\phi }_i\left(F;t\right)}{\partial F_j}}}{\displaystyle \underset{t_0}{\overset{t}{\int }}{\displaystyle \sum _{k=1}^{TH}\frac{\partial {\psi }_i\left[h\left(\tau \right);F;\tau \right]}{\partial h_k\left(\tau \right)}{v}_k\left(\tau \right)d\tau }}\\ +{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}\delta a_i^{\left(1\right)}\left(t\right){\phi }_i\left(F;t\right)dt}}{\displaystyle \underset{t_0}{\overset{t}{\int }}\frac{\partial {\psi }_i\left[h\left(\tau \right);F;\tau \right]}{\partial F_j}d\tau }\\ +{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(t\right){\phi }_i\left(F;t\right)dt{\displaystyle \underset{t_0}{\overset{t}{\int }}{\displaystyle \sum _{k=1}^{TH}\left[\frac{{\partial }^2{\psi }_i\left[h\left(\tau \right);F;\tau \right]}{\partial h_k\left(\tau \right)\partial F_j}{v}_k\left(\tau \right)\right]}d\tau }}};j=1,\mathrm{...},TF.\end{array}$$

(78)

The expressions in Eq. (78) are to be evaluated at the nominal values of the respective functions and parameters, but the respective indication (i.e., the superscript “zero”) has been omitted in order to simplify the notation.

The direct-effect term $\delta {\left(\partial R/\partial F_j\right)}_{dir}$ can be evaluated at this time for all variations $\delta F$, but the indirect-effect term $\delta {\left(\partial R/\partial F_j\right)}_{ind}$ can be evaluated only after having determined the variations $v^{\left(1\right)}\left(t\right)$ and $\delta a^{(1)}\left(t\right)$. The variation $v^{\left(1\right)}\left(t\right)$ is the solution of the 1^st-LVSS defined by Eq. (9). On the other hand, the variational function $\delta a^{(1)}\left(t\right)$ is the solution of the system of equations obtained by G-differentiating the 1^st-LASS. By definition, the G-differential of Eq. (15) is obtained as follows, for $i=1,\mathrm{...},TH$:

$$\begin{array}{l}{\left\{{\displaystyle \frac{d}{d\epsilon }}\left[a_i^{\left(1,0\right)}\left(t\right)+\epsilon \delta a_i^{\left(1\right)}\left(t\right)\right]\right\}}_{\epsilon =0}-\left\{{\displaystyle \frac{d}{d\epsilon }}\left[{\displaystyle \sum _{k=1}^{TH}\frac{\partial {\psi }_k\left[h^0\left(t\right)+\epsilon v^{\left(1\right)}\left(t\right);F^0+\epsilon \delta F;t\right]}{\partial h_i\left(t\right)}}\right]\right.\\ \times {\left.{\displaystyle \underset{t}{\overset{t_f}{\int }}\left[a_k^{\left(1,0\right)}\left(\tau \right)+\epsilon \delta a_k^{\left(1\right)}\left(\tau \right)\right]{\phi }_k\left(F^0+\epsilon \delta F;\tau \right)d\tau }\right\}}_{\epsilon =0}\\ ={\left\{{\displaystyle \frac{d}{d\epsilon }}{\displaystyle \frac{\partial D\left[h^0\left(t\right)+\epsilon v^{\left(1\right)}\left(t\right);F^0+\epsilon \delta F;t\right]}{\partial h_i\left(t\right)}}\right\}}_{\epsilon =0}.\end{array}$$

(79)

Performing the operations indicated in Eq. (79) and rearranging the various terms yields the following relations, for $i=1,\mathrm{...},TH$:

$$\begin{array}{l}\delta a_i^{\left(1\right)}\left(t\right)-{\displaystyle \sum _{k=1}^{TH}\frac{\partial {\psi }_k\left[h\left(t\right);F;t\right]}{\partial h_i\left(t\right)}{\displaystyle \underset{t}{\overset{t_f}{\int }}\delta a_k^{\left(1\right)}\left(\tau \right){\phi }_k\left(F;\tau \right)d\tau }}-{\displaystyle \sum _{m=1}^{TH}\frac{{\partial }^{2}D\left[h\left(t\right);F;t\right]}{\partial h_m\left(t\right)\partial h_i\left(t\right)}{v}_m^{\left(1\right)}\left(t\right)}\\ -{\displaystyle \sum _{k=1}^{TH}\left[{\displaystyle \sum _{m=1}^{TH}\frac{{\partial }^2{\psi }_k\left[h\left(t\right);F;t\right]}{\partial h_m\left(t\right)\partial h_i\left(t\right)}{v}_m^{\left(1\right)}\left(t\right)}\right]{\displaystyle \underset{t}{\overset{t_f}{\int }}{a}_k^{\left(1\right)}\left(\tau \right){\phi }_k\left(F;\tau \right)d\tau }}={\displaystyle \sum _{n=1}^{TF}{S}_{in}\left(F;t\right)\delta F_n,}\end{array}$$

(80)

where:

$$\begin{array}{l}{\displaystyle \sum _{n=1}^{TF}{S}_{in}\left(F;t\right)\delta F_n}\triangleq {\displaystyle \sum _{k=1}^{TH}\left[{\displaystyle \sum _{n=1}^{TF}\frac{{\partial }^2{\psi }_k\left[h\left(t\right);F;t\right]}{\partial F_n\partial h_i\left(t\right)}\delta F_n}\right]{\displaystyle \underset{t}{\overset{t_f}{\int }}{a}_k^{\left(1\right)}\left(\tau \right){\phi }_k\left(F;\tau \right)d\tau }}\\ +{\displaystyle \sum _{k=1}^{TH}\frac{\partial {\psi }_k\left[h\left(t\right);F;t\right]}{\partial h_i\left(t\right)}{\displaystyle \sum _{n=1}^{TF}{\displaystyle \underset{t}{\overset{t_f}{\int }}\left[a_k^{\left(1\right)}\left(\tau \right)\frac{\partial {\phi }_k\left(F;\tau \right)}{\partial F_n}\delta F_n\right]d\tau }}}+{\displaystyle \sum _{n=1}^{TF}\frac{{\partial }^{2}D\left[h\left(t\right);F;t\right]}{\partial F_n\partial h_i\left(t\right)}}\delta F_n.\end{array}$$

(81)

As indicated by the result obtained in Eq. (80), the variations $\delta a^{(1)}\left(t\right)$are coupled to the variations $v^{\left(1\right)}\left(t\right)$. Therefore, they can be obtained by simultaneously solving Eqs. (80) and (9), which together will be called the “2^nd-Level Variational Sensitivity System (2^nd-LVSS).” The solution of the 2^nd-LVSS, namely the vector $v^{\left(2\right)}\left(t\right)\triangleq {\left[\delta a^{(1)}\left(t\right),v^{\left(1\right)}\left(t\right)\right]}^{†}$, will be called the “2^nd-level variational sensitivity function.” Since the 2^nd-LVSS depends on the variations $\delta F$ (stemming from variations in the model parameters), it would need to be solved anew for each such variation. The repeated solving of the 2^nd-LVSS can be avoided by following the general principles underlying the 2^nd-FASAM [16], which considers the function $v^{\left(2\right)}\left(t\right)\triangleq {\left[v^{\left(1\right)}\left(t\right),\delta a^{(1)}\left(t\right)\right]}^{†}$to be an element in a Hilbert space denoted as $H_2\left({\Omega }_t\right)$. The Hilbert space $H_2\left({\Omega }_t\right)$ is considered to be endowed with an inner product denoted as ${〈{\chi }^{\left(2\right)},{\eta }^{\left(2\right)}〉}_2$, between two vectors ${\chi }^{\left(2\right)}\left(t\right)\triangleq {\left[{\chi }_1^{\left(2\right)}\left(t\right),{\chi }_2^{\left(2\right)}\left(t\right)\right]}^{†}\in H_2\left({\Omega }_t\right)$ and ${\eta }^{\left(2\right)}\left(t\right)={\left[{\eta }_1^{\left(2\right)}\left(t\right),{\eta }_2^{\left(2\right)}\left(t\right)\right]}^{†}\in H_2\left({\Omega }_t\right)$, with ${\eta }_1^{\left(2\right)}\left(t\right)\triangleq {\left[{\eta }_{1,1}^{\left(2\right)}\left(t\right),\dots ,{\eta }_{1,TH}^{\left(2\right)}\left(t\right)\right]}^{†}$, ${\eta }_2^{\left(2\right)}\left(t\right)\triangleq {\left[{\eta }_{2,1}^{\left(2\right)}\left(t\right),\dots ,{\eta }_{2,TH}^{\left(2\right)}\left(t\right)\right]}^{†}$, ${\chi }_1^{\left(2\right)}\left(t\right)\triangleq {\left[{\chi }_{1,1}^{\left(2\right)}\left(t\right),\dots ,{\chi }_{1,TH}^{\left(2\right)}\left(t\right)\right]}^{†}$, ${\chi }_2^{\left(2\right)}\left(t\right)\triangleq {\left[{\chi }_{2,1}^{\left(2\right)}\left(t\right),\dots ,{\chi }_{2,TH}^{\left(2\right)}\left(t\right)\right]}^{†}$, which is defined as follows:

$$\begin{array}{l}{〈{\chi }^{\left(2\right)},{\eta }^{\left(2\right)}〉}_2\triangleq {\displaystyle \underset{t_0}{\overset{t_f}{\int }}{\left[{\chi }^{\left(2\right)}\left(t\right)\right]}^{†}{\eta }^{\left(2\right)}\left(t\right)dt}={〈{\chi }_1^{\left(2\right)},{\eta }_1^{\left(2\right)}〉}_1+{〈{\chi }_2^{\left(2\right)},{\eta }_2^{\left(2\right)}〉}_1\\ ={\displaystyle \sum _{j=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{\chi }_{1,j}^{\left(2\right)}\left(t\right){\eta }_{1,j}^{\left(2\right)}\left(t\right)dt}}+{\displaystyle \sum _{j=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{\chi }_{2,j}^{\left(2\right)}\left(t\right){\eta }_{2,j}^{\left(2\right)}\left(t\right)dt.}}\end{array}$$

(82)

Following the general principles underlying the 2^nd-FASAM [16], the function $v^{\left(2\right)}\left(t\right)\triangleq {\left[v^{\left(1\right)}\left(t\right),\delta a^{(1)}\left(t\right)\right]}^{†}$ will be eliminated from the expression of each indirect-effect terms $\delta {\left(\partial R/\partial F_j\right)}_{ind}$, $j=1,\mathrm{...},TF$, defined in Eq. (78). This elimination is achieved by considering, for each index $j=1,\mathrm{...},TF$, a vector-valued function denoted as $a^{\left(2\right)}\left(t;j\right)={\left[a_1^{\left(2\right)}\left(t;j\right),a_2^{\left(2\right)}\left(t;j\right)\right]}^{†}\in H_2\left({\Omega }_t\right)$, with $a_1^{\left(2\right)}\left(t;j\right)\triangleq {\left[a_{1,1}^{\left(2\right)}\left(t;j\right),\dots ,a_{1,TH}^{\left(2\right)}\left(t;j\right)\right]}^{†}$and $a_2^{\left(2\right)}\left(t;j\right)\triangleq {\left[a_{2,1}^{\left(2\right)}\left(t;j\right),\dots ,a_{2,TH}^{\left(2\right)}\left(t;j\right)\right]}^{†}$. Using the definition provided in Eq. (82), we construct the inner product of Eqs. (9) and (80) with the vector $a^{\left(2\right)}\left(t;j\right)={\left[a_1^{\left(2\right)}\left(t;j\right),a_2^{\left(2\right)}\left(t;j\right)\right]}^{†}$, to obtain the following relation:

$$\begin{array}{l}{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_{1,i}^{\left(2\right)}\left(t;j\right)v_i^{\left(1\right)}\left(t\right)dt}}-{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_{1,i}^{\left(2\right)}\left(t;j\right){\phi }_i\left(F;t\right)dt{\displaystyle \sum _{k=1}^{TH}{\displaystyle \underset{t_0}{\overset{t}{\int }}\frac{\partial {\psi }_i\left[h\left(\tau \right);F;\tau \right]}{\partial h_k\left(\tau \right)}{v}_k^{\left(1\right)}\left(\tau \right)d\tau }}}}\\ +{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_{2,i}^{\left(2\right)}\left(t;j\right)\delta a_i^{\left(1\right)}\left(t\right)dt}}-{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_{2,i}^{\left(2\right)}\left(t;j\right){\displaystyle \sum _{k=1}^{TH}\frac{\partial {\psi }_k\left[h\left(t\right);F;t\right]}{\partial h_i\left(t\right)}}dt}}{\displaystyle \underset{t}{\overset{t_f}{\int }}\delta a_k^{\left(1\right)}\left(\tau \right){\phi }_k\left(F;\tau \right)d\tau }\\ -{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_{2,i}^{\left(2\right)}\left(t;j\right){\displaystyle \sum _{m=1}^{TH}\frac{{\partial }^{2}D\left[h\left(t\right);F;t\right]}{\partial h_m\left(t\right)\partial h_i\left(t\right)}{v}_m^{\left(1\right)}\left(t\right)dt}}}\\ -{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_{2,i}^{\left(2\right)}\left(t;j\right)}}{\displaystyle \sum _{k=1}^{TH}\left[{\displaystyle \sum _{m=1}^{TH}\frac{{\partial }^2{\psi }_k\left[h\left(t\right);F;t\right]}{\partial h_m\left(t\right)\partial h_i\left(t\right)}{v}_m^{\left(1\right)}\left(t\right)}\right]dt{\displaystyle \underset{t}{\overset{t_f}{\int }}{a}_k^{\left(1\right)}\left(\tau \right){\phi }_k\left(F;\tau \right)d\tau }}=Q^{\left(2\right)},\\ \end{array}$$

(83)

where:

$$Q^{\left(2\right)}\triangleq {\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_{1,i}^{\left(2\right)}\left(t\right)dt}}\left[{\displaystyle \sum _{k=1}^{TF}{q}_{ik}\left(F;t\right)\delta F_k}\right]+{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_{2,i}^{\left(2\right)}\left(t\right)dt}}\left[{\displaystyle \sum _{k=1}^{TF}{S}_{ik}\left(F;t\right)\delta F_k}\right].$$

(84)

Following the principles of the 2^nd-CASAM [16], the left-side of Eq. (83) will be identified with the indirect-effect term defined in Eq. (78), thereby determining the (yet undetermined) functions $a^{\left(2\right)}\left(t;j\right)={\left[a_1^{\left(2\right)}\left(t;j\right),a_2^{\left(2\right)}\left(t;j\right)\right]}^{†}$. For this purpose, the right-side of Eq. (78) is cast in the form of the inner product ${〈v^{\left(2\right)}\left(t\right),\left[\right]〉}_2={〈v^{\left(1\right)}\left(t\right),\left[\right]〉}_1+{〈\delta a^{(1)}\left(t\right),\left[\right]〉}_1$. The terms on the right-side of Eq. (78) involving the components of the function $\delta a^{(1)}\left(t\right)$ are already in the desired format, but the terms involving the components of the function $v^{\left(1\right)}\left(t\right)$ must be re-arranged, as follows:

(i)

The fourth term on the right-side of Eq. (78), is recast by using “integration by parts” as follows:

$$\begin{array}{l}{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(t\right)\frac{\partial {\phi }_i\left(F;t\right)}{\partial F_j}}}{\displaystyle \underset{t_0}{\overset{t}{\int }}{\displaystyle \sum _{k=1}^{TH}\frac{\partial {\psi }_i\left[h\left(\tau \right);F;\tau \right]}{\partial h_k\left(\tau \right)}{v}_k\left(\tau \right)d\tau }}\\ ={{\displaystyle \sum _{i=1}^{TH}{\displaystyle \sum _{k=1}^{TH}\left\{\left[{\displaystyle \underset{t_0}{\overset{t}{\int }}{a}_i^{\left(1\right)}\left(t\right)\frac{\partial {\phi }_i\left(F;t\right)}{\partial F_j}dt}\right]\left[{\displaystyle \underset{t_0}{\overset{t}{\int }}\frac{\partial {\psi }_i\left[h\left(\tau \right);F;\tau \right]}{\partial h_k\left(\tau \right)}{v}_k\left(\tau \right)d\tau }\right]\right\}}}}_{t=t_f}\\ -{{\displaystyle \sum _{i=1}^{TH}{\displaystyle \sum _{k=1}^{TH}\left\{\left[{\displaystyle \underset{t_0}{\overset{t}{\int }}{a}_i^{\left(1\right)}\left(t\right)\frac{\partial {\phi }_i\left(F;t\right)}{\partial F_j}dt}\right]\left[{\displaystyle \underset{t_0}{\overset{t}{\int }}\frac{\partial {\psi }_i\left[h\left(\tau \right);F;\tau \right]}{\partial h_k\left(\tau \right)}{v}_k\left(\tau \right)d\tau }\right]\right\}}}}_{t=t_0}\\ -{\displaystyle \sum _{i=1}^{TH}{\displaystyle \sum _{k=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}\frac{\partial {\psi }_i\left[h\left(t\right);F;t\right]}{\partial h_k\left(t\right)}{v}_k\left(t\right)}{\displaystyle \underset{t_0}{\overset{t}{\int }}{a}_i^{\left(1\right)}\left(\tau \right)\frac{\partial {\phi }_i\left(F;\tau \right)}{\partial F_j}d\tau }}}\\ ={\displaystyle \sum _{i=1}^{TH}{\displaystyle \sum _{k=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}\frac{\partial {\psi }_i\left[h\left(t\right);F;t\right]}{\partial h_k\left(t\right)}{v}_k\left(t\right)dt}{\displaystyle \underset{t}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(\tau \right)\frac{\partial {\phi }_i\left(F;\tau \right)}{\partial F_j}d\tau }}}.\end{array}$$

(85)

(ii)

The sixth (last) term on the right-side of Eq. (78), is recast by using “integration by parts,” as above, to obtain the following relation:

$$\begin{array}{l}{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(t\right){\phi }_i\left(F;t\right)dt{\displaystyle \underset{t_0}{\overset{t}{\int }}{\displaystyle \sum _{k=1}^{TH}\left[\frac{{\partial }^2{\psi }_i\left[h\left(\tau \right);F;\tau \right]}{\partial h_k\left(\tau \right)\partial F_j}{v}_k\left(\tau \right)\right]}d\tau }}}\\ ={\displaystyle \sum _{i=1}^{TH}{\displaystyle \sum _{k=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}\frac{{\partial }^2{\psi }_i\left[h\left(t\right);F;t\right]}{\partial h_k\left(t\right)\partial F_j}{v}_k\left(t\right)dt}{\displaystyle \underset{t}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(\tau \right){\phi }_i\left(F;\tau \right)d\tau }}}.\end{array}$$

(86)

Using in Eq. (78) the results obtained in Eqs. (85) and (86) yields the following expression for the indirect-effect term, for $j=1,\mathrm{...},TF$:

$$\begin{array}{l}\delta {\left({\displaystyle \frac{\partial R}{\partial F_j}}\right)}_{ind}\triangleq {\displaystyle \sum _{k=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{v}_k\left(t\right)dt}\left\{\frac{{\partial }^{2}D\left[h\left(t\right);F\left(\theta \right);t\right]}{\partial h_k\left(t\right)\partial F_j}\right.}\\ +{\displaystyle \sum _{i=1}^{TH}\frac{\partial {\psi }_i\left[h\left(t\right);F;t\right]}{\partial h_k\left(t\right)}{\displaystyle \underset{t}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(\tau \right)\frac{\partial {\phi }_i\left(F;\tau \right)}{\partial F_j}d\tau }}\left.+{\displaystyle \sum _{i=1}^{TH}\frac{{\partial }^2{\psi }_i\left[h\left(t\right);F;t\right]}{\partial h_k\left(t\right)\partial F_j}{\displaystyle \underset{t}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(\tau \right){\phi }_i\left(F;\tau \right)d\tau }}\right\}\\ +{\displaystyle \sum _{k=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}\delta a_k^{\left(1\right)}\left(t\right)\left\{\frac{\partial g_k\left(F;t\right)}{\partial F_j}\right.}}+{\displaystyle \frac{\partial {\phi }_k\left(F;t\right)}{\partial F_j}}{\displaystyle \underset{t_0}{\overset{t}{\int }}{\psi }_k\left[h\left(\tau \right);F;\tau \right]d\tau }\left.+{\phi }_k\left(F;t\right){\displaystyle \underset{t_0}{\overset{t}{\int }}\frac{\partial {\psi }_k}{\partial F_j}d\tau }\right\}dt.\end{array}$$

(87)

The left-side of Eq. (83) is now recast in the form of the inner product ${〈v^{\left(2\right)}\left(t\right),\left[\right]〉}_2$ by performing the following operations:

(i)

The second term on the left-side of Eq. (83) is rearranged by using “integration by parts” as follows:

$$\begin{array}{l}{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_{1,i}^{\left(2\right)}\left(t;j\right){\phi }_i\left(F;t\right)dt{\displaystyle \sum _{k=1}^{TH}{\displaystyle \underset{t_0}{\overset{t}{\int }}\frac{\partial {\psi }_i\left[h\left(\tau \right);F;\tau \right]}{\partial h_k\left(\tau \right)}{v}_k^{\left(1\right)}\left(\tau \right)d\tau }}}}\\ ={\displaystyle \sum _{k=1}^{TH}{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{v}_k\left(t\right)\frac{\partial {\psi }_i\left[h\left(t\right);F;t\right]}{\partial h_k\left(t\right)}dt}{\displaystyle \underset{t}{\overset{t_f}{\int }}{a}_{1,i}^{\left(2\right)}\left(\tau ;j\right){\phi }_i\left(F;\tau \right)d\tau }}}.\end{array}$$

(88)

(ii)

The fourth term on the left-side of Eq. (83) is rearranged by using “integration by parts” as follows:

$$\begin{array}{l}{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_{2,i}^{\left(2\right)}\left(t;j\right){\displaystyle \sum _{k=1}^{TH}\frac{\partial {\psi }_k\left[h\left(t\right);F;t\right]}{\partial h_i\left(t\right)}}dt}}{\displaystyle \underset{t}{\overset{t_f}{\int }}\delta a_k^{\left(1\right)}\left(\tau \right){\phi }_k\left(F;\tau \right)d\tau }\\ ={{\displaystyle \sum _{i=1}^{TH}{\displaystyle \sum _{k=1}^{TH}\left\{\left[{\displaystyle \underset{t_0}{\overset{t}{\int }}{a}_{2,i}^{\left(2\right)}\left(t;j\right){\displaystyle \sum _{k=1}^{TH}\frac{\partial {\psi }_k\left[h\left(t\right);F;t\right]}{\partial h_i\left(t\right)}}dt}\right]\left[{\displaystyle \underset{t}{\overset{t_f}{\int }}\delta a_k^{\left(1\right)}\left(\tau \right){\phi }_k\left(F;\tau \right)d\tau }\right]\right\}}}}_{t=t_f}\\ -{{\displaystyle \sum _{i=1}^{TH}{\displaystyle \sum _{k=1}^{TH}\left\{\left[{\displaystyle \underset{t_0}{\overset{t}{\int }}{a}_{2,i}^{\left(2\right)}\left(t;j\right){\displaystyle \sum _{k=1}^{TH}\frac{\partial {\psi }_k\left[h\left(t\right);F;t\right]}{\partial h_i\left(t\right)}}dt}\right]\left[{\displaystyle \underset{t}{\overset{t_f}{\int }}\delta a_k^{\left(1\right)}\left(\tau \right){\phi }_k\left(F;\tau \right)d\tau }\right]\right\}}}}_{t=t_0}\\ -{\displaystyle \sum _{i=1}^{TH}{\displaystyle \sum _{k=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}dt\left[-\delta a_k^{\left(1\right)}\left(t\right){\phi }_k\left(F;t\right)\right]\left[{\displaystyle \underset{t_0}{\overset{t}{\int }}{a}_{2,i}^{\left(2\right)}\left(\tau ;j\right)\frac{\partial {\psi }_k\left[h\left(\tau \right);F;\tau \right]}{\partial h_i\left(\tau \right)}d\tau }\right]}}}\\ ={\displaystyle \sum _{k=1}^{TH}{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}\delta a_k^{\left(1\right)}\left(t\right){\phi }_k\left(F;t\right)dt{\displaystyle \underset{t_0}{\overset{t}{\int }}{a}_{2,i}^{\left(2\right)}\left(\tau ;j\right)\frac{\partial {\psi }_k\left[h\left(\tau \right);F;\tau \right]}{\partial h_i\left(\tau \right)}d\tau }}}}\\ \end{array}$$

(89)

(iii)

The fifth term on the left-side of Eq. (83) is rearranged as follows:

$$\begin{array}{l}{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_{2,i}^{\left(2\right)}\left(t;j\right){\displaystyle \sum _{k=1}^{TH}\frac{{\partial }^{2}D\left[h\left(t\right);F;t\right]}{\partial h_k\left(t\right)\partial h_i\left(t\right)}{v}_k^{\left(1\right)}\left(t\right)dt}}}\\ ={\displaystyle \sum _{k=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{v}_k^{\left(1\right)}\left(t\right)}\left[{\displaystyle \sum _{i=1}^{TH}\frac{{\partial }^{2}D\left[h\left(t\right);F;t\right]}{\partial h_k\left(t\right)\partial h_i\left(t\right)}{a}_{2,i}^{\left(2\right)}\left(t;j\right)}\right]dt}.\end{array}$$

(90)

(iv)

The sixth term on the left-side of Eq. (83) is rearranged as follows:

$$\begin{array}{l}{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_{2,i}^{\left(2\right)}\left(t\right)}}{\displaystyle \sum _{k=1}^{TH}\left[{\displaystyle \sum _{m=1}^{TH}\frac{{\partial }^2{\psi }_k\left[h\left(t\right);F;t\right]}{\partial h_m\left(t\right)\partial h_i\left(t\right)}{v}_m^{\left(1\right)}\left(t\right)}\right]dt{\displaystyle \underset{t}{\overset{t_f}{\int }}{a}_k^{\left(1\right)}\left(\tau \right){\phi }_k\left(F;\tau \right)d\tau }}\\ ={\displaystyle \sum _{m=1}^{TH}{\displaystyle \sum _{i=1}^{TH}{\displaystyle \sum _{k=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{v}_m^{\left(1\right)}\left(t\right)}\left[\frac{{\partial }^2{\psi }_k\left[h\left(t\right);F;t\right]}{\partial h_m\left(t\right)\partial h_i\left(t\right)}{a}_{2,i}^{\left(2\right)}\left(t\right){\displaystyle \underset{t}{\overset{t_f}{\int }}{a}_k^{\left(1\right)}\left(\tau \right){\phi }_k\left(F;\tau \right)d\tau }\right]dt}}}\\ ={\displaystyle \sum _{k=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{v}_k^{\left(1\right)}\left(t\right)}{\displaystyle \sum _{i=1}^{TH}{\displaystyle \sum _{n=1}^{TH}\left[\frac{{\partial }^2{\psi }_n\left[h\left(t\right);F;t\right]}{\partial h_k\left(t\right)\partial h_i\left(t\right)}{a}_{2,i}^{\left(2\right)}\left(t\right){\displaystyle \underset{t}{\overset{t_f}{\int }}{a}_n^{\left(1\right)}\left(\tau \right){\phi }_n\left(F;\tau \right)d\tau }\right]dt}}}\end{array}$$

(91)

Inserting the results obtained in Eqs. (88)‒(91) into the left-side of Eq. (83) yields the following relation:

$$\begin{array}{l}{Q}^{\left(2\right)}={\displaystyle \sum _{k=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}dt{v}_k^{\left(1\right)}\left(t\right)}}\left\{a_{1,k}^{\left(2\right)}\left(t;j\right)\right.-{\displaystyle \sum _{i=1}^{TH}\frac{\partial {\psi }_i\left[h\left(t\right);F;t\right]}{\partial h_k\left(t\right)}{\displaystyle \underset{t}{\overset{t_f}{\int }}{a}_{1,i}^{\left(2\right)}\left(\tau ;j\right){\phi }_i\left(F;\tau \right)d\tau }}\\ \left.-{\displaystyle \sum _{i=1}^{TH}\frac{{\partial }^{2}D\left[h\left(t\right);F;t\right]}{\partial h_k\left(t\right)\partial h_i\left(t\right)}{a}_{2,i}^{\left(2\right)}\left(t;j\right)}-{\displaystyle \sum _{i=1}^{TH}{\displaystyle \sum _{n=1}^{TH}\frac{{\partial }^2{\psi }_n\left[h\left(t\right);F;t\right]}{\partial h_k\left(t\right)\partial h_i\left(t\right)}{a}_{2,i}^{\left(2\right)}\left(t\right){\displaystyle \underset{t}{\overset{t_f}{\int }}{a}_n^{\left(1\right)}\left(\tau \right){\phi }_n\left(F;\tau \right)d\tau }}}\right\}\\ +{\displaystyle \sum _{k=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}dt\delta a_k^{\left(1\right)}\left(t\right)\left\{a_{2,k}^{\left(2\right)}\left(t;j\right)-{\displaystyle \sum _{i=1}^{TH}{\phi }_k\left(F;t\right)dt{\displaystyle \underset{t_0}{\overset{t}{\int }}{a}_{2,i}^{\left(2\right)}\left(\tau ;j\right)\frac{\partial {\psi }_k\left[h\left(\tau \right);F;\tau \right]}{\partial h_i\left(\tau \right)}d\tau }}\right\}}}.\end{array}$$

(92)

The right-side of Eq. (92) can now be required to represent the indirect-effect term defined in Eq. (87), by imposing the requirement that the hitherto arbitrary function $a^{\left(2\right)}\left(t;j\right)={\left[a_1^{\left(2\right)}\left(t;j\right),a_2^{\left(2\right)}\left(t;j\right)\right]}^{†}$ be the solution of the following NIE-Volterra equations, for $i=1,\mathrm{...},TH;\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}j=1,\mathrm{...},TF$:

$$\begin{array}{l}{a}_{1,k}^{\left(2\right)}\left(t;j\right)-{\displaystyle \sum _{i=1}^{TH}\frac{\partial {\psi }_i\left[h\left(t\right);F;t\right]}{\partial h_k\left(t\right)}{\displaystyle \underset{t}{\overset{t_f}{\int }}{a}_{1,i}^{\left(2\right)}\left(\tau ;j\right){\phi }_i\left(F;\tau \right)d\tau }}-{\displaystyle \sum _{i=1}^{TH}\frac{{\partial }^{2}D\left[h\left(t\right);F;t\right]}{\partial h_k\left(t\right)\partial h_i\left(t\right)}{a}_{2,i}^{\left(2\right)}\left(t;j\right)}\\ -{\displaystyle \sum _{i=1}^{TH}{\displaystyle \sum _{n=1}^{TH}\frac{{\partial }^2{\psi }_n\left[h\left(t\right);F;t\right]}{\partial h_k\left(t\right)\partial h_i\left(t\right)}{a}_{2,i}^{\left(2\right)}\left(t\right){\displaystyle \underset{t}{\overset{t_f}{\int }}{a}_n^{\left(1\right)}\left(\tau \right){\phi }_n\left(F;\tau \right)d\tau }}}={\displaystyle \frac{{\partial }^{2}D\left[h\left(t\right);F\left(\theta \right);t\right]}{\partial h_k\left(t\right)\partial F_j}}\\ +{\displaystyle \sum _{i=1}^{TH}\frac{\partial {\psi }_i\left[h\left(t\right);F;t\right]}{\partial h_k\left(t\right)}{\displaystyle \underset{t}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(\tau \right)\frac{\partial {\phi }_i\left(F;\tau \right)}{\partial F_j}d\tau }}+{\displaystyle \sum _{i=1}^{TH}\frac{{\partial }^2{\psi }_i\left[h\left(t\right);F;t\right]}{\partial h_k\left(t\right)\partial F_j}{\displaystyle \underset{t}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(\tau \right){\phi }_i\left(F;\tau \right)d\tau };}\end{array}$$

(93)

$$\begin{array}{l}{a}_{2,k}^{\left(2\right)}\left(t;j\right)-{\displaystyle \sum _{i=1}^{TH}{\phi }_k\left(F;t\right)dt{\displaystyle \underset{t_0}{\overset{t}{\int }}{a}_{2,i}^{\left(2\right)}\left(\tau ;j\right)\frac{\partial {\psi }_k\left[h\left(\tau \right);F;\tau \right]}{\partial h_i\left(\tau \right)}d\tau }}\\ ={\displaystyle \frac{\partial g_k\left(F;t\right)}{\partial F_j}}+{\displaystyle \frac{\partial {\phi }_k\left(F;t\right)}{\partial F_j}}{\displaystyle \underset{t_0}{\overset{t}{\int }}{\psi }_k\left[h\left(\tau \right);F;\tau \right]d\tau }+{\phi }_k\left(F;t\right){\displaystyle \underset{t_0}{\overset{t}{\int }}\frac{\partial {\psi }_k}{\partial F_j}d\tau }.\end{array}$$

(94)

It follows from Eqs. (92)‒(94) that the indirect-effect term $\delta {\left(\partial R/\partial F_j\right)}_{ind}$ defined by Eq. (78) or, equivalently, Eq. (87) can be expressed in terms of the function $a^{\left(2\right)}\left(t;j\right)={\left[a_1^{\left(2\right)}\left(t;j\right),a_2^{\left(2\right)}\left(t;j\right)\right]}^{†}$ as follows, for $j=1,\mathrm{...},TF$:

$$\begin{array}{l}\delta {\left({\displaystyle \frac{\partial R}{\partial F_j}}\right)}_{ind}={\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_{1,i}^{\left(2\right)}\left(t;j\right)dt}}\left[{\displaystyle \sum _{n=1}^{TF}{q}_{in}\left(F;t\right)\delta F_n}\right]+{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_{2,i}^{\left(2\right)}\left(t;j\right)dt}}\left[{\displaystyle \sum _{n=1}^{TF}{S}_{in}\left(F;t\right)\delta F_n}\right]\\ ={\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_{1,i}^{\left(2\right)}\left(t;j\right)dt}}{\displaystyle \sum _{n=1}^{TF}\delta F_n}\left\{{\displaystyle \frac{\partial g_i\left(F;t\right)}{\partial F_n}}+{\displaystyle \frac{\partial {\phi }_i\left(F;t\right)}{\partial F_n}}{\displaystyle \underset{t_0}{\overset{t}{\int }}{\psi }_i\left[h\left(\tau \right);F;\tau \right]d\tau }\right.\\ \left.+{\phi }_i\left(F;t\right){\displaystyle \underset{t_0}{\overset{t}{\int }}\frac{\partial {\psi }_i\left[h\left(\tau \right);F;\tau \right]}{\partial F_n}d\tau }\right\}+{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_{2,i}^{\left(2\right)}\left(t;j\right)dt}}\left\{{\displaystyle \sum _{n=1}^{TF}\frac{{\partial }^{2}D\left[h\left(t\right);F;t\right]}{\partial F_n\partial h_i\left(t\right)}}\delta F_n\right.\\ +{\displaystyle \sum _{k=1}^{TH}\left[{\displaystyle \sum _{n=1}^{TF}\frac{{\partial }^2{\psi }_k\left[h\left(t\right);F;t\right]}{\partial F_n\partial h_i\left(t\right)}\delta F_n}\right]{\displaystyle \underset{t}{\overset{t_f}{\int }}{a}_k^{\left(1\right)}\left(\tau \right){\phi }_k\left(F;\tau \right)d\tau }}\\ \left.+{\displaystyle \sum _{k=1}^{TH}\frac{\partial {\psi }_k\left[h\left(t\right);F;t\right]}{\partial h_i\left(t\right)}{\displaystyle \sum _{n=1}^{TF}{\displaystyle \underset{t}{\overset{t_f}{\int }}\left[a_k^{\left(1\right)}\left(\tau \right)\frac{\partial {\phi }_k\left(F;\tau \right)}{\partial F_n}\delta F_n\right]d\tau }}}\right\}.\end{array}$$

(95)

The second-order sensitivities ${\partial }^{2}R\left[h;F\left(\theta \right)\right]/\partial F_j\partial F_n$ of the decoder-response with respect to the components of the feature function are obtained by adding the expression of the indirect-effect term obtained in Eq. (95) to the expression for the direct-effect term obtained in Eq. (77) and subsequently identifying the expressions that multiply the variations $\delta F_n$. The expressions thus obtained for ${\partial }^{2}R\left[h;F\left(\theta \right)\right]/\partial F_j\partial F_n$, for $j,n=1,\mathrm{...},TF$, are as follows:

$$\begin{array}{l}{\displaystyle \frac{{\partial }^{2}R\left[h;F\left(\theta \right)\right]}{\partial F_j\partial F_n}}\triangleq {\displaystyle \frac{{\partial }^{2}D\left[h\left(t\right);F\left(\theta \right);t\right]}{\partial F_n\partial F_j}}+{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(t\right)\frac{{\partial }^2{g}_i\left(h;F\right)}{\partial F_n\partial F_j}dt}}\\ +{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(t\right)\frac{{\partial }^2{\phi }_i\left(F;t\right)}{\partial F_n\partial F_j}{\displaystyle \underset{t_0}{\overset{t}{\int }}{\psi }_i\left[h\left(\tau \right);F;\tau \right]d\tau }}}+{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(t\right)\frac{\partial {\phi }_i\left(F;t\right)}{\partial F_j}{\displaystyle \underset{t_0}{\overset{t}{\int }}\frac{\partial {\psi }_i\left[h\left(\tau \right);F;\tau \right]}{\partial F_n}d\tau }}}\\ +{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(t\right)\frac{\partial {\phi }_i\left(F;t\right)}{\partial F_n}dt}}{\displaystyle \underset{t_0}{\overset{t}{\int }}\frac{\partial {\psi }_i\left[h\left(\tau \right);F;\tau \right]}{\partial F_j}d\tau }+{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_i^{\left(1\right)}\left(t\right){\phi }_i\left(F;t\right)dt}}{\displaystyle \frac{{\partial }^2{\psi }_i\left[h\left(\tau \right);F;\tau \right]}{\partial F_n\partial F_j}}\\ +{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_{1,i}^{\left(2\right)}\left(t;j\right)dt}}\left\{{\displaystyle \frac{\partial g_i\left(F;t\right)}{\partial F_n}}+{\displaystyle \frac{\partial {\phi }_i\left(F;t\right)}{\partial F_n}}{\displaystyle \underset{t_0}{\overset{t}{\int }}{\psi }_i\left[h\left(\tau \right);F;\tau \right]d\tau }\right.\\ \left.+{\phi }_i\left(F;t\right){\displaystyle \underset{t_0}{\overset{t}{\int }}\frac{\partial {\psi }_i\left[h\left(\tau \right);F;\tau \right]}{\partial F_n}d\tau }\right\}+{\displaystyle \sum _{i=1}^{TH}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}{a}_{2,i}^{\left(2\right)}\left(t;j\right)dt}}\left\{{\displaystyle \frac{{\partial }^{2}D\left[h\left(t\right);F;t\right]}{\partial F_n\partial h_i\left(t\right)}}\right.\\ +{\displaystyle \sum _{k=1}^{TH}\frac{{\partial }^2{\psi }_k\left[h\left(t\right);F;t\right]}{\partial F_n\partial h_i\left(t\right)}{\displaystyle \underset{t}{\overset{t_f}{\int }}{a}_k^{\left(1\right)}\left(\tau \right){\phi }_k\left(F;\tau \right)d\tau }}\left.+{\displaystyle \sum _{k=1}^{TH}\frac{\partial {\psi }_k\left[h\left(t\right);F;t\right]}{\partial h_i\left(t\right)}{\displaystyle \underset{t}{\overset{t_f}{\int }}{a}_k^{\left(1\right)}\left(\tau \right)\frac{\partial {\phi }_k\left(F;\tau \right)}{\partial F_n}d\tau }}\right\}.\end{array}$$

(96)

The NIE-Volterra system presented in Eqs. (93) and (94) is called the “2^nd-Level Adjoint Sensitivity System (2^nd-LASS)” and its solution, $a^{\left(2\right)}\left(j;t\right)={\left[a_1^{\left(2\right)}\left(j;t\right),a_2^{\left(2\right)}\left(j;t\right)\right]}^{†}$, $j=1,\mathrm{...},TF$, is called the “2^nd-level adjoint sensitivity function.” Since the sources on the right-sides of Eqs. (93) and (94) stem from the first-order sensitivities $\partial R\left[h;F\left(\theta \right)\right]/\partial F_j$, $j=1,\mathrm{...},TF$, they are dependent on the index “j”, which implies that for each first-order sensitivity $\partial R\left[h;F\left(\theta \right)\right]/\partial F_j$, there will correspond a distinct 2^nd-LASS, having a distinct solution $a^{\left(2\right)}\left(j;t\right)={\left[a_1^{\left(2\right)}\left(j;t\right),a_2^{\left(2\right)}\left(j;t\right)\right]}^{†}$, a fact that has been emphasized by using the index “j” in the list of arguments of this 2^nd-level adjoint sensitivity function. Therefore, there will be as many 2^nd-level adjoint functions as there are distinct first-order sensitivities $\partial R\left[h;F\left(\theta \right)\right]/\partial F_j$, which is equivalent to the number of components $F_j$of the “feature-function” $F\left(\theta \right)$. Notably, the integral operators on the left-sides of Eqs. (93) and (94) do not depend on the index “j”, which means that the same left-hand side needs to be inverted for computing the 2^nd-level adjoint function, regardless of the source term on the right-side (which corresponds to the particular component of the feature-function) of Eqs. (93) and (94). Therefore, if the inverses of the operators appearing on the left-sides of Eqs. (93) and (94) could be stored, they would not need to be inverted repeatedly, so the various 2^nd-level adjoint functions would be computed most efficiently.

The second-order sensitivities of the decoder-response with respect to the optimal weights/parameters ${\theta }_k,k=1,\mathrm{...},TW$, are obtained analytically by using the chain rule in conjunction with the expressions obtained in Eqs. (96) and (18), as follows:

$${\displaystyle \frac{{\partial }^{2}R\left[F\left(\theta \right)\right]}{\partial {\theta }_k\partial {\theta }_j}}={\displaystyle \frac{\partial }{\partial {\theta }_k}}\left\{{\displaystyle \sum _{i=1}^{TF}\frac{\partial R\left[F\left(\theta \right)\right]}{\partial F_i\left(\theta \right)}\frac{\partial F_i\left(\theta \right)}{\partial {\theta }_j}}\right\},j,k=1,\mathrm{...},TW.$$

(97)

When there are no feature functions but only individual model parameters, i.e., when $F_i\left(\theta \right)\equiv {\theta }_i$ for all $i=1,\mathrm{...},TF\triangleq TW$, the expression obtained in Eq. (96) yields directly the second-order sensitivities ${\partial }^{2}R/\partial {\theta }_i\partial {\theta }_j$, for all $i,j=1,\mathrm{...},TW$. In this case, the 2^nd-LASS would need to be solved $TW$-times rather than just $TF$-times $\left(TF<TW\right)$, when $TF$ feature functions can be constructed.

5. Illustrative Application of the 2^nd-FASAM-NIE-V to Neutron Slowing Down in an Infinite Homogeneous Hydrogenous Medium

The 2^nd-FASAM-NIE-V methodology developed in Section 4 will be applied in this Section to the illustrative model (considered in Section 3) that describes the energy distribution of neutrons in a homogeneous hydrogenous medium. As discussed in Section 4, the second-order sensitivities of the decoder response, $R\left[C\left(E\right)\right]$, with respect to the feature functions $F\left(\theta \right)$, ${\Sigma }_d\left(\theta \right)$, and S, will be determined by considering the second-order sensitivities to be “the first-order sensitivities of the first-order sensitivities.” Thus, the first-order sensitivities obtained in Eqs. (67)‒(69) will play the role of ”responses” in the application of the 2^nd-FASAM-NIE-V methodology.

5.1. Computation of Second-Order Sensitivities Stemming from $\partial R/\partial F\left(\theta \right)$

The second-order sensitivities stemming from the first-order sensitivity $\partial R/\partial F\left(\theta \right)$ will be determined from the expression of the first-order G-variation $\delta \left\{\partial R/\partial F\left(\theta \right)\right\}$, which is by definition obtained from Eq. (67) as shown below:

$$\begin{array}{l}\delta \left\{{\displaystyle \frac{\partial R}{\partial F\left(\theta \right)}}\right\}={\displaystyle \frac{d}{d\epsilon }}{\left\{{\displaystyle \frac{S^0+\epsilon \delta S}{E_s}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}\left[a^{\left(1,0\right)}\left(E\right)+\epsilon \delta a^{\left(1\right)}\left(E\right)\right]{\left(\frac{E_s}{E}\right)}^{\left[F\left({\theta }^0\right)+\epsilon \delta F\left(\theta \right)\right]}dE}\right\}}_{\epsilon =0}\\ =\delta {\left\{{\displaystyle \frac{\partial R}{\partial F\left(\theta \right)}}\right\}}_{dir}+\delta {\left\{{\displaystyle \frac{\partial R}{\partial F\left(\theta \right)}}\right\}}_{ind},\end{array}$$

(98)

where the direct-effect term $\delta {\left\{\partial R/\partial F\left(\theta \right)\right\}}_{dir}$ is defined as follows:

$$\begin{array}{l}\delta {\left\{{\displaystyle \frac{\partial R}{\partial F\left(\theta \right)}}\right\}}_{dir}\triangleq {\left\{{\displaystyle \frac{\left(\delta S\right)}{E_s}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right){\left(\frac{E_s}{E}\right)}^{F\left(\theta \right)}dE}\right\}}_{\theta ={\theta }^0}\\ +{\left\{\left[\delta F\left(\theta \right)\right]{\displaystyle \frac{S}{E_s}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right){\left(\frac{E_s}{E}\right)}^{F\left(\theta \right)}\ln \left(\frac{E_s}{E}\right)dE}\right\}}_{\theta ={\theta }^0},\end{array}$$

(99)

while the indirect-effect term $\delta {\left\{\partial R/\partial F\left(\theta \right)\right\}}_{ind}$ is defined as follows:

$$\delta {\left\{{\displaystyle \frac{\partial R}{\partial F\left(\theta \right)}}\right\}}_{ind}\triangleq {\left\{{\displaystyle \frac{S}{E_s}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}\delta a^{\left(1\right)}\left(E\right){\left(\frac{E_s}{E}\right)}^{F\left(\theta \right)}dE}\right\}}_{\theta ={\theta }^0}.$$

(100)

The direct-effect term can be computed at this time. On the other hand, the indirect-effect term can be computed after determining the variational function $\delta a^{\left(1\right)}\left(E\right)$, which is the solution of the G-differentiated 1^st-LASS defined in Eq. (50). By definition, the G-differential of Eq. (50) is provided by the following equation:

$$\begin{array}{l}{\displaystyle \frac{d}{d\epsilon }}{\left\{\left[a^{\left(1,0\right)}\left(E\right)+\epsilon \delta a^{\left(1\right)}\left(E\right)\right]\right\}}_{\epsilon =0}-{\displaystyle \frac{d}{d\epsilon }}{\left\{{\displaystyle \frac{F\left({\theta }^0\right)+\epsilon \delta F\left(\theta \right)}{E}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E}{\int }}\left[a^{\left(1,0\right)}\left(e\right)+\epsilon \delta a^{\left(1\right)}\left(e\right)\right]de}\right\}}_{\epsilon =0}\\ ={\displaystyle \frac{d}{d\epsilon }}{\left\{{\Sigma }_d\left({\theta }^0\right)+\epsilon \delta {\Sigma }_d\left(\theta \right)\right\}}_{\epsilon =0}.\end{array}$$

(101)

Performing the operations indicated in Eq. (101) yields the following equation to be satisfied by the function $\delta a^{\left(1\right)}\left(E\right)$ at the nominal parameter values (the superscript “zero” will be omitted for notational simplicity):

$$\delta a^{\left(1\right)}\left(E\right)-{\displaystyle \frac{F\left(\theta \right)}{E}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E}{\int }}\delta a^{\left(1\right)}\left(e\right)de}={\displaystyle \frac{\delta F\left(\theta \right)}{E}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E}{\int }}{a}^{\left(1\right)}\left(e\right)de}+\delta {\Sigma }_d\left(\theta \right).$$

(102)

Evidently, Eq. (102) would need to be solved repeatedly, for every parameter variation, in order to compute the function $\delta a^{\left(1\right)}\left(E\right)$ that would correspond to the respective parameter variation. As was shown in Section 4, the need for repeatedly solving Eq. (102)can be avoided by deriving an alternative expression for the indirect-effect term which does not involve $\delta a^{\left(1\right)}\left(E\right)$. This alternative expression is derived by introducing the 2^nd-LASS, the solution of which would replace the function $\delta a^{\left(1\right)}\left(E\right)$ in the alternative expression for the indirect-effect term. Notably, Eq. (102) is independent of variations in the forward function, $\delta C\left(E\right)$; therefore, the 2^nd-level adjoint sensitivity function will comprise a single component, denoted as $a^{\left(2\right)}\left(E;1\right)$, and the 2^nd-LASS will comprise a single equation for this component.

Following the principles of the 2^nd-FASAM-NIE-V, we use Eq. (45) to construct the inner product of Eq. (102) with a function $a^{\left(2\right)}\left(E;1\right)$, where the argument “1” indicates that this 2^nd-level adjoint sensitivity function will correspond to the first-order sensitivity $\partial R\left[C\left(E\right)\right]/\partial F\left(\theta \right)$, to obtain the following relation:

$$\begin{array}{l}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(2\right)}\left(E;1\right)\delta a^{\left(1\right)}\left(E\right)dE}-F\left(\theta \right){\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(2\right)}\left(E;1\right)}{\displaystyle \frac{dE}{E}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E}{\int }}\delta a^{\left(1\right)}\left(e\right)de}\\ =\delta F\left(\theta \right){\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(2\right)}\left(E;1\right)}{\displaystyle \frac{dE}{E}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E}{\int }}{a}^{\left(1\right)}\left(e\right)de}+\delta {\Sigma }_d\left(\theta \right){\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(2\right)}\left(E;1\right)dE}.\end{array}$$

(103)

The function $a^{\left(2\right)}\left(E;1\right)$ will be determined by requiring the left-side of Eq. (103) to represent the indirect-effect term defined in Eq. (100). For this purpose, the left-side will be recast using “integration by parts” into the following form:

$$\begin{array}{l}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(2\right)}\left(E;1\right)\delta a^{\left(1\right)}\left(E\right)dE}-F\left(\theta \right){\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(2\right)}\left(E;1\right)}{\displaystyle \frac{dE}{E}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E}{\int }}\delta a^{\left(1\right)}\left(e\right)de}\\ ={\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}\delta a^{\left(1\right)}\left(E\right)\left[a^{\left(2\right)}\left(E;1\right)-F\left(\theta \right){\displaystyle \underset{E}{\overset{E_s}{\int }}{a}^{\left(2\right)}\left(e;1\right)}\frac{de}{e}\right]dE}.\end{array}$$

(104)

Requiring the right-side of Eq. (104) to represent the indirect-effect term defined in Eq. (100) yields the following Volterra-type 2^nd-LASS to be satisfied by the 2^nd-level adjoint sensitivity function $a^{\left(2\right)}\left(E;1\right)$:

$$a^{\left(2\right)}\left(E;1\right)-F\left(\theta \right){\displaystyle \underset{E}{\overset{E_s}{\int }}{a}^{\left(2\right)}\left(e;1\right)}{\displaystyle \frac{de}{e}}={\displaystyle \frac{S}{E_s}}{\left({\displaystyle \frac{E_s}{E}}\right)}^{F\left(\theta \right)}.$$

(105)

The above 2^nd-LASS is to be satisfied at the nominal parameter values ${\theta }^0$. Notably, Eq. (105) is an “final-value problem” since the computation of $a^{\left(2\right)}\left(E;1\right)$ commences at the highest energy, where $a^{\left(2\right)}\left(E_s;1\right)=S/E_s$, and proceeds towards the lowest energy value, $E_{\mathcal{l}}$. For subsequent verification purposes, the closed-form explicit expression of $a^{\left(2\right)}\left(E;1\right)$ obtained by solving Eq.(105) is as follows:

$$a^{\left(2\right)}\left(E;1\right)=E^{-F\left(\theta \right)}\left[-SF\left(\theta \right)E_s^{F\left(\theta \right)-1}\ln E+SF\left(\theta \right)E_s^{F\left(\theta \right)-1}\ln E_s+S{E}_s^{F\left(\theta \right)-1}\right].$$

(106)

It follows from Eqs. (103)‒(105) that the indirect-effect term $\delta {\left\{\partial R/\partial F\left(\theta \right)\right\}}_{ind}$ is given by the following expression involving the 2^nd-level adjoint sensitivity function $a^{\left(2\right)}\left(E;1\right)$:

$$\delta {\left\{{\displaystyle \frac{\partial R}{\partial F\left(\theta \right)}}\right\}}_{ind}=\delta F\left(\theta \right){\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(2\right)}\left(E;1\right)}{\displaystyle \frac{dE}{E}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E}{\int }}{a}^{\left(1\right)}\left(e\right)de}+\delta {\Sigma }_d\left(\theta \right){\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(2\right)}\left(E;1\right)dE}.$$

(107)

Adding the expression for the indirect-effect term $\delta {\left\{\partial R/\partial F\left(\theta \right)\right\}}_{ind}$ obtained in Eq. (107) to the expression for the direct-effect term $\delta {\left\{\partial R/\partial F\left(\theta \right)\right\}}_{dir}$ obtained in Eq. (99) yields the following expression for the G-differential $\delta \left\{\partial R/\partial F\left(\theta \right)\right\}$:

$$\begin{array}{l}\delta \left\{{\displaystyle \frac{\partial R}{\partial F\left(\theta \right)}}\right\}=\delta F\left(\theta \right){\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(2\right)}\left(E;1\right)}{\displaystyle \frac{dE}{E}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E}{\int }}{a}^{\left(1\right)}\left(e\right)de}+\delta {\Sigma }_d\left(\theta \right){\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(2\right)}\left(E;1\right)dE}\\ +{\displaystyle \frac{\left(\delta S\right)}{E_s}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right){\left(\frac{E_s}{E}\right)}^{F\left(\theta \right)}dE}+\left[\delta F\left(\theta \right)\right]{\displaystyle \frac{S}{E_s}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right){\left(\frac{E_s}{E}\right)}^{F\left(\theta \right)}\ln \left(\frac{E_s}{E}\right)dE}.\end{array}$$

(108)

The above expression is to be evaluated at the nominal parameter values ${\theta }^0$.

It follows from the expression obtained in Eq. (108) that the second-order sensitivities (of the decoder response with respect to the feature functions) which stem from the first-order sensitivity $\partial R\left[C\left(E\right)\right]/\partial F\left(\theta \right)$ have the following expressions to be evaluated at the nominal parameter values ${\theta }^0$:

$${\displaystyle \frac{{\partial }^{2}R}{{\partial }^{2}F\left(\theta \right)}}={\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(2\right)}\left(E;1\right)}{\displaystyle \frac{dE}{E}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E}{\int }}{a}^{\left(1\right)}\left(e\right)de}+{\displaystyle \frac{S}{E_s}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right){\left(\frac{E_s}{E}\right)}^{F\left(\theta \right)}\ln \left(\frac{E_s}{E}\right)dE};$$

(109)

$${\displaystyle \frac{{\partial }^{2}R}{\partial {\Sigma }_d\left(\theta \right)\partial F\left(\theta \right)}}={\displaystyle \frac{1}{E_s}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(2\right)}\left(E;1\right)dE;}$$

(110)

$${\displaystyle \frac{{\partial }^{2}R}{\partial S\partial F\left(\theta \right)}}={\displaystyle \frac{1}{E_s}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right){\left(\frac{E_s}{E}\right)}^{F\left(\theta \right)}dE}.$$

(111)

Since the 1^st-level adjoint sensitivity function $a^{\left(1\right)}\left(E\right)$ is already available, the sensitivity ${\partial }^{2}R/\partial S\partial F\left(\theta \right)$ can be computed. The closed-form explicit expressions for the above second-order sensitivities are obtained by inserting the expression obtained in Eq. (106) for $a^{\left(2\right)}\left(E;1\right)$ and the expression obtained in Eq. (51) for $a^{\left(1\right)}\left(E\right)$, and performing the respective integrations. Carrying out these operations yields the following expressions:

$$\begin{array}{l}{\displaystyle \frac{{\partial }^{2}R}{\partial F\left(\theta \right)\partial F\left(\theta \right)}}=\left[S{\Sigma }_d\left(\theta \right)\right]\left\{-2{\left[1-F\left(\theta \right)\right]}^{-3}\left[1-{\left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)}^{F\left(\theta \right)-1}\right]\right.\\ \left.-2{\left[1-F\left(\theta \right)\right]}^{-2}{\left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)}^{F\left(\theta \right)-1}\ln \left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)-{\displaystyle \frac{F\left(\theta \right)}{1-F\left(\theta \right)}}{\left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)}^{F\left(\theta \right)-1}{\left[\ln \left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)\right]}^2\right\};\end{array}$$

(112)

$${\displaystyle \frac{{\partial }^{2}R}{\partial {\Sigma }_d\left(\theta \right)\partial F\left(\theta \right)}}={\displaystyle \frac{S}{1-F\left(\theta \right)}}\left\{{\displaystyle \frac{1}{1-F\left(\theta \right)}}\left[1-{\left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)}^{F\left(\theta \right)-1}\right]-F\left(\theta \right){\left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)}^{F\left(\theta \right)-1}\ln \left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)\right\};$$

(113)

$${\displaystyle \frac{{\partial }^{2}R}{\partial S\partial F\left(\theta \right)}}={\displaystyle \frac{{\Sigma }_d\left(\theta \right)}{1-F\left(\theta \right)}}\left\{{\displaystyle \frac{1}{1-F\left(\theta \right)}}\left[1-{\left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)}^{F\left(\theta \right)-1}\right]-F\left(\theta \right){\left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)}^{F\left(\theta \right)-1}\ln \left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)\right\}.$$

(114)

5.2. Computation of Second-Order Sensitivities Stemming from $\partial R/\partial {\Sigma }_d\left(\theta \right)$

The second-order sensitivities stemming from the first-order sensitivity $\partial R/\partial {\Sigma }_d\left(\theta \right)$ will be determined from the expression of the first-order G-variation $\delta \left\{\partial R/\partial {\Sigma }_d\left(\theta \right)\right\}$, which is by definition obtained from Eq. (68) as shown below:

$$\delta \left\{{\displaystyle \frac{\partial R}{\partial {\Sigma }_d\left(\theta \right)}}\right\}={\displaystyle \frac{d}{d\epsilon }}{\left\{{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}\left[C^0\left(E\right)+\epsilon \delta C\left(E\right)\right]dE}\right\}}_{\epsilon =0}={\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}\delta C\left(E\right)dE}.$$

(115)

Comparing Eq. (115) with Eq. (35), it becomes apparent that the following relation holds:

$$\delta \left\{\partial R/\partial {\Sigma }_d\left(\theta \right)\right\}={\left\{\delta R\left(C;\theta ;\delta C\right)\right\}}_{ind}/{\Sigma }_d\left(\theta \right)$$

(116)

Replacing the expression obtained for ${\left\{\delta R\left(C;\theta ;\delta C\right)\right\}}_{ind}$ in Eq. (65) into Eq. (116) yields the following expression:

$$\delta \left\{{\displaystyle \frac{\partial R}{\partial {\Sigma }_d\left(\theta \right)}}\right\}={\displaystyle \frac{F\left(\delta S\right)}{E_s{\Sigma }_d}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right)dE}+{\displaystyle \frac{S\left(\delta F\right)}{E_s{\Sigma }_d}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right)dE}{\left({\displaystyle \frac{E_s}{E}}\right)}^{F\left(\theta \right)}.$$

(117)

It follows from Eq. (117) that the second-order sensitivities that stem from $\partial R/\partial {\Sigma }_d$ have the following expressions:

$${\displaystyle \frac{{\partial }^{2}R}{\partial F\left(\theta \right)\partial {\Sigma }_d\left(\theta \right)}}={\displaystyle \frac{S}{E_s{\Sigma }_d\left(\theta \right)}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right)dE}{\left({\displaystyle \frac{E_s}{E}}\right)}^{F\left(\theta \right)};$$

(118)

$${\displaystyle \frac{{\partial }^{2}R}{\delta {\Sigma }_d\left(\theta \right)\partial {\Sigma }_d\left(\theta \right)}}=0;$$

(119)

$${\displaystyle \frac{{\partial }^{2}R}{\delta S\partial {\Sigma }_d\left(\theta \right)}}={\displaystyle \frac{F\left(\theta \right)}{E_s{\Sigma }_d\left(\theta \right)}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right)dE}.$$

(120)

The explicit closed-form expressions of the above second-order sensitivities are obtained by replacing the expression of $a^{\left(1\right)}\left(E\right)$ in Eqs. (118) and (120), respectively, and performing the respective integrations to obtain:

$${\displaystyle \frac{{\partial }^{2}R}{\partial F\left(\theta \right)\partial {\Sigma }_d\left(\theta \right)}}={\displaystyle \frac{S}{1-F\left(\theta \right)}}\left\{{\displaystyle \frac{1}{1-F\left(\theta \right)}}\left[1-{\left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)}^{F\left(\theta \right)-1}\right]-F\left(\theta \right){\left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)}^{F\left(\theta \right)-1}\ln \left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)\right\};$$

(121)

$${\displaystyle \frac{{\partial }^{2}R}{\delta S\partial {\Sigma }_d\left(\theta \right)}}={\displaystyle \frac{F\left(\theta \right)}{1-F\left(\theta \right)}}\left[1-{\left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)}^{F\left(\theta \right)-1}\right].$$

(122)

The expression for ${\partial }^{2}R/\partial F\left(\theta \right)\partial {\Sigma }_d\left(\theta \right)$ in Eq. (118) must be equivalent to the expression for ${\partial }^{2}R/\partial {\Sigma }_d\left(\theta \right)\partial F\left(\theta \right)$ in Eq. (110). Notably, therefore, these mixed second-order sensitivities are computed twice, using distinct expressions involving distinct adjoint functions, which provides an intrinsic mechanism for the stringent verification of the accuracy of the computations used for obtaining the numerical values of the respective adjoint functions.

5.3. Computation of Second-Order Sensitivities Stemming from $\partial R/\partial S$

The second-order sensitivities stemming from the first-order sensitivity $\partial R/\partial S$ will be determined from the expression of the first-order G-variation $\delta \left\{\partial R/\partial S\right\}$, which is by definition obtained from Eq. (69) as shown below:

$$\begin{array}{l}\delta \left\{{\displaystyle \frac{\partial R}{\partial S}}\right\}={\displaystyle \frac{d}{d\epsilon }}{\left\{{\displaystyle \frac{\left[F\left({\theta }^0\right)+\epsilon \delta F\left(\theta \right)\right]}{E_s}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}\left[a^{\left(1,0\right)}\left(E\right)+\epsilon \delta a^{\left(1\right)}\left(E\right)\right]dE}\right\}}_{\epsilon =0}\\ =\delta {\left\{{\displaystyle \frac{\partial R}{\partial S}}\right\}}_{dir}+\delta {\left\{{\displaystyle \frac{\partial R}{\partial S}}\right\}}_{ind},\end{array}$$

(123)

where the direct-effect term $\delta {\left\{\partial R/\partial S\right\}}_{dir}$ is defined as follows:

$$\delta {\left\{{\displaystyle \frac{\partial R}{\partial S}}\right\}}_{dir}\triangleq {\left\{{\displaystyle \frac{\delta F\left(\theta \right)}{E_s}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right)dE}\right\}}_{\theta ={\theta }^0},$$

(124)

while the indirect-effect term $\delta {\left\{\partial R/\partial S\right\}}_{ind}$ is defined as follows:

$$\delta {\left\{{\displaystyle \frac{\partial R}{\partial S}}\right\}}_{ind}\triangleq {\left\{{\displaystyle \frac{F\left(\theta \right)}{E_s}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}\delta a^{\left(1\right)}\left(E\right)dE}\right\}}_{\theta ={\theta }^0}.$$

(125)

The direct-effect term $\delta {\left\{\partial R/\partial S\right\}}_{dir}$ can be computed at this time. To determine the indirect-effect term $\delta {\left\{\partial R/\partial S\right\}}_{ind}$ without needing to compute the variational function $\delta a^{\left(1\right)}\left(E\right)$ by solving repeatedly Eq. (102), the steps outlined in Section 5.1 are applied using a 2^nd-order adjoint sensitivity function denoted as $a^{\left(2\right)}\left(E;3\right)$, where the argument “3” indicates that this function will correspond to the first-order sensitivity $\partial R/\partial S$. Thus, following the conceptual steps outlined in Eqs. (103)‒(107) but for the function $a^{\left(2\right)}\left(E;3\right)$ as the counterpart of the function $a^{\left(2\right)}\left(E;1\right)$, and for the indirect-effect term $\delta {\left\{\partial R/\partial S\right\}}_{ind}$ as the counterpart of the indirect-effect term $\delta {\left\{\partial R/\partial F\left(\theta \right)\right\}}_{ind}$, leads to the following expression for the indirect-effect term $\delta {\left\{\partial R/\partial S\right\}}_{ind}$:

$$\delta {\left\{{\displaystyle \frac{\partial R}{\partial S}}\right\}}_{ind}\triangleq {\left\{\delta F\left(\theta \right){\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(2\right)}\left(E;3\right)}{\displaystyle \frac{dE}{E}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E}{\int }}{a}^{\left(1\right)}\left(e\right)de}+\delta {\Sigma }_d\left(\theta \right){\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(2\right)}\left(E;3\right)dE}\right\}}_{\theta ={\theta }^0},$$

(126)

where the 2^nd-order adjoint sensitivity function $a^{\left(2\right)}\left(E;3\right)$ is the solution of the following Volterra-type 2^nd-LASS:

$$a^{\left(2\right)}\left(E;3\right)-F\left(\theta \right){\displaystyle \underset{E}{\overset{E_s}{\int }}{a}^{\left(2\right)}\left(e;3\right)}{\displaystyle \frac{de}{e}}={\displaystyle \frac{F\left(\theta \right)}{E_s}}.$$

(127)

For verification purposes, the explicit closed-form expression of the solution of Eq. (127) is provided below.

$$a^{\left(2\right)}\left(E;3\right)={\displaystyle \frac{F\left(\theta \right)}{E_s}}{\left({\displaystyle \frac{E_s}{E}}\right)}^{F\left(\theta \right)}.$$

(128)

Adding the expressions obtained in Eqs. (126) and (124) yields the following expression for the first G-variation $\delta \left\{\partial R/\partial S\right\}$:

$$\begin{array}{l}\delta \left\{{\displaystyle \frac{\partial R}{\partial S}}\right\}={\displaystyle \frac{\delta F\left(\theta \right)}{E_s}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right)dE}+\delta F\left(\theta \right){\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(2\right)}\left(E;3\right)}{\displaystyle \frac{dE}{E}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E}{\int }}{a}^{\left(1\right)}\left(e\right)de}\\ +\delta {\Sigma }_d\left(\theta \right){\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(2\right)}\left(E;3\right)dE}.\end{array}$$

(129)

It follows from Eq. (129) that the second-order sensitivities stemming from the first-order sensitivity $\partial R/\partial S$ have the following expressions:

$${\displaystyle \frac{{\partial }^{2}R}{\partial F\left(\theta \right)\partial S}}={\displaystyle \frac{1}{E_s}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(1\right)}\left(E\right)dE}+{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(2\right)}\left(E;3\right)}{\displaystyle \frac{dE}{E}}{\displaystyle \underset{E_{\mathcal{l}}}{\overset{E}{\int }}{a}^{\left(1\right)}\left(e\right)de}$$

(130)

$${\displaystyle \frac{{\partial }^{2}R}{\partial {\Sigma }_d\left(\theta \right)\partial S}}={\displaystyle \underset{E_{\mathcal{l}}}{\overset{E_s}{\int }}{a}^{\left(2\right)}\left(E;3\right)dE}.$$

(131)

$${\displaystyle \frac{{\partial }^{2}R}{\partial S\partial S}}=0.$$

(132)

Inserting the expressions of $a^{\left(2\right)}\left(E;3\right)$ and $a^{\left(1\right)}\left(E\right)$, respectively, into Eqs. (130) and (131), and performing the respective integrations, yields the following closed-form explicit expressions for the respective second-order sensitivities:

$${\displaystyle \frac{\partial R}{\partial F\left(\theta \right)\partial S}}={\displaystyle \frac{{\Sigma }_d\left(\theta \right)}{1-F\left(\theta \right)}}\left\{{\displaystyle \frac{1}{1-F\left(\theta \right)}}\left[1-{\left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)}^{F\left(\theta \right)-1}\right]-F\left(\theta \right){\left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)}^{F\left(\theta \right)-1}\ln \left({\displaystyle \frac{E_s}{E_{\mathcal{l}}}}\right)\right\};$$

(133)

$${\displaystyle \frac{{\partial }^{2}R}{\partial {\Sigma }_d\left(\theta \right)\partial S}}={\displaystyle \frac{F\left(\theta \right)}{1-F\left(\theta \right)}}\left[1-{\left({\displaystyle \frac{E_s}{E}}\right)}^{F\left(\theta \right)-1}\right].$$

(134)

The expression for ${\partial }^{2}R/\partial SF\left(\theta \right)$ in Eq. (111) must be equivalent to the expression for ${\partial }^{2}R/\partial F\left(\theta \right)\partial S$ in Eq. (130). The expression for ${\partial }^{2}R/\partial S\partial {\Sigma }_d\left(\theta \right)$ in Eq. (122) must be equivalent to the expression for ${\partial }^{2}R/\partial {\Sigma }_d\left(\theta \right)\partial S$ in Eq. (131). The equivalences of these corresponding expressions provide stringent verification criteria for the accuracy of the computation of the respective adjoint functions.

5.4. Discussion: Direct Computation of Second-Order Sensitivities Versus Their Indirect Computation via Feature Functions

Notably, the 9 second-order sensitivities of the decoder response with respect to the 3 feature functions,$F\left(\theta \right)$, ${\Sigma }_d\left(\theta \right)$ and $S$ were computed using 3 adjoint computations; each of these adjoint computations corresponds to one of the 3 first-order sensitivities of the decoder response with respect to the feature functions. Only 6 of these 9 second-order sensitivities have distinct values; the mixed sensitivities were computed twice, using different adjoint functions thus providing stringent verification criteria for the numerical computations of these functions. The second-order sensitivities of the decoder response with respect to the primary model parameters are obtained by applying the “chain-rule of differentiation” provided in Eq. (97) to the second-order sensitivities with respect to the feature functions.

On the other hand, the computation of the second-order sensitivities of the decoder response directly with respect to the primary model parameters would be performed by treating each of the first-order sensitivities defined in Eqs. (54)‒(57) as a “decoder/model response.” A shown in Section 3.1, there would be $TW+1$ first-order sensitivities in this case, which means that there would be $TW+1$ “2^nd-level adjoint sensitivity systems” to be solved, each one having a source-term that would correspond to one of the $TW+1$ first-order sensitivities. Evidently, it is considerably more advantageous computationally to consider compute the second-order sensitivities via “feature functions,” whenever possible, rather than directly with respect to the primary model parameters.

6. Discussion and Conclusions

This work has introduced the general mathematical framework of the 2^nd-FASAM-NIE-V methodology. The acronym “2^nd-FASAM-NIE-V” stands for “Second-Order Features Adjoint Sensitivity Analysis Methodology for Neural Integral Equations of Volterra Type.” The 2^nd-FASAM-NIE-V encompasses the mathematical framework of the (first-order) 1^st-FASAM-NIE-V methodology, which enables the most efficient computation of the exact expressions of all first-order sensitivities with respect to the feature functions and also with respect to the optimal values of the NIE-net’s parameters/weights, after the respective NIE-Volterra-net was optimized to represent the underlying physical system. The 1^st-FASAM-NIE-V methodology requires a single large-scale computation for determining the first-level adjoint sensitivity function that is subsequently used for computing the sensitivities using conventional numerical quadrature formulas. The 2^nd-FASAM-NIE-V requires as many large-scale computations (to solve the 2^nd-Level Adjoint Sensitivity System) as there are first-order sensitivities of the decoder response with respect to the feature functions. Subsequently, the second-order sensitivities of the decoder response with respect to the primary model parameters are obtained by applying the “chain-rule of differentiation” to the second-order sensitivities with respect to the feature functions.

The application of the 1^st-FASAM-NIE-V and the 2^nd-FASAM-NIE-V methodologies has been illustrated by using a well-known model for neutron slowing down in a homogeneous hydrogenous medium. This model has been chosen because the application of the 1^st-FASAM-NIE-V and the 2^nd-FASAM-NIE-V yields tractable explicit exact expressions for all quantities of interest, including the various adjoint sensitivity functions and first- and second-order sensitivities of the decoder response with respect to all feature functions and also with respect to the primary model parameters. This illustrative application highlights the unsurpassed efficiency of the 1^st-FASAM-NIE-V and the 2^nd-FASAM-NIE-V for second-order sensitivity analysis of NIE-Volterra nets. Ongoing research aims at developing the “Second-Order Features Adjoint Sensitivity Analysis Methodology for Neural Integro-Differential Equations,” which will enable, in premiere, the exact computations of second-order sensitivities of decoder responses with respect to optimized weights/parameters, based on the NIDE-models introduced in [22].