Second-Order Adjoint Sensitivity Analysis Methodology for Computing Exactly Response Sensitivities to Uncertain Parameters and Boundaries of Linear Systems: Mathematical Framework

This work presents the “Second-Order Comprehensive Adjoint Sensitivity Analysis Methodology (2 nd -CASAM)” for the efficient and exact computation of 1 st - and 2 nd -order response sensitivities to uncertain parameters and domain boundaries of linear systems. The model’s response (i.e., model result of interest) is a generic nonlinear function of the model’s forward and adjoint state functions, and also depends on the imprecisely known boundaries and model parameters. In the practically important particular case when the response is a scalar-valued functional of the forward and adjoint state functions characterizing a model comprising N parameters, the 2 nd -CASAM requires a single large-scale computation using the First-Level Adjoint Sensitivity System (1 st -LASS) for obtaining all of the first-order response sensitivities, and at most N large-scale computations using the Second-Level Adjoint Sensitivity System (2 nd -LASS) for obtaining exactly all of the second-order response sensitivities. In contradistinction, forward other methods would require (N 2 /2 + 3 N/2) large-scale computations for obtaining all of the first- and second-order sensitivities. This work also shows that constructing and solving the 2 nd -LASS requires very little additional effort beyond the construction of the 1 st -LASS needed for computing the first-order sensitivities. Solving the equations underlying the 1 st -LASS and 2 nd -LASS requires the same computational solvers as needed for solving (i.e., be used for solving the 1 st -LASS and the 2 nd -LASS. Since neither the 1 st -LASS nor the 2 nd -LASS involves any differentials of the operators underlying the original system, the 1 st -LASS is designated as a “first-level” (as opposed to a “first-order”) adjoint sensitivity system, while the 2 nd -LASS is designated as a “second-level” (rather than a “second-order”) adjoint sensitivity system. Mixed second-order response sensitivities involving boundary parameters may arise from all source terms of the 2 nd -LASS that involve the imprecisely known boundary parameters. Notably, the 2 nd -LASS encompasses an auto-matic, inherent, and independent “solution verification” mechanism of the correctness and accuracy of the 2 nd -level adjoint functions needed for the efficient and exact computation of the second-order sensitivities.


Introduction
The earliest use of adjoint operators for computing exactly and efficiently the first-order sensitivities of responses of a large-scale linear system comprising many parameters has appeared in the report by Wigner [2] on his work on the "nuclear pile" (i.e., nuclear reactor) using the linear neutron transport or diffusion equations. Cacuci [2] [3] conceived the rigorous 1 st -order adjoint sensitivity analysis methodology for generic large-scale nonlinear systems involving operator responses, comprising functional-type responses as particular cases. Representative developments of the adjoint sensitivity analysis methodology were reviewed in the books by Cacuci [4] and Cacuci, Ionescu-Bujor and Navon [5].
The first-order adjoint sensitivity analysis methodology was extended by Cacuci [6] [7] [8] to the "Second-Order Adjoint Sensitivity Analysis Methodology (2 nd -ASAM)" for linear and nonlinear systems, to enable the exact and efficient computation of all of the second-order sensitivities (i.e., functional derivatives) of model responses to parameters. However, all of the above works, as well as the vast majority of the application of the first-and second-order adjoint sensitivity analysis methodology considered that the phase-space location of the physical operators was subsequently used [10]- [15], either formally or in conjunction with variational formulations, to obtain approximate first-order sensitivities to boundary parameters of responses that were linear functionals (or ratios thereof) of the neutron flux in the context of linear neutron diffusion or neutron transport problems. Furthermore, the works referenced in the foregoing have considered the computation of response sensitivities to model parameters separately (rather than simultaneously) from the computations of response sensitivities to imprecisely known domain boundary parameters. Also, none of the above works considered responses that are simultaneously functions of the forward and adjoint state functions. Noteworthy, Cacuci [16] [17] has presented methodologies for the exact computation of 1 st -and 2 nd -order sensitivities of responses that are functionals of both the forward and the adjoint fluxes in a multiplying nuclear system with source. However, these works [16] [17] have specifically considered the neutron transport equation (rather than generic linear problems) within precisely known domains, without considering response sensitivities to uncertain domain parameters.
The aim of this work is to present the novel "Second-Order Comprehensive Adjoint Sensitivity Analysis Methodology (2 nd -CASAM)," which has the following features that generalize and extend all previously published works on this topic: 2) The system response considered within the 2 nd -CASAM framework is an operator-valued response that depends on both the forward and adjoint state-functions. Functional-valued responses are subsumed as particular cases.
2) The 2 nd -CASAM framework enables the efficient and exact computation of the 1 st -and 2 nd -order response sensitivities not only to uncertain model parameters but also to uncertain domain boundaries of generic linear systems.
This work is structured as follows: Section 2 presents the mathematical framework of the generic linear physical system comprising imprecisely known parameters and boundaries. Section 3 presents the mathematical framework of the 2 nd -CASAM. Section 4 offers concluding remarks. An accompanying work [18] presents an illustrative application of the general 2 nd -CASAM methodology to a paradigm model of generic particle/radiation transmission and/or evolution which has fundamental applications in many fields, including nuclear engineering (e.g., radiation detection, chemical reprocessing of spent reactor fuel, radioactive decay, etc.).
The quantities appearing in the Equation (1) are defined as follows: 2) α denotes a N α -dimensional column vector whose components are the physical system's imprecisely known parameters, which are subject to uncertainties; subset of the set of real scalars. The vector N α ∈  α is considered to include any imprecisely known model parameters that may enter into defining the system's boundary in the phase-se space of independent variables. The symbol "  " will be used to denote "is defined as" or "is by definition equal to." 2) ( ) 1 , , denotes the x N -dimensional phase-space position vector, defined on a phase-space domain denoted as x Ω which is defined as follows: . The lower boundary-point of an independent variable is denoted as x  α , where the superscript "  " denotes "lower" (e.g., the inner radius of a sphere or cylinder, the lower range of an energy-variable, etc.), while the corresponding upper boundary-point is denoted as x α , where the superscript "u" denotes "upper" (e.g., the outer radius of a sphere or cylinder, the upper range of an energy-variable, etc.). A typical example of "boundaries that depend on imprecisely known parameters" is provided by the boundary conditions needed for models based on diffusion theory, in which the respective "flux and/or current conditions" for the "boundaries facing vacuum" are imposed on the "extrapolated boundary" of the respective spatial domain. As is well known, the "extrapolated boundary" depends not only on the imprecisely known physical dimensions of the problem's domain, but also on the medium's microscopic transport cross sections and atomic number densities. For subsequent derivations, it is convenient to define the "vectors of endpoints" denotes a N ϕ -dimensional column vector whose components represent the system's dependent variables (also called "state functions");  -component column   vector whose components are operators (including differential, difference, integral, distributions, and/or infinite matrices) acting linearly on ϕ and nonlinearly on α .
6) All of the equalities in this work are considered to hold in the weak ("distributional") sense, since the right-sides ("sources") of Equation (1) and of other various equations to be derived in this work may contain distributions ("generalized functions/functionals"), particularly Dirac-distributions and derivatives and/or integrals thereof.
When differential operators appear in Equation (1), a corresponding set of boundary and/or initial conditions, which are essential to define the domain of ( ) , L x α , must also be given. Since this work considers only systems that are linear in the state function ( ) , the accompanying boundary and/or initial conditions must also be linear in , which means that they can be represented in operator form as follows: In Equation (2), the operator is a matrix comprising, as components, operators that act linearly on ( ) x ϕ and nonlinearly on α ; the quantity B N denotes the total number of boundary and initial conditions. The operator Equations (3) and (4) represent the "base-case" or nominal state of the physical system. Throughout this work, the superscript "0" will be used to denote "nominal" or "expected" values.
The vector-valued function ( ) is considered to be the unique nontrivial solution of the physical problem described by Equations (1) and (2). The linear operator ( ) , L x α is considered in this work to admit an adjoint operator, which will be denoted as : In Equation (5), the "product" notation compactly denotes the respective multiple integrals, while the dot indicates the "scalar product of two vectors" defined as follows: Formally, the inner product introduced in Equation (5)  ; ; In Equation (7), the formal adjoint operator ( ) is the solution of an adjoint system which can be written as follows: where the source for the adjoint equation is usually related to the system response under consideration. The domain of ( ) * , L x α is determined by selecting appropriate adjoint boundary and/or initial conditions, represented here in operator form as follows: In Equation (10), the subscript "A" indicates "adjoint," and the letter "B" indicates "boundary and/or initial conditions." The nominal value of the adjoint state function ( ) x ψ will be denoted as ( ) 0 x ψ , and is obtained by solving the adjoint system at the nominal parameter values 0 α : The system response (i.e., result of interest) associated with the problem modeled by Equations (1) and (2) is typically a real-valued nonlinear operator (function) of the system's forward and adjoint state-functions (i.e., dependent variables) and parameters, denoted as , which can be generally represented in the form where ( ) , ; ; F x ϕ ψ α is a suitably differentiable function with respect to its arguments.

The Second-Order Comprehensive Adjoint Sensitivity Analysis Methodology (2 nd -CASAM)
As has been already mentioned in the foregoing, the model and boundary para-     R ϕ ψ α is presented in Section 3.2, while the mathematical methodology for computing the 2 nd -order sensitivities is presented in Section 3.2.

Derivation of the 1 st -Level Adjoint Sensitivity System (1 st -LASS) for Computing Exactly and Efficiently the 1 st -Order Response Sensitivities to Model and Boundary Parameters
As shown by Cacuci [2] [3], in front of the sentence, the most general definition of the 1 st -order total sensitivity of the operator-valued model response , which is defined as follows: for ε ∈F , where F denotes the field of real scalars, and all (i.e., arbitrary) vectors ( ) ; The existence of the G-variation

( )
; F e x must satisfy the following condition for 1 2 , Numerical methods (e.g., Newton's method and variants thereof) for solving Equations (1) and (2) also require the existence of the 1 st -order G-derivatives of original model equations, in which case the components of the operators which appear in these equations must also satisfy the conditions described in Equations (15) and (16). Of course, if the first-order G-derivatives of the system's response do not exist, the computation of higher-order response sensitivities (G-derivatives) would be moot. Therefore, the conditions described in Equations (15) and (16) will henceforth be considered to be satisfied by the operators underlying the physical system, in which case the partial G-derivatives of ( ) R e at 0 e with respect to ϕ , ψ and α exist. These derivatives are row vectors defined as .
The notation { } ( ) 0 e indicates that the quantity within the braces is to be Equation (18)  (1) and (2). Thus, applying the definition of the G-differential to Equations (1) and (2) yields the following equations: where the superscript "(1)" indicates "1 st -Level" and the letter "B" indicates "boundary and/or initial conditions," and where the following definitions were used: The partial G-derivatives appearing in Equation (22) are matrices defined as follows: The partial G-derivatives that appear in Equation (23) are also matrices, having structures and components similar to those defined in Equations (24) and (25).
The first-order relationship between the variations ( ) δ x ψ and δα is determined by taking the G-differentials of Equations (9) and (10), which yields the following system of equations: where the superscript "(1)" indicates "1 st -Level," the letter "A" indicates "adjoint," the letter "B" indicates "boundary and/or initial conditions," and where the following definitions were used: As has been first shown by Cacuci [3], the adjoint sensitivity analysis methodology cannot be applied directly to an operator-valued response, but only to responses that are functionals (i.e., scalar-valued operators) of the state functions. For this purpose, the operator-valued response ( ) , ; R ϕ ψ α is expressed through its generalized Fourier (spectral) expansion ; , where the Fourier (spectral) coefficients ( ) ; F f e are functionals of the parameters, the forward and adjoint state variables, defined as follows: and where the quantities The 1 st -order G-differential of the response ( ) , ; ; ; , ; ; ; ; d ,  ; The appearance of the variations δϕ and δψ in the indirect-effect term is defined as follows: In particular, the inner product defined in Equation (37) also holds at the nominal parameter values 0 α .
2) Using Equation (37), construct the inner product of Equations (20) and (26) with a vector to obtain the following relation: , . ; ; which implies that the following system of equations is to be satisfied in the weak (distributional) sense: , , where the superscript "(1)" indicates "1 st -Level." The boundary and/or initial conditions represented by Equation (43)  ; , to a residual quantity that contains boundary terms involving only known values of δα , the index f and the nominal parameter values 0 α . This residual quantity will be denoted by ; ; P f δ α α , which is linear in δα and can be therefore written in the form .
In general, the residual quantity ( ) ( ) ; ; P f δ α α does not automatically vanish, although it may do so in particular instances. By considering extensions of ( ) ; , to the residual quantity ; ; P f δ α α , yields the following form for Equation (39): . ind ind Adding the results obtained in Equation (47)  , ; , where the quantity ( )   .

Derivation of the 2 nd -Level Adjoint Sensitivity System (2 nd -LASS) for Computing Exactly and Efficiently the 2 nd -Order Response Sensitivities to Model and Boundary Parameters
The 2 nd -order total sensitivity, denoted as ( ) ( ) F f ϕ ψ α is provided by the 2 nd -order G-differential of ( ) ; , ; ; , F δ e h h is defined inductively as "the total first-order G-variation (or 1 st -order total G-differential) of the 1 st -order total response variation (or 1 st -order G-differential)," as follows: ; ϕ ψ satisfy the conditions expressed in Equations (15) and (16) in a neighborhood e of 0 e . In order to proceed, it will be assumed that this is indeed the case. For the derivations to follow, it is convenient to introduce the vectors ; ; ; , and the variations around 0 g will be denoted as Using the definition provided in Equation (50), the 1 st -order G-differential For notational simplicity, the arguments of the functional   ; , , ; ; , , . ; , , ; ; ij  ; Since the source-terms of the 2 nd -LFSS depend on the parameter variations i δα , the 2 nd -LFSS is computationally expensive to solve. To avoid the need for solving the 2 nd -LFSS, the indirect-effect term defined in Equation (53) will be expressed in terms of a 2 nd -Level Adjoint Sensitivity System (2 nd -LASS), which will be constructed by following the general principles introduced by Cacuci [ ∈ b x H , is defined as follows: In particular, the inner product defined in Equation (60) with Equation (57) to obtain the following relation: .
The relation expressed by Equation (60) 3) The first term on the right-side of Equation (62) is now required to represent the indirect-effect term defined in Equation (53), which is accomplished by requiring that the following relation be satisfied: in Equation (62) to vanish. Even when it does not vanish, however, this bilinear concomitant will be reduced to a residual quantity, which will be denoted as 6) The second term on the right-side of Equation (66) is actually the indirect-effect term, (53), so that Equation (66) can be re-written in the following detailed form:  ; ; It is important to note that the 2 nd -LASS is independent of the variations δα  be required by re-computations and/or forward methods. It is also important to note that by solving the 2 nd -LASS N α -times, the "off-diagonal" 2 nd -order mixed sensitivities 1 2 2 i i F α α ∂ ∂ ∂ , 1 2 i i ≠ , will be computed twice, in two different ways (i.e., using distinct 2 nd -level adjoint functions), thereby providing an independent intrinsic (numerical) verification that the respective sensitivities are computed accurately.
Another important characteristic of using the 2 nd -LASS is the flexibility it provides for prioritizing the computation of the 2 nd -order sensitivities. The 2 nd -order mixed G-derivatives (sensitivities) corresponding to the largest relative 1 st -order response sensitivity should be computed first; the second largest relative 1 st -order response sensitivity should be considered next, and so on. Computing 2 nd -order partial sensitivities that correspond to vanishing 1 st -order sensitivities may also be of interest, since vanishing 1 st -order sensitivities may indicate critical points of the response in the phase-space of model parameters. Thus, only the 2 nd -order partial sensitivities which are deemed important should be computed; the unimportant ones can be deliberately neglected while knowing the error incurred by neglecting them.
It is noteworthy that the solving the equations underlying the 1 st -LASS and 2 nd -LASS require computational solvers for solving (i.e., "inverting") either the forward linear operator ( ) L α or the adjoint linear operator ( ) * L α . Only the various source on the right-sides of the 1 st -LASS and 2 nd -LASS differ from one another. Therefore, the same computer program and "solvers" can be used for solving the 1 st -LASS and the 2 nd -LASS as were already used for solving the original system of equations. For this reason, the 2 nd -LASS was designated as a "second-level" rather than a "second-order" adjoint sensitivity system, since the 2 nd -LASS does not involve any explicit 2 nd -order G-derivatives of the operators underlying the original system, but involves the inversion of the same operators that are needed to be inverted for solving the 1 st -LASS and the original system of equations.

Conclusions
This work has presented the "Second-Order Comprehensive Adjoint Sensitivity Analysis Methodology (2 nd -CASAM)" for the efficient and exact computation of 1 st -and 2 nd -order response sensitivities to imprecisely known parameters and domain boundary for a generic/general model for linear physical systems. The model's response (i.e., model result of interest) is considered to be a generic function/operator that depends on the model's forward and adjoint state functions, as well as on the imprecisely known boundary and model parameters. It has been shown that the novel 2 nd -CASAM has the following features that generalize and extend all previously published works on this topic: 1) The system response considered within the 2 nd -CASAM framework is an operator-valued response that depends on both the forward and adjoint state-functions. Functional-valued responses are subsumed as particular cases.
2) The 2 nd -CASAM framework enables the efficient and exact computation of the 1st-and 2 nd -order response sensitivities not only to uncertain model parameters but also to uncertain domain boundaries of generic linear systems. In particular, contributions to the first-order response sensitivities with respect to the specific parameters that characterize the imprecisely known model's domain boundary can arise solely from boundary terms or directly from the response if its definition includes the domain's boundary. On the other hand, mixed second-order response sensitivities involving boundary parameters may arise from all source terms (of the 2 nd -LASS) that involve the imprecisely known boundary parameters.
Companion works [18] [19] present illustrative applications of the 2 nd -CASAM to various responses for a paradigm model of generic particle/radiation transmission and/or evolution which has fundamental applications in many fields, including nuclear engineering (e.g., radiation detection, chemical reprocessing of spent reactor fuel, radioactive decay, etc.).

Conflicts of Interest
The author declares no conflicts of interest regarding the publication of this paper.