On the utilization of the adjoint method in microwave tomography

In microwave imaging, the adjoint method is widely used for the efficient calculation of the update direction, which is then used to update the unknown model parameter. However, the utilization and the formulation of the adjoint method differ significantly depending on the imaging scenario and the applied optimization algorithm. Because of the problem‐specific nature of the adjoint formulations, the dissimilarities between the adjoint calculations may be overlooked. Here, we have classified the adjoint method formulations into two groups: the direct and indirect methods. The direct method involves calculating the derivative of the cost function, whereas, in the indirect method, the derivative of the predicted data is calculated. In this review, the direct and indirect adjoint methods are presented, compared, and discussed. The formulations are explicitly derived using the two‐dimensional wave equation in frequency and time domains. Finite‐difference time‐domain simulations are conducted to show the different uses of the adjoint methods for both single source‐multiple receiver, and multiple transceiver scenarios. This study demonstrated that an appropriate adjoint method selection is significant to achieve improved computational efficiency for the applied optimization algorithm.


| INTRODUCTION
Microwave imaging is an imaging technique that relies on the contrast of the dielectric properties (permittivity and conductivity) between the target and its surrounding medium.][13][14][15][16] Microwave imaging data can be either qualitative or quantitative, depending on the reconstruction method. 17Qualitative methods, such as radar based methods, focus on the presence and location of the target, if any exists.On the other hand, quantitative methods (microwave tomography) aim to reconstruct the dielectric distribution by solving nonlinear scattering equations.Linear approximations (e.g., Born approximation) can be applied in low-contrast cases (i.e., the dielectric variation is relatively small).However, when there is a significant variation in dielectric properties inside the solution domain, the problem is treated as a non-linear inversion problem.In such cases, the misfit between measured and predicted data is mostly defined as Euclidean norm and minimized through an iterative process in which the dielectric distribution is updated.
Mainly, two different approaches can be followed to find the update direction with respect to a change in the unknown model parameter.In this study, these approaches are denoted as the direct and indirect methods: (1) The direct method involves calculating the derivative of the cost function, which is then used directly in gradientdescent-based algorithms, (2) the indirect method involves calculating the derivative of the predicted data, known as sensitivity or Jacobian, which is then used in cost function minimization, as in the Gauss-Newton (GN) method.The direct method does not perform an explicit computation of the Jacobian matrix, but the Jacobian is implicitly included in the calculations.On the other hand, the indirect method calculates the Jacobian matrix explicitly.Once the Jacobian matrix is obtained, the derivative of the cost function can also be determined.
In the study by Fang et al., 28 the cost function was expressed in the frequency domain, and GN iterative approach was utilized for the inversion.A scalar three-dimensional (3-D) finite-element method was used as a forward solver.The Jacobian matrix for the dual-mesh configuration was constructed using the adjoint formulation and provided in terms of field mesh elements.The proposed approach was later generalized and introduced as a nodal adjoint method, allowing the Jacobian matrix to be easily computed, especially for the Finite-Difference Time-Domain (FDTD) solution. 29In a recent study, 32 a discrete dipole approximation-based solver was combined with the nodal adjoint technique.Using the same simulation grid for both the forward and imaging domain discretizations, the Jacobian calculations were reduced to a simple multiplication of the two forward solutions.In Rubaek et al., 33 the nonlinear optimization problem was solved by calculating the update values in the GN algorithm using the conjugate gradient least squares method.The proposed algorithm operated on the Jacobian matrix, which was obtained by incorporating the adjoint technique and resulting in a set of inner product operations.In the work of Islam et al., 35 the optimization problem was solved via a fast implementation of the GN method, 51 and the Jacobian matrix was calculated for fraction parameters using the adjoint method.In another study, 4 the elements of the Jacobian matrix for a coaxial antenna model were derived using the reciprocity principle for adjoint calculations.The elements were then expressed in terms of scattering parameters.
The adjoint method, in addition to the Jacobian calculations outlined above, is also utilized for gradient-based optimization algorithms.El-Shenawee et al. 36 introduced a shape reconstruction inversion algorithm that uses an adjoint scheme.The gradient directions for the cost function were obtained by solving the forward problem twice in each iteration, independent of the number of model parameters.The derivation of the adjoint calculations was based on the methodology described in Dorn et al.'s work. 37Furthermore, Tournier et al. 38 presented another gradient-based optimization algorithm to minimize the cost function through the solution of the state and adjoint problems.The derivations for the adjoint problem were given using frequency domain Maxwell equations and impedance boundary conditions.
Even though frequency-domain reconstructions are more common in the literature, time-domain inversion methods have also been utilized due to their ability to use a broad spectrum of frequencies.Most of the time domain inversion methods have used a gradient-based optimization referred to as the forward-backward time-stepping method (FBTS), in which the received residual signals are reversed in time and utilized as sources. 47,52It is worth noting that the FBTS method has found applications in acoustics 19 and geophysics, 25 where it is called full waveform inversion.In a recent study by Saraskanroud et al., 43 FBTS was used to minimize the time-domain cost function.The gradient coefficients for the relative permittivity and conductivity were calculated using the adjoint method and expressed in terms of discontinuous Galerkin method coefficients. 53In Rydholm et al.'s study, 54 a time-domain algorithm 44 was used to obtain gradients with respect to Debye parameters through a first-order perturbation analysis using both the original and adjoint Maxwell equations.In Johnson et al, 55 the gradient calculations were carried out for 3-D problems using the adjoint approach, 47,56 in which the expressions for the update direction using the adjoint method were provided using operators and vectorial formalism.Moreover, Rekanos et al. estimated, 57 the dielectric distribution inside the imaging domain by minimizing the augmented cost functional.The gradients were obtained via the adjoint method, and the derivations were provided in terms of Lagrange multipliers.
In the aforementioned studies, either the direct or the indirect methods have been used depending on the preferred algorithm.Although the calculations differed from each other in terms of the number of forward simulation runs and the utilized adjoint sources, the method was mostly denoted as an adjoint method.However, the dissimilarities between different adjoint calculations have remained disclosed since the problems have been addressed specifically to the type of the forward models (i.e., frequency or time domain methods) and applied optimization techniques (e.g., Gauss-Newton or gradient descent).This diversity of the adjoint calculations and complex formulations may cause key points to be missed and could make the topic confusing.
Efforts have been made to enhance the comprehensibility of the adjoint method in several areas, including acoustic scattering problems, 58 electrical resistance tomography, 59 general boundary value problems (acoustic, electromagnetic, and elastic), 60 heat conductivity identification, 61 and system sensitivity analysis. 62To the best of the authors' knowledge, there is currently no existing study in the microwave imaging literature that presents, compares, and discusses adjoint calculations in a straightforward manner.We believe that presenting the analytical calculation of the adjoint formulations, as well as implementation differences, in a compact, simple, and comparative manner will help the researchers understand underlying mathematical principles and finally lead them to the selection of efficient solving methodologies.
This study contributes in the following ways: (1) Detailed calculations for the derivative of the cost function and the Jacobian are presented using the adjoint method.For simplicity of the notations and clarity of the derivations, two-dimensional (2D) transverse magnetic (TM) illumination was used.(2) Regardless of the applied forward solver, a general solution methodology is given in both time and frequency domains.(3) It is demonstrated that the direct and indirect method calculations differ from each other in terms of the adjoint sources and the number of required forward problems.(4) The flow diagrams of the derivative calculations, along with the number of forward simulation runs, are provided for two simulation scenarios to show variances in terms of implementations and computational efficiencies.
The remainder of this paper is organized as follows: In Section 2, the analytical derivations for computing the update directions in both frequency and time domains, and the numerical models used in simulations are presented.Next, in Section 3, the analytical expressions calculated in the previous section are verified via simulations, and flow diagrams for derivative calculations in two simulation scenarios are provided with computational loads.In Section 4, the results are presented, and subsequently, in Section 5, the computational aspects for specific imaging configurations, efficient calculation of the update direction, comparison between frequency and time domain adjoint method, and a brief explanation of the applications of the adjoint method in other inverse problems are discussed.

| METHOD
In this section, the derivative calculations and numerical simulations are described.The derivative calculations are divided into two categories: frequency and time domain formulations, each with two methods -direct and indirect.In the direct method, the derivative of the cost function with respect to relative permittivity distribution was calculated directly.On the other hand, in the indirect method, the derivative of the calculated electric field with respect to relative permittivity distribution (also known as sensitivity or Jacobian in discrete form) was calculated.The calculations for the derivative of the cost function were adapted from the work by Norton 58 and repeated here for completeness, while the sensitivity calculations in the frequency domain followed the methods outlined in Aghasi et al., 59 which formulates the problem for the Poisson equation.We also present numerical simulation models for two different scenarios involving N antennas.In the first scenario, only one antenna transmits while the others serve as receivers, whereas in the second scenario, each antenna acts as both transmitter and receiver.

| Formulation in the Frequency Domain
For the transverse magnetic-to-z (TM-z) mode, the equation satisfied by the z-component of the electric field, E z , in a lossless, isotropic, non-magnetic medium can be written as is the Laplacian operator in the transverse plane, k 2 0 ¼ ω 2 ε 0 μ 0 is free-space wave number, ω is the angular frequency, ε 0 is the free space permittivity, μ 0 is the free space permeability, ε r r ð Þ is the relative permittivity distribution, s r ð Þ ¼ jωμ 0 J z r ð Þ is source, and J z is the z-directed electric current density.Considering an imaging scenario with a single transmitter N tx ¼ 1 ð Þfor simplicity and N rx receivers, the cost function can be defined as where are the calculated and measured complex-valued electric fields at r i , respectively.The derivative of the cost function with respect to the relative permittivity distribution ε r ð Þ can be expressed as where Ã denotes the complex conjugate.

| The Direct Method
Here, we show that the derivative of the cost function ∂F=∂ε r ð Þcan be calculated directly using the adjoint method.The aim of this approach is to determine the term ∂F=∂ε r that links the variation in the relative permittivity δε r ð Þ with the modification in the cost function δF ð Þ, which is in the form of for an imaging domain of S. Consider an infinitesimal change in the relative permittivity, ε r !ε r þ δε r , that causes a small change in the electric field, E z !E z þ δE z .Inserting these variations in (1) and neglecting the small order terms The position dependence of the relative permittivity ε r ð Þ was omitted for brevity.Inserting the variation in the cost function, Now, consider an adjoint field Ẽz as a solution to the wave equation where the adjoint source was defined as Note that the expression in (9) appears in the integral equation given in (7).By inserting the left-hand side of ( 8) to (7), and omitting the position dependence in the notations, the expression in (7) becomes Using the relation the first integral in ( 10) can be written as Expanding the integration domain to infinity, the first term at the right-hand side of (11) vanishes by using the divergence theorem.Note that the actual electric field E z ð Þ and the adjoint field Ẽz À Á satisfy the Sommerfield radiation conditions 63,64 and the surface integral vanishes when j r j! ∞.As a result, equation (10) can be simplified to Inserting r 2 δE z term from ( 5), ( 13) can be expressed as Using the equivalence between ( 14) and ( 4), the derivative of the cost function can be expressed as This result indicates that independent of the number of receiver antennas N rx ð Þ, the derivative of the cost function can be obtained by solving two forward problems.The first forward problem involves the calculation of the electric field E z ð Þ using (1), while the second forward problem (8) involves the calculation of the adjoint electric field Ẽz À Á .For the second forward problem, the adjoint source is defined in (9), where the source corresponds to the sum of residual fields (i.e., the difference between the calculated and measured electric field) at measurement points.In Aghasi et al.'s study, 59 this source is denoted as the composite-source.Similar terminology is adopted here, with a slight modification, and the source is referred to as the composite-adjoint-source.

| The Indirect Method
Unlike the direct method, the aim here is to find sensitivity at receiver locations ∂E z =∂ε r j r¼r i , rather than finding ∂F=∂ε r directly.The relation between the changes in the relative permittivity δε r ð Þ to the change in the calculated field at measurement points δE z j r¼r i can be formulated as Left hand side of ( 16) can be expressed as Now, consider an adjoint field Ẽz i as a solution to the wave equation where we define adjoint-source at the i th receiver as Þappears in the integral equation given in (17).By inserting the left hand side of ( 18) into (17), and following the steps from ( 10) to ( 14) by replacing Ẽz with Ẽz i , we obtain This result shows that the sensitivity at a point r i can be calculated by solving the forward problem (1) to get E z and the adjoint problem (18) to get Ẽz i for a point source located at r i .
By computing the residual field (i.e., the difference between the calculated and measured electric field), and the sensitivities for each receiver (corresponding to N rx forward solutions) the derivative of the cost function defined in (3) can be obtained.

| Formulation in Time Domain
The wave equation in the time domain for a lossless, isotropic, non-magnetic medium can be written as and the initial conditions at t ¼ 0 can be defined as where c 0 ¼ 1= ffiffiffiffiffiffiffiffi ffi μ 0 ε 0 p is the speed of light in free space, and the source is s r,t ð Þ¼μ 0 ∂J z r, t ð Þ=∂t ð Þ .Similar to the frequency domain case, a single transmitter N tx ¼ 1 ð Þand N rx receivers were utilized.Using the E z field measured at r i over a time interval 0, T ½ , the cost function is defined as and the derivative of the cost function is expressed as

| The Direct Method
Similar to the frequency domain solution, we will seek a formulation in the form of (4).Introducing an infinitesimal variation in the relative permittivity and the electric field, the wave equation ( 21) and the cost function ( 23) becomes The wave equation for an adjoint field Ẽz r, t ð Þ can be defined in terms of the reverse-time variable τ, where τ ¼ T À t: with the following initial conditions at τ ¼ 0 (corresponds to terminal conditions at t ¼ T): The adjoint-source s r, τ ð Þ is defined as Using ( 27) and ( 29), and omitting the position dependence in the notations, ( 26) can be expressed as Equation ( 30) can be divided into two integrals: one including the second-order spatial-derivative of the Ẽz , and the other one including the second-order time-derivative of the Ẽz .By following the steps given in (11) and (12), the first integral can be rearranged as follows For the second integral, the derivative with respect to reverse time variable τ is formulated in terms of the time variable is applied twice using the initial and terminal conditions given in ( 22) and (28).In the first one, u and v are defined as u ¼ δE z t ð Þ and v ¼ ∂ Ẽz τ ð Þ=∂t, and in the second one, u ¼ ∂δE z =∂t and v ¼ Ẽz τ ð Þ.Consequently, the second integral becomes Therefore, (30) can be written as Inserting r 2 δE z t ð Þ term from (25), δF ε r ð Þ can be expressed as Again, integrating by parts in the time domain Finally, comparing (35) and ( 4), and inserting τ ¼ T À t, the derivative of the cost function is found as Similar to the frequency domain results, the derivative of the cost function can be calculated by solving two forward problems: E z r, t ð Þ and adjoint field Ẽz r, τ ð Þ can be obtained by solving (21) and ( 27) using the sources s r, t ð Þ and s r, τ ð Þ, respectively.The adjoint-source s r, τ ð Þ given in (29) corresponds to the sum of received residual signals.To obtain the adjoint electric field Ẽz r, τ ð Þ, this source needs to be propagated backward in time from t ¼ T to t ¼ 0. Instead of running the model backward in time, the residual field can be reversed in the time domain and backpropagated. 19,65

| The Indirect Method
The time domain sensitivity at time t ¼ t k due to a point source at measurement point r i can be expressed as Similar to the frequency domain derivations, the sensitivity is computed explicitly.For a given transmitter-receiver configuration, the perturbation in the measured data due to the perturbation in the relative permittivity can be expressed as Left hand side of (38) can be written as Now, an adjoint electric field Ẽz i r, τ ð Þ is defined as a solution to the wave equation, where τ ¼ T À t: with initial conditions at τ ¼ 0 (terminal conditions at t ¼ T) imposed as The time-dependent adjoint source is defined as where the source is located at r ¼ r i .By multiplying both sides of (39) with h t ð Þ and integrating in time, we arrive at Substituting the left hand side of ( 40) into (43) gives which is in similar form of (30).Replacing Ẽz with Ẽz i , and following the steps (31) to (35), we obtain , the sensitivity at a specific time point t ¼ t k can be expressed as This result implies that sensitivity at a space point r ¼ r i and a time point t ¼ t k can be calculated by solving the forward problem (21), and the adjoint problem (40).Depending on the value of t k , forward and backward propagating fields are cross-correlated. 66][69]

| Numerical Simulations
For numerical simulations, TM-mode 2D FDTD method was utilized.The computational domain was terminated with a 12-grid perfectly matched layer (PML) in all directions.The source was implemented as a sinusoidal carrier with a Gaussian envelope.The operation frequency was 5 GHz and N ¼ 16 antennas were placed in a uniform circular arrangement.The geometry utilized in simulations is depicted in Figure 1.Excluding the PMLs, the size of the computational domain was 600 mm with a grid size of Δx ≤ λ min =20, where λ min ¼ λ 0 = ffiffiffiffiffiffiffiffiffiffiffi ffi ε r,max p is the minimum wavelength, and ε r, max is the maximum value of the relative permittivity.The imaging area was a square region with a side length of 200 mm, and antennas were placed on a radius of 210 mm.The received field was sampled at 2000 time points with time step Δt ≤ Δx= 2c 0 ð Þ.For frequency domain simulations, the received time domain electric field was transformed to the frequency domain, and the field component at 5 GHz was used.Two simulation scenarios were run: (1) single source and (NÀ1) receivers, and (2) N sources/receivers.The simulations were conducted in both the frequency domain and time domain.

| Case-1: Single source and (N À 1) receivers
In this scenario, one of the antennas was used as the source, and the remaining N À 1 ð Þ antennas were in receiving mode.In contrast to typical microwave imaging setups where antennas serve as both transmitters and receivers, this particular scenario has been chosen for several reasons.First, it allows for the validation of the derived formulations with minimal computational complexity.Second, it helps to gain a clear understanding of the implementation of the adjoint method before transitioning to a more realistic imaging configuration.Last, it demonstrates the computational effectiveness of the direct method when the number of sources is less than the number of receivers.

Frequency Domain Simulations
The derivative of the cost function for M grid cells, a single transmitter N tx ¼ 1 ð Þ, and N rx ¼ N À 1 receiving antennas can be expressed in the frequency domain as where f ℝ M is the derivative of the cost function (which corresponds to ∂F=∂ε r ), J C NÀ1 ð ÞÂM ð Þ is the complex-valued Jacobian matrix (i.e., sensitivity, ∂E z =∂ε r ) and e C NÀ1 ð Þ is the complex-valued residual vector, that is, the difference between the calculated and measured data.Note that (47) corresponds to the vector domain representation of (3).
To calculate the derivative of the cost function f ð Þ, three simulations were conducted in the frequency domain, and their results were compared: (1) finite-difference (FD) method, (2) direct method, and (3) indirect method.For all simulations, the initial model was selected as an empty background with a relative permittivity of ε r r ð Þ ¼ 1.
F I G U R E 1 Simulation model.The gray shaded area shows the perfectly matched layer (PML) region.A dielectric object was inserted at the center of the imaging domain.The transmitting and receiving antennas were located at the circle with a radius of 210 mm.The square imaging area has a side length of 200 mm.(Note that the drawing is not to scale and is for illustrative purposes only.) In the first simulation, the FD method was used to calculate the Jacobian matrix.The sensitivity for the j th cell in the grid was calculated as where δε j is the change in permittivity at j th grid cell.A single simulation was run for the initial model, and the received field at each cell in the grid was calculated.Then, the relative permittivity of a single cell was perturbed by 10%, and the simulation was run again to calculate the E z at receiver points.The amount of perturbation was selected to be small enough to ensure linear approximation, but large enough to avoid numerical errors.A perturbation amplitude of 10% was used since it remains in the validity domain of the linear approximation. 70Using (48) at receiving points r ¼ r i ð Þ gives the j th column of the Jacobian matrix.Simulations were repeated for each grid cell resulting in M simulations.After the calculation of the Jacobian matrix, the derivative vector f ð Þ was obtained using (47).In the second simulation, the direct method was used to calculate the derivative of the cost function directly instead of computing the Jacobian matrix explicitly.E z in the simulation domain was calculated using the initial model, and then, the residual data e i ¼ E z r i ð ÞÀE meas z r i ð Þ was calculated at each receiver position.To calculate the adjoint field Ẽz , the residual data at all receiver points was used as the composite-adjoint-source.Subsequently, the entries of the derivative of the cost function f ð Þ were calculated by using where S j corresponds to the area of the j th grid cell.Note that only two simulations were run to obtain the derivative of the cost function, and the method is independent of the number of grid cells and the number of receivers.
In the third simulation, the indirect method was used to calculate the entries of the Jacobian matrix.A single simulation was run for the initial model to calculate the E z component at each grid cell.Then, i th receiver was used as an adjoint source to obtain the adjoint field Ẽz i .i th row of the Jacobian matrix was calculated by evaluating Since the field, values were calculated at all grid cells, multiplication of the real and adjoint fields provided the i th row of the Jacobian matrix.The simulations were repeated for each receiver, resulting in N À 1 ð Þsimulations.Then, the computed Jacobian J ð Þ was utilized in (47) to calculate the derivative of the cost function f ð Þ.

Time Domain Simulations
Similar to the frequency domain case, the derivative of the cost function for M grid cells, a single transmitter N tx ¼ 1 ð Þ, and N rx ¼ N À 1 receiving antennas, and N t time points can be expressed in the time domain as where f ℝ M is the derivative of the cost function (which corresponds to ∂F=∂ε r ), J ℝ NÀ1 Þ is the Jacobian matrix (i.e., sensitivity, ∂E z r, t ð Þ=∂ε r ) and e ℝ NÀ1 ð ÞÂN t ½ is the residual vector, that is, the difference between the calculated and measured data.Note that (51) corresponds to the vector domain representation of (24).
Three simulations were conducted to calculate the derivative of the cost function: (1) FD method, (2) direct method, and (3) indirect method.The initial model was selected as an empty background with a relative permittivity of ε r r ð Þ ¼ 1.In the first simulation, the FD method was used.Time domain electric fields at receiver points were calculated for the initial and perturbed model.As in the frequency domain case, the relative permittivity at the j th grid cell was changed 10% to ensure linearity.The Jacobian matrix was calculated using (48) and repeating the simulations for M grid cells.Then, the functional derivative was obtained using (51).
In the second simulation, the derivative of the cost function was determined through the direct method rather than explicitly calculating the Jacobian matrix.First, E z r, t ð Þ was calculated, and then, the adjoint field Ẽz r, t ð Þ, was calculated by reversing and backpropagating the composite-adjoint-source (29) in the time domain.The derivative of the cost function was calculated by using (36).Analogous to the frequency domain, only two simulations were run to obtain the derivative of the cost function, and the method is independent of the number of grid cells and the number of receivers.
In the third simulation, the indirect method was used to calculate the functional derivative.First, the E z r, t ð Þ field was obtained at each grid cell.Then, the adjoint source was applied at each receiver to obtain the adjoint field Ẽz i r,t ð Þ.The Jacobian was calculated by convolving the actual and adjoint fields, and then, (51) was used to calculate the derivative of the cost function (f ).Note that the residual (e) was convolved with the source function used for the adjoint field before computing (51). 71,72A total of N À 1 ð Þforward simulations, which are independent of the number of grid cells, were required to compute the Jacobian matrix.

Simulation Setup
To compare the derivative of the cost function obtained with the aforementioned methods, a cylindrical dielectric object with a radius of 24 mm and the highest relative permittivity of ε r,max ¼ 2 was placed at the center of the imaging domain (refer to Figure 1).The relative permittivity distribution was where r ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi is the magnitude of the position vector r ¼ x,y ð Þ, and a is the radius of the object.The maps for the derivative of the cost function were calculated by taking the absolute value of the derivative of the cost function ( ∂F=∂ε r j j , or f j j), and normalizing it to its maximum value.In terms of implementations, we would like to point out that the source term in (1) includes jωμ 0 term whereas the adjoint-sources in ( 9) and ( 42) do not.Therefore, in cases where the same forward solver is used for both field and adjoint-field calculations, the sources need to be configured accurately for precise computations.The same is also valid for the time-domain simulations since the adjoint source does not include the time-derivative ∂=∂t ð Þ of the current source.
The flow diagrams for the calculation of the derivative of the cost function, together with the required number of forward simulation runs (FDTD runs), are presented in Figure 2.
F I G U R E 2 Flow diagrams for the calculation of the derivative of the cost functionsingle source and N À 1 ð Þreceivers case.s is the number of forward simulation runs (FDTD runs), M is the number of grid cells, and N is the number of antennas where M ) N. The total number of forward simulation runs was provided at the bottom of the diagrams.

| Case-2: N Sources/Receivers
In this scenario, a realistic microwave imaging configuration is considered in which the antennas were used as both transmitter and receiver.Each antenna sequentially transmits while the response is measured on the remaining N rx ¼ N À 1 antennas.For N tx ¼ N transmitter antennas, the cost function was defined as The flow diagrams for the calculation of the derivative of the cost function are shown in Figure 3, along with the number of forward simulation runs.
For the FD method, the initial and perturbed permittivity distributions were used to calculate the received fields, and then this process was repeated for each transmitter.Similarly, for the direct method, the received field due to the actual and the composite-adjoint sources was calculated for each transmitter.For the indirect method, on the other hand, additional calculations were not required as the transmitter antenna changes.Due to reciprocity, the previously calculated adjoint fields, obtained using point sources at the receivers, can also be utilized as the electric fields generated by the transmit antennas.Likewise, this approach can be applied in cases where field values can be acquired through the utilization of the linearity and scaling.
The maps for the derivative of the cost function were computed for two imaging configurations using the indirect method.The calculation of the functional derivative is performed in both frequency and time domains.For the first configuration, the same dielectric object (radius: 24 mm, ε r,max ¼ 2) was used as in the case of a single transmit antenna (See Figure 4A.In the second configuration, the derivative map was computed where two dielectric cylinders with diameters of 15 and 21 mm were placed in the imaging domain with 120 mm separation (See Figure 5A).The relative permittivity of the left and right cylinders was defined according to (52) with maximum values of 2 and 4, respectively.The imaging area was selected as a rectangular region of 200 Â 80 mm.
F I G U R E 3 Flow diagrams for the calculation of the derivative of the cost function -N sources and N receivers case.s is the number of forward simulation runs (FDTD runs), M is the number of grid cells, and N is the number of antennas where M ) N. The total number of forward simulation runs was provided at the bottom of the diagrams.

| Case 1: Single Source and N À 1 ð ÞReceivers
The resulting maps for the derivative of the cost function obtained with a cylindrical object (radius = 24 mm, ε r,max ¼ 2) using the FD, direct, and indirect methods in the frequency domain are shown in Figure 6.In addition, the normalized derivative values at y ¼ 0 mm cut (x ¼ [À100,100] mm) and x ¼ 0 mm cut (y ¼ [À100, 100] mm) are given in Figure 7. Similarly, the obtained results in the time domain are shown in Figure 8 and Figure 9.
The results show that the derivative values obtained from all methods are in good agreement, verifying that the adjoint method was formulated and implemented correctly in the both frequency and time domains.As a consequence of using only a single transmit antenna to illuminate the object, the information derived from the derivative maps is limited in terms of making inferences about the object.(Detailed results for a more realistic microwave imaging configuration are presented in the subsequent section.) In terms of computational load, a comparison of the utilized methods shows that the FD, direct, and indirect methods required M þ 1 ð Þ, N, and 2 forward simulation runs (FDTD runs), respectively.Since the number of unknowns (the grid cell dielectric constants) was greater than the number of antennas (M ) N), the computational burden of the FD method was enormous.Therefore, FD method is not practical for the iterative solution procedure.For a single source case, the direct method was advantageous in terms of the required forward simulation runs since it involved only two forward simulation runs.In this scenario, using all antennas sequentially as transmitters allowed information to be obtained about the shape and permittivity distributions of the objects through the derivative of the cost function.
For the first configuration, the relative permittivity map of the imaging area and the maps for the normalized derivative of the cost function are shown in Figure 4. Likewise, Figure 5 illustrates the relative permittivity map of the imaging area, the map for the normalized derivative of the cost function, the normalized derivative, and the relative permittivity along the center line for the second configuration.The derivative of the cost function has a smooth variation around the dielectric objects, and a peak around the object centers can be observed for both configurations.In addition, the objects with greater relative permittivity resulted in a higher gradient value (larger amount of change), as expected.
In terms of computational load, compared to the single source and N À 1 ð Þreceiver case, the number of total forward simulation runs was multiplied by the number of antennas N ð Þ for the FD and the direct methods, resulting and (2 Â N) forward simulation runs (FDTD runs), respectively (see Figure 3).On the other hand, the number of total simulations did not change when the indirect method was used.Due to the reciprocity of the antennas, the field distribution was the same (or scaled if the source amplitudes were different) when an antenna radiates the actual source or the adjoint source.Therefore, one forward solution for each antenna, resulting in N forward simulation runs, was adequate to obtain the derivative of the cost function.Consequently, the computational burden of the FD method is the greatest compared to the adjoint methods.Among the adjoint methods, the indirect method has higher computational efficiency than the direct method when the antennas were used as both transmitter and receiver.

| Computational Efficiency
The computational efficiency of the direct and indirect methods depends on the imaging configuration and the utilization of the antennas (or sensors).Note that the computational efficiency is evaluated in terms of the number of forward solution runs since it is the most time-consuming part of the iterative microwave tomography algorithms.For an imaging configuration, assume N s and N r correspond to the number of source and receiver antennas, respectively.To calculate the derivative of the cost function using the direct approach, 2 Â N s forward simulations are required.On the other hand, for the indirect approach, the number of the required forward simulation runs is N s þ N r .However, the transmitters and receivers can be used interchangeably in an imaging scenario where the reciprocity principle holds.Therefore, there is no need to recalculate the field when an antenna switches roles.As a result, the number of forward simulation runs for the indirect method reduces to where N s \ N r ð Þcorresponds to the number of antennas that can be used as both transmitter and receiver, i.e., transceivers.If we compare the direct and indirect methods, the indirect method is efficient when This expression can be simplified as 2][73] Accordingly, even if one is not interested in the sensitivity or Jacobian, the indirect method may still be chosen for calculating the derivative of the cost function due to the smaller number of forward solutions.
In microwave imaging, antennas can function as both transmitters and receivers, with N s ¼ N r ¼ N, and N s \ N r ð Þ¼N.As a result, the number of forward solutions required for the direct and indirect methods can be expressed as 2 Â N and N, respectively.It can be observed that the indirect method requires half the forward solutions compared to the direct method to obtain the functional gradient.Therefore, the indirect method provides advantages in certain scenarios even when sensitivity is not required, particularly in scenarios where computational efficiency and reduced forward solutions are crucial factors.

| Calculation of the update direction for gradient descent and Gauss-Newton algorithms
Comparing the obtained results for commonly used Gauss-Newton and gradient descent (GD)-based algorithms in microwave imaging would provide valuable perspectives.Note that our aim is not to compare the reconstruction performance of these methods since other reconstruction parameters, such as convergence, regularization, image quality, and signal-to-noise ratio analysis, need to be analyzed, and they are out of scope.Instead, we compared the utilization of the different adjoint methods for calculating the update direction, u, that both GD and GN methods require.The update equation for the unknown model parameter m (in our case, relative permittivity ε r ) can be expressed as where n is the iteration number, α is the step length, and u is the update direction.For gradient-descent algorithm, the update direction u can be written as whereas, for Gauss-Newton algorithm, the update direction u can be expressed as For gradient descent algorithms, the update direction, u, can be obtained using both the direct and indirect methods.The direct method provides the derivative of the cost function rF ð Þ directly, while the indirect method obtains the functional derivative through the Jacobian and residual J T e À Á .The selection between the direct and indirect methods depends on the imaging configuration and computational efficiency requirements.For example, in cases where forward solution computation costs are the most substantial, as in microwave imaging, and the antennas can be used as both transmitter and receiver, the indirect approach becomes advantageous since it requires a smaller number of forward solutions compared to the direct method.Therefore, calculating the Jacobian first, then the gradient of the cost function, that is, the indirect method, becomes more efficient.However, it should be noted that this method requires a large-scale and memory-inefficient Jacobian matrix.This problem can be solved by dividing the Jacobian matrix into sub-matrices. 71n the other hand, for Gauss-Newton inversion, the update direction, u, requires Jacobian calculation; therefore, a comparison can be made between the indirect and FD methods.As expected, the indirect method leads to better computational efficiency since the number of the grid cells in the imaging domain M ð Þ was greater than the number of antennas M ) N ð Þ .The calculation of the Jacobian matrix offers another noteworthy observation.The Jacobian matrix is calculated row-by-row, where each row is computed by multiplying the actual electric field generated by the transmitter antenna with the adjoint field caused by the receiver antenna, see (50).Since the adjoint field is equivalent to the forward solution of the receiver antenna (with a scaling factor), each row of the Jacobian matrix corresponds to the multiplication of the two forward solution vectors associated with the selected pair of transmit and receive antennas. 749][80][81] The low adoption of quasi-Newton approaches might be due to the slower convergence rate compared to Newton-like methods. 80,813][84][85] In quasi-Newton methods, the Hessian matrix, which includes the computationally expensive second order derivative terms, is not used directly, but an approximate Hessian matrix is calculated using the gradient information obtained from previous iterations. 78,81The gradient information (update direction) can be calculated using both the direct and indirect adjoint approaches, similar to gradient-descent methods.As previously stated, the most computationally efficient method depends on the addressed problem and the imaging configuration.

| The comparison of frequency domain and time domain adjoint method
As mentioned in the introduction, frequency domain algorithms are often used more than time domain algorithms in microwave imaging.However, time-domain reconstruction brings the advantage of improved resolution due to its broadband information.Although the same outcome can be achieved through multifrequency or frequency hopping reconstruction techniques, 86,87 they bring the challenge of selecting appropriate frequencies and determining the number of frequencies to be used.Depending on the frequency selection strategy, non-linearity-related inaccuracies may arise in the resulting reconstructions.Therefore, the time-domain reconstruction is considered more robust than frequency-domain algorithms. 43,53On the other hand, in time-domain algorithms, the initial frequency domain data needs to be calibrated and transformed into the time domain.This requirement arises because microwave tomography experiments are usually conducted using a vector network analyzer, and the broadband data is acquired in the frequency domain.
In terms of computation time, time domain reconstruction methods in microwave tomography literature have been observed to be slower than their frequency domain counterparts. 53For time-domain reconstructions, the FBTS method, which utilizes the direct adjoint method-where received residual signals are reversed in time and used as sources-is usually preferred.Therefore, the forward problem needs to be solved twice for each transmitter at each iteration.In addition, the FBTS method uses gradient-based algorithms to minimize the cost function iteratively.Instead of using the direct method in the time domain, the indirect method can provide comparable computational efficiency to the frequency domain.First, the indirect method requires half the forward solutions, and second, the indirect method allows the use of the Gauss-Newton inversion, which converges faster than gradient descent when the initial model is sufficiently close to a local or global minimum.Therefore, using the indirect method with Gauss-Newton inversion in the time domain may provide an additional advantage in terms of computation time compared to the frequency domain.
The simulation setup used in this study has been configured to provide a functional gradient map that can serve as a good initial guess in the frequency domain.Therefore, the results of frequency and time domain reconstructions are quite close to each other.When the results for Case-1 (one source, NÀ1 receivers) are compared between the frequency domain and time domain, the functional gradient maps appear smoother in the time domain due to the utilization of wideband information.
Lastly, we would like to emphasize that for both frequency domain and time domain algorithms, the map of the derivative of the cost function strongly depends on the imaging configuration (i.e., operation frequency, number, and location of the antennas, size of the computational domain, grid size, etc.) and object parameters (location, size, proximity, and relative permittivity).Therefore, the location of the objects may not be clearly visible as in the examples given here.Nevertheless, the resulting derivative maps provide an initial guess and intuition for inverse problems solved by iterative methods.

| Applications of the adjoint method in inverse problems
The determination of the update direction using the adjoint method is frequently utilized in inverse problems.Therefore, the concepts presented here are not limited to microwave imaging and apply to scenarios where physical models are characterized by a mathematical framework that involves partial differential equations along with boundary and initial conditions.We would like to briefly provide some examples of inverse problems that utilize the adjoint method for optimization.
In the field of biomedical imaging, the adjoint method is widely used in ultrasound computed tomography. 19,88In this technique, the tissue is imaged using ultrasound transducers, and the acquired signals are used to reconstruct acoustic properties.The ultrasound transducers are reciprocal, that is, the sources (transmitters) and receivers (detectors) can be used interchangeably.Depending on the optimization algorithm, both the direct and indirect methods can be utilized.Usually, the direct adjoint method is used for the full waveform inversion.Using the indirect adjoint method for the optimization may provide a computational advantage in terms of the number of forward simulations.In addition, the indirect method allows the calculation of the Jacobian matrix, which can be used in fast converging algorithms.However, it should be noted that this method requires the development of a memory-efficient implementation of the Jacobian matrix.
In addition, the adjoint method finds application in many other disciplines, including seismic 71,72,89 and geophysics. 59,66,90In the seismic waveform inversion, seismic waves are generated through a controlled source and detected by seismic sensors placed on the earth's surface.Typically, the number of receivers exceeds the number of sources. 91,92Similarly, in ground-penetrating radar tomography, permittivity and conductivity distributions are usually investigated in the shallow subsurface using a few tens to a few hundred receivers. 66,93In these methods, the direct method is commonly utilized due to the greater number of receivers compared to sources.Typically, the cost function is minimized through the direct method (e.g., gradient-based methods) without calculating the Jacobian matrix.However, the computation of the Jacobian matrix can be preferred, as it provides additional information and allows for the performance of cumulative sensitivity and formal model resolution analyses. 66

| CONCLUSION
In this study, we focused on investigating the different applications of the adjoint method in microwave imaging.We have provided clear and concise analytical expressions for the adjoint formulations, filling a gap in existing literature where such formulations are predominantly presented based on specific applied problems.We classified the update direction calculations using adjoint method into two groups, the direct and the indirect methods, and presented the derivations independent of the forward solution method.The resulting analytical expressions were verified through numerical simulations, and the computational loads were measured in terms of the number of forward simulation runs.The results have demonstrated the differences between the adjoint methods and have shown the computational efficiency that can be achieved through the proper strategy.Moreover, the study has provided observations into how the method is utilized, thereby enhancing its practical utility.

F I G U R E 4 5
Simulation results of the first configuration for a cylindrical object (radius = 24 mm, ε r,max ¼ 2) inside the imaging area of 200 Â 200 mm.(A) Relative permittivity, (B) Normalized derivative of the cost function in the frequency domain, (C) Normalized derivative of the cost function in the time domain.The results are obtained for the N sources and N receivers case using the indirect method.Simulation results for the second configuration.(A) Relative permittivity map inside the imaging area of 200 Â 80 mm.The cylinders had diameters of 15 and 21 mm, and their maximum relative permittivity values were 2 and 4, respectively.(B) Relative permittivity values along y ¼ 0 mm cut (x = [À100, 100] mm).(C) Map for the normalized derivative of the cost function in the frequency domain.(D) Normalized derivative values along y ¼ 0 mm cut (x = [À100, 100] mm) in the frequency domain.(E) Map for the normalized derivative of the cost function in the time domain.(F) Normalized derivative values along y ¼ 0 mm cut (x = [À100, 100] mm) in the time domain.The results are obtained for the N sources/receivers case using the indirect method.

F I G U R E 6 8
Maps for the normalized derivative of the cost function obtained using the finite-difference method, direct method, and indirect method inside the imaging area of 200 Â 200 mm for a cylindrical object (radius = 24 mm, ε r,max ¼ 2).The results are obtained in the frequency domain for the single source and (NÀ1) receiver case.The color bar shows the normalized value of the absolute of the derivatives.(A) F I G U R E 7 Normalized values for the derivative of the cost function at (A) y ¼ 0 mm cut (x ¼ [À100, 100] mm).(B) x ¼ 0 mm cut (x ¼ [À100, 100] mm) for a cylindrical object (radius = 24 mm, ε r,max ¼ 2) using the finite-difference, direct, and indirect methods.The results are obtained in the frequency domain for the single source and (NÀ1) receiver case.Maps for the normalized derivative of the cost function obtained using the finite-difference method, direct method, and indirect method inside the imaging area of 200 Â 200 mm for a cylindrical object (radius = 24 mm, ε r,max ¼ 2).The results are obtained in the time domain for the single source and (NÀ1) receiver case.The color bar shows the normalized value of the absolute of the derivatives.

9
Normalized values for the derivative of the cost function at (A) y ¼ 0 mm cut (x ¼ [À100, 100] mm).(B) x ¼ 0 mm cut (x ¼ [À100, 100] mm) for a cylindrical object (radius = 24 mm, ε r,max ¼ 2) using the finite-difference, direct, and indirect methods.The results are obtained in the time domain for the single source and (NÀ1) receiver case.