Kurisummoottil Christo

Dirk Thomas

Slock
email: slock@eurecom.fr

BP-VB-EP BASED STATIC AND DYNAMIC SPARSE BAYESIAN LEARNING WITH KRONECKER STRUCTURED DICTIONARIES

Keywords: Sparse Bayesian Learning, Variational Bayes, Tensor Decomposition, Kronecker Structured Dictionary Learning, Belief Propagation

In many applications such as massive multi-input multi-output (MIMO) radar, massive MIMO channel estimation, speech processing, image and video processing, the received signals are tensors. In such applications, utilizing techniques from tensor algebra can be beneficial since it retains the tensorial structure in the received signal compared to processing on the matricized version of the same signal. Furthermore, the underlying parameters or states to be estimated are sparse in many of the above-said applications compared to the large system dimensions. In this paper, we propose techniques which allow handling the extension of sparse Bayesian learning (SBL) to time-varying states. Adding the parameters of the autoregressive process which is used to the model the time-varyings of the state leads to a non-linear (at least bilinear) state-space model. Belief propagation (BP) is a promising method to compute the minimum mean squared error (MMSE) or maximum a posteriori (MAP) estimates, but at the the expense of a high computational burden. However, inspired by a previous work on a combined BP and variational Bayes (VB) technique, we noted that using a combination of BP, VB, and expectation propagation (EP) can help to alleviate the computational complexity.

INTRODUCTION

The signal model for the recovery of a time varying sparse signal under Kronecker structured (KS) [START_REF] Duarte | Kronecker compressive sensing[END_REF][START_REF] Shakeri | Identifiability of Kronecker-structured dictionaries for tensor data[END_REF] dictionary matrix can be formulated as Observation: yt = (A

(t) 1 ⊗ A (t) 2 .... ⊗ A (t) N )xt + vt, State Update: xt = Fxt-1 + wt, (1)
where yt = vec(Yt), ⊗ represents the Kronecker product between two matrices, vec(•) representing the vectorized version of the tensor or matrix (•), Yt ∈ C I 1 ×I 2 ...×I N is the observations or data at time t, A t j,i ∈ C I j , the factor matrix A (t) T j,i ] T , justification for which can be seen in our previous work [START_REF] Thomas | Space alternating variational estimation and kronecker structured dictionary learning[END_REF]. We also define the measurement of dictionary matrix as,

A (t) = A (t) 1 ⊗ A (t) 2 .... ⊗ A (t) N = N j=1 A (t)
j . Sparse signal xt is modeled using an AR [START_REF] Duarte | Kronecker compressive sensing[END_REF] process with a diagonal correlation coefficient matrix F. xt contains only K non-zero entries (or significant EURECOM's research is partially supported by its industrial members: ORANGE, BMW, Symantec, SAP, Monaco Telecom, iABG, and by the projects DUPLEX (French ANR), MASS-START (French FUI) and EU ITN project SPOTLIGHT. coefficients), with K << M and thus the dictionary matrix to be learned allows a low rank representation. vt is assumed to be a white Gaussian noise, vt ∼ N (0, γ -1 I). The Sparse Bayesian Learning (SBL) algorithm was first introduced by [START_REF] Tipping | Sparse Bayesian learning and the relevance vector machine[END_REF] in the machine learning context and then proposed for the first time for sparse signal recovery by [START_REF] Wipf | Sparse Bayesian Learning for Basis Selection[END_REF]. The efficiency of the SBL approach revolves around the hierarchical prior modeling and the Bayesian LASSO [START_REF] Park | The Bayesian Lasso[END_REF] is a special case of SBL with Gaussian-Exponential hierarchical prior (equivalent to Laplacian marginal for xt). Dynamic autoregressive SBL (DAR-SBL) considered here is a case of joint Kalman filtering (KF) with a linear time-invariant diagonal state-space model, and parameter estimation, which can be considered an instance of nonlinear filtering.

Prior work on dictionary learning (DL) [START_REF] Lewicki | Learning overcomplete respresentations[END_REF][START_REF] Skretting | Recursive least squares dictionary learning algorithm[END_REF][START_REF] Aharon | K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF][START_REF] Roemer | Tensor-based algorithms for learning multidimensional separable dictionaries[END_REF][START_REF] Ding | Joint sensing matrix and sparsifying dictionary optimization for tensor compressive sensing[END_REF] focus on either maximum likelihood based schemes, LS or K-SVD algorithms. Due to space limitations, we refer the readers to a more detailed discussion on the state of the art in DL to our previous work [START_REF] Thomas | Space alternating variational estimation and kronecker structured dictionary learning[END_REF].

Contributions of this paper

• This paper is an extension of our previous work based on KS dictionaries [START_REF] Thomas | Space alternating variational estimation and kronecker structured dictionary learning[END_REF] to the case of time varying sparse state vector. Moreover, we observed that variational Bayes (VB) at the scalar level (for xt) and at the column level for the factor matrices A (t) j is quite suboptimal [START_REF] Thomas | Gaussian variational Bayes Kalman filtering for dynamic sparse Bayesian learning[END_REF] since the posterior covariance computed does not converge to the true posterior. Hence, this requires more accurate approximations like belief propagation (BP) which we propose here.

• Building on the framework of [START_REF] Riegler | Merging belief propagation and the mean field approximation: a free energy approach[END_REF], we combine BP and MF approximations in such a way as to optimize the message passing (MP) framework. The main focus being on getting reduced complexity (VB accounts for the low complexity) algorithms without sacrificing much on the performance (BP for performance improvement).

2. SBL DATA MODEL 1 Starting from the system model (1), diagonal matrices F (assumed to be deterministic) and Γ are defined with its elements, Fi,i = fi, fi ∈ (-1, 1) and

Γ -1 = diag (α), α = [α1, ...αM ].
Here αi represents the inverse variance of xi,t ∼ N (0,

1 α i ). Further, wt ∼ N (0, Λ -1 ), where Λ -1 =Γ(I -FF T )= diag ( 1 λ 1 , ..., 1 λ
M ). wt are the complex Gaussian mutually uncorrelated state innovation sequences. Hence we sparsify the prediction error variance wt also, with the same support as x0 and henceforth enforces the same support set for xt, ∀t. vt is independent of the wt process. The above signal model finds numerous applications including but not limited to 1) Bayesian adaptive filtering [START_REF] Sadiki | Bayesian adaptive filtering: principles and practical approaches[END_REF][START_REF] Ciochina | A family of optimized LMS-based algorithms for system identification[END_REF], 1 In this paper, boldface lower-case and upper-case characters denote vectors and matrices respectively. The operators tr(•), (•) T , (•) * , (•) H , • represent trace, transpose, conjugate, conjugate transpose and Frobenius norm, respectively. A complex Gaussian random vector with mean µ and covariance matrix Θ is denoted as x ∼ N (µ, Θ). diag (•) represents the diagonal matrix created by elements of a row or column vector. The operator < x > or E(•) represents the expectation of x. I N represents the identity matrix with dimension N . δ(•) represents the Kronecker delta function.

2) Wireless channel estimation: multipath parameter estimation as in [START_REF] Fleury | Channel Parameter Estimation in Mobile Radio Environments Using the SAGE Algorithm[END_REF][START_REF] Bajwa | Compressed channel sensing: A new approach to estimating sparse multipath channels[END_REF][START_REF] Cheng | Algebraic Channel Estimation Algorithms for FDD Massive MIMO systems[END_REF]. In Bayesian compressive sensing, a two-layer hierarchical prior is assumed for the xt as in [START_REF] Tipping | Sparse Bayesian learning and the relevance vector machine[END_REF]. The hierarchical prior is chosen such that it encourages the sparsity property of xt or of the innovation sequences vt. The state update gets represented as, p(xt/xt-1,

F, Γ)= M i=1 N (fixi,t-1, 1 α i ).
For the convenience of analysis, we reparameterize αi in terms of λi. Since λi ∝ αi, 1/λi's are also sparse and hence we can assume a Gamma prior for Λ, p(Λ) = e -bλ i , such that the marginal pdf of xt (student-t distribution) becomes more sparsity inducing than e.g., a Laplacian prior. The inverse of noise variance γ is distributed as, Gamma(γ; c, d). The advantage is that the whole machinery of linear MMSE estimation can be exploited, such as e.g., the KF. But this is embedded in other layers making things eventually non-Gaussian. Now the likelihood distribution can be written as, p(yt/xt, γ) ∝ γ N e -γ y t -A (t) x t 2 . To make these priors non-informative (Jeffrey's prior), we choose them to be small values a = c = b = d = 10 -5 . We define the unknown parameter vector θt = {xt, Λ, γ, F, A (t) }.

We also remark that the techniques proposed here (with KS DL) are an instance of gridless compressed sensing, for e.g. [START_REF] Ardah | A Gridless CS Approach for Channel Estimation in Hybrid Massive MIMO Systems[END_REF]. To further elucidate the application of the system model in (1), we consider the wireless channel estimation in the case of millimeter wave or massive MIMO systems. The time varying channel impulse response H (t) has per path a rank one contribution in four dimensions (Tx and Rx spatial multi-antenna dimensions, delay spread and Doppler spread). In the frequency domain

vec(H (t) ) = Np i=1 xi,t ht(ψi) ⊗ hr(φi) ⊗ v f (τi(t)) ⊗ vt(fi) = A(θ) (t) xt
where v f (.), vt(.) are appropriate Vandermonde vectors (possibly subsampled in the case of v f (.)). Hence we get a sum of rank one 4D tensors. We emphasize that the presented algorithm do not exploit parametric forms, because those parametric forms are uncertain. For eg., considering the massive MIMO channel estimation problem [START_REF] Thomas | Variational Bayesian learning for channel estimation and transceiver determination[END_REF], the array response at the mobile station (MS) is not exploitable in practice. Even the array response at the base station (BS) will typically require calibration to be exploitable or may have mutual coupling effects which are unknown. Doppler shifts can be clearly represented as Vandermonde vectors. Delays could be more or less clear, if one goes to frequency domain in OFDM, and one only takes into account the range of subcarriers for which the Tx/Rx filters can be considered frequency-flat. Then over those subcarriers, it's also Vandermonde. We do limit our discussions in this paper to the case of known A (t) . However, note that here we do consider the extension of the SAVE and BP based algorithm discussed herein to the case of unknown A (t) under the dynamic case (static SBL is a special case).

COMBINED BP/MF APPROXIMATION USING VARIATIONAL FREE ENERGY OPTIMIZATION

The time index t is omitted here for simplicity. In this section, we develop an MP framework combining BP and VB which are derived from variational free energy (VFE) which is a fundamental quantity in statistical physics. We refer the readers to [START_REF]Low Complexity Static and Dynamic Sparse Bayesian Learning Combining BP, VB and EP Message Passing[END_REF] for a detailed technical discussion on VFE and the theoretical formulation behind the derivation of MP expressions. We instead brief the results here. Assume that the posterior be factorized as, p(θ) = Fig. 1. Factor Graph (FG) for the dynamic SBL (at time t). Note that messages from the smoothing stage is not shown here.

N (i), N (a)
represent the number of neighbouring nodes of any variable node i or factor node a, respectively. Let ma→i represents the message passed from any factor node a to variable node i and ni→a represents the message passed from any variable node i to factor node a. The fixed point equations corresponding to the constrained optimization VFE can be written as follows [START_REF] Riegler | Merging belief propagation and the mean field approximation: a free energy approach[END_REF] (<>q represents the expectation w.r.t distribution q) qi(θi) = zi

a∈N BP (i) m BP a→i (θi) a∈N M F (i) m M F a→i (θi), ni→a(θi) = b∈N BP (i)\a m b→i (θi) b∈N M F (i) m b→i (θi), m M F a→i (θi) = exp(< ln fa(θa) > j∈N (a)\i n j→a (θ j ) ), m BP a→i (θi) = ( j∈N (a)\i nj→a(θj)fa(θa) j =i dθj), (2)

Dynamic BP-MF-EP based SBL

The figure 1 represents the FG (note that static case is a special case with the state update nodes being not present), where it is divided into two disjoint subsets ABP = f δ n,t ∀n, l, t and AMF represents rest of the factor or variable nodes. To combine BP and MF, we introduce the new variables hn,t = A (t) n,:xt, s l,t = f l x l,t-1 and the hard constraint factor nodes, f δ n,t = δ(hn,t -A (t) n,:xt), ∀n ∈ [1 : N ], t, and f∆ l,t = δ(s l,t -f l x l,t-1 ), ∀l ∈ [1 : M ], t. We can compute m f δ n,t →x l,t (x l,t ) = f δ n,t n h n,t →f δ n,t (hn,t)

l =l n x l ,t →f δ n,t (x l ,t ) l =l
dx l ,t . For notational brevity, we denote subscript (l, n) or (n, l) to represent the messages passed from l to n or viceversa. All the messages (beliefs or continuous pdfs) passed between them can be shown to be Gaussian [START_REF] Tan | Computationally efficient sparse Bayesian learning via belief propagation[END_REF] and thus it suffices to represent them by the mean and variance of the beliefs. The joint distribution of all the observations and parameters can be written as, p(yt, θt/y1:t-1) = p(yt/θt)p(θt/y1:t-1), where p(θt/y1:t-1) denotes the predictive distribution. Similar as in KF, first we compute the posterior distribution of θi,t given the observations till (t -1), which is called as the prediction stage.

Diagonal AR(1) ( DAR(1) ) Prediction Stage

Since there is no coupling between the scalars in the state update (1), it is enough to update the prediction stage using MF. However, the interation between x l,t and f l requires Gaussian projection, using expectation propagation (EP). For more detailed derivation, we refer to our previous work [START_REF] Thomas | Gaussian variational Bayes Kalman filtering for dynamic sparse Bayesian learning[END_REF] due to space limitations.

Measurement Update (Filtering) Stage

For the measurement update stage, the posterior for xt is inferred using BP. Note that we represent the mean of the messages by

x (t) n,l , ν (t)
n,l . The mean and variance of the beliefs computed at x l,t are denoted by x l,t|t , σ 2 l,t|t . In the measurement stage, the prior for x l,t gets replaced by the belief from the prediction stage. We refer to our previous work [START_REF]Low Complexity Static and Dynamic Sparse Bayesian Learning Combining BP, VB and EP Message Passing[END_REF] for detailed derivations and expressions for the messages. We define d l,t = (

N n=1 ν (t) -1 n,l ) -1 , r l,t = d l,t ( N n=1 x (t) n,l ν (t) n,l + x l,t|t-1 σ 2 l,t|t-1
). Given the messages, m f δ n,t →x l,t (x l,t ), the belief q(x l,t ) can be obtained as (

f λ i (λi) = p(λ k /a, b)), q(x l,t ) ∝ f λ i (λi) N n=1 m f δ n,t →x l,t ∝ N (x l,t ; x l,t|t , σ 2 l,t|t ), where σ -2 l,t|t = λ l,t + d -1 l,t , x l,t|t = r l,t 1+d l,t σ -2 l,t|t . (3)
One remark here is that compared to our previous work using VB [START_REF] Thomas | SAVE -Space alternating variational estimation for sparse Bayesian learning[END_REF], combining BP and MF gives a more accurate approximation of the error variance as shown in [START_REF] Thomas | Space alternating variational estimation and kronecker structured dictionary learning[END_REF], where σ 2 l,t|t incorporates the effect of all σ 2 l ,t|t , l = l.

Lag-1 Smoothing Stage

We show in [21, Lemma 1] that KF is not enough to adapt the hyperparameters, instead we need at least a lag 1 smoothing (i.e. the computation of x k,t-1|t , σ 2 k,t-1|t through BP). For the smoothing stage, we use BP with Gaussian Markov Random Fields (GMRF) based factorization. GMRF refers to the representation of BP [START_REF] Yedidia | Constructing Free Energy Approximations and Generalized Belief Propagation Algorithms[END_REF], when the underlying Gaussian distribution is expressed in terms of pairwise connections between scalar variables xi,t. We skip the detailed derivation and instead refer to our paper [START_REF]Low Complexity Static and Dynamic Sparse Bayesian Learning Combining BP, VB and EP Message Passing[END_REF]. Applying the MF rule from (2), the resulting Gaussian distribution has mean, σ -2 f i |t and variance, f i|t , the detailed derivations for which are in [21, Section 3.2.3]. The entire algorithm (a combination of BP, MF and EP, we call it as Combined BP-MF-EP DAR-SBL) is described in Algorithm 1. Also we remark that for the estimation of λ l , γ, we follow the same approach as in our paper [START_REF] Thomas | Gaussian variational Bayes Kalman filtering for dynamic sparse Bayesian learning[END_REF] and we refer to it for more details. One remark here is that another version called as Combined Vector BP-MF-EP DAR-SBL follows immediately from the derivations for Algorithm 1, where all the components of xt are considered jointly in the FG. Even though the performance will be higher (as observed in the simulations) for the vector case, it comes at the cost of a higher complexity due to the matrix inversion involved. Note that in Algorithm 1, we introduce temporal averaging for certain quantities (represented by <> |t ) in hyperparameter estimates and β being the temporal weighting coefficient which is less than one, see [START_REF] Thomas | Gaussian variational Bayes Kalman filtering for dynamic sparse Bayesian learning[END_REF] for more details. For the KS DL, the algorithm remains same as in our previous work [START_REF] Thomas | Space alternating variational estimation and kronecker structured dictionary learning[END_REF], which is denoted as space alternating variational estimation with Kronecker structured DL (SAVED-KS DL).

Suboptimality of SAVED-KS DL and Joint VB

First, we define the unfolding operation on an N th order tensor

Yt = [[A (t) 1 , ..., A (t) N ; x]] as [25] (Y (n) t is of size In × N i=1,i =n Ii) Y (n) t = A (t) n X (n) t (A (t) N ⊗ A (t) N -1 ...A (t) n+1 ⊗ A t n-1 ... ⊗ A (t) 1 ) T . ( 4
) From the expression for the error covariance in the estimation of the factor aji (tr{(

1 k=N,k =j < A T k A * k >) < X (j) T X (j) >}I), it is
clear that it does not take into account the estimation error in the other columns of Aj. The columns of Aj can be correlated, for e.g. if we consider two paths (say i, j) with same DoA but with different delays, the delay responses v f (τi(t)) and v f (τi(t)) may be corre-

Algorithm 1 Combined BP-MF-EP DAR-SBL with KS DL

Initialization: f l|0 , λ l|0 = a b , γ0 = c d , x l,0|0 = 0, σ 2 l,0|0 = 0, ∀l. Define Σ t-1|t-1 = diag (σ 2 l,t|t-1 ). for t = 1 : T do Prediction Stage: 1. From [12], x l,t|t-1 = f l|t-1 x l,t-1|t-1 , σ 2 l,t|t-1 = | f l|t-1 | 2 σ 2 l,t-1|t-1 + σ 2 f l |t-1 (| x l,t-1|t-1 | 2 + σ 2 l,t-1|t-1 ) + λ -1 l|t-1 . Filtering Stage: 1. Compute x (t) n,l , ν (t)
n,l from [21, eq. ( 5)] and update x l,t|t , σ -2 l,t|t from (3).

Compute ν

(t) l,n , x (t)
l,n from [21, eq. ( 7)]. 3. Continue steps 1) to 2) until convergence.

Smoothing Stage: Initialization

: Σ (0) t-1|t = Σ t-1|t-1 , x (0) t-1|t = x t-1|t-1 . Define B (t) =< F T A (t) T R -1 t A (t) F > +Σ t-1|t-1 , ht =< F T A (t) T > R -1 t yt. 1. Pi,j = -B (t) 2 i,j B (t) i,i + k∈N (i)\j P k,i
, µi,j = (hi,t +

k∈N (i)\j P k,i µ k,i ), ∀i, j. 2. σ -2 i,t-1|t = B (t) i,i + k∈N (i) P k,i , x i,t-1|t = σ 2 i,t-1|t (hi,t + k∈N (i) P k,i µ k,i )
Estimation of hyperparameters (Define:

x k,t = x k,t -f k x k,t-1 , ζt = βζt-1+ (1 -β) < yt -A (t) xt 2 >), b t = (< x k,t 2 > |t +b). 1. Compute f l|t , σ 2 f l|t from [21, eq. (11)], γt = c+N (ζ t +d) and λ l|t = (a+1) b t . SAVED-KS DL: aji = (bj ) 1 , bj = (Y (j) < X (j) >< ( 1 k=N,k =j A k ) T >)i, Υj,i = βj,iI, βj,i = tr{( 1 k=N,k =j < A T k A * k >) < X (j) T X (j) >}.
lated. However, since it is not clear how to model this dependency, we indeed keep it as a future work. This suboptimality in the error covariance estimate using SAVED-KS resulting from the correlation between the columns, can be avoided by using a joint VB [START_REF] Thomas | Space alternating variational estimation and kronecker structured dictionary learning[END_REF]. The joint VB estimates (mean and covariance) can be obtained as

M T j = A T 1,j =< γ > Ψ -1 j B T j , Ψj = (< γ >< X (j) ( 1 k=N,k =j < A T k A * k >)X (j) T >), (5)
where Vj =< X (j) >< (

1 k=N,k =j A k ) T >
and Bj is defined as with the first row of (Y (j) V T j ) removed. However, the joint VB involves a matrix inversion and is not recommended for large system dimensions. Nevertheless, it is possible to estimate each columns of Aj by BP, since each column estimate can be expressed as the solution of a linear system of equation from ( 5), a T j,i = Ψ -1 j bj,i. bj,i represents the i th column of B T j . The message passing expressions under BP (using a GMRF based FG) for the factor matrices can be written as

ζm,n = -(Ψj) 2 m,n /(ζm,m + k∈N (m)\j ζ k,m ), κm,n = (ζm,mκm,m)+ k∈N (m)\j ζ k,m κ k,m (Ψ j )m,n , (6)
where we initialize ζm,m = (Ψj)m,m, κm,m = (b j,i )m (Ψ j )m,m,ζm,n = 0, κm,n = 0. Finally the mean (κm) and variance (ζm) of the posterior belief can be computed as

κm = ζm,mκm,m+ k∈N (m) ζ k,m κ k,m ζm,m+ k∈N (m) ζ k,m , ζm = ζm,m + k∈N (m) ζ k,m . (7)
We remark that, the above BP based low complexity scheme for KS DL represents a major innovation compared to our previous work [START_REF] Thomas | Space alternating variational estimation and kronecker structured dictionary learning[END_REF], apart from the extension to the dynamic SBL case.

OPTIMAL PARTITIONING OF THE MEASUREMENT STAGE AND KS DL

In [START_REF] Thomas | Space alternating variational estimation and kronecker structured dictionary learning[END_REF], we derived the Fisher Information Matrix (FIM) for the KS DL where the sparse vector is static. Here, we reuse the FIM expressions to derive the optimal partitioning of the variables in the measurement stage. We refer to our paper [26, Lemma 1], where the main message was that if the parameter partitioning in VB is such that the different parameter blocks are decoupled at the level of FIM, then VB is not suboptimal in terms of (mismatched) Cramer-Rao Bound (mCRB). More detailed overview on mCRB can be found in [START_REF] Fortunati | Performance Bounds for Parameter Estimation under Misspecified Models [Fundamental findings and applications[END_REF]. If a finer partitioning granularity is used (such as up to scalar level as in MF), then VB becomes quite suboptimal, which can be alleviated by using BP instead.

Lemma 1. For the measurement stage, an optimal partitioning is to apply BP for the sparse vector xt and VB (SAVED-KS) for the columns of the factor matrices A (t) j,i assuming the vectors A (t) j,i are independent and have zero mean. However, if the columns of A (t) j are correlated, then a joint VB, with the posteriors of the factor matrices assumed independent, should be done for an optimal performance.

Proof: Let F (i) j = [0 I j ×I j (i-1) II j 0 I j ×I j (P j -i) ], Φj,t = vec(A (t) j )
. We observe that we can separate the contributions of J(Φt, xt) = [J(Φt) J(xt)], J(Φt) = [J(Φ1,t) ..... J(Φ N,t )] where, J(Φ j,t ) = F(xt)(Φ1,t ⊗ ...I I j P j .... ⊗ Φ N,t ),

A (t
J(xt) = [ F1( N j=1 Φ j,t ), ...., F M ( N j=1 Φ j,t )) ].
Further, the FIM for the case of SBL can be derived as [START_REF] Thomas | Space alternating variational estimation and kronecker structured dictionary learning[END_REF]

F IM =      E(γ)J(Φ t ) T J(Φ t ) 0 0 0 0 E(γ)J(x t ) T J(x t ) + E(Γ -1 ) 0 0 0 0 a E(Γ -2 ) 0 0 0 0 (N + c -1) E(γ -2 )      (8)
Here, γJ(xt) T J(Φt) = 0, since xt is zero mean. If the all columns of A (t) j are independent and zero mean, then E(J(Φt) T J(Φt)) becomes a diagonal matrix with no coupling between the free variables of any two different columns of the factor matrices. However, if any factor matrix is A (t) j is correlated, it is suboptimal to factorize the columns of A (t) j independently in the approximate posterior. Hence, in this case, a joint VB method (which has higher complexity) would be optimal to estimate the posterior distributions and this indeed justify the superior performance of joint VB approach described in Section 3.2.

SIMULATION RESULTS

For the observation model, the parameters chosen are N = 256, M = 200. For the simulations, we consider a 3 -D tensor with dimensions [START_REF] Tipping | Sparse Bayesian learning and the relevance vector machine[END_REF][START_REF] Skretting | Recursive least squares dictionary learning algorithm[END_REF][START_REF] Skretting | Recursive least squares dictionary learning algorithm[END_REF] and the number of non-zero elements of xt or the rank of the tensor (no of non-zero elements of xt) is fixed to be K = 16. All signals are considered to be real in the simulation. All the elements of the factor matrix A (t) j (time varying) are generated i.i.d. from a Gaussian distribution with mean 0 and variance 1. The rows of A (t) are scaled by √ 16 so that the signal part of any scalar observation has unit variance. Taking the SNR to be 20dB, the variance of each element of vt (Gaussian with mean 0) is computed as 0.01.

Consider the state update, xt = Fxt-1 + wt. To generate x0, the first 16 elements are chosen as Gaussian (mean 0 and variance 1) and then the remaining elements of the vector x0 are put to zero. Then the elements of x0 are randomly permuted to distribute the 30 non-zero elements across the whole vector. The diagonal elements of F are chosen uniformly in [0.9, 1). Then the covariance of wt can be computed as Γ(I -FF T ). Note that Γ contains the variances of the elements of xt (including t = 0), where for the non-zero elements of x0 the variance is 1. Following observations can be made from the simulations. In Figure 2, which is for static SBL case with DL, there is substantial improvement in NMSE compared to our previous work [START_REF] Thomas | Space alternating variational estimation and kronecker structured dictionary learning[END_REF]. Our proposed low complexity algorithm using BP has similar performance as that of joint VB which has higher complexity. In Figure 3, we evaluate the performance of the proposed BP-MF-EP DAR SBL and show that the parameter estimation benefits from BP. "MF DAR-SBL" refers to the sub-optimal version with no BP and only MF for filtering or smoothing of xt. Also we show the drastic improvement in performance with lag-1 smoothing for hyperparameter estimation compared to just using filtering.

CONCLUSION

N

  j ] which is unknown and the tensor product is repre-sented by [[A ; xt]], xt is the M (= N j=1 Pj)-dimensionalsparse signal and wt is the additive noise. We consider A

  ; a, b), where, Gamma(λi; a, b) = M i=1 Γ -1 (a)b a λ a-1 i

  f b (θ b ), where ABP , AMF represent the set of nodes belonging to the BP part and MF part, respectively. NBP (i) represents the number of neighbouring nodes of i which belong to the BP part, similarly NMF (i) is defined.

  ) and xt in (1) as, yt = (

Fig. 2 .

 2 Fig. 2. Static SBL: NMSE as a function of N .

Fig. 3 .

 3 Fig. 3. DAR-SBL: NMSE as a function of time.

We have presented here a low complexity algorithm for KS DL using a combination of BP, VB and EP. The motivation behind the proposed algorithm is to circumvent the suboptimality associated with incorrect posterior covariance computation for the columns of the factor matrices in a previous paper of ours which was entirely based on MF approximations. However, we are still unclear whether the proposed algorithm is robust enough to perform comparably with the variations in the model of KS A (t) , which is left as a future work.


Leonid Ryvkin

Focusing on gap lengths in the noisy violinist chip-firing problem

come

This description does not contain the information of the absolute position of the leftmost violinist, (but we don't care about this information anyways). I.e. a state for us is what is referred to as the shadow of a state in [START_REF] Khovanova | Ending States of a Special Variant of the Chip-Firing Algorithm[END_REF].

A (flat) N -clusteron is a state of N -1 consecutive zeroes, i.e. a • = 1 N -1 (where we use powers to express the concatenation of identical digits). This corresponds to N violinists in consecutive rooms.

We will now describe a move from one state of the system a t

• to the next one a t+1

• .

• Pick T with a T = 0 call T the triggering position.

• Set a T = 2,

• find the biggest j < T with a j > 0 and reduce it by one. If j does not exist skip this step

• find the smallest j > T with a j > 0 and reduce it by one. If j does not exist skip this step If no moves are possible (i.e. the sequence contains no more zeroes) we have reached a final state. We will denote a sequence of moves by a 1 • , ..., a n • and the corresponding triggering positions by T 1 , ..., T n-1 .

I.e., a move is a 0-gap triggering by turning into a 2 and two minus ones propagating left resp. right until they meet a non-zero number (which they would reduce by one) or disappear at the boundary of the state.

Here are all possible sequences of moves for N = 4 violinists, as described in Figure 1 of [START_REF] Khovanova | Ending States of a Special Variant of the Chip-Firing Algorithm[END_REF], in our notation:

Reaching a final state

To show that we always reach a final state, we will use the fact that there is a non-decreasing quantity:

Definition 1. We call |a • | = N -1
i=1 a i the norm of the state a • . We have the following property for norms of consecutive states a t

• , a t+1

• :

|a t+1 • | -|a t • | ∈ {0, 1, 2}, where the value 2 is only possible if |a t • | = 0 N -1
. This means that the norm is non-decreasing and can only stay the same or increase by one (unless the initial states is a flat N -clusteron, then it will increase by 2).

Theorem 2. For any initial state a 1

• ∈ N N -1 , the final state is reached after finitely many steps. Proof. Let k = max 1≤j≤N -1 a j and K = max(k, 2). Since the components of a either get smaller or get increased from 0 to 2, the number of possible states starting from a 1

• is finite (a subset of {0, 1, ..., K} n-1 ). Hence, either the assertion of the theorem is true or there is a sequence of n > 1 states a 1

• , ..., a n

• where the first one equals the last one. Let i be the leftmost changing value, i.e. the smallest value such that a t j = a t+1 j for all j < i and all t. Let now t be the step in which i changes, then

• either there are only zeros left to a t i , then |a t+1 • | > |a t • |, which implies that the sequence can not be a cycle, because the norm is non-decreasing.

• or there is a non-zero number left of i, which then changes (reduces by 1) in step t , which is a contradiction to our choice of i.

Hence there can be no cycle, i.e. we eventually run into a final state.

3 Characterization of possible final states [START_REF] Khovanova | Ending States of a Special Variant of the Chip-Firing Algorithm[END_REF] also gives a characterization of all possible final states when starting with a flat N -clusteron. Those are given by states containing just ones and one two (in an arbitrary position), i.e. they are of the form 1 k 21 N -2-k for k ∈ 0, ..., N -2. We will reprove this statement in this section. First we observe that if we start with a state in the space {0, 1, 2} N -1 , then we will always stay in this space (gaps of length > 2 will never appear). On the other hand, the norms of the states are bounded from above (even locally):

Theorem 3. Starting with a 1 • = 0 N -1 , there can never appear a state a t • containing a sub-sequence of length l with norm higher than l + 1. I.e. there is no l and i such that l-1 i=0 a t j+i > l + 1. Proof. We do induction by length l. For l = 1 this is just the statement that no numbers higher than 2 can appear from a state in {0, 1, 2} N -1 . Assume we have proven the statement for lengths up to l. Further assume t is the first step in which l j=0 a t i+j > l + 2 We can further assume that the value of l j=0 a t-1 i+j is l + 2. It is at most l + 2 by assumption, and the only case where it would increase by more than one during step t -1 is if the whole sub-sequence (a t-1 i , ..., a t-1 i+l ) was zero, in which case l j=0 a t i+j = 2 ≤ l -2. For the sum l j=0 a t-1 i+j to increase, either a t-1 i or a t-1 i+l has to be zero in which case l-1 j=0 a t-1 i+1+j or l-1 j=0 a t-1 i+j had value l + 2 , contradicting the induction hypothesis.

Figure 1 :

 1 Figure 1: Possible move sequences for 4 violinists.

Final states never contain zeroes and any (non-initial) state contains at least one 2. Since by the above the norm of the state can never surpass N -1 + 1 = N , we obtain: Corollary 4. Starting with a 1

• = 0 N -1 , the final position always consists of one two and the rest ones, i.e. is of the form 1 k 21 N -2-k for k ∈ 0, ..., N -2.

We still need to show that all positions of the two are possible, i.e. that 1 k 21 N -2-k is a possible final state for all k ∈ 0, ..., N -2. We will prove this by induction and need the following statement: Lemma 5. Assume there is a sequence of moves to obtain a state (a 0 , ..., a N -2 ) from 0 N -2 . Then there is a sequence of moves to obtain (a 0 , ..., a i-1 , 0, a i , ...., a N -2 ) from 0 N -1 for any i ∈ 0, ..., N -1.

Proof. If T 1 , ..., T t is the sequence of trigger points for (a 0 , ..., a N -2 ), then we set

and T1 , ..., Tt become a sequence of trigger points leading to (a 0 , ..., a i-1 , 0, a i , ...., a N -2 ) from 0 N -1 .

We can now prove that Theorem 6. Starting with a 1 • = 0 N -1 all states of the form 1 k 21 N -2-k for k ∈ 0, ..., N -2 can be final.

Proof. Again we proceed by induction. For N -1 = 1 (i.e. a single gap) there is only one final state (2) possible. Assume we have proven the statement for N -2 gaps. We now want to construct a sequence of moves turning 0 N -1 into 1 k 21 N -2-k .

• For k = 0, we use the above Lemma and the induction hypothesis to arrive at the state 021 N -3 . Then we trigger the leftmost spot.

• For k = N -2, we mirror the above construction, i.e. we generate 1 N -3 20 and trigger the rightmost position.

• For the other k we use the induction hypothesis and Lemma to arrive at the state 01 k-1 21 N -2-k . By triggering the only possible zero repeatedly we will move the zero further and further to the right until we arrive at the desired state.


Roberto Boccagna

Maurizio Bottini

Massimo Petracca

Mirko Rossi

Alessia Amelio

Guido Camata

Inte

Alan O'connor

Integrating Data-driven and Model-based Algorithms for Structural Health Monitoring

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

In the last few decades, Structural Health Monitoring (SHM) has benefited greatly from the outstanding progress made in the field of outlier detection-oriented algorithmic theory and from the considerable increase in CPU velocity and GPU performance in parallel computing. Such progress paved the way for anomaly detection in the context of Civil Engineering, making available the possibility of fast data collection, transmission and pro-cessing methods for a near real-time response of the overall structural state. Over the years, new, automated techniques stood alongside traditional, i.e., in-person, facility surveillance, and gradually supplanted visual inspection and monitoring as soon as networks of synchronized sensors were integrated in "smart structures". This paradigm change resulted in a substantial cost reduction for maintenance and, above all, in increased reliability when assessing structural conditions. Data-driven methods, i.e., detection techniques entirely based on the analysis of real data, enable the recognition of hidden patterns embedded in the data flow, and report any anomalous streams of information that indicate possible damage that would be otherwise undetectable. This reflects the underlying, foundational assumption of automated monitoring methods: damage, whether cracks, collapses, local mechanisms, or gradual deteriorations, are present as data anomalies that can be recognized. Whenever data analysis detects possible outliers, dedicated algorithms may be employed for damage localization and for evaluating damage type and extent. Then, in the decision-making phase, the remaining lifetime of the structure is assessed, and expert engineers can promptly establish whether or not to restrict access to the facility as a precautionary measure.

Artificial Intelligence (AI) has a prominent role in this context since it provides various solutions to the problems posed by SHM. In particular, Machine Learning (ML) methods are capable of recognizing atypical patterns in the datasets via statistical comparison of near real-time signals and an adequately large set of training data that the algorithms use to infer the processes that generate "healthy" and, if available, "damage" scenarios.

In the SHM arena, model-based algorithms are often adopted on their own or combined with datadriven algorithms for damage detection. Such methods provide for the construction of a Digital Twin (DT) of the examined structure, so that data collected by sensors can be used as forcing terms for Finite Element Method (FEM) dynamics numerical simulations. The more accurate the model in terms of geometrical, environmental and physical properties, the more reliable the predictions made by the algorithm.

Although the idea of combining data-driven and model-based approaches is definitely not new, such hybridization still represents a challenging issue for Civil and Computer Engineering as, without a unified framework, there is nothing to act as a mortar between the disciplines.

In this work, we propose a comprehensive solution for SHM, which has the advantage of being conceived as a chain of tools and algorithms designed and optimized specifically to perform their unique tasks. The Scientific ToolKit for Opensees (STKO) GUI for OpenSees [START_REF] Petracca | Stko user manual[END_REF]] is employed for the DT modeling and dynamics, as it offers a Python interface for customization, which conveniently connects the model-based and data-driven approaches, of which the data-driven approach is entirely developed in Python. Moreover, a unique data format (the HDF5 format) is adopted for I/O operations to avoid data incompatibilities.

THE MONITORING PROCESS

Currently, the main issues that arise in the automated monitoring process regard the establishment of a robust framework that covers all the intermediate steps from data acquisition to output production and interpretation. The Civil Engineering community offers many suggestions that may fit very specific cases [START_REF] Tibaduiza Burgos | Damage identification in structural health monitoring: A brief review from its implementation to the use of data-driven applications[END_REF]], but an objective standard for SHM is still lacking. In particular, the pre-processing sequence of manipulations that transform rough data time-series into a suitably small set of features to be processed via AI tools represents the most controversial part of the entire monitoring problem. Differences among the various solutions may already emerge during data collection, as different devices register different physical quantities at a given sampling rate.

The information embedded in the time series recorded for a fixed time window is typically extracted as a reduced set of representative features, as the whole time series would be intractable from the AI point of view. Such features can be divided into global, if extracted through multivariate analysis by considering data of a given time window as a whole, with the typical aim of estimating modal characteristics, or local, if data are unpacked before calculating temporal and/or spectral indicators for each time series. Further pre-processing steps, including data standardization and projection onto a subspace of lower dimension, may eventually be performed to eliminate redundant information, taking into account various evironmental effects and preparing the input for the ML algorithm.

Conversely, there is a general agreement about the ML techniques recommended for supervised and unsupervised learning; we note that we refer to unsupervised learning whenever the set used for AI training contains only data collected from structures in full, healthy, operational conditions, while supervised learning requires the training set to contain both regular and anomalous data. We present a broad, automated solution that self adapts in a smart way to the widest set of structures equipped with sensors, covering all possible environmental and external conditions and providing unsupervised as well as supervised solutions. In our scheme, data-driven and model-based algorithms are designed to coexist; moreover, the two strategies are mutually reinforcing, meaning that each of them compensates for the limitations of the other if considered alone for monitoring purposes. Nevertheless, the data-driven part keeps its supremacy within the proposed framework, as it can be run regardless of the realization of a digital twin of the monitored structure.

As we intend for our algorithm to catch all possible damages according to the broadest definition given by [START_REF] Farrar | An introduction to structural health monitoring[END_REF], we constructed two different AI solutions that work in parallel, to independently detect both local and global damages, if any. In the following sections, we outline our monitoring solutions. Our very first input consists of 3-dimensional time series of linear accelerations, angular velocities and inclinations, all sampled at the frequency of 1 KHz 1 , while the temperature value is appended to the features list at a later stage in order to mitigate its effects on the data. Optimal sensor placement is taken for granted in this context; optionally, filters for noise removal can be applied to time series before any other data processing method.

THE DATA-DRIVEN ALGORITHM

Here we introduce the portion of the framework dedicated to data analysis and anomaly detection. In what follows, when referring to data we intend the time series of monitored quantities collected within a fixed time window. The length of a single 1 The unit adopted for data acquisition is the MonStr device provided by ASDEA Hardware. time window is a user-defined parameter, although it is recommended to adopt a value greater than or equal to 100 times the fundamental period of the reference structure for optimal modal parameter estimation [START_REF] Ubertini | Automated modal identification in operational conditions and its application to bridges[END_REF]]. Time windows can be immediately consecutive or partially overlapping.

Univariate analysis. Local detectors

For a given time window: [START_REF] Peeters | Stochastic system identification for operational modal analysis: a review[END_REF]] to extract natural frequencies, damping ratios and modal shapes. Multivariate statistics is also applied to get classical indicators, again using the dedicated functions the TSFEL library makes available to perform this task.

Feature reduction

Regardless of feature arrangement, feature projection on a subspace of lower dimension is required to get rid of data redundancy, and account for environmental effects. Linear Principal Component Analysis (PCA) is adopted to accomplish this task after performing feature rescaling (through simple standardization or min-max scaling). The resulting features, which retain a fixed percentage of total variance, consist of linear combinations of the original ones and are mutually independent. Once computed for the training dataset, the same mathematical objects that perform data rescaling and projection are employed for test-data pre-processing to maintain data consistency.

Unsupervised training

In this case, training data are acquired from physical structures in an entirely healthy state so that AI has no chance to learn how to directly identify possible anomalies embedded in the data. The user can choose from among various algorithms belonging to the family of artificial neural network autoencoders (AE) best suited to their needs [Bank et al. (2020)]. The logic of such deep learning methods lies in the capacity to suitably reconstruct the inputs produced by the same process which generates the training instances, and badly reconstructing any instance whose underlying, production process differs from the "healthy" one. This is realized by inferring the statistics of reconstruction error from training data, that is the normalized ℓ 2 distance between the input and its reconstructed counterpart, then introducing a threshold to distinguish regular from anomalous instances (see Figure 1). Also, anomalous trends can be identified by keeping track of reconstruction errors over time, to monitor slow parameter variations that correspond to structure deteriorations.

Supervised training

In this case, training data can be split into two or more subsets, one corresponding to data coming from the structure in the healthy state, while the others concern various damaged scenarios. In its more basic form, we only distinguish between damaged/undamaged instances; however, if more accurate information is provided, the damaged subset may be furtherly split into smaller subsets according to damage location and severity (see e.g. the approach used in [START_REF] Parisi | Automated location of steel truss bridge damage using machine learning and raw strain sensor data[END_REF]). The production of damaged datasets typically requires the construction of the digital twin of the structure, which The input data x is mapped to the lower dimensional vector z through an artificial neural network (the encoder); then, x is obtained from z through the action of another artificial neural network that is symmetric with respect to the encoder (the decoder). The training phase aims to let the parameters of the network converge to the values that make x similar to x in some suitable metrics.

enables the user to generate dynamics for regular and custom damaged settings. Supervised ML algorithms are preferred in this case to exploit the flagged information contained in the data, bypassing any statistical study that is needed a posteriori in the unsupervised case. Variations of the multilinear perceptron (see Figure 2) are advised with a softmax function as an activator for the output layer in order for the output to directly classify the input instance as regular or anomalous. If in the training phase the data has been further refined, i.e., the output is non-binary, damage location and severity can also be inferred [START_REF] Parisi | Automated location of steel truss bridge damage using machine learning and raw strain sensor data[END_REF]].

Alongside Deep learning algorithms, ensemble learning methods like random forests [START_REF] Ho | Random decision forests[END_REF]] are successfully used in multiclass classification problems with improved accuracy with respect to artificial neural networks [Breiman (2001)], although the training phase is usually costly from a computational point of view (see Figure 3). However, the computational time taken by the two strategies for the inference phase is generally comparable. The use of random forests may also prevent algorithm overfitting [START_REF] Zhou | Structural damage detection method based on random forests and data fusion[END_REF]].

THE MODEL-BASED ALGORITHM

When considering the model-based situation alone, we refer to the model-based solution as the algorithm that uses dynamic data acquired by sensors as input forces for the discrete-time evolution of the FEM model corresponding to the structure under examination (see Figure 4). The idea is to let the model evolve according to the input data, periodically checking its overall conditions. This approach has the advantage of providing a direct, clear response compared to a data-driven solution that needs statistical interpretation of the re-ceived output, particularly when additional information beyond a rough reporting of anomalies is requested. The user can also visually inspect the behavior of the structure under given loads, and interpretation of results is thus immediate. The potential and limitations of this method, i.e., the reliability of the outputs, are totally related to the accuracy with which the numerical model reproduces the physical structure in terms of geometry, material, physical properties, and meshing. Moreover, as such properties are mutable over time and affect structural dynamics, a periodic updating of the model is requested in order for the actual facility and its digital twin to dynamically match. Our solution uses genetic algorithms to minimize the appropriately chosen metric that accounts for differences among natural frequencies and modal shapes (see [START_REF] Levin | Dynamic finite element model updating using simulated annealing and genetic algorithms[END_REF]). The tunable parameters for minimization are the elastic moduli of FEM elements.

We adopted STKO for the model-based portion as it is a straightforward solution to meet our need to maximize software and I/O interaction, as STKO's Python scripting interface permits the user to customize and program the pre and postprocessors as needed. Moreover, STKO includes all materials, elements, conditions and interactions offered in OpenSees, and also hosts all the codes needed for the AI software to perform.

MAIN ADVANTAGES OF OUR HYBRID

SOLUTION Our proposal for completely automated monitoring represents a step towards the establishment of a well-functioning standard for SHM, without reinventing the wheel of digital-twinning. In the following bulleted list, we highlight the pros of adopting our hybrid method, considering that our machinery becomes full data-driven in the absence of a digital replica of the structure. Furthermore, as with every digital twin framework, the numerical model can be used for "what-if" analysis and, mainly, to enrich the training set with custom damaged datasets, thus unlocking access to supervised methods.

Our improvements

• Double-check anomalies detected by the datadriven solution: in order to prevent false alarms, potential anomalies are verified by running the corresponding numerical dynamics and noticing if damage is really produced. • Damage type, localization and extension: if a binary classifier is employed as an AI detector, a numerical simulation provides missing information regarding damage characterization, suggesting possible actions to the user (closure of the facility, possible repairs). • Software compatibility: the entire workflow that realizes the sequence of operations described in the previous sections is perfectly embedded in the STKO program. The hardware and software solutions are designed by the same developers to accomplish a precise, specific task and be integrated one with the other. Hence, there is no need to integrate software programs produced by different companies. • Data format compatibility: in the various steps, all data are organized and exchanged in the HDF5 format. The central node itself rearranges measurements collected by distributed devices in an HDF5 file. This way, there is a complete inner coherence that prevents any data loss that often occurs in the process of conversion to different data formats, that is also time-wasting. • Smart use of CPUs and GPUs: dedicated, powerful graphic cards perform parallel computing for the AI part as the TensorFlow Python library for deep learning distributes by default independent processes on the GPU. The TSFEL library for features extraction gives the user the opportunity to split processes on available threads for a quicker response.

Concerning the model-based algorithm, STKO supports multi-processing solvers of OpenSees. Figure 5 depicts the workflow described in this section.

CONCLUSIONS

Digital twin implementation stands at the frontier of research in Civil and Computer Engineering, thus far representing the more reliable solution for automating the damage detection process for physical structures. We proposed our framework for letting data-driven and model-based approaches coexist and mutually reinforce each other, avoiding typical issues encountered in the realization of the scheme and motivating our algorithmic choices. The use of specifically designated software improves our solution's stability and strength, while the adoption of a unique data format for the whole process eliminates any possible incompatibilities or data loss. In the future, we will spend more time challenging our algorithm with simulated and real data and checking the accuracy of the outputs obtained. Focus will be placed on developing more user options so that the machine can adapt perfectly even to very unique structures requiring monitoring services.

Figure 1 :

 1 Figure 1: Schematic representation of an autoencoder.The input data x is mapped to the lower dimensional vector z through an artificial neural network (the encoder); then, x is obtained from z through the action of another artificial neural network that is symmetric with respect to the encoder (the decoder). The training phase aims to let the parameters of the network converge to the values that make x similar to x in some suitable metrics.

Figure 2 :Figure 3 :

 23 Figure2: Schematic representation of a multilayer perceptron. The input instance is mapped to a binary output for a classifier that distinguishes between regular and anomalous inputs, and has more than two neuronal components for a finer anomaly classification.Training Data

Figure 4 :

 4 Figure 4: Example of a FEM model of a structure modeled within STKO. The figure shows different mode shapes.

  ., Koenigstein, N., and Giryes, R. (2020). "Autoencoders." arXiv preprint arXiv:2003.05991. Barandas, M., Folgado, D., Fernandes, L., Santos, S., Abreu, M., Bota, P., Liu, H., Schultz, T., and Gamboa, H. (2020). "Tsfel: Time series feature extraction library." SoftwareX, 11. Breiman, L. (2001). "Random forests." Machine learning, 45. Buckley, T., Bidisha, G., and Pakrashi, V. (2022). "A feature extraction & selection benchmark for structural health monitoring." Structural Health Monitoring.

Figure 5 :

 5 Figure5: Flowchart of the process from data acquisition to output emission and decision making. Dotted lines refer to the flow of data acquired from the sensor network, dashed lines refer to connections among algorithms related to the data-driven portion, continuous lines to model-based processes. Continuously recorded data feeds data-driven methods, and numerical dynamics come into play whenever an alarm is emitted. If damage signaled by the alarm is confirmed, its extent and location is also evaluated. Periodically, modal parameters are extracted and employed to update the digital twin, which can be also used to perform a "what-if" analysis.


Reda Kasri

Eric Klumperink

Philippe Cathelin

Eric Tournier

Bram Nauta

A Digital Sine-Weighted Switched-Gm mixer for Single-Clock Power-Scalable Parallel Receivers

Keywords: Parallel receivers, Mixer-DAC, Direct Digital Frequency Synthesis, Harmonic rejection, N-path filter, low power CMOS design, scalable power consumption

This paper presents a mixed A/D architecture for parallel channelized RF receiver applications. Its power consumption scales with the number of active receivers and hence with the available overall data rate. A digital sine-weighted switched-Gm mixer with a DDFS per channel is proposed as a zero-IF mixer. The DDFS of all channels are programmable via a Look up Table and are driven by a single central clock. Each channel also exploits a 2-path filter to increase selectivity and interference robustness. To demonstrate the concept two parallel receivers were implemented in a 28nm UTBB FD-SOI CMOS, with 9.5mW/receiver, achieving 40dB of dynamic range, 13dB NF and better than 75 dB inter-receiver isolation.

I. INTRODUCTION

During the last decade, a trend has been to integrate ever more concurrently working parallel RF receivers on the same chip. This has been initiated with TV cable tuners in the 50-1000MHz range, and is now a key feature for several applications, such as carrier aggregation for 4G cellular communications, channelized satellite receivers and upcoming 5G applications. When implemented using traditional analog TV tuners, the two key challenges related to the wide bandwidth are: clean LO-generation for multiple simultaneous frequencies and harmonic mixing. In [START_REF] Gatta | An Embedded 65 nm CMOS Baseband IQ 48 MHz-1 GHz Dual Tuner for DOCSIS 3.0[END_REF], the LO-generation problem was addressed by combining Direct Digital Frequency Synthesis (DDFS) exploiting a DAC and band-pass filter in the reference path of each PLL. For rejecting harmonics, a Harmonic Rejection Mixer (HRM) along with an RF tracking filter were implemented. Each tuner requires its own PLL and inductor VCO, so that many parallel analog tuners would take a lot of chip area, while also introducing interference issues between the different uncorrelated frequencies.

Digitizing the whole spectrum, also known as Full Spectrum Capture (FSC) [START_REF] Wu | A 2.7 mW/Channel 48-1000 MHz Direct Sampling Full-Band Cable Receiver[END_REF], moves the channel selection and processing to digital. It allows a high number of received channels on the same chip, while only one clock is required. However, this solution puts a heavy burden on the ADC and digital signal processing, leading to a high power consumption. Moreover, the ADC power consumption is fixed and high as it always samples the entire band, even if only a few channels are needed, which is the case in many use scenarios. If the spectrum is continuous and at a rather low frequency (such as cable TV in the 50-1000MHz band), FSC may still be effective, but it becomes unpractical and inefficient when the spectrum becomes sparse or when only a few channels are needed. A more power efficient and power-scalable solution is wanted, fitting to a range of data demands, while also addressing the problem of multiple clock generation and interference between the channels.

The architecture we propose in this paper relies on a single central clock that drives the DDFS of each channel (see fig. 1), similar to [START_REF] Maxim | A DDFS Driven Mixing-DAC with Image and Harmonic Rejection Capabilities[END_REF]. However, our focus is on upcoming multichannel moderate performance (~40dB dynamic range) parallel reception applications, whereas [START_REF] Maxim | A DDFS Driven Mixing-DAC with Image and Harmonic Rejection Capabilities[END_REF] targets a traditional >60dB dynamic range TV application (>60dB worst case Harmonic rejection with a 10 bits DDFS). Our goal is to demonstrate that the solution in fig. 1 is efficient compared to FSC, as its power dissipation scales with the number of received channels. Moreover, it is different from [START_REF] Maxim | A DDFS Driven Mixing-DAC with Image and Harmonic Rejection Capabilities[END_REF] at circuit level to make it more power efficient. Due to the relaxed dynamic range and harmonic rejection constraints, an 8-bit DDFS plus mixer-DAC is sufficient, which can be realized compactly at low power (relaxed component mismatch requirements). The spectrum is split in channels of 5-50MHz, so that power efficient ADCs with <10fJ/conversion can be used, allowing for a few times oversampling and relaxed baseband filtering (<3rd order), while the ADC consumption can still be below 1mW for a 100MS/s 10bit ADC [START_REF] Murmann | ADC Performance Survey 1997-2016[END_REF].

The target of this paper is to demonstrate that multiple I/Q receivers including DDFSs and mixer-DACs can be realized with a single master clock, at <10mW/Channel with a channel dynamic range >40dB, which will be sufficient for many upcoming applications. To demonstrate the parallel receiver potential, we have integrated two receiver chains to test mutual interference, and have exploited FDSOI 28nm technology for the design, as it has low power signal processing capabilities and good MOSFET characteristics even at short channel lengths. Fig. 2(a) presents the proposed implementation of the mixer-DAC. It exploits 7 binary scaled transconductances that are all driven by the same RF-voltage (node In), while summing their currents at the output. An 8th bit directs the output current to either the 1 I+ or 1 I-output, i.e. as a sign bit (see below). The output current is thus the product of the analog input voltage and the digital code (DDFS output) that varies over time as a sampled walking sinewave, realizing a multiplying-DAC behavior. The DDFS is implemented here just by an SRAM and an address counter clocked by a single central clock. There is no regular output DAC, as in a traditional DDFS, but instead this DAC is merged with the mixer, i.e. it acts a multiplying DAC.

II. MIXER-DAC IMPLEMENTATION

The mixer essentially is single ended simplified version of the switched transconductor mixer proposed in [START_REF] Klumperink | A CMOS switched transconductor mixer[END_REF], with binary weighting added. This weighting is implemented putting identical switched-Gm mixer slices in parallel, one Gm slice is shown in fig. 2(b). If enabled, it operates as a linear and low noise CMOS invertertransconductor, which has favorable properties like high linearity, class-AB behavior and current re-use. bit and bit are switched between VDD and ground, allowing for low switchresistance regardless of the RF-input signal. To minimize power, a single ended switched transconductor is used, but a differential output voltage is created, similar to what happens in a single balanced analog mixer. The differential load impedances consist of two capacitors, C I+ and C I-connected to the 1 I+ and 1 I- baseband outputs. Intuitively, bit 8 can be seen as a sign-bit but also as the LO-signal of a mixer or chopper, see fig. 2 (c). Indeed, this switch acts like a chopper for the flicker noise of the Gmblocks. Moreover, the two switched capacitors act like an N-path band pass filter [START_REF] Ghaffari | Tunable High-Q N-Path Band-Pass Filters: Modeling and Verification[END_REF], i.e. they present a high impedance to the Gm at the channel-dependent switching frequency of bit 8 and a low impedance elsewhere. Only the wanted signal therefore develops significant voltage swing at the output of the Gm blocks. This improves the robustness of the receiver for other interfering channels non-linearity. To improve the out-of-band signal shunting to ground, a low switch resistance is wanted, thus the bit 8-transistors are large. However, as their output is at baseband, and the parasitic capacitance is in parallel to big baseband capacitors, this is not problematic. The number of DAC bits affects the theoretical maximal output SNR, i.e. SNR out 6.02 1.76 dB, 49.9 dB for an 8 bits DAC. Two parallel tuners receiving two uncorrelated frequencies were implemented. The DDFS concept allows to create any LO frequency compliant with Nyquist theorem from a single central clock: fs>2fmax (fs=2GHz in our case). To have experimental freedom to try different waveforms, we implemented a rather big 64Kb SRAM in this demonstrator. To reduce mismatch, we used a common centroid layout for the Gm bank (fig. 3), Gm dummies and a set of decoupling capacitor rows between each two Gm rows. To reduce clock phase mismatch between the different DAC cells, a balanced clock tree is used, For I/Q phase mismatch reduction, one counter was used to address both the I and Q SRAMs, and only one driver is used to drive both the I and Q clocks. A serial to parallel interface was used to set SRAMs data. Input matching was realized through an external matching resistance.

III. MEASUREMENTS

The presented measurement results are at PCB connector level and thus include losses and noise coming from PCB and off-chip matching-resistance and using an external clock of 2GHz. Fig. 4 shows the measured gain versus LO-frequency. It drops gradually from 33dB at 0.1GHz to 23dB at 1GHz due to charge sharing between the baseband capacitors and the parasitic output capacitance of the Gm-blocks, the NF increases accordingly from 7dB at 0.1GHz to 13dB at 1GHz. In fig. 5 the NF and gain versus IF-frequency is depicted. Thanks to the 8th bit chopping, the flicker noise corner frequency is measured at 20kHz instead of the usual hundreds of MHz in short length technologies such as 28nm FDSOI. The 1st order Gm-C low pass filtering shape is also depicted in the same figure, with a 3MHz cut off frequency, as expected for a 100pF baseband capacitor and 3.7mS average Gm. Harmonic Rejection Ratio (HRR) remains above 41dB for the 100MHz-1GHz spectrum band, as shown in fig. 6. Fig. 7 shows the IIP3 versus offset frequency: it varies from +2 dBm for 6MHz offset, to +16 dBm for 900MHz offset, which bears out the Gm linearity enhancement due to N-path filtering, mentioned earlier. The inter-receiver LO coupling rejection, versus the difference between the two receivers' bit8 frequencies is shown in fig. 8, Rx1 LO (LO1) was set to 100MHz and Rx2 LO (LO2) to 100MHz+Δf. Two tones were set at RF input: fin1 at LO1+80kHz and fin2 at LO2+90kHz. The rejection ratio is then the ratio between the output signal of Rx1 at 80kHz and the output of the same receiver at 90kHz that is down-converted due to the coupling of LO2 to Rx1.

The inter receiver LO coupling rejection shows a better than 75dB rejection. That, combined to the use of a single central clock, makes the proposed solution potentially suitable for parallel reception, and demonstrator performance/consumption can be extrapolated to an N-parallel tuner. Table .1 compares the demonstrator to the state of the art receivers for the same frequency range. The achieved noise figure is in the same range than for others, but as expected/targeted the clean dynamic range is less (roughly >40dB versus >60dB for others) to reduce power consumption. Each I/Q receiver consumes 9.5mW: 5mW for analog (Mixer DAC plus drivers) and 4.5mW for digital (SRAM plus counter), which is 26x-34x lower than for the high performance designs. Note that this power dissipation scales with the number of active channels, while for FSC the 300mW is always needed (for the ADC power alone!). Hence, our proposed solution is certainly more power up to a certain limit (here 26 channels). However, the exact break-even point depends on factors like required dynamic range, ADC energy efficiency and peak-to average ratio considerations, and is subject of a further work. IV. Conclusion

This paper motivates and demonstrates an architecture to realize power efficient reception in channelized receivers using only one central clock. Unlike Full Spectrum Capture, the power consumption is not constant, but scalable with the number of channels and hence with total data-rate. The architecture exploits a digital sine-weighted switched-Gm Mixer-DAC for power efficiency, followed by current mixer with embedded 2-path filter for improved interference robustness and low 1/f noise. The digital sine is generated via a look-up table that drives the Mixer-DAC. Two 6MHz receiver channels are realized on a demonstrator chip covering 100-1000MHz, achieving 23-33dB Gain, 7-13dB NF, >42 Image Rejection, >41dB worst case Harmonic Rejection, and better than <-75dB cross-talk spurs between the channels. Each channel consumes <9.5mW, 5mW for analog and 4.5mW for the SRAM and counter (digital sine generation). This makes the solution suitable for power critical applications with variable data-rate requirements.

Figure 1 :

 1 Figure 1: DDFS-based power-scalable and duplicable receiver.

Figure 2 :

 2 Figure 2: implemented parallel two-channel receiver (a); unit Gm slice with bias resistance (b); SRAM data (c)

Figure 3 :

 3 Figure 3: 7 bits Gm-DAC common centroid Gm bank layout including 127 unit gm blcok and decoupling capacitors

Figure 4 :

 4 Figure 4: measured NF and conversion gain versus equivalent LO-frequency (bit 8), for a constant IF-frequency=500kHz

Figure 5 :

 5 Figure 5: measured NF and conversion gain versus IF frequency at constant equivalent LO-frequency(bit 8) =100MHz

Figure 6 :Figure 7 :

 67 Figure 6: measured harmonic rejection ratio: 3rd, 5th, 7th and 9th harmonics

Figure 9 : 1 mm²

 91 Figure 9: 1 mm² Die micrograph of the implemented two parallel tuners.

Table 1 :

 1 Comparison table, with multi-channel receivers

		This work	[1]	[2]	[3]
	Architecture	low IF mixer DAC	low IF mixer	Full spectrum Sampling	low IF mixer DAC
	Frequency range [MHz] 100-1000	48-1000	48-1000	48-1000
	Number of channel	2	2	158	1
	Power Consumption [mW]	19*	650	300**	450
	Power Consumption /Channel [mW]	9.5*	325	1.9**	450
	Scalable power consumption	yes	yes	No	yes
	Gain [dB]	23-33	-	-	-
	IIP3 [dBm]	2-16	-	-	20
				25.4(W/o
	NF [dB]	7-13	5-7	LNA, W/	13
				digital)
			65 (with RF
	HRR [dB]	41<….<51	tracking	-	60
			Filter)
	Image Rejection [dB]	42	62	70	65
	Second receiver spurs
	on receiver1 rejection	<-75	-	-	-
	[dB]
	Total Die area [mm²]	1	10	1	1.2
	IC Technology	28nm FDSOI	65nm CMOS	28nm CMOS	0.13um CMOS
	* including digital DDFS

**only ADC power consumption (W/o LNA nor digital)

ACKNOWLEDGMENT

The authors thank Andreia Cathelin for technical discussions, Sylvain Clerc for precious help in digital synthesis and Julien Morel for help in measurements. Authors would like also to thank ST packaging and PCB teams.