Emilie Chouzenoux
email: emilie.chouzenoux@inria.fr

Víctor Elvira
email: victor.elvira@ed.ac.uk

Sparse Graphical Linear Dynamical Systems

Keywords: State-space models, graph inference, sparsity, graphical lasso, expectationmaximization, proximal algorithm ©2023 Emilie Chouzenoux

Introduction

Time series appear in applications of most fields of science and engineering, ranging from biomedicine, Earth observation, or network analysis, to name a few. The literature of timeseries analysis in the last century is vast, mostly dominated by auto-regressive movingaveraging (ARMA) models and their extensions [START_REF] Hamilton | Time series analysis[END_REF].

With the surge of the Bayesian paradigm, state-space models (SSMs) have became very popular in the last decades [START_REF] Hamilton | State-space models[END_REF][START_REF] Kim | State-space models with regime switching: classical and Gibbssampling approaches with applications[END_REF][START_REF] Sarkka | Bayesian Filtering and Smoothing[END_REF] for time series modeling. SSMs characterize complex systems through discrete-time models composed of a hidden (or latent) state that evolves in a Markovian manner. SSMs are composed of a state model, which can mimic realistically the complicated dynamics of the system, and the observation model, which links the temporal observations to the hidden state at the same time step. In SSMs, the Bayesian filtering task consists on computing the posterior probability density function (pdf) of the hidden state at a given time step, given all observations from the time series available up to that time. However, in most SSMs the posterior pdf is intractable and must be approximated, generally through particle filters [START_REF] Djuric | Particle filtering[END_REF][START_REF] Doucet | A tutorial on particle filtering and smoothing: Fifteen years later[END_REF][START_REF] Naesseth | Elements of sequential monte carlo[END_REF]. One relevant exception is the linear-Gaussian state-space model (LG-SSM), which allows for exact inference, when the model parameters are known, through the Kalman filter and the Rauch-Tung-Striebel (RTS) smoother (Sarkka, 2013, Chapter 8). The LG-SSM is arguably the most popular model and is still subject of intense research (Sarkka, 2013, Chapter 4). The literature of parameter estimation for LG-SSMs considers both probabilistic and point-wise approaches. The former methods can be applied for a wider class of SSMs (i.e., beyond LG-SSMs) and include for instance particle Markov chain Monte Carlo [START_REF] Andrieu | Particle Markov chain Monte Carlo methods[END_REF], sequential Monte Carlo squared [START_REF] Chopin | Smc2: an efficient algorithm for sequential analysis of state space models[END_REF], and nested particle filters [START_REF] Crisan | Nested particle filters for online parameter estimation in discretetime state-space markov models[END_REF][START_REF] Pérez-Vieites | Nested gaussian filters for recursive bayesian inference and nonlinear tracking in state space models[END_REF] (see [START_REF] Kantas | On particle methods for parameter estimation in state-space models[END_REF] for a review).

One of the main challenges in SSMs is the estimation of the model parameters. In its canonical version, the LG-SSM is composed of four parameters. The state model is parametrized by the so-called transition matrix, which linearly links consecutive states, and the covariance matrix of the additive Gaussian state noise. The observation model is parametrized by the observation matrix, which connects (also linearly) the state with the observation at the same time-step, and the covariance matrix of the additive Gaussian observation noise. Point-wise approaches are generally based on maximum likelihood estimation of the model parameters in LG-SSMs (Sarkka, 2013, Chap. 12). A first family of methods implement recursive schemes to compute first and second derivatives of the likelihood with respect to the unknown parameters (Segal andWeinstein, 1988, 1989;[START_REF] Nagakura | Computing exact score vectors for linear Gaussian state space models[END_REF][START_REF] Gupta | Computational aspects of maximum likelihood estimation and reduction in sensitivity function calculations[END_REF] (see the discussion in (Cappe et al., 2005, Sec. 10.2.4)). Iterative optimizers such as quasi-Newton [START_REF] Olsson | State-space models: from the EM algorithm to a gradient approach[END_REF] or Newton-Raphson [START_REF] Gupta | Computational aspects of maximum likelihood estimation and reduction in sensitivity function calculations[END_REF] are then employed to get a maximum likelihood estimate. However, due to the intricate form of the likelihood function, these optimization solvers must be used in a black-box setting, which prevents from the implementation of efficient linesearch procedures, and limits considerably the convergence guarantees. The second family of point-wise methods is based on expectation-minimization (EM) algorithms [START_REF] Shumway | An approach to time series smoothing and forecasting using the EM algorithm[END_REF] (Cappe et al., 2005, Sec. 10.4) (Sarkka, 2013, Sec. 12.2.3). EM algorithm for parameter estimation in LG-SSM exhibits a good convergence stability [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF][START_REF] Wu | On the convergence properties of the EM algorithm[END_REF], and a simplicity of implementation, which might explain its wide use on practical applied fields, such as finance, electrical engineering, and radar [START_REF] Sharma | Blind Kalman filtering for shortterm load forecasting[END_REF][START_REF] Sharma | Recurrent dictionary learning for state-space models with an application in stock forecasting[END_REF][START_REF] Frenkel | Recursive expectation-maximization (EM) algorithms for timevarying parameters with applications to multiple target tracking[END_REF]. The benefits and drawbacks of these algorithms are discussed in (Shumway and Stoffer, 1982, Sec. 1).

Another important family of approaches for time series analysis relies on graphical modeling. Static models for representing correlation properties within multivariate times series (assuming the series element as independant samples from a random variable), have been explored for decades [START_REF] Maathuis | Handbook of Graphical Models[END_REF]. Let us mention the famous Graphical Lasso algorithm [START_REF] Friedman | Sparse inverse covariance estimation with the graphical LASSO[END_REF], and its multiple variants [START_REF] Chandrasekaran | Latent variable graphical model selection via convex optimization[END_REF][START_REF] Benfenati | Proximal approaches for matrix optimization problems: Application to robust precision matrix estimation[END_REF][START_REF] Belilovsky | Learning to discover sparse graphical models[END_REF][START_REF] Bach | Learning graphical models for stationary time series[END_REF][START_REF] Ying | Does the ℓ 1 -norm learn a sparse graph under Laplacian constrained graphical models?[END_REF], for inference of sparse Gaussian graphical models. Modeling (non-independant) dynamical behaviors through (directed) graph representations, and in particular, causality, has also been explored [START_REF] Eichler | Graphical modelling of multivariate time series[END_REF][START_REF] Ioannidis | Semi-blind inference of topologies and dynamical processes over dynamic graphs[END_REF][START_REF] Giannakis | Topology identification and learning over graphs: Accounting for nonlinearities and dynamics[END_REF].

There is actually a tight link between state-space modeling and dynamical graphical modeling, as emphasized in [START_REF] Barber | Graphical models for time-series[END_REF]. For instance, in [START_REF] Sagi | Extended Kalman filter for graph signals in nonlinear dynamic systems[END_REF] the extended Kalman filter (EKF) is enhanced by taking a graphical perspective. [START_REF] Ioannidis | Inference of spatio-temporal functions over graphs via multikernel kriged Kalman filtering[END_REF] takes an approach of jointly estimating latent processes and the topology of a related graph. [START_REF] Alippi | Graph Kalman filters[END_REF]) also adopts a graphical perspective in linear-Gaussian models (or linearized versions of them through the EKF) in order to learn model parameters through deep neural networks. This approach bears links with differentiable particle filters (DPFs) [START_REF] Corenflos | Differentiable particle filtering via entropy-regularized optimal transport[END_REF][START_REF] Chen | Differentiable particle filters through conditional normalizing flow[END_REF]. In particular, the causality (in the Granger sense) within the hidden state of the state space model is directly related to the transition matrix of the model. As emphasized in [START_REF] Elvira | Graphical inference in linear-Gaussian state-space models[END_REF], such matrix has an interpretation as the adjacency matrix of a directed graph.

For the above reason, the approaches for graphical model inference and model parameter estimation are actually sometimes very similar, algorithmically speaking, the difference being mostly related to the interpretation of the parameters and their representation. The graphical modeling brings new insights such as the choice of specific priors (typically, sparsity) and algorithms (e.g., binary edge selection tools). In most works of graphical modeling, sparsity usually plays a key role [START_REF] Brendan | Discovering structure by learning sparse graphs[END_REF]. Indeed, the interpretation of graphical model requires an edge selection procedure, that can be enforced by imposing a sparsity prior at the graph inference stage. A typical choice is the ℓ 1 (i.e., Lasso) penalty [START_REF] Friedman | Sparse inverse covariance estimation with the graphical LASSO[END_REF][START_REF] Meinshausen | High-dimensional graphs and variable selection with the lasso[END_REF][START_REF] Chouzenoux | GraphEM: EM algorithm for blind Kalman filtering under graphical sparsity constraints[END_REF], which has the advantage of being convex 1 . Other types of priors have also been explored in [START_REF] Ying | Does the ℓ 1 -norm learn a sparse graph under Laplacian constrained graphical models?[END_REF][START_REF] Benfenati | Proximal approaches for matrix optimization problems: Application to robust precision matrix estimation[END_REF][START_REF] Chandrasekaran | Latent variable graphical model selection via convex optimization[END_REF][START_REF] Chouzenoux | Graphit: Iterative reweighted l1 algorithm for sparse graph inference in state-space models[END_REF][START_REF] Elvira | Graphical inference in linear-Gaussian state-space models[END_REF][START_REF] Kumar | A unified framework for structured graph learning via spectral constraints[END_REF][START_REF] Hippert-Ferrer | Learning graphical factor models with Riemannian optimization[END_REF], with the aim of imposing specific graph structures (e.g., low rank, bounded spectrum, block sparsity). For static graph modeling, proximal algorithms [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] (also including augmented Lagrangian primal-dual methods [START_REF] Komodakis | Playing with duality: An overview of recent primaldual approaches for solving large-scale optimization problems[END_REF]) are typically the answer of choice for accounting for the various class of priors of interest in graphical inference. Discrete optimization techniques for edge selection have also been explored (see [START_REF] Benfenati | A review on graph optimization and algorithmic frameworks[END_REF] and references therein). Dynamic graph inference has also been performed using proximal tools [START_REF] Ioannidis | Semi-blind inference of topologies and dynamical processes over dynamic graphs[END_REF]. In the particular case of state-space models, as we explained earlier, several methods, for instance based on Newton or EM, are available for parameter inference. But a common drawback of most available methods is that none allow to introduce in a direct manner a graphical aware prior. Newton-based techniques are, by essence, limited to differentiable priors which prevents the use of Lasso and low-rank penalties. EM methods can be modified to incorporate a prior term on the surrogate majorant function. But this usually leads to an M-step without close form, which prevents a directly implementable algorithmic solution. In our recent works [START_REF] Chouzenoux | Graphit: Iterative reweighted l1 algorithm for sparse graph inference in state-space models[END_REF][START_REF] Elvira | Graphical inference in linear-Gaussian state-space models[END_REF], we introduced more sophisticated EM-based schemes to cope with generic convex regularizations. The goal was to estimate the transition matrix within an LG-SSM, adopting a Granger-based graphical interpretation of such matrix. A fully Bayesian approach to estimate the transition matrix was taken in (Cox andElvira, 2022, 2023) where the space of sparse matrices was explored via reversible jump Markov chain Monte Carlo [START_REF] Green | Reversible jump MCMC[END_REF].

Despite this prolific literature, there exists a gap that we propose to fill in this work. Graphical modeling tools are either dedicated to represent static information related to statistical dependence between independent realizations of random variables (e.g., graphical Lasso approach), or to represent dynamic information related to causality behavior between time series samples. For the latter case, time series are either processed directly [START_REF] Ioannidis | Semi-blind inference of topologies and dynamical processes over dynamic graphs[END_REF], or through a state-space modeling including hidden state [START_REF] Elvira | Graphical inference in linear-Gaussian state-space models[END_REF]. However, we are not aware of any joint approach which includes both static and dynamic graphical modeling (hence, two graphs with two distinct purposes), in the context of state-space models. Our contribution lies within this goal, as we describe in the next section.

Summary of our main contributions

In this work, we introduce a joint graphical modeling for representing static and dynamical behaviors in the hidden state of a linear Gaussian state-space model, bridging the gap between the static graphical Lasso model from [START_REF] Friedman | Sparse inverse covariance estimation with the graphical LASSO[END_REF] and the dynamic model from [START_REF] Elvira | Graphical inference in linear-Gaussian state-space models[END_REF]. We then present a novel Bayesian inference method, called DGLASSO (Dynamic Graphical Lasso), to estimate both graphs, under a sparsity prior, jointly with the construction of the filtering/smoothing distribution for the time series. The method relies on an original optimization algorithm, implementing a block alternating proximal algorithm with efficient majorization-minimization inner steps. We establish the convergence of our algorithm using recent tools from nonlinear analysis. We then perform an extensive experimental validation of the proposed model and inference algorithm by means of experiments on synthetic data and real data from weather variability analysis. For reproducibility purpose, the code for DGLASSO algorithm, is available at https://pages.saclay.inria.fr/emilie.chouzenoux/log/DGLASSO_toolbox.zip.

Problem setup

Notation

Bold symbols are used for matrix and vectors. We denote by x 2 = √ x ⊤ x the Euclidean norm of x ∈ R N , where ⊤ states from the transpose operation and R N is the N -dimensional Euclidean space. We also introduce X F , X 2 and tr(X), the Frobenius norm, the spectral norm (i.e., largest singular value), and the trace, respectively, of elements X = (X(n, ℓ)) 1≤n≤N,1≤ℓ≤M ∈ R N ×M . Id N is the identity matrix of R N . S N denotes symmetric matrices of R N ×N . Both R N ×N and S N are Hilbert spaces, endowed with the Frobenius norm • F and the trace scalar product X, Y = tr(XY). S + N (resp. S ++ N ) is the set of N × N symmetric positive semidefinite (resp. definite) matrices of R N . Given a sequence of elements {x k } K k=1 of length K ≥ 1 and size N , we denote each element as x k = (x k (n)) 1≤n≤N and we use the notation x k 1 :k 2 to refer to the subsequence

{x k } k 2 k=k 1 , for 1 ≤ k 1 < k 2 ≤ K.
For convex analysis concepts, we rely on the notation in the reference textbook by [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]. We denote Γ 0 (H) the set of proper, lower semi continuous convex functions from a Hilbert space H to (-∞, +∞] (Bauschke and Combettes, 2017, Chap. 9), and ∂f the subdifferential of f ∈ Γ 0 (H) (Bauschke and Combettes, 2017, Chap. 16). With a slight abuse of notation, we make use of the extended form of the minus logarithm determinant function, defined as

(∀P ∈ S N ) -log det(P) = -log |P|, if P ∈ S ++ N , +∞, otherwise, (1)
with |P| the product of the eigenvalues of P. According to [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF], the function in Eq. ( 1) belongs to Γ 0 (S N ).

Graphs

We introduce here our notation for describing graphs. Most of our notation is inherited from [START_REF] Bühlmann | Statistics for High-Dimensional Methods, Theory and Applications[END_REF].

Let us define a graph G made of set of N vertices V = {v (n) s.t. n ∈ {1, . . . , N }} and of a set of edges E = {e (n,ℓ) s.t. (n, ℓ) ∈ E}. The latter gathers ordered pairs of distinct vertices, and as such, E ⊂ {1, . . . , N } 2 . Undirected graphs are made of undirected edges, that is such that (n, ℓ) ∈ E and (m, n) ∈ E, for every (m, n) ∈ {1, . . . , N } 2 . In contrast, directed graphs consist of directed edges, where we say that some e (n,ℓ)

∈ E is directed (from n to ℓ) if (ℓ, n) / ∈ E.
We can also distinguish reflexive graphs if self-loops are allowed (i.e., one can have (n, n) ∈ E), and nonreflexive graphs otherwise. Given these definitions, one can simply bound the cardinality of E for each category. For instance, for an undirected nonreflexive graph, Card(E) ≤ N (N -1)/2, while a directed nonreflexive graph has Card(E) ≤ N 2 and a directed reflexive graph has Card(E) ≤ N (N -1). Such graph definitions are binary, as it only described presence/absence of edges between vertex pairs. In this work, we require the notion of a weighted graph, where the edges (e (n,ℓ) ) (n,ℓ)∈E are associated to real valued weights (ω (n,ℓ) ) (n,ℓ)∈E . The edge positions and weight values of a graph are summarized in a matrix M ∈ R N ×N , where, for every (n, ℓ) ∈ E, M (n, ℓ) = ω (n,ℓ) and, for every (n, ℓ) ∈ {1, . . . , N } 2 / ∈ E, M (n, ℓ) = 0. An undirected (resp. directed) graph is thus associated to a symmetric (non symmetric) matrix M. A reflexive (resp. non reflexive) graph has non zero (resp. zero) entries in the diagonal of M. The number of edges is simply obtained as the number of non-zero entries in M. We finally define the so-called binary support of the (possibly sparse) matrix M. For every

M ∈ R N ×N , supp(M) = S ∈ {0, 1} N ×N , with, for every (n, ℓ) ∈ {1, . . . , N } 2 , S(n, ℓ) = 0 if and only if M (n, ℓ) = 0 (i.e., (n, ℓ) ∈ E).

Dynamical modeling

We consider the LG-SSM described, for k = 1, . . . , K, as:

x k = Ax k-1 + q k , (2)
y k = H k x k + r k , (3)
where,

• {x k } K k=1 ∈ R Nx and {y k } K k=1 ∈ R Ny , are the hidden state and the observations at each time k respectively,

• A ∈ R Nx×Nx is the transition matrix that we aim at estimating,

• {H k } K k=1 ∈ R
Ny×Nx maps the observation model matrices, possibly varying with k, that are assumed to be known,

• {q k } K k=1 ∼ N (0, Q) is the i.i.d
. state noise process, assumed to follow a zero-mean Gaussian model with covariance matrix Q ∈ S ++ Nx that we also aim at estimating,

• {r k } K k=1 ∼ N (0, R k ) is the i.i.d
. observation noise process, again zero-mean Gaussian with known covariance matrices R k ∈ S ++ Ny .

Throughout the paper, we denote P = Q -1 the precision matrix of the state noise. We furthermore assume an initial state distributed such that x 0 ∼ N (x 0 ; µ 0 , Σ 0 ) with known

µ 0 ∈ R Nx and Σ 0 ∈ S ++ Nx .
The state noises and the observation noises are mutually independent and also independent of the initial state x 0 .

In the next subsection, we provide some reminders about filtering/smoothing procedures for time series described by an LG-SSM.

Filtering and smoothing algorithms in linear dynamical systems

Both filtering and smoothing algorithms consist on the computation of a posterior probability density function (pdf) of the hidden state {x k } K k=1 . For every k ∈ {1, . . . , K}, the filtering distribution is p(x k |y 1:k ), where we denote as y 1:k = {y j } k j=1 the set of observations available up to the time step k, i.e., no future observations can be used to estimate x k . The filtering problem is suitable for online processing of the observations. The smoothing distribution is p(x k |y 1:K ), where K is the final time-step for which there is an available observation, i.e., for k ∈ {1, . . . , K -1} (note that it is also possible to condition on a subset of future observations, e.g., p(x k |y 1:k+τ ) with τ ∈ N ).

The filtering and smoothing distributions are in general intractable for most SSMs of interest. The LG-SSM is one of the exceptions that admit closed-form solutions.

Estimating the filtering and smoothing distributions is in general a challenging problem, since obtaining these distributions of interest is possible only in few models of interest. For instance, for the LG-SSM described in (2)-(3), it is possible to obtain the filtering and smoothing distributions, for k = 1, . . . , K, in the case where the model parameters A, Q, {H k } K k=1 , and {R k } K k=1 are known. Interestingly, these distributions can be obtained in an efficient sequential manner. In particular, the Kalman filter (KF) [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF] allows to obtain recursively (in a forward manner) the sequence of filtering distributions. Its smoothing counterpart, the Rauch-Tung-Striebel (RTS) smoother [START_REF] Briers | Smooting algorithms for state-space models[END_REF], runs backwards to obtain the sequence of smoothing distributions. We note that both algorithms require the model parameters to be known, which in the case of the LG-SSM presented in the previous section are A, Q = P -1 , {H k } K k=1 , and {R k } K k=1 . Algorithm 1 describes KF, which at each time step k ∈ {1, . . . , K} performs the (a) prediction/propagation step, where the mean µ k|k-1 and covariance Σ k|k-1 of the (state) predictive distribution are obtained; and the (b) update step, where the mean µ k and covariance Σ k of the filtering distribution Algorithm 1 Kalman Filter Input. Prior parameters µ 0 and Σ 0 ; model parameters A, P, {H k } K k=1 , and

{R k } K k=1 ; set of observations {y k } K k=1 .
Recursive step. For k = 1, . . . , K (a) Prediction/propagation step.

µ k|k-1 = Aµ k-1 (4) Σ k|k-1 = AΣ k-1 A ⊤ + P -1 (5) (b) Update step. ν k = H k µ k|k-1 (6) v k = y k -ν k (7) S k = H k Σ k|k-1 H ⊤ k + R k (8) K k = Σ k|k-1 H ⊤ k S -1 k (9) µ k = µ k|k-1 + K k v k (10)
Σ k = Σ k|k-1 -K k S k K ⊤ k (11) Output. {µ k , Σ k } K k=1 .
Then, for each k = 1, ..., K:

• state filtering pdf:

p(x k |y 1:k ) = N (x k ; µ k , Σ k )
• observation predictive pdf: p(y k |y 1:k-1 ) = N (y k ; ν k , S k ) are obtained. Algorithm 2 describes the RTS smoother. In this case, the iteration starts at k = K and runs backwards. It can be interpreted as a a refinement from the mean and covariance matrices of the filtering distribution, given by Kalman, updating them with information present in future observations. However, note that the observations are not reused in the RTS algorithm, i.e., all the required information in the observations is absorbed by the filtering distributions, which are used to produce the smoothing distributions.

Problem statement

Algorithms 1 and 2 are simple and efficient. However, they require the knowledge of the model parameters. In this paper, we assume {H k } K k=1 , and {R k } K k=1 to be known, and we address the problem of obtaining the filtering and smoothing distributions when matrices A and P are unknown, and must be estimated jointly with the filtering/smoothing step. To do so, we introduce a double graphical modeling of the state equations, where each matrix A and P now represents the weights of graphs with a specific, and complementary, statistical Algorithm 2 RTS Smoother Input. Filtering parameters {µ k , Σ k } K k=0 from the Kalman filter; model parameters A and P.

Initialization. Set µ s K = µ K and Σ s K = Σ K . Recursive step. For k = K, K -1, ..., 0 µ - k+1 = Aµ k (12) Σ - k+1 = AΣ k A ⊤ + P -1 (13) G k = Σ k A ⊤ Σ - k+1 -1 (14)
µ s k = µ k|k-1 + G k µ s k+1 -µ - k+1 (15)
Σ s k = Σ k|k-1 -G k Σ s k+1 -Σ - k+1 G ⊤ k (16) Output. {µ s k , Σ s k } K k=1 .
Then, for each k = 1, ..., K:

• state smoothing pdf:

p(x k |y 1:K ) = N (x k ; µ s k , Σ s k )
interpretation. We then propose an efficient and convergent inference approach to estimate both graphs under sparse priors given the observed sequence {y k } K k=1 , while also obtaining its filtering and smoothing distributions at every time steps.

Proposed model and inference method

We now introduce our graphical modeling perspective of the LG-SSM. This novel view gives a graph-based interpretation of the transition matrix A and precision matrix P. Then, we present the optimization methodology in order to estimate such matrices by exploiting the graphical perspective.

Sparse dynamical graphical model

Our novel approach interprets the transition matrix A and precision matrix P as a (weighted) directed and undirected graph, respectively. We now particularize this perspective for each matrix, deepening into the interpretability of such novel view, the complementarity of both graphs, and its benefits during the inference process.

State transition matrix

Matrix A governs the hidden process in (2) and can be seen as the matrix parameter of an order-one vector auto-regressive (VAR) unobserved process. For every n ∈ {1, . . . , N x } and ℓ ∈ {1, . . . , N x }, the entry A(n, ℓ) contains the information of how the n-th time series {x k (n)} K k=1 is affected by the ℓ-th time series {x k (ℓ)} K k=1 in consecutive time steps. More precisely, we can express the update n-th dimension of the latent state in the generative model as

x k (n) = Nx ℓ=1 A(n, ℓ)x k-1 (ℓ) + q k (n). ( 17
)
Thus, if A(n, ℓ) = 0, it implies that the ℓ-th time series of the latent state does not provide any information to predict the n-th time series one time step ahead conditioned to observing the information in all time series for which A(n, m) = 0, m ∈ {1, . . . , N x } \ ℓ.

We express this conditional independence by denoting x(n) ⊥ ⊥ x(ℓ)|X n . Here, for every n ∈ {1, . . . , N x }, we denote the set X n = {x(m)|A(n, m) = 0, m ∈ {1, . . . , N x }}, i.e., the prediction of x(n) observing the past of the time series in X n cannot be improved by observing the other time series. Note that in general, one can have A(n, ℓ) = A(n, ℓ) and A(n, n) = 0. Our interpretation is clearly connected to Granger causality [START_REF] Granger | Investigating causal relations by econometric models and cross-spectral methods[END_REF], and more in particular to conditional Granger causality [START_REF] Luengo | Hierarchical algorithms for causality retrieval in atrial fibrillation intracavitary electrograms[END_REF]. In particular, A(ℓ, n) = 0 for some (n, ℓ) ∈ {1, . . . , N x } 2 , it means that x k-1 (n) causes (in conditional Granger sense) x k (ℓ), for every k ∈ {1, . . . , K}. The conditional Granger causality directional relationships within the entries of the multivariate latent state time series x can hence be represented as a graphical model made of a reflexive directed graph with N x vertices, whose weights are those of the transpose matrix A ⊤ . The perspective of matrix A interpreted as a (weighted) directed graph bears some links with the work of graphical Granger causality by [START_REF] Shojaie | Discovering graphical granger causality using the truncating lasso penalty[END_REF], although we here model the interactions in the latent space instead of directly on the observations. In our case the graphical modeling is simpler, since the work by Shojaie and Michailidis (2010) considers all possible interactions in the observation space across time (i.e., the interpreted graph size is K • N y ). The price to pay of our modeling is the difficulty in inferring A, since it governs the latent process, hence it is never observed. We propose an advanced methodology to address the inferential task in Section 3. Illustrative example. In Figure 1, we display an illustrative example of the graphical model associated to the following state equations for N x = 5, for every k ∈ 1, . . . , K, with q k ∈ R Nx the latent state noise:

               x k (1) = 0.9 x k-1 (1) + 0.7x k-1 (2) + q k (1), x k (2) = -0.3 x k-1 (3) + q k (2), x k (3) = 0.8 x k-1 (5) + q k (3), x k (4) = -0.1 x k-1 (2) + q k (4), x k (5) = 0.5 x k-1 (3) + q k (5). ( 18
)
We display the matrix A, the associated binary matrix supp(A) and the resulting reflexive directed graph under this interpretation.

A =      0.9 0.7 0 0 0 0 0 -0.3 0 0 0 0 0 0 0.8 0 -0.1 0 0 0 0 0 0.5 0 0      (a) Matrix A. supp(A) =      1 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 1 0 0      (b) Binary support. 1 2 3 5 4 A(1, 2) A(2, 3) A(4, 2) A(5, 3) A(3, 5) A(1, 1) (c) Associated directed graph.
Figure 1: Graphical model associated to (18). Matrix A (a), its binary support (b) and associated reflexive directed graph (c). The edges are defined as non-zero entries of A ⊤ . Their thickness is proportional to the absolute entries of A ⊤ .

State noise precision matrix

Matrix Q denotes the noise covariance in the state Eq. ( 2). Since the noise is assumed to be Gaussian, this matrix, and more precisely, the associated precision matrix P = Q -1 , also has a direct interpretation in terms of graphical modeling, using the notion of Gaussian graphical model (GGM) (Bühlmann and Van De Geer, 2011, Section 13.4) [START_REF] Uhler | Gaussian graphical models: An algebraic and geometric perspective[END_REF]. Since we consider Q constant during the whole time series, let us denote the multivariate state noise r.v. at any time step as q ∼ N (0, Q). The GGM consists in a graphical modeling of the independence (or not) between the scalar random variables q(1), . . . , q(N x ). It is easy to prove that

q(n) ⊥ ⊥ q(ℓ)|{q(j), j ∈ 1, . . . , N x \{n, ℓ}} ⇐⇒ P(n, ℓ) = P(ℓ, n) = 0, (19)
i.e., the dimensions n and ℓ of q are independent given all other dimensions if and only if the entry P(n, ℓ) is zero (and obviously also P(ℓ, n) since the precision matrix is symmetric). Note that it is possible to condition in the l.h.s. of ( 19) only to the dimensions q(j) for which P(n, j) = 0 and the equivalence would still hold. The GGM relies on an undirected graph associated to a symmetric weight matrix equals to the inverse of the covariance matrix P = Σ -1 . In particular, namely (n, ℓ) / ∈ E if and only if P (n, ℓ) = P (ℓ, n) = 0. This GGM construction is at the core of the famous GLASSO (Graphical Lasso) formulation [START_REF] Friedman | Sparse inverse covariance estimation with the graphical LASSO[END_REF] (Maathuis et al., 2019, Section 9.7), whose goal is to build the maximum a posteriori estimator of P given realizations of the random vector q under a sparsity assumption on matrix P. The sparsity is here interpreted as a way to eliminate spurious edges in the graph associated to P. Illustrative example. In Figure 2, we display an illustrative example on the GGM associated to a given precision matrix P for N x = 5. We show the associated binary support matrix supp(P) and the resulting undirected graph under this interpretation. Although self-loops (i.e., non-zero diagonal elements) occur, we removed them from the graphical representation for ease of readability.

Proposed unifying view

We now summarize the graphical perspective on both A and Q and describe an unifying approach, where sparsity plays a key role. Matrix A is interpreted as the weight matrix of a directed graph with N x vertices. Sparsity (i.e., absence of edge in the graph) in A is interpreted as pair-wise partial/conditional independence, given a subset of the remaining time series, for a one-step ahead prediction of the hidden state. Matrix P = Q -1 is interpreted as the weight matrix of an undirected graph, related to a GGM describing the noise in the latent space. Sparsity in P is interpreted as pair-wise partial/conditional independence of two dimensions of the additive state noise, given a subset of the remaining dimensions. Both graphs are reflexive, with N x nodes and a maximum of N 2

P =      2 0 -0.1 0 0 0 0.9 0.3 -0.2 0.5 -0.1 0.3 0.8 0 0 0 -0.2 0 2 0 0 0.5 0 0 1.5      (a) Matrix P. supp(P) =      1 0 1 0 0 0 1 1 1 1 1 1 1 0 0 0 1 0 1 0 0 1 0 0 1      (b) Support matrix.
x edges for A (N x (N x -1) for P) associated to N 2

x weights (i.e., the N 2 x entries of A or P). Our perspective in the state process of the LG-SSM in ( 2) is that A encodes the way the information flows in consecutive time-steps between the nodes (state dimensions) of the network (vector state). Thus, its properties shape how the energy/information is transferred and dissipated (under the noise). In contrast, P = Q -1 encodes how information that is not in the system at time k -1 enters in the system at time k. In that respect, the interpreted graph with weight matrix P encodes the dependency of the new information across the nodes of the network.

We adopt the above perspective to estimate both A and Q by promoting properties in both graphs. Specifically, we introduce sparsity priors on the matrices, as the sparsity property is key to reach interpretability and compactness of the whole model. In particular, it allows to understand the inner structure of the latent space. Moreover, it can be helpful to speed up computations as the sparsity level is increased, e.g., when running the Kalman filter and RTS smoother. Our proposed method DGLASSO (Dynamic Graphical Lasso) hence aims at providing the maximum a posteriori (MAP) estimator of A and P (i.e., the weight matrices related to the graphical modeling of the latent state correlation and causality) under Lasso sparsity regularization on both matrices, given the observed sequence y 1:K . A visual representation of DGLASSO graphical model is given in Figure 3. The figure summarizes the relationships among the state entries of an LG-SSM using matrices (A, P) from Figures 1 and2.

Related works: Our approach DGLASSO generalizes important existing sparse graphical inference ones. For instance, our model with A = 0 (degenerate case) has no memory, and all the energy/information of the system is lost at each time step, thus the state dimensions only incorporate exogenous energy/information through the additive noises. This degenerate case is the same model than GLASSO [START_REF] Friedman | Sparse inverse covariance estimation with the graphical LASSO[END_REF] in the case when R k ≡ 0, and same than the robust GLASSO model (Benfenati et al., 2020, Sec.5.2) when R k ≡ σ 2 R Id. In contrast, if the state noise covariance matrix Q is known, DGLASSO coincides with our recent GraphEM framework [START_REF] Elvira | Graphical inference in linear-Gaussian state-space models[END_REF]. Probably the closer related work is [START_REF] Ioannidis | Semi-blind inference of topologies and dynamical processes over dynamic graphs[END_REF], which also introduces a joint graph modeling within an LG-SSM, capturing order-one causal relationships and instantaneous influence (i.e., order zero), through two sparse graphs. Their proposed inference method is an alternating optimization technique, that infers the two graphs under Lasso prior, jointly with the estimation of hidden state. In contrast with DGLASSO, in [START_REF] Ioannidis | Semi-blind inference of topologies and dynamical processes over dynamic graphs[END_REF], (i) the state model follows a structural vector autoregressive model (SVAR) where instantaneous causality and noise are distinguished, while DGLASSO assumes an order-one VAR in the hidden state; and (ii) the cost function does not result from a Bayesian modeling, and as such it is not related to a maximum a posteriori loss for the graph variables, (iii) the state estimation is point wise defined as the solution of an handcrafted optimization problem, while DGLASSO preserves a full Bayesian interpretation and hence allows the complete characterization of the filtering/smoothing state distributions. In particular, [START_REF] Ioannidis | Semi-blind inference of topologies and dynamical processes over dynamic graphs[END_REF] model does not recover GLASSO as a particular case.

Figure 3: Summary representation of the DGLASSO graphical model, for the example graphs presented in Figs. 1 and2. Blue (oriented) edges represent Granger causality between state entries among consecutive time steps, encoded in matrix A (Fig. 1). Magenta edges represent static (i.e., instantaneous) relationships between the state entries, at every time step, due to correlated state noise described by matrix P (Fig. 2).

Optimization problem

The considered MAP inference problem reads as an optimization problem that we formulate hereafter. More specifically, let us denote the posterior of the unknown parameter, p(A, P|y 1:K ), where the hidden states have been marginalized. It is direct to show, using Bayes rule and composition with the (strictly increasing) logarithmic function, that the maximum of p(A, P|y 1:K ) ∝ p(A, P)p(y 1:K |A, P), with p(A, P) some prior on the parameters A and P, coincides with the minimum of the following loss function:

(∀A ∈ R Nx×Nx )(∀P ∈ S Nx ) L(A, P) L 1:K (A, P) + L 0 (A, P). ( 20
)
with L 1:K (A, P)log p(y 1:K |A, P) and L 0 (A, P) =log p(A, P). According to (Sarkka, 2013, Chap. 12),

(∀A ∈ R Nx×Nx )(∀P ∈ S Nx ) L 1:K (A, P) = K k=1 1 2 log det(2πS k ) + 1 2 z ⊤ k S -1 k z k , (21)
where z k = y k -H k Aµ k-1 and S k is the covariance matrix of the predictive distribution p(y k |y 1:k-1 , A, P) = N (y k ; HAµ k-1 , S k ), both being obtained by the KF in Alg. 1 run for given (A, P) (see (Sarkka, 2013, Section 4.3)).

As already mentioned, the introduction of priors that induce sparsity is advantageous due to several reasons. First, it generally reduces over-fitting, particularly when K is low compared to the number of parameters to be estimated. Also, it enhances the interpretability. The zero elements in A and P have a clear interpretation of conditional independence between the time series. Ideally, we would use the ℓ 0 norm of A and P, i.e., penalizing the number of non-zero entries of the matrix. However, this norm is known to have undesirable properties such as being non-convex, non continuous, and associated to a improper law p(A). We thus propose instead the regularization term

(∀A ∈ R Nx×Nx )(∀P ∈ S Nx ) L 0 (A, P) = λ A A 1 + λ P P 1 . ( 22
)
The ℓ 1 norm in ( 22), defined as the sum of absolute values of the matrix entries, is a proper convex function, that leads to the so-called Lasso regularization [START_REF] Bach | Optimization with sparsity-inducing penalties[END_REF]. Note that this penalty, used in numerous works of signal processing and machine learning [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF][START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF], including graph signal processing [START_REF] Friedman | Sparse inverse covariance estimation with the graphical LASSO[END_REF][START_REF] Benfenati | Proximal approaches for matrix optimization problems: Application to robust precision matrix estimation[END_REF], is associated with a joint Laplace prior distribution on A and P. Such joint distribution factorizes (i.e., the prior assumes independence on both parameters), the means are zero, and the scale parameters are proportional to, respectively, λ A and λ P . The larger the regularization parameter λ A (or λ P ), the higher sparsity of A (or P), with the degenerate case of a null A (and P) when the regularization parameter grows.

General minimization procedure

The expressions ( 21)-( 22) provide an effective way to evaluate (20). However, due to the recursive form in (21), it is challenging to derive direct quantities (e.g., gradient) for L 1:K . Moreover, despite its simple expression, the regularization term ( 22) is non differentiable. For both reasons, the minimization of ( 20) is a challenging question.

We propose a block alternating majorization-minimization (MM) technique to jointly infer the probabilistic filtering/smoothing distributions of the SSM, along with MAP estimates of (A, P). Our method presents the advantage of sound convergence guarantees and the ability to incorporate sparsity priors on both A and P. The general idea of MM is to replace a complicated optimization problem by a sequence of more tractable ones [START_REF] Sun | Majorization-minimization algorithms in signal processing, communications, and machine learning[END_REF][START_REF] Hunter | A tutorial on MM algorithms[END_REF]. Surrogate approximations of the cost function are built iteratively, following a majorization principle. For any estimates A ∈ R Nx×Nx and P ∈ S ++ Nx (i.e., in the interior domain of definition of L) of (A, P), a majorizing approximation is constructed for the likelihood term. It is required to satisfy both

(∀A ∈ R Nx×Nx )(∀P ∈ S Nx ) Q(A, P; A, P) ≥ L 1:K (A, P), (23)
and also the so-called tangency condition

Q( A, P; A, P) = L 1:K ( A, P). (24)
The MM algorithm then alternates between Majorization step to build function Q + L 0 satisfying conditions ( 23) and ( 24), and Minimization step to minimize this majorizing approximation. In our proposed approach, we adopt a block alternating implementation of the MM method, where only one variable (i.e., A or P) is updated at each iteration, following a cyclic rule [START_REF] Jacobson | An expanded theoretical treatment of iteration-dependent majorize-minimize algorithms[END_REF][START_REF] Hong | A unified algorithmic framework for block-structured optimization involving big data: With applications in machine learning and signal processing[END_REF]. This alternating strategy considerably simplifies the majorizing function construction as well as its minimization. We furthermore add so-called proximal terms in both updates. Let us recall that, for a function φ : H → (-∞, +∞] ∈ Γ 0 (H), with H a Hilbert space with endowed norm • , the proximity operator 2 of function f at V ∈ H is defined as [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] prox

φ ( V) = argmin V∈H φ(V) + 1 2 V -V 2 . ( 25
)
In our context, the considered Hilbert space is either R Nx×Nx for the update of the transition matrix or S Nx for the update of the precision matrix, and the endowed norm is in both cases the Frobenius norm • F . Introducing proximity terms thus amounts to adding to each majorizing function a quadratic distance to the previous iterate, weighted by a positive factor. This modification preserves the MM interpretation of the method, while ensuring improved stability and convergence guarantees. As we show below, the iterates belong to the interior of domain of the loss function by construction. Namely for every i ∈ N, (A (i) , P (i) ) ∈ R Nx×Nx × S ++ Nx , so the precision matrix remains invertible along the iterations and the algorithm is well defined. The resulting DGLASSO approach is summarized in Algorithm 3. DGLASSO aims at providing ultimately the MAP estimates for the matrix parameters (A, P) of the considered LG-SSM, through the minimization of (20). The covariance state noise MAP estimate is straightforwardly obtained by inversion of the precision state noise matrix provided as an output of DGLASSO. The state filtering/smoothing pdf associated to each estimates are obtained by running KF/RTS loops when setting (A, P) equal to DGLASSO outputs. Next sections are dedicated to the (i) construction of the majorizing function, (ii) discussion about the resolution of each inner step, and (iii) convergence analysis.

Building the majorizing function

In this section, we derive a theorem regarding the expression of the loss function L 1:K and a valid majorant function for it.

Theorem 1 The loss function can be expressed as

(∀A ∈ R Nx×Nx )(∀P ∈ S Nx ) L 1:K (A, P) = 1 2 K k=1 (x k -Ax k-1 ) ⊤ P (x k -Ax k-1 ) - K 2 log det(2πP) + log p(x 0:K |y 1:K , A, P) + log p(x 0 ) - K k=1 p(y k |x k ). ( 28
)
2. See also http://proximity-operator.net/ Algorithm 3 DGLASSO algorithm Inputs. Prior parameters µ 0 and Σ 0 ; model parameters {H k } K k=1 and {R k } K k=1 ; set of observations {y k } K k=1 ; hyper-parameters (λ A , λ P ) > 0; stepsizes (θ A , θ P ) > 0; precisions (ε, ξ) > 0.

Initialization. Set (A (0) , P (0) ) ∈ R Nx×Nx × S ++ Nx .
Recursive step. For i = 0, 1, . . .:

(a) //Update transition matrix (i) Build surrogate function satisfying (23) and (24) at A = A (i) and P = P (i) .

(ii) Run Algorithm 4 with precision ξ to solve

A (i+1) = prox A→θ A Q(A,P (i) ;A (i) ,P (i) )+θ A λ A A 1 A (i) , = argmin A∈R Nx×Nx Q(A, P (i) ; A (i) , P (i) ) + λ A A 1 + 1 2θ A A -A (i) 2 F . (26)
(b) //Update noise precision matrix (i) Build surrogate function satisfying (23) and (24) at A = A (i+1) and P = P (i) .

(ii) Run Algorithm 5 with precision ξ to solve P (i+1) = prox P→θ P Q(A (i+1) ,P;A (i+1) ,P (i) )+θ P λ P P 1 P (i) , = argmin i+1) , P; A (i+1) , P (i) ) + λ P P 1 + 1 2θ P P -

P∈S Nx Q(A (
P (i) 2 F . (27)
If A (i+1) -A (i) F ≤ ε A (i) F and P (i+1) -P (i) F ≤ ε P (i) F , stop the recursion by returning (A (i+1) , P (i+1) ).

Output. MAP estimates of the transition and state noise precision matrices.

Moreover, consider the outputs of Algs. 1 and 2 for a given A ∈ R Nx×Nx and P ∈ S N ++ x . Denote

     Ψ = 1 K K k=1 Σ s k + µ s k (µ s k ) ⊤ , ∆ = 1 K K k=1 Σ s k G ⊤ k-1 + µ s k (µ s k-1 ) ⊤ , Φ = 1 K K k=1 Σ s k-1 + µ s k-1 (µ s k-1 ) ⊤ , (29)
where, for every k ∈ {1, . . . , K}, G k = Σ k A ⊤ ( AΣ k A ⊤ + P -1 ) -1 (see Alg. 2). Then, conditions (23) and (24) hold with

(∀A ∈ R Nx×Nx )(∀P ∈ S Nx ) Q(A, P; A, P) = K 2 tr P(Ψ -∆A ⊤ -A∆ ⊤ + AΦA ⊤ ) - K 2 log det(2πP). ( 30
)
As a consequence, for every (θ A , θ P ) > 0, (∀A ∈ R Nx×Nx )(∀P ∈ S Nx ) L(A, P) ≤ Q(A, P; A, P) + L 1:K ( A, P) -Q( A, P; A, P)

+ λ A A 1 + λ P P 1 + 1 2θ A A -A 2 F + 1 2θ P P -P 2 F , (31)
with equality holding for A = A and P = P.

Proof See Appendix A.

Theorem 1 allows to build, for any tangent point ( A, P), a majorizing approximation (31) for L. DGLASSO method leverages this property designing a block alternating MM scheme. At each iteration i ∈ N, two steps are conducted, namely the update of (a) the transition matrix, (b) the noise precision matrix. Each steps follows an MM structure, that is it first builds a majorizing approximation for L at the current estimates, using Theorem 1, and then minimizes it with respect to the active variable (A in step (a), or P in step (b)). Processing the variables in two distinct steps allows to build upon the desirable convex structure of (31) with respect to one of the variable, the other being fixed. The good performance of MM approaches combined with block alternating steps have been illustrated in [START_REF] Hong | A unified algorithmic framework for block-structured optimization involving big data: With applications in machine learning and signal processing[END_REF][START_REF] Chouzenoux | A block coordinate variable metric forwardbackward algorithm[END_REF][START_REF] Hien | Inertial block proximal method for non-convex nonsmooth optimization[END_REF]. In particular, convergence guarantees are at reach, as we will show in Section 4.

Remark 2 From Theorem 1, the construction of the majorizing function of L at each tangent point ( A, P) requires to run the KF Algorithm 1 and the RTS Algorithm 2. This feature was already present in previous EM-based approaches for LG-SSM parameter estimations (see for instance (Sarkka, 2013, Chap. 12) and [START_REF] Elvira | Graphical inference in linear-Gaussian state-space models[END_REF]). Due to its block alternating form, the proposed DGLASSO algorithm requires to build twice per iteration a majorizing approximation, which means that the data must be processed by KF/RTS twice per iteration. In contrast, gradient-based methods relying on sensitivity equations [START_REF] Gupta | Computational aspects of maximum likelihood estimation and reduction in sensitivity function calculations[END_REF] only require the KF recursions for building their updates, but not the RTS recursions, being thus more computationally efficient. However, such approaches cannot easily encompass non-smooth sparsity priors. Moreover, if no prior is assumed (i.e., λ A = λ P = 0), the aforementioned methods exhibit very slow convergence in this problem as we will show in our numerical experiments. This might be because the gradient of the neg-log-likelihood is extremely ill-conditioned in this context.

Resolution of the inner problems

We now discuss the structure and resolution of the inner problems ( 26) and ( 27) arising in Algorithm 3.

Minimization with respect to A: Let A ∈ R Nx×Nx and P ∈ S ++ Nx . By definition of the proximity operator (25) and the majorant expression in (30), Eq. ( 26) amounts to minimizing the following function

minimize A∈R Nx×Nx C 1 (A) θ A K 2 tr P(Ψ -∆A ⊤ -A∆ ⊤ + AΦA ⊤ ) + θ A λ A A 1 + 1 2 A -A 2 F . (32)
Remarkably, the problem above is a special instance of a multivariate Lasso regression problem [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], for which many efficient iterative solvers are available. The specificity here is that the problem is strongly convex thanks to the proximal term. We thus suggest the use of the Dykstra-like algorithm by [START_REF] Bauschke | A Dykstra-like algorithm for two monotone operators[END_REF], whose iterations are recalled in the Appendix B. This method presents the advantage of fast convergence rate, ease of implementation and no parameter tuning.

Minimization with respect to P: Let A ∈ R Nx×Nx and P ∈ S ++ Nx . The update of Eq. ( 27) solves a minimization problem with generic form minimize P∈S Nx C 2 (P)

θ P K 2 tr (PΠ) - θ P K 2 log det P + θ P λ P P 1 + 1 2 P -P 2 F , (33)
where we denote

Π = Ψ -∆ A ⊤ -A∆ ⊤ + AΦ A ⊤ . ( 34
)
Here we have used the definition of the proximity operator (25) and the majorant expression in (30) (ignoring the constant multiplicative term in the logarithm). Matrix Π ∈ S Nx in (34) is constant with respect to variable P. Remarkably, (33) reads as a regularized form of the famous GLASSO problem [START_REF] Friedman | Sparse inverse covariance estimation with the graphical LASSO[END_REF], and gets actually equivalent to it when θ P → ∞. Matrix Π plays the role of the empirical covariance matrix in GLASSO, and 2λ P K acts as the weight on the ℓ 1 term. The proximal term works as a Tikhonovlike regularizer, ensuring the strong convexity of the problem and thus the unicity of its minimizer. Moreover, by the definition of the log-determinant in Eq. ( 1), the solution of (33) belongs to P ∈ S ++ Nx , i.e., the precision matrix is invertible and thus a covariance matrix can be deduced from it by inversion. Standard GLASSO solvers can be easily modified to solve (33). We present in the Appendix B the complete derivations, when applying the Dykstralike algorithm from [START_REF] Bauschke | A Dykstra-like algorithm for two monotone operators[END_REF], which presents the advantage of fast convergence and ease of implementation.

Convergence analysis

We now present our convergence proof for the proposed DGLASSO algorithm presented in Algorithm 3. Our analysis is made assuming that the inner steps ( 26) and ( 27) are solved in an exact manner. The extension of the result to the case of an inexact resolution of the subproblems is discussed at the end of the section.

Descent property

Lemma 1 Assuming exact resolution of (26) and (27), the sequence {A (i) , P (i) } i∈N produced by DGLASSO algorithm satisfies

(∀i ∈ N) L(A (i+1) , P (i+1) ) ≤ L(A (i) , P (i) ). ( 35
)
Proof Let i ∈ N. First, let us Theorem 1 at ( A, P) = (A (i+1) , P (i) ) (Inequality (a)), and the definition of

P (i+1) in (27) (Inequality (b)), L(A (i+1) , P (i+1) ) (a) ≤ Q(A (i+1) , P (i+1) ; A (i+1) , P (i) ) + L 1:K (A (i+1) , P (i) ) -Q(A (i+1) , P (i) ; A (i+1) , P (i) ) + λ A A (i+1) 1 + λ P P (i+1) 1 + 1 2θ A A (i+1) -A (i+1) 2 F + 1 2θ P P (i+1) -P (i) 2 F ( 36
) (b) ≤ Q(A (i+1) , P (i) ; A (i+1) , P (i) ) + L 1:K (A (i+1) , P (i) ) -Q(A (i+1) , P (i) ; A (i+1) , P (i) ) + λ A A (i+1) 1 + λ P P (i) 1 + 1 2θ A A (i+1) -A (i+1) 2 F + 1 2θ P P (i) -P (i) 2 F . (37) Inequality (37) simplifies into L(A (i+1) , P (i+1) ) ≤ L 1:K (A (i+1) , P (i) ) + λ A A (i+1) 1 + λ P P (i) 1 = L(A (i+1) , P (i) ). ( 38
)
Applying Theorem 1 now at ( A, P) = (A (i) , P (i) ) (Inequality (a)), and the definition of A (i+1) in (26) (Inequality (b)), leads to

L(A (i+1) , P (i) ) (a)
≤ Q(A (i+1) , P (i) ; A (i) , P (i) ) + L 1:K (A (i) , P (i) ) -Q(A (i) , P (i) ; A (i) , P (i) )

+ λ A A (i+1) 1 + λ P P (i) 1 + 1 2θ A A (i+1) -A (i) 2 F + 1 2θ P P (i) -P (i) 2 F ( 39
) (b) ≤ Q(A (i) , P (i) ; A (i) , P (i) ) + L 1:K (A (i) , P (i) ) -Q(A (i) , P (i) ; A (i) , P (i) ) + λ A A (i) 1 + λ P P (i) 1 + 1 2θ A A (i) -A (i) 2 F + 1 2θ P P (i) -P (i) 2 F , (40)
which simplifies into

L(A (i+1) , P (i) ) ≤ L 1:K (A (i) , P (i) ) + λ A A (i) 1 + λ P P (i) 1 = L(A (i) , P (i) ), (41)
and concludes the proof.

If the cost function L is lower bounded (e.g., if it is coercive), Lemma 1 implies the convergence of sequence L(A (i) , P (i) ) i∈N to a finite value and, as such, the existence of cluster points in {A (i) , P (i) } i∈N . This is however a rather weak convergence result and we propose hereafter a thorough analysis relying on recent tools of nonlinear analysis [START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems:an approach based on the kurdyka-lojasiewicz inequality[END_REF][START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF] combined with the works [START_REF] Tien | An inertial block majorization minimization framework for nonsmooth nonconvex optimization[END_REF][START_REF] Hien | Inertial block proximal method for non-convex nonsmooth optimization[END_REF] on the convergence of block alternating MM schemes.

Convergence guarantees

Theorem 3 Consider the sequence {A (i) , P (i) } i∈N generated by DGLASSO, assuming exact resolution of both inner steps (26) and (27). If the sequence {A (i) , P (i) } i∈N is bounded, then {A (i) , P (i) } i∈N converges to a critical point of L.

Proof The convergence analysis relies in proving that the exact form of DGLASSO algorithm is a special instance of TITAN algorithm from [START_REF] Tien | An inertial block majorization minimization framework for nonsmooth nonconvex optimization[END_REF], and as such, inherits (Tien et al., 2020, Theorem 7) and (Tien et al., 2020, Theorem 8), under our assumptions.

• Let us introduce the following notations.

(∀A ∈ R Nx×Nx )(∀P ∈ S Nx ) f (A, P) = 1 2 K k=1 (x k -Ax k-1 ) ⊤ P (x k -Ax k-1 ) + log p(x 0:K |y 1:K , A, P) -log p(x 0 ) - K k=1 log p(y k |x k ), ( 42
) (∀A ∈ R Nx×Nx ) g 1 (A) = λ A A 1 , (43)
(∀P ∈ S Nx ) g 2 (P) = - K 2 log det(2πP) + λ P P 1 , (44) so that
(∀A ∈ R Nx×Nx )(∀P ∈ S Nx ) L(A, P) = f (A, P) + g 1 (A) + g 2 (P), (45)
with f lower semi-continuous function, g 1 ∈ Γ 0 (R Nx×Nx ) and g 2 ∈ Γ 0 (S Nx ). Moreover, let us denote

(∀A ∈ R Nx×Nx )(∀P ∈ S Nx ) J (A, P) = K 2 tr P(Ψ -∆A ⊤ -A∆ ⊤ + AΦA ⊤ ) , (46)
and, for every A ∈ R Nx×Nx and P ∈ S ++ Nx ,

(∀A ∈ R Nx×Nx ) u 1 (A; A, P) = J (A, P) + 1 2θ A A -A 2 F + f ( A, P) -J ( A, P), (47)
(∀P ∈ S Nx ) u 2 (P; A, P) = J ( A, P) + 1 2θ P P -P 2 F + f ( A, P) -J ( A, P).

By Theorem 1, the following majorization properties hold for every A ∈ R Nx×Nx and

P ∈ S ++ Nx , (∀A ∈ R Nx×Nx ) u 1 (A; A, P) ≥ f (A, P), (49)
(∀P ∈ S Nx ) u 2 (P; A, P) ≥ f ( A, P).

(50)

and we have the tangency condition, for every A ∈ R Nx×Nx and P ∈ S ++ Nx ,

u 1 ( A; A, P) = f ( A, P), (51)
u 2 ( P; A, P) = f ( A, P).

(52)

Then, straightforward computations allow to rewrite the iterates of Algorithm 3 as follows: i) , P (i) ) + g 1 (A),

(∀i ∈ N)      A (i+1) = argmin A∈R Nx×Nx u 1 (A; A (
P (i+1) = argmin P∈S Nx u 2 (P; A (i+1) , P (i) ) + g 2 (P), (53)
which identifies with the iterative scheme TITAN from [START_REF] Tien | An inertial block majorization minimization framework for nonsmooth nonconvex optimization[END_REF], in the case of two blocks and setting the extrapolation step to zero. The rest of the proof amounts to check the fulfillment of the assumptions required for (Tien et al., 2020, Theorem 7) and (Tien et al., 2020, Theorem 8).

• Let us denote, for every A ∈ R Nx×Nx and P ∈ S ++ Nx ,

(∀A ∈ R Nx×Nx ) ũ1 (A; A, P) = u 1 (A; A, P) + g 1 (A), (54)
(∀P ∈ S Nx ) ũ2 (P; A, P) = u 2 (P; A, P) + g 2 (P). ( 55)

Functions u 1 and u 2 are quadratic and strongly convex with respective strong convexity constants θ -1 A and θ -1 P . Since both g 1 and g 2 are convex, functions ũ1 and ũ2 are also strongly convex, with respective strong convexity constants θ -1

A and θ -1 P . Let i ∈ N. According to the optimality conditions of both equations in (53), there exists

T (i+1) 1 ∈ ∂ ũ1 (A (i+1) ; A (i) , P (i) ) ∈ R Nx×Nx and T (i+1) 2 ∈ ∂ ũ2 (P (i+1) ; A (i+1) , P (i) ) ∈ S Nx such that    tr T (i+1) 1 (A (i) -A (i+1) ) ≥ 0, tr T (i+1) 2 (P (i) -P (i+1) ) ≥ 0. ( 56
)
Moreover, by strong convexity of both ũ1 and ũ2 ,    ũ1 (A (i) ; A (i) , P (i) ) ≥ ũ1 (A (i+1) ; A (i) , P (i) ) + tr

T (i+1) 1 (A (i) -A (i+1) ) + 1 2θ A A (i+1) -A (i+1) 2
F , ũ2 (P (i) ; A (i+1) , P (i) ) ≥ ũ2 (P (i+1) ; A (i+1) , P (i) ) + tr T (i+1) 2 (P (i) -P (i+1) ) + 1 2θ P P (i+1) -P (i+1) 2 F . (57) Hence, using (56), ũ1 (A (i) ; A (i) ,

P (i) ) ≥ ũ1 (A (i+1) ; A (i) , P (i) ) + 1 2θ A A (i+1) -A (i+1) 2 F , ũ2 (P (i) ; A (i+1) , P (i) ) ≥ ũ2 (P (i+1) ; A (i+1) , P (i) ) + 1 2θ P P (i+1) -P (i+1) 2 F , (58)
and, using ( 49)-( 50)-( 51)-(52),

L(A (i) , P (i) ) ≥ L(A (i+1) , P (i) ) + 1 2θ A A (i+1) -A (i+1) 2 F , L(A (i+1) , P (i) ) ≥ L(A (i+1) , P (i+1) ) + 1 2θ P P (i+1) -P (i+1) 2 F . (59)
It means that the so-called NDSP (nearly sufficiently decreasing property) condition from [START_REF] Tien | An inertial block majorization minimization framework for nonsmooth nonconvex optimization[END_REF] holds with, for every i ∈ N, (γ

(i) 1 , γ (i) 2 , η (i) 1 , η (i)
2 ) ≡ (0, 0, θ -1 A , θ -1 P ). Remark that our proof for (59) constitutes an alternative proof for Lemma 1.

• According to [START_REF] Gupta | Computational aspects of maximum likelihood estimation and reduction in sensitivity function calculations[END_REF], function f is continuously differentiable. As g 1 ∈ Γ 0 (R Nx×Nx ) and g 2 ∈ Γ 0 (S Nx ), we have for every A ∈ R Nx×Nx and P ∈ S Nx ,

∂ A (f (A, P) + g 1 (A)) = ∇ A f (A, P) + ∂g 1 (A), ∂ P (f (A, P) + g 2 (P)) = ∇ P f (A, P) + ∂g 2 (P), (60)
so that [START_REF] Tien | An inertial block majorization minimization framework for nonsmooth nonconvex optimization[END_REF], Assumption 3(i)) holds. Moreover, by construction of the majorizing function J , it is also continuously differentiable and we have the coincidence of the gradient at the majorization point, namely for every A ∈ R Nx×Nx and P ∈ S ++ Nx ,

∇ A f ( A, P) = ∇ũ 1 ( A; A, P), (61)
∇ P f ( A, P) = ∇ũ 2 ( P; A, P).

Thus, according to (Tien et al., 2020, Lem.2), [START_REF] Tien | An inertial block majorization minimization framework for nonsmooth nonconvex optimization[END_REF], Assumption 2) is verified. Moreover, as the restriction of J to both of its variables is quadratic, the partial gradients ∇u 1 and ∇u 2 are linear and thus Lipschitz continuous on any bounded subset of R Nx×Nx and S Nx , respectively, which yields the fulfillment of (Tien et al., 2020, Assumption 3).

• Since NSDP and (Tien et al., 2020, Assumption 2) hold, the result (Tien et al., 2020, Theorem 6) allows to show that, if the sequence {A (i) , P (i) } i∈N is bounded, then every limit point (A * , P * ) of it is a critical point of L. Moreover, (Tien et al., 2020, Assumption 3) is satisfied and (Tien et al., 2020, Cond.1) and (Tien et al., 2020, Cond.4) are trivially met in our case. We can also show that the loss function L satisfies the Kurdyka-Lojasiewicz property [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF]) (using a similar proof than (Chouzenoux et al., 2019, Lemma 2)). Thus, (Tien et al., 2020, Theorem 7) holds which concludes our proof.

Discussion

The proof of Theorem 3 is grounded on the recent works by [START_REF] Tien | An inertial block majorization minimization framework for nonsmooth nonconvex optimization[END_REF] and [START_REF] Hien | Inertial block proximal method for non-convex nonsmooth optimization[END_REF], generalizing the works by [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF] and [START_REF] Chouzenoux | A block coordinate variable metric forwardbackward algorithm[END_REF], for establishing the convergence of block alternating MM schemes under the Kurdyka-Lojasiewicz (KL) inequality assumption [START_REF] Lojasiewicz | Unepropriété topologique des sous-ensembles analytiques réels[END_REF]. The latter assumption is a powerful tool of nonlinear functional analysis that has been popularized in the seminal paper by [START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems:an approach based on the kurdyka-lojasiewicz inequality[END_REF]. In their paper, the authors show how to prove the convergence of the iterates generated by a large family of minimization schemes, under the sole requirement that the function to minimize satisfies the KL inequality. The latter requirement is very mild, as it holds for a large class of functions, non necessarily convex, as soon as they can be embedded within an o-minimal structure. Semi-algebraic and analytical functions are example of such functions. In our analysis, for the sake of conciseness, we skipped the proof showing that L satisfies KL as it is identical to the one of (Chouzenoux et al., 2019, Lem.2) .

Theorem 3 assumes that both inner steps ( 26) and ( 27) are solved in an exact manner. Due to strong convexity, each problem has a unique solution, which ensures the well posedness of the study. However, in both case, the solution does not take a closed form (except when λ A or λ P equals zero, in such case we retrieve the MLEM updates from (Sarkka, 2013, Chap. 12)). Iterative inner solvers are thus necessary, and we proposed some efficient ones in Section 3.5. The extension of our convergence study to the case of inexact resolution of ( 26) and ( 27) is not straightforward, up to our knowledge. Convergence of inexact proximal schemes for KL losses has been studied in various works, such as [START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms,forward-backward splitting[END_REF][START_REF] Cherni | SPOQ lp-over-lq regularization for sparse signal recovery applied to mass spectrometry[END_REF][START_REF] Chouzenoux | Variable metric forward-backward algorithm for minimizing the sum of a differentiable function and a convex function[END_REF] to name a few. But we are not aware of any study covering the block alternating MM scheme considered here, and thus left the convergence study of the inexact implementation of DGLASSO as future work. In practice, we impose a rather demanding condition on the stopping rule for the inner solvers of ( 26) and ( 27) and did not observe any numerical instabilities of the proposed algorithm.

Experiments

We perform a thorough evaluation study on various controlled scenarios where the ground truth matrices denoted (A * , P * ) (as well as Q * = (P * ) -1 ) are predefined, and the time series {y k , x k } K k=1 are build directly from the LG-SSM model ( 2)-( 3) using such matrices. The goal is then to provide estimates ( A, P, Q) of (A * , P * , Q * ), given the observation of {y k } K k=1 . Except otherwise stated, all the compared methods have moreover access to a perfect knowledge of (R k , H k , µ 0 , Σ 0 ), in top of accessing the time series {y k } K k=1 . The hidden state {x k } K k=1 is, by definition, not assumed to be known, and is only used to compute quality metrics on test sets. We first work on synthetic data in Section 5.1, where the structure, the sparsity level, the conditioning, of the sought matrices (A * , P * ) is controlled. This allows us to evaluate our method on multiple cases, to discuss its parameter tuning, and to compare it to benchmarks in terms of inference quality and complexity. We then address in Section 5.2, a set of realistic problems of graph inference arising in weather variability analysis, using four datasets built upon the Neurips data challenge [START_REF] Runge | The causality for climate competition[END_REF]. This second set of experiments aims at evaluating the ability of DGLASSO to model and estimate a large class of graph structures (here, 200 different graphs per dataset), in comparison with other state-of-the-art graph detection methods.

All codes are run on a Desktop Dell Latitude computer, with 11th Gen Intel(R) Core(TM) i7-1185G7 at 3.00GHz, equipped with 32Go Ram, using Matlab R2021a software. The code is publicly available at https://pages.saclay.inria.fr/emilie.chouzenoux/log/DGLASSO_toolbox.zip, for reproducibility purpose.

Synthetic data

Datasets

We set K = 10 3 , R k = σ 2 R Id Ny for every k ∈ {1, . . . , K}, µ 0 ∈ R Nx is a vector of ones, Σ 0 = σ 2 0 Id Nx with (σ R , σ 0 ) = (10 -1 , 10 -4 ). Matrix H k is set to identity matrix for every k ∈ {1, . . . , K}, so that N x = N y . This setting models a one-to-one correspondence between states and observations. Identifiability issues are hence avoided.

We set N x = N y = 9, and we rely on block-diagonal matrices (A * , P * ), made of 3 blocks with dimensions 3 × 3. The diagonal blocks of A * are randomly set as matrices of auto-regressive processes of order one, AR(1), rescaled to have a spectral norm equals 0.9 to ensure the process stability. The diagonal blocks of P * are randomly set following the ground truth matrices (A * , P * , Q * ) and the estimated ( A, P, Q) (when available). For instance,

RMSE(A

* , A) = A * -A 2 F A * 2 F . (63)
RMSE(P * , P) and RMSE(Q * , Q) are defined in a similar fashion. We also compute areaunder-curve (AUC) and F1 score comparing the non-zero entries (that is, the graph edges positions) of the sparse matrices (A * , P * ) and their estimates ( A, P). A threshold value of 10 -10 (in absolute value) is used for the detection hypothesis. We furthermore evaluate the ability of the estimated model parameters to actually describe and infer the time series (both observed and hidden states). To that end, we build test time series (x test , y test ), not seen by the algorithms, constructed by running the ground truth LG-SSM (i.e., with ground truth matrix parameters (A * , P * , Q * )). We then run KF and RTS algorithms 1 and 2, respectively, using either the ground truth matrices (A * , P * , Q * ) or their estimations ( A, P, Q), to build, for every k ∈ {1, . . . , K}, the predictive distribution means

(µ * k , ν * k , µ s * k ) and ( µ k , ν k , µ s k )
, respectively. This allows in particular to compute the cumulative normalized mean squared error (cNMSE) between the predictive distribution means using either the ground truth model matrices or the estimated ones. Namely, we calculate

cNMSE(ν * , ν) = K k=1 ν * k -ν 2 K k=1 ν * k 2 , ( 64
)
as well as cNMSE(µ * , µ), and cNMSE(µ s * , µ s ). Finally, we evaluate the negative logarithm of the marginal likelihood L 1:K ( A, P) as defined in (21), on the test time series.

Tuning the regularization parameters

Let us first discuss the setting of the DGLASSO hyper-parameters (λ A , λ P ), accounting for the sparsity priors on A and P, respectively. To that aim, we ran DGLASSO for 100 values of hyperparameters (λ A , λ P ), regularly spaced on a log-scale grid between 1 and 10 2 , and repeated the experience on 50 randomly generated time series, for dataset A. We display in Figure 4 the averaged values, over the random runs, of several quantitative metrics, as a function of hyperparameters (λ A , λ P ). We also report in the caption the averaged RMSE scores obtained when running DGLASSO with (λ A , λ P ) = (0, 0) (i.e., MLEM result). As it can be observed, both F1/RMSE for the estimation of the transition matrix A are mostly governed by the value of λ A , while the quality scores for the state noise precision matrix P are mostly influenced by λ P . Note that the RMSE scores on Q, not shown here, follow similar evolution than RMSE on P. Just inspecting F1/RMSE appears not informative enough to set parameters (λ A , λ P ), as each metric and each parameter seems to push towards a different goal. The maps of cNMSE(ν * , ν) and of the marginal likelihood logloss L 1:K show very similar behavior. Note that the later is a practically interesting metric, because it does not require the knowledge of the ground truth matrices. On this example, it however appears not discriminative enough. Typically, it stays almost constant for a (too) wide value range of (λ A , λ P ), which confirms again the ill-posedness of the minimization problem. The maps for cNMSE(µ * , µ) and cNMSE(µ s * , µ s ) are very similar. Interestingly, the minimization of both these quantities, related to the state mean distributions, seems to provide a meaningful value range for the regularization parameters, narrowing down around values that achieve an optimal compromise between (i) good estimation of the sought matrices, (ii) good estimation of the sparse matrices support, and (iii) good predictive behavior for time series inference by KF/RTS techniques. Same conclusions were reached for the other three datasets. Note additionally that the DGLASSO results outperform for a wide range of (λ A , λ P ) those obtained with MLEM, confirming the validity of the proposed sparsity prior. This will be more deeply examined in an upcoming section. In our forthcoming experiments, except otherwise stated, we fixed (λ A , λ P ) through a rough grid search to minimize cNMSE(µ * , µ) averaged on few (typically 5) runs. On real datasets, we advocate the use of the marginal likelihood on test trials, as a control variate.

Influence of conditioning number

We report in Table 1 the results, averaged over 50 realizations of the time series generation. We do not provide the scores for the estimation of A (resp. (P, Q)) for GLASSO/rGLASSO (resp. GRAPHEM), as those methods are not designed for such task. The general trend in the results is a slight decrease of the quality of estimation for all methods, when c Q increases (i.e., from dataset A to D). This is expected, as an ill-conditioned matrix Q complicates the mathematical structure of the likelihood term.

Regarding the estimation of A, GRAPHEM method presents the best results for the three considered metrics. DGLASSO is second best, while MLEM follows. As already stated, GLASSO/rGLASSO do not estimate A (i.e., they assume this matrix to be zero, that is the process is i.i.d. without any Markov structure). Regarding the estimation of (P, Q), DGLASSO is outperforming the benchmarks in terms of RMSE score. The second best, in terms of RMSE is MLEM. GLASSO and rGLASSO gets very bad RMSE scores, probably because of the model mismatch induced by assuming A to be zero. The edge detection performance of DGLASSO are excellent, except in few cases where rGLASSO gets slightly better results. MLEM gets poorer results than DGLASSO in all metrics. In particular it exhibits a bad F1 score, as it does not include sparse prior, and thus does not succeed to eliminate any spurious edges.

Regarding the distribution mean calculation, it is remarkable that DGLASSO outperforms in all examples the benchmarks, which shows that the proposed formulation allows to provide model parameters that are best suited for inference using KF/RTS at test phase. The marginal likelihood log-loss is also minimized with the quantities provided by our method. This could appear as counter-intuitive, as the method is not specifically designed to minimize this loss (while MLEM is). The advantage of DGLASSO is that it accounts for prior knowledge, making the inferred model more robust, and less prone to overfitting, which is translated into better behavior on an independent test time series.

We display in Figure 5 box plots assessing the variability to the RMSE and F1 scores, for both MLEM and DGLASSO methods, on 50 runs on dataset A and dataset D. The RMSE values are in most runs lower (i.e., better) for the proposed method. Both methods are quite stable with respect to the time series generation, as the plots show few outliers. For dataset D, corresponding to a more challenging case with an ill-conditioned matrix P * , and dataset D (bottom). DGLASSO outperforms MLEM on all runs, in terms of accuracy metrics (the higher, the better), while its RMSE scores are better (i.e., lower) for DGLASSO for most runs. Dataset D is more challenging, in terms of inference, thus yielding to more spread results for both methods.

the results are slightly more spread, especially for the metrics related to the recovery quality of this matrix. Remark that the F1 scores of MLEM are constant, equal to 0.5. As already pointing out, this is expected as MLEM is not designed to perform an edge detection task. We refer the reader to Appendix C.3 for extra experiments on the synthetic datasets, assessing the performance of DGLASSO when the sparsity level of A * increases.

Complexity analysis

We now examine the time complexity of the method, as well as of MLEM, when processing dataset A, for various values of the time series length K, namely K ∈ {100, 200, 500, 1000, 2000, 5000} (we recall that our previous experiments were all done with K = 1000). We display in Figure 6(left) the computing time in seconds for computing ( A, P), for DGLASSO and MLEM, averaged over 50 realizations. We also display, in Figure 6(middle) for the same experiment, the RMSE between A * and A, and in Figure 6(right) the metric cNMSE(µ * , µ). One can notice that our method is slightly more demanding in terms of computational time, but it scales similarly than its non regularized version MLEM. Moreover, as already observed in our previous experiments, DGLASSO outperforms MLEM on both metrics shown here, in all tested values of K. As expected, the results are better for higher values of K, at Table 1: Results for the four considered datasets A to D, with an increasing conditioning number of P * equals to log 10 (c) ∈ {0.1, 0.2, 0.1, 1}, respectively. We evaluate the methods in terms of estimation quality for (A, P, Q), using either RMSE as defined in (63), and edge detection scores (AUC, F1), as well as in terms of inference quality using cNMSE metrics defined in (64). the price of an increased computational time. Interestingly, the regularization still yields improved results for very large K.

Weather data

Experimental settings

We now evaluate our method on realistic graph datasets arising from causal discovery studies in the field of weather variability tracking. Specifically, we consider two sets of 200 sparse matrices A * ∈ R Nx , with N x = 5 or 10 respectively, representing the ground truth causal graphs used to produce WEATH datasets in the Neurips 2019 data challenge [START_REF] Runge | The causality for climate competition[END_REF] 5 . For each A * , we create times series (x k , y k ) K k=1 using ( 2)-(3), with K = 10 3 , H * = Id (i.e., N y = N x ), and (σ R , σ 0 ) = (10 -1 , 10 -4 ). We set Q * as a block diagonal matrix of J blocks with dimensions B j 1≤j≤J , with J j=1 B j = N x . Here, we consider two settings for the conditioning of Q * , namely one with a condition number close to one (i.e., Q * is close to identity matrix) and another with high condition number. The main characteristics of the datasets and their names are summarized in Table 2.

We evaluate our results in terms on quality assessment metrics of the estimated A when compared to its ground truth A * , as this is the quantity of interest for these datasets. We compute RMSE( A, A * ), as well as the precision, recall, specificity, accuracy, and F1 score for detecting the non-zero entries of the transition matrix (that is, the graph edges positions). A threshold value of 10 -10 on the absolute entries is used for the detection hypothesis.

As for comparison, we also provide the results obtained with MLEM and GRAPHEM approaches, as in our previous experiments. The same stopping criterion than in our previous set of experiments is used. The hyperparameters are finetuned using the strategy depicted on the previous section. We used rough grids for the finetuning and keep the parameters fixed for all the dataset, to avoid any overfitting behavior. In addition to these EM-based methods, we provide comparisons with two Granger-causality approaches for graphical modeling, namely pairwise Granger Causality (PGC) and conditional Granger Causality (CGC). Those approaches are based on sparse vector autoregressive (VAR) models. We allow the order of the VAR process to be up to p = 11, CGC is run with a maximum distance of 5 causal links, and in each experiment we display the best results (in F1 score) for the statistical tests performed with the significance level α ∈ {0.01, 0.05, 0.1} (see more details in [START_REF] Luengo | Hierarchical algorithms for causality retrieval in atrial fibrillation intracavitary electrograms[END_REF]). As PGC and CGC do not provide a weighted graph estimation, no RMSE score is computed in those cases.

Results

We summarize our results in Table 3. We also report the averaged inference time for all methods. Remarkably, the proposed method outperforms all benchmarks in all the considered metrics, related to the inference of the graph weights (RMSE metric) and the edge detection task (binary classification metrics). As expected, the result quality tends to slightly degrade when a more complicated structure is assumed for Q * (see, for instance, WeathN5a versus WeathN5b), and when the dataset size increases (see, for instance, WeathN5a versus WeathN10a). MLEM has very poor edge detection metrics, since it does not exploit any sparsity prior on the sought matrix. GRAPHEM has better behavior, but, in these datasets, it stays way behind the quality of DGLASSO, by all metrics (which was not the case for the synthetic data). This shows that our method really makes the change when dealing with complex graphs and correlated noise in the observations. In terms of computational time, the proposed method has similar complexity than the two other EM-based approaches. The results of CGC and PGC are satisfactory in the first two examples although the binary metrics are far from the performance of DGLASSO. PGC still gives meaningful results in the last two examples, but CGC has a low performance due to a high number of false negatives (i.e., many existing links are not discovered). We also note that while PGC and CGC have significantly lower running times for the case with N y = 5, the computational cost for N y = 10 is closer to DGLASSO while the performance gap is higher. Thus, in this examples we show that as the dimension grows, the computational cost of all methods is comparable while the proposed method has clearly a better performance. We also display some examples of inferred graphs in Figures 7 and8. We can observe that DGLASSO is able to capture in a very accurate manner the structure of the graphs, despite the wide variability of the dataset. MLEM and GRAPHEM capture the main edges, but their graphs are perturbed by several spurious edges, which shows how important it is to adopt a joint graph modeling, with sparsity priors on each.

Conclusion

This paper proposes a joint graphical modeling approach that incorporates a graphical perspective on the dynamics of the hidden state of a linear Gaussian state-space model. From top to bottom: Original graph representation of A * , and of its estimation A, using DGLASSO, MLEM, GRAPHEM, respectively.

In particular, we propose a joint approach that considers a sparse undirected graph as a prior on the precision matrix of the hidden state noise and a sparse directed graph for the transition matrix that models the state dynamics. By bridging the gap between the static graphical Lasso model and the dynamic state-space model, we provide a novel comprehensive framework for interpretable inference with time-series data. The presented inference method, based on an efficient block alternating majorization-minimization algorithm, enables simultaneous estimation of both graphs and the construction of the filtering/smoothing From top to bottom: Original graph representation of A * , and of its estimation A, using DGLASSO, MLEM, GRAPHEM, respectively. distribution for the time series. The established convergence of our algorithm, departing from recent nonlinear analysis tools, enhances the reliability and practicality of our approach. Through extensive experimental validation on synthetic and real weather variability data, we have demonstrated the effectiveness and potential of our proposed model and inference algorithm. The results showcase its ability to uncover meaningful insights from time-series data, contributing not only to better forecasting performance but also to a better understanding of complex phenomena in various scientific and engineering domains. Fu-ture research can further explore and extend the presented framework to tackle even more challenging state-space models, and more complex graphical structures.

A. Proof of Theorem 1

For any initial state x 0:K , with non zero probability, the neg-log-likelihood (A, P) → (L 1:K (A, P) =log p(y 1:K |A, P)) reads, according to Bayes' rule, (∀A ∈ R Nx×Nx )(∀P ∈ S Nx )

L 1:K (A, P) =log p(x 0:K , y 1:K |A, P) + log p(x 0:K |y 1:K , A, P).

Moreover, using again Eq. ( 2) and the statistical model of the state noise,

(∀A ∈ R Nx×Nx )(∀P ∈ S Nx ) (67) K k=1 log p(x k |x k-1 , A, P) = - 1 2 K k=1 (x k -Ax k-1 ) ⊤ P (x k -Ax k-1 ) -log det(2πP) , (68)
= - 1 2 K k=1 (x k -Ax k-1 ) ⊤ P (x k -Ax k-1 ) + K 2 log det(2πP), (69)
which concludes the first part of the proof.

Let us now consider some A ∈ R Nx×Nx and P ∈ S ++ Nx . We start by recalling some known results arising from the EM methodology [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF][START_REF] Wu | On the convergence properties of the EM algorithm[END_REF], that we specify here for our context for the sake of clarity. First, we notice that (∀A ∈ R Nx×Nx )(∀P ∈ S Nx ) L 1:K (A, P) =log p(y 1:K |A, P)p(x 0:K |y 1:K , A, P)dx 0:K , (70) since log p(y 1:K |A, P) is constant with respect to the integration variable, and the distribution p(x 0:K |y 1:K , A, P) integrates to one. Then, according to (65) and ( 70), the expectation of the neg-log-likelihood multiplied by p(x 0:K |y 1:K , A, P) over all possible values of the unknown state reads:

(∀A ∈ R Nx×Nx )(∀P ∈ S Nx ) L 1:K (A, P) = q(A,P; A, P) p(x 0:K |y 1:K , A, P) log p(x 0:K , y 1:K |A, P)dx 0:K + h(A,P; A, P) p(x 0:K |y 1:K , A, P) log p(x 0:K |y 1:K , A, P)dx 0:K . (71)

In particular, for (A, P) = ( A, P), L 1:K ( A, P) = q( A, P; A, P) + h( A, P; A, P) so that (∀A ∈ R Nx×Nx )(∀P ∈ S Nx ) L 1:K (A, P) -L 1:K ( A, P) = q(A, P; A, P)q( A, P; A, P) + h(A, P; A, P)h( A, P; A, P). (72) Using Gibbs's inequality, h(A, P; A, P) ≤ h( A, P; A, P), with equality at (A, P) = ( A, P). Thus, using (72), that is (∀A ∈ R Nx×Nx )(∀P ∈ S Nx ) L 1:K (A, P) ≤ q(A, P; A, P) + L 1:K ( A, P)q( A, P; A, P), (73) where equality holds at (A, P) = ( A, P). Notice that, for any function reading (∀A ∈ R Nx×Nx )(∀P ∈ S Nx ) Q(A, P; A, P) = q(A, P; A, P) + ct /A,P ,

we obviously still have

(∀A ∈ R Nx×Nx )(∀P ∈ S Nx ) L 1:K (A, P) ≤ Q(A, P; A, P) + L 1:K ( A, P) -Q( A, P; A, P), (75)
where equality hereagain holds at (A, P) = ( A, P). Using ( 75), ( 20), ( 22) and noticing that, for every (θ A , θ P ) > 0,

(∀A ∈ R Nx×Nx )(∀P ∈ S Nx ) 1 2θ A A -A 2 F ≥ 0, 1 2θ P P -P 2 F ≥ 0, (76)
with equality holding at (A, P) = ( A, P), we deduce the desired majorizing property (∀A ∈ R Nx×Nx )(∀P ∈ S Nx ) L(A, P) ≤ Q(A, P; A, P) + L 1:K ( A, P) -Q( A, P; A, P)

+ λ A A 1 + λ P P 1 + 1 2θ A A -A 2 F + 1 2θ P P -P 2 F , (77)
with equality holding at (A, P) = ( A, P). The remaining of the proof amounts to expliciting the expression for (A, P) → Q(A, P; A, P) satisfying ( 74) with function q defined as in ( 71):

(∀A ∈ R Nx×Nx )(∀P ∈ S Nx ) q(A, P; A, P) =p(x 0:K |y 1:K , A, P) log p(x 0:K , y 1:K |A, P)dx 0:K . (78)

Following (Sarkka, 2013, Theorem 12.4) (see also an alternative proof in [START_REF] Elvira | Graphical inference in linear-Gaussian state-space models[END_REF], Sec.III-B)), ( 71

Using (Bauschke and Combettes, 2017, Chap. 24) (see also [START_REF] Benfenati | Proximal approaches for matrix optimization problems: Application to robust precision matrix estimation[END_REF]), for every α > 0,

Z = UDiag 1 2 (ω(n) + ω(n) 2 + 4γ) 1≤n≤Nx U ⊤ (93)
where ω = (ω(n)) 1≤n≤Nx gathers the eigenvalues of W -γΠ ∈ S Nx and U ∈ R Nx×Nx is an orthogonal matrix such that W -γΠ = UDiag(ω)U ⊤ .

(94)

C. Additional experiments

We present here additional experimental results completing Section 5.1.

C.1 Maximum likelihood calculation

We compare MLEM approach from [START_REF] Shumway | An approach to time series smoothing and forecasting using the EM algorithm[END_REF] with the sensitivity equation approach from [START_REF] Segal | A new method for evaluating the log-likelihood gradient, the hessian, and the fisher information matrix for linear dynamic systems[END_REF], for computing maximum likelihood estimates of the sought matrices. For the latter, we use the gradient expressions from [START_REF] Nagakura | Computing exact score vectors for linear Gaussian state space models[END_REF], and performed the minimization using a Quasi-Newton routine from the Optimization Toolbox from Matlab, and used the notation MLQN. We ran both MLEM and MLQN methods on the four synthetic datasets A to D. Figure 9 displays the evolution of the loss function along iterations for both MLEM and MLQN algorithms. One can notice that both reach a similar value at convergence. Moreover, MLEM requires very few iterations (typically, less than 10) to stabilize, while MLQN displays a slower convergence profile despite the quasi-Newton (i.e., second-order) acceleration strategy. This might be due to the ill conditioning of the minimization problem. In terms of time complexity, both methods are comparable, requiring a full pass on the time series, through a Kalman update, at each iteration. Given these results, we decided to keep only MLEM as a comparative benchmark for our experiments. (i) , Q (i) ) along iterations using either MLEM or MLQN approach.

Figure 10: Box plots for the RMSE and F1 scores for retrieving (A, P, Q) matrices, when running DGLASSO on one single LG-SSM time series using dataset A, and 50 random initializations (A (0) , P (0) ). Noticeably, low variability is observed for all metrics.

C.2 Robustness to initialization

DGLASSO algorithm amounts to minimizing a non-convex loss function. As such, its results might be sensitive to the initialization of the algorithm. To evaluate this aspect, we consider the computation of ( A, P) given the observation of a single time series generated by the ground truth LG-SSM, when using 50 different initializations of DGLASSO algorithm. To do so, we use the same initialization strategy as discussed above, now with (a, p) randomly selected as p ∼ U ([0, 1]) and a ∼ U ([0, 1]). Figure 10 displays the box plots for RMSE and F1 scores obtained for dataset A. One can notice that the box plots are very concentrated, showing a good robustness of the method to its initialization, with the wider spreading observed for the F1 score on A. Similar behavior was observed for the other three datasets.

C.3 Influence of sparsity level

We evaluate here the performance of DGLASSO, as well as the benchmarks, when varying the sparsity level of the ground truth matrices. To do so, we perform slight changes in the dataset, to vary the sparsity pattern of the matrices (i.e., the edge structure of the graphs). First, we modify the ground truth matrix A * by keeping s A ∈ {27, 15, 10, 5} entries of it, within the 27 block diagonal ones, to be non-zero, the others being set to zero. We then rescaled the matrix to keep a spectral norm equals 0.99. Matrix Q * is taken from dataset A. The results are reported in Table 4.

As we can observe, the performance of MLEM, in terms of RMSE and F1 score, drop dramatically when the sparsity level on A increases. This is expected as this approach does not promote any sparsity prior on matrix A. The best AUC are either obtained by MLEM, DGLASSO or rGLASSO (for P), depending on the test cases. GLASSO/rGLASSO metrics slightly improve when A * gets sparser, which is expected, as their assumption of a zero transition matrix gets more realistic. Hereagain, DGLASSO outperforms the benchmarks in most cases.

Figure 2 :

 2 Figure 2: Matrix P (a), its binary support (b), and the associated undirected graph (c) with edge thickness proportional to the absolute entries of P.

Figure 4 :

 4 Figure4: Evolution of RMSE, F1, cNMSE and loss scores on the estimation of A (left) and P (right) by DGLASSO, as a function of hyperparameters (λ A , λ P ), for dataset A (averaged on 10 runs). As a comparison, the averaged RMSE scores for (λ A , λ P ) = (0, 0) (i.e., MLEM) on this example were (0.077, 0.106) for (A, P), respectively.

Figure 5 :

 5 Figure 5: Box plots for the RMSE and F1 scores when running MLEM (left) and DGLASSO (right) on 50 randomly generated LG-SSM time series, for dataset A (top)and dataset D (bottom). DGLASSO outperforms MLEM on all runs, in terms of accuracy metrics (the higher, the better), while its RMSE scores are better (i.e., lower) for DGLASSO for most runs. Dataset D is more challenging, in terms of inference, thus yielding to more spread results for both methods.

Figure 6 :

 6 Figure 6: Evolution of the complexity time (left), RMSE(A * , A) (middle) and cNMSE(µ * , µ) (right) metrics, as a function of the time series length K, for experiments on dataset A averaged over 50 runs.

Figure 7 :

 7 Figure 7: Graph inference results on 4 examples extracted from the dataset WeathN5a.From top to bottom: Original graph representation of A * , and of its estimation A, using DGLASSO, MLEM, GRAPHEM, respectively.

Figure 8 :

 8 Figure 8: Graph inference results on 4 examples extracted from the dataset WeathN10a.From top to bottom: Original graph representation of A * , and of its estimation A, using DGLASSO, MLEM, GRAPHEM, respectively.

  x 0:K , y 1:K |A, P) = log p(x 0 ) + K k=1 log p(x k |x k-1 , A, P) + K k=1 log p(y k |x k ).

  )-(74) hold for(∀A ∈ R Nx×Nx )(∀P ∈ S Nx ) Q(A, P; A, P) = K 2 tr P(Ψ -∆A ⊤ -A∆ ⊤ + AΦA ⊤ ) -K 2 log det(2πP),(79)Since P ∈ S ++ Nx , and Ψ ∈ S Nx (by construction), we have equivalently,-γ∆γ P -1 ∆ ⊤ P + 2γ ZΨ + P -1 Z -P -1 W = 0.(89)Thus,Z = lyapunov P -1 , 2γΨ, γ(∆ + P -1 ∆ ⊤ P) + P -1 W , (90)where A = lyapunov(X, Y, Z) provides the solution to the Lyapunov equation XA + AY = Z.Proximity of log-determinant term. Let γ > 0 and W ∈ S ++ Nx . Then, by definition, Z = prox W→γ(-log det(W)+tr(WΠ))

Figure 9 :

 9 Figure 9: Evolution of L 1:K(A (i) , Q (i)) along iterations using either MLEM or MLQN approach.

Table 2 :

 2 5. https://causeme.uv.es/static/datasets/TestWEATH/ Dataset N x = N y (B j ) 1≤j≤J log 10 (cond(Q * )) Details about the datasets. Each dataset is associated with 200 examples for matrix A * .

	WeathN5a	5	(2, 3)	0.1
	WeathN5b	5	(2, 3)	1
	WeathN10a	10	(5, 5)	0.1
	WeathN10b	10	(5, 5)	1

Table 3 :

 3 Results for climate datasets along with computing times. All the metrics are averaged over the 200 examples of the dataset.

		method	RMSE accur. prec. recall spec.	F1	Time (s.)
		DGLASSO	0.108	0.937	0.894 0.998 0.894 0.937	0.608
		MLEM	0.140	0.413	0.413 1.000 0.000 0.584	0.596
	WeathN5a	GRAPHEM	0.127	0.703	0.595 1.000 0.496 0.742	0.606
		PGC	-	0.772	0.902 0.515 0.953 0.652	0.019
		CGC	-	0.672	0.828 0.285 0.945 0.415	0.026
		DGLASSO	0.166	0.773	0.668 0.992 0.619 0.788	0.630
		MLEM	0.197	0.413	0.413 1.000 0.000 0.584	0.376
	WeathN5b	GRAPHEM	0.186	0.629	0.536 1.000 0.368 0.694	0.470
		PGC	-	0.675	0.677 0.469 0.819 0.544	0.017
		CGC	-	0.634	0.659 0.263 0.895 0.369	0.023
		DGLASSO	0.202	0.948	0.898 0.925 0.954 0.890	1.363
		MLEM	0.264	0.219	0.219 1.000 0.000 0.359	0.834
	WeathN10a	GRAPHEM	0.224	0.511	0.311 1.000 0.374 0.473	1.445
		PGC	-	0.879	0.904 0.504 0.983 0.644	0.232
		CGC	-	0.773	0.539 0.211 0.932 0.278	0.358
		DGLASSO	0.192	0.866	0.633 0.994 0.829 0.769	0.557
		MLEM	0.342	0.219	0.219 1.000 0.000 0.359	0.989
	WeathN10b	GRAPHEM	0.219	0.855	0.620 0.994 0.816 0.757	0.655
		PGC	-	0.799	0.558 0.473 0.890 0.506	0.154
		CGC	-	0.750	0.407 0.218 0.900 0.265	0.178

Table 4 :

 4 Results for sparsity levels s A ∈ {27, 15, 10, 5} of A * , using (P * , Q * ) from dataset A. The first set of results, with s A = 27 identifies with the first row block of Tab. 1.

			Estimation of A	Estimation of P	Estim. Q	State distrib.	Predictive distrib.
	Method cNMSE(µ s RMSE AUC F1 RMSE AUC F1 RMSE DGLASSO 0.061 0.843 0.641 0.082 0.778 0.698 0.083 6.394×10 -8 MLEM 0.076 0.817 0.500 0.105 0.857 0.500 0.102 1.095×10 -7 = 27 GLASSO NA NA NA 0.818 0.804 0.496 1 073.510 4.485×10 -6 A rGLASSO NA NA NA 0.764 0.924 0.598 31.689 2.826×10 -6	1.050×10 -7 1.803×10 -7 7.180×10 -6 5.492×10 -6	2.984×10 -4 4.843×10 -4 1.000 1.000	12 307.169 12 341.205 28 459.294 22 957.693
		GRAPHEM	0.045	0.895	0.847	NA	NA	NA	NA	4.364×10 -6	6.944×10 -6	2.980×10 -4	29 035.030
	s = 15 A	DGLASSO MLEM GLASSO rGLASSO	0.108 0.159 NA NA	0.916 0.933 NA NA	0.781 0.313 NA NA	0.077 0.105 0.463 0.425	0.883 0.899 0.702 0.723	0.855 0.500 0.576 0.598	0.075 0.101 7.088 4.562	8.268×10 -7 1.298×10 -6 2.064×10 -5 1.379×10 -5	1.157×10 -6 2.110×10 -6 3.723×10 -5 2.768×10 -5	5.183×10 -3 6.052×10 -3 1.000 1.000	11 801.220 11 806.483 14 258.010 14 432.765
		GRAPHEM	0.127	0.928	0.721	NA	NA	NA	NA	5.077×10 -5	6.455×10 -5	3.802×10 -3	29 034.715
	s = 10 A	DGLASSO MLEM GLASSO rGLASSO	0.089 0.167 NA NA	0.981 1.000 NA NA	0.938 0.220 NA NA	0.084 0.105 0.365 0.339	0.768 0.900 0.667 0.749	0.686 0.500 0.500 0.699	0.084 0.101 0.469 0.632	9.386×10 -7 2.383×10 -6 3.026×10 -5 1.970×10 -5	1.415×10 -6 3.918×10 -6 5.958×10 -5 4.375×10 -5	9.748×10 -3 2.592×10 -2 1.000 1.000	12 290.789 11 819.739 13 173.462 13 497.061
		GRAPHEM	0.117	1.000	0.745	NA	NA	NA	NA	8.956×10 -5	1.124×10 -4	1.226×10 -2	29 034.727
	s = 5 A	DGLASSO MLEM GLASSO rGLASSO	0.100 0.234 NA NA	0.823 0.893 NA NA	0.774 0.116 NA NA	0.084 0.106 0.242 0.262	0.764 0.856 0.667 0.704	0.682 0.500 0.500 0.566	0.085 0.102 0.306 0.451	9.217×10 -7 2.232×10 -6 1.406×10 -5 1.195×10 -5	1.206×10 -6 3.538×10 -6 3.095×10 -5 2.884×10 -5	9.981×10 -3 5.610×10 -2 1.000 1.000	12 268.667 12 312.749 12 892.144 12 834.737
		GRAPHEM	0.127	0.931	0.573	NA	NA	NA	NA	1.060×10 -4	1.225×10 -4	1.477×10 -2	29 034.456

* , µ) cNMSE(µ s * , µ s ) cNMSE(ν * , ν) L 1:K ( A, P)

x k-1 (1)x k-1 (2) x k=1 (3) x k-1 (4) x k-1 (5) x k (1) x k (2) x k (3) x k (4) x k (5) x k+1 (1)x k+1 (2)x k+1 (3)x k+1 (4)x k+1 (5)

http://www-syscom.univ-mlv.fr/ benfenat/Software.html

https://pages.saclay.inria.fr/emilie.chouzenoux/LogicielEN.html

Acknowledgments É.C. acknowledges support from the European Research Council Starting Grant MAJORIS ERC-2019-STG-850925. The work of V. E. is supported by the Agence Nationale de la Recherche of France under PISCES (ANR-17-CE40-0031-01), the Leverhulme Research Fellowship (RF-2021-593), and by ARL/ARO under grant W911NF-22-1-0235.

The authors thank Gustau Camps-Valls and his team, for providing the ground truth data for (Runge et al., 2020).

the procedure from [START_REF] Moré | Algorithms for bound constrained quadratic programming problems[END_REF], with predefined conditioning number c, with log 10 (c) ∈ {0.1, 0.2, 0.5, 1}, leading to datasets A, B, C and D, respectively.

DGLASSO provides estimates ( A, P) as a direct output. The estimate Q is simply deduced as Q = ( P) -1 . Precision parameters in DGLASSO are set to (ε, ξ) = (10 -3 , 10 -4 ). We initialize DGLASSO with P (0) = 10 -1 Id Nx , and A (0) equals to a stable auto-regressive order one matrix with entries A (0) (n, m) = (10 -1 ) |n-m| , rescaled to have spectral norm equals to 0.99. The setting of regularization parameters (λ A , λ P ) is discussed in the dedicated Section 5.1.4. Performance of the method for varying initializations is discussed in Appendix C.2.

Comparisons to other methods

In our experiments, we also compare DGLASSO with other model inference techniques for LG-SSM.

We consider the EM method from [START_REF] Shumway | An approach to time series smoothing and forecasting using the EM algorithm[END_REF] (denoted after by MLEM) to compute ( A, Q) as maximum likelihood estimates (i.e., no regularization is employed) of matrices (A * , Q * ), the estimation P being defined here as the inverse of Q. Comparison with another state-of-the-art approach for maximum likelihood computation, relying on sensitivity function minimization from [START_REF] Gupta | Computational aspects of maximum likelihood estimation and reduction in sensitivity function calculations[END_REF], is presented in Appendix C.1.

Second, we consider three model inference techniques that explicitly incorporate a sparse graphical prior knowledge on the sought matrices. Namely, we compare to GLASSO [START_REF] Friedman | Sparse inverse covariance estimation with the graphical LASSO[END_REF], that considers a simple static and noiseless version of the LG-SSM (i.e., A is empirically set to zero and observation noise is neglected). Matrix P is then obtained through a maximum a posteriori formulation under an ℓ 1 prior. The convex GLASSO loss is minimized using the proximal splitting method described in (Benfenati et al., 2020, Alg.1). Matrix Q is deduced by inversion of the resulting P. We also compare with the robust GLASSO (rGLASSO) approach introduced in [START_REF] Benfenati | Proximal approaches for matrix optimization problems: Application to robust precision matrix estimation[END_REF], that explicitly accounts for the expression of R k in the maximum a posteriori loss expression. For the sake of fair comparison, we use hereagain an ℓ 1 penalty to obtain P although more sophisticated priors could be encompassed by rGLASSO. For both aforementioned methods, we rely on the Matlab toolbox provided by the authors 3 . Finally, we provide the results obtained with the GRAPHEM method we recently introduced in [START_REF] Elvira | Graphical inference in linear-Gaussian state-space models[END_REF]. GRAPHEM provides a maximum a posteriori estimate A under an ℓ 1 prior, while Q is empirically set to σ 2 Q Id Ny , with a finetuned σ Q > 0. The Matlab toolbox provided by the authors 4 is used to produce the results for this method.

All the comparative methods are programmed in the same language, they are initialized using a similar strategy as our DGLASSO method, and similar stopping criteria and hyperparameter tuning approach are employed, for fair comparisons.

Evaluation metrics

We first evaluate the results of our method, as well as the comparative benchmarks, through quality assessment metrics. Namely, we use the relative mean square error (RMSE) between and

Hereabove, (µ s k , Σ s k ) 0≤k≤K-1 denotes the mean and covariance of the smoothing distribution obtained when running Algs. 1-2 using ( A, P). Moreover, the matrix

is defined in Alg. 2. This concludes the proof.

B. Proximal algorithms to solve the inner steps

We present Algorithms 4 and 5, that are proximal splitting algorithms to solve, respectively, the inner problems ( 32) and ( 33). Specifically, both algorithms are special instances of the Dykstra-like splitting algorithm from [START_REF] Bauschke | A Dykstra-like algorithm for two monotone operators[END_REF] (see also (Combettes and Pesquet, 2011, Sec.5)), for the minimization of the sum of two convex but functions. Sequence (A n ) n∈N (resp. (P n ) n∈N ) is guaranteed to converge to the solution of problem (32) (resp. ( 33)). The proximity steps involved in Algorithms 4 and 5 have closed form expressions that can be found for instance in the reference book [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]. We explicit them hereafter for the sake of completeness.

Proximity of ℓ 1 . Let γ > 0 and V ∈ R Nx×Nx . Then,

which amounts to applying the soft thresholding operator with weight γ on every entry of the matrix input V.

Proximity of quadratic term. Let γ > 0 and W ∈ R Nx×Nx . Then, by definition,

The optimality condition for (87) gives

Algorithm 4 Proximal splitting method to solve (32)

Inputs. A, P, Ψ, ∆, Φ. Precision ξ > 0.

1. Setting. Set stepsize ϑ ∈ (0, 2).

2.

Initialization. Set V 0 = A.

3. Recursive step. For n = 1, 2, . . .:

Output. Transition matrix A = A n .

Algorithm 5 Proximal splitting method to solve (33)

Inputs. A, P, Ψ, ∆, Φ. Precision ξ > 0.

1. Setting. Set stepsize ϑ ∈ (0, 2) and Π as in (34).

2. Initialization. Set V 0 = P.

3. Recursive step. For n = 1, 2, . . .:

If |C 2 (A n ) -C 2 (A n-1 )| ≤ ξ, stop the recursion.

Output. Precision matrix P = P n .