Sergiy Pereverzyev
email: sergiy.pereverzyev@uibk.ac.at

Kateřina Hlaváčková-Schindler

Graphical Lasso Granger method with 2-levels-thresholding for recovering causality networks

Keywords: causality network, gene causality network, Granger causality, graphical Lasso method, 2-levels-thresholding

The recovery of the causality networks with a number of variables is an important problem that arises in various scientific contexts. For detecting the causal relationships in the network with a big number of variables, the so called Graphical Lasso Granger (GLG) method was proposed. It is widely believed that the GLG-method tends to overselect causal relationships. In this paper, we propose a thresholding strategy for the GLG-method, which we call 2-levels-thresholding, and we show that with this strategy the variable overselection of the GLG-method may be overcomed. Moreover, we demonstrate that the GLG-method with the proposed thresholding strategy may become superior to other methods that were proposed for the recovery of the causality networks.

Introduction

Causality is a relationship between a cause and its effect (its consequence). One can say that inverse problems solving, where one would like to discover unobservable features of the cause from the observable features of an effect [START_REF] Engl | Regularization of inverse problems[END_REF], i.e., searching for the cause of an effect, can in general be seen as a causality problem.

A causality network is a directed graph with nodes, which are variables {x j , j = 1, . . . , p}, and directed edges, which are the causal influences between the variables. We write x i ← x j if the variable x j has a causal influence on the variable x i . Causality networks arise in various scientific contexts.

For example, In Cell Biology one considers causality networks which involve sets of active genes of a cell. An active gene produces a protein. It has been observed that the amount of the protein which is produced by a given gene may depend on, or may be causally influenced by, the amount of the proteins which are produced by other genes. In this way, causal relationships between genes and the corresponding causality network arise. These causality networks are also called gene regulatory networks. An example of such a network is presented in Fig. 1. This network is achieved from the biological experiments in [START_REF] Li | Discovery of time-delayed gene regulatory networks based on temporal gene expression profiling[END_REF], and it can be found in the BioGRID database. This network has been used in several works [START_REF] Sambo | CNET: an algorithm for reverse engineering of causal gene networks[END_REF][START_REF] Lozano | Grouped graphical Granger modeling for gene expression regulatory networks discovery[END_REF][START_REF] Shojaie | Discovering graphical Granger causality using the truncating lasso penalty[END_REF] as a test network. Knowledge of the correct causality networks is important for changing them. In Cell Biology, these networks are used in the research of the causes of genetic diseases. For example, the network in Fig. 1 consists of genes that are active in the human cancer cell HeLa [START_REF] Whitfield | Identification of genes periodically expressed in the human cell cycle and their expression in tumors[END_REF]. If one wants to suppress the genes expression in this network, then the primary focus of the suppression therapy should be on the causing genes. For the use of the causality networks in other sciences see, for example, [START_REF] Pereverzyev | Graphical Lasso Granger method with 2-levels-thresholding for recovering causality networks[END_REF].

How can causality network be recovered? In practice, the first information that can be known about the network is the time evolution (time series) of the involved variables {x j t , t = 1, . . . , T }. How can this information be used for inferring causal relationships between the variables?

The statistical approach to the derivation of the causal relationships between a variable y and variables {z j , j = 1, . . . , p} using the known time evolution of their values {y t , z j t , t = 1, . . . , T, j = 1, . . . , p} consists in considering a model of the relationship between y and {z j , j = 1, . . . , p}. As a first step, one can consider a linear model of this relationship: y t ≈ p j=1 β j z j t , t = 1, . . . , T. The coefficients {β j , j = 1, . . . , p} can be specified using the least-squares method. Then, in Statistics [START_REF] Wikipedia | Causality -Wikipedia, The Free Encyclopedia[END_REF] by fixing the value of a threshold parameter β tr > 0, one says that there is a causal relationship y ← z j if |β j | > β tr .

For detecting causal relationships between variables {x j , j = 1, . . . , p} the concept of the so called multivariate Granger causality has been proposed. This concept originated in the work of Clive Granger [START_REF] Granger | Investigating causal relations by econometric models and crossspectral methods[END_REF], who was awarded the Nobel Prize in Economic Sciences in 2003. Based on the intuition that the cause should precede its effect, in Granger causality one says that a variable x i can be potentially caused by the past versions of the involved variables {x j , j = 1, . . . , p}.

Then, in the spirit of the statistical approach and using a linear model for the causal relationship, we consider the following approximation problem:

x i t ≈ p j=1 L l=1 β j l x j t-l , t = L + 1, . . . , T, (1)
where L is the so called maximal lag, which is the maximal number of the considered past versions of the variables. The coefficients {β j l } can be determined by the least-squares method. As in the statistical approach, one can now fix the value of the threshold parameter β tr > 0 and say that

x i ← x j if L l=1 |β j l | > β tr . (2)
It is well known that for a big number of genes p, as it is pointed out for example in [START_REF] Lozano | Grouped graphical Granger modeling for gene expression regulatory networks discovery[END_REF], the causality network, which is obtained from the approximation problem (1), is not satisfactory. First of all, it cannot be guaranteed that the solution of the corresponding minimization problem is unique. Another issue is connected with the number of the causality relationships that is obtained from (1). This number is typically very big, while one expects to have a few causality relationships with a given gene. To address this issue, various variable selection procedures can be employed. The Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF] is a well known example of such a procedure. In the regularization theory, this approach is known as the l 1 -Tikhonov regularization. It has been extensively used for reconstructing the sparse structure of an unknown signal. We refer the interesting reader to [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF][START_REF] Ramlau | A Tikhonov-based projection iteration for nonlinear ill-posed problems with sparsity constraints[END_REF][START_REF] Grasmair | Sparse regularization with l q penalty term[END_REF][START_REF] Fornasier | Theoretical foundations and numerical methods for sparse recovery[END_REF][START_REF] Lorenz | Gradient descent for Tikhonov functionals with sparsity constraints: theory and numerical comparison of step size rules[END_REF] and the references therein.

The causality concept that is based on the Lasso was proposed in [START_REF] Arnold | Temporal causal modeling with graphical Granger methods[END_REF] and is named Graphical Lasso Granger (GLG) method. However, it is stated in the literature that the Lasso suffers from the variable overselection. And therefore, in the context of the gene causality networks several Lasso modifications were proposed. In [START_REF] Lozano | Grouped graphical Granger modeling for gene expression regulatory networks discovery[END_REF], the so called group Lasso method was considered for recovering gene causality networks using the multivariate Granger causality. The corresponding method can be named Graphical group Lasso Granger (GgrLG) method. And in [START_REF] Shojaie | Discovering graphical Granger causality using the truncating lasso penalty[END_REF], the truncating Lasso method was proposed. The resulting method can be named Graphical truncating Lasso Granger (GtrLG) method.

Nevertheless, it seems that an important tuning possibility of the Lasso, namely an appropriate choice of the threshold parameter β tr , has been overlooked in the literature devoted to the recovery of the gene causality networks. In this paper, we are going to show that the GLG-method, which is equipped with an appropriate thresholding strategy and an appropriate regularization parameter choice rule, may become a superior method in comparison to other methods that were proposed for the recovery of the gene causality networks.

The paper is organized as follows. In Sect. 2, we recall the GLG-method. The quality measures of the graphical methods are presented in Sect. 3. In Sect. 4, we use the network from Fig. 1 to compare the performance of the known graphical methods with the ideal version of the GLG-method, which we call the optimal GLG-estimator. Such a comparison demonstrates the potential of the GLGapproach. In Sect. 5, we propose a thresholding strategy for the GLG-method that allows its automatic realization, which we describe in Sect. 6. Then again we use the network from Fig. 1 to compare the performance of the proposed version of the GLG-method with other graphical methods. It turns out that the proposed method has a superior quality compared to the known methods. The paper is finished with the conclusion and outlook in Sect. 7.

Graphical Lasso Granger method

Let us specify the application of the least-squares method to the approximation problem [START_REF] Arnold | Temporal causal modeling with graphical Granger methods[END_REF]. For this purpose, let us define the vectors

Y i = (x i L+1 , x i L+2 , . . . , x i T ) , β = (β 1 1 , . . . , β 1 L , β 2 1 , . . . , β 2 L , . . . , β p 1 , . . . , β p L )
, and the matrix

X = (x 1 t-1 , . . . , x 1 t-L , x 2 t-1 , . . . , x 2 t-L , . . . , x p t-1 , . . . , x p t-L ); t = L + 1, . . . , T .
Then, in the least-squares method, one considers the following minimization problem:

Y i -Xβ 2 → min β , (3)
where • denotes the l 2 -norm.

As it was mentioned in the introduction, the solution of (3) defines unsatisfactory causal relationships and various variable selection procedures should be employed instead. A well-known example of such procedures is the Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF]. In this procedure, one considers the following minimization problem:

Y i -Xβ 2 + λ β 1 → min β . (4)
Solution of (4) for each variable {x i , i = 1, . . . , p} with the causality rule (2) defines an estimator of the causality network between the variables {x i }, and in this way one obtains the Graphical Lasso Granger (GLG) method [START_REF] Arnold | Temporal causal modeling with graphical Granger methods[END_REF].

Quality measures of the graphical methods

A graphical method is a method that reconstructs the causality network, which is a directed graph, with the variables {x j }. The quality of a graphical method can be estimated from its performance on a known causality network. The network in Fig. 1 has been used for testing methods' quality in several publications [START_REF] Sambo | CNET: an algorithm for reverse engineering of causal gene networks[END_REF][START_REF] Lozano | Grouped graphical Granger modeling for gene expression regulatory networks discovery[END_REF][START_REF] Shojaie | Discovering graphical Granger causality using the truncating lasso penalty[END_REF]. What measures can be used for estimating the quality of a graphical method?

First of all, let us note that a causality network can be characterized by the so called adjacency matrix A = {A i,j | {i, j} ⊂ {1, . . . , p} } with the following elements:

A i,j = 1 if x i ← x j ; A i,j = 0 otherwise.

The adjacency matrix A true for the causality network in Fig. 1 is presented in Fig. 2. There, the white squares correspond to A i,j = 1, and the black squares -to the zero-elements. The genes are numbered in the following order: CDC2, CDC6, CDKN3, E2F1, PCNA, RFC4, CCNA2, CCNB1, CCNE1. Now, imagine that there is a true adjacency matrix A true of the true causality network, and there is its estimator A estim , which is produced by a graphical method. The quality of the estimator A estim can be characterized by the following quality measures: precision (P), recall (R), F 1 -score (F 1 ). See, for example, [START_REF] Pereverzyev | Graphical Lasso Granger method with 2-levels-thresholding for recovering causality networks[END_REF] for the detailed definition of these measures.

As it was already mentioned, the causality network in Fig. 1 has been used for testing quality of graphical methods. In particular, in [START_REF] Shojaie | Discovering graphical Granger causality using the truncating lasso penalty[END_REF] one finds the above mentioned quality measures for the following methods: GgrLG, GtrLG and CNET. CNET is a graph search-based algorithm that was introduced in [START_REF] Sambo | CNET: an algorithm for reverse engineering of causal gene networks[END_REF]. The data {x j t } is taken from the third experiment of [START_REF] Whitfield | Identification of genes periodically expressed in the human cell cycle and their expression in tumors[END_REF] consisting of 47 time points, and the maximal lag L is taken to be equal to 3. The quality measures from [START_REF] Shojaie | Discovering graphical Granger causality using the truncating lasso penalty[END_REF] are presented in Table 1. As it is seen from the table, CNET has the highest F 1 -score. However, CNET is the most computationally expensive among the considered methods that does not allow its application to large networks. GgrLG has a good recall but a poor precision, and thus, GtrLG can be considered as a better method among the considered methods.

Optimal GLG-estimator

As we have seen in the previous section, the graphical methods, which are based on the Lasso modifications, were tested on the network in Fig. 1 (Table 1). However, the application of the graphical method that is based on the pure Lasso (GLG) to the network in Fig. 1 has not been reported. Moreover, it seems that the possibility of varying the threshold parameter β tr in GLG also has not been considered in the literature devoted to the reconstruction of the causality networks.

Assume that the true causality network with the variables {x j } is given by the adjacency matrix A true . Assume further that the observation data {x j t } is given. What is the best reconstruction of A true that can be achieved by the GLG-method? The answer to this question is given by, what we call, the optimal GLG-estimator. Let us specify its construction.

First of all, let us define the following quality measure, which we call F smeasure: F s = 1 p 2 A true -A estim 1 , 0 ≤ F s ≤ 1. F s-measure represents the number of false elements in the estimator A estim that is scaled with the total number of elements in A estim . Now, let β i (λ) denote the solution of the minimization problem (4) in the GLG-method, and β j i (λ) = (β j 1,i , . . . , β j L,i ). Then, the GLG-estimator A GLG (λ, β tr ) of the adjacency matrix A true is defined as follows:

A GLG i,j (λ, β tr ) = 1 if β j i (λ) 1 > β tr ; A GLG i,j (λ, β tr ) = 0 otherwise.
The optimal GLG-estimator A GLG,opt of the true adjacency matrix A true is the GLG-estimator A GLG (λ, β tr ) with the parameters λ, β tr such that the corresponding F s-measure is minimal.

The optimal GLG-estimator of the adjacency matrix for the causality network in Fig. 1 is presented in Fig. 2. Its quality measures can be found in Table 2. We used the same data {x j t } as in [START_REF] Sambo | CNET: an algorithm for reverse engineering of causal gene networks[END_REF][START_REF] Lozano | Grouped graphical Granger modeling for gene expression regulatory networks discovery[END_REF][START_REF] Shojaie | Discovering graphical Granger causality using the truncating lasso penalty[END_REF]. Also, as in [START_REF] Lozano | Grouped graphical Granger modeling for gene expression regulatory networks discovery[END_REF][START_REF] Shojaie | Discovering graphical Granger causality using the truncating lasso penalty[END_REF], we take the maximal lag L = 3. As one can see, the optimal GLG-estimator reconstructs almost completely the causing genes of the most caused gene in the network. The recall of A GLG,opt is equal to the highest recall in Table 1, but precision and F 1 -score are considerably higher.

Of course, A GLG,opt is given by the ideal version of the GLG-method, where we essentially use the knowledge of A true . How close can we come to A GLG,opt without such a knowledge? To answer this question, let us first decide about the choice of the threshold parameter β tr .

Thresholding strategy

The purpose of the threshold parameter β tr is to cancel the causal relationships x i ← x j with small β j i (λ) 1 . When can we say that β j i (λ) 1 is small? We propose to consider the following guideindicators of smallness:

β i min (λ) = min{ β j i (λ) 1 , j = 1, . . . , p | β j i (λ) 1 = 0 }, β i max (λ) = max{ β j i (λ) 1 , j = 1, . . . , p }. (5)
In particular, we propose to consider the threshold parameter of the following form:

β i tr,α (λ) = β i min (λ) + α( β i max (λ) -β i min (λ) ). (6)
As a default value we take α = 1/2. In the optimal GLG-estimator A GLG,opt tr,1/2 with the threshold parameter

β i tr,1/2
we choose λ such that the corresponding F s-measure is minimal. For the causality network in Fig. 1, this estimator is presented in Fig. 2. Its quality measures can be found in Table 2. One observes that although there is some quality decrease in comparison to A GLG,opt , the quality measures are still higher than for the methods in Table 1. However, can this quality be improved?

The choice of the threshold parameter β i tr,1/2 rises the following issue. With such a choice we always assign a causal relationship, unless the solution of (4) β i (λ) is identically zero. But how strong are these causal relationships compared to each other? The norm β j i (λ) 1 can be seen as a strongness indicator of the causal relationship x i ← x j .

Let us now construct a matrix A GLG,opt;β tr,1/2

, similarly to the adjacency matrix

A GLG,opt tr,1/2
, where instead of the element 1 we put the norm β j i (λ) 1 , i.e.

A GLG,opt;β tr,1/2

(i, j) = β j i (λ tr,1/2 opt,i ) 1 if β j i (λ tr,1/2 opt,i ) 1 > β i tr,1/2 , A GLG,opt;β tr,1/2 (i, j) =0 otherwise.
This matrix is presented in Fig. 2. One observes that the false causal relationships of the estimator A GLG,opt tr,1/2 are actually weak. This observation suggests to use a second thresholding that is done on the network, or adjacency matrix, level.

We propose to do the thresholding on the network level similarly to the thresholding on the gene level. Namely, let us define the guideindicators of smallness on the network level similarly to [START_REF] Fornasier | Theoretical foundations and numerical methods for sparse recovery[END_REF]:

A min = min{ A GLG,opt;β tr,1/2 (i, j) = 0 }, A max = max{ A GLG,opt;β tr,1/2 (i, j) }.
And similarly to [START_REF] Granger | Investigating causal relations by econometric models and crossspectral methods[END_REF], define the threshold on the network level as follows:

A tr,α1 = A min + α 1 ( A max -A min ). (7)
We find it suitable to call the described combination of the two thresholdings on the gene and network levels as 2-levels-thresholding. The adjacency matrix obtained by this thresholding strategy is the following:

A GLG,opt
tr,1/2;α1 (i, j) =1 if A GLG,opt;β tr,1/2 (i, j) > A tr,α1 , A GLG,opt tr,1/2;α1 (i, j) =0 otherwise.

It turns out that with α 1 = 1/4 in [START_REF] Grasmair | Sparse regularization with l q penalty term[END_REF] the optimal GLG-estimator can be fully recovered.

An automatic realization of the GLG-method

For an automatic realization of the GLG-method, i.e. a realization that does not rely on the knowledge of the true adjacency matrix A true , in addition to a thresholding strategy one needs a choice rule for the regularization parameter λ in (4). For such a choice, we propose to use the so called quasi-optimality criterion [START_REF] Tikhonov | Use of the regularization method in non-linear problems[END_REF][START_REF] Bauer | Regularization independent of the noise level: an analysis of quasi-optimality[END_REF][START_REF] Kindermann | On the convergence of the quasioptimality criterion for (iterated) Tikhonov regularization[END_REF]. Some details of the application of this criterion can be found in [START_REF] Pereverzyev | Graphical Lasso Granger method with 2-levels-thresholding for recovering causality networks[END_REF].

The reconstruction obtained by the GLG-method with the 2-levels-thresholding and quasi-optimality criterion A GLG,qo tr,1/2;1/4 is presented in Fig. 2. Its quality measures can be found in Table 2. One observes that there is a little decrease in recall in comparison to the optimal GLG-method; however, this recall is the same as for the GtrLG-method (Table 1). But due to the highest precision, the F 1 -score remains to be higher than for the methods in Table 1. Thus, one may say that the proposed realization of the GLG-method outperforms the methods in Table 1.

Nevertheless, one may still wonder, why the proposed realization of the GLGmethod captures only the causal relationships of the most caused gene. It appears that the value of the maximal lag L plays an important role in the selection of the causal relationships.

In the modifications of the GLG-method the authors of [START_REF] Lozano | Grouped graphical Granger modeling for gene expression regulatory networks discovery[END_REF][START_REF] Shojaie | Discovering graphical Granger causality using the truncating lasso penalty[END_REF] considered L = 3. All results presented so far were also obtained with L = 3. It turns out that for L = 4 the optimal GLG-estimator (see Fig. 2) delivers a much better reconstruction of the causality network. In particular, two more caused genes are recovered.

The proposed automatic realization of the GLG-method with L = 4 (Fig. 2) recovers an additional caused gene in comparison to the realization with L = 3. Also, all considered quality measures for our automatic realization of the GLGmethod with L = 4 (Table 2) are considerably higher than for the methods in Table 1. We would like to stress that no use of the knowledge of A true is needed for obtaining A GLG,qo tr,1/2;1/4 , and no readjustment of the design parameters α, α 1 is necessary.

Conclusion and outlook

The proposed realization of the Graphical Lasso Granger method with 2-levelsthresholding and quasi-optimality criterion for the choice of the regularization parameter shows a considerable improvement of the reconstruction quality in comparison to other graphical methods. So, the proposed realization is a very promising method for recovering causality networks. Further tests and developments of the proposed realization are worthwhile. In particular, applications to larger causality networks are of interest. A GLG,opt tr,1/2;1/4

A GLG,qo tr,1/2;1/4 As an open problem for the future, one could consider a study of the choice of the maximal lag and its possible variation with respect to the caused and causing genes.

Fig. 1 .

 1 Fig. 1. Causality network of the human cancer cell HeLa genes from the BioGRID database (www.thebiogrid.org).

9 AFig. 2 .

 92 Fig.2. The adjacency matrix A true for the causality network in Fig.1and its various GLG-estimators.

Table 1 .

 1 Quality measures of the known graphical methods.

		P	R	F1
	GgrLG	0.24	0.44	0.3
	GtrLG	0.3	0.33	0.32
	CNET	0.36	0.44	0.4

Table 2 .

 2 Quality measures of the various GLG-estimators.

		Fs	P	R	F1
	GLG-opt	6.2%	1	0.44	0.62
	GLG-opt; tr, 1/2	14.8%	0.38	0.56	0.45
	GLG-qo; tr, 1/2; 1/4	7.4%	1	0.33	0.5
	GLG-opt, L = 4	3.7%	0.88	0.78	0.82
	GLG-qo;L = 4;tr, 1/2; 1/4	7.4%	0.71	0.56	0.63


François Léonard

Noisy data clusters are hollow

Keywords: Shell manifold, hollow clusters, high dimensional data, Big Data, measurement, multidimensional statistic, multivariate analysis, equipment monitoring

A new vision in multidimensional statistics is proposed impacting several areas of application. In these applications, a set of noisy measurements characterizing the repeatable response of a process is known as a realization and can be seen as a single point in . The projections of this point on the N axes correspond to the N measurements. The contemporary vision of a diffuse cloud of realizations distributed in is replaced by a cloud in the shape of a shell surrounding a topological manifold. This manifold corresponds to the process's stabilized-response domain observed without the measurement noise. The measurement noise, which accumulates over several dimensions, distances each realization from the manifold. The probability density function (PDF) of the realization-to-manifold distance creates the shell. Considering the central limit theorem as the number of dimensions increases, the PDF tends toward the

where µ fixes the center shell location and  fixes the shell thickness. In the proposed vision, the likelihood of a realization is a function of the realization-to-shell distance rather than the realization-to-manifold distance. The demonstration begins with the work of Claude Shannon followed by the introduction of the shell manifold and ends with practical applications to monitoring equipment.

Introduction

This paper focuses on a link between the measurement noise of a repeatable process and hollow clusters in and proposes a method to exploit this link in different applications. The proposed link has been only partially revealed by papers about fuzzy clustering, shell clustering, noise clusters, hollow clusters in imaging, vacuum clusters, data topology, topological noise and persistent homology.

In a generalization of the fuzzy k-means clustering algorithm, Coray [START_REF] Coray | Clustering algorithms with prototype selection[END_REF] suggested the use of non-linear circular prototypes for clusters. Some shell-clustering algorithms and their application to boundary detection and surface approximation have been proposed by Krishnapuram et al. [START_REF] Krishnapuram | Fuzzy and Possibilistic Shell Clustering Algorithms and their Application to Boundary Detection and Surface Approximation-Part I[END_REF] with the concept of expected thickness of the shell clusters. Some authors [START_REF] Ester | A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise[END_REF][START_REF] Ertöz | Finding Clusters of Different Sizes, Shapes, and Densities in Noisy, High Dimensional Data[END_REF] propose clustering tools adapted to noisy data. For example, a noise cluster prototype is introduced in the hope that all the noisy points can be part of the noise cluster [START_REF] Ohashi | Fuzzy Clustering and Robust Estimation[END_REF][START_REF] Davé | Characterization and Detection of Noise in Clustering[END_REF]. Hollow clusters have been observed in imaging [START_REF] Lee | The nonlinear statistics of highcontrast patches in natural images[END_REF]. In patternrecognition techniques, vacuum shell detection is used in clustering methods [START_REF] Tso | Classification Methods for Remotely Sensed Data[END_REF].

The topic of data topology is linked to many theoretical aspects and application fields. Hollow shapes can be detected by homology, which is a computable topological invariant of shapes. Niyogi et al. [START_REF] Niyogi | Finding the homology of submanifolds with high confidence from random samples[END_REF] show that if the data is noisy, in the sense that it is drawn from a probability distribution concentrated around a manifold, the homology of that manifold can still be computed with high confidence from that data. In topology, the noisy data is handled more systematically by the persistent homology methods developed recently by several authors [START_REF] Edelsbrunner | Topological persistence and simplification[END_REF][START_REF] Hanson | On the Strengthening of Topological Signals in Persistent Homology through Vector Bundle Based Maps[END_REF][START_REF] Chazal | Weak feature size and persistent homology: computing homology of solids in R n from noisy data samples[END_REF]. Pekelis and Holmes [START_REF] Pekelis | Are the data hollow?[END_REF] present statistical approaches to calculating the homology of data and testing whether data lie on a closed manifold with a hollow interior. The reader is referred to Carlsson [START_REF] Carlsson | Topology and Data[END_REF] for a comprehensive review of this topic.

The link between measurement noise and hollow clusters starts with the definition of repeatable process. A repeatable process offers an opportunity for pattern recognition, monitoring, trending, diagnosis and prognosis. The time dimension is implicit in the context of monitoring but not necessary in pattern recognition. Repeatable processes are found in many fields: medical, biological, chemical, mechanical, electrical, quantum mechanics, social, stock market, manufacturing, communication and imaging. A process can be defined as a transfer function having internal states, external inputs (independent variables) and external measured responses (dependent variables). A process can be observed from the measurements of independent and dependent variables. Measurements contain both information and noise. With regard to information, in some specific cases, e.g. communications, the information set is known and defined by a finite set of discrete values. In more common cases, the information is unknown and can be defined by a multidimensional function corresponding to a unique internal state (Fig. 1a). As for noise, this is a zero-mean random process, correlated or not between measurements. Noise can be seen as the sum of numerous sources of irrecoverable information. Hidden state processes with a small impact contribute to noise. When the time dimension is implicit, a realization is a set of close time measurements recorded on a same channel or on many channels. In practical terms, the operating condition of the process is time-variant and, as in photography, the realization is impaired by a blur when the time spread of the measurement set is too large. In the following development, a realization corresponds to a picture of a stabilized process and can be taken discontinuously with an irregular sampling rate. Note that the "time" dimension can be replaced by the appropriate dimension related to the context (e.g. spatial). The measurement set corresponding to a realization can share a common engineering unit or not. The normalization of measurements having different engineering units is an item of great concern: the selected normalization method contributes to the success or failure of the application algorithm (monitoring, pattern recognition…). Turning lastly to internal states, a response function can contain more than one internal state which can be visible or not from the dependent variables and can be determined by a sequence of independent variable values (Markov chain) or hidden. Figure 1 illustrates the response function corresponding to a single, several and multiple internal states. A hidden state can be defined as a state not linked to the input variables with causality. In the latter case, an unmeasured independent variable which contributes inputs to the process can be transposed as a hidden state of the process. The additive effect of numerous small contributions of unmeasured independent variables and hidden states to the process yields a diffused response function (Fig. 1d).

In multidimensional statistics, a realization corresponds to a single point in . A set of realizations corresponding to the same state process yields a cluster in . This paper starts with the Shannon vision of the cluster constellation corresponding to the received signals: a stack of n-spheres with n=N-1, each n-sphere centered on its corresponding transmitted symbol before the addition of noise. From this vision of multiple n-sphere clusters associated with already known symbols having no dimension (a point in for each state process), we will jump to a single cluster generated by an unknown multivariate one-state process. In both cases, the clusters are hollow.

Back to Shannon

Hollow clusters of measured data were first introduced by Claude Shannon [START_REF] Shannon | Communication in the Presence of Noise[END_REF] in the mid 20th century. His demonstration of the limit capacity of a communication channel is based on the limit of the stacking of n-spheres in a volume delimited by a maximum amplitude signal. The n-sphere radius corresponds to the average length of the noise vector added to a transmitted symbol. The error rate appears as a function of the overlap between the shells, each shell corresponding to a transmitted symbol. In Shannon's demonstration, a communication channel comprises an information source (the already known symbol set), a transmitter, a noise source, a receiver and a destination. Consider only one symbol

  1 , , N x x  x  (1)
of the set transmitted using N dimensions. The corresponding received signal is

  1 1 1 , , , , n N N x x x          y x ε   (2)
where n  is the noise added for the dimension "n" and ε the noise vector added to the transmitted symbol in . The average distance from the signal received to the symbol transmitted

        y x (3)
corresponds to the radius of the n-sphere and becomes

2 1 N n n               (4)
considering the use of the Euclidian metric. The measurement dispersion

  2 2           y x (5)
is called the hardness of the n-sphere. The increasing µ/ ratio with the increasing number of dimensions is called the "sphere-hardening phenomenon" in the communication theory [START_REF] Willems | Information and communication theory: communication theory[END_REF]. In the present mathematical development,  is defined as the half-thickness of the shell. Indeed, when the added noise   n  is a set of N zero-mean Gaussian PDF independent random variables, the PDF of a realization is the normal distribution   In a practical application where x is unknown and m y is the realization "m", eqs. 4 and 5

are replaced by the estimates

1 1 ˆM M m M m M      y y (6)
and

  2 2 1 1 ˆM M m M M m M        y y (7)
where

1 1 M M m m M    y y (8)
is the estimated value of x calculated from M realizations.

The PDF of a realization is defined by

  2 ˆ, M M    when     : 1, m m M  y and by (9)   2 ˆ, M M M       when     : 1, m m M  y where (10)   1 ˆ1 M M M M       (11

Shell manifold

Let   1 2
, , , , ,

n n l L x f w w w w    ,   1, n N  , (12)
be the N-coordinates system of the manifold X in defined by L variables l w . Note that L=0 for the specific case of the n-sphere. In the general context of an application, the N dependent variables are the response of the process and the L independent variables are the process inputs.

Consider the m th realization 1,

m m m nm nm Nm Nm m m x x x              y xε   1, , , , , , , , ,
of the realization set Y with 1, , ,

, , , ,

m m n m Nm y y y      y   , m  y Y , (14)
and where ,

n m
 is the noise added for the dimension n and m ε is the noise vector in which contributes to the realization m.

Consider the cluster corresponding to repeatable realizations point

m m   y x
ε for a unique point point x of the manifold, the average distance from the realization to the point is

point point m         y x ( 15
)
and corresponds to the radius of the n-sphere eq. 3.

In the case of realizations set Y corresponding to different points on the manifold, considering the special case of a equal amplitude additive Gaussian noise   2 0 0,  for all dimensions and a local small surface patch S having a negligible curvature with respect to the noise vector length with N-L>>1, we can apply some axis rotations to in order to align the local L-surface patch S of the manifold in the bearing of L axis of the new N coordinates system of . The mean distance of the realization normal to the surface

0 0 2 S 0 point 1 n n N L n n N L N L N                             ( 16
)
is shorter than the noise vector since only N-L coordinates contribute to this distance.

Conversely the corresponding variance

    0 0 2 2 2 2 2 S S 0 point 1 1 2 2 n n N L n n N N L N L                                        . ( 17
)
increases when the number of manifold dimensions L increases. This simplest case of a Gaussian noise having a uniform amplitude 0   on the N dimensions and independent of the manifold coordinates yields a hollow cluster of realizations where the PDF of a realization is

    2 0 0 , 2 N L N N L       .
In the general case of the realizations set Y built from dimensions having different noise sources, we conclude that locally

S point    
in the vicinity of a surface patch S. For the general case of M realizations locally available near the surface patch S, we define the average realization-to-manifold distance as

S 1 1 M m m m M       y x  (18)
where m x  is the nearest surface point of the estimated manifold in the vicinity of the surface patch S with respect to the realization m y . This average distance fixes the shell position with respect to the manifold. Since N>L, the shell surrounds the manifold. The estimation of realization dispersion relative to the surface

  2 2 S S 1 1 M m m m M          y x  (19)
is called the hardness of the shell. The S   value is defined as the half-shell thickness.

Note that the average distance between two realizations y i and y j corresponding to the same set of independent variables (same location on the manifold) is S 0

2 2 i j N E N N L              y y (20)
for the simplest case of a Gaussian noise having a uniform amplitude 0  on the N dimensions and independent of the manifold coordinates. The nearest neighbor distance of a realization is on average less than the previous average distance (eq. 20) and diminishes with the increasing measurement population. When the noise is a function of the manifold coordinates, the shell position and hardness also vary with the manifold coordinates. In the more general case of a diffuse manifold (Fig. 1d), which is not defined by eq. 12, the PDF of the realizations corresponds to the multidimensional convolution of the manifold density with the  

2 S S ,      function.
The cluster appears hollow only if the total noise amplitude is significantly greater than the manifold width.

In the case of multiple internal states, the resulting PDFs of the realization can be superimposed. To sum up, Fig. 2 shows an artistic 3D illustration of a hollow cluster and the PDF of realizations when the manifold is known (b) and when it is estimated (c).

When the manifold is estimated, the quadratic form  

A sample application with a generic comparison algorithm

The algorithm compares two vector quantities X and Y in . The vector may be a realization or an estimate (bin average or manifold interpolation). The average of the distance ||X-Y|| is assumed constant until the repeatable response of the process is no longer there. In the following algorithm sequence, the "match" threshold level is set at four times the standard deviation These errors correspond to the estimation error vector length in . If X or Y is a realization, no error is estimated for it. For the case of a clustering followed by an interpolation (Fig. 3), the Ex or Ey errors correspond to the manifold dispersion  M . This dispersion can be estimated from error propagation in the clustering process taking account of the weight of each cluster in the interpolation and the distance between an interpolated point on the manifold with the surrounding clusters. More computationally expensive, the manifold dispersion estimation can be obtained by bootstrapping [START_REF] Efron | Bootstrap Methods: Another Look at the Jackknife[END_REF]. Whatever the method used, the computational complexity deployed in these error estimations substantially exceeds that needed for process lines 6-11.

Lines 9 and 10 allow the detection of an abrupt distance change between X and Y. To detect a slow trend, the shelldist and shellvar values must be frozen when the population m is large enough to be representative of the average initial response. As a coding option, the shelldist and shellvar values can be estimated considering the new realization only when match is true.

Practical monitoring applications

A monitoring application software includes measuring, database management, numerical signal processing, estimation, comparison and alarm management. Signal processing transforms the raw measurements into a distilled form in a representation space (envelope, correlogram, cyclogram, Fourier, time-frequency, scale-frequency, high-order spectrum…) where the signal-to-noise ratio appears at its maximum in the context of the application task. The estimation builds the initial history of a process and refreshes the actual process average response. With the "match" logic case replaced by a "no-alarm" logic case (lines 4 and 8), the previous generic comparison algorithm can be used for monitoring. Usually two comparisons are done, one between the initial history and the actualized average response for slow trend detection, and the second between the actualized average response and the last realization for the detection of a fast growing defect [START_REF] Léonard | Trending and pattern recognition using the spherehardening phenomenon of signal multidimensional projection[END_REF].

In vibro-acoustic monitoring of commutating equipment, the signal envelope is compared to a reference envelope built from data history [START_REF] Léonard | Trending and pattern recognition using the spherehardening phenomenon of signal multidimensional projection[END_REF]. The equipment response is found in the N time samples of the vibro-acoustic envelope [START_REF] Léonard | Vibro-acoustic signature treatment process in high-voltage electromechanical switching system[END_REF], typically 1000 samples. The n-sphere method described in eqs. 1-11 has been applied with limited success since the equipment response is significantly a function of the temperature. The binning of the reference envelope as a function of a set of temperature ranges was a first step toward increasing the accuracy of the method. The trajectory manifold (L=1) is modeled by a set of bins, each bin containing an average location in , an average temperature and the corresponding population. In a second step, we added a krigging algorithm [START_REF] Krige | A review of the development of geostatistics in South Africa. Advanced Geostatistics in the Mining Industry[END_REF][START_REF] Matheron | Les concepts de base et l'évolution de la géostatistique minière. Advanced Geostatistics in the Mining Industry[END_REF][START_REF] Journel | Mining Geostatistics[END_REF] to interpolate the reference value from the bins. The last realization is then compared to the reference interpolated at the same temperature as the realization. The bin values, populations, temperatures and dispersions are used as krigging inputs. Each of the 1000 time samples is interpolated. In monitoring a hydroelectric generating unit [START_REF] Léonard | Application de la statistique et du krigeage à la surveillance de groupes hydroélectriques à Hydro-Québec[END_REF][START_REF] Léonard | Hydro Equipment: Learning from Hydro-Québec's Monitoring Experience[END_REF], the operating condition is determined by many independent variables: gate opening, water levels, rotor excitation current, cooling water and air temperatures,… If we include vibrational data (peak-peak, RMS, mean and spectral line values) as part of the dependent variables, the total number of measured variables reaches over 100. In the proposed method, a dynamic clustering estimates the manifold corresponding to the behavior of the unit. The clustering in is a dynamic clustering inspired from the one described in [START_REF] Léonard | Dynamic clustering of transient signals, patent applications WO2012162825 and US20140100821 A[END_REF][START_REF] Bellemare | Partial discharge sniffer to ensure worker safety in underground vaults[END_REF] and described in the addendum. Since the measured variables are expressed in different engineering units, we added a dynamic normalization of variables as a function of the maximum observed range of each variable. The clustering is usually applied using the information from the L independent variables (see mask = 1 in the addendum). The clustering can also use the information about some dependent variables in order to distinguish between a limited number of internal states (Fig. 1b) but this feature has not been tested. The average value of each dependent variable and the average square value are calculated for each cluster.

When enough clusters with enough populations are found from the data history, krigging is used to interpolate the N dependent variables averaged in clusters. Krigging uses the cluster values, dispersions and populations to estimate the interpolated (or extrapolated) value and its accuracy. The interpolated value and its accuracy  M are both used in the comparison with the last realization.

For both applications, the accuracy of the interpolation is added quadratically (Fig. 3, lower) to the shell realization dispersion to fix the alarm bound (typically ±4 ). Before the comparison with the alarm bound, the mean distance S   is subtracted from the distance between the interpolation and the last realization. Note that every change in the equipment response pushes the resulting distance to increase toward the positive bound (+4 ). The newly developed method dramatically reduces the rate of false alarms (type 1 error) for a same missing-defect rate (type 2 error).

Clustering (or binning) is an unavoidable part of the method. "Clustering should be thought of as the statistical counterpart to the geometric construction of the pathconnected components of a space, which is the fundamental building block upon which algebraic topology is based" [START_REF] Coray | Clustering algorithms with prototype selection[END_REF]. In the proposed applications, the clustering is done function of the L independent variables assuming a single continuous and monotonous manifold. Moreover, the measurement noise of independent variables is not considered: this noise is implicitly transferred to dependent variables.

Conclusion

For a repeatable response process, data clusters built from noisy data are hollow when the noise amplitude is greater than the dispersion generated by the hidden internal states of the process. When the clusters are built from repeatable realizations, clustering and interpolation can be used to estimate the manifold corresponding to denoised information. The probability density function of realizations generates a shell enveloping the manifold and the corresponding model can be estimated from realizations. Big data subsets related to measurements are shell manifolds.

The proposed vision suggests a new way to process information about noise in multidimensional statistics. The PDF contribution of a single measurement channel in the context of a multichannel measurement has been explained and exploited. The likelihood of a realization appears as function of the realization-to-shell distance rather than the realization-to-manifold distance.

Many of the current tools of multivariate analysis are affected by this vison since they do not target hollow clusters or shell clusters. For each of these tools, the question is: can we find an algorithm improvement to adapt the tool to noisy data clusters? This new proposed vision is a call to everyone looking for new tools adapted to noisy data clusters including large multidimensional datasets.

Addendum -Clustering algorithm

The normalization of input variables having different engineering units is not addressed in the following description of a clustering algorithm. Dynamic clustering is a compromise between optimality and computational complexity. In the proposed clustering algorithm, the final number of clusters kmax is predetermined and fixes the memory dimension. The cluster centroid locations [START_REF] Zheng | Dynamic clustering analysis of documents based on cluster centroids[END_REF] are adjusted dynamically at each step considering the last clustering result and the new realization. The decision rule is based on a Euclidian distance calculated from a subset of available dimensions. In a realization, L measurements correspond to independent variables and N others correspond to dependent variables. The function  

1 : independent variable 0 : dependent variable mask i     with   1, i N L   (A-1)
defines which subset is used in the distance calculation

      1/ 2 2 1 ( ) N L i i masked i k m a s ki C k X              C X . (A-2)
The average realization-to-cluster centroid variance cvar(k) is calculated for each cluster based on a masked distance calculation. The square root of this variance is averaged in order to share a unique shell distance value, shelldist, for all clusters. The shelldist value adjusts the cluster-to-cluster sampling distance on the manifold dynamically. Note that the increasing kmax value reduces the sampling distance and, also, the average cluster population as well as the accuracy of the centroid estimate. The number of clusters kmax fixes the tradeoff between the sampling distance and the cluster population. Clustering is successful when the clusters are uniformly distributed on the manifold. For the same realization set, different sampling orders of realizations yield similar results.

In order to reduce the computational expense, the maximum allowed distance dmax is defined as cdist shelldist  and corresponds to the decision level to merge a realization to the nearest cluster without the computational task of considering a possible fusion of two clusters. One cluster fusion 1 k k   decision implies the use of a k-Nearest Neighbors (k-NN) algorithm [START_REF] Fix | Discriminatory analysis, nonparametric discrimination[END_REF][START_REF] Cover | Nearest Neighbor Pattern Classification[END_REF]

realizing   1 2 k k  
cluster-to-cluster comparisons without optimization. cdist is set typically at about 1 to 3, a function of L and the desired accuracy. Sharing the same dmax value for all clusters forces the cluster variances to converge into a similar range which corresponds to a uniform sampling of the explored part of the manifold. The algorithm stores the first kmax realizations as new clusters and starts the subsequent clustering. Considering X as a new realization, C(k) as a matrix of the kmax cluster centroids, V(k) as a matrix of cluster dispersion, p(k) as a vector of cluster populations and shelldist as the average shell cluster distance, the algorithm steps are: At the beginning of the clustering process, the shelldist is underestimated (line 9) since every new cluster has zero variance. An underestimated value of dmax increases the probability of cluster inter-distance matrix calculation, increasing the accuracy and also the computational expense. When a second member X j is merged with a first member X k of a new cluster, the variance   As demonstrated in monitoring a hydroelectric generating unit in Fig. 4, the algorithm converges to minimize calls to cluster-to-cluster matrix calculation. Note the "0" count during initialization and the asymptotic pattern of the curve which converges to add new calls only when the sampling touches an unexplored part of the manifold.

            1 k p k k p k     C X C ;                     1 k k k p k k p k       V C X C X V  ; '  : Hadamard product'            

Figure 1 :

 1 Figure 1: 2D illustration of the response function, a dependent variable (DV) function of an independent variable (IV), for a) a single internal state, b) several internal states, c) a hysteresis internal state and d) an additive effect of numerous internal states.

.

  However, when the added noise   n  is a set of N zero-mean random variables uncorrelated, the central limit theorem stipulates that the PDF tends also toward the normal distribution   number of dimensions.

1 M y and 2 My

 12 ) is an additional distance explained by a new realization not taken into account in the M y estimation. The case of a realization as an average number of realizations and that of the comparison of two estimations calculated for different sets of realizations are discussed in[START_REF] Léonard | Quantitative analysis of signal-related measurements for trending and pattern recognition[END_REF].

  to estimate the likelihood of a realization where  M is the manifold dispersion estimated in the vicinity of the realization.

Figure 2 :

 2 Figure 2: Cutaway view of a cluster showing a manifold (left). Realizations on the outside shell are in black, on the inside in grey and in the cut view in white. On the right of the figure: a) a cutaway view through the manifold surface, b) the corresponding PDF for a known manifold and c) the general case of an unknown manifold estimated from realizations averaging.

Figure 3 :

 3 Figure 3: Sequence of operations in estimating the likelihood of the last realization. Realizations are clustered in order to estimate the inner manifold. An interpolation of this manifold is then compared to the last realization. In the upper right illustration, the cluster disk is proportional in size to the cluster population. Clustering estimates the values S   and

  the Smallest Distance Realization-to-Cluster and the corresponding cluster k' 8. If sdrc > dmax Then 'Calculates the cluster-to-cluster distances triangle matrix and finds dscc, the smallest distance, and the corresponding cluster couple (j,k)'; If sdrc > sdcc Then 'merge the two corresponding clusters' Merges the realization to the nearest cluster'

Figure 4 :

 4 Figure 4: Counts of cluster-tocluster matrix calculation as a function of the number of processed realizations. In this example, N=14, L=10, cdist=1.5, kmax=200, three years monitoring.

  at the initialization step. The contributing variances Ex and Ey of the estimates X and Y can be included in the line 8 as:

	where	alpha 	  0,1	is the smoothing factor of the moving average with a weighting
	factor	alphaweight alpha 	 1  	 alpha alphaweigh t 	calculated before line 9 and with
	0 d shelldist alphaweight  8. If 	4 (  	2 shellvar Ex Ey  	2 1 / 2 )	then match=. false. ; else match=.
		true.
											1/2 shellvar :
	1. m =0; 'Number of processed comparisons'
	2. shelldist =0 : 'Mean manifold-to-shell distance'
	3. shellvar =0 ; 'Variance of calculated distances to the mean distance'
	4. match = .true. ;
	5. Wait, calculate or find in memory the X and Y vectors;
	6. d=||X-Y|| ; 'Calculate the distance'
	7. If m=0 Then Goto 10 ; Endif ; If m=1 Then Goto 9 ; Endif;
	8. If	d shelldist 	4  	shellvar	1/ 2	then match = . false. ; else match = . true. ; Endif ;
	9.	shellvar		((	d shelldist 	)	2	     1 m	) / shellvar m	;
	10.	 shelldist d m shelldist   	  /	m	 1  ;
	11. m =m+1 ;
	12. goto 5 ;
	Lines 8 to 10 are the kernel of the application algorithm and contain the essential part of
	the innovation. Lines 9 and 10 can be replaced by moving averages as:
	9.	( shellvar alpha d shelldist (   	2 ) ( 1  	) alpha shellvar alphaweight ) 
	10.	 1    shelldist alpha d  	 alpha shelldist alphaweight  

Acknowledgements

I would like to thank Dr. Tomasz Kaczynski, professor at Sherbrooke university and Dr. Odile Marcotte, Deputy Director at the Centre de Recherches Mathématiques, for their comments and suggestions during the writing of this article. I should also acknowledge the support provided by my administration which gave me all its trust and the budgets needed to develop and promote this new statistical vision. Mention must also be made of my colleagues who are responsible for implementing the algorithms in industrial systems, in particular engineer-researcher Michel Gauvin, who was the first to implement these algorithms in a monitoring system which is currently operational on our equipment.


A Cozza

F Monsef

Layered Electric-Current Approximations of Cylindrical Sources

Keywords: Cylindrical sources, equivalent currents, modal theory, antennas

Formulations of the equivalence theorem only involving electric currents either imply non-null inner fields, or the inclusion of a bulk perfect magnetic conductor. The possibility of maintaining a homogeneous background medium across the equivalent-current surface and null inner fields without recurring to magnetic currents is here investigated for the case of cylindrical sources. The proposed procedure keeps using the definition of equivalent electric currents as found in A. Love's theorem, introducing further layers of auxiliary electric currents to mimic the role of the missing magnetic ones. Optimal scalar coefficients applied as weights on each layer current are found, so to minimize the mean quadratic error in the fields radiated outside the auxiliary shell region, while at the same time minimizing internal radiation. The procedure is inherently well-conditioned, thanks to the reduction in the number of degrees of freedom involved, equal to the number of layers, as forced by a regression approach. A thorough numerical analysis is presented, where residual errors are found to be in the range of a few percent points.

Introduction

Equivalence theorems are a workhorse of antenna theory, widely applied in order to simplify the analysis and the computation of fields generated by a source of radiation. Initially derived by Love and MacDonald more than a century ago [START_REF] Love | The integration of the equations of propagation of electric waves[END_REF][START_REF] Macdonald | The integration of the equations of propagation of electric waves[END_REF], they can be flexibly declined into versions that can suite a large number of practical configurations [START_REF] Smythe | The double current sheet in diffraction[END_REF][START_REF] Harrington | Time-Harmonic Electromagnetic Fields[END_REF][START_REF] Rengarajan | The field equivalence principle: Illustration of the establishment of the non-intuitive null fields[END_REF]. Independently from the specific version, they all share the same basic idea, illustrated in Fig. 1: given an original source Γ, identify a spatial distribution of equivalent currents over a closed surface Σ containing Γ, such that these currents generate a perfect reproduction of the field distribution expected for Γ over Ω e , outside Σ. Analogue versions exist for a source Γ outside Σ, a case not treated in this paper .

Referring to the nowadays standard formulation due to Schelkunoff [START_REF] Schelkunoff | Some equivalence theorems of electromagnetics and their application to radiation problems[END_REF], equivalent electric, J eq (r ), and magnetic, M eq (r ), currents are required, if the equivalent problem in Fig. 1(b) is to work in a homogeneous medium across the surface Σ. In this case, it can be shown that [START_REF] Love | The integration of the equations of propagation of electric waves[END_REF] J eq (r ) = δ(rr Σ ) n(r Σ ) × H (r )

(1a) M eq (r ) = δ(r -r Σ )E (r ) × n(r Σ ), (1b)
where E (r ) and H (r ) are, respectively, the electric and magnetic field generated by the original source at the position r ; δ(•) is Dirac's delta distribution, here used in order to represent a singular current sheet limited to Σ; n(r Σ ) is the outward unit vector normal to Σ, at the position r Σ . A time convention exp(jωt ) is assumed, but only phasors are shown. All quantities are frequency dependent, but the pulsation ω is dropped for the sake of brevity. In case the fields radiated by a source were obtained through experiments, both electric and magnetic fields would need to be accessed; in practice, only probes sensitive to one kind of field could be available. Equivalent electric currents suffice if Ω i is filled with an ideal perfect magnetic conductor (see Fig. 1(c)), as described in [START_REF] Schelkunoff | Some equivalence theorems of electromagnetics and their application to radiation problems[END_REF], ensuring a perfect field reproduction over Ω e . Unfortunately, reproducing this kind of scenario in a numerical model, a bulk material in Ω i would behave as a scat-

W i R 1 R L Figure 2:
The paradigm considered in this paper: L multiple-layer auxiliary electric-current distributions are found in a shell shown in shade, each occupying a surface Σ i . Internal and external test surfaces Σ i T and Σ e T , as introduced in sec. [START_REF] Harrington | Time-Harmonic Electromagnetic Fields[END_REF], are used for the optimization of the layer weights. The new inner and outer regions are referred to as primed quantities, to differentiate them from those defined in standard equivalent-theorem formulations (see Fig. 1).

terer in case waves impinged onto it. Formulations using only electric currents while keeping a homogeneous medium across Σ would be required in this case, as discussed in [START_REF] Martini | An equivalence theorem based on the use of electric currents radiating in free space[END_REF].

Solutions for this problem can be defined, as shown in [START_REF] Harrington | Time-Harmonic Electromagnetic Fields[END_REF][START_REF] Las-Heras | Radial field retrieval in spherical scanning for current reconstruction and nf-ff transformation, Antennas and Propagation[END_REF][START_REF] Colinas | Application of the Tikhonov regularization technique to the equivalent magnetic currents near-field technique, Antennas and Propagation[END_REF][START_REF] Barriere | Mapping of equivalent currents on high-speed digital printed circuit boards based on near-field measurements, Electromagnetic Compatibility[END_REF]11], but imply non-null fields in Ω i . This can be a problem since inward radiation inevitably leads to outward waves after propagating through Ω i . In this paper we are interested in equivalent currents that radiate only outside Σ, as a continuation of the fields radiated by the original source Γ; therefore, null fields are required in Ω i .

To this end, magnetic equivalent currents could be converted into an additional electric contribution J M (r ) = ∇×M eq (r )/jωµ [12]. But in case of an experimental characterization of the original source, this operation is ill-posed, since it involves estimating derivatives from measurements. Critically, normal derivatives would require sampling fields at different distances from the original source. Moreover, this further electric-current contribution now involves spatial derivatives of a Dirac distribution, which does not correspond to a double-sheet of currents.

In order to overcome these difficulties, this paper suggests using multiple layers of electric currents, defining a shell region Ω ML as shown in Fig. 2. The rationale is to make the fields radiated by the different layers of currents interfere, in such a way as to emulate the interplay between fields radiated by electric and magnetic currents in Love's theorem.

Our proposal is not to set up an inverse problem where the equivalent current distributions over each layer would act as unknowns, since this kind of approach is often ill-conditioned, particularly in the case of measured data. We rather attempt extending the use of distributions defined as in (1a) over at least two surfaces, by defining optimal scalar coefficients acting as weights on each of these current distributions. As it will be shown in the rest of the paper, the proposed method, though effective, does not allow a perfect equivalence, hence our referring to these currents as auxiliary rather than equivalent.

Scalar representation of the problem

In the following we will consider the case of original sources described by means of electric-current distributions confined within a shell region of space Ω ML (Fig. 2). For ease of derivation, a cylindrical geometry will be assumed, with translation invariance to axial displacements.

Dealing with cylindrical sources, it is convenient to proceed by separating the original currents modelling the source over the region Γ (Fig. 1) into axial and transversal components

J s (r ) = J s (r ) + J ⊥ s (r ) = ẑ J s (r ) + φJ ⊥ s (r ), (2)
where bold variables represent vector quantities while hatted ones are unit vectors. In order to simplify our derivation, we introduce magnetic currents [12, section 7.12]

M s (r ) = ẑ M s (r ) = - 1 jωǫ ∇ × J ⊥ s (r ). (3)
In the following, instead of the original current distribution (2) we will consider a source described by axial currents J s (r ) and

M s (r ).
Maxwell's equations for cylindrical configurations can be expressed as [12, section 14.1]

∇ 2 + ω 2 ǫµ 2 E (r ) = -jωµJ s (r ) -∇ • (ẑ × M ⊥ s (r )) (4)
∇ 2 + ω 2 ǫµ 2 H (r ) = -jωǫM s (r ) -∇ • (ẑ × J ⊥ s (r )), (5)
hence describing the source through only axial electric and magnetic currents yields scalar decoupled equations

∇ 2 + ω 2 ǫµ 2 E (r ) = -jωµJ s (r ) (6)
∇ 2 + ω 2 ǫµ 2 H (r ) = -jωǫM s (r ), (7)
corresponding to a TM (transverse-magnetic) and TE (transverse-electric) problem, respectively. Hence, ( 6) is solved by

E (r ) = ωµ 4 Γ dr ′ H 0 (k 0 r -r ′ )J s (r ′ ), (8)
where H 0 (•) is Hankel's function of zero-th order of the first kind (i.e., for an outgoing wave). Equation ( 8) is also a solution to [START_REF] Martini | An equivalence theorem based on the use of electric currents radiating in free space[END_REF], by substituting electric-field related quantities with magnetic ones, as according to the duality principle; for this reason, the rest of the paper will only focus on the case of axial electric sources. Solutions found for this problem can be readily applied to the case of azimuthal electric sources thanks to (3).

Fields radiated by a single layer of electric currents

We first seek a modal expansion of the electromagnetic field radiated by the original source in Γ, by applying the addition theorem to [START_REF] Las-Heras | Radial field retrieval in spherical scanning for current reconstruction and nf-ff transformation, Antennas and Propagation[END_REF], here recalled for the sake of completeness, for the special case of the outgoing Hankel function of the zeroth order

H 0 (k o r -r ′ ) =        ∞ m=-∞ J m (k o r ′ )H m (k o r )e jm(ϕ-ϕ ′ ) r ′ < r ∞ m=-∞ H m (k o r ′ )J m (k o r )e jm(ϕ-ϕ ′ ) r ′ > r (9)
where J m (•) and H m (•) are Bessel's and Hankel's functions of the mth order of the first kind. From ( 8) and ( 9) we straightforwardly obtain

E (r ) = ωµ 4 ∞ m=-∞ γ m ψ h m (r ) r ∈ Ω e , (10)
having introduced a shorthand representation of cylindrical modes

ψ h m (r ) = H m (k o r )e jmϕ (11a) ψ j m (r ) = J m (k o r )e jmϕ (11b)
and the associated modal weights

γ m = Γ dr ′ ψ j m (r ′ ) * J e (r ′ ), (12)
with * the complex conjugate.

The magnetic field can also be represented in a similar manner; making use of the electric-field curl equation,

H ⊥ (r ) = j 4 ∞ m=-∞ γ m ∇ × ẑψ h m (r ) = = ∞ m=-∞ γ m r jm r -φ∂ r ψ h m (r ), (13)
where ∂ r is the first-order derivative taken with respect to the radial distance r . In order to understand what is the effect of removing equivalent magnetic currents in a free-space configuration, we consider a cylindrical surface Σ, of radius R, over which

J eq (r ) = -δ(r -R) ∞ m=-∞ γ m e jmϕ ∂ r H m (k o r ), (14)
as given by (1) and (13).

The electric field generated by these currents can be computed by applying [START_REF] Las-Heras | Radial field retrieval in spherical scanning for current reconstruction and nf-ff transformation, Antennas and Propagation[END_REF] to the current distribution (14). The application of the addition theorem (9) yields

Ê (r ) = ωµ 4 ∞ m=-∞ γ m η e m (R)ψ h m (r ) r ∈ Ω e η i m (R)ψ j m (r ) r ∈ Ω i (15) with η e m (R) = -jπR J m (k o R)∂ R H m (k o R) (16a) η i m (R) = -jπRH m (k o R)∂ R H m (k o R). (16b)
The corresponding magnetic field is

Ĥ⊥ (r ) = j 4 ∞ m=-∞ γ m r jm r -φ∂ r η e m (R)ψ h m (r ) r ∈ Ω e η i m (R)ψ j m (r ) r ∈ Ω i . ( 17
)
Fig. 3 shows an example of the electric and magnetic field distributions for the ideal case expected for the equivalence theorem, and those observed in the case where only equivalent electric currents are present, i.e., a partial application of the equivalence theorem. These results refer to the test-case described in section 8, dealing with a directive linear source four wavelengths long. This source has a directive leftward radiation.

In the case of partial application of the equivalence theorem, the leftward radiation is well reproduced, whereas a spurious rightward radiation appears. This last is due to the fact that equivalent electric currents that are at the origin of the leftward radiation do also radiate towards the inner region Ω i ; this inward radiation would be counterbalanced by equivalent magnetic currents in the equivalence-theorem free-space formulation. Left alone, the inward contribution may focus inside the inner region Ω i and subsequently diverge rightwards, causing the spurious radiation that did not appear in the original configuration. Hence the need for an alternative mechanism of control of the inward radiation.

Multiple-layer auxiliary currents

Reproducing the results of Love's formulation would require, apart for a common constant factor, to ensure that η e m (R) = 1, ∀m and η i m (R) = 0, ∀m. In other words, the η i,e m (R) coefficients act as distortions, in the modal expansion of the field distribution generated by the original electric currents J e (r ), by passing, e.g., for the region external to Σ, from modal coefficients {γ m } to {γ m η e m (R)}.

Our proposal is to introduce multiple layers of auxiliary current distributions, defined as in Fig. 2 over L concentric cylindrical surfaces {Σ l }, with respective radii {R l }. The rational behind this procedure is that by applying different weights to the individual electric current distributions, the resulting overall modal coefficients can be forced to match more closely those expected for [START_REF] Love | The integration of the equations of propagation of electric waves[END_REF]. The electric currents over each layer is defined as in Love's equivalence theorem, but are further weighted by coefficients {A l }, leading to overall fields (15) and (17) with a new set of coefficients

η i,e m = L l =1 A l η i,e m (R l ) (18)
As discussed in the next section, the problem is now to identify the optimal weights {A l } that reproduce as closely as possible the electromagnetic field distribution generated by the original source over Ω ′ e , and null fields over Ω ′ i , while the region of space Ω ML occupied by the layers will act as a buffer region of finite thickness.

Residual errors are minimized over test surfaces Σ e T and Σ i T , respectively for external and internal fields. Clearly, the equivalence theorem implies that any choice of these surfaces is equivalent. Therefore, in the following, we assess the accuracy of the approximate field over cylindrical test surfaces Σ e T and Σ i T of radii R e T > R L and R i T < R 1 , respectively, as illustrated in Fig. 2. Only errors on tangential fields will be evaluated. Even though the proposed procedure can be extended to any number of layers, we will hereafter set L = 2, in order to look for solutions offered by the simplest possible case.

The rms errors are here defined as

e 2 E (R T ) = (ωµ/4) -2 Σ T dr Ê (r ) -E (r ) 2 (19) e 2 H (R T ) = (1/4) -2 Σ T dr φ • Ĥ⊥ (r ) -H ⊥ (r ) 2 , (20)
perfect reproduction of the original field distributions over Ω ′ e would now imply the condition η e m = 1, ∀m. Hence, from ( 10), (13), ( 15) and ( 17)

e 2 E ,H (R T ) =        ∞ m=-∞ |η e m -1| 2 |β E ,H m (R T )| 2 R T > R L ∞ m=-∞ |η i m | 2 |β E ,H m (R T )| 2 R T < R 1 (21)
where

β E m (R T ) = γ m H m (k o R T ) R T > R L γ m J m (k o R T ) R T < R 1 (22)
and

β H m (R T ) = γ m ∂ r H m (k o R T ) R T > R L γ m ∂ r J m (k o R T ) R T < R 1 (23)
are proportional to the individual modal contributions to the original field distributions required by Love's formulation.

Optimal layer weights

The errors in (21) are ultimately functions of the layer weights {A l } and should be minimized in order to find an optimal excitation of the auxiliary-current layers. In this respect, all summations will be limited to terms |m| < M, where M can be chosen from one of the heuristic criteria routinely applied in electromagnetic theory [13,14].

The simplest option for minimizing the rms errors consists in observing that e E can be interpreted as the length of the difference vector found between two vectors whose scalar elements are respectively {β E m η e m } and {β E m }. In this respect, 21) can be minimized by means of a leastsquare (LS) approach, by introducing the normal equation [15]

e 2 E (R T ), R T > R L in (
β E (R e T )η e A = β E (R e T )c 1 , (24)
where

β E = diag{β E -M , . . . , β E M }, η e is an 2M + 1 × L ma- trix whose columns are composed by the {η e m (R l )}, A = [A 1 , . . . , A L ]
T are the L layer weights to be optimized and c a is a 2M +1 column vector of elements equal to a, here a = 1. In practice, (24) translates the problem of finding the optimal

{A l } such that {η e m } ≃ 1, ∀ m ∈ [-M, M],
A = β E (R e T )η e † β E (R e T )c 1 (25)
with † the Penrose-Moore inverse. In general, this solution will depend on the set of LS-weights β and therefore on whether one intends to optimize the reproduction of the electric or of the magnetic field. In practice, the solution is weakly dependent on which of the two sets (22) or (23) are chosen, since it can be shown that

β H m (r ) β E m (r ) = H ′ m (k o r ) H m (k o r ) ≃ j ∀r W, (26)
with W the maximum transversal dimensions of the source. Hence, the optimization weights {β E ,H m } are practically proportional one to the other for a given m. As a consequence, both electric and magnetic fields will be asymptotically optimized at the same time. This trend is indeed confirmed by the results presented in sec. 8 for external test surfaces, while it is not necessarily the case for inner test surfaces, where choosing between the two sets of {β m } can have a measurable effect on the inner error, as shown in sec. 8.4.

Internal & external optimization

The above results are divided into internal and external regions defined by the auxiliary-source surfaces. It is therefore natural to wonder if layer weights could be optimized for both regions at the same time, or if any joint optimization can only be obtained by degrading the best performance attainable when dealing only with one region.

In order to gain some insight into this issue, we recall that for

k o r ≫ 1 H m (k o r ) ≃ 2 πk o r e j(k o r -mπ/2-π/4) , (27)
while

H ′ m (k o r ) ≃ jH m (k o r ). Since J m (k o r ) = ReH m (k o r ), in case k o R l ≫ 1, ∀l ∈ [1, L], (16) is asymptotically approxi- mated by η e m (R) ≃ 1 -j(-1) m e 2jk o R (28a) η i m (R) ≃ -2j(-1) m e 2jk o R , (28b)
for a single layer of auxiliary electric currents. When applying these results to (18), for L layers,

η e m ≃ L l =1 A l -j(-1) m L l =1 A l e 2jk o R l (29a) η i m ≃ -2j(-1) m L l =1
A l e 2jk o R l .

(29b)

A close look at (28) shows that the {η e m } appear to be given by an offset version of the {η i m }, by a fixed quantity not depending on the modal index m. As a result, when the condition of null internal field η i m = 0, ∀m, is fulfilled, η e m = 1, ∀m, hence implying a perfect reproduction of the field topography over the external region. Therefore, internal and external optimization can go hand in hand as the above condition is sufficient to ensure an exact reproduction of the field topographies dictated by Love's formulation of the equivalence theorem, at least within an asymptotic framework. Moreover, (29) also explains how the L layer weights can control M ≫ L modal distortions, thanks to the asymptotic factorization of the portions impacted by the l and m indexes.

In practice, the minimization of the metrics in (21) leads to sub-optimal solutions, resulting in inevitably m-dependent modal distortions {η i,e m } . In fact, it should be recalled that the LS problem discussed in sec. 5 is typically overdetermined, since L ≪ M in most configurations: in the present paper we will consider L = 2, while M 1 even for weakly directive sources. Under these conditions, the only possibility for acceptable solutions must pass through the minimization of the portion shared by the external and internal distortions, as seen in ( 29). In any other case, the number of degrees of freedom would be too small with respect to the modal weights to control them independently.

More generally, a joint optimization over the internal and external surfaces can be considered, by extending (24) as

β E (R e T )η e β E (R i T )η i A = β E (R e T )c 1 c 0 , (30)
which corresponds to a Tikhonov regularization term in the original normal equation (24). Indeed, the squared-norm of the residual error between the left-and the right-hand side terms in (30) reads as

β E (R e T )η e A -β E (R e T )c 1 2 + β E (R i T )η i A 2 , (31)
whence a Tikhonov matrix Γ = β E (R i T )η i , i.e., Γ ∈ C 2M+1,L . The definition of Γ implies that the regularization is not operating in a isotropic way along the 2M + 1 modal dimensions, and is rather controlled by the modal coefficients expected for the radiation over the internal test surface Σ i T .

Approximate solutions

As opposed to the above procedure, alternative approximate approaches can be defined for the computation of layer weights. Two such solutions are here considered, one obtained from the asymptotic development of Hankel's functions, while the second is based on a local reasoning looking at current layers as a collection of directive sources.

The simplicity of these approximate solutions makes them attractive alternatives to the optimal solution presented in the previous section. A direct comparison of their respective accuracy is presented in sec. 8.

Asymptotic solution

Starting from the asymptotic expressions of η i,e m (R) in ( 29), the problem defined by setting the constraints η i m = 0 and η e m =

1 can be solved in closed form. For the special case of L = 2, the result is given by

A 1 = -A 2 e -j2k o d (32a) A 2 = 1 1 -e -j2k o d . (32b)
Since (32) were obtained by requiring minimal modal distortions, the above procedure does not aim at reproducing the field distribution on a global scale as (25), but rather attempts to optimize the modal weights on an individual basis. As shown in sec. 8, while this procedure provides results similar to those obtained with the optimal procedure over the inner region Ω ′ i , its performance in the outer region is definitely worse.

End-fire approach

An alternative approximate solution can be proposed by looking at Love's formulation from a different viewpoint. Consider the electric and magnetic equivalent currents it defines at a generic position r ∈ Σ; these two elementary sources make up the smallest directive source, sometimes referred to as a Huygens' source. Inward radiation along -r can be demonstrated to be nil.

This observation provides a clue as to how multiple layers of electric currents can emulate the results of the equivalence theorem, by forcing the radiation from a second layer to interfere destructively with the first one along inward directions, while summing up in phase for outward directions. Such a configuration can be interpreted as a (continuous) collection of end-fire arrays, oriented radius-wise; within this framework, the layer weights A 1 and A 2 are which are similar to those found in (32). The main difference, apart the way the results were derived, is that in the case of the end-fire approach the excitation of the two radiating elements is the same apart for a phase shift, as usually done for this kind of arrays, whereas in the asymptotic case we still compute two different auxiliary current distributions. The possibility of generating the information for the second layer from only a single one can be an asset when these data are measured rather than computed, e.g., from near-field facilities, as it would cut the measurement time by half. The end-fire approach falls short of the global approach employed in the optimal procedure as it does not aim to control the modal distortions. As shown in sec. 8, the residual error increases not only in the external region, but also in the internal one, by giving raise to a spurious focusing roughly covering the space occupied by the original source. Still, its overall ability in controlling the modal distortions is not much worse than that displayed by the asymptotic approximation.

A 1 = -A 2 e -jk o d (33a) A 2 = 1 1 -e -j2k o d , (33b)

Numerical results

The accuracy and effectiveness of the proposed approaches are assessed for the case of directive sources described by a linear current distribution J s (r ) of length W s

J s (r ) = ẑe -jβx sin(πx/W s )δ(y)δ(z) |x| ≤ W s /2, (34)
tapered toward its ends by a sinusoidal profile in order to produce a dominant main lobe. Fig. 4 shows the electric field distribution they produce for W s = 2λ and W s = 4λ. The amplitudes of their modal coefficients are shown in Fig. 5. These configurations will be used throughout this section as references in evaluating the performance of layer coefficients derived by means of the optimal procedure (25), the asymptotic approximation (32) or the end-fire approximation (33).

All of the following results related to the optimal procedure involve the {β E m } as LS weights. Differences observed when using the {β H m } are discussed in 8.4.

Modal distortions

The question of how accurately two layers of electric currents can emulate these configurations can be assessed by means of several indicators. The first here considered, perhaps the most straightforward, is the set of modal distortions {η m } for the external and internal modal coefficients ( 16), which account for the imperfect reproduction of the original modal coefficients. Some examples of results for the external distortions are given in Fig. 6, for three values of the inner radius R 1 ; the case R 1 = W s /2 corresponds to the boundary of the minimum cylinder containing the original source. Unless otherwise mentioned in the rest of the paper, the second layer is always found at a distance d = λ/4 outside the first one.

Fig. 6 shows that the end-fire approximation, though yielding a limited error over the modal indexes m where most of the radiated energy is found (see Fig. 5), actually presents a rapidly diverging error for |m| > 10. More fundamentally, η e m for the end-fire case present a systematic error (bias) that will be shown in sec. 8.3 to have a major impact on the field topographies.

On the contrary, the optimal procedure has a very good performance for low-order modal indexes, where the {β m } coefficients enforce the most stringent conditions for the LS algorithm. It can be noticed how this performance comes to the expenses of a distortion that peaks for |m| = 12, though the distortions are always bounded in these examples and have a minor impact on the overall results, since they apply to modes that are very weakly excited.

Increasing R 1 external distortions flatten out over a wider range of modes for the three approaches, but the end-fire approximation still displays a flat residual offset, as shown in Fig. 6(c). The differences in the performances of the three solutions are highlighted in Fig. 7, where |η e m -1| is shown to be only minimized in the case of optimal layer weights, while the asymptotic solution is only slightly better than the end-fire one.

Concerning the internal distortions, Fig. 8 confirms that the results obtained with the optimal and asymptotic solutions are very similar, while the end-fire solution is again characterized by a flat offset that explains the appearance of a relatively important inner radiation, as shown in sec. 8.3. Again, setting auxiliary currents further away from the original source region

(a) R 1 /W s = 1/2; (b) R 1 /W s = 1 and (c) R 1 /W s = 3/2.
strongly reduces the absolute value of the internal distortions; at the same time, the internal distortions also increase at a lower rate. It is worth observing that this trend is not shared by the end-fire solution, as its performance over low-order modes remains virtually unchanged.

Field errors on test surfaces

The modal distortions translate into errors on the electromagnetic field generated by the auxiliary sources, which can be measured on test surfaces as discussed in sec. 4. An example of typical results is given in Figs. 9 and 10, for the case R 1 /W s = 1/2 for end-fire and optimal procedures. In these two cases, internal (external) test surfaces were defined at half a wavelength inside (outside) the innermost (outermost) layer of auxiliary electric currents. These results show that although the optimal procedure is expected to ensure a better performance,

i m }, computed in the case W s = 4λ with ∆ R = λ/4, for R 1 /W s = 1/2 (black lines) and R 1 /W s = 3/2 (red lines).
the residual errors over the internal test surface can be locally stronger than by applying the end-fire approach. Moreover, the optimal procedure ensures that these errors are on average equal to zero, as opposed to the end-fire approximation; as it will be shown in sec. 8.3, these coherent deviations lead to residual focusing over the inner region.

Residual errors can be summarized by defining relative errors for the electric and magnetic fields

ǫ 2 E ,H (R T ) = e 2 E ,H (R T ) E E ,H , (35)
where E E ,H stands for the square of the L 2 norm of the electric (magnetic) field over the external test surface Σ e T . Recalling (15), ( 17) and ( 21) where E E ,H are also taken as references for the internal errors, in order to assess the relative intensity of the residual field generated over the inner region Ω ′ i . These figures of merit are shown in Figs. 11 to 13. Several trends can be pointed out. First of all, as ∆R approaches λ/2, residual errors increase. This result can be understood by noticing that, under an asymptotic approximation, the field distributions observed over two layers associated to each mode will be proportional one to the other if they are half a wavelength apart. Hence, rather than conveying complementary information that could be made to interfere, the two sets of auxiliary currents become increasingly redundant.

ǫ 2 E ,H (R e T ) = ∞ m=-∞ |η e m -1| 2 |β E ,H m (R e T )| 2 ∞ m=-∞ |β E ,H m (R e T )| 2 (36a) ǫ 2 E ,H (R i T ) = ∞ m=-∞ |η i m | 2 |β E ,H m (R i T )| 2 ∞ m=-∞ |β E ,H m (R e T )| 2 , (36b)
At the same time, as the ratio R 1 /W s increases, the errors systematically reduce. This behavior is easily understood within the framework introduced in sec. 6: in the limit for k 0 R 1 ≫ 1, all modal coefficients can be controlled at the same time. Fundamentally for the same reason, errors on the magnetic and electric fields over the external test surface are practically identical, while this is not true over the internal test surface, where the asymptotic expressions do not hold.

As expected, all errors are consistently lower for the optimal procedure. More importantly, the optimal procedure takes better advantage of a larger R 1 /W s , with a faster rate of reduction in the errors, when compared with the two approximate cases. The performance of the asymptotic approximation approaches that of the optimal procedure for larger R 1 /W s , since the former is correct in the limit for k 0 R 1 → ∞.

Field topographies

An example of the field topographies generated by two-layer auxiliary currents is shown in Fig. 14, emulating the linear source (34) when W s = 4λ. The three sets of electric and magnetic field distributions are to be compared to the original case reproduced in Fig. 4(b).

Differences between the performance of the three solutions are shown in Figs. 15 and16, as percentages of the peak field amplitude observed over the external region Ω ′ e ; two values of R 1 /W s are considered, namely 0.5 and 1.5. Independently from the dimensions of the shell-region, it is found that the end-fire approach always leads to a residual (coherent) focusing within the inner region Ω ′ i , as opposed to an incoherent interference pattern observed with the other two solutions. Errors in the external region drop when applying the optimal procedure.

Cross comparisons

This section draws some cross comparisons between the performance as already presented and the layer weights that generated them. The first such comparison is presented in Table 1 where the layer weights applied for the three methods are summarized. The weights {A l } are presented in a normalized form, as A 2 /A 1 .

We first need to recall that the weights of the approximate solutions only depend on the relative distance between the layers. By the same token, the approximate solutions are by definition only operating by phase-shifting the individual radiations from the auxiliary layers. Conversely, for the optimal procedure, this property is never explicitly enforced; the fact that all the results in Table 1 involve |A 1 /A 2 | = 1 is therefore worth noticing, as it implies that there is no need for amplitude-compensation schemes. and the internal test surface Σ i T (lower graph), with R i T = R -λ/2, for a varying distance ∆R between two current layers. Four values of R/W s are considered, equal to 0.5, 0.75, 1.0 and 1.5. Results for a source width W s = 4λ, optimal layer weights.

As a consequence, it appears that the major differences in the performance between the approximate and optimal procedures are in fact due to minor differences in the phase-shift angle enforced between the layers. While the asymptotic solution requires an in-phase radiation of the two layers, the optimal one spans phase-shift angles between 0.26 and 2.2 degrees, for R 1 /W s passing from 1.5 to 0.5. Correspondingly, the external error decreases more steeply for the optimal procedure, with a relative improvement of about a factor 2 to 10 with respect to the asymptotic solution.

Table 2 seeks to point out how the choice of electric or magnetic weights {β m } affects the performance of the optimal procedure. To this end, the layer weights and errors were computed for the two choices of {β m }: all the quantities in Table 2 noted by the letter 'R' stand for the ratio of results obtained by choosing {β H m } and those obtained with {β E m }; e.g., R A is the ratio of A 2 /A 1 obtained in the two cases, of which only the marginal phase-shift angle is shown, as expressed in degrees. All other ratios are self-explaining.

It can be concluded that the kind of {β m } chosen has little impact on the error metrics on the external test surface, with variations of the order of the percent point; this trend is reinforced when passing from W s = 2λ to W s = 4λ. Moreover, the same conclusions hold for both electric and magnetic field external errors. As opposed to these conclusions, the internal errors are more sensitive to the choice of the {β m } weights.

The asymptotic proportionality between {β E m } and {β H m }, discussed in sec. 5, is confirmed, with residual phase-shift angle of R A passing from 1.8 to 0.23 degrees for W s = 2λ.

A last comparison is needed between the performance of the Table 1: Absolute value and phase shift of the relative weight A 2 /A 1 , for the three methods of computation of the layer weights, for the case W s = 4λ and ∆R = λ/4. The corresponding relative errors ǫ 2 E (R T ) give a measure of how such small adjustments in the layer weights lead to large differences in the quality of reconstruction of the field distribution. Shorthand notations are used for the end-fire (ef), asymptotic (as) and optimal (opt) solutions. external-field optimization and the joint internal-external one. It was suggested in sec. 6 that, asymptotically, reproducing the original radiation over the external region must pass for a null radiation over the internal region. As such, a joint optimization can be regarded as redundant, as long as L ≪ M. In this respect, the results shown in Fig. 17 support this analysis, since the internal modal distortions are practically unmodified when switching between the two optimization strategies. The joint optimization should therefore be considered only in those cases where L ≃ M.

R 1 /W s Type |A 2 /A 1 | ∠A 2 /A 1 ǫ 2 E (R e T ) ǫ 2 E (R i T ) 0

Conclusions

This paper analyzed the possibility of defining auxiliary current distributions, extending the idea of Love's equivalent theorem to the case of electric-only currents. Rather than proceeding through the definition of an inverse problem, our proposal is based on a hybrid solution, keeping the electric current distribution defined by Love's theorem. The role of magnetic currents is taken over by further layers of electric currents weighted in such a way as to approximate the results expected for the equivalence theorem. The resulting procedure ensures a good performance, with residual errors of the order of a few percent points.

Figure 1 :

 1 Figure1: Original source of radiation (a) and two equivalent-current configurations, with a homogeneous background medium using electric and magnetic currents (b) and with only electric currents with the internal region filled by a perfect magnetic conductor (c).

Figure 3 :

 3 Figure 3: An example of electric (b) and magnetic (d) fields radiated by a single layer of equivalent electric currents, here represented by the dashed circle. Reference field distributions expected for the case of electric and magnetic equivalent currents are shown on the left column, respectively for the electric (a) and magnetic (c) field. Corresponding field distributions generated by the electric currents alone are given on the right column. Data represent the amplitude of the fields measured in dB, normalized to the peak amplitude observed for the reference distribution.

Figure 4 :

 4 Figure 4: Longitudinal electric-field distribution generated by original sources with: (a) W s = 2λ and (b) W s = 4λ. The inner regions with null radiation marks the minimum-sized cylinders containing the original sources.

Figure 5 :

 5 Figure 5: Moduli of the modal coefficients {γ m } for the case of W s = 2λ and W s = 4λ, serving as references.

Figure 6 :

 6 Figure 6: Real and imaginary parts of the overall external-field distortions {η e m }, computed in the case W s = 4λ and ∆ R = λ/4 for: (a) R 1 /W s = 1/2; (b) R 1 /W s = 1 and (c) R 1 /W s = 3/2.

Figure 7 :

 7 Figure 7: The distance |η e m -1|, computed in the case W s = 4λ and ∆ R = λ/4 for R 1 /W s = 1.

Figure 8 :

 8 Figure 8: Absolute value of the overall internal-field biases {η i m }, computed in the case W s = 4λ with ∆ R = λ/4, for R 1 /W s = 1/2 (black lines) and R 1 /W s = 3/2 (red lines).

Figure 9 :

 9 Figure 9: Real and imaginary parts of the longitudinal electric field for W s = 4λ, R 1 /W s = 1/2 and ∆R = λ/4, computed over cylindrical external (top figure) and internal (bottom figure) test surfaces, with R i T = R 1 -λ/2 and R e T = R 2 + λ/2, respectively. Results obtained under the end-fire approximation are shown as solid lines, reference results as dots.

Figure 10 :

 10 Figure10: Same configuration as in Fig.9. Results obtained with optimal layer weights.

5 Figure 11 :

 511 Figure 11: Comparison between relative errors on the electric and magnetic fields over the external test surface Σ e T (upper graph), with R e T = R + ∆R + λ/2

5 Figure 12 :

 512 Figure12: Same configuration as in Fig.11; results for the asymptotic solution.

Figure 13 :Table 2 :

 132 Figure13: Same configuration as in Fig.11; results for the end-fire approximation.

Figure 14 :Figure 15 :

 1415 Figure 14: Field distributions generated by two-layer auxiliary currents with R 1 = W s /2 and spacing ∆R = λ/4, for the source (34) with W s = 4λ. The modulus of the longitudinal electric (left column) and transverse magnetic (right column) fields, expressed in dB and normalized to their peak value, for the (a)-(b) end-fire, (c)-(d) asymptotic and (e)-(f) optimal solutions. The auxiliarycurrent layers are shown as dashed circles.

Figure 16 :

 16 Figure 16: Same as in Fig. 15, for the case R 1 /W s = 3/2.

Figure 17 :

 17 Figure 17: Absolute value of the overall internal-field modal distortions {η i m }, computed in the case W s = 4λ with ∆ R = λ/4, for R 1 /W s = 1/2 (black lines) and R 1 /W s = 3/2 (red lines).

  but weighting the individual errors according to the energy contributed by each mode, as assessed by |β E m | 2 . The LS solution to (24) is thus given by

Acknowledgement

Part of the work here reported was funded by the French National Research Agency through the grant ANR-12-ASTR-0005, within the framework of the 2012 ASTRID program.