Thomas Fromenteze
email: thomas.fromenteze@unilim.fr

Ettien L Kpré

Cyril Decroze

David Carsenat

Okan Yurduseven

Mohammadreza Imani

Jonah Gollub

David R Smith

Unification of compressed imaging techniques in the microwave range and deconvolution strategy

In recent years, several techniques of imaging have emerged based on compressive sensing. Thus, using passive devices able to code and sum the information measured on a radiating aperture, complex UWB beamforming hardwares can be replaced by inverse problems to solve using deconvolution techniques. Two similar compressive techniques have been developed independently, based on a passive multiplexer and on a metamaterial aperture. This paper proposes an unification of these approaches, using a common mathematical formulation and comparing the deconvolution strategy adopted by each research team.

I. INTRODUCTION

The domain of radar imaging in the microwave range has received a substantial attention in recent years thanks to its application in various domains such as medical diagnosis and concealed objects detection. However, UWB signals and large radiating apertures are required to reach high resolutions in both range and cross-range axis, imposing heavy technical constraints on the active devices required for the generation and reception of waves. In recent years, the development of the compressive sensing theory allowed new imaging techniques to emerge. This new principle is based on the use of a coding aperture, made of a set of radiating elements distributed randomly each with a different resonant frequency. Through experimental characterization of the coded aperture, its complex radiation pattern at each frequency can be predicted, making it possible to reconstruct an image by solving an inverse problem using the aperture's transfer functions and the measured compressed signal at its output. Fig. [START_REF] Carsenat | UWB antennas beamforming using passive time-reversal device[END_REF] depicts a generic compressive system.

The waves radiated by an antenna are reflected on a target of spatial reflectivity M (r) and measured by the coding aperture's elements. The received waves are coded depending on the radiating elements' positions and complex weights, corresponding to the compressive device's transfer functions C(r R X , f ) and summed at the unique output. The expression of the measured compressed signal g(f ) is, assuming the first Born approximation:

g(f ) = r R X r G(rT X , r, f ) M (r) G(r, rR X , f ) C(rR X , f ) ∂r ∂rR X (1)
with:

Compressed output g(f )

C o d i n g a p e r t u r e Radiating elements at f0

C(r R X , f 0) Tra nsf er fun ctio ns • Target M (r) + 0 r T X rR X + r G( rTX , r, f0 ) G(r, r RX , f 0)
Fig. 1: Imaging system composed of a transmitting antenna and a receiving compressive aperture.

• r R X : the positions of the coding aperture's surface.

• r T X : the position of the transmission antenna.

• G(r 1 , r 2 , f ): the Green function between r 1 and r 2 .

Thus, the aim of the compressed imaging technique is to compute an estimation of M (r) depending on the measured signal g(f ) and on the pre-characterized compressive device's transfer functions. This inverse problem is often confronted to ill-conditioning issues, source of inaccuracy in the reconstruction process. Thus, regularization techniques must be applied to improve the estimation of the spatial reflectivity of the scene.

In the following sections, two techniques that have been developed independently are presented and compared. A first one is based on a passive multiplexer connected to a conventional antenna array while the second one is a metamaterial radiating aperture.

II. PASSIVE MULTIPLEXER

This technique has been introduced in [START_REF] Carsenat | UWB antennas beamforming using passive time-reversal device[END_REF] and is based on a passive multiplexer, which is a 1-to-N ports highly dispersive filter with lowly correlated transfer functions. These properties are reached using the multipath propagation is a closed cavity, able to supply different frequency selectivity depending on the input's port position. Thus, the modal diversity represents here the support of the compression and requires an important volume of cavity. However, compact devices has been introduced in [START_REF] Fromenteze | Miniaturized Device for Passive Microwave UWB Beamforming[END_REF][3] using 2D chaotic cavities and in [START_REF] Fromenteze | Waveform coding for passive multiplexing: Application to microwave imaging[END_REF] with the conception of a 2D photonic crystal. Fig. (2) depicts a compressive device connected to an array of UWB Vivaldi antennas. Waves are radiated by a transmission antenna, reflected by a target and received by the compressive device. Thus, a unique signal is measured for the reconstruction of the scene.

Compressed output rRX r TX r • Target G (r T X , r, f ) G (r , rRX , f )
Fig. 2: Illustration of a passive connected to an array of UWB Vivaldi antennas.

The characterization of the compressive device's transfer function is performed between the access ports of the antennas and the unique output, using cabled measurement tools. For this technique, the first step is to reconstruct the waveforms received on the compressive device's aperture to make them repropagate in the space. The expression of these reconstructed signals is:

S rec (r R X , f ) = [C(r R X , f )] + g(f ) (2)
where [C(r R X , f )] + is a pseudo-inverse of the matrix [C(r R X , f )] computed for each frequency independently. In [START_REF] Carsenat | UWB antennas beamforming using passive time-reversal device[END_REF], matched filtering was used, leading to [C(r R X , f )] = [C(r R X , f )] * . This pseudo-inverse corresponds to a simple phase compensation of the device's impact, and is the most stable solution since no magnitude inversion, source of the illconditioning problems, is involved. In [START_REF] Fromenteze | Waveform coding for passive multiplexing: Application to microwave imaging[END_REF], the Tikhonov regularization has been adapted to this application by computing the following pseudo-inverse:

[C(rR X , f )] + = ([C(rR X , f )] † [C(rR X , f )] + β[I]) -1 [C(rR X , f )] † (3)
with β, a positive regularization parameter added to a least-square pseudo-inverse to prevent the inversion of too low values. This computation is performed on a frequency by frequency basis, meaning that the dimension of

[C(r R X , f )] is (N R X × 1) in this calculus.
In [START_REF] Fromenteze | Miniaturized Device for Passive Microwave UWB Beamforming[END_REF], the truncated singular value decomposition (tSVD) has been introduced as an alternative to improve the conditioning of the inverse problem. A matrix factorization is used:

[C(r R X , f )] = [U (r R X )] f [Σ(r R X )] f [V (r R X )] † f (4)
where [U (r R X )] f and [V (r R X )] f are unitary matrix made of orthonormal basis weighted by the singular values of the quasi-diagonal matrix [Σ(r R X )] f {σ 1 , σ 2 , ..., σ N }, arranged in descending order. The conditioning is improved by setting the lowest singular values to zero and computing the following pseudo-inverse:

[C(r R X , f )] + = [V (r R X )] f [Σ(r R X )] -1 f [U (r R X )] † f (5) with [Σ(r R X )] -1
f a quasi-diagonal matrix filled with { 1 σ1 , ..., 1 σ β , 0, ..., 0}. All the remaining singular values must satisfy σ > β, with β the regularization parameter. The tSVD and the Tikhonov pseudo-inverse are closely related, the singular value decomposition of the inversed matrix introduced in eq. ( 3) leads to:

diag([Σ(r R X )] Tikhonov ) = { σ 1 σ 2 1 + β , σ 2 σ 2 2 + β , ...} (6)
Thus, a similar truncation of the lowest singular values is applied depending on the regularization parameter's value. The main advantage of the tSVD is the possibility to display the useful elements of the linear problem to apply the truncation, but it requires a more expansive computation to perform the matrix factorization at each frequency. Once the signals received by the aperture have been estimated, conventional imaging algorithm can be computed with them to reconstruct an image. The MIMO range migration algorithm introduced in [START_REF] Zhuge | Three-Dimensional Near-Field MIMO Array Imaging Using Range Migration Techniques[END_REF] is based on fast Fourier transforms and has been shown efficient, even with 3D reconstructions [START_REF] Fromenteze | Passive UWB multiplexing device for a single-port MIMO RADAR[END_REF].

III. METAMATERIAL APERTURE

An other technique introduced in [START_REF] Hunt | Metamaterial apertures for computational imaging[END_REF] and developed in [START_REF] Lipworth | Metamaterial apertures for coherent computational imaging on the physical layer[END_REF] uses a radiating aperture made of coupled resonant cELC unit cells able to radiate the required coded patterns depending on the frequency. This metamaterial aperture is associated with a transmission horn matched in the same bandwidth (Fig. 3). The study of the compressive device's radiation has been first developed by modeling the field from feeding coaxial probes within parallel-plate waveguide and their couplings to radiating cELC unit cells considered as magnetic dipoles. The measurement of the compressive device's transfer functions can also be performed by a near field scan. A cabled measurement is not possible here since the response of the coupled radiating elements is the main source of dispersion, but the near field scan allows for a more accurate measurement of the complete system's radiation. The deconvolution strategy is slightly different from the one presented in the previous section. The radiation of the compressive aperture is discretized over a set of voxels in a region of interest:

Compressed output r R X r T X G ( r T X , r , f ) G(r , r R X , f ) r • Target
[H(r, f )] = r R X [G(r T X , r, f )].[G(r, r R X , f )].[C(r R X , f )] (7)
where "." stands for the element-wise multiplication operator. Thus, the reconstructed image M rec (r) is evaluated in the following way:

M rec (r) = [H(r, f )] + g(f ) (8)
For this technique, several types of pseudo-inverse has also been studied for the computation of [H(r, f )] + . In [START_REF] Hunt | Metamaterial apertures for computational imaging[END_REF], the iterative algorithm TwIST is adapted to this problem using a-priori information on the sparsity of the imaging scene to optimize the estimation. Matched filtering and l 2 -norm leastsquares methods have also been applied for the imaging of continuous targets. This deconvolution method differs from the one presented in the previous section since the complete system's radiation is compensated by the pseudo-inversion, while the previous approach was split in two steps to prevent multiple pseudo-inversions of the compressive device's transfer functions that are independent of the scene to image. Thus, more accurate results can be expected with this second technique, at the cost of a longer computation time. In the following section, a Matlab simulation is presented to compare the two approaches in terms of accuracy and computation time. Thus, a generic compressive device is simulated, associated with an isotropic transmission antenna.

IV. COMPARISON OF THE DECONVOLUTION METHODS

An aperture made of 30 × 30 radiating elements regularly spaced on a 1m 2 square are randomly weighted depending on the frequency is studied in the 2-10 GHz range. The compressive aperture is centered on the xz-axis and a transmission antenna is set on the side of the array at the coordinate (x = -0.6m, z = 0m). For the same scenario, a comparison of point spread functions (PSF) is presented depending on the deconvolution approach. The frequency vector is sampled by 1001 points and the imaging scene is discretized by 90×13×90 voxels in a volume delimited by x ∈ [-0.5, 0.5]m, y ∈ [0.31, 0.76]m, z ∈ [-0.5, 0.5]m. In each simulation, the target is set in the center of the y-axis, at y = 0.5m. Two 3D PSF are presented in Fig. (4) and Fig. ( 5) using the deconvolution techniques developed for the passive multiplexer and for the metamaterial aperture respectively. For both cases, the reconstructed field is normalized and the y = 0.5m plane is extracted to compare the level of sidelobes. This simulation is performed applying the first Born approximation with a punctual target set at (x = 0m, y = 0.5m, z = 0m). Both reconstruction are performed using Tikhonov pseudo-inversions with optimized regularization parameters. The results are obtained running Matlab on a computer with a 3.5 GHz 6-core processor and 16 Gb of RAM. The average level of sidelobes is lower in the case of a pseudo-inversion of the radiation compared to the case of a reconstruction of the aperture's field re-propagated with a range migration algorithm. Indeed, the RMA is based on a stationary phase approximation allowing for a reduction of the required computing time at the price of slight distortion of the images. 1D cuts for x = 0m and z = 0m are extracted from the y = 0.5m plane and compared in Fig. [START_REF] Lipworth | Metamaterial apertures for coherent computational imaging on the physical layer[END_REF]. For the studied scenario, the computation time required with each deconvolution method are gathered in table (I). In each case, a part of the algorithm is independent of the imaging scene and is pre-computed. In the case of a spatial deconvolution, a 1.5 Gb matrix [H(r, f )] is pre-computed for the depicted scenario and saved, working with double-precision floats. The longer computation time for this case is explained by the matrix loading required before the reconstruction of a scene. To conclude, a simulation of a gun shape is presented with the same radiating aperture to show the imaging capacity of this compressive system. The image obtained for conventional digital beamforming using the signals measured from each receiving element is compared with the one obtained with only one compressed signal measured at the output of the imaging device (Fig. 7). For this case, a random matrix of 10001 frequency points has been used, associated with the range migration algorithm. The estimation of the 30 × 30 received signals can be performed accurately enough to reconstruct a 3D gun shape close to the non-compressed case.

V. CONCLUSION

In this paper, the similarities between two techniques of compressed imaging in the microwave range have been presented. Although independently developed, this two techniques are basically the same and share a common mathematical formulation. The use of a near field scan seems to be the best way to characterize the imaging system since it allows for Fig. 7: Imaging of a gun using conventional digital beamforming (left), compared to compressed imaging (right).

an accurate measurement of the radiation patterns, including the inter-element coupling that is not compensated for the case of the passive multiplexer. Furthermore, the use of range migration algorithm should be adapted to the metamaterial apertures when fast imaging is required, subject to a slight loss of accuracy. This new technology represents a promising way to simplify the complex hardware architectures required for high resolution microwave imaging. Future efforts in each research group will be focused on improving the imaging capabilities of these techniques, implementing compressive aperture to aperture imaging [START_REF] Fromenteze | Passive UWB multiplexing device for a single-port MIMO RADAR[END_REF], GPU acceleration, and polarization imaging.

Fig. 3 :

 3 Fig. 3: Metamaterial aperture made of coupled cELC unit cells, associated with a radiating horn.

Fig. 4 :

 4 Fig. 4: 3D point spread function for a reconstruction of aperture's field and a range migration algorithm and 2D extraction from the y = 0.5m plane

Fig. 6 :

 6 Fig. 6: Comparison of the PSF for 1D cut along the x-axis and z-axis.

TABLE I :

 I Comparison of the computation time depending on the chosen deconvolution strategy

		Pre-computation time	Computation time
	Aperture reconstruction & RMA	4.8 s	2.4 s
	Spatial deconvolution	38 min	66 s


Thomas Fromenteze
email: thomas.fromenteze@unilim.fr

Ettien L Kpré

Cyril Decroze

David Carsenat

Passive UWB Beamforming: a N to M Compression Study

Recent works have demonstrated the feasibility of microwave imaging using compressive techniques, exempting the use the of active delay lines, phase shifters, or moving parts to achieve beamforming. With this method, waves are coded in a passive way by a compressive device to reduce the complexity of the transmitter and/or receiver chains of the telecommunication and radar systems requiring beamsteering. Since this technique is based on frequency diversity, the reduction of the compressive device's volume imposes a diminution of the amount of driven antennas. In this article, the improvement brought by simultaneous excitations of the compressive device is presented. Adapting a new mathematical formulation, it is shown that M inputs can send independent waveforms allowing the beamsteering of an N-elements antenna array, while maintaining N > M .

I. INTRODUCTION

UWB radar techniques are used in a wide range of applications such as medical diagnosis and concealed objects detection. These techniques often take advantage of the possible transmission of microwaves in materials that are totally opaque in the optical domain. The capability of an imaging system to resolve targets depends of two parameters:

• The down range resolution δ r = c

2B

• The cross range resolution δ cr = min[ c 0.7079B , λ c ] R L [START_REF] Zhuge | Short-range ultra-wideband imaging with multiple-input multiple-output arrays[END_REF] with c the wave celerity, B the absolute bandwidth, λ c the central wavelength, R the range of the focal point, and L the size of the radiating aperture. Thus, an high resolution imaging system requires the use of UWB hardwares connected to a large amount of antennas, representing an heavy burden regarding its cost and complexity. In the time domain, the aim of these active devices is to generate delays applied to each radiating elements and achieve beamsteering in the desired direction, as depicted in Fig. [START_REF] Zhuge | Short-range ultra-wideband imaging with multiple-input multiple-output arrays[END_REF]. It is worth noting that this principle can be applied in transmission as well as in reception.

In [START_REF] Carsenat | UWB antennas beamforming using passive time-reversal device[END_REF], a compressive technique has been presented to achieve beamforming with UWB signals in a passive way. This technique is based on the use of a passive multiplexer controlled by its input waveform, connected to an antenna array. This principle has been extended in [START_REF] Fromenteze | Waveform coding for passive multiplexing : Application to microwave imaging[END_REF], developing a compact compressive device and a new formulation of signal reconstruction. An other technique based on compressed sensing has been introduced in [START_REF] Hunt | Metamaterial apertures for computational imaging[END_REF] and extended in [START_REF] Lipworth | Metamaterial apertures for coherent computational imaging on the physical layer[END_REF], using a radiating device made of coupled resonators able to transmit signals in different directions of the space simultaneously depending on the frequency. Both these techniques are basically the same, exploiting the frequency diversity of a passive device to radiate orthogonal patterns at each frequency. A new study is presented in the next section, focused on the improvements

1 2 3 N d x y θ 0 τ 1 τ 2 τ 3
Fig. 1: Illustration of a time-domain beamforming: the delay applied to each radiated pulses depends on the direction of radiation and of the inter-element spacing.

brought by using several receiving chains connected to the compact compressive device introduced in [START_REF] Fromenteze | UWB passive beamforming for large antenna arrays[END_REF]. The comparison between conventional and compressed beamforming is depicted to demonstrate the compressive capability of this technique and the modal diversity gain provided by using simultaneous measurements of independent coded waveforms.

II. MULTIPLE-PORTS COMPRESSIVE BEAMFORMING

A. Compressed acquisition theory

The study presented in this article is based on a compression from a set of N signals received by an antenna array, and compressed into M signals after the propagation in a coding device, with M < N (Fig. 2) The expression of Y k (f ), the signal compressed at the device's output k is defined in eq. ( 1):

Passive multiplexer S 1 (f ) . . . S n (f ) . . . S N (f ) H( 1, 1, f ) Y 1 (f ) H(n, m, f ) Y m (f ) H( N, M, f ) Y M (f )
Y k (f ) = N n=1 H(n, k, f ) S(n, f ) (1)
with H(k, m, f ) the compressive device transfer function between the ports k and m and S(m, f ) the signal received by the m th antenna. This compression corresponds to a generalization of the experiments presented in [START_REF] Carsenat | UWB antennas beamforming using passive time-reversal device[END_REF][3][6] [START_REF] Fromenteze | UWB passive beamforming for large antenna arrays[END_REF], where the output was unique. In this article, the use of multiple outputs is studied to evaluate the gain for beamforming performed with a large number of antennas. Following the compressive sensing theory, the sensing matrix

[S(f )] rec = [H(f )] + [Y (f )] (2)
[H(f )] + is a pseudo-inverse computed by truncated singular value decomposition (tSVD) of [H(f )]. In this purpose, the matrix [H(f )] is factorized into three matrix for each frequency:

[H(f )] N ×M = [U (f )] N ×N [Σ(f )] N ×M [V (f )] † M ×M (3)
with N the number of antennas and M the number of receiving chains. [U (f )] and [V (f )] are unitary matrices made of orthonormal set of basis, weighted by the singular values ordered in the quasi-diagonal matrix [S(f )]. The conditioning of this linear application is defined by the following equation:

cond([H(f ]) = σ 1 σ N (4)
σ 1 and σ N stand for the largest and the smallest singular values of [H(f )]. The condition number cond([H(f ]) relates the accuracy of the solution given by the pseudo-inversion. An ideally conditioned problem would involve a condition number of 1, meaning that all the singular values have the same significance and are not sensitive to the inversion. For an ill-conditioned problem, the lowest singular values can be set to zero to improve the accuracy of the pseudo-inversion, expressed as :

[H(f )] + = [V (f )][Σ(f )] -1 t [U (f )] † (5)
where

[Σ(f )] -1
t is a quasi-diagonal matrix filled by { 1 σ1 , ..., 1 σ β , 0, ..., 0}. All the remaining singular values must satisfy σ > β, with β a regularization parameter optimized empirically depending on the intrinsic conditioning of the problem. By computing this pseudo-inversion frequency by frequency, the inverse matrix [H(f )] + allows for a reconstruction of the waveforms received by the antenna array from the measured signals. Finally, radiation patterns C(θ, f ) can be computed with the reconstructed signals [S(f )]. In the case of a far field approximation, the beamforming is computed in a conventional way by correlating the waveforms with plane waves :

C(θ, f ) = N n S rec,n (f ) exp j 2πf c d(n -1) sin(θ) (6)
with d the inter-antenna spacing and θ the beamsteering direction. In the next section, this principle is applied to a compact device connected to an antenna array to illustrate the beamforming improvement compared to conventional beamforming techniques applied with the same amount of receivers. [START_REF] Fromenteze | UWB passive beamforming for large antenna arrays[END_REF] is used as a passive multiplexer (Fig. 3). This device is a microwave substrate RT/duroid 6006 with a hole engraved on the its top conductor. 20 SMA ports are regularly connected on its edges, with small distance variations to break periodicity. Despite its simple aspect, this microwave device presents very useful properties linked to the association of regular and convex boundaries delimiting the wave propagation [START_REF] Draeger | One-channel time -reversal in chaotic cavities: Experimental results[END_REF] [START_REF] De Rosny | Coherent backscattering of an elastic wave in a chaotic cavity[END_REF]. These properties can be illustrated by replacing the field propagating in a cavity by its equivalent plane wave spectrum. The path followed by a plane wave can be represented as a ray bouncing on the cavity's boundary (Fig. 4). Starting from the same initial conditions, a ray in a chaotic cavity tends to fill the entire space after several reflexions compared to a regular cavity where the ray are focused in particular areas. Furthermore, a slight modification of the conditions of a ray in a chaotic cavity can have an important impact in the way it propagates, unlike the case of the regular cavity. Thus, even with a limited frequency diversity which represent the support of the signals compression, these properties help to decrease the correlation level between the device's transfer functions. The aim of this study is to demonstrate the interest of using several outputs simultaneously to achieve a better reconstruction of the signals when a large amount of antennas are connected to the compressive device. Thus, the transfer function between 16 inputs and 4 outputs are measured in the 2-4 GHz frequency range (300 points) to evaluate the singular value spectrum depending on the amount of connected outputs (Fig. 5). More information is measured when new outputs are connected, leading to higher singular values and to an improved conditioning for a given regularization parameter.

B. Application to a compressive device 1) Study of the singular values:

A 2D chaotic cavity introduced in [6]
2) Compressed beamforming: In conventional digital beamforming, as many receivers as antennas are required. With the same amount of receiving chains, the compressive beamforming can be performed with a larger number of antennas connected to the passive multiplexer. For this example, 16 antennas are connected to the chaotic device and compressed into 4 waveforms (Fig. 6). The beamforming capabilities of this system are compared with a conventional digital beamforming array made of four elements. These experiments are simulated by Matlab programs using the device's measured transfer functions in the frequency range 2-4 GHz. For both cases, the inter-antenna space is set to d = 0.7 * λ c = 7 cm . The figure [START_REF] Fromenteze | UWB passive beamforming for large antenna arrays[END_REF] depicts the results obtained for a beam-steering in the 25 • direction, for a conventional system made of four elements and for a compressive system connected to 16 antennas. In both cases, the antennas are isotropically radiating and the mutualcoupling is not considered. The radiation patterns computed with eq. ( 6) are converted from frequency to time domain data using an inverse Fourier transform. This figure shows that using the same number of receiving chains and inter-antenna space, the compressive technique allows for a thinner formed beam since the equivalent radiation aperture is wider and appropriately sampled. Furthermore, the grating lobes observed for the beamsteering with the conventional system are mitigate by the compressive device, due to a statistical cancellation of the off-axis frequency components [START_REF] Fromenteze | UWB passive beamforming for large antenna arrays[END_REF]. The time required for the computation of the tSVD pseudoinverse was less than 3.5 × 10 -2 s and can be done only once and for all since the compressive device's transfer functions are stationary. The estimation of the received signals is calculated in less than 2 × 10 -3 s. Radiation patterns are extracted at the focusing time t = 0 ns and presented in Fig. [START_REF] Draeger | One-channel time -reversal in chaotic cavities: Experimental results[END_REF]. 3) Simultaneous measurements and beamforming: The last part of this study is focused on the impact of the amount of connected outputs to achieve beamforming. An example is depicted in Fig. [START_REF] De Rosny | Coherent backscattering of an elastic wave in a chaotic cavity[END_REF] to show the improvement provided by increasing the amount of connected outputs for a beamsteering in the -60 • direction. Fig. 9: Improvement of beamforming depending on the amount of connected outputs, for a beamsteering in the -60 • direction.

In this example, the level of noise tends to go down when the amount of connected outputs increase, meaning that new useful information are measured and decompressed to get better estimations of the waveforms received by the antennas. Thus, the peak signal-to-noise ratio (PSNR) is defined to determine the efficiency of the reconstruction for each beamsteering direction Φ, depending on the amount of connected outputs:

PSNR(Φ) = -10 log 10   θ t [|C Φ (θ, t)|-|C Φ,ideal (θ, t)|] n θ n t   (7)
with C Φ (θ, t) the time domain radiation pattern computed with the compressive technique and C Φ,ideal (θ, t) the equivalent ideal radiation pattern calculated by plane wave projections over the position of the 16 antennas, both for a beamsteering in the direction Φ. Depending on the application, the trade-off between quality of the reconstruction and complexity of the associated hardware can be chosen. The computed results are gathered in Fig. (10). As expected, the quality of the beamforming directly depends on the amount of measured signals from the compressive device, allowing for a better of the received signals.

III. CONCLUSION

A compressive multiplexing technique has been studied in this paper, starting from previous works. The use of several outputs connected to a coding device has been shown efficient to increase the modal diversity required to achieve the compressed measurement of waveforms received by an antenna array made of 16 elements, over 4 receiving chains. Thus, the improvement compared to conventional digital beamforming using the same number of receivers and the same inter-element spacing has demonstrated the possible implementation of this technique to reduce the hardware complexity of the UWB radars.

Fig. 2 :

 2 Fig. 2: Compression of a set of N signals received by an antenna array into M signals measured to compute digital beamforming.

  [H] is sparse and must present a low level of correlation between its transfer functions, ensuring reduced inter-signals projections during the compressing process. Starting from the measured signals [Y (f )], the received waveforms [S(f )] are reconstructed to achieve digital beamforming using a pseudo-inversion of the matrix [H(f )] at each frequency. Thus, the estimation of the signals [S(f )] is:

Fig. 3 :

 3 Fig. 3: Planar compressive device connected to 20 SMA ports.

Fig. 4 :

 4 Fig. 4: Comparison between a regular and a chaotic cavity : propagation of a ray, starting from the same initial conditions

Fig. 5 :

 5 Fig. 5: Normalized magnitude of the singular values depending on the number of outputs connected to the compressive device.

Fig. 6 :

 6 Fig.6: Compressive beamforming system connected to 16 antennas, connected to four independent receiving chains. In the study of conventional beamforming, 4 antennas are directly connected to the receiver's outputs, while maintaining the same inter-antenna spacing.

Fig. 7 :

 7 Fig. 7: Comparison of the time domain radiation patterns computed for a conventional beamforming with 4 antennas compared to a compressed beamforming of 16 antennas connected to four receivers via the compressive device. The beamsteering is performed in the 25 • direction.

Fig. 8 :

 8 Fig. 8: Radiation patterns extracted at the focusing time from Fig. (7).

Fig. 10 :

 10 Fig. 10: Peak Signal to Noise ratio depending on the amount of connected outputs and of the steering angle.


Thomas Fromenteze
email: thomas.fromenteze@unilim.fr

Cyril Decroze

David Carsenat

Ettien L Kpré

Imagerie MIMO ULB à haute résolution par un composant passif à port unique

Une technique de compression permettant l'acquisition d'images RADAR 3D à haute résolution est présentée dans cet article. Basée sur un principe d'acquisition comprimée, les informations reçues par un réseau d'antennes MIMO ULB sont mesurées sur un port unique. Après reconstruction des formes d'ondes reçues en postprocessing, il est alors possible de réaliser des images 3D de façon rapide grâce aux développements récents en terme d'algorithmes MIMO.

Introduction

Les récentes améliorations des RADAR microondes à haute-résolution ont été rendus possibles grâce aux développements de systèmes MIMO (Multiple-Inputs Multiple-Outputs) fonctionnant en ultra-large bande, permettant ainsi un accroissement conjoint des quantités d'informations spatiales et temporelles mesurées [START_REF] Bliss | Multiple-input multiple-output (MIMO) radar and imaging : degrees of freedom and resolution[END_REF]. Les applications de tels systèmes sont multiples, allant du domaine médical à la détection d'armes enfouies. Cependant, l'implémentation de ces techniques à grande échelle est fortement limitée par le coût et la complexité des systèmes actifs ULB multipliés par le nombre de voies d'émission et de réception. Ainsi, une technique a été introduite dans [START_REF] Carsenat | UWB Antennas beamforming using passive time reversal device[END_REF] afin de réduire la complexité des systèmes actifs nécessaires au fonctionnement des RADAR ULB. Cette méthode de multiplexage passif utilise la propagation au sein d'un composant dispersif afin de compresser les formes d'ondes émises ou reçues. Ce principe a été étendu dans [START_REF] Fromenteze | Passive UWB multiplexing device for a single-port MIMO RADAR[END_REF] au cas des RADARs MIMO, réalisant un prototype fonctionnant de façon passive dans le domaine ultra large bande. La théorie associée à cette technique est présentée dans la section suivante, suivie de résultats expérimentaux.

Principe théorique

Un montage illustratif est présenté en figure [START_REF] Bliss | Multiple-input multiple-output (MIMO) radar and imaging : degrees of freedom and resolution[END_REF]. Une mesure en réflexion est réalisée à l'entrée du composant compressif de canaux [H(f )], reliant le port d'accès aux N antennes. Le signal mesuré S(f ) correspond alors au produit matriciel suivant : Afin d'extraire la réponse des cibles du signal comprimé mesuré en réflexion du composant passif, des matrices de pseudo-inverses sont d'abord calculées par une de Tikhonov :

S(f ) = [H(f )] 1×N [M (f )] N ×N [H(f )] T N ×1 , ( 1
[H(f )] + = ([H(f )] † [H(f )] + β[I N ]) -1 [H(f )] † (2)
où β est un paramètre de régularisation permettant d'améliorer le conditionnement du terme inversé et † est l'opérateur transposé-conjugué. Si les canaux [H(f )] sont orthogonaux entre eux, on a alors :

[H(f )] + [H(f )] = [H(f )][H(f )] + = [I N ] (3) Ainsi, une estimation [M r (f )] de [M (f )] peut être cal- culée : [M r (f )] = [H(f )] + [S(f )]([H(f )] + ) T (4)
A partir de l'équation (3), le développement suivant est réalisé :

[M r (f )] = [H(f )] -1 [H(f )] [I N ] [M (f )] [H(f )] T ([H(f )] -1 ) T [I N ]
(5) menant à une parfaite reconstruction de la réponse des cibles. En pratique, cette déconvolution ne peut pas être idéale pour une fréquence unique, mais l'effet constructif en temporel permet de s'en approcher en fonction du niveau de corrélation des canaux. Cette technique est donc d'autant plus efficace que la bande passante utilisée est importante, permettant d'augmenter la diversité fréquentielle des canaux et de limiter leurs niveaux de corrélation. Dans la prochaine section, la faisabilité de ce principe est démontrée par un montage expérimental.

Montage experimental et mesures

Le composant utilisé pour cette preuve de concept est une cavité reverberante de 0.6m 3 dans laquelle une port d'entrée et 20 ports de sortie sont reliés à des sondes ULB (Fig. 2). Le niveau de corrélation entre chaque couple de canaux est calculé (Fig. 4). Un rapport moyen de 17 relevé entre les niveaux d'autocorrélations et d'intercorrelations permet d'assurer une projection inter-signaux limitées lors de l'étape de reconstruction. L'image 3D de la sphère métallique est reconstruite, malgré la présence quelques résidus spatiaux probablement liés à de légères modifications des canaux [H(f )] depuis l'étape de caractérisation.

Dans la mesure où l'approche développée reste équivalente aux systèmes MIMO conventionnels, il est possible de disposer d'une ouverture synthétique plus grande en utilisant le même nombre d'antennes au moyen d'un réseau en croix [START_REF] Zhuge | Three-Dimensional Near-Field MIMO Array Imaging Using Range Migration Techniques[END_REF]. L'ouverture équivalente ainsi crée correspond à la convolution spatiale de la position des antennes d'émission et de réception et permet de mettre en évidence la faible redondance spatiale d'informations qui tend à apparaitre dans le cas du réseau rectangle utilisé en MIMO. La figure (8) présente le montage du réseau utilisé pour lequel l'espacement entre antennes de 6.25cm sur chaque branche est utilisé.

Cette mesure reproduit une expérience présentée dans [START_REF] Zhuge | Three-Dimensional Near-Field MIMO Array Imaging Using Range Migration Techniques[END_REF] où un algorithme rapide de Range Migration basé sur une décomposition en ondes planes est adapté au cas MIMO. Le but de cette expérience est ainsi d'exploiter la rapidité de cette approche conventionnelle et de l'adapter au cas compressif étudié. Après capture d'un unique signal mesuré par l'analyseur de réseau, la matrice d'interaction des antennes avec la scène est reconstruite puis traitée par l'algorithme d'imagerie présenté dans [START_REF] Zhuge | Three-Dimensional Near-Field MIMO Array Imaging Using Range Migration Techniques[END_REF]. Une forme d'arme à feu est utilisée pour cette expérience, afin de présenter les performances de cette technique pour des cibles plus complexes et d'introduire sa potentielle implémentation au domaine de la détection d'objets enfouis. L'image obtenue en 3D ainsi qu'une coupe dans le plan de la cible sont présentées en figure (9). Cette image a été calculée en moins de 2s grâce à l'implémentation de l'algorithme de Range Migration.

La reconstruction obtenue permet de retrouver la forme de l'arme à feu dans l'espace. Une rotation de 45 • de la scène est obtenue lors de l'application de l'algorithme rapide introduit dans [START_REF] Zhuge | Three-Dimensional Near-Field MIMO Array Imaging Using Range Migration Techniques[END_REF], dû à une reconstruction orientée selon les branches du réseau en croix.

Conclusion

Une méthode d'imagerie à haute résolution ne nécessitant que la mesure d'un unique signal a été réalisée. Cette preuve de concept permet de démontrer la possible réduction de complexité des architectures actives mises en jeux dans les systèmes d'imagerie microonde MIMO. Dans la mesure où la qualité de la reconstruction des signaux dépend de la diversité fréquentielle des canaux du composant compressif, cette technique est particulièrement bien adaptée au domaine ULB, voyant ses performances augmenter avec la bande passante utilisée contrairement à l'ensemble des systèmes existants. Le principe d'acquisition "one shot" permettrait l'implémentation de cette technique dans de nombreux domaines allant de la détection d'objets enfouis à l'imagerie médicale. Une fois les signaux reconstruits par pseudo-inversion, cette approche est compatible avec l'ensemble des algorithmes d'imagerie développés pour les systèmes conventionnels, permettant de s'adapter aux contraintes de temps de calcul et de complexité de cible, impliquant par exemple des formes complexes et des matériaux hétérogènes.

Un nouveau composant a été développé pour de futures expériences (Fig. 10). Les câbles reliant les antennes au composant compressif ont été supprimés dans la mesure où ils représentaient la plus grande source d'instabilité. De prochaines études théoriques et expérimentales seront ainsi menées afin d'améliorer la reconstruction des signaux compressés et la qualité des images hauterésolution obtenues par cette technique.

Figure 1 .

 1 Figure 1. Illustration d'un RADAR MIMO ULB à port d'accès unique. Les formes d'ondes émises et reçues par le réseau d'antennes sont codées par la propagation dans un composant compressif, permettant la mesure d'une unique forme d'onde.

Figure 2 .

 2 Figure 2. Représentation schématique de la cavité reverbérante utilisée comme composant dispersif.

Figure 3 .

 3 Figure 3. Caractérisation des canaux de la cellule réverbérante et vue des sondes ULB utilisées pour les créer.

Figure 4 .

 4 Figure 4. Niveau de corrélation maximum entre chaque couple possible de canaux et de pseudoinverses.

Figure 5 .

 5 Figure 5. Montage du RADAR MIMO à port unique

Figure 6 .

 6 Figure 6. Cible sphérique placée en face du réseau rectangulaire.

Figure 7 .

 7 Figure 7. Gauche : Image 3D reconstruite à partir du signal comprimé -Droite : Coupe à une distance de 50 cm.

Figure 8 .

 8 Figure 8. Réseau en croix utilisant le même nombre d'éléments que précédemment, connecté à la cellule réverbérante.

Figure 9 .

 9 Figure 9. Reconstruction d'une cible en forme d'arme à feu en 3D et extraction d'une image 2D dans le plan de la cible. La scène a subi une rotation de 45 • du fait de l'algorithme rapide d'imagerie utilisé [5].

Figure 10 .

 10 Figure 10. Nouveau système compressif appliqué à l'imagerie microonde MIMO ULB "One shot".

Références

XIX èmes Journées Nationales Microondes, 3-4-5 Juin 2015 -Bordeaux


Van-Tinh Tran
email: van-tinh.tran@univ-lyon1.fr

Alex Aussem
email: aaussem@univ-lyon1.fr

A Practical Approach to Reduce the Learning Bias under Covariate Shift

Covariate shift is a specific class of selection bias that arises when the marginal distributions of the input features X are different in the source and the target domains while the conditional distributions of the target Y given X are the same. A common technique to deal with this problem, called importance weighting, amounts to reweighting the training instances in order to make them resemble the test distribution. However this usually comes at the expense of a reduction of the effective sample size. In this paper, we show analytically that, while the unweighted model is globally more biased than the weighted one, it may locally be less biased on low importance instances. In view of this result, we then discuss a manner to optimally combine the weighted and the unweighted models in order to improve the predictive performance in the target domain. We conduct a series of experiments on synthetic and real-world data to demonstrate the efficiency of this approach.

Introduction

Selection bias, also termed dataset shift or domain adaptation in the literature [START_REF] Moreno-Torres | A unifying view on dataset shift in classification[END_REF], occurs when the training distribution P (x, y) and the test distribution P (x, y) are different. It is pervasive in almost all empirical studies, including Machine Learning, Statistics, Social Sciences, Economics, Bioinformatics, Biostatistics, Epidemiology, Medicine, etc. Selection bias is prevalent in many realworld machine learning problems because the common assumption in machine learning is that the training and the test data are drawn independently and identically from the same distribution. The term "domain adaptation" is used when one builds a model from some fixed source domain, but wishes to deploy it across one or more different target domains. The term "selection bias" is slightly more specific as it assumes implicitly that there exists a binary variable S that controls the selection of examples in the training set, in other words we only have access to the examples that have S = 1. For instance, case-control studies in Epidemiology are particularly susceptible to selection bias, including bias resulting from inappropriate selection of controls in case-control studies, bias resulting from differential loss-to-follow-up, incidence-prevalence bias, volunteer bias, healthy-worker bias, and nonresponse bias [START_REF] Hernán | A structural approach to selection bias[END_REF].

It is well known that one may account for the difference between P (x, y) and P (x, y) by re-weighting the training points using the so-called importance weight, denoted as β(x, y) = P (x, y)/P (x, y). Formally, let {h θ * } θ∈Θ be a model family from which we want to select an optimal model h θ * (x) = h(x, θ * ) for our learning task and let l(y, h(x, θ)) be the loss function we would like to minimize, the optimal model we are searching for is the one that minimizes the expected loss over the test (or target) distribution:

θ * = argmin θ∈Θ (x,y)∼P
β(x, y)P (x, y)l(y, h(θ, x))

So in practice, weighting the empirical loss of the training instances by β(x, y) provides a well-justified solution to the selection bias problem.

In general, the estimation of β(x, y) with two different distributions P (x, y) and P (x, y) is unsolvable, as the two terms could be arbitrarily far apart. One simple assumption we can make about the connection between the distributions of the source and the target domains is that P (x, y) and P (x, y) differ only in P (x) and P (x) while their conditional distribution P (y|x) remains unchanged. This specific selection bias is known as covariate shift in the literature [START_REF] Shimodaira | Improving predictive inference under covariate shift by weighting the log-likelihood function[END_REF]. In this case, the weighting term reduces to β(x) = P (x)/P (x) and effective adaptation is possible. At first glance, it may appear that covariate shift is not a problem because, for classification, we are only interested in P (Y |X) which remains unchanged. In fact, Shimodaira [START_REF] Shimodaira | Improving predictive inference under covariate shift by weighting the log-likelihood function[END_REF] showed that there are circumstances under which the predictive performance is jeopardized by covariate shift. This happens typically when the parametric model family {P (Y |X, θ)} θ∈Θ is misspecified, that is, there does not exist any θ ∈ Θ such that P (Y |X = x, θ) = P (Y |X = x) for all x ∈ X , so none of the models in the model family can exactly match the true relation between X and Y .

The intuitive reason why covariate shift under model misspecification is a problem is that the optimal (misspecified) model performs better in dense regions of the input space than in sparse regions, because the dense regions dominate the average classification error, which is what we want to minimize. If the dense regions of X are different in the training and test sets, the optimal model on the training set will no longer be optimal on the test set. In other words, the optimal model depends on P (x), and if P (x) = P (x), then the optimal model for the target domain differs from that for the source domain. It was proven that, if the support of P (x) (the set of x for which P (x) > 0) is contained in the support of P (x), then the optimal model that maximizes this re-weighted log likelihood function asymptotically converges to the optimal model for the target domain [START_REF] Shimodaira | Improving predictive inference under covariate shift by weighting the log-likelihood function[END_REF] and a large body of research has been devoted to the estimation of P (x)/P (x) e.g. [START_REF] Zadrozny | Learning and evaluating classifiers under sample selection bias[END_REF], [START_REF] Huang | Correcting sample selection bias by unlabeled data[END_REF], [START_REF] Sugiyama | Direct importance estimation with model selection and its application to covariate shift adaptation[END_REF], [START_REF] Cortes | Learning bounds for importance weighting[END_REF], [START_REF] Bickel | Discriminative learning under covariate shift[END_REF], [START_REF] Kanamori | A least-squares approach to direct importance estimation[END_REF], [START_REF] Kanamori | Statistical analysis of kernel-based least-squares density-ratio estimation[END_REF], [START_REF] Nguyen | Estimating divergence functionals and the likelihood ratio by convex risk minimization[END_REF]. However, reweighting methods do not necessarily improve the prediction accuracy as they also dependent on the extent to which the model is misspecified [START_REF] Wen | Robust learning under uncertain test distributions: Relating covariate shift to model misspecification[END_REF].

In this paper, we show analytically that, despite the fact that the unweighted model is globally more biased than the weighted one, the former may locally be less biased on low importance instances. In view of this result, we design a simple algorithm that combines the weighted and the unweighted models in order to improve the predictive performance in the target domain. More specifically, we prove that an optimal B always exists such that, in the region where β(x) ≤ B , the biased model trained on the unweighted sample should be preferred to the unbiased one, and vice-versa. We propose a practical procedure to estimate this threshold value from training data.

The remainder of this paper is structured as follows. In Section 2, we define some key concepts used along the paper and state some results that will support our analysis. Then in Section 3, we conduct a theoretical analysis to prove that an optimal (but not necessarily unique) B always exists and discuss a manner to optimally combine the weighted and the unweighted models in order to improve the predictive performance in the target domain. In section 4, a series of experiments are carried out on toy problems and real-world data sets to assess the effectiveness of this approach.

Preliminaries

In this section, we define some key concepts used along the paper and state some results that will support our analysis. Consider the supervised learning problem where we observed n training samples, denoted by ((x t ; y t ) : t = 1, ..., n), where x t ∈ X ⊂ R d are i.i.d training input points drawn from some probability distribution p(x) and y t ∈ Y ⊂ R are the corresponding training output values drawn from a conditional probability distribution p(y|x). We are interested in predicting the output value y at an input point x using a model h θ (x) = h(x, θ) parameterized by θ ∈ Θ ⊂ R m . Under covariate shift assumption, the test inputs follow a different probability distribution p (x) while the conditional probability distribution of test output p(y|x) remains unchanged. The ratio β(x) = p (x) p(x) is called the importance of x. Given a loss function l(y, h(x, θ)) : X × Y × Y → [0, ∞), we shall consider throughout this paper, the following loss functions:

-EL-Tr: Expectation of loss over training distribution p(x, y) = p(x)p(y|x) Loss 0 (h θ ) = E x,y∼p [l(y, h(x, θ))] = p(x) p(y|x)l(y, h(x, θ))dydx -EL-Te: Expectation of loss over test distribution p (x, y) = p (x)p(y|x)

Loss 1 (h θ ) = E x,y∼p [l(y, h(x, θ))] = p (x) p(y|x)l(y, h(x, θ))dydx -EL-IWTr: Expectation of Importance-weighted loss over training distribu- tion Loss β (h θ ) = E x,y∼p [β(x)l(y, h(x, θ))]
-B-LEL-Te: We then define Local Expectation of loss over test distribution given β(x) ≤ B of any given hypothesis h θ :

loss(h θ , β(x) ≤ B) = β(x)≤B p (x) Y p(y|x)l(y, h(x, θ))dydx
We also define the optimal parameters of EL-Tr, EL-Te and EL-IWTr:

     θ 0 = argmin θ Loss 0 (h θ ) θ 1 = argmin θ Loss 1 (h θ ) θ β = argmin θ Loss β (h θ ).
It may easily be shown that EL-IWTr is equal to EL-Te,

E x,y∼p [β(x)l(y, h(x, θ))] = p(x) p(y|x) p (x) p(x) l(y, h(x, θ))dydx = p (x) p(y|x)l(y, h(x, θ))dydx
Therefore, minimizing EL-IWTr is equivalent to minimizing EL-Te. Nonetheless, while h θ β is globally less biased than h θ0 , we will show next that it is more biased than h θ0 on low-importance instances. Note that B-LEL-Te can be rewritten as:

loss(h θ , β(x) ≤ B) = β(x)≤B β(x) Y p(x)p(y|x)l(y, h(x, θ))dydx Suppose β(x) takes on continuous value in [b 0 , b M ]
where b 0 > 0, we may rewrite B-LEL-Te as following:

loss(h θ , β(x) ≤ B) = B b0 b β(x)=b Y p(x)p(y|x)l(y, h(x, θ))dydxdb Let L(h θ , β(x) = b) = β(x)=b Y p(x)p(y|x)l(y, h(x, θ))dydx, then: loss(h θ , β(x) ≤ B) = B b0 bL(h θ , β(x) = b)db Similarly, if β(x) takes on discrete values in {b i } M i=0 such that b 0 < b 1 < ... < b M , we rewrite B-LEL-IWTr as: loss(h θ , β(x) ≤ B) = k(B) i=0 b i L(h θ , β(x) = b i )
where k(B) is the largest integer such that b k(B) ≤ B. From the definitions above, we may write

     Loss 1 (h θ ) = loss(h θ , β(x) ≤ b M ), Loss 0 (h θ ) = ∞ b0 L(h θ , β(x) = b)db, for continuous β(x), Loss 0 (h θ ) = M i=0 L(h θ , β(x) = b i ), for discrete β(x).
As aforementioned, a model h(x, θ) is said to be correctly specified if there exist parameter θ * ∈ Θ such that h(x, θ * ) = f (x), otherwise it is said to be misspecified. It is obvious that if a model is correctly specified, the optimal parameter θ of EL-Tr, EL-Te, and any B-LEL-Te coincide. Therefore, the model that minimizes EL-Tr will perform well on the test data globally (i.e., minimizing EL-Te) as well as locally (i.e., B-LEL-Te) in any region of the form β(x) < B. Yet, in practice, almost all models are more or less misspecified. So minimizing EL-Tr θ 0 is not necessarily equivalent minimizing EL-Te. Since EL-Te is equal to EL-IWTr, the parameter minimizing of EL-IWTr θ β , which can be estimated from data, will also minimize EL-Te as shown in [START_REF] Shimodaira | Improving predictive inference under covariate shift by weighting the log-likelihood function[END_REF], [START_REF] Zadrozny | Learning and evaluating classifiers under sample selection bias[END_REF]. However, due to the model misspecification, θ β does not necessarily minimize B-LEL-Te. In fact, we will prove that there exist some

B * (h θ β ) ∈ [b 0 , b M ] such that B-LEL-Te of θ β
exceeds that of θ 0 by proving a stronger conclusion that for all model h θ , with θ ∈ Θ, there exist some B * (h θ ) ∈ [b 0 , b M ] such that B*-LEL-Te of h θ exceeds that of h θ0 , in other words any h θ is locally more biased than h θ0 when predicting instance with β(x) ≤ B * .

In addition, the estimation of θ β may subject to high variance since it involves instance weighting, which is known to reduce the effective samples size [START_REF] Cortes | Learning bounds for importance weighting[END_REF], [START_REF] Gretton | Covariate shift by kernel mean matching[END_REF]. Hence the idea to use h θ0 of instead of h θ β to predict the test instances with β(x) ≤ B * .

Problem analysis

In this section, we conduct theoretical analyses for a simple and then a more general selection bias mechanism. Those analyses will be used to derive a practical procedure aiming at reducing the bias due to covariate shift with misspecified regression or classification learning models.

We first show how EL-Tr is related to B-LEL-Te,

Lemma 1. Suppose β(x) takes on continuous value in [b 0 , b M ] with b M > b 0 > 0, then: Loss 0 (h θ ) = 1 b M loss(h θ , β(x) ≤ b M ) + b M b0 1 B 2 loss(h θ , β(x) ≤ B)dB Proof. For continuous β(x): b M b0 1 B 2 loss(h θ , β(x) ≤ B)dB = b M b0 loss(h θ , β(x) ≤ B)d -1 B = loss(h θ , β(x) ≤ B) -1 B | b M b0 - b M b0 -1 B d(loss(h θ , β(x) ≤ B)) By definition, loss(h θ , β(x) ≤ B) = B b0 bL(b, h θ )db, so loss(h θ , β(x) ≤ b 0 ) = 0 and d(loss(h θ , β(x) ≤ B)) = BL(h θ , β(x) = B)dB. Thus: b M b0 1 B 2 loss(h θ , β(x) ≤ B)dB = -1 b M loss(h θ , β(x) ≤ b M ) + b M b0 1 B (BL(h θ , β(x) = B)dB)
By definition, we have Loss

0 (h θ ) = b M b0 L(h θ , β(x) = B)dB, so: b M b0 1 B 2 loss(h θ , β(x) ≤ B)dB = - 1 b M loss(h θ , b M ) + Loss 0 (h θ )
which concludes the proof A similar results holds in the discrete case.

Corollary 1. Suppose β(x) takes on discrete values {b i } M i=0 such that b 0 < b 1 < ... < b M , then: Loss 0 (h θ ) = 1 b M loss(h θ , β(x) ≤ b M ) + M -1 k=0 1 b k - 1 b k+1 loss(h θ , β(x) ≤ b k ) Proof. M -1 k=0 1 b k - 1 b k + 1 loss(h θ , β(x) ≤ b k ) + 1 b M loss(h θ , β(x) ≤ b M ) = 1 b 0 - 1 b 1 [b 0 L(h θ , β(x) = b 0 )] + 1 b 1 - 1 b 2 [b 0 L(h θ , β(x) = b 0 ) + b 1 L(h θ , β(x) = b 1 )] + ... + 1 b M -1 - 1 b M [b 0 L(h θ , β(x) = b 0 ) + ... + b M -1 L(h θ , β(x) = b M -1 )] + 1 b M [b 0 L(h θ , β(x) = b 0 ) + b 1 L(h θ , β(x) = b 1 ) + .. + b M L(h θ , β(x) = b M )] = b 0 L(h θ , β(x) = b 0 ) 1 b 0 - 1 b 1 + 1 b 1 - 1 b 2 + ... + 1 b M -1 - 1 b M + 1 b M + ... + b M -1 L(h θ , β(x) = b M -1 ) 1 b M -1 - 1 b M + 1 b M + b M L(h θ , β(x) = b M ) 1 b M = M i=0 L(h θ , β(x) = b i ) = Loss 0 (h θ )
In view of Corollary 1, we may now state the following theorem, Theorem 1. Suppose there exists two real values, b 0 and b 1 , such that b 0 < 1 < b 1 and a subset X 0 ⊂ X such that

β(x) = b 0 if x ∈ X 0 b 1 if x / ∈ X 0 ,
then there exists a threshold B * such that:

loss(h θ1 , β(x) ≤ B * ) ≥ loss 1 (h θ0 , β(x) ≤ B * ).
In fact, B * can take any value in [b 0 , b 1 ).

Proof. By definition, Loss 0 (h θ0 ) ≤ Loss 0 (h θ1 ), using Lemma 1, we may write: When the assumptions of Theorem 1 holds, we say that the covariate shift scheme follows a simple step distribution. The equality in Theorem 1 only occurs when θ 0 minimizes EL-Te and θ 1 minimizes EL-Tr. Such condition indicates that covariate shift does not have an effect on searching for optimal θ, which is a rare case as shown by other studies. Theorem 1 shows that for simple step distribution where inclusion in the training sample is either proportional to b -1 0 (over-sampled instances), or to b -1 1 (under-sampled instances), h θ0 exhibits a lower bias compared to h θ1 on the low importance test instances. This type of selection bias mechanism is actually quite common. For instance, prospective cohort studies in epidemiology are by design prone to covariate shift because selection criteria are associated with the exposure to potential risk factors.

Loss 0 (h θ0 ) = 1 b 1 loss(h θ0 , β(x) ≤ b 1 ) + 1 b 0 - 1 b 1 loss(h θ0 , β(x) ≤ b 0 ) = 1 b 1 Loss 1 (h θ0 ) + 1 b 0 - 1 b 1 loss(h θ0 , β (
Theorem 2. For all θ ∈ Θ, there exists a threshold B * (h θ ) such that

loss(h θ , β(x) ≤ B * (h θ )) ≥ loss(h θ0 , β(x) ≤ B * (h θ )) (1)
B * (h θ ) could take any value in the set below:

B * (h θ ) = argmax B (loss(h θ , β(x) ≤ B) -loss(h θ0 , β(x) ≤ B))
The equality occurs whenever θ 1 is also a minimum for EL-Tr.

Proof. We prove by contradiction that Theorem 2 holds. Assume that inequality 1 does not hold for B * (h θ ) defined above:

loss(h θ , β(x) ≤ B * (h θ )) -loss(h θ0 , β(x) ≤ B * (h θ )) < 0 (2)
By definition of B * (h θ ), we may show that, for all

B ∈ [b 0 , b M ], loss(h θ , β(x) ≤ B) -loss(h θ0 , β(x) ≤ B) < 0 Thus, for all B ∈ [b 0 , b M ] loss(h θ0 , β(x) ≤ B) > loss(h θ , β(x) ≤ B)
Now, using Lemma 1 for continuous β(x), we have:

Loss 0 (h θ0 ) = 1 b M loss(h θ0 , β(x) ≤ b M ) + b M b0 1 B 2 loss(h θ0 , β(x) ≤ B)dB > 1 b M loss(h θ , β(x) ≤ b M ) + b M b0 1 B 2 loss(h θ , β(x) ≤ B)dB = Loss 0 (h θ )
Hence, Loss 0 (h θ0 ) > Loss 0 (h θ ), contradicts the fact that θ 0 = argmin θ Loss 0 (h θ ) is the optimal hypothesis under the unweighting scheme and θ = argmin θ Loss 0 (h θ ).

If the two terms in inequality1 are equal, then we can prove similarly that Loss 0 (h θ0 ) = Loss 0 (h θ ), which implies that θ 1 is also a minimal solution of EL-Tr. The demonstration for discrete β(x) values follows similarly.

Theorem 2 states that any model h θ with θ ∈ Θ is outperformed by h θ0 learned from the unweighted training samples in terms of bias when predicting examples with β(x) ≤ B * (h θ ). This is also applied to model h θ β which minimizes EL-IWTr. In addition, the estimation of θ β may exhibit a higher variance due to the effective sample size reduction as discussed in [START_REF] Cortes | Learning bounds for importance weighting[END_REF][START_REF] Gretton | Covariate shift by kernel mean matching[END_REF]. These results altogether suggest that h θ0 should be preferred to h θ β for predicting the instance's outputs in the region β(x) ≤ B * (h θ ), termed low-importance region. Therefore, for any learning task with covariate shift, we shall train two distinct models, one with and the other without the importance weighting scheme. Then, we shall use the latter to predict instances satisfying β(x) ≤ B * (h θ ) and use the former to predict the remaining instances. The optimal value for B * (h θ ) may be estimated from the training data. The set of all possible empirical threshold B * (h θ β ) can be obtained empirically by solving the following problem :

B * (h θ ) = argmax B 1 n i∈{1,..,n} β(xi)≤B β(x i )[l(y i , h(x i , θ β )) -l(y i , h(y i , θ 0 ))] (3)
As n grows to infinity, it follows from the law of large numbers that,

B * (h θ ) → B * (h θ )
Therefore, B * (h θ β ) could be estimated empirically either from training data or by cross validation. In this study, we use a 5-fold importance weighted cross validation to estimate B * (h θ β ) as suggested in [START_REF] Sugiyama | Direct importance estimation with model selection and its application to covariate shift adaptation[END_REF]. It should be emphasized that B * (h θ β ) is not necessarily unique. For instance, any value between b 0 and b 1 in Theorem 1 is admissible as mentioned earlier.

Experiments

In this section, we assess the ability of our "hybrid approach" to reduce the learning bias under covariate shift based on Theorem 2. We first discuss the strategies employed to estimate the importance weights: one is based explicitly on the true bias mechanism, the other is based on linear density-ratio model.

We emphasize that the latter does not require any prior knowledge of the true sampling probabilities to estimate the β(x) values, and uses the test input features instead. In fact, the estimation of distribution is a hard problem, thus it is more appealing to directly estimate β(x). Indeed, a large body of work has been devoted to this line of research e.g. [START_REF] Zadrozny | Learning and evaluating classifiers under sample selection bias[END_REF], [START_REF] Huang | Correcting sample selection bias by unlabeled data[END_REF], [START_REF] Sugiyama | Direct importance estimation with model selection and its application to covariate shift adaptation[END_REF], [START_REF] Nguyen | Estimating divergence functionals and the likelihood ratio by convex risk minimization[END_REF], [START_REF] Cortes | Learning bounds for importance weighting[END_REF], [START_REF] Bickel | Discriminative learning under covariate shift[END_REF], [START_REF] Kanamori | A least-squares approach to direct importance estimation[END_REF]. From the many references, we choose the Unconstrained Least-Square Importance Fitting (uLSIF) estimator for β(x) that was proved to be successful with covariate shift. We then study a toy regression problem to show if covariate shift corrections based on our method reduces prediction error on the test set when the learning model is misspecified. We then test our approach on real world benchmark data sets, from which the training examples are selected according to various biased sampling schemes as suggested in [START_REF] Kanamori | A least-squares approach to direct importance estimation[END_REF].

Importance ratio estimation

As aforementioned, we use two weighting schemes in ours experiments, one is derived from the true selection bias mechanism an one is Unconstrained Least-Square Importance Fitting (uLSIF), a method based on linear density-ratio models [START_REF] Kanamori | A least-squares approach to direct importance estimation[END_REF]. Formally, it assumes that the density ratio β(x) can be approximated by a linear model β(x) = M i=1 α i h i (x) where the basis functions h i , i = 1, ..., M are chosen so that h i (x) ≥ 0 for all input value x. The coefficients α 1 , ..., α M are parameters of the linear model and are estimated from data by minimizing the empirical square error between weighted biased distribution (from training data) and the bias-free distribution of x:

min α 1 2n n i=1 ( β(x i )) 2 - 1 n n i=1 β(x i ) + λ.Reg(α)
where {x i } n i=1 and {x i } n i=1 , are the training and test inputs, Reg(α) is the regularization term, introduced to avoid overfitting. A heuristic choice of h i (x) proposed in [START_REF] Kanamori | A least-squares approach to direct importance estimation[END_REF] is a Gaussian kernel centered at the test points {x i } n i=1 when the number of test points is small (less than 100) or at template points {x i } 100 i=1 , which is a random subset of test set when the number of test points is large for computation advantage. The kernel width and the regularization term Reg(α) are optimized by cross-validation with grid search.

Toy regression problem

Consider the following training data generating process: x ∼ N (µ 0 , σ 0 ) and y = f (x) + , where µ 0 = 0, σ 0 = 0.5, f (x) = -x + x 3 , and ∼ N (0, 0.3). In the test data, we have the same relationship between x and y but the distribution of the covariate x is shifted to x ∼ N (µ 1 , σ 1 ), where µ 1 = 0, σ 1 = 1. The training and test distributions, along with their ratio are depicted in Fig. 1a and1b. The minimization of EL-Tr is obtained using the unweighted Least Square Regression (uLSR) method for the normal regression while minimization of EL-Te is performed by the weighted Least Square Regression (wLSR). As shown in (b) True importance weights q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -2 -1 0 Fig. 1: An illustrative example of fitting a function f(x) using a linear model with/without the weight importance scheme (wLSR/uLSR) and a combination of both (termed "Mix").

[10], wLSR is unbiased thus it should perform better than uLSR, which is biased, on test data. However, as can be seen in Fig. 1c, uLSR (red dashed line) seems to better approximate the y = f (x) curve (in blue) than wLSR (black dashed line) on instances in the interval (-1, 1). As may be seen in Fig. 1d, the hybrid model that optimally combines wLSR and uLSR, based on Theorem 1, achieves a lower Mean Square Error (MSE) compared to wLSR. The experiment was repeated 30 times for each number of sample size. It should be noted that the hybrid model always outperforms the weighted model and the gain in performance on the test set is more noticeable for larger training sizes.

Simple step sample selection distribution

In this second experiment, we consider a simple step distribution with known or estimated selection probabilities and we apply this selection scheme on a variety of UCI data sets in order to assess the efficiency of our bias correction procedure in more realistic scenarios. We use a SVM classifier for both classification and regression tasks. Experiments are repeated 100 times for each data set. In each trial, we randomly select an input feature x c to control the bias along with 300 training samples. We then apply the following single step probability distribution as discussed in Theorem 1,

P (s = 1|x = x c i ) = p s = p1 = 0.9 if x c i ≤ mean(x c ) p2 = 0.9 1+exp(r) otherwise
where r is a parameter that controls the strength of the selection bias. In each trial r takes a random value from a normal distribution N (2, 0.1). With these parameters, the selection probability for instances having an x c value (e.g. a degree of exposure to some risk factor) above the mean is between 7 to 10 times smaller than for those having of a lower value. This is a scenario that typically arises in epidemiological cohort studies when subjects are includes in the study according to some exposure factor. Consider the two following weighting schemes. The first one: β = p (x)/p(x) = p(s = 1)/p(s = 1|x) ∼ 1/p s assumes that the bias mechanism is known exactly.

β(x) ∼ p -1 s ∼ b1 = 1 if x c i ≤ mean(x c ) b2 = 1 + exp(r)otherwise
In practice, however, the selection probability is rarely known exactly. So let us assume that the estimation of β is subject to some error and let us consider the following approximate weighting scheme:

β(x) ∼ p -1 s ∼ b1 = 1 if x c i ≤ mean(x c ) b2 = 1 + exp(r) if otherwise
where r = r + N (0, 0.1) is our noisy estimate of r. For each weighting scheme, we fit a true weighted model (denoted as P in Table 1) and an approximated weighted model (denoted as P in Table 1). As p 1 < 1 and p 2 > 1, our weighting mechanism satisfies the assumptions of Theorem 1, so we set B * = 1. We report the mean square errors (MSE) in Tab.1. All values are normalized by the MSE of the unweighted model (our gold standard). As may be seen from the plots in Fig. 2a and2b, the combined models outperform the weighted ones. That is, when using either exact probability ratio, the results obtained with P mix are better than that of P . The same observation can be made when the estimated probability ratios are used instead (i.e., Pmix versus P ) and except on the Banknote data set. The gain is significant at the significance level 5% using the Wilcoxon signed rank test.

General covariate selection mechanisms

In this last experiment, we use the same setting as above but we use a more general distribution:

P (s = 1|x = x c i ) = ps =      p1 = 0.9 if x c i ≤ mean(x c ) p2 = 0.1 if x c i > mean(x c ) + 0.8 × 2σ(x c ) p3 = 0.9 - x c i -mean(x c ) 2σ(x c )
otherwise.

where σ(x c ) denotes the standard deviation of x c . As may be observed, the assumptions required in Theorem 1 do no hold anymore with this more general sample selection distribution. According to Eq.3, we need to estimate B * (h θ ) empirically from data. We consider again two importance weighting schemes: one is based on the true underlying probability and is referred to as P , while the other is based on the uLSIF estimator. As may be observed from Table 2 and Figures 2c and2d that performances of the hybrid models are significantly improved with respect to the weighted models, except with the Congressional Voting and Banknote data sets. q q q q q q q q q q q q q q 0.80 0.90 1.00 1.10 0.80 0.95 1.10 weighted mix (a) q q q q q q q q q q q q q q 0.80 0.90 1.00 1.10 0.80 0.95 1.10 weighted mix (b) q q q q q q q q q q q q q 0.96 1.00 1.04 0.96 1.02 weighted mix (c) q q q q q q q q q q q q q q 0.96

Conclusions

In this paper, we showed that the standard importance weighting approach used to reduce the bias due to covariate shift can easily be improved when misspecified Even for more general covariate shift scenarios, we proved that there always exist a threshold for the importance weight below which the test instances should be predicted by the globally biased model. In view of this result, we proposed a practical procedure to estimate this threshold and we discussed a simple procedure to combine the weighted and unweighted prediction models. The method was shown to be effective in reducing the bias on several UCI data sets.

  x) ≤ b 0 ) Similarly, Loss 0 (h θ1 ) = It is easily shown that the right hand side of inequality above is non-negative due to the definition of θ 1 . It follows that loss(h θ1 , β(x) ≤ b 0 ) -loss(h θ0 , β(x) ≤ b 0 ) ≤ 0 which, given the assumption about β(x), is equivalent to, loss(h θ1 , β(x) = b 0 ) -loss(h θ0 , β(x) = b 0 ) ≤ 0 Thus the Theorem is true when B * = b 0 . It is also true for any other B * ∈ [b 0 , b 1 ) as a consequence.

  examples High-importance examples (c) True function, uSLR and wSLR on test data.

Fig. 2 :

 2 Fig. 2: MSE gain of the mix model vs. MSE gain of the weighted model. Points below the diagonal line indicate that the mix model outperforms the weighted model. Figures (a) and (b): simple step distribution covariate shift used in the first experiment with the weighted model based on (a) the true selection probability and (b) based on the estimated selection probability. Figures (c) and (d): covariate shift in used in the second experiment when the weighted model based was based on: (c) true selection probability; (d) on uLSIF.

Table 1 :

 1 Mean test error averaged over 200 trials with different weighting schemes on 15 UCI data sets. Data sets marked with '*' are regression problems. P denotes the weighting scheme using the true selection probability and P denotes the weighting scheme using a noisy selection probability. For each pair of weighted and mix models, the better prediction value is highlighted in boldface

	Data set	No weighting P	P mix	P	P mix
	India diabetes	1.000 ± 0.020 0.966 ± 0.019 0.960 ± 0.018 0.968 ± 0.019 0.962 ± 0.018
	Ionosphere	1.000 ± 0.128 0.915 ± 0.105 0.902 ± 0.107 0.911 ± 0.104 0.897 ± 0.106
	BreastCancer	1.000 ± 0.039 1.020 ± 0.044 1.013 ± 0.044 1.020 ± 0.044 1.013 ± 0.043
	GermanCredit	1.000 ± 0.008 1.000 ± 0.007 0.996 ± 0.008 1.000 ± 0.008 0.996 ± 0.008
	Australian credit	1.000 ± 0.006 0.963 ± 0.008 0.947 ± 0.010 0.964 ± 0.008 0.947 ± 0.010
	Mushroom	1.000 ± 0.068 0.090 ± 0.057 0.872 ± 0.060 0.888 ± 0.058 0.874 ± 0.056
	Congressional Voting 1.000 ± 0.033 1.026 ± 0.039 0.993 ± 0.038 1.030 ± 0.038 1.000 ± 0.037
	Banknote	1.000 ± 0.040 0.970 ± 0.043 0.978 ± 0.038 0.969 ± 0.042 0.975 ± 0.039
	Airfoil self noise*	1.000 ± 0.023 0.997 ± 0.015 0.961 ± 0.012 0.993 ± 0.015 0.958 ± 0.012
	Abanlone*	1.000 ± 0.032 0.984 ± 0.020 0.960 ± 0.020 0.985 ± 0.021 0.961 ± 0.020
	Auto MGP*	1.000 ± 0.084 0.939 ± 0.066 0.933 ± 0.067 0.939 ± 0.066 0.930 ± 0.067
	Boston Housing*	1.000 ± 0.057 1.037 ± 0.053 0.994 ± 0.050 1.037 ± 0.053 0.994 ± 0.050
	Space GA*	1.000 ± 0.009 1.021 ± 0.007 0.962 ± 0.008 1.018 ± 0.008 0.961 ± 0.008
	Cadata*	1.000 ± 0.013 1.038 ± 0.022 1.029 ± 0.017 1.037 ± 0.022 1.029 ± 0.017

Table 2 :

 2 Mean test error averaged over 200 trials for different weighting schemes on UCI data set. Data sets marked with * are for regression problems. P denotes the weighting scheme based on the true selection probability and uLSIF denotes the weighting scheme using the uLSIF estimator. For each pair of weighted and mix models, the better prediction value is highlighted in boldface. ± 0.021 0.980 ± 0.018 0.975 ± 0.018 1.016 ± 0.021 1.006 ± 0.021 Ionosphere 1.000 ± 0.087 1.006 ± 0.087 0.988 ± 0.085 1.028 ± 0.093 1.007 ± 0.087 BreastCancer 1.000 ± 0.019 1.004 ± 0.018 0.993 ± 0.019 1.000 ± 0.018 0.993 ± 0.019 GermanCredit 1.000 ± 0.008 1.003 ± 0.008 0.999 ± 0.008 1.009 ± 0.008 1.001 ± 0.008 Australian credit 1.000 ± 0.009 0.972 ± 0.007 0.967 ± 0.007 1.007 ± 0.008 1.005 ± 0.008 Mushroom 1.000 ± 0.558 1.011 ± 0.054 0.963 ± 0.051 0.991 ± 0.054 0.989 ± 0.054 Congressional Voting 1.000 ± 0.037 1.023 ± 0.036 1.010 ± 0.037 0.987 ± 0.036 0.997 ± 0.036 Banknote 1.000 ± 0.060 1.083 ± 0.057 0.962 ± 0.062 0.962 ± 0.061 0.979 ± 0.058 Airfoil self noise* 1.000 ± 0.007 0.995 ± 0.007 0.995 ± 0.007 1.011 ± 0.008 1.001 ± 0.008 Abanlone* 1.000 ± 0.007 1.001 ± 0.008 1.001 ± 0.007 1.005 ± 0.007 0.998 ± 0.006 Auto MGP* 1.000 ± 0.026 0.990 ± 0.025 0.970 ± 0.025 1.015 ± 0.027 0.994 ± 0.026 Boston Housing* 1.000 ± 0.043 0.984 ± 0.031 0.940 ± 0.032 1.036 ± 0.040 0.989 ± 0.042 Space GA* 1.000 ± 0.006 1.005 ± 0.005 0.980 ± 0.006 1.000 ± 0.005 0.996 ± 0.005 Cadata* 1.000 ± 0.012 1.008 ± 0.013 1.006 ± 0.012 1.023 ± 0.013 1.010 ± 0.012 training models are used. Considering a simple class of selection bias mechanisms, we proved analytically that the unweighted model exhibits a lower prediction bias compared to the globally unbiased model in the low importance input subspace.

	Data set	No weighting P	P mix	uLSIF	uLSIF mix
	India diabetes	1.000

Acknowledgments: This work was partially supported by a grant from the European ENIAC Joint Undertaking (INTEGRATE project).