Vicente Zarzoso

Parameter Estimation in Block Term Decomposition for Noninvasive Atrial Fibrillation Analysis

Keywords: Curaçao, Dutch Antilles, Dec. 10-13, 2017

Atrial fibrillation (AF) is the most common sustained cardiac arrhythmia encountered in clinical practice. Recently, a tensor decomposition approach has been put forward for noninvasive analysis of AF from surface electrocardiogram (ECG) records. Multilead ECG data are stored in tensor form and factorized via the block term decomposition (BTD). An accurate selection of parameters, including the number of block terms and the rank of the Hankel matrix factors, is necessary to guarantee physiologically significant results by this approach. The present work proposes to estimate the matrix rank by exploiting the characteristics of atrial activity during AF, which can be approximated by an autoregressive (AR) model in short records. To test this idea, three AR model order estimates are considered: Akaike's information criterion, minimum description length and partial autocorrelation function. The quality of the resulting tensor decompositions is evaluated in terms of both computational and physiologically related indices. Numerical experiments demonstrate that these model order estimation methods can find matrix rank values leading to accurate BTD approximations of the AF ECG tensor and physiologically plausible results.

I. INTRODUCTION

Held responsible of up to 25% of strokes, atrial fibrillation (AF) is the most prevalent cardiac arrhythmia, especially affecting the elderly [START_REF] January | AHA/ACC/HRS guideline for the management of patients with atrial fibrillation: A report of the American College of Cardiology/American Heart Association Task Force on practice guidelines and the Heart Rhythm Society[END_REF]. The mechanisms underlying its genesis and self-perpetuation are not yet fully understood, making the management of this challenging condition unsatisfactory. For these reasons AF is considered as the last great frontier of cardiac electrophysiology. The noninvasive analysis of AF can be performed with the aid of the surface electrocardiogram (ECG), a widespread clinical tool due to its safety and costeffectiveness. During AF the electrical wavefront propagation across the atria becomes disorganized, so that the P wave of normal atrial activation is replaced by rapidly oscillating f waves occurring all throughout the recording but masked by the QRST complex of ventricular activity (VA) at each heartbeat. Features such as amplitude and frequency of the atrial activity (AA) signal in the ECG reflect AF complexity and are linked to the probability of spontaneous cardioversion (return to normal sinus rhythm) and therapeutic success [START_REF] Bollmann | Frequency analysis of human atrial fibrillation using the surface electrocardiogram and its response to Ibutilide[END_REF]- [START_REF] Zarzoso | Non-invasive prediction of catheter ablation outcome in persistent atrial fibrillation by fibrillatory wave amplitude computation in multiple electrocardiogram leads[END_REF]. An accurate computation of such features requires the estimation of AA and the suppression of VA in surface ECG records.

The multi-sensor nature of the standard ECG, composed of 12 leads, can be exploited to perform AA extraction. The ECG data can be arranged in a space-time matrix and then decomposed using suitable matrix methods for blind source separation, e.g., principal and independent component analysis [START_REF] Rieta | Atrial activity extraction for atrial fibrillation analysis using blind source separation[END_REF]- [START_REF] Zarzoso | Extraction of ECG characteristics using source separation techniques: exploiting statistical independence and beyond[END_REF]. For matrix decompositions to be unique, constraints are necessary such as orthogonality of spatial and/or temporal factors (sources) or statistical independence between the sources. These constraints, however, may not always be coherent with physiological considerations.

To surmount this limitation, a tensor approach has recently been put forward for AF analysis [START_REF] Ribeiro | Atrial signal extraction in atrial fibrillation electrocardiograms using a tensor decomposition approach[END_REF]- [START_REF] Ribeiro | Enhanced block term decomposition for atrial activity extraction in atrial fibrillation ECG[END_REF]. The main advantage of tensor decompositions is their uniqueness under quite mild structural conditions, thus avoiding additional constraints that may lack physiological grounds. The block term decomposition (BTD) presented in [START_REF] Lathauwer | Blind separation of exponential polynomials and the decomposition of a tensor in rank-(lr, lr, 1) terms[END_REF] is chosen for its adequacy with the characteristics of AA during AF. Indeed, atrial signals may be approximated by all-pole models and can be mapped onto Hankel matrices with ranks matching the number of poles. Accordingly, the Hankel matrices associated with the ECG leads are stored in a third-order tensor, which is then processed by BTD. Experimental results show the pertinence of this approach and its potential superiority to matrix decompositions, especially in short data records and limited spatial diversity [START_REF] Ribeiro | Atrial signal extraction in atrial fibrillation electrocardiograms using a tensor decomposition approach[END_REF]- [START_REF] Ribeiro | Enhanced block term decomposition for atrial activity extraction in atrial fibrillation ECG[END_REF]. However, an accurate selection of parameters including the number of block terms and the rank of the Hankel matrix factors is crucial to BTD's satisfactory performance. The optimal choice of model parameters is a challenging open issue not only in BTD [START_REF] Lathauwer | Blind separation of exponential polynomials and the decomposition of a tensor in rank-(lr, lr, 1) terms[END_REF], [START_REF] Hunyadi | Block term decomposition for modelling epileptic seizures[END_REF] but also in other tensor decompositions.

The present contribution proposes to take advantage of the atrial signal characteristics to estimate the rank of the matrix factors in the BTD of AF ECG tensors. In short time intervals, the atrial signal admits an autoregressive (AR) model, whose order is linked to the number of poles and the rank of the associated Hankel matrix. Three popular AR model order estimation methods, namely, Akaike's information criterion, minimum description length and partial autocorrelation function, are applied to real AF ECG data in observation segments free from VA. Their model order estimates are used as matrix rank parameter in BTD. Parameter selection quality is assessed by computational as well as physiological indices measuring tensor approximation and atrial signal estimation performance, showing the appropriateness of the proposed approach.

II. BLOCK TERM DECOMPOSITION

AA extraction in AF ECG records can be formulated as the blind source separation problem:

X = MS (1)
where each row of X 2 IR K⇥N contains the N time samples of one ECG lead, S 2 IR R⇥N contains the source signals and M 2 IR K⇥R denotes the mixing matrix. Matrix S contains the atrial and ventricular sources as well as other artifacts and interference contributing to the ECG. Only matrix X is observed in the above model, and the goal is to estimate the two matrix factors in eqn. (1) from the observation of X. The value of K depends on the number of leads considered. When processing a standard ECG record, we have K = 12 leads.

In [START_REF] Lathauwer | Blind separation of exponential polynomials and the decomposition of a tensor in rank-(lr, lr, 1) terms[END_REF], BTD is put forward to solve the BSS problem using a tensor decomposition approach. The idea behind BTD is to map the kth row of X, denoted x k , onto a Hankel matrix, denoted H

x k , with entries:

[H x k ] i,j = [x k ] i+j 1 (2)
which is then stacked along the mode-3 of a third-order tensor X 2 IR I⇥J⇥K :

X :,:,k = H x k k = 1, 2, . . . , K. (3)
Hankel matrix H

x k is of dimensions (I ⇥ J), where I = N/2, J = N/2 + 1 if N is even, and I = J = (N + 1)/2 if N is odd. According to model [START_REF] January | AHA/ACC/HRS guideline for the management of patients with atrial fibrillation: A report of the American College of Cardiology/American Heart Association Task Force on practice guidelines and the Heart Rhythm Society[END_REF], tensor X admits the following decomposition in R block terms:

X = R X r=1 H sr m r (4)
where m r denotes the rth column of M and the source Hankel matrix H sr is defined as in eqn. [START_REF] Bollmann | Frequency analysis of human atrial fibrillation using the surface electrocardiogram and its response to Ibutilide[END_REF]. If the matrix factor H sr has rank L r , the rth block term in model ( 4) has multilinear rank-(L r , L r , 1). Sufficient conditions for the uniqueness of BTD model ( 4) are given in [START_REF] Lathauwer | Blind separation of exponential polynomials and the decomposition of a tensor in rank-(lr, lr, 1) terms[END_REF]:

C1) Matrices A = [A 1 , A 2 , . . . , A R ] and B = [B 1 , B 2 , . . . , B R ] with dimensions I ⇥ P R r=1 L r and J ⇥ P R r=1 L r , respectively, are full column rank, where H sr = A r B T r , r = 1, 2, . . . , R. C2) Mixing matrix M does not contain collinear columns.
Condition C1 requires in particular that I, J P R r=1 L r , which would impose a theoretical lower bound on the sample size N . This condition, however, is sufficient but not necessary.

AA in AF presents a narrowband power spectrum, typically concentrated around a peak frequency in the 3 to 9 Hz band. As a result, the atrial sources accept the all-pole or exponential model:

[S] r,n def = s r (n) = Lr X `=1 c r,`z (n 1) r,`n = 1, 2, . . . , N (5)
where z r,`i s the `th pole of the rth source and c r,`s cales the exponential term. A signal with L r distinct poles is characterized by a Hankel matrix with rank L r . If the poles are different, condition C1 holds for sufficient sample size and decomposition (4) is unique. Milder conditions guarantee uniqueness even in the case of common poles [START_REF] Lathauwer | Blind separation of exponential polynomials and the decomposition of a tensor in rank-(lr, lr, 1) terms[END_REF]. If these conditions are fulfilled, the BTD of tensor X defined in (3) will yield the mixing matrix columns as vector factors and the source Hankel matrices as matrix factors, like in eqn. (4). According to (2), the source signals can then be recovered by averaging along the antidiagonals of the associated Hankel matrices.

As recently shown in [START_REF] Ribeiro | Atrial signal extraction in atrial fibrillation electrocardiograms using a tensor decomposition approach[END_REF]- [START_REF] Ribeiro | Enhanced block term decomposition for atrial activity extraction in atrial fibrillation ECG[END_REF], BTD can provide adequate results and outperform matrix approaches in scenarios with limited observation length and spatial diversity. Nevertheless, its performance is highly dependent on the appropriate choice of model parameters, i.e., the number of block terms R and the rank of the matrix factors L r . In the sequel, we will assume for simplicity that all sources have the same number of poles, L r = L, r = 1, 2, . . . , R, a typical assumption in previous works. The remaining of the paper focuses on the estimation of the matrix factor rank L by means of AR analysis.

III. AR MODELING

As pointed out by all-pole equation ( 5), the sources accept the AR model:

s r (n) = L X `=1 a r,`sr (n `) + e r (n) (6)
where e r (n) is the error signal. AR coefficients a r,`d efine the Lth degree polynomial p r (z

) = z L P L `=1 a r,`z L `,
whose roots are the poles z r,`. The AR model coefficients can be estimated by solving the Yule-Walker equations via the efficient Levinson-Durbin recursion or well-known alternative techniques. As a result, parameter L may be obtained by suitable AR model order estimation methods. Among the most popular existing techniques, we consider three such estimators in this work: Akaike's information criterion, minimum description length and partial autocorrelation function, which are briefly recalled next.

Akaike's information criterion (AIC) estimates the AR model order as the value of L minimizing the cost function:

AIC(L) = N log " L + 2L
where " L represents the prediction error, i.e., the power of the error term e r (n) in eqn. [START_REF] Rieta | Atrial activity extraction for atrial fibrillation analysis using blind source separation[END_REF] when the model order is set to L. AIC tends to underestimate the model order when the signal under analysis is not an AR process, and tends to overestimate it as the sample size N increases [START_REF] Hayes | Statistical Digital Signal Processing and Modeling[END_REF]. The minimum description length (MDL) is based on the criterion:

MDL(L) = N log " L + L log N.
MDL is a consistent model order estimator, as it converges to the true order as N ! 1 when the signal follows an AR model [START_REF] Hayes | Statistical Digital Signal Processing and Modeling[END_REF]. Finally, the partial autocorrelation (PAC) function criterion is also considered. The PAC sequence can be directly computed from the reflection coefficients obtained in the solution of the Yule-Walker equations [START_REF] Hayes | Statistical Digital Signal Processing and Modeling[END_REF]. We estimate L as the maximum index of the PAC coefficients lying outside the statistically significant threshold of ±1.96/ p N . For its increased robustness and stability, we choose the modified covariance method (also known as forward-backward or least squares method) [START_REF] Hayes | Statistical Digital Signal Processing and Modeling[END_REF] to estimate the AR model coefficients and prediction error used in the AIC and MDL methods.

In practice, the source signals are unknown and their rank cannot be estimated as explained above. To circumvent this problem, we propose to compute as an initial approximation the AR model order estimates in the TQ segments of the ECG. TQ segments are the time intervals between two consecutive beats and therefore only contain atrial signal contributions, as illustrated Fig. 1.

IV. EXPERIMENTAL RESULTS

This section performs a series of numerical experiments on synthetic and real data to test the ideas presented above. In the first place, synthetic signals in controlled conditions are employed to assess the influence of the selected model parameters on BTD's computational performance in terms of tensor reconstruction quality and number of iterations for convergence. These indices are shown to improve when the selected parameters approach the ground truth. Next, we illustrate the application of the AR model based matrix rank estimators in a real AF ECG record, with the aim of optimizing the computational performance of BTD while providing adequate atrial signal estimation quality and physiologically reasonable results in this biomedical problem. In the following experiments, BTD is computed by the nonlinear least squares (NLS) implementation available in Tensorlab MATLAB toolbox [START_REF] Vervliet | [END_REF]. The maximum number of iterations is set to 1000.

A. Synthetic Data

A mixture of R = 5 sources and K = 5 observations is generated with a mixing matrix composed of normalized Gaussian random entries. Each source is a sum of 6 sinusoids with random frequencies and N = 250 samples. Since a sinusoid is defined by two complex conjugate poles, we have L = 12 as the true theoretical value of the Hankel matrix factors associated with each source signal. Figure 2 shows the normalized mean squared error (NMSE), defined as k X X k 2

Fro /kX k 2 Fro , of the tensor reconstructed for different values of the matrix factor rank L assumed by BTD. We plot the best tensor NMSE out of 10 independent random initializations of the BTD factors for each L. Also shown are the corresponding NMSE of the estimated sources and the number of iterations required by the NLS-BTD algorithm for convergence. Performance indices are computed after permutating and scaling the estimated temporal factors through a greedy algorithm based on the normalized crosscorrelation with the actual sources. One can observe that these computational performance indices drop significantly when the value of L used in the algorithm approaches the ground truth and, therefore, they can be used as a guide to estimate the BTD matrix factor rank. from the first patient of an AF signal database acquired during catheter ablation procedures at the Cardiology Department of Princess Grace, Monaco. The record, over a minute long and acquired at a 977 Hz sampling rate, is preprocessed by a zero-phase forward-backward type-II Chebyshev bandpass filter with 0.5 and 40 Hz cutoff frequencies to suppress highfrequency noise and baseline wandering. For BTD processing, the single beat-to-beat interval (QQ segment) with the longest TQ interval duration is considered, of which bipolar limb lead II is shown in Fig. 1. The full beat and the TQ segment are composed of 1468 and 1017 samples, respectively.

2) BTD factor rank estimation via AR modeling: The AR order estimation methods described in Sec. III are applied to the TQ segment in each ECG lead separately, yielding the values shown in Table I. For each method, the model orders averaged over the 12 ECG leads are considered as an estimate of the matrix factor rank of the atrial sources present in the observed AF ECG; hence: LAIC = 95, LMDL = 48 and LPAC = 17. 3) BTD computational performance: Next, the whole single-beat 12-lead ECG segment, including the QT interval containing VA, is downsampled by a factor of two and hankelized as explained in Sec. II, resulting in a third-order tensor with dimensions (367⇥368⇥12). Downsampling is necessary as the original tensor with dimensions (734 ⇥ 735 ⇥ 12) poses difficulties to Tensorlab for some block term configurations considered in these experiments. The tensor is processed by NLS-BTD assuming different values of the number of block terms R and the matrix factor rank L. The best reconstructed tensor NMSE over 10 independent random initializations and the corresponding number of iterations for convergence are shown in Fig. 3 (top). A noticeable decrease in the number of iterations can be observed from L = LMDL = 48. Tensor approximation accuracy also improves significantly after that value of L for a sufficient number of block terms.

4) Atrial signal analysis: Apart from these indices of computational performance, other physiologically motivated measures quantifying AA content are also considered such as dominant frequency (DF) and spectral concentration (SC), which are calculated as described in [START_REF] Castells | Spatiotemporal blind source separation approach to atrial activity estimation in atrial tachyarrhythmias[END_REF]. This experiment assumes, as in previous works, that all atrial activity can be represented by a single source, which is selected as the signal with the highest SC among the sources with DF in the typical [START_REF] Matsuo | Clinical predictors of termination and clinical outcome of catheter ablation for persistent atrial fibrillation[END_REF][START_REF] Ribeiro | Atrial signal extraction in atrial fibrillation electrocardiograms using a tensor decomposition approach[END_REF]-Hz frequency band. Figure 3 (bottom) plots the DF and SC indices of the atrial source in the best tensor approximation solution. The obtained DF does not show significant variations over the considered range of BTD parameters but remains physiologically plausible as it lies close to the DF of an intracardiac electrogram (EGM) acquired with the aid of a catheter located in the left atrial appendage, DF EGM = 6.8 Hz. Higher values of SC, reflecting clearer atrial signal components, are generally obtained as the number of block terms increases. Fig. 4 plots the AA signal reconstructed by BTD with R = 12 block terms and the different AR-based matrix rank estimates. Reconstructions are obtained with random initializations of the tensor factors different from those of Fig. 3. Also shown as ground truth is the result by the adaptive singular value cancellation (ASVC) method for QRST template matching and subtraction [START_REF] Alcaraz | Adaptive singular value cancelation of ventricular activity in single-lead atrial fibrillation electrocardiograms[END_REF]. Following closely the guidelines of [START_REF] Alcaraz | Adaptive singular value cancelation of ventricular activity in single-lead atrial fibrillation electrocardiograms[END_REF], ASVC is computed as the SVD-based synchronized average of the optimal number of QT intervals in the full observed recording (88-s long in this example) that are most similar in terms of correlation to the QT interval under analysis. Processing only the single-beat interval, the BTD estimates achieve a level of atrial signal quality comparable to ASVC's. BTD may even appear superior in the QT interval, where ASCV seems to underestimate the existing AA contribution as quantified by the QT-to-TQ segment power ratio P QT /P TQ . For reference, this ratio reaches 8.6 dB in lead V1.

V. CONCLUSIONS

The present contribution has proposed to estimate the rank of the matrix factors in the BTD of AF ECG records by computing the AR order of the observed TQ segments, containing atrial contributions only. In the illustrative example, AIC and MDL yield rank estimates leading to adequate performance of BTD in terms of tensor approximation error, computational complexity and atrial signal parameters, for a sufficient number of block terms. The proposed approach may arguably be considered as suboptimal since all block terms are assumed to have the same matrix rank and the rank of the observed mixture is expected to be higher than that of individual sources as soon as any two sources do not have the same poles. Despite these apparent limitations, experimental results are encouraging. Further work should aim at the estimation of the number of block terms, assess the impact of atrial signal temporal variability on the tensor estimates and experimentally evaluate their performance in a full AF ECG database.

ACKNOWLEDGMENT

The author is a member of the Institut Universitaire de France.

  Fig.1. Single-beat AF ECG segment processed by BTD in the real data experiments of Sec. IV-B. The segment is composed of an initial QT interval, containing a mixture of AA and VA, followed by a TQ segment, in which only AA is present. The rank of the BTD matrix factors is computed by several AR model order estimators from the TQ segment. Then the full interval (QT+TQ) is decomposed via BTD. All 12 leads of the standard ECG are considered in the analysis, but only bipolar limb lead II is shown here for clarity.

Fig. 2 .

 2 Fig. 2. Synthetic data experiments: influence of BTD matrix factor rank L. NMSE of reconstructed tensor for the best BTD initialization (solid line), corresponding source reconstruction NMSE (dashed line) and number of iterations for BTD convergence (dotted line, right axis). The dashed vertical line shows the actual theoretical value of the source Hankel matrix rank (L = 12).

Fig. 3 .Fig. 4 .

 34 Fig. 3. Experiments on real AF ECG data: influence of number of BTD block terms R and matrix factor rank L. (Top) NMSE of reconstructed tensor for best BTD initialization (solid lines, left axis) and corresponding number of iterations for BTD convergence (dotted lines, right axis). (Bottom) Dominant frequency (solid lines, left axis) and spectral concentration (dotted lines, right axis) of the corresponding atrial source.

  1. Single-beat AF ECG segment processed by BTD in the real data experiments of Sec. IV-B. The segment is composed of an initial QT interval, containing a mixture of AA and VA, followed by a TQ segment, in which only AA is present. The rank of the BTD matrix factors is computed by several AR model order estimators from the TQ segment. Then the full interval (QT+TQ) is decomposed via BTD. All 12 leads of the standard ECG are considered in the analysis, but only bipolar limb lead II is shown here for clarity.

TABLE I .

 I AR MODEL ORDER ESTIMATES FOR THE PROCESSED TQ SEGMENT OF THE AF ECG RECORD. VALUES REPRESENT MEAN ± STANDARD DEVIATION COMPUTED OVER THE 12 LEADS.

	AIC	MDL	PAC
	95 ± 3.5	48± 5	17± 5.2


Rubén Martín-Clemente

Vicente Zarzoso
email: zarzoso@i3s.unice.fr

L1-PCA Signal Subspace Identification for Non-sphered Data Under the ICA Model

Principal component analysis (PCA) is an ubiquitous data compression and feature extraction technique in signal processing and machine learning. As compared with the classical L2-norm PCA, its L1-norm version offers increased robustness to outliers that are usually present in faulty data. Recently, L1-PCA was shown to perform source recovery when the observed data follow an independent component analysis (ICA) model. However, proof of this result requires the data to be sphered, i.e., to be preprocessed to constrain their covariance matrix to be the identity. The present contribution extends this result by relaxing the sphering assumption and allowing the data to have arbitrary covariance matrix. We prove that L1-PCA is indeed able to identify the mixing matrix columns associated with the strongest independent sources, thus performing signal subspace identification with improved robustness to outliers. Numerical experiments illustrate and confirm the theoretical findings.

I. INTRODUCTION

Principal Component Analysis (PCA) is probably the most popular multivariate data analysis technique [START_REF] Jolliffe | Principal Component Analysis[END_REF], as it constitutes the basis for a variety of dimensionality reduction and noise filtering approaches. Multivariate data can be considered as vectors in a high-dimensional vector space, say R p . The primary goal of PCA, as defined originally by Pearson, is to find the direction of the line that best fits the data in the p-dimensional space. In the traditional approach, this goal is obtained by maximizing the sum of the squares of the data projections onto the sought direction. Actually, several best fitting directions, or principal axes, can be computed by solving repeatedly the above optimization problem under the constraint that the nth principal axis has to be orthogonal to the previous (n 1) axes. The subspace V ⇢ R p spanned by the orthogonal basis vectors pointing in the main principal directions is called the signal subspace, another fundamental concept in multivariate signal processing.

It is intuitive that the principal axes reveal much of the structure of the data. In the context of the processing of faulty data, PCA is mostly useful when the data of interest is concentrated around only a few principal axes, so that the dimension of the signal subspace V is much lower than p, whereas the noise is distributed isotropically through all R p . In this case, a significant reduction of the noise, while preserving the information content of the original data, can be obtained by projecting the faulty data onto the signal subspace.

Recently, L1-norm PCA has been proposed in an attempt to increase robustness to outliers of classical PCA. Experi-mentally, the L1-PCA approach has proven to be effective for restoration of faulty data in, e.g., the reconstruction of occluded images [START_REF] Kwak | Principal component analysis based on L1-norm maximization[END_REF], [START_REF] Markopoulos | Optimal algorithms for L1-subspace signal processing[END_REF], [START_REF] Markopoulos | Efficient L1-Norm Principal-Component Analysis via Bit Flipping[END_REF], as well as in pattern recognition and dimensionality reduction [START_REF] Brooks | A pure L 1 -norm principal component analysis[END_REF], [START_REF] Kwak | Principal component analysis based on L1-norm maximization[END_REF]- [START_REF] Markopoulos | Optimal algorithms for L1-subspace signal processing[END_REF], [START_REF] Meng | Improve robustness of sparse PCA by L1 -norm maximization[END_REF], [START_REF] Pang | Robust tensor analysis with L1-norm[END_REF], [START_REF] Yu | An efficient algorithm for L1-norm principal component analysis[END_REF]- [START_REF] Wang | Block principal component analysis with L1-norm for image analysis[END_REF]. Although fast yet suboptimal algorithms for L1-PCA exist [START_REF] Kundu | Fast computation of the L1-principal component of real-valued data[END_REF], [START_REF] Kwak | Principal component analysis based on L1-norm maximization[END_REF], one of the most attractive features of this approach is that it can be also performed by optimal algorithms with guaranteed global convergence, as recently shown in [START_REF] Markopoulos | Optimal algorithms for L1-subspace signal processing[END_REF].

The link between L1-PCA and independent component analysis (ICA), another popular multivariate data processing technique [START_REF] Comon | Handbook of Blind Source Separation: Independent component analysis and applications[END_REF], [START_REF]Independent Component Analysis: Principles and Practice[END_REF], is carried out in [START_REF] Martín-Clemente | On the link between L1-PCA and ICA[END_REF], thus opening the possibility to perform ICA using globally optimal algorithms for L1-PCA. L1-PCA is shown to extract the independent sources under the ICA model in the case where the data are sphered, i.e., they have been whitened to present an identity covariance matrix. The present contribution extends this result by relaxing the assumption that the data be sphered, which defines a more general framework. A theoretical analysis supported by experimental results shows that L1-PCA is able to identify the dominant signal subspace under the ICA model, thus enabling source recovery via beamforming in a further processing stage.

II. PROBLEM FORMULATION

Let us suppose that we observe N samples {x 1 , x 2 , . . . , x N } of a p-dimensional zero-mean random vector x, which can be arranged into the (p ⇥ N ) matrix

X = [x 1 , . . . , x N ].
Let w 2 R p , kwk 2 = 1, define a projection direction. As mentioned above, the goal of classical PCA is to find the random variable y = w T x with the largest empirical variance. The first principal component vector solves the following problem:

max kwk 2 =1 N X n=1 (w T x n ) 2
or, in matrix form:

max kwk 2 =1
w T X 2 .

(

) 1
Though mathematically appealing, the L2-norm is rather sensitive to impulsive noise or outliers since squaring the projections overemphasizes the effects of far-off data points, which is a serious drawback.

To increase the robustness of PCA when processing faulty data, we can try to eliminate the outliers from the data or, alternatively, replace the L2-norm cost function by a more robust criterion. The following L1-norm based variant of traditional PCA is proposed in [START_REF] Kwak | Principal component analysis based on L1-norm maximization[END_REF]:

max kwk 2 =1 w T X 1 (2)
where kyk 1 = P N n=1 |y i | represents the sum of the absolute entries of vector y = w T X. The above two PCA criteria are generally not equivalent, even when the data contain no outliers. A key problem is to determine under which conditions the solutions of ( 1) and ( 2) are the same, or at least similar, in the absence of outliers. When the data are sphered, i.e., the covariance matrix R

x = 1 N XX T equals the identity matrix, it follows that w T X 2 = kwk 2 and therefore no principal axis can be defined or identified according to [START_REF] Brooks | A pure L 1 -norm principal component analysis[END_REF]. This is no surprising as sphering is often performed by classical PCA, so that the information that the L2 criterion can extract from sphered data is somewhat exhausted. By contrast, L1-PCA can yield interesting directions even for sphered data. A generative model that can provide useful insights is the case where each data point x i is made up of the linear superposition

x i = q X j=1 a j s ij (3)
where s ij represents the degree to which a j 2 R p is present in observation x i and q is the actual dimension of the data.

The assumption that coefficients s ij are mutually independent for different values of j is physically plausible in many realworld problems, and gives rise to the independent component analysis (ICA) model [START_REF] Comon | Handbook of Blind Source Separation: Independent component analysis and applications[END_REF]. The solutions of L1-PCA under the ICA model when p = q and the data are sphered are analyzed in depth in [START_REF] Martín-Clemente | On the link between L1-PCA and ICA[END_REF], showing that projections associated with independent components are locally stable stationary points of the L1-PCA criterion. The remaining of this paper explores the extension of this result to the more general case of non-sphered data, where R

x 6 = I, and proves that L1-PCA yields principal directions in a similar sense to L2-PCA but with the benefit of increased robustness to outliers.

III. THEORETICAL ANALYSIS

We conduct a preliminary study for the case q = 2, implying that the data points can be described by only two independent components. This case is simple enough to allow a theoretical analysis, while still retaining the important features of the problem. Results are stated without proof for lack of space.

A. General Analysis for q = 2

We start by adopting a distributional viewpoint so that, according to model (3), the observed data vector x can be expressed as

x = a 1 s 1 + a 2 s 2 = As (4)
with A = [a 1 , a 2 ] and s = [s 1 , s 2 ] T . Components s 1 and s 2 represent zero-mean, mutually statistically independent random variables, which are also assumed to have unit variance. If this is not the case, the power of the variables can always be incorporated into the length of the vectors a i without loss of generality. Given a projection unit-norm vector w, the projected data are defined by y = w T x = g 1 s 1 + g 2 s 2 , where g i = w T a i . Denoting by a ji the jth entry of a i , we can also write

g i = p X j=1 w j a ji .
Observe that, assuming ergodicity conditions, the L1-norm criterion in [START_REF] Chong | An Introduction to Optimization[END_REF] becomes (up to an irrelevant scale factor) proportional to E{|y|} for large enough sample size:

kyk 1 ! N !+1
N E{|y|} where y = w T X. For this reason, instead of (2), let us consider here the equivalent problem

max w E{|y|} subject to kwk 2 = 1
whose Lagrangian is given by

L(w, ) = 1 2 E{|y|} + (kwk 2 2 1) (5)
where is the Lagrange multiplier. Differentiating with respect to the components of w by the chain rule, and using some formulas in [START_REF] Martín-Clemente | On the link between L1-PCA and ICA[END_REF], it can be proven that

@L @w n = 2 X i=1 a ni G i + 2 w n
where

G i = Z 1 1 sf i (s)F j ✓ g i g j s ◆ ds (6)
with i, j 2 {1, 2}, i 6 = j,. Symbols f i (•) and F j (•) denote the probability density function and the distribution function, respectively, of s j . This equation can be rewritten in matrix form as @L @w = Ab + 2 w

where b = [G 1 , G 2 ] T .
Finally, the stationary points of the criterion must verify:

@L @w = 0 ) Ab = 2 w. ( 7
)
This equation is difficult to solve because b depends nonlinearly on w through g in [START_REF] Kwak | Principal component analysis based on L1-norm maximization[END_REF]. Hence, rather than searching for a general solution, we will consider a simple yet significant particular problem for now.

B. Particular Case: Data Clustered Around 1D Subspace Now let us focus on the important situation where the data appear approximately clustered around a one-dimensional subspace, a straight line. Without loss of generality, we assume that this line is in the direction of a 1 , the case of a 2 being treated in a totally analogous manner. A necessary condition for the observed data clustering is that ka 1 k 2 ka 2 k 2 under the unit-variance source assumption. Vector a 1 defines a principal axis in the sense of the dominant signal subspace, while the term a 2 s 2 is interpreted here as a small amount of added noise defining a noise subspace. We are going to study whether L1-PCA is able to identify the signal subspace. This is tantamount to verifying whether

w ⇡ ± a 1 ka 1 k 2 (8)
is a solution to [START_REF] Kwak | Feature extraction for one-class classification problems: enhancements to biased discriminant analysis[END_REF]. Since, by definition,

g = w T A ⇡ ± h 1, a T 1 a 2 /ka 1 k 2 i
we can assume that g2 g1 ' 0. Then, the following Maclaurinbased approximation of ( 6) holds for i = 2:

G 2 ⇡ Z 1 1 sf 2 (s)  F 1 (0) s f 1 (0) g 2 g 1 ds = f 1 (0) g 2 g 1
where we have exploited the assumption that s 2 has zero mean and unit variance. Clearly, g2 g1 ' 0 implies G 2 ' 0. Substituting G 2 = 0 in (7), we readily arrive at [START_REF] Li | L1-norm based 2DPCA[END_REF]: this proves that a 1 is a stationary point of the L1-norm criterion. We know check its stability.

C. Local Stability Analysis

Let L w be the Hessian matrix associated to [START_REF] Kundu | Fast computation of the L1-principal component of real-valued data[END_REF], L w = F w + H w , where

H w = " @ 2 @w i @w j p X n=1 w 2 n # i,j i, j = 1, . . . , p
is the Lagrange multiplier and F w is the Hessian of 1 2 E{|y|},

F w def =  @ 2 @w i @w j 1 2 E{|y|} i,j i, j = 1, . . . , p.
Assuming the solution [START_REF] Li | L1-norm based 2DPCA[END_REF], it can be shown that

L w = f 1 (0) ka 1 k 2 a 2 a T 2 ka 1 k 2 ✓Z 1 0 sf 1 (s) ds ◆ I.
Let v be any vector orthogonal to w. Geometry shows that

v T a 2 = ⌥kvk 2 ka 2 k 2 sin(\a 1 a 2 )
where \a 1 a 2 is the angle between a 1 and a 2 . Then,

v T L w v kvk 2 2 = ka 2 k 2 2 ka 1 k 2 sin 2 (\a 1 a 2 )f 1 (0) ka 1 k 2 Z 1 0 sf 1 (s) ds.
When ka 1 k 2 is greater enough than ka 2 k 2 , it follows that v T L w v < 0. According to [START_REF] Chong | An Introduction to Optimization[END_REF], this ensures that (8) maximizes the criterion. Otherwise, if f 1 (0) R 1 0 sf 1 (s) ds, as in a supergaussian distribution [START_REF] Martín-Clemente | On the link between L1-PCA and ICA[END_REF], then v T L w v may be positive:

in this case, we have to minimize, not maximize, the L1-norm to recover a 1 .

According to this theoretical analysis, L1-PCA determines the signal subspace in a way similar to L2-PCA. This result, combined with the fact that L1-PCA is more robust against outliers than L2-PCA, provides additional support for the use of L1-PCA in the analysis of faulty data.

IV. NUMERICAL EXPERIMENTS

Some computer experiments are carried out to support the validity and generality of the theoretical development of the previous section. For simplicity we will assume that p = q, i.e., the observations have the same dimensionality as the underlying independent components generating the data. [START_REF] Brooks | A pure L 1 -norm principal component analysis[END_REF] T . These vectors make an angle of 45 and 30 , respectively, with the horizontal axis so that R

A. Case q = 2 Let a 1 = ↵ p 2 [1, 1] T and a 2 = 1 2 [ p 3,
x = AA T 6 = I. Assuming that w ⇡ ± a1 ka1k2 as in [START_REF] Li | L1-norm based 2DPCA[END_REF], it follows that ka 1 k 2 ka 2 k 2 , as required, when ↵ ! 1. We generate N = 200 samples of the random variable x using model (4), where s 1 and s 2 are drawn from independent zero-mean unit-variance uniform distributions. L1-PCA is then applied on the data using the iterative algorithm presented in [START_REF] Kwak | Principal component analysis based on L1-norm maximization[END_REF]. First, Figure 1 (top) shows vector a 1 for ↵ = 1.5 (solid-line arrow), the vector w L1 obtained by L1-PCA (bold solid line), the first principal component w L2 generated by traditional L2-PCA (dashed line) and the scatter plot of the observations of x. In the figure, vectors a 1 , w L1 and w L2 have been shifted and scaled differently for clarity. In this particular experiment, observe that L1-PCA estimates a 1 , the vector pointing in the main direction of the data, better than L2-PCA: the angle between w L1 and a 1 equals 3.2 , whereas the angle between w L2 and a 1 is 6.3 . Next, we repeat the experiment after adding several outliers to the same data. Results are shown in Figure 1 (bottom). As expected, the estimation is negatively affected by the presence of outliers: the angle between w L1 and a 1 now becomes 3.6 . Yet L1-PCA proves more robust than L2-PCA, since the angle between w L2 and a 1 is now equal to 12.5 .

Figure 2 (left) shows the angle of vector w L1 determined by L1-PCA with the horizontal axis as a function of ↵ after convergence of the L1-PCA algorithm. As ↵ increases, this angle tends to 45 so that w L1 is parallel to a 1 as expected. On the other hand, if ↵ tends to zero, the angle approaches 30 and w L1 coincides with a 2 , which is also consistent with our analysis as a 2 now becomes dominant. Finally, when ↵ = 1, so that a 1 and a 2 have the same length, w L1 points in the direction of the sum of both vectors. In Figure 2 (right), as in a previous experiment, data contain 5% of Gaussian outliers and therefore the figure also compares the robustness of L1-PCA and L2-PCA.

B. Case q > 2

We now generate N = 500 independent samples of x = P 7 j=1 a j s j , x 2 R 7 . The coordinates of vectors a j are drawn from a standardized Gaussian distribution, except for a 4 , where the variance is set to 10, so that a 4 defines the main axis of the data. The variables s j are zero-mean and uniformly distributed. Figure 3 shows the angle of vector w L1 recovered by L1-PCA with each of the basis vectors a i , simply computed as arccos ⇣ |w T ai| kwkkaik

⌘

, averaged over 100 independent Monte Carlo realizations. A confidence interval is also depicted: the total length of the vertical line equals twice the standard deviation of the angle. As expected, w L1 is more aligned with a 4 than with any other vector a i , i 6 = 4. This result demonstrates that the theoretical analysis of Sec. III can be generalized to the case of data with more than two dimensions (q > 2). For comparison, we also show the results obtained by L2-PCA. V. CONCLUSIONS This work has found a relationship between L1-PCA and the principal axes defining the dominant signal subspace when the data are not sphered. This result is consistent with the fact that the criterion seeks for projections with large amplitude. The dominant source subspace identified by L1-PCA proves more robust in the presence of outliers than classical L2-PCA. Further work should aim at the theoretical characterization of the influence of outliers and the extension of our analysis to more than two underlying dimensions.

Column number

Fig. 1 .Fig. 2 .

 12 Fig.1. Scatter plot of a two-dimensional observation (p = 2) composed of two independent components (q = 2) and N = 200 samples. The independent component direction a 1 is represented by a solid line arrow. The L1-PCA solution w L1 (bold solid-line arrow) is aligned with the main axis of the data, as predicted by our analysis. For the purposes of comparison, the L2-PCA solution w L2 is also shown (dashed-line arrow). Vectors are not drawn to scale. Top plot: No outliers. Bottom plot: 5% of outliers. Outliers are drawn at random from a normalized Gaussian distribution. It is apparent in this figure that L1-PCA is more robust to outliers than L2-PCA.

Fig. 3 .

 3 Fig.3. Results with q = 7 random directions, where a 4 has the strongest magnitude, and N = 500 samples. Each pair of vertical lines represents the confidence intervals of the angles formed by L1-PCA vector w L1 (left line) and the L2-PCA vector w L2 (right line) with each data direction over 100 independent Monte Carlo realizations. L1-PCA succeeds in finding the principal direction naturally defined by the dominant axis a 4 . The results are also similar to those obtained by applying L2-PCA.

  Proc. CAMSAP, 7th IEEE Intl. Workshop on Computational Advances in Multi-Sensor Adaptive Processing, Curaçao, Dutch Antilles, Dec.10-13, 2017.

in:

ACKNOWLEDGMENTS

This work is partially funded by the Spanish Ministry of Economy and Competitiveness under project TEC2014-53103-P CMANS.

V. Zarzoso is a member of the Institut Universitaire de

France.


Alexandru Dimca
email: dimca@unice.fr

Polynomial

S Barland

P Reynaud-Bouret

F Cauneau

; K Guillouzouic

U Kuhl

T Passot

F Planchon

Polynomial interpolation in higher dimensions

Keywords:

f (p i ) = P(p i ) for all i = 1, 2, ..., n.

The hope is that this polynomial will in fact satisfy f (t) = P(t) for any real number t, and hence our discovered formula would allow us to make predictions as well. Fig. 1 The n = 4 points (p i , c i ) are on the parabola y = P(x), with P(x) = x 2 -3x + 1.

The following result was first published by Waring in 1779, rediscovered by Euler in 1783, and published by Lagrange in 1795.

Theorem 1. For any n distinct real numbers p i , i = 1, 2, ..., n, and any given n values c i , i = 1, 2, ..., n, there is a unique polynomial P(x) of degree at most n -1, such that P(p i ) = c i for all i = 1, 2, ..., n.

To give a formula for P(x), consider, for any i = 1, 2, ..., n, the degree n -1 polynomial

Q i (x) = Π j=1,n;j =i (x -p j ) Π j=1,n;j =i (p i -p j ) ,
and note that Q i (p i ) = 1 and Q i (p k ) = 0 for any k = i. With this notation one has

P(x) = ∑ i=1,n c i Q i (x).
Consider the vector space of polynomials of degree at most d, denoted by R[x] ≤d , and the linear map given by evaluation

d : R[x] ≤d → R S = R n , (Q)(p i ) = Q(p i ),
for any i = 1, 2, ..., n, where R S denotes the vector space of all functions S → R.

Interpolation in dimension 2

Let now p i = (x i , y i ) for i = 1, 2, ..., n be n points in the plane R 2 . Assume each point has an associated value c i ∈ R, thought of as the result of a measurement at the point p i . Let S = {p 1 , p 2 , ..., p n } be the set of these points, consider the associated function f : S → R given by f (p i ) = c i and look for the minimal degree d such that there is a polynomial

Q ∈ R[x, y] of degree d satisfying f (p i ) = Q(x i , y i ) for all i = 1, 2, ..., n.
This is the old question, but the setting is new: the answer now depends on the position of the points p i in the plane.

Example 1 (3 points in the plane). If the 3 points p 1 , p 2 and p 3 are not collinear, then we can take d = 1. Indeed, as above, we can construct 3 polynomials Q 1 , Q 2 and Q 3 , by taking the equations of lines passing through two of the points p i . When the 3 points p 1 , p 2 and p 3 are collinear, then the minimal degree is d = 2. Indeed, there are conics passing through two of these points and avoiding the remaining one.

We discuss now a special type of interpolation node, i.e., a special class of choices for the points p i 's. Consider a finite family of lines L j : j (x, y) = a j x + b j y + c j = 0 in the plane R 2 , for j = 1, 2, ..., m. If these lines are generic, i.e., no two are parallel and no three are concurrent, then we get precisely

N = m 2
intersection points, which will play the role of our points p i .

Theorem 2. For any m generic lines in the plane and any given N values c i associated to their intersection points p i , there is a unique polynomial P(x, y) of degree at most m -2, such that P(p i ) = c i for all i = 1, 2, ..., N. More precisely, the evaluation map Note that the degree of P is much smaller than the number of interpolation points p i , namely m -2 < N = m(m -1)/2. The injectivity claim is easy, using Bezout Theorem about the intersection of two plane curves. The surjectivity is subtle, the proof uses Hodge theory, see [START_REF] Dimca | Koszul complexes and pole order filtrations[END_REF]. The case when C is a Chebyshev curve is particularly interesting, see [START_REF] Dimca | Chebyshev curves, free resolutions and rational curve arrangements[END_REF].

d : R[x, y] ≤d → R S = R N , (Q)(p i ) = Q(x i , y i ), p i = (x i , y i )

Projective Duality and Interpolation

From now on we move from the real field R to the complex field C, and from the affine plane C 2 to the projective plane P 2 , with coordinates (x : y : z). A point p in P 2 is given by 3 homogeneous coordinates

p = (a : b : c),
where a, b, c ∈ C are not all zero. To such a point we can associate a line L p in P 2 , given by the equation

L p : ax + by + cz = 0.
Hence to a set of points S = {p i : i = 1, 2, ..., n} in P 2 , we can associate a line arrangement A S = {L p i : i = 1, 2, ..., n} in P 2 . The multiplicity of a point p in a line arrangement A is the number of lines of the arrangement A passing through p. For more on line arrangements we refer to [START_REF] Dimca | Hyperplane Arrangements: An Introduction[END_REF].

Splitting type of a line arrangement

Why to pass from points to lines ? Because line arrangements have a lot of geometry. In particular, for any line arrangement A in P 2 one can define a rank two vector bundle E = E(A) on the projective plane P 2 , the bundle of logarithmic vector fields along A. If L is a generic line in P 2 , then the restriction E|L splits as a direct sum of two line bundles on L = P 1 , with first Chern classes given by two negative integers, say (-a, -b), with 0 ≤ a ≤ b.

The pair (a, b) is called the splitting type of E and satisfies a + b = |A| -1.

For details, see [START_REF] Abe | On the splitting types of bundles of logarithmic vector fields along plane curves[END_REF][START_REF] Dimca | Syzygies and logarithmic vector fields along plane curves[END_REF].

A new look at the 1-dimensional case: a refinement

The fact that the evaluation map

d : R[x] ≤d → R S = R n , (Q)(p i ) = Q(p i ), is surjective for d ≥ n -1 is equivalent to the claim that dim ker d = dim{Q ∈ R[x] ≤d : Q(p i ) = 0 for any i} = d + 1 -n,
for d ≥ n -1. Now fix an integer k ≥ 1, consider a new point q ∈ R, but q / ∈ S, and define a new evaluation map

d,q,k : R[x] ≤d → R n+k = R n × R[x] ≤k-1 , where (Q)(p i ) = Q(p i ) ∈ R and (Q)(q) = T k-1 Q(q) ∈ R[x] ≤k-1 is the (k - 1)
-st Taylor expansion of the polynomial Q at the point q.

In particular, (Q)(q) = T k-1 Q(q) = 0 if and only if the first (k -1) derivatives of Q vanish at q, namely Q (j) (q) = 0 for all 0 ≤ j ≤ k -1. It is easy to show that this new evaluation map is surjective for d ≥ n + k -1, and hence

dim ker d,q,k = d + 1 -n -k, for d ≥ n + k -1.
The practical interest of this refinement is that, for instance, a zero Taylor expansion of high order means very small values for the polynomial Q in the neighborhood of the given point q.

A new 2-dimensional interpolation problem

Starting with a given set of points S = {p i : i = 1, 2, ..., n} in P 2 , we fix an integer k ≥ 1, and define I(S) d to be the vector space of homogeneous polynomials Q ∈ C[x, y, z] d such that Q(p i ) = 0 for any i = 1, 2, ..., n. Choose then a generic point q ∈ P 2 and consider the vector space V(d, S, k, q) of homogeneous polynomials Q ∈ I(S) d such that T k-1 Q(q) = 0. The expected dimension of this vector space is dim e V(d, S, k, q) := dim

I(S) d - k + 1 2 ,
when this number is positive. An important special case is when d = k, which is also the simplest case to consider. In this setting, we introduce the following notion, see [START_REF] Cook | Line arrangements and configurations of points with an unexpected geometric property[END_REF].

Definition 1. We say that the set S admits an unexpected curve of degree k if dim V(k, S, k, q) > dim e V(k, S, k, q) ≥ 0.

The main result in this direction is the following, see [START_REF] Cook | Line arrangements and configurations of points with an unexpected geometric property[END_REF].

Theorem 3. Let S be a finite set of n points in P 2 and let (a S , b S ) be the splitting type of the dual line arrangement A S . Then S admits an unexpected curve of degree k if and only if the following hold.

• a S + 1 ≤ k ≤ b S -1;
• the multiplicity of any intersection point in A S is at most a S + 1.

An example: the complete polygonal arrangements

Consider a regular polygon with N ≥ 3 edges, and the associated line arrangement A consisting of the following 2N + 1 lines

• the N lines determined by the N edges of the polygon,

• the N symmetry axes of the polygon, and

• the line at infinity. This type of line arrangement occurs in the following result, see [START_REF] Di Marca | Unexpected curves arising from special line arrangements[END_REF].

Theorem 4. For N even, the complete N-polygonal arrangement has an unexpected curve of degree N. The proof uses the theory of supersolvable line arrangements to show that a = N -1 and b = N + 1. As an example, for N = 4 we get a = 3 and b = 5. Hence A admits an unexpected curve of degree 4 by the result in see [START_REF] Cook | Line arrangements and configurations of points with an unexpected geometric property[END_REF].

In the case N = 4, the complete N-polygonal arrangement is the dual arrangement A S , where the set of points S consists of the points (0 : 0 : 1), (0 : 1 : 0), (1 : 0 : 0), (1 : 1 : 1), (0 : 1 : 1), (1 : 0 : 1), (1 : 1 : 0), (-1 : 1 : 0), (1 : 1 : 2). This situation was considered first by B. Harbourne in [START_REF] Harbourne | Asymptotics of linear systems, with connections to line arrangements[END_REF].

(p 4 ,

 4 c 4 ) = (3, 1)

Fig. 2 m

 2 Fig. 2 m = 4 generic lines in the plane

Fig. 3

 3 Fig. 3 Complete 4-polygonal arrangement; the line at infinity is drawn as the red circle.

  The above results say that d is injective if and only if d ≤ n -1, and d is surjective if and only if d ≥ n -1. Such results are true over any field: R, C or even finite fields.