P M R De Oliveira

J H De

M Goulart

C A R Fernandes

Senior Member, IEEE V Zarzoso

Blind Source Separation in Persistent Atrial Fibrillation Electrocardiograms Using Block-Term Tensor Decomposition with Löwner Constraints

Keywords: Atrial Fibrillation, Block Term Decomposition, Constrained Alternating Group Lasso, Electrocardiogram, Löwner Matrices. I

The estimation of the atrial activity (AA) signal from electrocardiogram (ECG) recordings is an important step in the noninvasive analysis of atrial fibrillation (AF), the most common sustained arrhythmia encountered in clinical practice. This problem admits a blind source separation (BSS) formulation that has been recently posed as a tensor factorization, using the Hankelbased block term decomposition (BTD), which is particularly well-suited to the estimation of exponential models like AA during AF. However, persistent forms of AF are characterized by short R-R intervals and very disorganized (or weak) AA, making it difficult to model AA directly and perform its successful extraction through Hankel-BTD. To overcome this drawback, the present work proposes a tensor approach to estimate QRS complexes and subtract them from the ECG, resulting in a signal that, ideally, only contains the AA component. Such an approach tackles the problem of blind separation of rational functions, which models QRS complexes explicitly. The data tensor admitting a BTD is built from Löwner matrices generated from each lead of the observed ECG. To this end, this paper formulates a variant of the recently proposed constrained alternating group lasso (CAGL) algorithm that imposes Löwner structure on the decomposition blocks. This is done by performing an orthogonal projection, which we explicitly derive, at each iteration of CAGL. Results from experiments with synthetic data show the consistency of the proposed Löwner-constrained AGL (LCAGL) in extracting the desired sources. Experimental results obtained on a population of 20 patients suffering from persistent AF show that the proposed variant outperforms other tensor-based methods in terms of atrial signal estimation quality from ECG records as short as a single heartbeat.

distinguished by an uncoordinated and irregular atrial depolarization [START_REF] January | AHA/ACC/HRS Focused Update of the 2014 AHA/ACC/HRS Guideline for the Management of Patients With Atrial Fibrillation[END_REF]. Persistent AF represents a particularly complex form of this arrhythmia, where extensive atrial remodelling has taken place due to sustained AF, significantly affecting atrial activity (AA) and causing AF perpetuation. This cardiac rhythm disturbance is considered an economical burden, as a patient suffering from AF costs, anually, approximately $8 700 more to healthcare providers than a patient without AF in the USA. It is also associated with a higher morbidity, representing a major health concern, since about half million hospitalizations with AF as the primary diagnosis are registered every year, also in USA [START_REF] January | AHA/ACC/HRS Focused Update of the 2014 AHA/ACC/HRS Guideline for the Management of Patients With Atrial Fibrillation[END_REF]. If the number of new diagnosis stays stable, it is estimated that AF will increase from 2.3 to at least 12.1 million patients in the USA by 2050 [START_REF] Mainardi | Understanding Atrial Fibrillation: The Signal Processing Contribution[END_REF], and from 8.8 to about 17.9 million in the European Union by 2060 [START_REF] Krijthe | Projections on the number of individuals with atrial fibrillation in the European Union, from 2000 to 2060[END_REF], becoming then a new epidemic. Worldwide, the number of AF patients was estimated to be around 33.5 million in 2010 [START_REF] Morillo | Atrial fibrillation: the current epidemic[END_REF]. Considering these worrisome statistics and the fact that the electrophysiological mechanisms of AF are not totally understood, it is not surprising that this challenging cardiac condition has been a topic of intense research in the past few years, and is expected to be investigated even more deeply.

During AF episodes, the atrial depolarization, represented by the P wave during normal sinus rhythm, is replaced by fibrillatory waves, or f waves, which are present continuously throughout the electrocardiogram (ECG) recording. The f waves are masked by the QRS-T complex that corresponds to the ventricular depolarization and repolarization, i.e., the ventricular activity (VA), in each heartbeat. Figure 1 illustrates an example of a persistent AF ECG recording with short R-R intervals and fine AF, i.e., when the AA has a very low amplitude. These two characteristics are quite common in persistent AF ECGs.

In order to provide a better understanding of the complex mechanisms behind AF, a detailed analysis of the AA signal is necessary. A cost-effective noninvasive method to perform this task is to extract the AA from the standard 12-lead ECG. During AF, the AA and the VA signals are typically assumed uncoupled, which makes the extraction of AA from the ECG admit a blind source separation (BSS) formulation [START_REF] Rieta | Atrial activity extraction for atrial fibrillation analysis using blind source separation[END_REF]. The block-term decomposition (BTD) built from Hankel matrices, proposed as a technique to solve BSS problems in [START_REF] Lathauwer | Blind separation of exponential polynomials and the decomposition of a tensor in rank-(lr, lr, 1) terms[END_REF], was used to noninvasively extract the AA signal from AF ECG recordings, and shown to outperform the matrix-based techniques in this particular application [START_REF] Zarzoso | Parameter estimation in block term decomposition for noninvasive atrial fibrillation analysis[END_REF]- [START_REF] De Oliveira | Block term decomposition of ECG recordings for atrial fibrillation analysis: temporal and inter-patient variability[END_REF]. Also, tensorbased techniques do not require orthogonality or mutual inde- pendence between the sources, enjoying essential uniqueness under mild constraints, and can perform BSS from rather short signal records, on a beat-to-beat basis.

The BTD tensor model is flexible, as it can accommodate sources with different characteristics, i.e., with diferent number of poles (or degrees, in the case of rational functions). Two commonly used methods for computing this decomposition are the nonlinear least squares (NLS) method implemented in the Tensorlab MATLAB Toolbox [12] and the alternating least squares method with enhanced line search (ALS-ELS) [START_REF] Lathauwer | Decompositions of a higher-order tensor in block terms -Part III: Alternating least squares algorithms[END_REF]. However, their performance depends considerably on their initialization and model structure, i.e., the rank of the matrix factors and the number of blocks of the tensor.

To overcome these limitations, an algorithm to compute the BTD called alternating group lasso (AGL) and a variant called constrained AGL (CAGL) that deals with linear constraints over block matrices were recently proposed [START_REF] De | Alternating group lasso for block-term tensor decomposition with application to ECG source separation[END_REF]. The principle of AGL is to compute BTD as the minimization of a least squares (LS) fit term plus a regularization term in order to penalize the rank of the blocks. AGL is able to find an approximate BTD with weak assumptions on its model parameters (number of block terms and multilinear ranks), since only an upper bound is imposed on these parameters. In CAGL, approximate projections are performed at each iteration to ensure that the block matrices have low rank and belong to a specified subspace, a problem known as structured low-rank approximation (SLRA) [START_REF] Markovsky | Structured low-rank approximation and its applications[END_REF]. Indeed, the estimation of low-rank block matrices is a crucial condition to extract the desired information contained in the (unique) model parameters. CAGL was also assessed as an AA extraction tool in [START_REF] De | Alternating group lasso for block-term tensor decomposition with application to ECG source separation[END_REF], showing its ability to extract more information of the AF ECG signal than other existing techniques, as this method could extract two possible atrial source signals, probably having their origin in different parts of the atria.

The Hankel-BTD has proven to be a useful and powerful AA extraction tool in AF analysis, but a satisfactory performance is only achieved for segments with long R-R intervals and with well defined AA, visible in most of the segment. However, as stated earlier in this section, recordings with short R-R intervals (< 0.75 s) and fine AF (amplitude of the f waves lower than 0.1 mV) are quite more common during persistent AF episodes. When assessed in these challenging cases, the Hankel-BTD as well as the matrix-based techniques do not provide satisfactory results [START_REF] De Oliveira | Löwner-based tensor decomposition for blind source separation in atrial fibrillation ECGs[END_REF]. Automatically selecting the AA signal among the estimated sources after performing BSS is also an unsolved issue, since no systematic method is reported in the literature for this purpose.

Aiming to avoid such limitations in these common and challenging persistent AF scenarios, the present work proposes the BTD built from Löwner matrices as a solution for BSS of rational functions [START_REF] Debals | Löwner-based blind signal separation of rational functions with applications[END_REF] to model the VA and separate it from the AA. This strategy suits the characteristics of VA in ECG recordings, since the QRS complex can be well approximated by rational functions [START_REF] Lócsi | Approximating poles of complex rational functions[END_REF], [START_REF] Fridli | Rational functions systems in ECG processing[END_REF], and when mapped onto Löwner matrices, the degree of the rational function matches the rank of the Löwner matrix [START_REF] Debals | Löwner-based blind signal separation of rational functions with applications[END_REF]. Modeling VA instead of AA is a reasonable strategy in these difficult cases (weak AA and/or short R-R intervals), common in persistent AF episodes, as the AA signal becomes very difficult to model [START_REF] De Oliveira | Löwner-based tensor decomposition for blind source separation in atrial fibrillation ECGs[END_REF]. The VA estimated by the Löwner-BTD is then subtracted from the ECG, resulting in a signal that contains mainly AA. Another advantage of the proposed method is that no technique for atrial source selection is needed, since the resulting signal is already the desired signal.

The main contributions of the paper can be summarized as follows:

• We put forward a new tensor-based method for noninvasive AA estimation in the challenging scenarios of low amplitude f-waves or short R-R intervals, typical of persistent AF episodes. This work is a significant extension of the preliminary results reported in the conference paper [START_REF] De Oliveira | Löwner-based tensor decomposition for blind source separation in atrial fibrillation ECGs[END_REF]. Indeed, while the NLS method is used to compute the Löwner-BTD and tested in a reduced database in [START_REF] De Oliveira | Löwner-based tensor decomposition for blind source separation in atrial fibrillation ECGs[END_REF], we provide herein an improved method for Löwner-BTD computation and statistically more significant results by doubling the size of the AF patient database. This paper formulates a variant of CAGL that deals with Löwner constraints, imposing a low-rank Löwner structure. We show how this constraint can be enforced by an alternating projection approach, for which we derive an explicit expression of the projection onto the subspace of Löwner matrices. The resulting LCAGL is compared to the Löwner-BTD computed with the NLS method (used in the experiments reported in [START_REF] De Oliveira | Löwner-based tensor decomposition for blind source separation in atrial fibrillation ECGs[END_REF]), as well as the recently proposed CAGL dealing with Hankel constraints. Experimental results with synthetic and real AF ECG segments show the superiority of the proposed technique in challenging cases of this complex arrhythmia. Successful AA estimation is achieved from observation records as short as a single heartbeat, making the new method potentially very attractive in real-time AF analysis.

The rest of this paper is organized as follows. Section II presents the notation and some basic definitions of tensor algebra. Section III recalls BTD as a tensor approach to solve BSS problems and how to build a tensor from an ECG data matrix. Section IV describes AGL for computing BTD, while Section V formulates its variant with Löwner constraints and derives the linear projection onto the Löwner subspace. A numerical evaluation of the proposed approach is shown in Section VI with synthetic signals. Section VII presents the database setup and results of experiments with real AF data. Section VIII discusses the advantages and limitations of the present work and, finally, Section IX concludes the paper and provides some prospects of future research.

II. NOTATION AND ALGEBRA PREREQUISITES

To ease the comprehension of the paper, this section summarizes the notation used in the present work and introduces some definitions that are not common in research fields outside tensor algebra.

Scalars, vectors, matrices and tensors are represented by lower-case (a, b, ...), boldface lower-case (a, b, ...), boldface capital (A, B, ...) and calligraphic (A, B, ...) letters, respectively. The matrix transpose operator is represented by (•)

T , ||•|| F is the Frobenius norm and • is the outer product. Diag(a) builds a diagonal matrix by placing the elements of a along the diagonal and I N represents an identity matrix of size N ×N . A third-order tensor A ∈ C I1×I2×I3 , with scalar entries a i1,i2,i3 , has its frontal slices represented by

A ..i3 ∈ C I1×I2 . A matrix A ∈ C I1×I2 ,
||A|| 2,1 = I2 i2=1 ||a .i2 || . (1)
Given a vector y of length N , one can build a Hankel matrix H ∈ R I×J , where

I = J = N +1 2 if N is odd or I = N 2 and J = N 2 + 1 if N is even, with entries: h i,j y i+j-1 (2)
with i = 1, ..., I, j = 1, ..., J.

One can also build a Löwner matrix L ∈ R I×J from y, where I + J = N , with entries:

i,j y xi -y zj x i -z j (3)
where i = 1, ..., I, j = 1, ..., J and y tn is the t th n time sample of y, sampled in the time set T = {t 1 , t 2 , ..., t N }, which is partitioned into two different non-overlapping time sets: X = {x 1 , x 2 , ..., x I } and Z = {z 1 , z 2 , ..., z J }. Two simple partitioning methods are the interleaved partition, i.e., X = {t 1 , t 3 , t 5 , ...}, and Z = {t 2 , t 4 , t 6 , ...}, and the block partition, i.e., X = {t 1 , t 2 , ..., t I }, and Z = {t I+1 , t I+2 , ..., t N }. As stated in [START_REF] Debals | Löwner-based blind signal separation of rational functions with applications[END_REF], both partitioning methods give similar performance in a BSS scenario. III. TENSOR-BASED BSS OF ECG DATA This section discusses the tensor-based BSS formulation and its connection with AF ECGs. First, Section III.A describes the matrix-based BSS formulation with its constraints to guarantee uniqueness, motivating the formulation of a tensor-based BSS problem with uniqueness under mild conditions. The Hankel-BTD for modeling the AA is recalled in Section III.B, whereas its Löwner-constrained version to model the VA is introduced in Section III.C.

A. BSS Formulation

An ECG recording composed of K leads and N time samples can be modeled as a matrix factorization:

Y = MS ∈ R K×N (4)
where M ∈ R K×R is the mixing matrix, modeling the propagation of the cardiac electrical sources from the heart to the body surface, S ∈ R R×N is the source matrix that contains the atrial and ventricular sources and R is the number of sources. Since the goal is to estimate M and S from matrix Y only, one can see that source estimation in ECG recordings is a BSS problem.

As previously stated, matrix-based methods require some constraints in order to guarantee the uniqueness of the factorization in (4). For example, principal component analysis (PCA) [START_REF] Jolliffe | Principal Component Analysis[END_REF] requires the columns of the mixing matrix M to be orthogonal, whereas independent component analysis (ICA) [START_REF] Comon | Independent component analysis, A new concept?[END_REF] methods require statistical independence between the sources and non-Gaussianity. Although mathematically coherent, such constraints may lack physiological grounds, hindering the interpretation of the results. By contrast, tensorbased techniques do not require such constraints and can ensure uniqueness under relaxed conditions.

As it will be detailed later, the BTD with constrained matrix factors can well suit the characteristics of ECG signals, which makes this tensor factorization technique a powerful tool for AA extraction from AF ECGs. The BTD of an arbitrary thirdorder tensor T ∈ R I×J×K is written as:

T = R r=1 E r • c r ( 5
)
where R is the number of blocks, c r ∈ R K is nonzero and E r ∈ R I×J is a structured matrix that has rank L r , admitting the factorization E r = A r B T r , where A r ∈ R I×Lr and B r ∈ R J×Lr have rank L r . We may then rewrite (5) as: One can see that the BTD is a decomposition of T in multilinear rank-(L r ,L r ,1) terms [START_REF] Lathauwer | Blind separation of exponential polynomials and the decomposition of a tensor in rank-(lr, lr, 1) terms[END_REF], represented by a sum of the outer product of its matrix and vector factors, as illustrated in Figure 1. Note that the number of sources R equals the number of blocks in the BTD model ( 5)-( 6). Several results can be invoked to guarantee the essential uniqueness of this tensor decomposition. For example, in [6, Theorem 2.2], it is shown that the BTD is essentially unique if the following conditions are satisfied:

T = R r=1 A r B T r • c r . (6)
1) The matrix factors

A = A 1 A 2 . . . A R ∈ R I× R r=1 Lr and B = B 1 B 2 . . . B R ∈ R J× R r=1 Lr are full-column rank, which requires R r=1 L r ≤ min(I, J). 2) Matrix C = c 1 c 2 . . . c R ∈ R K×R does not
contain proportional columns. For the processed ECG segments, we typically have that min(I, J) ≈ N +1 2 ≥ R r=1 L r . Also, matrix C corresponds to the mixing matrix M, whose columns represent the contribution of each source to the ECG leads, which are unlikely to be proportional. Essentially unique means that the structure of the decomposition does not change if the r th 1 and r th 2 terms are permuted, as long as they have equal ranks, and if E r is rescaled, as long as c r is counterscaled. Milder uniqueness conditions can be found in [START_REF] Lathauwer | Blind separation of exponential polynomials and the decomposition of a tensor in rank-(lr, lr, 1) terms[END_REF].

If the matrices E r are linearly structured in the sense that they belong to a given subspace U S , then the frontal slices of T , given by T ..k , also belong to U S . For the BSS application of this decomposition, two commonly used structures are the Hankel structure, that focuses on separating complex exponential functions, and the Löwner structure, focusing on the separation of rational functions. In the following, we will describe the principles behind these choices.

B. AA Modeling via Exponential Functions

The Hankel-BTD is especially adapted to source signals composed of a sum of complex exponentials:

s(t) = F f =1 λ f z t f (7)
where z f ∈ C is a possibly complex pole, λ f the associated scale factor (that can also be complex-valued), 1 ≤ f ≤ F , and t denotes continuous time. This model is indeed a good representation of the quasi-periodic behaviour of the AA signal in AF episodes [START_REF] Stridh | Spatiotemporal QRST cancellation techniques for analysis of atrial fibrillation[END_REF].

It is shown in [START_REF] Lathauwer | Blind separation of exponential polynomials and the decomposition of a tensor in rank-(lr, lr, 1) terms[END_REF] that a BTD built from Hankel matrices solves the BSS problem (4) when the underlying sources follow the exponential model [START_REF] Comon | Independent component analysis, A new concept?[END_REF]. This results relies on the fact that a Hankel matrix associated with an F -pole exponential model has a rank bounded by F , which follows from the Vandermonde decomposition of the Hankel matrix. To exploit this idea to perform BSS, a tensor is first obtained from the data matrix Y by mapping each of its K rows into a Hankel matrix

H (k) ∈ R I×J , 1 ≤ k ≤ K, according to (2).
Next, the tensor is built by stacking each Hankel matrix along the third dimension, as frontal slices, of a third-order tensor Y H ∈ R I×J×K , that is:

Y ..k = H (k) = R r=1 m k,r H (r) S (8)
where

H (r)
S is a Hankel matrix built from the r th source of S. Because the described mapping is linear and each row of Y is a linear combination of the sources with weights determined by one row of M, it follows that the third-order tensor Y H admits a BTD model and can be written as:

Y H = R r=1 H (r) S • m .r . (9)
The Hankel-BTD suits the characteristics of AA in AF episodes since atrial signals can be approximated by a sum of complex exponentials like [START_REF] Comon | Independent component analysis, A new concept?[END_REF] and mapped into Hankel matrices with rank equal to the number of poles [START_REF] Zarzoso | Parameter estimation in block term decomposition for noninvasive atrial fibrillation analysis[END_REF]. However, in persistent AF, ECG recordings present short R-R intervals, significantly masking the AA in the recording. Also, f-waves during AF are characterized for being disorganized with very low amplitude, sometimes lower than the noise. Such issues make it difficult to fit the AA signals to this exponential model, compromising the AA extraction via Hankel-BTD.

C. VA Modeling via Rational Functions

To circumvent the aforementioned difficulty in the persistent AF scenarios considered in this work, we turn to another alternative: VA modeling. Indeed, the QRS morphology is not significantly changed during AF, as compared to normal sinus rhythm (NSR), and hence it is reasonable to focus on VA estimation. When subtracting the VA estimate from the ECG, the resulting signal will ideally contain the AA signal.

Previous works notice that QRS complexes can be approximated by rational functions of low degree [START_REF] Debals | Löwner-based blind signal separation of rational functions with applications[END_REF], [START_REF] Lócsi | Approximating poles of complex rational functions[END_REF], [START_REF] Fridli | Rational functions systems in ECG processing[END_REF]:

s(t) = a(t) + F f =1 D f d=1 c f,d (t -p f ) d (10)
where a(t) is a polynomial of degree A, F is the number of complex poles, denoted p f , D f is their multiplicity, and

c f,d = 1/(u f -v d )
are the scalar entries of a Cauchy matrix based on the vectors u ∈ C F and v ∈ C D f , with

u f = v d , ∀f, d.
Interestingly, a second kind of BTD suits the characteristics of signal model [START_REF] De Oliveira | Temporal stability of block term decomposition in noninvasive atrial fibrillation analysis[END_REF]. Instead of Hankel matrices, this alternative BTD is based on Löwner matrices built as in [START_REF] Krijthe | Projections on the number of individuals with atrial fibrillation in the European Union, from 2000 to 2060[END_REF]. The key result underying this decomposition is that when a signal following model ( 10) is mapped into a Löwner matrix, the degree of the rational function matches the rank of the matrix [START_REF] Debals | Löwner-based blind signal separation of rational functions with applications[END_REF]. Hence, the idea is to map each row of the data matrix Y into a Löwner matrix L (k) ∈ R I×J as in ( 6), 1 ≤ k ≤ K. Next, the tensor is built by stacking each Löwner matrix along the third dimension, as frontal slices, of a third-order tensor Y L ∈ R I×J×K . Similarly as with Hankel matrices, one can show that:

Y ..k = L (k) = R r=1 m k,r L (r) S (11)
where

L (r) S
∈ R I×J is a Löwner matrix built from the r th row of S. One can see that the procedure for constructing Y ..k in ( 11) is a linear mapping and, for each r, the outer product between matrix L (r) S and the r th column of M, i.e., m .r , is performed to build a third-order tensor containing the contribution of the r th source to the ECG tensor. Putting together the contribution of all sources, the third-order tensor Y L admits a BTD model and can be written as:

Y L = R r=1 L (r) S • m .r . ( 12
)
In short, the Löwner-BTD tensor model ( 12) is directly related to the rational function model [START_REF] De Oliveira | Temporal stability of block term decomposition in noninvasive atrial fibrillation analysis[END_REF], and is therefore well-suited to VA estimation in AF ECGs. To do so, one first needs to estimate the BTD model factors from third-order tensor Y L , an important step that is addressed next.

IV. ALTERNATING GROUP LASSO ALGORITHM FOR BTD

In this section, the classical approach for the computation of the BTD model factors and its limitations are recalled. Then, we review the AGL algorithm, which was conceived to overcome such limitations.

In general, an approximate BTD is computed by minimizing the Euclidean distance between the observed data tensor Y ∈ C I×J×K and a model of fixed structure with respect to the model components:

f (A, B, C) Y - R r=1 A r B T r • c r 2 F . (13)
In the special case of interest shown in eqn. (12), L r must belong to the subspace of Löwner matrices with dimensions

( N 2 × N 2 )
, denoted S L . The mode-3 slices Y ..k of the observed tensor are Löwner by construction. However, a solution ( Â, B, Ĉ) of ( 13) may not satisfy Âr BT r ∈ S L , due to noise and modeling imperfections. Note that even if the sum R r=1 ( Âr BT r )ĉ k,r is Löwner, there is no guarantee that so are matrices Âr BT r , r = 1, 2, ..., R. Also, algorithms based on [START_REF] Lathauwer | Decompositions of a higher-order tensor in block terms -Part III: Alternating least squares algorithms[END_REF] are strongly dependent of the initialization of their matrix factors and do not estimate the model parameters, i.e., the number of blocks and their ranks.

To overcome such limitations, instead of using a fixed BTD structure as [START_REF] Lathauwer | Decompositions of a higher-order tensor in block terms -Part III: Alternating least squares algorithms[END_REF], an algorithm called AGL and its constrained version for Hankel matrices called CAGL, are introduced in [START_REF] De | Alternating group lasso for block-term tensor decomposition with application to ECG source separation[END_REF]. This method includes penalization terms promoting lowrank blocks and controlling the number of blocks as:

F (A, B, C) f (A, B, C) + γ g(A, B, C) (14)
where f (•, •, •) is the same as in [START_REF] Lathauwer | Decompositions of a higher-order tensor in block terms -Part III: Alternating least squares algorithms[END_REF], γ > 0 is a regularization parameter and g is a regularization function of the form:

g(A, B, C) A 2,1 + B 2,1 + C 2,1 . (15)
Due to the geometric properties of the mixed 2,1 -norm, solutions where A, B and C have null columns (for sufficiently high γ values) will be induced, allowing one to select the relevant low-rank blocks. This method is called group lasso and is a generalization of the the lasso estimator principle [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF].

Since problem ( 14) is nonconvex (and nonsmooth), but convex by blocks, a block coordinate descent (BCD) approach is employed in [START_REF] De | Alternating group lasso for block-term tensor decomposition with application to ECG source separation[END_REF]. BCD consists in partitioning the set of optimization variables and sequentially solving convex subproblems in each subset of variables, fixing the others.

Consider Â(ν) , B(ν) and Ĉ(ν) the estimates of A, B and C, respectively, obtained at iteration ν. Fixing B = B(ν) and C = Ĉ(ν) in ( 14) the subproblem in A of iteration (ν + 1) can be written as [START_REF] De | Alternating group lasso for block-term tensor decomposition with application to ECG source separation[END_REF]:

min A 1 2 Y-W (ν+1) A (A) 2 F +γ A 2,1 + τ 2 A-Â(ν) 2 F ( 16
)
where W (ν+1)

A is a linear map that depends on B(ν) and Ĉ(ν) and τ is positive.

Analogously, strictly convex subproblems can be derived for B and C. AGL solves these subproblems alternatively, through the updates of Â(ν) , B(ν) and Ĉ(ν) , in this order, at each iteration ν. When all subsets are updated, one iteration of the algorithm is completed. Iterations are repeated until convergence.

V. L ÖWNER-CONSTRAINED ALTERNATING GROUP LASSO

In order to ensure the Löwner structure of the matrix factors at the end of iterations, this section explains how Löwner constraints can be imposed in CAGL, and then details the linear projection that can be implemented in SLRA algorithms.

A. Imposing Löwner Constraints in L r

We propose a formulation to deal with Löwner constraints in CAGL, called here Löwner-constrained AGL (LCAGL). We highlight that LCAGL follows the same principle of CAGL, which is to ensure a specified structure over the block matrices. However, instead of using a projection on the Hankel subspace as described in [START_REF] De | Alternating group lasso for block-term tensor decomposition with application to ECG source separation[END_REF], LCAGL uses an optimal projection on the Löwner subspace that will be explicitly given later.

Löwner constraints can be enforced by solving a SLRA problem for each block matrix at each iteration of LCAGL. For this purpose, Cadzow's Algorithm (CA) [START_REF] Cadzow | Signal enhancement-a composite property mapping algorithm[END_REF] is used at the end of each iteration. CA performs alternating projections, where the one that leads to the set of matrices with rank up to L r is performed by the truncated singular value decomposition (SVD), whereas the one that leads to the Löwner subspace is formulated below. It is important to highlight that, in practice, the rank used in the SVD truncation step is computed by counting the current number of simultaneously nonzero columns of A r and B r . The algorithm is summarized in Table I.

Inputs:

Data tensor Y, penalty parameter γ, proximal term weight τ , sampling period ∆, initialization ( Â(0) , B(0) , Ĉ(0) ) Outputs:

Approximate BTD factors ( Â, B, Ĉ)

1: ν ← 1 2: while stopping criteria not met do 3: Solve ( 16) to obtain Â(ν) 4:

Solve ( 16) for B to obtain B(ν) 5:

for r = 1, . . . , R do 6:

( Â(ν) r , B(ν) r ) ← Cadzow( Â(ν) r ( B(ν) r ) T , ∆) 7:
Solve ( 16) for C to obtain Ĉ(ν) 8:

ν ← ν + 1

B. Orthogonal Projection onto the Löwner Subspace

In this work, we consider for simplicity signals with an even number of samples, i.e., I = J = N/2, and the interleaved partitioning method. In this way, we can rewrite the matrix 3) as:

L (k) ∈ R N 2 × N 2 of Equation (
L (k) =     y k,t 1 -y k,t 2 t1-t2
. . .

y k,t 1 -y k,t N t1-t N . . . . . . . . . y k,t N -1 -y k,t 2 t N -1 -t2
. . .

y k,t N -1 -y k,t N t N -1 -t N     . (17)
Considering that the signals are regularly sampled with a sampling period ∆ = t n -t n-1 , for n = 2, ..., N , (17) becomes:

L (k) =     y k,t 1 -y k,t 2 -∆ . . . y k,t 1 -y k,t N -(N -1)∆ . . . . . . . . . y k,t N -1 -y k,t 2 (N -3)∆ . . . y k,t N -1 -y k,t N -∆     (18)
where each element is given by:

(k) i,j = y k,t2i-1 -y k,t2j [2(i -j) -1]∆ . ( 19
)
Note that this formulation can be easily adapted to signals that are irregularly sampled. For this other scenario, the denominator [2(i -j) -1]∆ in [START_REF] Jolliffe | Principal Component Analysis[END_REF] is replaced by ∆ ij = t 2i-1 -t 2j . Since in practice the signals are regularly sampled, the formulation of ( 19) will be kept in the sequel.

For a given matrix E (k) ∈ R I×J , its projection onto the Löwner subspace can be obtained by minimizing the LS cost function given by:

F(a, b) = N/2 i=1 N/2 j=1 e (k) i,j - a i -b j ∆[2(i -j) -1] 2 . ( 20
)
where a i = ŷk,t2i-1 and b j = ŷk,t2j , for i, j = 1, ..., N/2, with ŷk. being the estimate of y k. . Taking the derivative of F with respect to an arbitrary a p we have:

N/2 j=1 2 e (k) p,j - a p -b j ∆[2(p -j) -1] -1 ∆[2(p -j) -1] = 0 . (21)
Then,

N/2 j=1 e (k) p,j + b j ∆[2(p -j) -1] -1 ∆[2(p -j) -1] = N/2 j=1 -a p (∆[2(p -j) -1]) 2 . (22)
In this way, the solution for a p , for p = 1, ..., N/2, can be obtained as:

âp = N/2 j=1 e (k) p,j + bj ∆[2(p-j)-1] 1 ∆[2(p-j)-1] N/2 j=1 1 (∆[2(p-j)-1]) 2 . ( 23
)
Similarly, for b p we have:

bp = N/2 i=1 e (k) i,p - ai ∆[2(i-p)-1] 1 ∆[2(i-p)-1] - N/2 i=1 1 (∆[2(i-p)-1]) 2 . ( 24
)
The necessary conditions for optimality can be expressed as:

âp -α -1 p N/2 j=1 u p,j b j = z p bp + β -1 p N/2 i=1 v i,p a i = w p (25)
for p = 1, ..., N/2, where:

α p = N/2 j=1 u p,j , for u p,j = 1/{∆[2(p -j) -1]} 2 (26)
β p = N/2 i=1 v i,p , for v i,p = -1/{∆[2(i -p) -1]} 2 (27)
and

z p = N/2 j=1 y p,j ∆[2(p -j) -1] α -1 p ( 28
)
w p = N/2 i=1 y i,p ∆[2(i -p) -1] β -1 p . (29)
Hence, a system of linear equations can be solved in order to estimate â = [â 1 , ..., âN/2 ] and b = [ b1 , ..., bN/2 ]:

I N/2 -Diag(α -1 )U Diag(β -1 )V T I N/2 â b = z w (30)
where matrices U and V, and vectors α, β, z and w are constructed from u p,j , v i,p , α p , β p , z p and w p in (29), (30), (31), (32) and (33), respectively. Note that this linear system has an infinite number of solutions because the block matrix of (30) has rank N -1, as a consequence of the vertical shift invariance of the Löwner structure. Indeed, the Löwner construction is invariant to a vertical shift of the signal because the Löwner matrix built from signal y k. is equal to the one built from y k. = y k. + c, for any constant c ∈ R. To fix this ambiguity, without loss of generality, we can assume that the first sample is null, i.e., a 1 = 0, so that the linear system has now N -1 linearly independent equations in N -1 variables, thus providing now a unique solution. A matrix with Löwner structure can then be reconstructed from estimated vectors â and b.

VI. NUMERICAL EVALUATION ON SYNTHETIC DATA

This section presents some numerical results on synthetic data, which are convenient for assessing the performance of LCAGL since in these scenarios the ground-truth is known. This assessment is carried out using two different synthetic signals that simulate persistent AF recordings. Both signals contain 12 leads, generated by a random mixing matrix and 3 sources: AA, VA, and noise.

To simulate the VA signals, two synthetic QRS complexes modeled by rational functions were generated according to the model proposed in [START_REF] Bognár | Geometric interpretation of QRS complexes in ECG signals by rational functions[END_REF]:

S V A [n] = R(e -iΘ r m a [e in ]) (31)
where Θ is the coefficient parameter, responsible for the symmetric/asymmetric behavior of the model, R(•) is the real part of its complex argument, and r m a [•] is a basic normalized rational function given by:

r m a [z] := 1 -|a| 1 -āz m ( 32
)
where ā is the reciprocal of the pole a and m is a multiplicity parameter. Two synthetic QRS complex models with parameters a = 0.8, m = 2, and Θ = {π/2, 0} are generated.

To simulate the AA signal during AF, the model proposed in [START_REF] Stridh | Spatiotemporal QRST cancellation techniques for analysis of atrial fibrillation[END_REF] that mimics the f waves is used. This model with P harmonics is given by:

S AA [n] = - P p=1 x p [n] sin (p θ[n]) (33)
with modulated amplitude and phase respectively given by:

x p [n] = 2 pπ x + ∆x sin 2π fx fs n (34) and θ[n] = 2π f0 fs n + ∆f f f sin 2π f f fs n (35)
where x is the sawtooth amplitude, ∆x is the modulation peak amplitude, f x is the amplitude modulation frequency, f s is the sampling frequency, f 0 is the frequency value upon which θ[n] depends linearly, ∆f is the maximum frequency deviation and f f is the modulation frequency. Two synthetic f-wave models were generated with the parameters presented in Table II To simulate the interference typically present in an ECG, additive white Gaussian noise (AWGN) with zero-mean and variance σ 2 , respectively, is introduced. Finally, the mixing matrix is also generated according to a Gaussian distribution, with scaling factors chosen to obtain an average power ratio between the components consistent with clinical ECGs.

One generated ECG signal mimics the challenging scenario where the AA is weak, while the other synthetic signal mimics the scenario where the R-R interval is short. Both synthetic AF ECG signals are illustrated in Figure 3.

The normalized mean square error (NMSE) between the estimated and original VA source signals is computed to evaluate estimation quality. As LCAGL is relatively robust to initialization, no Monte Carlo runs were performed and the value of γ between 8 × 10 -4 and 10 -2 that provided the best solution was chosen.

Figure 4 shows the VA estimates (blue dashed line) and the respectives ground-truths (gray solid line) for the two synthetic signals along with the computed NMSE described previously and the rank of the block. The low NMSE shows that a successful source extraction is achieved while the Löwner structure of the block is guaranteed. Also, one can see that the rank of the block, rank B , increases with the number of QRS complexes in the segment, as expected, since the signal has more transient components, needing more poles to be well modeled.

VII. EXPERIMENTAL RESULTS IN REAL AF ECGS

In order to validate the proposed tensor approach, experimental results using real AF ECG recordings, comparing LCAGL with the Löwner-BTD computed using the NLS method, and with CAGL imposing Hankel constraints, are reported and discussed in this section. First, the indices used to measure AA quality and the observed database are summarized. Then, the experimental results from the two challenging scenarios of this arrhythmia, which were previously described in Section I, are shown and discussed.

A. Signal Quality Measurement in AF Episodes

Measuring estimation quality (or the AA content) of real signals is a difficult task. Since in practice there is no ground truth for comparison, one needs to take advantage of some features which are typically present in AA during AF episodes and use them as surrogate performance criteria. For example, in the frequency domain, the AA during AF has a peak between 3 and 9 Hz. The position of this peak is called dominant frequency (DF). In this section, two parameters used to measure AA estimation quality are presented. The first one, used to measure AA extraction quality, is the spectral concentration (SC), defined as the relative amount of energy around the DF. The SC is computed as in [START_REF] Castells | Spatiotemporal blind source separation approach to atrial activity estimation in atrial tachyarrhythmias[END_REF]:

SC = 1.17fp fi=0.82fp P S (f i ) Fs/2 fi=0 P S (f i ) (36)
where f p is the value of the DF, F s is the sampling frequency, f i is the discrete frequency and P S is the power spectrum of the source signal computed using Welch's method, as in [START_REF] Castells | Spatiotemporal blind source separation approach to atrial activity estimation in atrial tachyarrhythmias[END_REF].

The second parameter, that also provides a quality measurement of AA extraction, is the power contribution to lead V1, denoted P (r), given by [START_REF] De Oliveira | Source analysis and selection using block term decomposition in atrial fibrillation[END_REF]:

P (r) = 1 N ||m V 1,r s r. || 2 (37)
in mV 2 , where m V 1,r is the contribution of the r th source to lead V1 and s r. is the r th source in time domain. The P (r) of an AA source is expected to be relatively strong (> 10 -4 mV 2 ), since lead V1 is the one that typically best reflects AA in AF ECGs, as this lead strongly correlates with the AA from the right atrium and moderately correlates with that from the left atrium [START_REF] Raine | Principal component analysis of atrial fibrillation: Inclusion of posterior ECG leads does not improve correlation with left atrial activity[END_REF].

B. AF Database and Experimental Setup

Real ECG data from 20 patients suffering from persistent AF are used in this study. All the recordings belong to a database provided by the Cardiology Department of Princess Grace Hospital Center, Monaco. The recordings are acquired The 10 segments with disorganized and/or weak AA, from patients P1 to P10 are composed of one heartbeat, i.e., the QRS complex followed by the T wave and the visible f waves, and have between 0.8 and 1.7 seconds. The other 10 segments with short R-R intervals from patients P11 to P20 have 1.5 seconds of duration with at least 2 QRS complexes. All the 20 segments are downsampled by a factor of 10, since this significantly reduces the computational complexity without any noticeable loss of quality, as the signals' energy is concentrated in lower frequencies and no spectral aliasing is generated [START_REF] De | Alternating group lasso for block-term tensor decomposition with application to ECG source separation[END_REF].

The NLS method used for comparison in the experiments is available in Tensorlab MATLAB Toolbox and is set up as in [START_REF] De Oliveira | Löwner-based tensor decomposition for blind source separation in atrial fibrillation ECGs[END_REF]. As this method is often very sensitive to initialization, 10 Monte Carlo runs are performed for each segment and the best one is chosen. For CAGL and LCAGL, no Monte Carlo runs are needed and the value of γ between 8×10 -4 and 10 -2 In Figure 6, we can see the observed recording, the estimated VA and AA of the segment of Patient P2 processed by LCAGL. The segments are shown now on lead V1 for a better clarity of the estimated AA, in the time domain on the left and in the frequency domain on the right, with a DF equal to 6.67 Hz, typical of AA. It can be seen that the SC provided by LCAGL is 58.73% and P (r) value equal to 2.02 × 10 -4 mV 2 , while the ones provided using the Löwner-BTD computed by the NLS method, choosing the best match (R, L r ), are 5.36% and 1.59 × 10 -4 mV 2 , respectively, for the same segment. It should be noted that the Hankel-BTD computed by CAGL was not able to extract the AA from this recording, nor from any other used recording characterized by a disorganized and/or weak AA. Indeed, Hankel-BTD attempts to model the AA present in the ECG, which is a difficult task under these conditions. Due to the weak AA amplitude, its estimate is basically a straight line as shown in Figure 5 (Top).

Expanding the previous result to the population of 10 patients with segments presenting this particular characteristic, Figure 7 (Top) shows the SC values of the AA source estimated by the Löwner-BTD computed by the NLS method and LCAGL, for each patient. Regarding this AA quality index, LCAGL outperforms the NLS method in all but patient P4, where the SC of NLS is 1% higher than the one provided by LCAGL. In addition, the NLS method was not able to successfully extract the AA from patient P7, and CAGL could not successfully extract the AA from any of the 10 observed patients, as stated before. Also, the average SC over the 10 patients is 45.48% for LCAGL and 27.12% for the NLS method. Figure 8 shows how the P (r) values vary through the population of patients for both Löwner-based methods. No significant difference between the methods can be observed. However, a significant difference is observed between the two patient populations. The P (r) of the patients with a disorganized and/or weak AA signal is lower than the ones characterized by short R-R intervals and varies in a shorter range. These factors are indicative of fine AF. and LCAGL, for each patient whose segments are characterized by short R-R intervals. In 6 out of the 10 patients of this population, LCAGL provided an AA estimate with significantly higher SC than NLS, while for the other 4 patients the SCs of both methods are practically the same. More precisely, for those 4 patients LCAGL still produces a slightly higher SC value (around 1%). In addition, the mean SC over these 10 patients is 45.72% for LCAGL and 39.25% for the NLS method. Figure 9 shows the VA and AA estimates on lead V1 of the processed segment of patient 20. One can see that the several QRS complexes are well estimated, resulting in a clear AA signal when subtracted from the original ECG. The SC provided by LCAGL is 52.48%, practically equal to the 52.45% provided using the NLS method, manually choosing the best match (R, L r ) for the same segment. LCAGL provides a slightly stronger P (r), than the NLS method, which yields 3.11 × 10 -3 mV 2 for this segment.

D. Segments With Short R-R Intervals

Figure 10 shows the observed segment of Patient P18 on lead V1 and compares AA source estimations by the tensorbased methods. The AA estimate by CAGL (blue line) is deformed, probably due to the presence of 2 QRS complexes and the shape of the AA that may not be well approximated by an all-pole model, hampering the AA extraction. This loss of physiological shape also occurs in Figure 5 (bottom). Although the AA source estimated by NLS is better than the one produced by imposing the Hankel structure, it contains some ventricular residuals, due to an inaccurate VA source estimation. Finally, it can be seen that for this segment, LCAGL is the method that provides the best VA cancellation from the original recording, yielding the best AA estimate (green line).

We note that, for all but patients P8 and P17, both Löwnerbased methods (NLS and LCAGL) provide the same DF for the AA source estimate. Before concluding this experimental assessment, it is important to highlight that classical techiques that estimate the VA and subtract it from the ECG to retrieve the AA, such as average beat subtraction and adaptive singular value cancellation [START_REF] Alcaraz | Adaptive singular value cancellation of ventricular activity in single-lead atrial fibrillation electrocardiograms[END_REF] will certainly fail in theses short recordings, as they need sufficiently long ECGs to perform an accurate VA estimation.

VIII. DISCUSSION

Extracting the AA from AF ECGs is a particularly challenging task when the AA presents a very low amplitude and short R-R intervals appear often along the recording, as typically occurs in persistent forms of AF. In these cases, methods that focus on f-waves estimation do not provide a clean AA extraction or even fail. The proposed LCAGL method overcomes such limitations by taking into account that QRS complexes can be well modeled by rational functions. Estimating the VA can be simpler in such scenarios as the QRS morphology is not significantly changed during AF episodes as compared to NSR.

When compared to other tensor-based methods, LCAGL is more robust to the initialization of model parameters, as observed in the experiments on real and synthetic data. Also, the method guarantees the Löwner structure of BTD matrix factors. As compared to classical techniques for AA extraction based on QRS cancellation, LCAGL can operate successfully in very short ECG records, about one heartbeat long, whereas competing techniques require recordings typically composed of tens of seconds to estimate the QRS template well. This interesting advantage is a direct consequence of the deterministic signal model assumed by the BTD tensor approach, and renders LCAGL very attractive for real-time AF analysis. Another advantage lies in the fact that, after source separation, no additional processing is needed to automatically identify the AA source of interest, since the AA contribution results from the subtraction of the estimated VA from the original ECG. In consequence, potential errors resulting from AA source misidentification are avoided.

Limitations of the present work include the relative short sample size (20 patients) of the persistent AF population used in the experimental assessment, from which statistically solid conclusions are difficult to draw. The instability of LCAGL in estimating the T-wave should also be mentioned. Indeed, Fig. 6 shows an example where the T-wave is also captured with the QRS complex, whereas in Fig. 9 the T-waves are missing in the VA signal estimated by the same method. A possible explanation is that it is quite hard to distinguish the Twave from the AA signal in segments with short R-R intervals, as one can see in Fig. 9. It should be remarked, however, that T-waves pose difficulties to most AA extraction methods in the literature.

IX. CONCLUSION

Persistent AF is an advanced stage of this arrhythmia where the ECG typically presents a fast heart rate, resulting in short R-R intervals and an AA signal that is difficult to identify, making the AA extraction a challenging task. The present work has put forward and investigated a tensor-based approach to model the VA and separate it from the AA in persistent AF ECGs, showing a number of benefits in these challenging scenarios. This approach consists in explicitly modeling the VA with rational functions and mapping them into Löwner matrices in order to built a BTD tensor model. The factors of this third-order tensor model are computed using the LCAGL algorithm introduced in this work. The Löwner structure of the matrices that compose the tensor is guaranteed by a linear projection onto the Löwner subspace derived here.

After an evaluation on synthetic signals, the Löwner-BTD computed by LCAGL was applied to a population of 20 patients suffering from persistent AF, whose segments are characterized by short R-R intervals and a disorganized and/or weak AA signal contribution. Experimental results with real AF ECGs in very short observation windows (around 1-1.5 seconds) showed LCAGL's better AA extraction performance than two other tensor-based methods in this challenging clinical scenario, while offering the possibility of beat-to-beat online processing.

More robust tensor-based techniques should be studied in future works, in order to derive methods for the computation of an approximate BTD model that requires less parameter tuning and is able to better capture the T-wave when hard to discern from the f-waves. Also, experiments in a larger population of AF patients are envisaged in order to confirm the clinical relevance of these results.

Fig. 1 :

 1 Fig.1: An example of a persistent AF ECG recording, with short R-R intervals and fine AF.

Fig. 2 :

 2 Fig. 2: Graphical representation of the BTD of an arbitrary third-order tensor.

Fig. 3 :

 3 Fig.3: Synthetic AF ECG models simulating the challenging cases where the AA signal is weak (Model 1, top) and the R-R intervals are short (Model 2, bottom). Solid blue, dashed red and dash-dotted yellow lines represent VA, AA and noise, respectively.

Fig. 4 :

 4 Fig. 4: Synthetic AF ECG models. Original and estimated VA signals in blue dashed and gray solid lines, respectively, along with the estimation NMSE (dB) and the rank of the associated Löwner matrix factor, rank B .

Fig. 5 :

 5 Fig.5: Real AF ECG records. Top plots: an example of a single heartbeat segment of an AF ECG recording with a weak AA signal (solid green) and its respective AA signal estimate (dashed red) by the Hankel-BTD method CAGL. Bottom plots: an example of an AF ECG recording used in these experiments with 3 QRS complexes (solid blue) and its respective AA signal estimate (dashed orange), also by Hankel-BTD CAGL. The segments are shown on lead V1.

Fig. 6 :

 6 Fig. 6: Real AF ECG records. Left: Original recording, VA and AA estimates on lead V1 by LCAGL of Patient P2 in the time domain. The AA signal is vertically shifted for clarity. Right: Power spectrum of AA estimate.

Fig. 7 :Fig. 8 :

 78 Fig. 7: Real AF ECG records. SC values (%) of the AA source estimated by NLS and LCAGL for each patient of the observed database. Top: patients P1-P10 with disorganized and/or weak AA signals. Bottom: patients P11-P20 with short R-R intervals.

Figure 7 (Fig. 9 :

 79 Figure 7 (bottom) shows the SC values of the AA source estimated by the Löwner-BTDs computed by the NLS method

Fig. 10 :

 10 Fig. 10: Real AF ECG records. Estimated atrial sources contribution to lead V1 in a segment of short R-R intervals of Patient P18 by CAGL and the Löwner-BTDs computed by the NLS method and LCAGL. AA signal estimates are vertically shifted for clarity.

  with scalar entries a i1,i2 , has its i th 1 row and i th 2 column represented by a i1. and a .i2 , respectively. Symbol || • || denotes the l 2 -norm and || • || 2,1 denotes the mixed l 2,1 -norm, defined for an arbitrary matrix A with I 2 columns as:

TABLE I :

 I Löwner-Constrained Alternating Group Lasso (LCAGL) algorithm.

TABLE II :

 II Parameters of the synthetic AA signal models.

	Model P	x	∆x	fx	fs	f 0	∆f	f f
	1	5	150	50	0.08 1000	6	0.2	0.10
	2	3	60	18	0.50 1000	8	0.3	0.23

ACKNOWLEDGMENT P. M. R. de Oliveira was funded by a PhD scholarship from the IT Doctoral School (ED STIC) of Université Côte d'Azur. J. H. de M. Goulart's work was partially supported by the MIAI chair "LargeData" of Université Grenoble Alpes.

V. Zarzoso holds the 3IA chair "IAblation" of Université Côte d'Azur.

This work has been supported in part by the French government, through the 3IA Côte d'Azur Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-19-P3IA-0002.