Constrained Pseudoinverses for the Electromagnetic Inverse Source Problem

Inverse source strategies have proven to be quite relevant for several applications in advanced electromagnetics. These schemes are based on the solution of ill-posed problems in which current or near-field distributions are reconstructed from far-field (or from less informative field) information. Standard strategies, that can include physical constraints such as Love conditions, often rely on standard pseudoinverse definitions and yield solutions that are, at times, far from the physical ones. This work proposes a different approach focusing on defining and analyzing a new family of pseudoinverses that takes advantage of small-in-dimension subspaces containing a priori information. The new solutions returned by the new pseudoinverses will be a suitable average between a solution living entirely in the vector space containing the a priori information and a solution obtained via norm-minimizing approaches. The contribution presents both theoretical analyses and numerical experiments showing the practical effectiveness of the novel mathematical tool.

that describes the radiation [5], [6].This null-space, rooted in the physical and mathematical nature of the problem, is present independently of the choice of discretization; in fact, inverse approaches based either on the multipole expansion of the fields [7], [8] or on a nonorthogonal basis [9], [10], [11] all have to select a specific solution to the ill-posed problem.Several approaches have been adopted in the literature for selecting the appropriate solution.They include the choice of a specific solution norm (e.g., the popular Moore-Penrose (MP) pseudoinverse approaches, see [12] and references therein), criteria embedded in an iterative procedure [13], [14], or approaches leveraging on explicit conditions such as the Love conditions [15], [16], [17] and related [18], [19].
Moreover, the dimension of the (numerical) nullspace of the radiation matrix is further increased by the exponentially decaying contribution of the evanescent modes on the total field, when observed at far-field distances from the source.Because only a finite number of electromagnetic modes actually propagate [20], even if the presence of the null space and the consequent rank deficient matrix can be tackled via pseudoinversion, the loss of information associated with the evanescent modes is irreversible, degrades the accuracy of the solution [21], [22], and can only be handled through the injection of a priori information into the mathematical problem.This strategy has been used widely in the past to reduce the required number of samples below the Nyquist limit or, equivalently, to improve field reconstruction given a fixed number of samples (see [23], [24] and references therein).In [23], the number of samples and their position are found iteratively through the use of an orthonormal basis constructed by the simulation of the antenna under test with different design parameters.A similar construction of an overcomplete basis is presented in [24] and a given set of real samples is used to guide an iterative algorithm in the selection of the optimal subbasis.Historically, the task of finding a set of bases on which to project and approximate a given datum has also been accomplished by the compressive sampling community [25] covering a wide range of applications and integrating signal processing strategies [26].
The above-mentioned approaches are effective when the amount of information to inject is abundant enough to delineate an optimal basis to describe the observation space.Our approach will instead deal with a different scenario in which the information available is not sufficient to completely eliminate the ill-posedness of the problem, a pseudoinversion is still required, and the question of how the information at hand can be used to choose a better pseudoinversion definition with respect to a standard, a priori-information agnostic choice remains.
The purpose of this work is to propose, define, and analyze a new family of pseudoinverses that will be able to take advantage of a small-in-dimension subspace containing a priori information.The new operators can be seen as a properly chosen average between a solution living entirely in the vector space containing the a priori information and a solution obtained via a norm minimizing approach.This average is built so that the operator remains a pseudoinverse, although not of MP type.Theoretical developments are alternated with numerical studies to show the effectiveness and practical relevance of our approach when applied to real case scenarios.
The article is organized as follows.The notation and the relevant background on the inverse source problem are set in Section II.In Section III, we present a novel analysis of the radiation operator via vector spherical harmonics (VSH) basis in which the effects of the nonradiating sources and the evanescent fields on the conditioning is highlighted.The new constrained pseudoinverse is presented in Section IV while the numerical results are delineated in Section V. Finally, Section VI presents our conclusions and venues for future investigations.Very preliminary results of this article have been presented as a conference contribution [27].

II. BACKGROUND AND NOTATION
Let Ω − be a closed domain in R 3 with boundary Γ and complement Ω + = R 3 \Ω − .Consider a source located in Ω − which radiates in R 3 the electric and magnetic fields E and H.According to the equivalence theorem, the electromagnetic problem in Ω − can be changed (medium, fields, and sources [28]) by placing on Γ equivalent magnetic and electric current densities M, J that radiate E, H in Ω + and arbitrary new fields E ′ , H ′ in Ω − .These currents must satisfy the Maxwellian conditions at the interface where n is the surface-normal unit vector in r directed from Ω − to Ω + , E + and H + are the original electric and magnetic fields evaluated in r ∈ Ω + , and E ′− and H ′− are the new fields in r ∈ Ω − .In (1) and ( 2), the external and internal fields are evaluated in the limit r → r ∈ Γ taken from Ω + and Ω − , respectively.The e −iωt time dependence of the physical quantities is assumed throughout the article.The inverse source problem finds M, J given field observations on a surface Γ m in Ω + .In this context, the electric field integral operator (EFIO) and the magnetic field integral operator (MFIO) applied to a vector function, e.g., J, are defined respectively as where r ∈ Γ m is the observation point, k is the wavenumber, and g(∥r − r ′ ∥) = (e ik∥r−r ′ ∥ /(4π∥r − r ′ ∥)).If r ∈ Γ , then T r and K r are denoted by T and K, respectively.These integral operators are involved in a variety of surface and radiation electromagnetic problems.Among them, the electric field integral equation (EFIE) relates the electric current J on the surface Γ of a perfect electric conductor (PEC) to an incident electric field E i .The electromagnetic field generated by a given source in Ω − can be evaluated in r ∈ Γ m , by means of equivalent currents as where η = √ µ/ϵ is the wave impedance of the considered homogeneous medium, chosen to be the same in Ω − and Ω + .Finally, n × E + and n × η H + are obtained as In the following, we will consider three different discretizations of these continuous problems that are described in the three following paragraphs.
As first discretization scheme, the equivalent currents J and M are approximated with div-conforming Rao-Wilton-Glisson (RWG) basis functions { f n } n [29] defined on a triangular mesh of Γ , so that M ≈ It is then possible to describe the fields they radiate in a homogeneous medium through where R is the discretization of the continuous operator R in ( 9) and the elements of C r,δ , S r,δ are obtained, in a pointmatching fashion, as being the number of observation points.The linear system (11) samples the electric and the magnetic field in r with e δ and h δ containing the Cartesian coordinates of the vectors ⟨δ(r − r m ), E + ⟩ and ⟨δ(r − r m ), H + ⟩, respectively.An alternative to (11) can be obtained by using curl-conforming rotated RWG functions { n × f n } n in place of the Dirac delta functions as test functions.The discretized boundary problem (8) then reads with Finally, the tangential fields in (10) are approximated with where Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Finally, a third discretization scheme uses div-conforming Buffa-Christiansen (BC) dual basis functions {g n } n , defined on the dual mesh obtained through barycentric refinement of the original mesh [30] in place of the RWG basis function.In particular, the BC-discretized radiation operators S r,δ , C r,δ , are defined with [S r,δ ] mn = ⟨δ(r − r m ), S r g n ⟩ and [C r,δ ] mn = ⟨δ(r − r m ), C r g n ⟩.The RWG-and BC-based discretization schemes are linked through the well-conditioned change of basis matrix G mix , with In the following, the different discretization strategies delineated above will be combined to form mappings between equivalent surfaces and measurement surfaces that take into account the mapping properties of the operators and the constraints in chaining the different bases.

III. SPECTRAL ANALYSIS
We present here a VSH analysis that highlights the sources of ill-posedness of (10) that describes the radiation.Consider the VSH basis [31] with where P lm are the associated Legendre polynomials and the integers l and m are such that 0 ≤ l and −l ≤ m ≤ l.Applying (3) and (4) to a vector function expanded with ( 14) and (15) yields T r X lm = − a r J l (ka)H (1)  l (kr )U lm (18) where a is the radius of Γ , r = ∥r∥, and J l , H (1)   l are the Riccati-Bessel and the first-kind Riccati-Hankel functions, respectively.Expressions (18) to (21) are used in (10) and R t is discretized with the VSH basis, truncated to l ≤ L, as where the subscript ϕ is used to distinguish the VSH and RWG discretizations.Then, In (23), the superscripts X X , XU , U X , and UU denote the test and source functions couples used to discretize the associated block; for example, where the subscripts on X lm , U lm span l ≤ L, −l ≤ m ≤ l, and ⟨a, b⟩ t,ϕ R t,ϕ are diagonal, their singular values cannot be obtained by direct inspection.Instead, a singular value decomposition (SVD) is computed and the singular values of R t,ϕ are shown in Fig. 1 The second half of the singular values (i.e., i > 1350) belongs to the nullspace of the operator, while the singular values of the first half decay more rapidly as the observation distance increases [20], [21].These singular values are associated with the evanescent fields and extend the numerical nullspace of the operator.For this reason, in practical applications, the condition number has to be limited by truncation of the singular values lower than the measurement noise floor (NF), causing an overall loss of near-field information [22].
In this work, an operator tailored for the inverse source problem uses a pseudoinverse which acts on a small subspace of a priori solutions and is able, if built properly, to recover more accurately the weakly radiating part of the solution, starting from far-field samples.

IV. NEW CONSTRAINED PSEUDOINVERSE A. Insights and Derivation of the Proposed Pseudoinverse
First, we recall that a generic pseudoinverse A † of A ∈ C M×N satisfies any of the properties and the well-known MP pseudoinverse A + of A is the unique operator which satisfies i), ii), iii), and iv).If one is only interested in finding one solution x 0 (potentially among many) of the linear system Ax = b, one can choose x 0 = A (1) b, where A (1) is a generic pseudoinverse satisfying i) since [32] Ax 0 = AA (1) b = AA (1) When A (1) varies in the set of all pseudoinverses satisfying i), denoted by A{1}, x 0 = A (1) b can be shown to deliver all possible solutions of the linear system.Since in many practical cases, there are solutions of rank-deficient matrices that have Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
little to no physical meaning, the condition i) alone does not ensure that the pseudoinverse and the corresponding solution have any link with the physical problem.
If we start from the opposite standpoint, we could build a pseudoinverse from existing (physical) solutions of the problem x i , i = 1, . . ., P and the corresponding right-handsides b i = Ax i .By defining the matrices X ∈ C N ×P as X j,i = [x i ] j and B ∈ C M×P as B j,i = [b i ] j , we could now choose as pseudoinverse of A the matrix XB † , which is a pseudoinverse of type ii) if the pseudoinverse B † of B is of type ii).This pseudoinverse provides a vector x 0 = XB † b which is always a linear combination of solutions of the system for selected RHSs and, if those are properly chosen, such a linear combination will exhibit physical properties.However, XB † / ∈ A{1} unless BB † = I, i.e., this choice of pseudoinverse does not provide in general an exact solution satisfying (25).In this work, we propose a pseudoinverse that would somehow be in between the above described approaches.In particular, the pseudoinverse we propose to adopt reads which is combining the previously defined approaches.Differently from XB † , A ‡ provides solutions which satisfy ( 25) since A ‡ ∈ A{1}: with B = AX.If B † B = I, the proposed pseudoinverse yields an exact reconstruction if the RHS is entirely described by a linear combination of the columns of B, i.e., there exists v such that b = Bv .In this case, we have where Xv is the exact solution associated with b = Bv = AXv .

B. On a Possible Strategy in Selecting Parameters for the Proposed Pseudoinverse
In the context of the inverse source problem-corresponding to A = R-the proposed pseudoinverse (26) may enhance the solution of (11) with near-field information encoded in X.The constrained pseudoinverse ( 26) is flexible and can be used in different inversion contexts and applications.In particular, different inversion scenarios may lead to selecting different pseudoinverse types for B † and R (1) [the latter can be of other types additionally to i)] and to different choices for the coupled matrices X and B = RX.In the following, we will consider the MP pseudoinverses B + and R + as particular choices for B † and R (1) , respectively.
The columns of X can be filled with solutions found either in real or simulated environments.In the following, we will opt for the latter.The choice of the vectors x i is heavily dependent on the scenario.In general, however, two distinct sources that belong to the same design space up to a given perturbation of one or more of the design parameters (position and shape) will generate the same field up to a given precision.Thus, the space of a priori solutions x i can be generated by solving several radiation problems in which these unknown parameters differ.In the following, we will provide initial application examples that will help fixing ideas.

V. NUMERICAL RESULTS
The proposed pseudoinverse is first numerically tested in a customary inverse source problem: the electric and magnetic fields generated by three electric dipoles, E + and H + , are sampled on a finite number of points and reconstructed at arbitrary positions.The dipoles, oscillating at 5 GHz, are placed in a sphere of radius a = 4 cm (≈0.67λ) and surface Γ discretized with triangular elements of average edge length h = λ/10, where λ is the free-space wavelength.The reference magnetic and electric equivalent currents on Γ are approximated with RWG functions whose coefficients m and j are The observation points on which field samples are evaluated with (11) are located on a spherical surface Γ m , concentric to Γ , and are Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Fig. 3. Reconstruction error evaluated on spherical surfaces concentric to Γ starting from samples generated 1λ away from Γ .The MP and the proposed pseudoinverse are compared, and the latter is tested when varying the number of a priori sources (with ∆α) and their distance from the original source (with d).A white NF of −50 dB is applied to the far-field observations.distributed on uniform polar and azimuthal angle grids.The field observations are perturbed with white Gaussian noise, generating a NF of −50 dB; Γ m is set 1λ away from Γ along its outward radial direction, and thus in a far-field region with respect to the dipole source.The injected white noise limits the reconstruction accuracy throughout all of the spatial band of radiation, in particular of the weakly radiating modes [22].As benchmark of (26), A (1) -here chosen to be A (1)  = R + as a particular case-is regularized using a truncated-SVD.The truncation index is chosen to minimize the reconstruction error at Γ m for the specific NF and right-hand-side under consideration.Then, a set of a priori solutions forming the columns of X is obtained by evaluating the equivalent currents coefficients of additional electric dipoles displaced inside Γ .These additional dipoles [the white dots in Fig. 2(b)] are displaced on spheres centered on the real sources [red diamonds in Fig. 2(b)]; the electric moment of all dipoles is identical.This test was designed to account for uncertainty in the knowledge of the position of the electromagnetic source.The goal is then to better reconstruct the fields everywhere with such a priori knowledge.The effectiveness of the choices of X, B in ( 26) is investigated when varying the angular resolution ∆α, which controls the distance between the constraining sources (and therefore their number).Once the solution of ( 26) is obtained, the fields are radiated on spherical surfaces concentric to Γ , and the relative error ϵ is evaluated on each sphere as where v and v ref are the reconstructed and reference vectors, respectively.The reconstruction errors of the new pseudoinverse (26) and of the MP pseudoinverse are shown in Fig. 3, the latter regularized through the same optimum criterion of R + in (26) while keeping the same observation points; as expected, reducing ∆α lowers the reconstruction error, since the dimension of the column space of X is increased in a meaningful way (auxiliary sources clustering around the source dipoles).We also note that at a distance of 1λ (on Γ m ), where the reconstruction points coincide with the observation points, the fields are not perfectly reconstructed and differ slightly between the different reconstruction methods.These effects can in part be explained by numerical error propagation through ill-conditioned matrices and by the usage of different pseudoinverses.Finally, and to better appreciate the displacement of the a priori sources and the near-field reconstruction improvement, in Fig. 2, the error ( 32) is evaluated on a square box with edge length 2.2a, centered at the origin: the reconstruction error for the proposed method is lower than that obtained via the MP pseudoinverse, everywhere on the sides of the box and, in particular, in the near-field regions around their centers.The constrained pseudoinverse is now tested on a more complex scenario in which a deformed and scaled-down PEC shuttle model is illuminated by a plane wave propagating in the direction r = [1, 0, 1] and oscillating at 5 GHz.The reference scattered field is sampled in the far-field region, and ( 26) is used to reconstruct more accurate scattered nearfields.The test could be representative of the case where the scatterer geometry is only approximately known and a priori assumptions can be made on the structure rather than on its displacement.First, the EFIE ( 12) is solved on Γ , a triangular mesh of the shuttle illustrated in Fig. 4, and the resulting electric current is denoted by j i .Then, the equivalent magnetic and electric currents are obtained on an equivalent surface Γ wrapping the scatterer following the appropriate discretization scheme as m = −G −1 mix e, j = ηG −1 mix h, with e = T r j i , h = −K r j i .These currents are the reference solution used to scatter the reference electric field in Ω + , being E + = −C r,δ m+S r,δ j .The constraining space of the proposed pseudoinverse makes use of two differently deformed shuttle meshes, denoted by Γ 1 and Γ 2 in Fig. 4. To obtain the matrix X of a priori solutions, the additional scatterers are illuminated with the same plane-wave applied to the reference Γ , the corresponding EFIE are solved, and the equivalent currents are found on Γ .Thus, X and B = RX are composed of two columns only.Finally, the reconstruction is made on concentric spheres, the smallest just big enough to enclose Γ , starting from far-field measurements Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.(1λ from a sphere of 8 cm radius containing Γ ) corrupted with NF of −60 dB.The reconstruction error is shown in Fig. 4, where the MP pseudoinverse is compared with the constrained pseudoinverse through ϵ(E + ), and there is clearly an improvement at near-field distances.Finally, in Fig. 5, ϵ(E + ) is evaluated on a cut of a surface wrapping Γ at distance ≈0.1λ; the constrained pseudoinverse yields a lower reconstruction error everywhere on the surface under consideration.

VI. CONCLUSION
This work has presented a new constrained pseudoinverse which exploits a priori information in its definition to guide the selection of a specific solution of the ill-posed inverse source problem.By its very definition, the proposed pseudoinverse should be used when a priori information is available about the source under investigation.Once this preliminary condition has been assessed, the a priori vectors can be accommodated according to a variety of strategies.We proposed and analyzed two lines of investigation: the first uses a priori information in the form of spatial knowledge of the source; the second uses a priori geometrical assumptions.In both scenarios, the proposed constrained pseudoinverse performed better than the MP pseudoinverse in reconstructing the near-fields, at the cost of generating the a priori solutions.

Fig. 1 .
Fig. 1.Singular values of the radiation operator R t,ϕ discretized with the VSH basis, evaluated at different distances; the σ i are ordered with the SVD index.

Fig. 2 .
Fig. 2. Reconstruction error on a square box of side length 2.2λ centered in the origin.(a) MP reconstruction.(b) Proposed pseudoinverse ∆α = π/2.The real sources are represented with red diamonds, the constraining dipoles with white dots; the background shows ϵ(E + ).

Fig. 4 .
Fig. 4. Reconstruction error evaluated on concentric spherical surfaces: the MP and the proposed pseudoinverse are compared, the latter uses Γ 1 andΓ 2 to generate a priori vectors while the mesh of the scatterer is Γ (the colors of the geometries only highlight the deformations).On the x-axis, the distance 0λ coincides with the sphere of radius 8 cm enclosing Γ ; the distance 1λ denotes the observation surface on which the electric field is sampled and corrupted with a Gaussian white NF of −60 dB.

Fig. 5 .
Fig. 5. Near-field reconstruction error.(a) MP reconstruction.(b) Proposed pseudoinverse given the a priori meshes of Fig. 4. The scatterer surface Γ is in red, the equivalent surface Γ in black and ϵ(E + ) is evaluated in the near-field of Γ (distance ≈ 0.1λ).