Efficient Computation of the Magnetic Polarizability Tensor Spectral Signature using POD

Our interest lies in the identification of hidden conducting permeable objects from measurements of the perturbed magnetic field in metal detection taken over range of low frequencies. The magnetic polarizability tensor (MPT) provides a characterisation of a conducting permeable object using a small number of coefficients, has explicit formula for their calculation and a well understood frequency behaviour, which we call its spectral signature. However, to compute such signatures, and build a library of them for object classification, requires repeated solution of a direct (full order) problem, which is typically accomplished using a finite element discretisation. To overcome this issue, we propose an efficient reduced order model (ROM) using a proper orthogonal decomposition (POD) for the rapid computation of MPT spectral signatures. Our ROM benefits from output certificates, which give bounds on the accuracy of the predicted outputs with respect to the full order model solutions. To further increase the efficiency of the computation of the MPT spectral signature, we provide scaling results, which enable an immediate calculation of the signature under changes in the object size or conductivity. We illustrate our approach by application to a range of homogenous and inhomogeneous conducting permeable objects.


Introduction
There is considerable interest in using the magnetic polarizability tensor (MPT) characterisation of conducting permeable objects to classify and identify hidden targets in metal detection. The MPT is a complex symmetric rank 2 tensor, which has 6 independent coefficients, although the number of independent coefficients for objects with rotational or reflectional symmetries is smaller [14]. Its coefficients are a function of the exciting frequency, the object's size, its shape as well as its conductivity and permeability. Explicit formulae for computing the tensor coefficients have been derived [4,14,16,17] and validated against exact solutions and measurements [15,16]. Also, the way in which the tensor coefficients vary with the exciting frequency is theoretically well understood [17] offering improved object classification.
The frequency (or spectral) behaviour of the MPT, henceforth called its spectral signature, has been exploited in a range of different classification algorithms including simple library classification for homogeneous [5] and inhomogeneous objects [18], a k nearest neighbours (KNN) classification algorithm [21] and machine learning approaches [33]. The MPT classification of objects has already been applied to a range of different applications including airport security screening [22,21], waste sorting [12] and anti-personal landmine detection [27]. The aforementioned supervised classification techniques rely on a library of MPT spectral signatures to learn how to classify the objects. The purpose of this paper is to describe an efficient computational tool for computing this library.
One approach to obtaining a library of spectral signatures is to use a metal detector or dedicated measurement device to obtain MPT coefficients of different objects [38,37,26], however, to do so, over a range of frequencies for a large number of objects, is time consuming and will result in unavoidable measurement errors and noise. Therefore, there is considerable interest in their automated computation. By post processing finite element method (FEM) solutions to eddy current problems using commercial packages (e.g with ANSYS as in [26]) MPT coefficients can be obtained, however, improved accuracy, and a better understanding, can be gained by using the available explicit expressions for MPT coefficients, which rely on computing finite element (FE) approximations to a transmission problem [14,16,17]. Nevertheless, to produce an accurate MPT spectral signature, the process must be repeated for a large number of excitation frequencies leading to potentially expensive computations for fine discretisations (with small mesh spacing and high order elements). The present paper addresses this issue by proposing a reduced order model, in the form of a (projected) proper orthogonal decomposition (POD) scheme, that relies on full order model solutions computed using the established open source FE package, NGSolve, and the recently derived alternative explicit expressions formulae for the MPT coefficients [17]. The use of NGSolve [29,28] ensures that the solutions to underlying (eddy current type) transmission problems are accurately computed using high order Hpcurlq conforming (high order edge element) discretisations (see [19,30,36] and references therein) and the POD technique ensures their rapid computation over sweeps of frequency.
Reduced order models (ROMs) based on POD have been successfully applied to efficiently generate solutions for new problem parameters using a small number full order model snapshots in a range of engineering applications including mechanics [23,25], thermal problems [35,6], fluid flow [20,24] as well as electromagnetic problems with application to integrated circuits [13]. However, they have not been applied to the computation of MPT spectral signatures. A review of current POD techniques is provided in [11,9].
The main novelty of the work is the application of a POD approach for the efficient and accurate computation of the MPT spectral signature and the derivation of output certificates that ensure accuracy of the reduced order predictions. This ROM approach is motivated by the previous success of POD approaches and the theoretical study [17], which shows the spectral behaviour of the MPT is characterised by a small number of functions and, hence, has a sparse representation. The practical computation requires only computing full order model solution snapshots at a small number of frequencies and the evaluation of the MPT spectral signature follows from solving a series of extremely small linear systems. A second novelty is the presentation of simple scaling results, which enable the MPT spectral signature to easily be computed from an existing set of coefficients under the scaling of an object's conductivity or object size.
The paper is organised as follows: In Section 2, the eddy current model, which applies in metal detection, and the asymptotic expansion of the perturbed magnetic field in the presence of a conducting permeable object, which leads to the explicit expression of the MPT, is briefly reviewed. Then, in Section 3, the FE model used for full order model problem is described. Section 4 presents the POD reduced order model scheme. This is followed, in Section 5, by the derivation of results that describe the scaling of the MPT under parameter changes. Sections 6 and 7 present numerical examples of the POD scheme for computing the frequency behaviour of the MPT and examples of the scaling of the MPT under parameter changes, respectively.

The eddy-current model and asymptotic expansion
We briefly discuss the eddy-current model along with stating the asymptotic expansion that forms the basis of the magnetic polarizability description of conducting objects in metal detection.

Eddy-current model
The eddy current model is a low frequency approximation of the Maxwell system that neglects the displacement currents, which is valid when the frequency is small and the conductivity of the body is high. A rigorous justification of the model involves the topology of the conducting body [3]. The eddy current model is described by the system where E α and H α are the electric and magnetic interaction fields, respectively, J 0 is an external current source, i :" ?´1 , ω is the angular frequency, µ is the magnetic permeability and σ is the electric conductivity. We will use the eddy current model for describing the forward and inverse problems associated with metal detection.

Forward problem
In the forward (or direct) problem, the position and materials of the conducting body B α are known. The object has a high conductivity, σ " σ˚, and a permeability, µ " µ˚. For the purpose of this study, the conducting body is assumed to be buried in soil, which is assumed to be of a much lower conductivity so that σ « 0 and have a permeability µ " µ 0 :" 4πˆ10´7H/m. A background field is generated by a solenodial current source J 0 with support in the air above the soil, which also has σ " 0 and µ " µ 0 . The region around the object is B c α :" R 3 zB α as shown in Figure 1. Note that a similar model also applies in the situation of identifying hidden targets in security screening [22,21] and waste sorting [12] amongst others. Figure 1: A diagram showing a hidden conducting object Bα, buried in soil, with a current source located in the air above.
The forward model is described by the system (1), which holds in R 3 , with µpxq " and the regions B α and B c α are coupled by the transmission conditions which hold on Γ α :" BB α . In the above, rus Γα :" u|`´u|´denotes the jump, the`refers to just outside of B α and the´to just inside and n denotes a unit outward normal to Γ α . The electric interaction field is non-physical in B c α and, to ensure uniqueness of this field, the condition ∇¨E α " 0 is imposed in this region. Furthermore, we also require that E α " Op1{|x|q and H α " Op1{|x|q as |x| Ñ 8, denoting that the fields go to zero at least as fast as 1{|x|, although, in practice, this rate can be faster.

Inverse problem
In metal detection, the inverse problem is to determine the location, shape and material properties (σå nd µ˚) of the conducting object B α from measurements of pH α´H0 qpxq taken at a range of locations x in the air. As described in the introduction, there are considerable advantages in using spectral data, i.e. additionally measuring pH α´H0 qpxq over a range of frequencies ω, within the limit of the eddy current model. Here, H 0 denotes the background magnetic field and E 0 and H 0 are the solutions of (1) with σ " 0 and µ " µ 0 in R 3 . Similar to above, we also require the decay conditions E 0 " Op1{|x|q and H 0 " Op1{|x|q as |x| Ñ 8. Note that practical metal detectors measure a voltage perturbation, which corresponds to ş S n¨pH α´H0 qpxqdx over an appropriate surface S [16]. For very small coils, this voltage perturbation is approximated by m¨pH α´H0 qpxq where m is the magnetic dipole moment of the coil [16].
A traditional approach to the solution of this inverse problem involves creating a discrete set of voxels, each with unknown σ and µ, and posing the solution to the inverse problem as an optimisation process in which σ and µ are found through minimisation of an appropriate functional e.g. [32]. From the resulting images of σ and µ one then attempts to infer the shape and position of the object. However, this problem is highly ill-posed [8] and presents considerable challenges mathematically and computationally in the case of limited noisy measurement data.
Instead, we seek an approximation of the perturbation pH α´H0 qpxq at some point x exterior to B α , which allows objects to be characterised by a small number of coefficients in a MPT that are easily obtained from the measurements of pH α´H0 qpxq once the object position is known, which can be found from a MUSIC algorithm for example [4]. The object identification then reduces to a classification problem, as discussed in the introduction.

The asymptotic expansion and MPT description
Following [4,14] we define B α :" αB`z where B is a unit size object, α is the object size and z is the object's translation from the origin as shown in Figure 2. Then, using the asymptotic formula obtained by Ammari, Chen, Chen, Garnier and Volkov [4], Ledger and Lionheart [14] have derived the simplified form pH α´H0 qpxq i " pD 2 x Gpx, zqq ij pMq jk pH 0 pzqq k`O pα 4 q, which holds as α Ñ 0 and makes the MPT explicit. The relationship between the leading order term in the above to the dipole expansion of pH α´H0 qpxq is discussed in [16]. In the above, Gpx, zq :" 1{4π|x´z| is the free space Laplace Green's function, D 2 x G denotes the Hessian of G and Einstein summation convention of the indices is implied. In addition, M " pMq jk e j be k , where e i denotes the ith orthonormal unit vector, is the symmetric rank 2 MPT, which describes the shape and material properties of the object B α and is frequency dependent, but is independent of the object's position z. We will sometimes write MrαB, ωs to emphasise this. The above formulation, and the definition of M below, are presented for the case of a single homogenous object B, the extension to multiple inhomogeneous objects can be found in [18,17].
Using the derivation in [17], we state the explicit formulae for the computation of the coefficients of M, which are particularly well suited to a FEM discretisation. The earlier explicit expressions in [14,15,16] are equivalent for exact fields. We use the splitting pMq ij :" pN 0 q ij`p Rq ij`i pIq ij obtained in [17] with pRrαB, ωsq ij :"´α where N 0 rαBs, RrαB, ωs and IrαB, ωs are real symmetric rank 2 tensors, which each have real eigenvalues.
In the above,μ r pξq :" " µ r :" µ˚{µ 0 ξ P B 1 ξ P B c , and ν :" α 2 ωµ 0 σ˚, δ ij is the Kronecker delta and the overbar denotes the complex conjugate. The computation of (5) rely on the solution of the transmission problems [17] ∇ˆμ´1 r ∇ˆθ where Γ :" BB and Note also that we choose to introduceθ i´e iˆξ , which can be shown to satisfy the same transmission problem as (6) except with a non-zero jump condition for rnˆμ´1 r ∇ˆθ p0q i s Γ and the decay conditionθ

Full order model
To approximate the solutions to the transmission problems (6) and (7) we truncate the unbounded domain B c at a finite distance from the object B and create a bounded domain Ω containing B. On BΩ, we approximate the decay conditions (6e) and (7f) by nˆθ p0q i " nˆpθ p0q i´e iˆξ q " 0 and nˆθ p1q i " 0, respectively. On this finite domain, we approximate the associated weak variational statements to these problems using FEM with a Hpcurlq conforming discretisation with mesh spacing h and order elements p where Hpcurlq :" u : u P pL 2 pΩqq 3 , ∇ˆu P pL 2 pΩqq 3 ( , and L 2 pΩq denotes the standard space of square integrable functions. In Section 3.1 we provide their weak formulations and provide their discretisation in Section 3.2. Henceforth, we call this discrete approximation the full order model.

Weak formulation of the problem
Following the approach advocated in [19] for magnetostatic and eddy current problems, we add a regularisation term ε ş Ωθ p0q i¨ψ dξ, where ε is a small regularisation parameter, to the weak variational statement of (6), written in terms ofθ p0q i , in order to circumvent the Coulomb gauge ∇¨θ p0q i " 0. For details of the small error induced by this approximation see [19,36]. Then, by choosing an appropriate set of Hpcurlq conforming finite element functions in W phpq Ă Hpcurlq, we obtain the following discrete regularised weak form for (6) : Find real solutionsθ In a similar manner, the discrete weak variational statement of (7) is: Find complex solutions θ for all ψ phpq P Y ε X W phpq where the overbar denotes the complex conjugate.
For what follows it is beneficial to restate (10) in the following form: Find θ ş Ω u¨vdξ denotes the L 2 inner product over Ω and ω indicates the list of the problem parameters tω, σ˚, µ r , αu that one might wish to vary. Note that r`¨;¨,¨˘is a function of µ r as θ p0,hpq i depends on µ r .

Finite element discretisation
For the implementation of the full order model, we use NGSolve [1,29,28,36] along with the hierarchic set of Hpcurlq conforming basis functions proposed by Schöberl and Zaglmayr [30], which are available in this software. In the following, for simplicity, we focus on the treatment of θ p1,hpq i and drop the index i as each direction can be computed in a similar way (as canθ p0,hpq i ). We denote these basis functions with N pkq pξq P W phpq leading to the expression of the solution function along with the weighting functions θ p1,hpq pξ, ωq :" ψ phpq pξ, ωq :" where N d is the number of degrees of freedom. Here, and in the following, the bold italic font denotes a vector field and the bold non-italic Roman font represents a matrix (upper case) or column vector (lower case). With this distinction, we rewrite (13) in matrix form as where Npξq is the matrix constructed with the basis vectors N k pξq as its columns, i.e.
With this, we may also rewrite (11) as follows and, with a suitable choice of l i pωq, we may rewrite (15) as the linear system of equations where the coefficients of Apωq and rpθ p0,hpq , ωq are defined to be NGSolve offers efficient approaches for the computational solution to (16) using preconditioned iterative solvers [36,19], which we exploit. Following the solution of (16), we can obtain θ p1,hpq pξ, ωq using (14) and, by repeating the process for i " 1, 2, 3, we get θ p1,hpq i pξ, ωq. Then pMrαB, ω, µ r , σ˚sq ij , for the full order model, is found by using (5).

Reduced order model (ROM)
A traditional approach for the computation of the MPT spectral signature, i.e. the variation of the coefficients MrαB, ωs with frequency, would involve the repeated solution of the N d sized system (16) for different ω. To reduce the computational cost of this, we wish to apply a ROM in which the solution of (16) is replaced by a surrogate problem of reduced size. Thus, reducing both the computation cost and time to produce a solution for each new ω. In particular, in Section 4.1, we describe a ROM based on the POD method [9,2,11,31] and, in Section 4.2, apply the variant called projection based POD (which we denote by PODP), which has already been shown to work well in the analysis of magneto-mechanical coupling applied to MRI scanners [31]. To emphasise the generality of the approach, the formulation is presented for an arbitrary list of problem parameters denoted by ω. In Section 4.3 we derive a procedure for computing certificates of accuracy on the ROM solutions with negligible additional cost.

Proper orthogonal decomposition
Following the solution of (16) for qpωq for different values of the set of parameters, ω, we construct a matrix D P C N dˆN with the vector of solution coefficients as its columns in the form where N ! N d denote the number of snapshots. Application of a singular value decomposition (SVD) e.g. [7,10] gives where U P C N dˆNd and V˚P C NˆN are unitary matrices and Σ P R N dˆN is a diagonal matrix enlarged by zeros so that it becomes rectangular. In the above, V˚" V T is the hermitian of V.
The diagonal entries pΣq ii " σ i 1 are the singular values of D and they are arranged as σ 1 ą σ 2 ą ... ą σ N . Based on the sparse representation of the solutions to (7) as function of ν, and hence ω, (and hence also the sparse representation of the MPT) found in [17], we expect these to decay rapidly towards zero, which motivates the introduction of the truncated singular value decomposition (TSVD) e.g. [7,10] where U M P C N dˆM are the first M columns of U, Σ M P R MˆM is a diagonal matrix containing the first M singular values and pV M q˚P C MˆN are the first M rows of V˚. The computation of (20) constitutes the off-line stage of the POD. Using (20) we can recover an approximate representation for each of our solution snapshots as follows where ppV M q˚q j refers to the jth column of pV M q˚.

Projection based proper orthogonal decomposition (PODP)
In the online stage of PODP, q P ODP pωq « qpωq is obtained by taking a linear combination of the columns of U M where the coefficients of this projection are contained in the vector p M . We choose to also approximate lpωq in a similar way so that θ p1,hpq pξ, ωq « pθ p1,hpq q PODP pξ, ωq :"Npξqq P ODP pωq " where P Y pP ODP q Ă Y ε X W phpq . Substituting these lower dimensional representations in to (15) we obtain the following Then, if we choose o M pωq appropriately, we obtain the linear system which is of size MˆM where A M pωq :" pU M q˚ApωqU M and r M pθ p0,hpq , ωq :" pU M q˚rpθ p0,hpq , ωq. Note, since M ă N ! N d , this is significantly smaller than (16) and, therefore, substantially computationally cheaper to solve. After solving this reduced system, and obtaining p M pωq, we obtain an approximate solution for θ p1,hpq pξ, ωq using (22). Focusing on the particular case where ω " ω, from (12) we observe that we can express A and r as the simple sums where the definitions of A p0q , A p1q and r p1q pθ p0,hpq q are obvious from (17), (12) and the definition of ν. Then, by computing and storing pU M q˚A p0q U M , pU M q˚A p1q U M , pU M q˚r p1q pθ p0,hpq q, which are independent of ω, it follows that A M pωq and r M pθ p0,hpq , ωq can be efficiently calculated for each new ω from the stored data. In a similar manner, by precomputing appropriate data, the MPT coefficients in (5) can also be rapidly evaluated for each new ω using the PODP solutions. This leads to further considerable computational savings. We emphasise that the PODP is only applied to obtain ROM solutions for θ p1q pξ, ωq and not to θ p0q pξq, which does not depend on ω.

PODP output certification
We follow the approach described in [11], which enables us to derive and compute certificates of accuracy on the MPT coefficients obtained with PODP, with respect to those obtained with full order model, as a function of ω. To do this, we set i pωq " θ p1,hpq i pωq´pθ p1,hpq i q PODP pωq P Y phpq , where we have reintroduced the subscript i, as we need to distinguish between the cases i " 1, 2, 3. Although i also depends on ξ, we have chosen here, and in the following, to only emphasise its dependence on ω. We have also introduced Y phpq " Y ε X W phpq for simplicity of notation, and note that this error satisfies which is called the error equation [11] and which is called Galerkin orthogonality [11]. The Riesz representation [11] of rp¨; θ p0,hpq i , ωq denoted bŷ r i pωq P Y phpq is such that so that Then, by using the alternative set of formulae for the tensor coefficients [17] pRrαB, ωsq ij "´α written in terms of the full order solutions, we obtain the certificates for the tensor entries computed using PODP stated in the lemma below. Note that the formulae stated in (5) are used for the actual POD computation of pR P ODP rαB, ωsq ij and pI P ODP rαB, ωsq ij , but the form in (29) is useful for obtaining certificates. Also, as pN rαBsq ij is independent of ω we have pN 0,P ODP rαBsq ij " pN 0 rαBsq ij and we write M P ODP rαB, ωs " N 0,P ODP rαBs`R P ODP rαB, ωs`iI P ODP rαB, ωs for the MPT obtained by PODP.
Lemma 4.1. An error certificate for the tensor coefficients computed using PODP išˇp where p∆rωsq ij :" and α LB is a lower bound on a stability constant.
Proof. We concentrate on the proof forˇˇpRrαB, ωsq ij´p R P ODP rαB, ωsq ijˇa s the proof for the second bound is similar and leads to the same result. Recalling the symmetry of RrαB, ωs, we have pRrαB, ωsq ij " 1 2 ppRrαB, ωsq ij`p RrαB, ωsq ji q so that which follows since ν and θ p0,hpq i are real valued and where we have dropped the dependence of ω on i for simplicity of presentation. Thus, Next, using (25), we make the observation thaťˇˇˇ Following similar steps to [11, pg47-50], and introducing a Riesz representation in (27), we can find that and, combining this with (31), completes the proof.
The efficient evaluation of (30) follows the approach presented in [11, pg52-54], adapted to complex matrices and with the simplification that we compute a Riesz representationr i pωq P Y ph0q using lowest order elements for computational efficiency. The computations are split in to those performed in the off-line stage and those in the on-line stage as follows.
In the off-line stage, the following p2M`1qˆp2M`1q hermitian matrices are computed where, since G pj,iq " pG pi,jq q H , it follows that, in practice, only the 3 matrices G p1,1q , G p2,2q and G p3,3q are required for computing the certificates on the diagonal entries of the tensors, and the further 3 matrices G p1,2q , G p1,3q and G p2,3q are needed for the off-diagonal terms. In the above, pM 0 q ij " N piq , N pjq are the coefficient of a real symmetric FEM mass matrix for the lowest order, with N piq , N pjq P W ph0q being typical lowest order basis functions, and where P p 0 is a projection matrix of the FEM basis functions from order p to the lowest order 0, U pM,iq is the U M obtained in (20) for the ith direction. The stability constant α LB " λ min minp1, ω ω 1 q is obtained from the smallest eigenvalue of an eigenvalue problem [11, pg56], which, in practice, is only performed once for smallest frequency of interest ω 1 .
In the on-line stage, we evaluate for each ω by updating the vector We then apply (30) to obtain the output certificates.

Scaling of the MPT under parameter changes
Two results that aid the computation of the frequency sweep of an MPT for an object with scaled conductivity and an object with a scaled object size from an already known frequency sweep of the MPT for the same shaped object are stated below.
Lemma 5.1. Given the MPT coefficients for an object αB with material parameters µ r and σ˚at frequency sω, the coefficients of the MPT for an object, which has the same B, α and µ r , but with conductivity sσ˚, at frequency ω, are given by pMrαB, ω, µ r , sσ˚sq ij "pMrαB, sω, µ r , σ˚sq ij , where pMrαB, sω, µ r , σ˚sq ij denote the coefficients of the original MPT at frequency sω.
Lemma 5.2. Given the MPT coefficients for an object αB with material parameters µ r and σ˚at frequency s 2 ω, the coefficients of the MPT for an object sαB, which is the same as B apart from having size sα, at frequency ω, are given by pMrsαB, ω, µ r , σ˚sq ij "s 3 pMrαB, s 2 ω, µ r , σ˚sq ij , where pMrαB, s 2 ω, µ r , σ˚sq ij denote the coefficients of the original MPT at frequency s 2 ω.

Numerical examples of PODP
The PODP algorithm has been implemented in the Python interface to the high order finite element solver NGSolve package led by the group of Schöberl [29,36,28] available at https://ngsolve.org. The snapshots are computed by solving (9) and (11) using NGSolve and their Hpcurlq conforming tetrahedral finite element basis functions of order p on meshes of spacing h [30]. Following the solution of (16), and the application of (13), the coefficients of MrαB, ωs 2 follow by simple post-processing of (5). If desired, PODP output certificates can also be efficiently computed using the approach described in Section 4.3. The Python scripts for the computations presented can be accessed at https://github.com/BAWilson94/MPT-Calculator.

Conducting permeable sphere
We begin with the case where B α " αB is a permeable conducting sphere of radius α " 0.01 m and B is the unit sphere centred at the origin. The sphere is chosen to have a relative permeability µ r " 1.5 and conductivity σ˚" 5.96ˆ10 6 S/m. To produce the snapshots of the full order model, we set Ω to be a ball 100 times the radius of B and discretize it using a mesh of 26 385 unstructured tetrahedral elements, refined towards the object, and a polynomial order of p " 3. We have chosen this discretization since it has already been found to produce an accurate representation of MrαB, ωs for 10 2 ă ω ă 10 8 rad/s by comparing with exact solution of the MPT for a sphere [34,16]. Indeed, provided that the geometry discretisation error is under control, performing p-refinement of the full order model solution results in exponential convergence to the true solution [14]. We follow two different schemes for choosing frequencies ω for generating the solution vectors qpωq required for D in (18). Firstly, we consider linearly spaced frequencies ω min ď ω n ď ω max , n " 1, 2, . . . , N , where, as in Section 4.1, N is the number of snapshots, and denote this choice of samples by "Lin" in the results. Secondly, we consider logarithmically spaced frequencies ω min ď ω n ď ω max and denote this regime by "Log" in the results.
Considering both linearly and logarithmically spaced frequencies with ω min " 1ˆ10 2 rad/s, ω max " 1ˆ10 8 rad/s and N " 9, 13, 17, in turn, to generate the snapshots, the application of an SVD to D in (19) leads to the results shown in Figure 3 where the values have been scaled by σ 1 and are strictly decreasing. We observe that "Log" case produces singular values σ i {σ 1 , which tend to 0 with increasing i, while the "Lin" case produces σ i {σ 1 , which tend to a finite constant with increasing i. Also shown is the tolerance T OL " 1ˆ10´3, i.e. we define M such that σ The superior performance of logarithmically spaced frequency snapshots over those linearly spaced is illustrated in Figure 4 paq where the variation of the condition number κpA M pωqq with ω for the frequency range and snapshots presented in Figure 3 is shown. Included in Figure 4 pbq is the corresponding error measure |epΛ i pωqq| :" |Λ exact i pωq´Λ P ODP i pωq|{|Λ exact i pωq| with ω, where Λ i pωq " λ i pRrαB, ωs`N 0 rαBsqì λ i pIrαB, ωsq, λ i p¨q indicates the ith eigenvalue. Note that, since the results for i " 1, 2, 3 are identical on this scale, only i " 1 is shown. From the results shown in this figure, we see that, with the exception of N " 17, the logarithmically spaced frequency snapshots results in a lower condition number compared to the linear ones and that all the logarithmically spaced snapshots result in a smaller error. Further tests reveal that the accuracy of the PODP using N " 9, 13, 17 and logarithmically spaced snapshots remains similar to that shown in Figure 4 pbq for T OL ď 1ˆ10´3 for this problem. We complete the discussion of the sphere by showing a comparison of λ i pRrαB, ωsq and λ i pIrαB, ωsq, each with ω, for the full order model, PODP using N " 9 and the exact solution in Figure 5. Again, the results for i " 1, 2, 3 are identical and, hence, only i " 1 is shown. In this figure, we observe excellent agreement between PODP, the full order model solution and exact solution. In Figure 6, we show the output certificates pR P ODP rαB, ωs`N 0,P ODP rαBsq ii˘p ∆rωsq ii (summation of repeated indices is not implied) and pI P ODP rαB, ωsq ii˘p ∆rωsq ii , each with ω, obtained by applying the technique described in Section 4.3 for the case where i " 1 and with N " 17, 21 and T OL " 11 0´6. Similar certificates can be obtained for the other tensor coefficients. We observe that certificates are indistinguishable from the the MPT coefficients obtained with PODP for low frequencies in both cases and the certificates rapidly tend to the MPT coefficients for all ω as N is increased. Note that T OL " 1ˆ10´6 is chosen as larger tolerances lead to larger certificates, however, this reduction in tolerance does not substantially affect the computational cost of the ROM. Although the effectivity indices p∆rωsq 11 {|pRrαB, ωs´R P ODP rαB, ωsq 11 | and p∆rωsq 11 {|pIrαB, ωs´I P ODP rαB, ωsq 11 |) of the PODP with respect to the full order model are clearly larger at higher frequencies, we emphasise that they are computed at negligible additional cost, they converge rapidly to the MPT coefficients obtained with PODP as N is increased and give credibility in the PODP solution without the need of performing additional full order model solutions to validate the ROM. for paq pN 0 rαBs`RrαB, ωsq 11 using N " 17, pbq pIrαB, ωsq 11 using N " 17, pcq pN 0 rαBs`RrαB, ωsq 11 using N " 21 and pdq pIrαB, ωsq 11 using N " 21, each with ω.
The computational speed-ups offered by using the PODP compared to a frequency sweep performed with the full order model are shown in Figure 7 where N " 9, 13, 17 and logarithmically spaced snapshots are chosen with ω min " 1ˆ10 2 rad/s, ω max " 1ˆ10 8 rad/s, as before. For the comparison, we vary the number of output points N 0 produced in a frequency sweep and measure the time taken to produce each of these frequency sweeps using a 2.9 GHz quad core Intel i5 processor and also show the percentage speed up offered by each of these PODP sweeps. Also shown is the break down of the computational time for the offline and online stages of the PODP for the case where N " 13. Note, in particular, that the computational cost increases very slowly with N 0 and that the additional cost involved in computing the output certificates is negligible. The breakdown of computational costs for other N is similar.

Conducing permeable torus
Next, we consider B α " αB to be a torus where B has major and minor radii a " 2 and b " 1, respectively, α " 0.01 m and the object is permeable and conducting with µ r " 1.5, σ˚" 5ˆ10 5 S/m. The object is centred at the origin so that it has rotational symmetry around the e 1 axis and hence MrαB, ωs has independent coefficients pMrαB, ωsq 11 and pMrαB, ωsq 22 " pMrαB, ωsq 33 , and thus N 0 rαBs, RrαB, ωs, IrαB, ωs each have 2 independent eigenvalues. To compute the full order model, we set Ω to be a sphere of radius 100, centred at the origin, such that it contains B and discretise it by a mesh of 26142 unstructured tetrahedral elements, refined towards the object, and a polynomial order of p " 3. This discretisation has already been found to produce an accurate representation of MrαB, ωs for the frequency range with ω min " 1ˆ10 2 rad{s and ω max " 1ˆ10 8 rad{s with the full order model. The reduced order model is constructed using N " 13 snapshots at logarithmically spaced frequencies with T OL " 1ˆ10´4. Figure 8 shows the results for λ i pN 0 rαBs`RrαB, ωsq and λ i pIrαB, ωsq, each with ω, for both the full order model and the PODP. The agreement is excellent in both cases. In Figure 9 we show the output certificates pR P ODP rαB, ωs`N 0,P ODP rαBsq ii˘p ∆rωsq ii (no summation over repeated indices implied) and pI P ODP rαB, ωsq ii˘p ∆rωsq ii , each with ω, obtained by applying the technique described in Section 4.3 for the case where N " 17 and T OL " 1ˆ10´6. Note that we increased the number of snapshots from N " 13 to N " 17 and have reduced the tolerance to ensure tight certificates bounds.

Conducting permeable tetrahedron
The third object considered is where B α " αB is a conducting permeable tetrahedron. The vertices of the tetrahedron B are chosen to be at the locations v 1 "¨0 0 0‚ , v 2 "¨7 0 0‚ , v 3 "¨5 .5 4.6 0‚ and v 4 "¨3 , the object size is α " 0.01 m and the tetrahedron is permeable and conducting with µ r " 2 and σ˚" 5.96ˆ10 6 S/m. The object does not have rotational or reflectional symmetries and, hence, MrαB, ωs has 6 independent coefficients and, thus, N 0 rαBs, RrαB, ωs, IrαB, ωs each have 3 independent eigenvalues. To compute the full order model, we set Ω to be a cube with sides of length 200 centred about the origin and discretise it with a mesh of 21427 unstructured tetrahedral elements, refined towards the object, and a polynomial order of p " 3. This discretisation has already been found to produce an accurate representation of MrαB, ωs for the frequency range with ω min " 1ˆ10 2 rad{s and ω max " 1ˆ10 8 rad{s.
The reduced order model is constructed using N " 13 snapshots at logarithmically spaced frequencies with T OL " 1ˆ10´4. Figure 10 shows the results for λ i pN 0 rαBs`RrαB, ωsq and λ i pIrαB, ωsq, each with ω, for both the full order model and the PODP. The agreement is excellent in both cases. In Figure 11 we show the output certificates pR P ODP rαB, ωs`N 0,P ODP rαBsq ij˘p ∆rωsq ij and pI P ODP rαB, ωsq ij˘p ∆rωsq ij , both with ω, for i " j and i ‰ j obtained by applying the technique described in Section 4.3 for the case where N " 21 and T OL " 1ˆ10´6. Once again, we increased the number of snapshots from N " 13 to N " 21 and have reduced the tolerance to ensure tight certificates bounds, except at large frequencies. i " j pbq pIrαB, ωsq ij , with i " j, pcq pN 0 rαBs`RrαB, ωsq ij , with i ‰ j, pdq pIrαB, ωsq ij with i ‰ j, each with ω.

Inhomogeneous conducting bar
As a final example we consider B α " αB to the inhomogeneous conducting bar made up from two different conducting materials. The size, shape and materials of this object are the same as those presented in Section 6.1.3 of [18]. This object has rotational and reflectional symmetries such that MrαB, ωs has independent coefficients pMrαB, ωsq 11 , pMrαB, ωsq 22 " pMrαB, ωsq 33 and, thus, N 0 rαBs, RrαB, ωs, IrαB, ωs each have 2 independent eigenvalues. To compute the full order model, we set Ω to be a sphere of radius 100 centred about the origin and discretise it with a mesh of 30209 unstructured tetrahedral elements, refined towards the object, and a polynomial order of p " 3. This discretisation has already been found to produce an accurate representation of MrαB, ωs for the frequency range with ω min " 1ˆ10 2 rad{s and ω max " 1ˆ10 8 rad{s. The reduced order model is constructed using N " 13 snapshots at logarithmically spaced frequencies with T OL " 1ˆ10´4. Figure 10 shows the results for λ i pN 0 rαBs`RrαB, ωsq and λ i pIrαB, ωsq, each with ω, for both the full order model and the PODP. The agreement is excellent in both cases. In Figure 13 we show the output certificates pR P ODP rαB, ωs`N 0,P ODP rαBsq ii˘p ∆rωsq ii (no summation over repeated indices implied) and pI P ODP rαB, ωsq ii˘p ∆rωsq ii , both with ω, obtained by applying the technique described in Section 4.3 for the case where N " 23 and T OL " 1ˆ10´6. Note that we increased the number of snapshots from N " 13 to N " 23 and have reduced the tolerance to ensure tight certificates bounds, except at large frequencies.

Numerical examples of scaling
In this section we illustrate the application of the results presented in Section 5.

Scaling of conductivity
As an illustration of Lemma 5.1, we consider a conducting permeable sphere B α " αB where α " 0.01 m with materials properties µ r " 1.5 and σ p1q " 1ˆ10 7 S/m and a second object, which is the same as the first except that σ p2q " sσ p1q " 10σ p1q . In Figure 14, we compare the full order computations of MrαB, ω, µ r , σ p1q s and MrαB, ω, µ r , σ p2q s with that obtained from (32). We observe that the translation predicted by (32) is in excellent agreement with the full order model solution for MrαB, ω, µ r , σ p2q s.

Scaling of object size
To illustrate Lemma 5.2, we consider a conducting permeable tetrahedron B p1q α " α p1q B " 0.01B with vertices as described in Section 6.3 and material properties µ r " 1.5 and σ˚" 1ˆ10 6 S/m. Then, we consider a second object B p2q α " α p2q B " sα p1q B " 0.015B, which, apart from its size, is otherwise the same as B p1q α . In Figure 15, we compare the full order computations of Mrα p1q B, ω, µ r , σ˚s and Mrα p2q B, ω, µ r , σ˚s with that obtained from (33). We observe that the translation and scaling predicted by (33) is in excellent agreement with the full order model solution for Mrα p2q B, ω, µ r , σ˚s. " α p1q B " 0.01B with µr " 1.5 and σ˚" 1ˆ10 6 S/m, α " 0.01 m and a second tetrahedron, which is the same as the first except that B p2q α " α p2q B " sα p1q B " 0.015B: showing the translation and scaling predicted by (33) compared with the full order model solutions for paq λ i pN 0 rαB, µrs`RrαB, ω, µr, σ˚sq and pbq λ i pIrαB, ω, µr, σ˚sq.

Conclusions
An application of a ROM using PODP for the efficient computation of the spectral signature of the MPT has been studied in this paper. The full order model has been approximated by Hpcurlq conforming discretisation using the NGSolve finite element package. The offline stage of the ROM involves computing a small number of snapshots of the full order model at logarithmically spaced frequencies, then in the online stage, the spectral signature of the MPT is rapidly and accurately predicted to arbitrarily fine fidelity using PODP. Output certificates have been derived and can be computed in the online stage at negligible computational cost and ensure accuracy of the ROM prediction. If desired, these output certificates could be used to drive an adaptive procedure for choosing new snapshots, in a similar manner to the approach presented in [11]. However, by choosing the frequency snapshots logarithmically, accurate spectral signatures of the MPT were already obtained with tight certificate bounds. In addition, simple scaling results, which enable the MPT spectral signature to be easily computed from an existing set of coefficients under the scaling of an object's conductivity or object size, have been derived. A series of numerical examples have been presented to demonstrate the accuracy and efficiency of our approach for homogeneous and inhomogeneous conducting permeable objects. Future work involves applying the presented approach to generate a dictionary of MPT spectral signatures for different objects for the purpose of metallic object identification using a classifier.