Student Member, IEEE C Nsiala Nzéza

R Gautier

Member, IEEE G Burel

Blind Multiuser Identification in Multirate CDMA Transmissions: A New Approach

A new blind sequences synchronization method in multiuser multirate CDMA context is proposed in this paper. Instead of the approach based on the FROBENIUS Square Norm Behaviour (FSNB) previously proposed in [1]-[3], we develop a new blind method based on the Maximum Eigenvalue Behaviour (MEVB) according to analysis window shifts. Theoretical analysis shows that the MEVB-based criterion provides a significant improvement of performances, and proves that it is a powerful tool for blind synchronization. We show that the improvement is mainly due to suppression of synchronisation peaks masking that occurred with the previous method. Simulation results confirm that, performing the synchronization process using this new criterion, allows one to achieve very good performances at the receiver side in term of chip error rate (CER) and bit error rate (BER), even at very low SNRs.

I. INTRODUCTION

Spread spectrum signals have been used in the military domain for long time for secure communications [START_REF] Picholtz | Theory of spreadspectrum communications -a tutorial[END_REF]. Nowadays their field of application includes civilian transmissions, especially CDMA transmissions [START_REF] Magill | Spread-spectrum tehcnology for comercial applications[END_REF]. Thanks to the properties of the pseudo-random sequences used, the CDMA technique allows to solve the problem of the increasing number of users in a frequency band. Moreover, these signals are difficult to detect, especially in a non cooperative context (e.g. spectrum surveillance), because they are often below the noise level, due to low Signal-to-Noise Ratios (SNRs).

Several blind approaches (i.e. when the process of recovering data from multiple simultaneously transmitting users without access to training sequences) have been addressed in the literature [START_REF] Wang | Blind multiuser detection: A subspace approach[END_REF]- [START_REF] Buzzi | Iterative cyclic subspace tracking for blind adaptative multiuser detection in multirate cdma systems[END_REF]. In another way of research, many partially blind multiuser schemes, that exploit some known channel properties, have been proposed [START_REF] Wang | Blind adaptive interference suppression for CDMA communication based on eigenspace tracking[END_REF], [START_REF] Honig | Blind estimation of the pseudo-random sequence of a direct spread spectrum signal[END_REF]. But, all these methods are not blind in the sense used in our non cooperative context. Indeed, they require some prior knowledge about users parameters. This knowledge is not available in a non cooperative context.

Hence, this article is organized as follows: Section II will introduce the signal model and assumptions made. Section III will present the new approach. Section IV will deal with the theoretical analysis of the MEVB-based criterion and compare it with the FSNB-based one. The simulations results will be detailed in Section VI, and our conclusions will be drawn in Section VII.

II. SIGNAL MODELING AND ASSUMPTIONS

Let us consider the uplink multirate CDMA transmissions using the variable spreading length (VSL) technique (i.e. sequences have the same chip period, and data rates are tied to sequences length), the downlink can be viewed as a particular case of the uplink. Let us denote S, a set of available data rates

R 0 < R 1 < • • • < R S-1 . By denoting N i
u the number of active users transmitting at R i and N u the total number of users such that S-1 i=0 N i u = N u , the received signal can be expressed as:

y(t) = S-1 i=0 N i u -1 n=0 +∞ k=-∞ a n,i (k)h n,i (t -kT si -T dn,i ) + b(t) (1) where h n,i (t) = Li-1 k=0 c n,i (k)p i (t -kT c ).
In [START_REF] Burel | Blind estimation of the pseudo-random sequence of a direct spread spectrum signal[END_REF], the subscript (•) n,i refers to the n th user transmitting at R i , denoted throughout this article as the (n, i) th user. Accordingly:

• a n,i (k) are the baseband symbols of variance σ 2 an,i for the (n, i) th user, whereas p i (t) is the convolution of the transmission filter, channel filter (which takes into account channel echoes, fading, multipaths and jammers) and receiver filter for each rate.

• The term h n,i (t), defined in [0 T si [, is a virtual filter corresponding to the convolution of all filters of the transmission chain with the spreading sequence

{c n,i (k)} k=0•••Li-1
, where L i is the spreading factor for the (n, i) th user.

• Because of the VSL technique, the symbol period T si for the users transmitting at the rate R i is tied to the common chip period T c : T si = L i T c , and s n,i stands for the (n, i) th signal.

• The term T dn,i is the corresponding transmission delay for the (n, i) th user; it is assumed to satisfy: 0 T dn,i < T si and to remain constant during the observation.

• b(t) is a centered white Gaussian noise of variance σ 2 b . • Signals are assumed to be independent, centered, noiseuncorrelated and received with the same power: σ σ 2 s0,0 , for all (n, i), n = 0, 1,

• • • , N i u -1, i = 0, 1, • • • , S -1.
• Finally, the SN R (in dB) at the detector intput is negative (signal hidden in the noise).

III. PROPOSED APPROACH

The blind multiusers detection scheme proposed in [START_REF] Nsiala Nzéza | Blind multiuser detection in multirate cdma transmissions using fluctuations of correlation estimators[END_REF] allowed us to estimate symbol periods T si , thus data rates. Then, the synchronization process is performed within each group i of users transmitting at the same data rate. The received signal is sampled and divided into N non-overlapping temporal windows of duration T F = T si = M T e , M ∈ N ⋆ , where T e is the sampling period, T si is the symbol period estimated in [START_REF] Nsiala Nzéza | Blind multiuser detection in multirate cdma transmissions using fluctuations of correlation estimators[END_REF], and denoted within a group T s for more clearness.

Denoting Y the sampled received signal matrix which each column contains M signal y(t) samples, leads us to get the following (M × N )-matrix:

Y =    y(t) • • • y(t + (N -1)T s ) . . . • • • . . . y(t + T s -T e ) • • • y(t + N T s -T e )    (2)
Thanks to theoretical analysis, (2) can be expressed as:

Y = N i u -1 n=0 h 0 n + h -1 n a T n (3)
where the vector a

T n = [• • • , a n (m) • • • ]
contains all the symbols of the n th user. Since filters h n (t) are defined in [0 T s [ vectors h n 0 and h -1 n are defined, for each interfering user whithin the group, as follows :

• the vector h -1 n contains the end of the corresponding spreading waveform during T st n , followed by zeros during t n ;

• the vector h 0 n contains zeros during t n , followed by the beginning of the corresponding spreading waveform during T st n . where t n , n = 0, 1, • • • , N i u -1, represents the temporal shift between an analysis window and the beginning of a whole symbol of each corresponding signal, as illustrated on Fig. 1 for two users, T d represents the transmission delay of the second signal, and D f the analysis window translation.

Vectors h 0 n and h -1 n allow one to take into account the temporal shifts t n between an analysis window and the beginning of a whole symbol of the n th interfering signal in the correlation matrix defined as: R = YY * , where (•) * stands for the transpose conjugate of (•). Then, like in [START_REF] Burel | Blind estimation of the pseudo-random sequence of a direct spread spectrum signal[END_REF]- [START_REF] Nsiala Nzéza | Blind synchronization and sequences identification in CDMA transmissions[END_REF], and under the assumptions of independent, centered and noiseuncorrelated signals, and since the signals are assumed to be received with the same power, the correlation matrix R can be written as:

R = σ 2 b    ρ T s T e N i u -1 n=0 (1 - t n T s )v 0 n (v 0 n ) * + t n T s v -1 n (v -1 n ) * + I ( 4
)
where ρ is the signal to noise and interference ratio (SNIR) compared to one of the other users, I is the identity matrix, and vectors v 0 n and v -1 n are normalized versions of h 0 n and h -1 n . Then, we set:

β = ρ Ts Te , d f = D f Ts α n = tn Ts , τ n = T dn Ts (5) Thus, α n , τ n , d f ∈ [0 1[. The (M × M ) correlation matrix R can be rewritten as: R = σ 2 b   β N i u -1 n=0 (1 -α n )v 0 n (v 0 n ) * + α n v -1 n (v -1 n ) * + I   (6)
Since sequences are assumed to be slightly correlated, the matrix R eigenvalues in an unspecified order can be expressed as:

     λ 0 n = σ 2 b {β (1 -α n ) + 1} , n = 0, • • • , N i u -1 λ -1 n = σ 2 b {βα n + 1} , n = 0, • • • , N i u -1 λ n = σ 2 b , n = 2N i u , • • • , M -1 (7)
Equations ( 6) and ( 7) highlight 2N i u eigenvalues associated to the signal space and the M -2N i u ones associated to the noise space (all are assumed to be equal on average to the noise power). Fig. 1 and (7) evidence that each value α n will change according to the normalized analysis windows shifts d f in cyclic fashion. Thus, to take it into account, we set:

α = d f -τ n , n = 0, • • • , N i u -1 x ≡ x modulo 1, x ∈ R (8)
Then, using ( 6), ( 7) and ( 8) led us to derive a blind synchronization criterion based on the FROBENIUS Square Norm behaviour (FSNB) of ( 6) according to analysis windows shifts [START_REF] Burel | Blind estimation of the pseudo-random sequence of a direct spread spectrum signal[END_REF]- [START_REF] Nsiala Nzéza | Blind synchronization and sequences identification in CDMA transmissions[END_REF]. Nevertheless, theoretical analysis of the extended FSNB-based criterion to the multiuser case, evidence 40th IEEE -Asilomar Conference on Signals, Systems and Computers, Oct. 29 -Nov. 1, 2006, pp. 2162-2166, Pacific Grove, California ©2006 IEEE. Personal use of this material is permitted. However, permission to reprint/republish this material for advertising or promotional purposes or for creating new collective works for resale or redistribution to servers or lists, or to reuse any copyrighted component of this work in other works must be obtained from the IEEE. the performances degradation further to any increase in the number of users and/or decrease of SNRs. Consequently, we developed a new blind synchronization scheme based no longer on the study of the FROBENIUS square norm behaviour, but rather on each maximum eigennvalue behaviour (MEVB) with analysis window shifts, as detailed in section below.

IV. THEORETICAL ANALYSIS

The FROBENIUS Square Norm of ( 6), denoted R 2 F is defined as:

R 2 F = N i u -1 n=0 (λ 0 n ) 2 + (λ -1 n ) 2 + (M -N i u )σ 2 b (9)
Using ( 8) in ( 9), and only keeping the variable part of the result, allowed us to express the extended FSNB-based criterion previously proposed in [START_REF] Burel | Blind estimation of the pseudo-random sequence of a direct spread spectrum signal[END_REF], and denoted here F :

F (d f ) = 1 + N i u -1 n=0 { d f -τ n 2 -d f -τ n } (10)
Equation ( 10) is quadratic convex function in each interval [τ n τ n+1 [. The extended FSNB-based synchronization scheme consists in maximizing [START_REF] Wang | Blind adaptive interference suppression for CDMA communication based on eigenspace tracking[END_REF], which is equivalent to find the criterion (10) maxima (synchronization peaks). Let us recall that a peak is a curve point from which while moving by lower or higher values, the curve is always decreasing. From ( 7) and ( 8), synchronization peaks may be at points d f = τ n . However, we proved that according to the choice of τ n , e.g. when two consecutive τ n and τ n+1 are very close, points at d f = τ n are not synchronization peaks, as illustrated on Fig. 2, where: τ 0 = 0, τ 1 = 0.0448, τ 2 = 0.068, τ 3 = 0.3853.

Hence, this considerations led us to propose, as detailed here after, a new blind synchronization method based on the maximum largest eigenvalue behaviour according to analysis window shifts (MEVB).

Without loss of generality, let us set: σ 2 b = β = 1, and introduce (8) in [START_REF] Haghighat | A subspace scheme for blind user identification in multiuser ds-cdma[END_REF]. Thus, the consecutive largest eigenvalues of ( 6) can be rewritten as:

(a)

λ 0 n = 2 -d f -τn λ -1 n = 1 + d f -τn (b) λ 0 n+1 = 2 -d f -τn+1 λ -1 n+1 = 1 + d f -τn+1 (11
) For any analysis window shift d f such that τ n ≤ d f < τ n+1 and for any user n, (11) (a) shows that λ 0 n decreases while λ -1 n increases, and conversly. Moreover, (11) (a) and (b) also show that the n th user eigenvalue decreases starting from its maximum value (synchronization peak) in d f = τ n , while that of the (n + 1) th increases, and the values for which they are equal are local minima. Hence, the MEVB-based criterion is defined as the maximum value between two consecutive largest eigenvalues. Consequently, the MEVB-based criterion is equivalent to the function C(d f ) defined as:

C(d f ) = max τn≤d f <τn+1 λ 0 n , λ -1 n+1 , n = 0, • • • , N i u -1 (12)
Since eigenvalues are linear functions of d f , the MEVBbased criterion (12) always presents maxima, thus synchronization peaks in d f = τn , contrary to the extended FSNBbased one, as shown on Fig. 3, with the same parameters that those on Fig. 2. Let us note that the MEVB-based criterion can easily be derived for downlink transmissions by setting τ n = 0, n = 0, • • • , N i u -1 in all equations above. Hence, the discussion above highlights that, after synchronizing using the criterion (12) performances will be better than those that will be obtained after synchronizing using the criterion [START_REF] Honig | Blind estimation of the pseudo-random sequence of a direct spread spectrum signal[END_REF], due to the lack of synchronization peaks masking as it will be shown in section VI. Once the synchronization process has been performed, sequences can be identified and the numbers of interfering users within the group detreminated as detailed in the following section.

V. SEQUENCES IDENTIFICATION AND NUMBER OF USERS

DETERMINATION

Since it was set σ 2 b = β = 1, once one of the interfering users is synchronised, e.g. d f = τ 0 , which corresponds to set α 0 = 0, the correlation matrix becomes:

40th IEEE -Asilomar Conference on Signals, Systems and Computers, Oct. 29 -Nov. 1, 2006, pp. 2162-2166, Pacific Grove, California ©2006 IEEE. Personal use of this material is permitted. However, permission to reprint/republish this material for advertising or promotional purposes or for creating new collective works for resale or redistribution to servers or lists, or to reuse any copyrighted component of this work in other works must be obtained from the IEEE.

R =v 0 v * 0 + N i u -1 n=1 (1 -α n )v 0 n (v 0 n ) * + α n v -1 n (v -1 n ) * + I (13)
Hence, still under the assumptions of almost uncorrelated sequences, (13) highlights a maximum eigenvalue which associated vector contains the corresponding spreading sequence (apart from the effects the global transmission filter), 2(N i u -1) eigenvalues generated by the other users, and M -2N i u + 1 eigenvalues equal, on average, to the noise power. Then this process is performed in an iterative way so as to get the N i u largest eigenvalues, which the N i u associated vectors contain the corresponding spreading sequences.

Moreover, the number of interfering users N i u within a group is directly equal to the number of synchronization peaks. In addition, in the particular case of downlink, i.e. τ n = 0, α n = α, n = 0, 1, • • • , N i u -1, when users are synchronized, i.e. α = 0, the correlation matrix becomes:

R = N i u -1 n=0 v n v * n + I (14)
Hence, the correlation matrix (14) has N i u largest eigenvalues and M -N i u eigenvalues equal which associated vector contains the corresponding spreading sequence (apart from the effects the total transmission filter), and the others are on average equal to the noise power. In this case, the number of interfering users is equal to the number of the largest eigenvalues.

At last, linear algebra techniques described in [START_REF] Burel | Blind estimation of the pseudo-random sequence of a direct spread spectrum signal[END_REF]- [START_REF] Nsiala Nzéza | Blind synchronization and sequences identification in CDMA transmissions[END_REF] and applied to the estimated eigenvectors, allow one to identify sequences used at the transmiter side and to recover transmitted symbols. It is obvious that the synchronization process errors will be propagated to the sequences identification and symbols recovering process. Hence, Section VI will focuss on performances in term of mean chip error rate (MCER) and mean bit error rate (MBER) at the receiver after synchronizing using both MEVB and FSNB criteria. We did not compare both criteria with others, since all schemes adressed in the literature need prior knowledge of at least one transmission parameter.

VI. SIMULATION RESULTS

Simulations were carried out in uplink with 4 signals of 170µs, each of them spread by a complex GOLD sequence of length L = 127. The chip frequence was F c = 100 MHz, the initial sampling frequence was F e = 300 MHz, the SNR was -5dB at the detector input, i.e., only the MAI noise is considered. The initial number of windows was N = 58, with a duration of 550µs, and the number of samples was 32768. The symbols belong to a QP SK constellation. The symbol period is blindly estimed as described in [START_REF] Nsiala Nzéza | Blind multiuser detection in multirate cdma transmissions using fluctuations of correlation estimators[END_REF] and was T s = 1.27µs. So, it was set: T F = T s , and normalized arbitrary analysis windows shifts at the beginning of the process was set: α 0 = 0.3228, α 1 = 0.4226, α 2 = 0.7533, α 3 = 0.9423.

Fig. 4 illustrates the experimental MEVB criterion curve and evidences 4 synchronization peaks. Hence, there are 4 users transmitting at this data rate, as expected.

Fig. 5 shows the probability of synchronizing at least one user, i.e. the probability of detecting and estimating at least a synchronization peak. It cleary confirms that this probability is higher when using the MEVB-based criterion than that when using the FNSB-based one.

Fig. 6 highlights the performances in term of mean chip error rate (MCER), i.e. the average ratio of the number of erroneous sequences chips on the total number of sequences chips after synchronizing using either MEVB or FSNB criteria. It is shown that, after synchronizing using the MEVB-based criterion, the MCER obtained is lower than this obtained after synchronizing using the FSNB one, with a gain of 3 dB. Then, since the sequences length is 127, we can notice that in average, only one chip at most is erroneous. Those results are very important since a non cooperative context is assumed here. Hence, once a synchronization has been performed through the MEVB criterion or the FSNB one, the MCER is very low, even at very low SNRs, and it is the lowest after synchronizing using the MEVB criterion.

40th IEEE -Asilomar Conference on Signals, Systems and Computers, Oct. 29 -Nov. 1, 2006, pp. 2162-2166, Pacific Grove, California ©2006 IEEE. Personal use of this material is permitted. However, permission to reprint/republish this material for advertising or promotional purposes or for creating new collective works for resale or redistribution to servers or lists, or to reuse any copyrighted component of this work in other works must be obtained from the IEEE. Mean chip error rate after synchronizing using the FSNB criterion after synchronizing using the MEVB criterion Figure 6. Mean Chip Error Rate, N i u = 4

Fig. 7 shows the performances in term of mean bit error rate (MBER). It clearly evidences very good performances in term of MBER after synchronizing using both criteria and estimating transmitted symbols. In agreement with the results shown on Fig. 6, after synchronizing using the MEVB criteria, the MBER is the lowest, with gain of 2 dB. Typicaly, in some cooperative systems, the MBER is about 10 -8 which very close to one obtained here after synchronization using the MEVB criterion at -5 dB.

VII. CONCLUSION

In this paper, taking benefit of the blind multirate multiuser detection scheme previously proposed in [START_REF] Nsiala Nzéza | Blind multiuser detection in multirate cdma transmissions using fluctuations of correlation estimators[END_REF], led us to develop and perform a new blind synchronization method based no longer on the study of the FROBENIUS square norm behaviour (FSNB) [START_REF] Burel | Blind estimation of the pseudo-random sequence of a direct spread spectrum signal[END_REF]- [START_REF] Nsiala Nzéza | Blind synchronization and sequences identification in CDMA transmissions[END_REF], but rather on the evolution of eigenvalues maxima with analysis window shifts, denoted MEVB (Maximum EigenValue Behaviour), whithin each group of interfering users transmitting at the same data rate.

Indeed, theoretical analysis of the extended FSNB-based criterion [START_REF] Burel | Blind estimation of the pseudo-random sequence of a direct spread spectrum signal[END_REF]- [START_REF] Nsiala Nzéza | Blind synchronization and sequences identification in CDMA transmissions[END_REF] to the multiuser case, evidence the performances degradation due to an increase of synchronization peaks masking. We proved that the masking of synchronization peaks is due to the definition of the FSNB-based criterion. Then, we demonstrated that the MEVB-based citerion exhibits better performances than the FSNB-based criterion due to the lack of synchronization peaks masking. Moreover, this new method always allows the determination of the number of interfering users transmitting at the same data rate, since synchronization peaks are never masked, even at very low SNRs.

Then we showed that, synchronizing using both criteria allowed one to obtain very good peformances in terms of mean chip and bit error rates, but the lowest rates were obtained after synchronizing using the MEVB-based criterion, in agreement with theoretical results. However, computation of the FSNB criterion is faster. Consequently, the choice between those two blind synchronization criteria depends on expected -14 -12 -10 -9 -8 -5 10 Mean bit error rate after synchronizing using the FSNB criterion after synchronizing using the MEVB criterion performances and available computational power.

VIII. ACKNOWLEDGEMENT

This work was supported by the Brittany Region (France).

Figure 1 .

 1 Figure 1. Relative position of two signals and an analysis window in uplink.

  Figure 2.FSNB Criterion: case of the masking of synchronization peaks, N i u = 4

Figure 3 .

 3 Figure 3. MEVB Criterion in uplink, N i u = 4

Figure 4 .

 4 Figure 4. Experimental MEVB-based criterion, SNR =-5 dB, N i u = 4

Figure 5 .

 5 Figure 5. Probability of synchronizing at least one user, N i u = 4

Figure 7 .

 7 Figure 7. Mean bit error rate for differents SNRs, N i u = 4


Student Member, IEEE C Nsiala Nzéza

R Gautier

Member, IEEE G Burel

Blind Multiuser Detection in Multirate CDMA Transmissions Using Fluctuations of Correlation Estimators

This paper deals with the problem of blind multiuser detection in multirate direct-sequence code division multiple access (DS-CDMA). Direct-Sequence Spread Spectrum (DS-SS) signals are well-known for their low probability of interception: their statistics are similar to those of a noise; furthermore, they are usually transmitted below the noise level. Here, the method based on the fluctuations of autocorrelation estimators, previously described in a single user context [1], is extended to multiusers. From the variable spreading length technique, we will evidence an increase of these periodical fluctuations and, on average, a relationship between their peak amplitudes and sequence lengths. A theoretical analysis will show that detection is possible, even with a very low signal-to-noise ratio (SNR) at the detector input. This approach will be illustrated with experimental results. The assessment of symbol periods allows one to start the synchronization process proposed recently in [2], [3], which also permits the determination of the number of interfering users in each group of users transmitting at the same symbol period, thus, conversely at the same data rate.

I. INTRODUCTION

DS-CDMA systems are nowadays of increasing importance in wireless cellular communications because of their inclusion in most of the proposals about both terrestrial and satellitebased standards for third-generation (3G) wireless networks [START_REF]WCDMA-the radio interface for future mobile multimedia communications[END_REF], [START_REF]Satellite UMTS/IMT2000 W-CDMA air interfaces[END_REF]. One of the salient features of 3G cellular systems is the capability of supporting transmission data as diverse as voice, packet data, low-resolution video, compressed audio, etc... Since these heterogeneous services produce digital information streams with different data rates, their implementation requires the use of multirate CDMA systems where each user may transmit his data at one among a set of available data rates. An easy way to view the multirate CDMA transmission is to consider the variable spreading length (VSL) technique where all users employ sequences with the same chip period; moreover, the data rate is tied to the length of the spreading code of each user.

Many blind schemes and algorithms have been devised to either improve the performance of a CDMA receiver in a multirate multiuser context or reduce its complexity. Some prior knowledge of user parameters, e.g. the signature waveform [START_REF] Roy | Subspace blind adaptive detection for multiuser CDMA[END_REF], the processing gain, the code of a group of active users [START_REF] Wang | Group-blind multiuser detection for uplink CDMA[END_REF], the chip rate [START_REF] Haghighat | A subspace scheme for blind user identification in multiuser DS-CDMA[END_REF], [START_REF] Buzzi | Iterative cyclic subspace tracking for blind adaptive multiuser detection in multirate CDMA systems[END_REF] is always assumed, but its nature depends on the technique employed. Here, we propose a blind detection method requiring no prior knowledge about the transmitter. It is typically the case in blind signals interception in the military field or in spectrum surveillance.

These considerations led us to extend the detection scheme initially proposed in [START_REF] Burel | Detection of spread spectrum transmission using fluctuations of correlation estimators[END_REF] to multirate CDMA transmissions. This method should also allow a blind multi-standard detection through a differentiation of various standards based on an analysis of their set of available data rates [START_REF] Williams | Personal area technologies for internetworked services[END_REF].

This report is organized as follows: Section II will introduce the signal model and assumptions made. Section III will describe our blind-detection approach. Section IV will deal with the theoretical analysis of the proposed method. The simulations results will be detailed in Section V, and our conclusions will be drawn in Section VI.

II. SIGNAL MODELING AND ASSUMPTIONS

Let us consider the general case of asynchronous DS-CDMA system where each user can transmit at one out of S available data rates

R 0 < R 1 < • • • < R S-1 . By denoting N i
u the number of active users transmitting at R i and N u the total number of users such that S-1 i=0 N i u = N u , the complex equivalent of the received signal can be expressed as:

y(t) = S-1 i=0 N i u -1 n=0 +∞ k=-∞ a n,i (k)h n,i (t -kT si -T dn,i ) + b(t) (1) where h n,i (t) = Li-1 k=0 c n,i (k)p i (t-kT c ).
In (1), the subscript (•) n,i refers to the n th user transmitting at R i , denoted throughout this report as the (n, i) th user. Accordingly:

• a n,i (k) are the baseband symbols of variance σ 2 an,i for the (n, i) th user, whereas p i (t) is the convolution of the transmission filter, channel filter (which takes into account channel echoes, fading, multipaths and jammers) and receiver filter for each rate.

• The term h n,i (t) is a virtual filter corresponding to the convolution of all filters of the transmission chain with the spreading sequence {c n,i (k

)} k=0•••Li-1
, where L i is the spreading factor for the (n, i) th user.

• Because of the VSL technique, the symbol period T si for the users transmitting at the rate R i is tied to the common chip period T c : T si = L i T c , and s n,i stands for the (n, i) th signal.

• The term T dn,i is the corresponding transmission delay for the (n, i) th user; it is assumed to satisfy: 0 T dn,i < T si and to remain constant during the observation.

• b(t) is a centered white Gaussian noise of variance σ 2 b . • Signals are assumed to be independent, centered, noiseuncorrelated and received with the same power: σ 2 sn,i = σ 2 s0,0 , for all (n, i). • Finally, the SN R (in dB) at the detector intput is negative (signal hidden in the noise).

By denoting s(t) the global noise-unaffected received signal, (1) can be rewritten as:

y(t) = s(t) + b(t) = S-1 i=0 N i u -1 n=0 s n,i (t -T dn,i ) + b(t) (2)

III. PROPOSED APPROACH

The case under study is a non-cooperative context with no prior information available (unknown spreading sequence, symbol period,...). The basic principle of any detection method, whatever the application, is to take profit of the difference between the statistical properties of a signal and those of noise. Indeed, in some simple applications, the great difference between signal and noise frequencies permits the detection of a signal by filters. But, here, the application is far more complex: indeed, to reduce the probability of interception, a spread spectrum signal has to be very similar to noise, e.g. its autocorrelation and that of a white noise are both close to a Dirac function because of the pseudorandom sequence. The novelty of the proposed approach is to be based not on the autocorrelation, but on the fluctuations of its estimators. To compute the fluctuations, let us divide the received signal into temporal windows and denote T F the window duration and M the number of windows. Then, an autocorrelation estimator is applied to each window prior to the computation. By using the m th window, let us compute an estimation of the autocorrelation for any signal r(t):

R m rr (τ ) = 1 T F TF 0 r m (t)r * m (t -τ )dt (3)
where r m (t) is the signal sample over the m th window. The use of M windows allows us to estimate the second-order moment of the estimated correlation R m rr (τ ) as:

Φ(τ ) = E | R rr (τ )| 2 = 1 M M -1 m=0 | R m rr (τ )| 2 (4)
where E(•) is the estimated expectation of (•). Hence, Φ(τ ) is a measure of the fluctuations of R m rr (τ ). Applying (3) to (2) with the assumptions of independent, centered and noiseuncorrelated signals leads to the following equations for both uplink and downlink transmissions:

R ss (τ ) = S-1 i=0 N i u -1 n=0 R sn,isn,i (τ ) (5)
R yy (τ ) = R ss (τ ) + R bb (τ ). (6)
where R sn,isn,i (τ ) and R bb (τ ) are the estimates of the (n, i) th noise-unaffected signal and that of the noise autocorrelation fluctuations, respectively. Indeed, since the fluctuations are computed from many randomly-selected windows, they do not depend on the signals relative delays T dn,i .

IV. THEORETICAL ANALYSIS

To evidence the efficiency of the proposed method, even at very low SN R, let us successively investigate the contributions of noise and signal through the analysis of the secondorder moment of the autocorrelation estimator.

A. Contribution of the noise

At first, let us consider the additive noise alone, i.e. with no hidden spread spectrum signal. Indeed, since the analysis is focused on the fluctuations of autocorrelation estimators, but not on the autocorrelation itself, only the fluctuations due to the additive noise are uniformly distributed over all values of τ . Depending on the spreading sequence properties, the Multiple Access Interference (MAI) noise generates rather null, or very low, incoherent fluctuations. Since the noise is random, the fluctuations of the autocorrelation estimator, denoted Φ b (τ ), are also random. Let us characterize them by their mean (m Φ b ) and standard deviation (σ Φ b ). If the frequency response of the receiver filter is flat in [-W/2, +W/2] and zero outside, like in [START_REF] Burel | Detection of spread spectrum transmission using fluctuations of correlation estimators[END_REF], we get:

m Φ b = σ 4 b W T F (7)
where σ 2 b is the additive noise variance. The standard deviation of Φ b (τ ) is given by:

σ Φ b = 2 M σ 4 b W T F (8)

B. Contribution of the signal

Let us, now, consider only the global spread spectrum signal to show that, on average, high amplitudes of the fluctuations of the autocorrelation estimator, denoted Φ s (τ ), are obtained for each τ multiple of each symbol period

T si , i = 0, • • • , S -1.
Since the symbol periods are different, let us denote Φ i (τ ) the fluctuations of the autocorrelation estimator of the N i u signals s n,i transmitting at the same data rate T si . Let us also term m Φi the mean value of the fluctuations Φ i (τ ) for each value τ multiple of T si ; it leads to:

Conference IEEE-GLOBECOM'2006, San Francisco, California, Nov. 27 -Dec. [START_REF] Burel | Detection of spread spectrum transmission using fluctuations of correlation estimators[END_REF]2006 ©2006 IEEE. Personal use of this material is permitted. However, permission to reprint/republish this material for advertising or promotional purposes or for creating new collective works for resale or redistribution to servers or lists, or to reuse any copyrighted component of this work in other works must be obtained from the IEEE.

Φ i (τ ) = m Φi • pgn Ts i (τ ) (9)
where pgn Ts i (τ ) = +∞ k=-∞ δ(τ -kT si ), and the function δ(τ ) is defined as:

δ(τ ) = 1, if τ = 0 0, if τ = 0 (10)
Consequently, the signals being independent and centered by assumption, the fluctuations Φ s (τ ) of the global noiseunaffected signal, for each value τ multiple of T si , can be expressed as:

Φ s (τ ) = S-1 i=0 Φ i (τ ) = S-1 i=0 m Φi • pgn Ts i (τ ) (11)
Let us, now, extend the results about the single user context [START_REF] Burel | Detection of spread spectrum transmission using fluctuations of correlation estimators[END_REF] to the multirate CDMA one; since the signals are assumed to be received with the same power (equal to σ 2 s0,0 , for example) and T si = L i T c , for each value τ multiple of T si , the average value amplitude m Φi of the estimator fluctuations in each group of N i u users transmitting at the same data rate can be written as:

m Φi = N i u T si T F σ 4 s0,0 = N i u L i T c T F σ 4 s0,0 (12)

C. Remarks and analysis

At this point, it is worth noting that, on average, (12) shows an increase of fluctuation amplitudes concomitant with that of the number of users. Moreover, it indicates that, on average, the amplitude of fluctuations is proportional to the length of the sequences used.

Therefore, the longer the sequence is, the higher the peaks of autocorrelation fluctuations are, and thus the biggest amplitude is usually exhibited by the fluctuations of the N 0 u users transmitting at the lowest data rate. This approach is thus, a powerful tool to estimate symbol periods. it allows one to differentiate the various transmitted data rates and, thus, to distinguish between the different standards. Indeed, (12) shows that, for each multiple of each T si , high fluctuations of the autocorrelation estimator are obtained, and their amplitudes are tied to L i and N i u . So, in presence of a spread spectrum signal, the fluctuations curve highlights equispaced peaks whose average spacing corresponds to the estimated period symbol Tsi . It ensues that, in each group of N i u users, one can start the blind synchronization process proposed in [START_REF] Nzéza | Blind synchronization and sequences identification in CDMA transmissions[END_REF], [START_REF] Gautier | Synchronisation et estimation aveugle de séquences d'étalement pour une transmission de type CDMA en liaison descendante[END_REF] and determine the number of interfering users N i u in each group, as illustrated in section V.

Let us assume that:

T s0 = ω i T si , ∀i = 1, • • • , S -1, where {ω 1 , • • • , ω S-1
} is a set of increasing integers, i.e. all data rate are multiple of the lowest one. Therefore, whenever τ = k • ω i T s0 , k ∈ N * , a superimposition of the fluctuations peaks ensues from (11), but according to (12) their amplitudes are differnt and tied to each symbol period T si . Hence, the estimation of symbol periods is unaffected.

In addition, let us denote ρ, the SN R of a given signal such that σ 2 s0,0 = ρσ 2 b , and P t = N u σ 2 s0,0 + σ 2 b , the total received power. Their use in (12) leads to the normalized mean amplitudes of fluctuations:

m Φi = N i u T si T F P t N u + ρ -1 2 (13)
Thanks to this normalization in (13), in practice the calculation capacities are not overflowed. Moreover, on average the fluctuations amplitude is decreased, for each value τ multiple of T si , for fixed N i u and N u when the SN R at the detector input is very low. Since N i u < N u , for fixed ρ, (13) highlights a decrease of the fluctuations amplitude when the total number N u of users is increasing.

V. SIMULATION RESULTS

To determine whether a spread spectrum signal is hidden in the noise, let us compare Φ(τ ) (fluctuations of the global noise-affected signal) with theoretical mean m ©2006 IEEE. Personal use of this material is permitted. However, permission to reprint/republish this material for advertising or promotional purposes or for creating new collective works for resale or redistribution to servers or lists, or to reuse any copyrighted component of this work in other works must be obtained from the IEEE. and the simulation parameters, and thus:

m (th) Φ b = m Φ b σ (th) Φ b = σ Φ b (14)
Hence, the theoretical maximum of the noise fluctuations, i.e. the detection threshold, is deduced from (14):

m (th) Φ b + 3 • σ (th) Φ b (15)
The aim of theoretical calculation is to predict the average and the maximum of the fluctuations curve when there is no signal hidden in the noise. Whenever a spread spectrum signal is hidden in the noise, the average of the curve deviates from the theoretical average, and in particular, the maximum of the curve is above the theoretical maximum of noise fluctuations. It is precisely this difference between the theoretical predictions and the values of the obtained fluctuations that permits the detection of the signal hidden in the noise. When there is only noise, this difference is very weak.

A. Symbol periods estimation and standards differenciation

Simulations were, thus, carried out on assuming a processing gain of L 0 = 31 and L 1 = 127 (Complex GOLD sequences), and on fixing the common chip period, F c = 150 M Hz, the sampling period, F e = 300 M Hz, the window duration T F = 2 µs, and number of windows, M = 300. For the sake of simplicity, we set

N 0 u = 2, N 1 u = 2, N u = 4
, and SN R ≈ -5 dB (we assume only the MAI noise). The general case of uplink is considered in this report. Fig. 1 illustrates the fluctuations of the autocorrelation estimator Φ(τ ) versus τ (in µs). The curve clearly highlights two sets of equispaced peaks of different amplitudes. It means that two sets of spread spectrum signals transmitting at two different rates are hidden in the noise; the theoretical analysis, described above, demonstrated us that these peaks resulted from multiple symbol periods, which allowed the determination of these symbol periods (in µs): T s0 = 0.2061, T s1 = 0.8467. In this configuration, simulations showed that detection was possible up to a SN R = -25 dB at the detector input. It is worth underlining that detection remains possible for far much lower SN Rs on condition to use more numerous analysis windows, but the computation time is far much longer.

To re-estimate these values, for each set of autocorrelation fluctuations, let us look for a maximum near the farthest peak on the right side. For example, let us consider the 150 th peak and check that: 150 × T s0 ≈ 31µs and 150 × T s1 ≈ 127µs. Hence, there are maxima at the 150 th peak in agreement with Figs. 2 and3.

Table I gives the detector average relative MSE (Mean Square Error of estimated periods) for SN Rs (in dB) in the range [-14 -5] for different sequences (4 users for each type of sequences) as well as the standard deviation of the MSE. The MSE is usually negligible even at very low SN Rs, and the dispersion of estimated symbol periods around the real ones is very low. It proves the efficiency and power of our method for the blind detection of multiusers. Moreover, depending on sequence properties, this error tends to zero, especially at high SN Rs, which is the case for civilian communications.

B. Determination of the number N i u of interfering users in each group

As described in [START_REF] Nzéza | Blind synchronization and sequences identification in CDMA transmissions[END_REF], [START_REF] Gautier | Synchronisation et estimation aveugle de séquences d'étalement pour une transmission de type CDMA en liaison descendante[END_REF], the knowledge of the estimated symbol period, T si , allows one to perform the synchronization process in each group of users transmitting at the same data rate; one should note that the window duration T F is reajusted over the synchronization process so that: T F = T si .

The fast and efficient blind synchronization criterion developed in [START_REF] Nzéza | Blind synchronization and sequences identification in CDMA transmissions[END_REF], [START_REF] Gautier | Synchronisation et estimation aveugle de séquences d'étalement pour une transmission de type CDMA en liaison descendante[END_REF] In downlink, since the signals are transmitted simultanously, the criterion curve presents only one synchronization peak. Hence, the number of largest eigenvalues of one of the submatrices extracted from the global signal correlation matrix is equal to the number N i u of interfering users in the corresponding group.

In uplink, the criterion curve obtained in the interval [0, T si ] presents several peaks, and their number is equal to the number N i u of interfering users transmitting the same data rate:

• In the first group of users, i.e. T s0 = T s0 = 0.2061 µs, Fig. 4 highlights two synchronization peaks. Hence, the number of users is equal to N 0 u = N 0 u = 2, as expected.

• In the second group, i.e. T s1 = T s1 = 0.8467 µs, Fig. 5 also highlights two synchronization peaks. Consequently, there are two users transmitting in this group. So, the number of users is equal to

N 1 u = N 1 u = 2.
Thus, simulation results show that our method allows one to detect the different data rates, and, especially, to determine them via the symbol periods estimation. Then, combined with the blind synchronization method that we recently proposed in [START_REF] Nzéza | Blind synchronization and sequences identification in CDMA transmissions[END_REF], [START_REF] Gautier | Synchronisation et estimation aveugle de séquences d'étalement pour une transmission de type CDMA en liaison descendante[END_REF], our approach allows one to determine the number of interfering users transmitting at the same data rate, even with a very low SNR at the detector input.

VI. CONCLUSIONS

The scheme developed, here, for blind-detection of multiusers in multirate CDMA transmissions proved its efficiency and power both in uplink and downlink. We also highlighted its ability to estimate symbol periods. Indeed, despite the high similarity between the autocorrelation of a DS-SS signal and that of a noise, we demonstrated that the fluctuations of correlation estimators are higher when spread spectrum signals are hidden in the noise.

The average spacing evidenced on the fluctuations curves between the high, equispaced peaks allowed the determination of symbol periods. We also evidenced a relationship between fluctuation amplitudes, sequence lengths and the number of users transmitting at the same data rate. It allowed one to determine the number of user groups transmitting at the same data rate and thus to distinguish one standard from others as previously suggested in [START_REF] Williams | Personal area technologies for internetworked services[END_REF]. As proposed in [START_REF] Nzéza | Blind synchronization and sequences identification in CDMA transmissions[END_REF], [START_REF] Gautier | Synchronisation et estimation aveugle de séquences d'étalement pour une transmission de type CDMA en liaison descendante[END_REF], the number of users within a group can be found through the synchronization process.

Moreover, this blind detection scheme is applicable to passband transmissions with unknown carrier frequencies. In the single user context, it was already implemented in an operational system by scanning all the available frequency bands and recovering the carrier frequencies as done in classical transmissions.

VII. ACKNOWLEDGEMENT

This work was supported by the Brittany Region (France).

Figure 1 .

 1 Figure 1. Fluctuations Φ(τ ), N 0 u = 2, N 1 u = 2, SN R = -5 dB.

Figure 2 .

 2 Figure 2. Estimation of Ts 0 .

  (th) Φ b and maximum of the noise fluctuations Φ b (τ ). Together (7) and[START_REF] Haghighat | A subspace scheme for blind user identification in multiuser DS-CDMA[END_REF] show that m Φ b and σ Φ b are affected by only the noise variance Conference IEEE-GLOBECOM'2006, San Francisco, California, Nov. 27 -Dec.[START_REF] Burel | Detection of spread spectrum transmission using fluctuations of correlation estimators[END_REF] 2006

Figure 3 .

 3 Figure 3. Estimation of Ts 1 .

1 Figure 5 .

 15 Figure 5. Synchronization criterion in uplink, N 1 u = 2, SN R = -5 dB.

TABLE I DETECTOR

 I PERFORMANCES.

	Sequences (Length)	Average MSE	Standard deviation
	OVSF (32)	3.5685 10 -5	0.0059
	WALSH (16)	3.0162 10 -6	0.0017
	OVSF+WALSH (32, 16)	1.2 10 -2	0.1095
	OVSF+WALSH+GOLD (32, 16, 127)	8.1 10 -3	0.09

  is based on the behavior of the FROBENIUS square norm of submatrices that scan the intercepted signal correlation matrix.

					Conference IEEE-GLOBECOM'2006, San Francisco, California, Nov. 27 -Dec. 1, 2006
	2.35	x 10 6
	2.3
	2.25
	2.2
	2.15
	2.1
	2.05
	2	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8
						Time (in µs)
										1.85	x 10 6
										1.8
										1.75
										Criterion	T S
										1.7
										1.65
											synchronization criterion
											synchronization peaks
										0 1.6	0.03	0.05	0.08	0.11	0.15	0.16	0.19	0.22	0.24
											Time (in µs)

0 Figure 4. Synchronization criterion in uplink, N 0 u = 2, SN R = -5 dB.

©2006 IEEE. Personal use of this material is permitted. However, permission to reprint/republish this material for advertising or promotional purposes or for creating new collective works for resale or redistribution to servers or lists, or to reuse any copyrighted component of this work in other works must be obtained from the IEEE.


Adrien Chan-Hon-Tong

On the existence of two interior-point methods for linear programming following two different paths, and it implications

Recently, the central-path log-barrier methods (and more generally interiorpoint methods relying on central-path framework) have been proven not strongly polynomial on worse case linear program instances. A potential bypass could be to consider interior-point but not-central-path methods. In particular, this paper remarks the existence of two kind of interior-point based methods for linear feasibility which are equivalent from complexity point of view but with slightly different mechanisms. The main contribution of this paper is numerical experiments which proves that the paths followed by those two algorithms are not equal on all instances.

Introduction and motivation

The development of central-path log-barrier methods for linear programming dates back to 1980 with for example [START_REF] Renegar | A polynomial-time algorithm, based on newton's method, for linear programming[END_REF]. This framework has even been recently improved with an efficient data-structure [10] (which is a deterministic version of [START_REF] Michael B Cohen | Solving linear programs in the current matrix multiplication time[END_REF]). Using [10], a linear program related to a matrix A ∈ R M ×N with N variables, M = O(N ) constraints (matrix not too flat), and with total binary size L can be solved in less than O(N ω L) operations (where ω is the exponent of matrix multiplication/inversion).

Due to this low complexity, those algorithms are overwhelmingly the stateof-the-art of linear programming. However, this central-path framework could be challenged by recent negative result [START_REF] Xavier Allamigeon | Log-barrier interior point methods are not strongly polynomial[END_REF]: the classical central-path log-barrier methods are not strongly polynomial as 2 M steps could be required to deal with the curvatures of the central paths of worse case instances (while still being polynomial as those instances have L ≥ 2 M ). Even more, [START_REF] Xavier Allamigeon | No selfconcordant barrier interior point method is strongly polynomial[END_REF] proves that any algorithm whose current point makes line-move into the central-path of those instances is not strongly polynomial.

But, those negative results only holds for methods following the classical central-path. So, it could be interesting to consider as-efficient interior-point not-central-path algorithms for linear programming.

Indeed, there exists methods for linear program designed to deal with specific families of linear program instance in strongly polynomial times [9] for combinatorial matrix or Chubanov algorithm [START_REF] Chubanov | A polynomial algorithm for linear optimization which is strongly polynomial under certain conditions on optimal solutions[END_REF] or even the simplex for markov decision process [START_REF] Post | The simplex method is strongly polynomial for deterministic markov decision processes[END_REF]. But those methods are not related to interior point ones.

One potential example is [START_REF] Chan-Hon-Tong | A new algorithm for linear programming in critical systems[END_REF] which introduces the self-concordant Perceptron algorithm and makes the conjecture that it does not follow the central-path. Yet, this statement is not trivial as [START_REF] Chan-Hon-Tong | A new algorithm for linear programming in critical systems[END_REF] works with linear feasibility (rather than linear programming) without exploiting entirely the behavior of Newton descent on self-concordant barrier functions.

In this context, the contribution of this paper is to introduce and compare numerically two correctly equivalent methods for linear feasibility: one inspired by central-path log-barrier and one inspired by self-concordant Perceptron. The main result is that the paths followed by the two algorithms is not always the same in those numerical experiments.

Of course, there may exist worse case instances for both those algorithms simultaneously. But, the fact that those algorithms follows different paths (while having equivalent complexity) still opens a potential bypass to [START_REF] Xavier Allamigeon | Log-barrier interior point methods are not strongly polynomial[END_REF][START_REF] Xavier Allamigeon | No selfconcordant barrier interior point method is strongly polynomial[END_REF] negative results.

2 Two equivalent but different solvers for linear feasibility

Linear feasibility

Linear feasibility is the problem of finding x ∈ X A = {x ∈ R N , Ax > 0} given A ∈ R M ×N with the assumption that X A ̸ = ∅ (where 0 is the vector filled with 0). Let stress that if there exists x, Ax > 0, then, there exists χ, Aχ ≥ 1 (because one could form χ by dividing x by the minimal component of Ax which is strictly positive by assumption). As recalled in [START_REF] Chan-Hon-Tong | A new algorithm for linear programming in critical systems[END_REF], given a linear feasibility solver, one is able to solve general linear programming with same complexity (only strongly polynomial pre/post operations are required).

However, this requires in practice many steps of pre-processing. For example, linear feasibility solvers are allowed to have undefined behavior on input instances which does not satisfy core assumption X A ̸ = ∅. So, without assumption, the pre-processing should take the input (which can be unbounded or without admissible point) and cast it into a 16 times more constraint-variables matrix with four-times total binary size to ensure that the resulting A matrix verifies X A ̸ = ∅. Thus, despite it is theoretically true that linear feasibility is strongly polynomially equivalent to linear programming, this equivalent is not practical.

But, the only scope of this paper is to focus on [START_REF] Xavier Allamigeon | No selfconcordant barrier interior point method is strongly polynomial[END_REF] recent negative results. For this purpose, considering linear feasibility is not an issue.

Then, this paper offers two correctly-comparable algorithms based on selfconcordant theory for linear feasibility: one inspired by classical central-path log-barrier method and one inspired by self-concordant Perceptron.

Self-concordant theory

The two methods considered in this paper are based on self-concordant theory which describes a class of functions on which Newton descent performs efficiently. Due to the volume of this theory, this paper only states some key results, but a complete overview can be found in [START_REF] Nemirovski | Interior point polynomial time methods in convex programming[END_REF].

Notation: this paper uses log to represent log 2 i.e. log(2) = 1.

If Ψ(x) is a self-concordant function (mainly any sum of quadratic, linear, constant and -log term), with a minimum Ψ * , then, the Newton descent starting from x start allows to compute x ϵ such as Ψ(x ϵ ) -

Ψ * ≤ ε in O(Ψ(x start ) -Ψ * + log log( 1 ε
)) damped Newton steps. Each step of this descent is:

• λ Ψ (x) ← (∇ x Ψ) T (∇ 2 x Ψ) -1 (∇ x Ψ) • x ← x - 1 1+λΨ(x) (∇ 2 x Ψ) -1 (∇ x Ψ)
Even more, if Φ θ is a family of functions parameterized by θ which weights a linear term, with in addition of a common M self concordant function (e.g. a sum of M log terms). Then, the minimum of Φ (1± O(1) √ M )θ can be approximate from the minimum of Φ θ with O(1) Newton steps.

This way, the minimum of 2 L c T x- A m x). The point is that from the 2 L -related minimum one could form the solution of the linear program min

x,Ax≤b c T x. This is the core idea of the classical central-path log-barrier methods where the central-path is the set of minimums for the different values of the weight θ of the linear objective.

Theorem

Scope

The two algorithms considered in this paper are not original as it: there are direct application of self concordant theory [START_REF] Nemirovski | Interior point polynomial time methods in convex programming[END_REF] for linear feasibility. The contribution of this paper is to compare their paths (numerically) in next section. However as linear feasibility is not classical, this section describes precisely those algorithms. Let first recall that solving linear feasibility is the ability to produce x ∈ X A given a matrix A ∈ Z M ×N such that X A = {x ∈ R N , Ax ≥ 1} ̸ = ∅ with a number of steps depending on M and L the total binary size of A.

Claims

Let F, G be the two families of self-concordant functions:

F θ (v) = v T AA T v 2 1 T AA T 1 + θ 2 1 T v - m∈{1,...,M } log(v m ) (1)
G θ (x) = θ1 T Ax - m∈{1,...,M } log(A m x + 1) (2)
Claim 1: For all θ > 0, both F θ and G θ have a minimum which will be written ω(θ) and χ(θ)

Claim 2:

χ(1) = 0 and F 1 (1) -F 1 (ω(1)) ≤ M Claim 3: Both AA T ω(2 -L ) > 0 and Aχ(2 -L ) > 0 Remark: considering previous claims, the self concordant theory (see 2.2) allows to produce v, x such that AA T v > 0 and Ax > 0 with O( √ M L) Newton steps by decaying θ from 1 to 2 -L with factor (1-O(1) √ M ) performing O(1) Newton steps to approximate ω((1 -O(1) √ M )θ) and χ((1 -O(1) √ M )θ) from ω(θ) and χ(θ).

Lemmas required by the theorem

Lemma 1

For all θ ∈]0, 1], let f θ (t) = θt -log(t) and g θ (t) = θt -log(t + 1), then g θ (t) ≥ log(θ) and f θ (t) ≥ log(θ).
Proof f θ and g θ are trivially lower bounded as they goes to infinity on the extremity of their domains. As smooth convex functions, they have thus a minimum which is reached when θ -

1 t = f ′ θ (t) = 0 and θ -1 t+1 = g ′ θ (t) = 0 leading to f * θ = 1 + log(θ) ≥ log(θ) and g * θ = 1 -θ + log(θ) ≥ log(θ).
Lemma 2

∀A ∈ Z M ×N with total binary size L such that X A ̸ = ∅, the problems min p∈R N ,Ap≥1

1 T Ap (3) min q∈R N ,Aq≥1 q T q (4)
are well defined, and any optimal solutions p * , q * verify 1 T Ap * ≤ 2 L and (q * ) T (q * ) ≤ 2 L . Proof As X A ̸ = ∅ both problems admits feasible vectors. Then, both problem are bounded: by 4M for eq.3 and 0 for eq.4. Finally, some optimal solutions of these problems can be written as a linear systems extracted from A, thus, all the corresponding components can be written as a fraction of sub-determinant of A (due to Cramer rules) so 1 T Ap * and (q * ) T (q * ) are bounded by the maximal determinant extracted from A i.e. bounded by 2 L .

Proof of the theorem

2.5.1 Proof of claim 1 By construction, F θ (v) ≥ m∈{1,...,M } f θ (v m ) (as v T AA T v ≥ 0) and G θ (x) ≥ m∈{1,...,M } g θ (A m x). Then, seeing lemma 1, ∀θ ∈]0, 1], both F θ (v) ≥ M log(θ) and G θ (x) ≥ M log(θ).
Let just remark that A T 1 = 0 is impossible: let consider q * the solution of eq.4, then, (q 1)) ≤ M implying that only O(M ) Newton steps are required to initialize the algorithm (consisting of alternating Newton step and θ = 1 decaying) with a sufficient approximation of ω(1).

* ) T A T 1 = (Aq * ) T 1 ≥ 1 T 1 = M so F θ is well defined as (q * ) T A T 1 = (q * ) T 0 = 0 is impossible. 2.5.2 Proof of claim 2 Trivially, (∇G 1 ) 0 = 1 T A -1 T A = 0 so 0 = min x G 1 (x) = χ(1). Similarly, F 1 (v) ≥ m∈{1,...,M } v m -log(v m ) = m∈{1,...,M } f 1 (v m ) >= 0. So F 1 (ω(1)) ≥ 0, but, F (1) = 1 2 + M 2 . So F (1) -F (ω(

Proof of claim 3

For G θ :

Let consider a proof by contradiction: if there exists k,

A k χ(2 -L ) ≤ 0, then, let consider G 2 -L (χ(2 -L )+p * )-G 2 -L (χ(2 -L ))
where p * is any optimal solution from eq.3 problem.

Then,

G 2 -L (χ(2 -L ) + p * ) -G 2 -L (χ(2 -L )) = 1 T Ap * 2 L - m∈{1,...,M } [log(A m χ(2 -L ) + A m p * + 1) -log(A m χ(2 -L ) + 1)] < 1 T Ap * 2 L -log(A k χ(2 -L ) + A m p * + 1) -log(A k χ(2 -L ) + 1) (because Ap * ≥ 1,
so it strictly increase all log terms -but the increase could be small for m

̸ = k) < 1 -log( A k χ(2 -L )+1+1 A k χ(2 -L )+1 ) (because 1A T p ≤ 2 L seeing lemma 2) < 1 -1 = 0 (because 1+δ δ ≥ 2, ∀δ ∈]0, 1]) So G 2 -L (χ(2 -L )+p * )-G 2 -L (χ(2 -L )) < 0. This is a contradiction as χ(2 -L ) is expected to be the minimum of G 2 -L . So Aχ(2 -L ) > 0 For F θ : By definition of ω, (∇F θ ) ω(θ) = 0. Yet, ∀m ∈ {1, ..., M }, ((∇F θ ) ω(θ) ) m = A m A T ω(θ) + θ 2 -1 ω(θ)m
. So, if ω(θ) is higher bounded independently from θ, then, as soon as θ will be lower than this bound, one would have

A m A T ω(θ) = 1 ω(θ)m -θ 2 > 0. For bounding ω(θ), let first remark that ∀θ ∈]0, 1], ω(θ) T AA T ω(θ) 1 T AA T 1 ≤ 4M because otherwise F θ ( ω(θ)
2 ) ≤ F θ (ω(θ)) -3M + M : -3M because the quadratic term is higher than 4M and divided by 4 and +M because -

m log( ω(θ) 2 ) = M - m log(ω(θ)).
Independently ∀v ≥ 0, one could consider q * solution of eq.4 and the quantity (q * ) T A T v. Then:

• (q * ) T A T v ≤ ||q * || 2 v T AA T v (Cauchy inequality) • (q * ) T A T v = (Aq * ) T v ≥ 1 T v So, ∀v ≥ 0, (1 T v) 2 ≤ v T AA T v (q * ) T q * .
Combining, the bound on v and the fact that ω(θ) T AA T ω(θ)

1 T AA T 1
≤ 4M and using the lemma 2 ensuring (q * ) T q * ≤ 2 L , this leads to ∀θ ∈]0, 1], ω(θ) ≤ 2 -L 1.

In particular, for θ = 2 -L , (∇F

2 -L ) ω(2 -L ) = 0 leads to ∀m ∈ {1, ..., M }, A m A T ω(2 -L ) = 1 ω(2 -L )m -2 -L-1 > 1 2 -L-1 > 0.
3 Main results: numerical experiments on paths related to those two families Importantly, the previous section is only relevant for completeness and understanding: those F and G algorithms are not really new, and thus, not interesting by itself. But, the goal of this paper is that the paths induced by F and G have been found numerically different on some instances. This point is interesting as it is the beginning to bypass negative results from [START_REF] Xavier Allamigeon | No selfconcordant barrier interior point method is strongly polynomial[END_REF].

Comparing paths

The idea is to compare the set of A T ω(θ) and the set of χ(ν) to check if they are equivalent. If not, then the fact that there exists instance with high curvature of {χ(ν)} ν∈]0,1] has no impact on the complexity of an algorithm following

{A T ω(θ)} θ∈]0,1] .
Precisely, as it is hard to compare those two sets from geometrical point of view, an idea is to check if those two sets are close at some point. Indeed, a point z is in the area of quadratic convergence of a function

G ν iff ∃ν ∈]0, 1], (∇ z G ν ) T (∇ 2 z G ν )(∇ z G ν ) ≤ 1 4 (5)
So, this paper offers to check overlap between the F path and the G path by considering for all θ:

min ν∈]0,1] ∇ ||χ(ν)|| ||A T ω(θ)|| A T ω(θ) G ν ) T (∇ 2 ||χ(ν)|| ||A T ω(θ)|| A T ω(θ) G ν )(∇ ||χ(ν)|| ||A T ω(θ)|| A T ω(θ) G ν (6)
The idea of considering ||χ(ν)|| ||A T ω(θ)|| A T ω(θ) is to check for colinearity between χ(ν) and A T ω(θ) rather than equality. Indeed, A T ω(θ) tends to have O(M ) norm while χ(ν) tends to have a O( 1 θ ) norm. Yet, if A T ω(θ) is exactly χ(ν) with just a different scale, then, ||χ(ν)|| ||A T ω(θ)|| A T ω(θ)F= χ(ν). Code used for those experiments is in appendix. Importantly, this code is not designed for solving efficiently linear program or linear feasibility: the purpose is only to consider eq.6.

Results on random instances

This paper offers first to consider random instances. Those instance are built around two vector α, β such that 1 T A ||1 T A|| ≈ α, but, with the constraint that ∃k, A k ||A k || α ≈ -1, and, with the constraint that Aβ > 0.

First, there exists a small proportion (around 5%) of random instances for which the two paths are completely separated i.e. ∀θ,

min ν∈]0,1] (∇ ||χ(ν)|| ||A T ω(θ)|| A T ω(θ) G ν ) T (∇ 2 ||χ(ν)|| ||A T ω(θ)|| A T ω(θ) G ν )(∇ ||χ(ν)|| ||A T ω(θ)|| A T ω(θ) G ν ) > 1 4 .
For most instances, there exists θ values for which ||χ(ν)|| ||A T ω(θ)|| A T ω(θ) enters the G path. Yet, this is even more informative because, contrary to what could be though the fact that the F path overlap the G one at some point does not force the F path to continue to follow the G one.

This proves that F and G paths are both different and that the mechanism of the optimization of each family is different from the other. Thus, the curvature of the G path has no influence of the behavior of F minimization.

Results on hard instances

The code also implements the hard instance found in [START_REF] Xavier Allamigeon | Log-barrier interior point methods are not strongly polynomial[END_REF] i.e. b > 0, x 1 -t 2 b > 0, x 2 -tb > 0 x 2r > 0, x 2r-1 > 0 ∀1 ≤ j < r: t 1-1/2 j (x 2j-1 + x 2j -tx 2j+2 ) > 0, max(tx 2j-1 , tx 2j ) -x 2j+1 > 0 and the additional constraint ϵb -x 1 > 0 (to force x 1 to be close to 0 which is the optimal value).

However, even for small value of r (and moderate value of t), the difference between χ and ω is important leading to numerical issue when considering

(∇ ||χ(ν)|| ||A T ω(θ)|| A T ω(θ) G ν ) T (∇ 2 ||χ(ν)|| ||A T ω(θ)|| A T ω(θ) G ν )(∇ ||χ(ν)|| ||A T ω(θ)|| A T ω(θ) G ν ).
Yet, this shows again that χ and ω are significantly different.

Other feature highlighting potential difference

In addition to numerical results from previous section, it is also interesting to see that F and G behave differently regarding the number of small/large components of Aχ or AA T ω.

Indeed, the minimum of G θ is lower when a larger numbers of A k χ(θ) are large. Thus, considering uniquely the values of G, the convergence could be faster when the problem is highly constrained as most A k χ tend to be small during all the optimization. Inversely, the minimum of F θ is lower when a small number of A k A T ω(θ) are large. Indeed, ω(θ) k large implies

A k A T ω(θ) high as A k A T ω(θ) = 1 ω(θ) k -1 2θ
. Thus, considering uniquely the values of F , the convergence could be faster when the problem is low constrained as most A k χ tend to be high during all the optimization.

Off course, the optimization does not depend of the value of F, G (this will be the case only if each Newton step was providing the minimal improvement). Yet, this highlights that different mechanisms take place if F and G.

Perspectives

First, this paper does not claim that hard instances could not make both F and G minimization not strongly polynomial. Yet, it just remark that the paths induced by F and G are not always the same (numerically). Despite this claim could be strengthen by formal argument about the difference between F and G, this still invites to look for interior-point not-central-path algorithm.

This could be the beginning of a bypass of [START_REF] Xavier Allamigeon | No selfconcordant barrier interior point method is strongly polynomial[END_REF] negative result.

[9] Eva Tardos. A strongly polynomial algorithm to solve combinatorial linear programs. Operations Research, 1986.

[10] Jan van den Brand. A deterministic linear program solver in current matrix multiplication time. In Proceedings of the Fourteenth Annual ACM-SIAM Symposium on Discrete Algorithms, pages 259-278. SIAM, 2020. ) . sum ( ) l = t o r c h . s q r t ( l 2 ) x = x -1 / ( 3 + 3 * l ) * t o r c h . mv( d 2 g i n v , dg ) i f l < 0 . 0 5 : i f t h e t a < 0 . 9 : a l l t h e t a X . append ( ( t h e t a , x . c l o n e ( ) ) ) t h e t a = t h e t a * 0 . 9 5 i f k % 10 == 0 : p r i n t ( t h e t a ) p r i n t ( " s o l v i n g w i t h F " ) a l l t h e t a V i f l < 0 . 0 5 : i f t h e t a < 0 . 9 : a l l t h e t a V . append ( ( t h e t a , v . c l o n e ( ) ) ) t h e t a = t h e t a * 0 . 9 5 i f k % 10 == 0 : p r i n t ( t h e t a ) p r i n t ( " c h e c k d i s t a n c e b e t w e e n v and x " ) f o r t h e t a , v i n a l l t h e t a V : ATv = t o r c h . mv(AT, v ) p r o x i m i t y = [ ] f o r nu , x i n a l l t h e t a X : A T v s c a l e d = ATv * t o r c h . s q r t ( ( x * x ) . sum ( ) ) / t o r c h . s q r t ( ( ATv * ATv ) . sum ) . sum ( ) l = t o r c h . s q r t ( l 2 ) x = x -1 / ( 3 + 3 * l ) * t o r c h . mv( d 2 g i n v , dg ) i f l < 0 . 0 5 : i f t h e t a < 0 . 9 : a l l t h e t a X . append ( ( t h e t a , x . c l o n e ( ) ) ) t h e t a = t h e t a * 0 . 9 5 i f k % 10 == 0 : p r i n t ( t h e t a ) p r i n t ( " s o l v i n g w i t h F " ) a l l t h e t a V i f l < 0 . 0 5 : i f t h e t a < 0 . 9 : a l l t h e t a V . append ( ( t h e t a , v . c l o n e ( ) ) ) t h e t a = t h e t a * 0 . 9 5 i f k % 10 == 0 : p r i n t ( t h e t a ) p r i n t ( " c h e c k d i s t a n c e b e t w e e n v and x " ) f o r t h e t a , v i n a l l t h e t a V : ATv = t o r c h . mv(AT, v ) p r o x i m i t y = [ ] f o r nu , x i n a l l t h e t a X : A T v s c a l e d = ATv * t o r c h . s q r t ( ( x * x ) . sum ( ) ) / t o r c h . s q r t ( ( ATv * ATv ) . sum

Source code of numerical experiments

# x {2 j +1} < t * x {2 j -1} b e c o m e s t * x {2 j -1}-x {2 j +1}>0 A [ c o u n t e r , 2 * j -1 ] = t A [ c o u n t e r , 2 * j + 1 ] = -1 c o u n t e r += 1 # x {2 j +1} < t * x {2 j } b e c o m e s t * x {2 j -1}-x {2 j }>0 A [ c o u n t e r , 2 * j ] = t A [ c o u n t e r , 2 * j + 1 ] = -1 c o u n t e r += 1 # t ˆ( 1 -1 / 2 ˆj ) * ( x {2 j -1} + x {2 j } ) -x {2 j +2} >0 A [ c o u n t e r , 2 * j + 2 ] = 1 A [ c o u n t e r , 2 * j -1 ] = t * * ( 1 -1 / ( 2 * * j ) ) A [ c o u n t e

  m∈{1,...,M } log(b m -A m x) can be approximated in O( √ M L) Newton steps from the minimum of c T x-m∈{1,...,M } log(b m -

i m p o r t t o r c h i m p o r t numpy d e f

  s o l u t i o n (A, x ) : Ax = t o r c h . mv(A, x ) r e t u r n ( Ax > 0 ) . f l o a t ( ) . sum ( ) == A . s h a p e [ 0 ] d e f g e n e r a t e r a n d o m i n s t a n c e (M=100 , N=10 , m a r g i n = 0 . 0 0 5 ) : median = ( t o r c h . r a n d (N) * 2 -1 ) . d o u b l e ( ) median = median / t o r c h . s q r t ( ( median * median ) . sum ( ) ) mean = ( t o r c h . r a n d (N) * 2 -1 ) . d o u b l e ( ) mean -= ( mean * median ) . sum ( ) * median A = ( t o r c h . r a n d (M, N) * 2 -1 ) . d o u b l e ( ) f o r m i n r a n g e (M) : i f m < M * 4 // 5 : A [m] += mean e l s e : A [m] -= mean i f (A [m] * median ) . sum ( ) <= m a r g i n * m a r g i n : A [m] += m a r g i n * median -(A [m] * median ) . sum ( ) * median f o r m i n r a n g e (m) : a s s e r t (A [m] * median ) . sum ( ) > m a r g i n * m a r g i n r e t u r n A, median d e f F (A, t h e t a , v , AAT=None ) : i f AAT i s None : AT = t o r c h . t r a n s p o s e (A, 1 , 0 ) AAT = t o r c h . matmul (A, AT) q = ( v * t o r c h . mv(AAT, v ) ) . sum ( ) r e t u r n q / AAT. sum ( ) / 2 + t h e t a * v . sum ( )t o r c h . l o g ( v ) . sum ( ) d e f DF(A, t h e t a , v , AAT=None ) : i f AAT i s None : AT = t o r c h . t r a n s p o s e (A, 1 , 0 ) AAT = t o r c h . matmul (A, AT) AATv = t o r c h . mv(AAT, v ) f o r m i n r a n g e ( v . s h a p e [ 0 ] ) : AATv [m] = AATv [m] / AAT. sum ( ) + t h e t a -1 / v [m] r e t u r n AATv d e f D2F (A, v , AAT=None ) : i f AAT i s None : AT = t o r c h . t r a n s p o s e (A, 1 , 0 ) AAT = t o r c h . matmul (A, AT) e l s e : AAT = AAT. c l o n e ( ) AAT = AAT / AAT. sum ( ) f o r m i n r a n g e ( v . s h a p e [ 0 ] ) : AAT[m ] [ m] = AAT[m ] [ m] + 1 / v [m] * * 2 r e t u r n AAT d e f lambdaF (A, t h e t a , v , AAT=None ) : i f AAT i s None : AT = t o r c h . t r a n s p o s e (A, 1 , 0 ) AAT = t o r c h . matmul (A, AT) d f = DF(A, t h e t a , v , AAT) l 2 = ( d f * t o r c h . mv( t o r c h . i n v e r s e ( D2F (A, v , AAT) ) , d f ) ) . sum ( ) r e t u r n t o r c h . s q r t ( l 2 ) d e f b a r r i e r (A, x ) : r e t u r n t o r c h . l o g ( t o r c h . mv(A, x ) + 1 ) . sum ( ) d e f G(A, t h e t a , x ) : r e t u r n t h e t a * t o r c h . mv(A, x ) . sum ( )b a r r i e r (A, x ) d e f DG(A, t h e t a , x ) : tmp = t o r c h . z e r o s ( x . s h a p e [ 0 ] ) . d o u b l e ( ) f o r n i n r a n g e ( x . s h a p e [ 0 ] ) : tmp [ n ] = (A [ : , n ] / ( t o r c h . mv(A, x ) + 1 ) ) . sum ( ) r e t u r n t h e t a * A . sum ( 0 )tmp d e f D2G(A, x ) : N = x . s h a p e [ 0 ] tmp = t o r c h . z e r o s (N, N ) . d o u b l e ( ) f o r i i n r a n g e (N ) : f o r j i n r a n g e (N ) : tmp [ i ] [ j ] = (A [ : , i ] * A [ : , j ] / ( t o r c h . mv(A, x ) + 1 ) *

  d o u b l e ( ) f o r k i n r a n g e ( 1 0 0 0 ) : i f s o l u t i o n (A, x ) : p r i n t ( " done " ) b r e a k dg = DG(A, t h e t a , x ) d2g = D2G(A, x ) d 2 g i n v = t o r c h . i n v e r s e ( d2g ) l 2 = ( dg * t o r c h . mv( d 2 g i n v , dg )

  = [ ] t h e t a = 1 AT = t o r c h . t r a n s p o s e (A, 1 , 0 ) AAT = t o r c h . matmul (A, AT) v = t o r c h . o n e s (A . s h a p e [ 0 ] ) . d o u b l e ( ) * 0 . 1 f o r k i n r a n g e ( 1 0 0 0 ) : i f s o l u t i o n (A, t o r c h . mv(AT, v ) ) : p r i n t ( " done " ) b r e a k d f = DF(A, t h e t a , v , AAT) d 2 f = D2F (A, v , AAT) d 2 f i n v = t o r c h . i n v e r s e ( d 2 f ) l 2 = ( d f * t o r c h . mv( d 2 f i n v , d f ) ) . sum ( ) l = t o r c h . s q r t ( l 2 ) v = v -1 / ( 3 + 3 * l ) * t o r c h . mv( d 2 f i n v , d f )

1 #

 1 ( ) ) i f t o r c h . mv(A, A T v s c a l e d ) . min ( ) > -1: p r o x i m i t y . append ( ( lambdaG (A, nu , A T v s c a l e d ) , nu ) ) p r o x i m i t y = s o r t e d ( p r o x i m i t y ) p r i n t ( t h e t a , p r o x i m i t y [ 0 : 2 ] ) p r i n t ( " now w i t h h a r d i n s t a n c e s " ) r = 3 t = 1 0 0 0 n u m v a r s = 2 * r + 1 num eqns = 3 * r + 3 A = t o r c h . z e r o s ( ( num eqns , n u m v a r s ) ) . d o u b l e ( ) c o u n t e r = 0 # hack >0 A [ c o u n t e r , 0 ] = 1 c o u n t e r += -x1 > -t ˆ2 b e c o m e s -x1 + h a c k * t ˆ2 > 0 A [ c o u n t e r , 1 ] = -1 A [ c o u n t e r , 0 ] = t * * 2 c o u n t e r += 1 # x2 < t b e c o m e s -x1 + h a c k * t ˆ2 > 0 A [ c o u n t e r , 2 ] = -1 c o u n t e r += 1 # x {2 r -1} > 0 A [ c o u n t e r , 2 * r -1 ] = 1 c o u n t e r += 1 # x {2 r } > 0 A [ c o u n t e r , 2 * r ] = 1 c o u n t e r += 1 f o r j i n r a n g e ( 1 , r ) :

  r , 2 * j ] = t * * ( 1 -1 / ( 2 * * j ) ) c o u n t e r += 1 # x1 < 0 . 0 0 0 0 0 0 1 h a c k ( b e c a u s e 0 i s o b v i o u s l y t h e o p t i m a l v a l u e ) A [ c o u n t e r , 0 ] = 0 . 0 0 0 0 0 0 1 A [ c o u n t e r , 1 ] = -1 c o u n t e r += 1 p r i n t ( " s o l v i n g w i t h G" ) a l l t h e t a X = [ ] t h e t a = 1 x = t o r c h . z e r o s (A . s h a p e [ 1 ] ) . d o u b l e ( ) f o r k i n r a n g e ( 1 0 0 0 ) : i f s o l u t i o n (A, x ) : p r i n t ( " done " ) b r e a k dg = DG(A, t h e t a , x ) d2g = D2G(A, x ) d 2 g i n v = t o r c h . i n v e r s e ( d2g ) l 2 = ( dg * t o r c h . mv( d 2 g i n v , dg )

  = [ ] t h e t a = 1 AT = t o r c h . t r a n s p o s e (A, 1 , 0 ) AAT = t o r c h . matmul (A, AT) v = t o r c h . o n e s (A . s h a p e [ 0 ] ) . d o u b l e ( ) * 0 . 1 f o r k i n r a n g e ( 1 0 0 0 ) : i f s o l u t i o n (A, t o r c h . mv(AT, v ) ) : p r i n t ( " done " ) b r e a k d f = DF(A, t h e t a , v , AAT) d 2 f = D2F (A, v , AAT) d 2 f i n v = t o r c h . i n v e r s e ( d 2 f ) l 2 = ( d f * t o r c h . mv( d 2 f i n v , d f ) ) . sum ( ) l = t o r c h . s q r t ( l 2 ) v = v -1 / ( 3 + 3 * l ) * t o r c h . mv( d 2 f i n v , d f )

  ( ) ) i f t o r c h . mv(A, A T v s c a l e d ) . min ( ) > -1: p r o x i m i t y . append ( ( lambdaG (A, nu , A T v s c a l e d ) , nu ) ) p r o x i m i t y = s o r t e d ( p r o x i m i t y ) p r i n t ( t h e t a , p r o x i m i t y [ 0 : 2 ] )