Learning Network Granger causality using Graph Prior Knowledge
Authors: Lucas Zoroddu
Pierre Humbert
Laurent Oudre
Venue: Accepted by TMLR ()

Abstract: Understanding the relationships among multiple entities through Granger causality graphs
within multivariate time series data is crucial across various domains, including economics,
finance, neurosciences, and genetics. Despite its broad utility, accurately estimating Granger
causality graphs in high-dimensional scenarios with few samples remains a persistent chal-
lenge. In response, this study introduces a novel model that leverages prior knowledge in
the form of a noisy undirected graph to facilitate the learning of Granger causality graphs,
while assuming sparsity. In this study we introduce an optimization problem, we propose
to solve it with an alternative minimization approach and we proved the convergence of
our fitting algorithm, highlighting its effectiveness. Furthermore, we present experimental
results derived from both synthetic and real-world datasets. These results clearly illustrate
the advantages of our proposed method over existing alternatives, particularly in situations
where few samples are available. By incorporating prior knowledge and emphasizing spar-
sity, our approach offers a promising solution to the complex problem of estimating Granger
causality graphs in high-dimensional, data-scarce environments.

Published in Transactions on Machine Learning Research (06/2024)
Learning Network Granger causality Knowledge
using Graph
Prior
Lucas Zoroddu
Université Paris Saclay, Université Paris Cité, ENS Paris Saclay, CNRS, SSA, INSERM, Centre Borelli, F-91190, Gif-sur- Yvelle, France
Picrre Humbert
Université Paris-Saclay, CNRS, Inria, Laboratoive de mathématiques d' d'Orsay, F-91405, Orsay, France
Laurent Oudre
Université Paris Saclay, Université Paris Cité, ENS Paris Sactay, CNRS, SSA, INSERM, Centre Borelli, 7-91190, Gif-sur- Yuette, France
Reviewed on OpenReview: ht s: // op en ré net or um DN ut 5f yR
Abstract
Understanding the relationships among multiple entities through Granger causality graphs within multivariate time series data is crucial across various domains, including economics, finance, neurosciences, and genetics. Despite its broad utility, accurately estimating Granger causality graphs high-dimensional scenarios with few samples remains a persistent chal- lenge. In response, this study introduces novel model that leverages prior knowledge in the form of a noisy undirected graph to facilitate the learning of Granger causality graphs, while assuming spansily. In this study we introduce an optimization problem, propose to solve it with an alternative minimization approach and we proved the convergence of our fitting algorithm, highlighting its effectiveness. Furthermore, we present experimental results derived from both synthetic and real-world datasets. These results clearly illustrate the advantages of our proposed method over existing alternatives, particularly sations where few samples are available. By incorporating prior knowledge and emphasizing spar- sity, our approach offers promising solution to the complex problem of estimating Granger causality graphs in dimessiinll,,,al,, data-scarce environments.
1 Introduction
Multivariate time series analysis plays a crucial role in understanding and forecasting real-world phenomena involving multiple interdependent variables. Within such complex systems, Granger causality (Granger 1969) graph has emerged as a powerful tool to uncover directional relationships and causal influences. It has therefore been the basis for wide range of applications in fields such as economics (Stock & Watson 2016), neuroscience (Seth et al. 2015), gene regulation (Yao et al., 2013), protein-protein interactions (Zou et al. 2010) etc. Granger causality assesses whether the past values of one variable can help predict the future values of another variable. For instance, we say that random variable X Granger-causes variable Y if the past values of X contain information that enhances our ability to predict Y beyond what we could achieve using only the historical values of Y itsclf. While Granger causality is powerful tool for studying relationships between two variables, many real-world phenomena involve multiple interconnected variables. Network Granger causality extends the basic Granger causality concept by assuming linear causal relationships between p variables. Hence, allows to capture the complex causal relationships within multiple variables. In this context, the aim is not identify whether one variable inffuences another but to unravel the intricate causal structure that underlies entire system. Network Granger causality is particularly

---

Published in Transactions on Machine Learning Research (06/2024)
relevant in fields like neuroscience where researchers seek to understand how different regions of the brain interact and influence each other over time (Seth et al., 2015), in finance to capture market dynamics or in climatology to investigate complex interactions shaping Earth's climate system. It is also applicable to discover cause-and-effect relationships within genetic pathways and protein-protein intcractions, shedding light on the regulatory mechanisms that govern these biological processes (ao el al. 2013).
A conventional approach to learn Granger causality graphs involves estimating the parameters of a Vee- tor AutoRegressive model (VAR) based on observed multivariate time series [Lütkepohl] 2005). However, inferring parameters of VAR model can be non trivial, particularly when dealing with high-dimensional data and limited samples. To address this issue, rescarchers have studied various regularization techniques, approaching the problem from both frequentist and Bayesian perspectives. From a frequentist standpoint, Basu et al. (015)) explored the application of a Group Lasso penalty and its consistency properties. Since then, several sparsity-inducing penalties have been proposed in the works of Kock & Callot (2015). Lin &
Michailidis (2017). and Nicholson et al. (2020). On the Bayesian side, diverse prior distributions governing the parameters of the VAR model were investigated in the literature. These include the use of Gaussian priors (Sims 1993 Gaussian-invertee Wishart priors (Banbura et al. 20100) and hierarchical normal priors, as explored by Ghosh et al. (2018). However, as underlined by Duan ct al. (2023), the resulting Granger causality graphs from these methods often exhibit undesired characteristics, ranging from excessive density to pronounced disconnection, potentially conflicting with established scientifie knowledge.
To mitigate these issues and enhance the fidelity of Granger causality graphs, researchers have investigated for the incorporation of supplementary information in the form of network structures. For instance, in gene network analysis, genes are often grouped into distinct pathways, and it is common observation that interactions within pathway are more frequent than those bridging different pathways (Marlin et al. 2012) In response, Yao et al. (2013) have introduced penalization terms into the optimization process to ensure that the derived Granger graph aligns with this a priori knowledge, Another approach involves adopting a tree-rank prior distribution, enforcing the graph of Granger causalitics to be a subgraph within the union of spanning trees Duan et al. (2023). More recently, Lin el al. (2024) proposed to solve a Structural Equation Model incorporating constraint to obtain DAG, Furthermore, in scenarios characterized by signals from physical processes and recorded by sensors distributed across multiple spatial locations, the Euclidean k-NN (k-Nearest Neighbors) graph is often leveraged and finds widespread adoption within the domain of Graph Signal Processing (Ortega et al., 2018), where it serves as a foundational element for tasks such as signal filtering, denoising, and prediction. For instance, in (Isufi et al. 2019), a VAR(MA) model is fitted as a polynomial of the Laplacian matrix of the kNN graph. Nonetheless, most of works within the scientific literature operate under the idealized assumption of possessing a complete and precise prior graph an unrcalistic supposition in many real applications.
Contributions: In this paper, we propose an optimization problem for learning VAR model parameters by utilizing prior knowledge about relationships within time series data, represented in the form of an undirected graph (cf Figure D Our approach enables the incorporation of an incomplete and noisy prior undirected network when learning the Granger causality Network, and we jointly learn VAR parameters while denoising the initial graph. To do that, we propose a two-block coordinate descent algorithm and we prove that it converges to a set of stationary points, strengthening the reliability of our approach. Moreover, we show that the optimization problem corresponds the computation of a Maximum A Postcriori (MAP) of a particular statistical model. To validate the effectiveness of our method, we conduct. a series of experiments encompassing synthetic and real-world datasets. Our findings convincingly demonstrate the superiority of our proposed model over vanilla alternatives for varying noise levels over the prior network, particularly in scenarios marked by data scarcity. This underscores the potential of our method as a valuable tool for deriving Granger causality graphs in practical settings across a spectrum of applications.
2 Preliminaries
In this section we present the basics of (Network) Granger causality and its relationships with Vector Au- toregressive (VAR) models. We also give brief overview of existing methods for learning VAR parameters.
2

---

Published in Transactions on Machine Learning Research (06/2024)
Current Methods
Proposed Method
Inpat lime peries
Input time seur
W mmy
Am
mm
Prior weighted Graph
Nm
AG
N MMM MM
N
MM
Figure 1: Diagram of the proposed method. Left:
Existing methods mainly infer causality graphs from time series only. Right: Our proposed method leverage both time series and prior undirected graph (encoding 'likelihood' of link presence) to infer the causal graph
2.1 Network Granger causality
In order to define the Network Granger causality, we recall the concept of classic Granger causality. Definition 1 (Granger causality (Granger 1969 ). Let x and be two stationary time series such that:
y/t = Σ atylt - + + Σ
Aalit - i] + c(t]
1=1
where cit] o²) d is the order of the model, and the parometers ag, Bs fitted minimizing the least square error belween and ils reconstruction using 11. We say that Grunger-causes if one of the (Beb is non zero.
Remark 2. In Granger 1969),
the Granger causality introduced as a statistical test based on the null 3p == O. One can use a F-test to know whether the vull hypothesis if
hypothesis Ho 3== B2
rejected or not.
Remark 3 (Link between causality and Granger causality). While Granger causality can effectively capture causal relationships, its interpretation leans more towards predictive power than direct expression of true causality. For instance, if x and are driven by an unknown 2, il is possible that (x. y) rejects the null hypothesis whereas there is no causality between them.
To generalize the Granger causality p variables, introduce Vector Autoregressive (VAR) models. A VAR model of order d 2 1, denoted VAR(d), explains the values of multivariate time series at time t using linear combination of its d previously observed values.
Definition (Vector Autoregressive Model (VAR)). Given d malrices (CT)d in RexP, a VAR(d) is defined at each time 1,2,... .....p by:
d X = Σ C C -
(2)
7=1
3

---

Published in Transactions on Machine Learning Research (06/2024)
where X = [x. Xyl+)) is random p-dimensional time series and e|t] NV(0, 02), a 0, is some innovalion noise.
Traditionally, Granger causality has been based on the assumption of VAR model Lütkepohl (2005)
primarily conducting tests on the VAR coefficients in a bivariate context. However, in complex real-world systems with numerous time series, analyzing the relationship between only two series can result in misleading conclusions due to confounding factors [ūtkepohh] (2005). In practice, VAR models are often used to analyze
relationships between several variables of interest. Indeed, the matrices [C7)+=1 in Equation (2) capture specific temporal dependencies between the dimensions and are associated with the notion of Network Granger causality (NGC), which was then introduced a generalization of the Granger causality in Basu ct al. 2015 to adjust for possible confounders or jointly consider multiple scries.
Definition 5 (Network Granger causality (Basu et al., 2015 ). Let X = (X1lt) xp(t) E Rº for t = 1, 1,2.... be a p-dimensional lime series following the VAR model defined in Eq. 2 Then, X5[-] is called Granger cause of X,I:] if at least one matriz entry {C] T T= -......) d} is non-zero.
The intuition is that if C # for some T, like the bivariate case, means that past values of X, are useful predict future values of X,, independently from the other variables.
Indeed, we can reduce the model to a bivariate case for each couple (i,j) writing:
X
Σ Σ
X 1 +J++
Σ C x,t
r=l k
and then use the definition of the classical Granger causality.
Note that, like for the classical Granger causality (p 2). NGC does not necessarily capture true causal relationships, but rather indicates the power of prediction of some variables to others. Nevertheless, NGC remains a powerful tool for understanding interactions between random time series, and its estimation is of practical interest.
2.2 Estimation of the VAR parameters
Given N samples [t) E R" n = 1...,, N, at each time t E [1, d] stored in X E RPXN one possibility to estimate the VAR parameters (CT) is to maximize the likelihood Model 2 This leads to the following Least Square problem:
arg min
1 2N
X
Eox
7
(3)
C
C4...
This minimisation can be divided into p independent sub-problems, which are casier to handle. For j 1 ly,,, py they are given by:
1
CJ =|(
arg
min Cd
(4)
2N
where C refers to the j-th row of C Each sub-problent be solved efficiently.
is standard linear problem and can therefore
In high dimensional settings, can still be challenging to fit VAR(d) models because the model has d x p2 degrees of freedom. Some regularization approaches were developed to overcome this limitation, both from frequentist and Bayesian point of view. Several authors suggest adding a penalty term to learn the VAR parameters. Thus, the optimisation problem becomes:
arg
min Od 31
[4
Σ
2 XP, j...... C
(5)
2N
4

---

Published in Transactions on Machine Learning Research (06/2024)
Method
Penalization term
Guarantees
Statistical Model
Ridge
A||u|[2
10)
Lasso
Selection consistency
d, rid
Laplace (0, 1) Laplace(O, me
Adaptive Lasso
A Σ wx/ea po/e
Selection consistency Zou. 2006
ind
&
Group Lasso
A M wy [
Direction consistency Basu
2015
dy g o² ind N (Mrs rg
+ Gamma m1t11 X²
Table 1: Main regularizations for linear regression.
where (P+)-= are some penalization terms.
The Lasso shrinkage technique (I.1 penalization), as introduced by Tibshirani (TTibshirani 1996)), has been integrated into various models to account for specific data properties. In the study by Basu et al. (Basu et al. 2015)), ; for example, the Group Lasso is proposed and explored to leverage knowledge of grouping structures. Lasso shrinkage proves effective in enhancing lag order selection, as demonstrated in the work of Shojaic and Michailidis (Shojaic & Michailidis 2010 using truncated Lasso estimator and in Nicholson et al.'s ( ((Nicholson el al. 2020 hierarchical model.
The Adaptive Lasso (or weighted Lasso), introduced by Zou ((Zou 2006)), finds applications in improving
variable selection (ooo 2006g), addressing non-stationary time series, and proposing online algorithms (Messner & Pinson 2019 by leveraging previous parameters to estimate the next ones. It's worth noting that the conventional approach for utilizing Adaptive Lasso involves initially solving the least square problem (without constraints) to obtain an initial estimator ČOLS, followed by setting the weights wi,j |C,,y
Additionally, Adaptive Lasso has proven to be a relevant method for incorporating prior knowledge about parameters, exemplified by the consideration of past values in ((Mesnne & Pinson 2019)). Another commonly used estimator, the Ridge estimator, is also examined in Giovanni's study (Ballarin 2022
within the context of VAR models, where the authors discuss the implications of anisotropic penalization. It's important to note that for some of these models, theoretical guarantees have been established regarding the selection consistency of Lasso- based methods (refer to Table D
As highlighted in the introduction, some Bayesian models have been introduced learn VAR paramcters by incorporating prior distributions. These priors aim to prioritize certain lags, with lower lags commonly assumed to be more informative than higher ones, or to impose specific structures the graph, such
spanning trees structure ((Sims) 1993 Banbura et al., 2010 Ghosh et al. 2018 Duan et al. 2023 ) A
comprehensive survey of existing Bayesian VAR models and their associated sampling algorithms is presented in the work of Miranda Agrippino and Ricco ((Miranda Agrippino & Rieco 2018
However, note that in some cases, frequentist penalizations correspond to compute the Maximum Likelihood Estimator (MLE) of a particular Bayesian model, so it can be useful to consider a model from both point of view to have deeper understanding of the solution. We refer to Table 1 for the correspondences between conventional penalization terms and their associated Bayesian models.
3 Model and Framework
In this section, we present our main contributions. In Section 3.1 we first introduce an optimization problem allowing to incorporate graph prior knowledge in the estimation of the parameters of VAR model. Then, in Section
3.4 we describe an algorithm based on alternating minimization to solve this problem. In section 3.5 we prove that this algorithm converges to a stationary point and that its time complexity for fixed number of iterations is of the same order as a Lasso estimator. Finally, we link this op- timization problem to a maximum a posteriori of a certain statistical model, which presented in Section 3.6
In the following, we make the two classical assumptions:
5

---

Published in Transactions on Machine Learning Research (06/2024)
Assumption 1. The time series are generated by VAR(1) model:
Xit| = C. - 1] + =(t]
(6)
so we only need to learn one matriz C. However, note that the assumption d = 1 is not limiting since a VAR(d) model can always be written as a VAR(1) model (Lütkepohl, 2005). This generalization is detailed Section 3.7
Assumption 2. Homdlton (1994) showed that obtaining reliable estimates for VAR parameters necessitates
stationary of [ Subsequently, in Littkepohl (2005) work, was established that the model G is stable if and only (CC) 1, where p is the spectral radius. This assumption will be maintained in our analysis.
3.1 Problem formulation
In general, the cstimation of the parameters of the VAR model (2), i.c. the matrix C, requires the observation of a long stationary realization ion of the p-dimensional time series. However, in many applications, we observe only short replicas of the time series, and additional information must be incorporated into the model to obtain accurate estimates.
We propose to leverage prior knowledge on the structure the matrix C by assuming that is statistically related to a given matrix A prror which is the adjacency matrix of an undirected graph encoding a priori relationships between the individual dimensions. More specifically, the idea is to encode the following property: if two nodes (i,5) are not likely to be linked, Aprior is small and the value of the associated coefficient Cag is more likely to be closed to zero (meaning that there is no Granger causality between the two time scries X.[ and x)(.)). While forbidding certain links in causal discovery is standard approach and it used in recent methods like PCMCI Runge et al. (2019), ; this assumption may be loo strong and this issue will be addressed in the paper (ef Remark Sh. In addition, since the prior information is rarely perfectly accurate for most applications, we assume that A prrio is not necessarily the optimal prior knowlodge. So we want to allow the model to refine this prior through iterations, i.e. computing a matrix A close to A but not equal. To understand the relationships between these matrices, statistical point of view is derived Apior and A.
in 3.6 and provides the mathematical relationships between C,A prior
Formally, let consider N independent wide-sense stationary multivariate time series (Assumption 2 x( dl.... (1 d]. In addition, suppose that we have access to a matrix prior, summarizing our prior knowledge on the relationships between each pair of variables. To estimate the VAR parameters, we introduce the following optimization problem:
min A C
1
Σ n=1
N
x - c - 1]
+
Σ
Ca,l + CCaal
N
Aij
ISisjšP
+ 2X
Σ
log(2A,,) ( + 1A A prior 12
(7)
ISi<jšP
subject to Ai, 2 0, Ai, = AJi isi<jšp.
Equation
contains 3 terms:
1 x cx I]ll² corresponds to the Least Square problem objective: this term al- N
ii)
M
#=1
lows to measure the difference between the original signals and their reconstructions. Recall that we only consider in this soction VAR(1) model, hence
Cayl+ IC +2 Σ log(2A.,)) the penalization term that takes into account the ISi<j5P
(ii) E
ISKKSS A
graph prior knowledge. This penalization is inspired by the one used in Adaptive Lasso models, where the terms 1/A,j act as weights. The higher the A the closer and j are in the graph, and the lower the penalty for C. and C,,-- It should be noted that the additional term composed of the sum of log is normalization term linked to the associated statistical model (see Section 3.6

---

Published in Transactions on Machine Learning Research (06/2024)
(ii) ||A - A is regularization term to take into account that A is not necessarily the optimal prior knowledge and could therefore be further refined. Actually, adding this term in the optimization problem is equivalent to imposing a normal prior distribution to the coefficient of Ai,j (see Section
3.6 In practice, this term allows to increases the robustness to the prior knowledge noisc.
The symmetry constraint A., = A3, for 1 j S p applied as represents the adjacency matrix of the denoised undirected prior graph knowledge.
Remark 6. This work presents a model leveraging prior undirected weighted graph. This model choice
is based on the idea that for some applications, we have ideas about which variables are linked together but we do not know the directions of causality (i or j - i?). However, it is possible to slightly modify the model to leverage directed graph prior knowledge (the optimization problem remains the same except that we split. the part with Aij and that we remove the symmelry constrainl).
Finally, the problem 7 depends on two hyper parameters and 7.. À controls the sparsity of the learned graph and y controls the confidence the prior graph.
3.1.1 Interpretation of the optimization problem
A prior is pXp symmetric matrix, seen as the adjacency matrix of an undirected graph, corresponding to the prior knowledge of the Network Granger causality structure.
C is p X p matrix corresponding to the VAR(1) paramcters and seen as the adjacency matrix of the directed graph of Granger causality.
A is p p symmetric matrix, seen as the adjacency matrix of an undirected graph, corresponding to the denoised prior graph A prior during the optimization process.
x (t = 1,2) are the multivariate time scries in Rp' generated by the VAR(1) model with param- eters C.
The idea of this optimization problem is to find a sparse matrix C allowing to forecast X using A as a first estimation of the parameters model. The main advantage here is that the prior network is taking into account during the learning process and allows to perform a relevant variable selection. Indeed, even though the variable selection consistency proved for estimator like Lasso or Adaptive Lasso (cf (Zoul (0066), this property holds for an infinite number of samples. Here, by assigning weights pairwise relationships, our method allows to guide the variable selection using (noisy) prior knowledge in order be efficient in settings with few samples available.
3.2 Motivating examples
Before entering in the details of the solving method, we illustrate here some cases for which our method could be used, and how the graph A prior can be constructed depending on the problem.
Multi-sensor networks: In the context of multi-sensor networks, A can be computed by focusing on the distances between the sensors. For cortain applications, the recordings are associated with physical processes, suggesting that causal links are more likely to exist between nearby sensors than distant ones: for instance, temperature data are governed by complex equations but exhibit spatial continuity. Consequently, in the experiment with the Molene dataset (Sec. 4.5.2), prrir is computed applying Gaussian kernel to the pairwise distances: Prior K2 (d(x,,,))) in order to increase the penalization of causal links between
far away sensors (cf Figure] 2
Expert knowledge from the literature can be used to construct a score that reflects how likely two proteins interact based on previous experiments. In this case, prior is the score the couple (i.j). This can lead to prior knowledge graph about pairwise proteins interaction allowing to lead the cansal discovery process. An interesting example the one presented in Carlin et al. (017) by aggregating various established pathways from Pathway Common Cerami al. 2010)
Protein interactions:

---

Published in Transactions on Machine Learning Research (06/2024)
Landeriia
Doufrment
Pontig
Cars arnind
Aarient
Pleemo
Guer
Nazaire
Baule
E EECUUbI IN
Figure 2: Prior graph for the multisensor network example.
Intracardiac electrocardiograns: In this scenario, we are measuring electrical activities in the heart using multiple electrodes. While this example falls within the case of multi-sensor networks, and the prior graph can be constructed as previously explained, A prior can also be computed by focusing on the similar- ity between signals recorded by the electrodes, using cross-correlation for example. allows guiding the algorithm to find causality links between the most similar signals (ef Figure 3
Compute pairwise cross correlations
Input time series
Prior Graph
Ca.
Csa
X5
X3
Cya
CLa
C2.4
X4
X5
Figure 3: Construction of the graph prior for intracardiac signals example.
3.3
Links with existing models
The work of Yao al. (013) is the closest one to ours. They propose leverage prior knowledge A prior p solving:
Č argmin [X[t] CX 1]13 + I -
(8)
with N M + || or N M |M] However this model works under the assumption that we have prior knowledge about the exact values of the matrix C, which can be unrcalistic in some applications. Consequently, this framework impose to know the sign of values in C since the optimization problem penalizes the component wise difference between C and prior The authors propose to address this issue by hand, computing Ridge estimator of the VAR model chose the signs of A prior On the contrary, by exploiting weighted Lasso, our model does not suffer from these limitations and only require relative prior knowledge about relationships
8

---

Published in Transactions on Machine Learning Research (06/2024)
Our method is also closed to an Adaptive Lasso estimator obtained with:
Cayl+ C
min Σ - C
Σ ISKCsSP
(9)
C
A
Indeed, considering only the two first terms of our problem, we obtain exactly the Adaptive Lasso one. However the Adaptive Lasso estimator is not robust to the errors in the weights, since the asymptotic consistency is achieved taking for weights unbiased estimators of the true solution Louu,0000 Thus, the last term in the optimization problem allows to refine the prior through iterations using information in the time series resulting in a more robust version of Adaptive Lasso as will be seen the experiments (Section
3.4 Solving method:
A-AdaptiveLasse (AALasso)
The function to minimize (7) is not convex in (A, C). Indeed, we can show that even with fixed, the function in A is not convex (c.f. Appendix Aj. However, the function is convex in C with A fixed (Adaptive Lasso problem) and have closed form for the roots of the derivative in A (with C fixed). Thus, alternating minimization algorithm to solve this problem.
3.4.1 C uuddate.
For fixed A, the optimization problem
7 with respect to C is:
N
Ciyl + Cjal
min
x 60 tl
C
+
Σ
(10)
N
A
ISKcjS
From Eq. 10 we see that the optimization step in C is an Adaptive Lasso problem with weights equal to 1/A.. Zou 2006). A common way to solve Adaptive Lasso regression is to remark that it can be written like Lasso problem. Indeed, considering Ci2 Cia we can write the problem Lasso one transforming X.
A.,
While there is no closed formula compute the Lasso estimator in the general case, common ways to solve it are to use: (i) the least-angle regression (LARS) algorithm Efron et al. (2004), (ii) algorithms based on coordinate descent Wu & Lange (2008), (iii) the ADMM algorithm Boyd (2010). A recent survey presents
some of these algorithms and provides their convergence rates (Zhao & Huo 2023)
Details of the ADMM
updates solve Problem 10 are presented below. Let E [1,p],
we first rewrile 10
as p subproblems:
p
Ca,ll
min
|x|||] - C,Xlt 11l1
+
i= i...,..,
(11)
N
A
with X,[t] and X,[t- 1] the vectors containing the N samples of the variable i. Then we rewrite Lasso problem:
as a
min C::
N -
- 1111 + Σ Casl i= 1=1,..,
(12)
where C.2- C.a and X|t =1]= x Xt.j)ssssp- In order to use the ADMM algorithm, we rewrite the A.
problem as follows:
| - Č, - 1]]2 + + subject to C, = 0.
min C.
N
Finally, the ADMM updates are:
C arg min 1 ~ [|x[[t] C - 1112++ _IC, - - w|2
C,
9

---

Published in Transactions on Machine Learning Research (06/2024)
(X(t -
p 1X + P w))
U+ arg Al|111 -
SN + W)
-soft-thresooddin C + W)
W W + + -
Bj -t if x
where [Sc(
0 if -i5s,, j j
at <
and W,U are auxiliary variables specific to the ADMM algorithm.
3.4.2 A update.
For fixed C, the optimization problem 7 with respect to A is:
min
|
+
22
Σ log(2A,.)) + 7 || A prior
(13)
A
Aij
F
ISKCSSP
ISKCSSP
subject to Arj 2 0, A.j A,s isi<jsp.
To address the symmetry constraint, straightforward way is to optimize over the upper diagonal values and to set Aja Aij for i < j. The minimisation can then be carried out by directly calculating the exact minimum, which is given in the next proposition.
Proposition 7. The roots of the derivative with respect to Alm of the objective function
13 are:
prior m.
1
2k
+ e 2ikr/3
+ 23/T
-9 q
k
0,1,2 (14)
2
27
(A 1m
p A
2y
where
pelor
8y 2
22 C[(Ci, + |Cm,,l)
9
A = 4p² + 27q²
Furthermore, there erists al least one posilive root for the derivalive with respect to Alm, and the global minimum on the interval J0, +oo is attained at one of these roots.
Proof of Finding the roots of the objective function derivative with respect to A equivalent to finding the roots of 3-degree polynomial, allowing us to apply the Cardan formula. Subsequently, the objective function in As,m diverges towards infinity at the boundaries of j0, +o0[, there exists at least positive root, and the minimum is allained at one of these roots.
3.4.3 Alternating Minimization algorithm.
The final AALasso algorithm is presented in Algorithm AALasso depends on four input parameters:
(i) À controls the sparsity of the learned graph. The larger X, the more sparse the solution.
(ii) 7 controls the confidence in the prior graph. A large y will constrain A to stay close to A prior whereas small value allows A to deviate from Aprior
(iii) The choice of (0) has an impact on the solution since we start by solving the problem C with fixed A. The straightforward idea is to take A (0) prior but note that other initializations can be chosen fc.f. Appendix B.3 For instance, AATasso be seen generalization of Lasso and
10

---

Published in Transactions on Machine Learning Research (06/2024)
Algorithm 1: Fitting algorithm.
input: Nittrr X, A prior
output: C, A
A (0) + Aprior
for to Nuer do
C fc(C, (r-1)) where fc denotes the update in (3.4.1
A (A(O)), A) where SA denotes the update in 3.4.2)
return C (Ntwr) A (Nnor)
ALssL since these two estimators correspond to the first iteration of AALasso for particular choices of A ( Indeed, taking ( equals to the one matrix (1)sa.sp the C update corresponds to the Lasso algorithm, while taking A ( (wa,))sujsp where way are the weights computed by the Least Square Estimator, the first step corresponds to the LS+ALasso estimator.
(iv) The number of iterations Miter of the alternating minimization algorithm impact directly the runtime and the performances, and in practice we will chose a relatively small number of iterations (around 10), according to synthetic experiments conducted in 4
Remark 8. Note that if it erists on iteration to such that A O, the coefficients C and C will be zero for >> Tu- Nezt, in order to allow the algorithm to add edge that is not originally present in the prior graph, we set the minimum values of the adjacency matrix to 0.
Remark As explained in Remark G is possible leverage directed prior graph with slights modifi- cations of the model. Regarding the solving method, the update remains the same, and for the A update, we compute the solutions of p² 3-degree polynomial rather than p(p 1)/2).
3.5 Theoretical properties
In this section, we prove the convergence towards set of stationary points of our algorithm to solve the optimization problem 7 and we show that the time complexity asymptotically in the same order than the one of Lasso estimator.
3.5.1 Convergence
Classical results of alternating minimization convergence assume that the objective function is differentiable (c.f. Grippo & Sciandrone (2000)). However it is not applicable to our case since objective function in is not differentiable in C.2 O. We now introduce the lower directional derivatives to address this issue. Definition 10 (Lower directional derivative). For any x IPP and any U RP, we denote the (lower) directional derivative of f at x in the direction by f(( lim inf ((=+x))1(x)) XTO
Definition 11 (Stationary points). We say that z is a stationary point of f if 2 E dom and f'(z;)) O. Ve.
We say that is a coordinate-wise wise minimum point of f if E domf and f(z + (0..., Uk...., .)....,)) 2 f(2),
Yug R for all h: = 1,..,, N.
Using this framework, the following theorem holds.
Theorem 12. (Convergence of A ALasso)
The sequence ((C(r),A A())) generated by the two-blocks alternating minimization is well defined and bounded. Moreover, every cluster point is a stationary point of the problem
Proof. The proof is given in Appendix A.2
11

---

Published in Transactions on Machine Learning Research (06/2024)
3.5.2 Computational complexity
Recall that when fitting a Lasso estimator for learning VAR parameters in p dimension, we fit p independent Lasso estimators dimension p. In the following, p denotes the dimension and N the number of samples. Recall that our algorithm iteratively solves two steps, A and C. First, as detailed in [.4.1] the C update is solved by transforming the Adaptive Lasso problem into a Lasso problem in O(p N) and using ADMM to fit the Lasso parameters in Opp² N x p²) (see A.3 for more details). Thus, the time complexity of the C update is in Opp² + N x p²). Then the A update is done computing closed formula in ((1) for each value Aij, 1 I < j S p, so this step is in O(p²), which is negligible compared to the C update. Finally, considering a fixed number of iterations Niter, the AALasso algorithm has a time complexity in (Ne x p² + Niter x N x p²), which is approximately Niter times the one of the Lasso estimator solved with the ADMM algorithm.
3.6 Link with Statistical Model
It interesting to adopt the Bayesian point of view to deeper understand the statistical hypothesis behind the model and how the algorithm works. Here, it can be shown that the optimization problem 7 presented in the previous section is equivalent to the Maximum A Postcriori (MAP) of the probabilistic graphical model presented above in Eq. 15
A²
C
X
A
with
Aij ~ N(A i.j = 1, .....
Ci ~ Laplaco 0, Ai i,j = 1, 1.... X X(t 1},03)
(15)
Figure Observable variables are in grey and latent variables white.
The solulions of Problem 7 correspond to the Masimum A Posteriori (MAP) of the statis- tical model 15 under the assumption that A follows the improper distribution 1s,(B).
Theorem 13.
Proof of 13 The proof is given in Appendix A.1
Then, the normalization term log(2AA.) introduced in corresponds to the normalization of the Laplace distribution. Intuitively, it allows A to remain meaningful regarding the Laplace distribution and it avoids parameters to go to infinity.
This probabilistic graphical model allows a good understanding of the model introduced. An unknown graph g" first generates parameter matrix (Ct from Laplace distribution with variances equal to (A (note that Laplace distribution encourages sparsity). and C* allows generate the multivariate time series X following VAR(1) model. On the other side, we observe noisy (Gaussian noise) version of * To
12

---

Published in Transactions on Machine Learning Research (06/2024)
leverage the information provided by A to learn C* we propose to jointly infer the two latent variables A and C
3.7 Generalization to VAR(d) models
In the proceding section, we discussed VAR(1) models for simplicity. However, certain applications require the consideration of VAR(d) models and we detail the generalization of our model in this section. A VAR(d) model can be expressed as a VAR(1) model as follows. Let X a process defined by a VAR(d) model:
x = Σ C C2 =|t],
where X[t] (Xi Xplt]) is random p-dimensional time series and e[t]
N ox0,)x4p). ox 0. For
C
C
Cd'
1
0
0
2 d, let X (X X[t 1)..... X d + 1})² and C
then the process X(t)
0
satisfies the VAR(1) model:
X CX(t - 1|+ (clt),....cd dt 11)².
Note that the stability assumption of the VAR model is now satisfied if p CC) < 1 (cf (Lütkepohl 2005
Using this last point, our model is still applicable for some d E N assuming we have a prior matrix defined by:
A prior
A prior
A prior
1
0
prior
0
0
Thus, solving this problem with number of lags d 1 is equivalent to solve a problem with lag in new variables.
dimension p d, since we can consider the lags Note that straightforward idea is to set Aptior A prior for [ = 1..... d but this generalization allows to leverage prior knowledge about relationships likelihood at each order.
4 Experiments
In this section, large series of experiments is carried out to assess the effectiveness and applicability of AALasso to learn Granger causality graphs. Our algorithm is tested using both synthetic and real-world data sets. To utilize AALasso effectively, the input requirements include a multivariate time scries and the availability of an adjacency matrix representing a prior network structure. By employing synthetic and real datasets, we evaluate AALasso's robustness across various scenarios, including limited number of samples and several levels of noise, and real-world datasets. These experiments provide valuable insights into the algorithm's capabilities and its potential applications in diverse domains. The code reproduce our experiments with synthetic data will be available.
4.1 Task and evaluation metrics
In these experiments, the objective is to learn Network Granger causality (NGC) from given multivariato time series and a prior nelwork. We suppose that these series follow a VAR(1) model, hence, learning the NGC is equivalent to fit the VAR parameters, ie., given X R and prior x we want to estimate the matrix C.
Since VAR models are usually employed for forecasting tasks, a standard metric evaluate estimators is the normalized Root Mean Square Error (nRMSE) of the step predictions. Let X be multivariate
13

---

Published in Transactions on Machine Learning Research (06/2024)
time series and [t] be the reconstruction using the fitted VAR model time t, the nRMSE is defined by:
Σ, [ X Σ, [
nRMSE(XX)
Note that, the main objective of this paper is to learn the underlying NGC, so are more interested in learning relevant graph than allowing a good reconstruction (even though the two tasks are correlated). To evaluate the quality of the learned graph, we compute the F1-score between C and C (sce c.g. (Pasdcloup) et al. 2016 for more in information). This metric is defined follows:
2 precision recall
F1-score =
precision + recall
where
True Positives
precision
True Positives False Posilives
and
True Positives
recall True Positives False Negatives
Here, the precision measures the proportion of correctly identified causality links, while recall measures the ability to capture all causality links. Note that the F1-score is only available for synthetic data need to have access to the true graph (the true VAR model).
Although these two measures are related (a good graph should lead to good reconstruction), it should be noted that a good reconstruction can be achieved by relatively dense graph. Given that sparsity is a desired property, the F1 score is used to understand whether the learned graph can efficiently reconstruct time series while avoiding irrelevant edges.
4.2 Methods
The aim these experiments is to show that AALasso can exploit prior knowledge to improve its performance compared to existing methods. We compare our estimator to the classical estimators: the Lasso and the Adaptive Lasso with weights equal to the least squares estimates (noted LS + ALasso, cf (Zou] 2006)), ). Moreover, since the first step of our algorithm is equivalent to solve an Adaptive Lasso problem assuming that the weights are given by W.j A compare method with this first step (denoted 1-AALasso)
to demonstrate the usefulness perform several steps. Note that do not show the performances of the Least Squares estimator since the results poor in settings with only few samples. Finally, note that the Lasso and ALaso algorithms do not take into account the prior matrix, so the prior noise will not impact their results.
Implementation details. For all experiments, we used the package asgl (Álvroo Méndez Civicta et al. 2021 implemented using cvxpy (Diamond & Boyd 2016 to solve Lasso and Adaptive Lasso regression problems.
4.3
Datasets
4.3.1
Synthetic Data
Synthetic data are generated with respect to the statistical model 15 To define the matrix A', first
generate 40 points in [o. 1]2 uniformly at random. Then, we construct a matrix D E R.*P by applying the Gaussian kernel to the pairwise Euclidean distance between the points (the standard deviation of the
Gaussian kernel is taken equal to the median values all pairwise distances). A is obtained by randomly setting to 0 ratio 7m 0.5 of values of D (mispecified edges) and cutting 0 values smaller than 0.7 to promote sparsity. Finally, the VAR parameters are drawn from Laplace A?,) and A prior D + E, where & is a symmetric matrix whose subdiagonal values sampled from i.i.d. Gaussian distribution with
14

---

Published in Transactions on Machine Learning Research (06/2024)
Algorithm 2: Data generation.
input : :p: Tmx T, T, TA
output: A Aprior X
Generate randomly p points in [0, 1]3
Compute the euclidean distance matrix D.
A' is obtained setting to 0 ratio of Tm values of D (mispecified edges). Generate VAR parameters C. ~ Laplace(O,, A.))-
for 1 S i j S p do Randomly set Ci or Cj to (directed graph).
A prior- D & where & is a symmetric matrix where subdiagonal values are sampled from i.i.d.
gaussian distribution with variance 2
Sample X(t) ~ WN(CC' 1]
return A prior X
variance o'A (varying in {0.02, 0.1, 2.1,0... 0.25.,....... 0.35) to test several level of noises). Note that Aprior is computed
from D and not from A Indeed, A is obtained by performing sparse variable selection from D; therefore, we do not incorporate this variable selection into A prior and the objective to assess the effectiveness of AALssso in accurately retrieving this selection. This data generation is summarized in Algorithm (2) At the end, for each experiment, sample N = {80, 200...., 500} different time series X[t] ~ X((C' 1]. ox),. xx 0,1, which we split into training and test sets of equal sizes. We repeat this procedure 20 times for each value of N.
4.3.2 Breast Cancer Network
Dataset The Heritage Provider Network DREAM & Breast Cancer Network Prediction dataset focuses on prodicting causal protein networks using time scries data from reverse phase protein array (RPPA) experiments. The aim to advance breast cancer understanding by using complex time series data and computational modeling to uncover causal relationships within protein networks responding to various stimuli and inhibitors across different cell lines. It involves examining four cell lines (BT549, BT20, MCF7, and UACC812) under four inhibitor conditions (AKT, AKT + MET, FGFR1 + FGFR3, and DMSO control) and eight ligand stimuli (Serum, PSB, FGF, Insulin, FGF1, HGF, NRG1, and IGF1) at multiple time points (t tt=00 0, 5 min, 15 min, 30 min, 1 hr, 2 hr, and 4 hr). Here, the task is to create 32 causality networks, one for each combination of stimulus ligand and cell line, from protein probe across the time points. Formally, given a 3-uple (cell line, ligand, inhibitor), we observe a multivariate times series X RpxN with N = and p (41, 45, 48) and want to learn a graph of Granger causalities g.
Choice of prior In (Carlin al. 2017). the authors successfully utilized a prior network derived from the Pathway Commons database version 3 Cerami et al.) (2010) to enhance their analysis. In essence, this network prior was developed by aggregating various established pathways from Pathway Commons, where protein interactions closely aligned with the concept of causal influence. It assumed that interactions declared in these pathway databases implied that perturbations to upstream regulatory proteins could lead to either direct or indirect perturbations of downstream target proteins connected via directed paths. The adjacency matrix was obtained performing a heat diffusion over the initial graph (ef (Carlin et al. 2017 for more dctails). Importantly, the network prior was independent of training data and could be reused across experiments. The utilization of this network prior improved performances, highlighting the importance of taking into account prior knowledge for the causality inference task. However, in (Carlin et al. 2017) the authors did not take into account this network in their inference algorithm: they averaged a posteriori the output of the GENIE3 (Huynh-Thie Thu et al.] 2010) algorithm with the adjacency matrix of this prior. We therefore chose A prior equals to this adjacency matrix.
Note that the objective of this challenge to infer causality but not necessarily cemporal/Granger Granger causality. Thus, we do not expect achieve better results than the challengers but the objective is to demonstrate that taking into account a prior can be relevant.
15

---

Published in Transactions on Machine Learning Research (06/2024)
Pre- processessn of time series We normalize the time series pre-processing step to satisfy the first-order stationary assumption:
X 7 + x(1: x
(16)
ox
where is the mean of X and ox its standard deviation.
4.3.3 Molène Dataset
Dataset The Molène datase contains hourly temperatures recorded by sensors at p = 32 locations Britany (France) during N = 746 hours. Here the objective is to understand the spatio-lemporal dynamics of the temperature and assess the extent to which the model can describe complex phenomena (such weather) by considering only data and geographical information (sensor positions).
Choice of Aprior The prior adjacency matrix is A exp(-dist(i, ji, where dist(i,j) is the Euclidean distance between the stations i and and dist is the median of all computed distances. Note that this case good example of what could be a 'noisy" prior knowledge, since the euclidean geometry is isotropic contrarily to the weather dynamics.
Pre- processing of time series We consider the first derivative of the signals rather than original signals in order to verify as much as possible the widc-sense stationary property.
4.4 Results synthetic data
4.4.1 Comparison with classical algorithms
In this part, we compare AALasso with standard method compute Network Granger Causality. Thus we focus on methods based solving the classic NGC optimization problem with some regularization (Lasso, Adaptive Lasso). The aim of these experiments is to show that our method is able to leverage prior knowledge to increase the accuracy of the retrieved NGC.
For all of the 20 experiments, we performed Niter = 10 (see 4.4.4 iterations in the alternating minimization algorithm ( using half of the training set. The parameters and were selected by cross-validation minimizing the nRMSE over the second half of the training set (the validation set).
0.200
1-AALasso
0.8
AALasso
0.175
Prediction Error
LStALasso
F1 SCORE
0.150
Lasso
1-AALasso AALasso
0.5
LS LLasoo
0.4
Lasso
0.025
50
100
150
200
250
100
150
200
250
Samples
Samples
Figure 5: rNMSE and Fi-score in function of the number of samples used for training using 4 = 0.1. We plotted the 90% confidence intervals.
The results in Figure 5 exhibit better reconstruction and greater Fi-score when utilizing AALasso rather than vanilla methods when the number of samples is lower than 140. From 40 to 140 samples, AALasso
16

---

Published in Transactions on Machine Learning Research (06/2024)
0.7
1-AALasso
AALasso
ALssoo
0.5
Lasso
F1 score
GT sparsity
0.4
0.2
0.0
Sparsity
Figure 6: Fl-score in function of the sparsity of the learned graph a number of edges p(p-1)
(obtained using several values for A) with σл 0.1 and N = 40. A sparsity of 0 means an empty graph while sparsity of mcans a dense graph. We plotted the 90% confidence intervals.
returns Fi-scores from 0.6 to 0.8 while LS+ALass Fi-score ranges from 0.5 and 0.73 (an average gain of 0.1). Thus, our algorithm effectively leverages the additional information in settings with few samples, and approach enables fine-tuning of the graph while remaining good forecasting power. Results presented in Table g provides more details regarding the computation of F1-score (precision and recall). We observe that the gain in the F1-score is a consequence of a important gain in precision. Indeed, whereas the first iteration of the algorithm results in relatively good recall and precision, the next iterations lead to important improvements of precision with limited loss in recall. Moreover, the results shown in Figure D indicate that AALasso outperforms Lasso and ALaso estimators for sparsity levels near the ground truth one (with manual selection of X). These results confirm that the model allows to increase the performances of the variable selection and that the gains in Figure 2 are not just consequence of a better selection of X.
When the number of samples increases, the LS+ALasso estimator provides better reconstruction than AALasso. This behavior can be explained by the fact that the prior knowledge given is noisy so the AATasso estimator is biased. In practice allows to roduce the variance when few samples are available, but when the number samples becores large enough to perform statistical inference directly from time series data, this bias lead slightly less accurate results. Thus, the question is know whether the time series are informative to know if adding biased prior knowledge will improve not the results. However, recall that are interested in graph learning task the Fi-score is more informative than the reconstruction error, and shows satisfying results even for relatively large number of samples (until certain threshold, here 250 samples). This can be explained by the fact that the algorithm prefers precision to recall, and will focus on returning true causal relations, sacrifying reconstruction and recall. Finally, the difference rcsults between the first ileration and the complete optimization process of A ALasso points out the interest. of the alternating minimization.
N Samples
100
180
Metrics
P
F1
R
Lasso
0.28 (0.04)
0.61 (0.09)
0.38 (0.04)
0.51 (0.1) 0.75 (0.08)
0.6 (0.09)
0.57 (0.08)
0.83 (0.06)
0.68 (0.07)
ALasoo
0.54 (0.09)
0.54 (0.11)
0.53 (0.07)
0.69 (0.1)
0.68 (0.1)
0.68 (0.09)
0.64 (0.08)
0.83 (0.06)
0.72 (0.07)
A-AALasso
0.44 (0.06) 0.73 (0.1)
0.55 (0.06)
0.65 (0.06)
0.82 (0.09)
0.72 (0.06)
0.83 (0.06)
0.84 (0.06)
0.83 (0.05)
AALasso
0.55 (0.07) 0.71 (0.11) 0.62 (0.07)
0.76 (0.06)
0.79 (0.1)
0.77 (0.07)
0.87 (0.05) 0.81 (0.07) 0.84 (0.05)
Table 2: Procision, Recall and F1-score in function of the number of samples. We took noise over the prior
network with a 0.1.
17

---

Published in Transactions on Machine Learning Research (06/2024)
Experiments in higher dimension Finally, to understand whether our method was able to deal with high-dimensional environment (i.e N << p), we conducted the same experiments taking p E {60, 100, 160 and N = 40 (still 20 samples for training and 20 for validation).
Dimension
60
100
160
Metrics
P
F1
P
R
F1
R
F1
Lasso
0.31 (0.07)
0.48 (0.09)
0.37 (0.07)
0.26 (0.04)
0.32 (0.06)
0.28 (0.04)
0.26 (0.03)
0.24 (0.04)
0.25 (0.03)
LS+ALasso
0.58 (0.06)
0.38 (0.10) 0.45 (0.08)
0.52 (0.06)
0.20 (0.05)
0.29 (0.05)
0.48 (0.05)
0.13 (0.04)
0.20 (0.05)
1-AALass
0.45 (0.05)
0.60 (0.10) 0.51 (0.06)
0.40 (0.04)
0.47 (0.06)
0.43 (0.04)
0.39 (0.03)
0.38 (0.05)
0.39 (0.04)
AALasso
0.56 (0.06) 0.60 (0.11) 0.57 (0.08)
0.47 (0.05) 0.45 (0.07) 0.46 (0.05)
0.45 (0.02)
0.37 (0.06) 0.40 (0.04)
Table 3: Precision, Recall and F1-score in function of the dimension. We took noise over the prior network
with A 0.1 and N = 40 samples.
From Table we see that in settings where N << p. AALasso remains belter than the other standard methods regarding the Fl score. More of that, the higher the dimension, the higher the gap between LS+ALasso and AALasso regarding the Fl score (from less than 1.2 times better in dimension 40 to twice better in dimension 160). Thus, our method efficiently leverages prior knowledge, and it is of high interest in this kind of settings.
4.4.2 Influence of the prior network
In this section, we present results for various configurations of prior noises. Similar to the experiments in the previous section, for each configuration, we conducted 20 experiments, we took Niter = 10, we utilized half of the training set and the parameters À and yy were selected via cross-vale validation,,,, minimizing the normalized Root Mean Square Error (nRMSE) over the second half. Figure 7 presents the results for Prediciton error and Fl-score with N 40, 100, 140 samples for varying noise levels. It is important assossing the robustness of our method to prior matrix noise, where the noise specifically corresponds to A prior noise and not to the noise of the VAR model. Additionally, since Lasso and LS+ do not leverage prior knowledge, variations in this noise do not impact the results. The findings demonstrate that AALasso exhibits robustness to noise, displaying a comparable prediction error than the Lasso or LS+ ALaso one (even better than AALsso in few samples settings) and a consistently better Fl-score for all tested configurations. Furthermore, these results indicate that our model effectively generalizes Adaptive Lasso (corresponding to the first iteration) and enables refinement of results over subsoquent iterations.
note that we are
Tables 4 and 5 provide further insights into the high Fi-score achieved with AALasso for both N 40 and N 140 scenarios. The AALasso's behavior remains consistent across all tested noise levels. While the first iteration yiclds a good recall but limited precision, subsequent iterations lead to significant increase in precision (an average gain of 0.15), especially in high noise level settings (gain of 0.2) while maintaining a good recall (an average loss of 0.05).
4.4.3 Comparison with other Causal Discovery algorithms
While we focused in this paper on learning Network Granger Causality, is interesting to compare our method with other causal network discovery models. A recent survey Assaad et al. (2022) presents some state of the art methods to infer causality from time series. The main ones are based the conditional independence framework (Spitess et al.] (2000) which consists in running conditional independence tests to check whether variable X is a parent of another variable Y conditioned to some other variables Zi,.., 2n- Based on the PC algorithm Spirtes & Glymour (1991), in Runge et al. (2019). the authors introduced a new model (PCMCT) specific to causal discovery from times series, using what they call 'the momentary conditional independence (MCI)* to remove false positive returned by the PC algorithm. Some variants of PCMCI were then designed, like PCMCT+ Runge (2020) which includes contemporaneous links and improves
18

---

Published in Transactions on Machine Learning Research (06/2024)
0.20
1-AALasso
1-AALasso AALasso
AALasso
Prediction Error
6.10
Lasso
Lasso
SCORE
0.55
0.34
0.32
0.30
0.65
0.06
0.25
0.10
0.15
0.20
0.25
0.30
0.35
2.55
0.10
0.15
0.20
0.25
0.30
0.25
O.A.
0.08
1-AALasso
AALasso
Prediction Error
LS+ALasso
0.07
SCORE
1-AALass0
0.08
AALasso
0.70
LS
Lasso
0.65
0.60
0.15
0.20
0.35
5.05
0.10
0.20
0.25
0.50
A
A
0.055
Prediction Error
0.050
F1 SCORE
1-AALa550
0.040
Lasso
1-AALasso
LS+ALasso
Lasso
9.05
0.10
0.20
025
OAA
OAA
Figure 7: Prediction error and F1-score in function of the noise
(Top) N = 40 samples, (Middle)
N = 100 samples, (Bottom) N = 140 samples.
the reliability of tests, or LPCMCI Gerhardus & Runge 2020 which designed to address the issue of the presence of latent confounders.
Note that is possible include some prior knowledge when using methods based on conditional indepen- dence. These prior information can be expressed there is causal link j, there is not any link between and j, i is leaf, is root is ancestor of j. This kind of prior knowledge do not exactly match with the one used in this paper (the information is binary in this framework) but we tested PCMCI leveraging prior (denoled by Prior PCMCT) by encoding forbidding edges (i.,i) or (j.i) when A < (here we tested T E (0.2, 0.3)). All of these methods are implemented in hhtpps:/// jase://jakbbuuggg wwabbbuung..giuub..............e github. io/tigramiiiit t and we used this package perform the same experiments than the one Soction 4.4.1 We also add the results obtained with CGC-2SPR Yao et al. (2013), the method presented in Section 3.3 that leverage the same kind of prior knowledge than method but using L.2 penalty.
19

---

Published in Transactions on Machine Learning Research (06/2024)
A
0.1
0.25
0.35
Metrics
P
P
R
P
Lasso
0.28 (0.04)
0.61 (0.09)
0.38 (0.04)
0.28 (0.04)
0.61 (0.09) 0.38 (0.04)
0.28 (0.04)
0.61 (0.09)
0.38 (0.04)
AL
0.54 (0.09) 0.54 (0.11) 0.53 (0.07)
0.54 (0.09)
0.54 (0.11) 0.53 (0.07)
0.54 (0.09)
0.54 (0.11)
0.53 (0.07)
1-AALasso
0.44 (0.06) 0.73 (0.1)
0.55 (0.06)
0.39 (0.06)
0.71 (0.11) 0.49 (0.07)
0.37 (0.08)
0.7 (0.11)
0.47 (0.07)
AALasso
0.58 (0.07) 0.69 (0.11) 0.63 (0.08)
0.5 (0.07) 0.67 (0.12) 0.57 (0.07) 0.45 (0.06) 0.65 (0.12) 0.53 (0.07)
Table Precision, Recall and F1-score in function of 0A. We used N = 40 samples for training.
0.1
0.25
0.35
Metrics
P
R
F1
R.
F1
P
R
F1
Lasso
0.54 (0.13)
0.78 (0.08)
0.63 (0.1)
0.54 (0.13)
0.78 (0.08)
0.63 (0.1)
0.54 (0.13)
0.78 (0.08)
0.63 (0.1)
LS ALasso
0.67 (0.25)
0.77 (0.1)
0.69 (0.14)
0.67 (0.25)
0.77 (0.1) 0.69 (0.14)
0.67 (0.25)
0.77 (0.1)
0.69 (0.14)
Aasso
0.75 (0.05)
0.81 (0.08)
0.78 (0.06)
0.71 (0.06)
0.8 (0.08)
0.75 (0.06)
0.66 (0.07) 0.78 (0.1)
0.71 (0.07)
AALass
0.82 (0.06) 0.79 (0.09)
0.8 (0.07)
0.76 (0.06)
0.76 (0.1) 0.76 (0.07)
0.71 (0.09) 0.74 (0.11) 0.72 (0.09)
Table 5: Precision, Recall and F F1-score in function of σ We used N = 140 samples for training.
Samples
40
100
180
Metries
R
F1
F1
F1
Lasso
0.28 (0.04)
0.61 (0.09)
0.38 (0.04)
0.51 (0.1)
0.75 (0.08)
0.6 (0.09)
0.37 (0.08)
0.83 (0.06)
0.68 (0.07)
ALssO
0.54 (0.09)
0.54 (0.11)
0.53 (0.07)
0.69 (0.1)
0.68 (0.1)
0.68 (0.09)
0.64 (0.08)
0.83 (0.06)
0.72 (0.07)
CGC-2SPR
0.45 (0:22)
0.42 (0.16)
0.36 (0.07)
0.66 (0.22)
0.57 (0.11)
0.58 (0.10)
0.84 (0.12)
0.64 (0.10)
0.72 (0.06)
PCMCI+
0.89 (0.08)
0.38 (0.07)
0.53 (0.07)
0.92 (0.05)
0.56 (0.06)
0.69 (0.05)
0.92 (0.04)
0.67 (0.05)
0.77 (0.04)
LPCMCI
0.64 (0.08)
0.5 (0.06)
0.56 (0.06)
0.64 (0.07)
0.63 (0.06)
0.63 (0.06)
0.64 (0.08)
0.72 (0.06)
0.68 (0.06)
Prior PCMC++
0.68 (0.07)
0.5 (0.06)
0.58 (0.05)
0.7 (0.06)
0.64 (0.06)
0.66 (0.05)
0.71 (0.05)
0.74 (0.05)
0.72 (0.04)
1-AALass
0.44 (0.06)
0.73 (0.1)
0.55 (0.06)
0.63 (0.06)
0.82 (0.09)
0.72 (0.06)
0.83 (0.06)
0.84 (0.06)
0.83 (0.05)
AATasso
0.55 (0.07)
0.71 (0.11)
0.62 (0.07)
0.76 (0.06)
0.79 (0.1)
0.77 (0.07)
0.87 (0.05)
0.81 (0.07)
0.84 (0.05)
Table 6: Precision, Recall and Fl-score for all causal discovery methods tested, p = 40, σ = 0.1.
Compared to methods that do not leverage prior knowledge (Lasso, PCMCI and its variants), Table 6 show that AALasso allows for better network reconstruction, highlighting its ability to efficiently leverage prior information. Moreover, AALasso remains superior to Prior PCMCI+, showing that the way we introduce the prior information in our model is relevant and efficient. In details, we remark that PCMCI+ allows for a very good Precision for the retrieved causal links (from 0.89 to 0.92), but the Recall remains very low (from 0.38 to 0.67). This method achieves the best Precision but the worst Recall, leading to a Fl-score lower than the one of AALasso.
4.4.4 Empirical time complexity
Concerning time complexity, discussed theoretically in Section 3.5.2 we observed that the runtime of AALasso is approximately proportional to Niter times that of Lasso estimator. Therefore, it is useful to analyze the convergence rate of our method to find a balance between precision and computational efficiency. Figure [ shows the behavior of the Prediction error and F1-score averaged on all experiments conducted to visualize the convergence over the successive iterations. The convergence seems to be achieved for N = 10 in average (soe Appendix B.1 for more experiments), and that motivated the choice of Niter for comparing the performances of the algorithms. Then we can compare runtimes of Lasso, ALsso and AALasso with Nuer 10 iterations. To well understand the gain and the tradeoff between runtime and performances,
20

---

Published in Transactions on Machine Learning Research (06/2024)
plot the F1 score in function of the runtimes for Nuter {0,5, 10, with N 40 samples and σ 0.1
in Figure 9 and for other scenarios in Appendix B.5
We remark that AALasso with 5 iterations is a good trade off to optimize both the runtime and the f'I-score.
0.075
0.74
Lasso
Reconstruction Error
Reconstruction Error
LS
0.070
0.72
AALasso
0.065
0.70
Lasso
0,000
ALasso AALasso
0.055
0.050
10
17 16
iterations
N iterations
Figure 8: Prediction error and Fl-score in function of the number of iterations used for training using synthetic data with σд E (0.02, 0.1, 0.25, 0.35] and N E {40, 100, 140, 180, 250}.
0.65
0.60
F1 score
AALasso (1 iterations)
0.50
AALasso (2 iterations)
AALasso (5 iterations)
0.45
AALasso (10 iterations)
AALasso (15 iterations)
0.40
Lasso
LS+ ALasso
Runtime (in s)
Figure 9: Fl score in function of runtimes for Lasso, LS ALasso and AALasso for Niter {0,5, 10, 15) with
N = 40 samples, OA 0.1.
4.5 Experiments on real-world data
4.5.1 HPN DREAM8 Challenge
For this dataset, the choice of A prior is detailed in Section 4.3.2 and we compare the results with Lasso and Least Squares methods (note that we replace the ALaso method here because the results were better when fitting the Least Squares estimator).
The hyperparameters are selected using data for the cells BT549, MCF7, and UACC812 and the test set contains the data for all pairs (ligand, inhibitor) for the cell BT20. We use two common metrics compare the performances of the algorithms: the overall accuracy, measuring the relevance of variable selection, and
21

---

Published in Transactions on Machine Learning Research (06/2024)
0.7
0.6
0.6
ROCMJC
Accuracy
0.0
04
0.3
0.4
Leac Squares
ASSO
KALASSO
Least Squares.
LASSO
A6,A500
Mathod
Milthod
Figure 10: Results obtained on the DREAM8 dataset.
the Receiver Operating Characteristic Area Under the Curve (ROCAUC), allowing design threshold- independent mentioned in (Bradley] 1997] and measuring the ability of model to discriminate the two classes. Formally, let isijss E RP*P the parameters computed and (Yiz)ısussp {0, 1 the ground truth values (causality or not), then
Σ (,
I if |x] >
The accuracy is computed by Acc
I553SP
where T, x
taking
g
0 otherwise
T == 0.05. This metrics measures the relevance of variable selection.
The ROCAUC is computed by ROCAUC fo TPR(FPR- dt where TPR(t) stands for True Positive Rate (sensitivity), and FPR(() stands for False Positive Rate (1 specificiity with a threshold
As shown in Figure I0 our method AALasso demonstrates better performance when compared to traditional estimators such as Lasso and Ordinary Least Squares (OLS) regarding both accuracy and ROCAUC. Thus, by incorporating the adjacency matrix derived from the Pathway Commons database as a prior network, A.Lasso effectively harnessed valuable prior knowledge about protein interactions. Our results showcase that LAssso outperform the baseline methods, emphasizing the significance of considering such prior information in causality inference tasks. The Figure 11 presents example of graph learned, and remark that the Lasso estimator explain all the variables only with few number of variables (presence of columns in the matrix). On the contrary, the AALasso estimates are more homogeneous, following a directed version of the prior structure.
4.5.2 Molene Dataset
VAR(1) model For this dataset, we train the models with points, still selecting and y by cross- validation. Figure (12) compares the resulting graphs of Granger causalities using our method AALasso and Lasso. We observe that the graph resulting of the AALasso is sparse while remaining connected and allows good visualization of the physical process. Moreover, contrarily to the Lasso one, it is consistent with the Euclidean structure, confirming that the algorithm leverages the given prior matrix. This simple example encourages the using of AALasso to learn Granger causality for dynamic or physic system by taking into account the physics of the problem.
VAR(3) model For this example, we applied the generalization of our model with the order d We computed the same matrix prroo than for d (cf Section 4.3.3) and then considered the generalized
22

---

Published in Transactions on Machine Learning Research (06/2024)
0
10
20
10
20
30
(a)
0-
10
20
10
10
20
30
(c)
(d)
Figure 11: Example of graphs learned. (a) Least Squares method. (b) Lasso method. (c) AALasso. (d) Prior network.
prior matrix given by
A prior
A prior
A prior
1
0
Arior
0
0
The comparison between the ALasoo and the AALasso estimates for d = 3 is presented if Figure 13 Like for the case d = 1, the AALasso graph is sparser than the LS+ALasso one, and allows better visualization of the physical process. Moreover, using the first order (only 3 and 1 edges for order 2 and 3) which is consistent with diffusion process. Finally, the size of the edges for AALasso seems to be proportional to the order (edges for order 2 and 3 are longer than the ones for order 1) which is again consistent with the physics (information take more time to travel longer distances).
remark that AATasso explains the main part of the signal only by
5 Discussion and conclusion
In conclusion, this paper has introduced novel method designed to efficiently learn Granger causalities settings with limited samples. Our approach stands out by effectively incorporating prior knowledge through the utilization of noisy adjacency matrix. We demonstrated the convergence of our algorithm AALase and
23

---

Published in Transactions on Machine Learning Research (06/2024)
Laniin
Mo
SHi
Brieuc
R
Landernea
Landeria
Dowimente
Boulmmet
Quimpx
Ponthor
Qaima
Prip-Sy
Concameeuu
Concatioon
Canert
Plarmeg!
Vorness
Guérande
Saint-
Gueran
Bauto
La Bault E'coutLA
Scoublac
AALaso
(b)
Figure 12: C Results on the Molène dataset for (a) AALasso and (b) Lasso. Darker colors indicate larger weight.
Order
Order 2 (a) LS ALasso
Order 3
Order 1
Order 2 (b) AALasso
Order 3
Figure 13: C Results on the Molène dataset. with for (a) L. + ALasso and (b) AALasso. Darker colors indicate larger weight.
showed that the time complexity was in the same order of magnitude as that of the Lasso. To empirically validate our method and demonstrate its efficacy, we conducted series of experiments. We selected variety of datasets, including synthetic data and real-world examples like the Breast Cancer wotwork and the Molène
24

---

Published in Transactions on Machine Learning Research (06/2024)
Dataset. In these experiments, we employed rigorous evaluation metrics to assess the performance of our method across different scenarios, showcasing its versatility and applicability. Thus the incorporation of prior information has proven instrumental in achieving superior accuracy and robustness when compared to classical algorithms in this domain.
While our method allows to incorporate prior knowledge in the learning process, it could be interesting to add structure the learned graph and go beyond sparsity. Indoed, the framework we have presented here can be readily extended to learn graphs with specific structural constraints, such spectral and adjacency constraints, similar to those outlined in (Kumar et al. 2020). This flexibility arises from the fact that our model operates with symmetric matrix containing positive values, making it amenable to range of applications that necessitate tailored graph structures.
Finally, it could be interesting lo study whether this framework could be used to learn lime-varying graphs of Granger causality (ef (Gao & Yang 2022 ), for example by considering the graph learned at time as a prior for the graph at time 1.
References
Charles K Assaad, Emilie Devijver, and Eric Gaussier. Survey and evaluation of causal discovery methods for time series. Journal of Artificial Intelligence Research, 73:767-819, 2022.
Giovanni Ballarin. Ridge regularized estimation of var models for inference. arxie, 06 2022. URL https: //doi.org/10 10. 48550/arXiv. 2105.00860
Marta Banbura, Domenico Giannone, and Lucrezia Reichlin. Large bayesian vector auto regressions. Jowrnal of Applied Econometrics, 25(1):71-92, 2010. URL https:///cconappert repeaeaaaaeeeeeee.oooo.e.ree.. repec. jaa japmet:
v:25:y:2010:1:1:pp::
Sumanta Basu, Ali Shojaie, and George Michailidis. Network granger causality with inherent grouping structure. Journal of Machine Learning Research, 16(13):417-453, 2015. URL http:///mml..oor org/ //aperr
s/v16/basu16a.htt html
Stephen Boyd. Distributed optimization and statistical learning via the alternating direction method of multúpliers. Foundalions and Trendso in Machine Leurning, 3(1):1-12) 2010. ISSN 1935-8245. doi:
URL http://dx. doi. org.10.0100.. 1561/2200000016
Andrew P. Bradley. The use of the area under the roc curve in the evaluation of machine learning algorithms. Pattern Recognition, 30(7):1145 1159, 1997. ISSN 0031-3203. doi: http:://bbbi.rr/110111////// 3(96)001..-2. URL http:://wwwaaoo/oooc////////////
Daniel E. Carlin, Evan O. Paull, Kiley Graim, Christopher K. Wong, Adrian Bivol, Peter Ryabinin, Kyle Ellrott, Artem Sokolov, and Joshua M. Stuart. Prophetic granger causality to infer gene regulatory networks. PLOS ONE, 12(12):1 21, 12 2017. doi: 10.1371 //oournal..oo..011000.. URL httpps://oi...... g/10. 10t.117//jouunall.mmmm 0077000030000
Ethan G Cerami, Benjamin F. Gross, Emnk Demir, Tgor Rodchenkov, Özgün Babur, Nadia Anwar, Nikolaus Schultz, Gary D Bader, and Chris Sander. Pathway commons, web resource for biological pathway data. Nucleic acids research, (suuppl 1):D685 D690, 2010.
Steven Diamond and Stephen P. Boyd. Cvxpy: A python-embedddd modeling language for convex optimiza- tion. Journal of machine learning research JMLR, 17, 2016. URL sttp:://ppi.sematticcchoo (CorppssID:629800)
Leo L Duan, Zeyu Yuwen, George Michailidis, and Zhengwu Zhang. Low tree-rank bayesian vector autore- gression models. Journal of Machine Learning Research, 24(286):1 35, 2023. URL http://jmlr.org//
apers/v24/22-0366 html
Bradlcy Efron, Trevor Hastie, Iain Johnstone, and Robert Tibshirani. Least angle regression. The Annals of Statistics, 32(2):407~ 499, 2004. doi: 10. 10.1214400000051 1214/0000533000000000 URL httpp:://dot wtt.oi..or/........ 44/09055364000000087
25

---

Published in Transactions on Machine Learning Research (06/2024)
Wei Gao and Haizhong Yang. Time-varying group lasso granger causality graph for high dimensional dynamic system. Paltern Recognilion, 130:108789 2022. ISSN 0031-3203. doi: http::///doo.rrg//.000///..aaaa 022.108789.
Andreas Gerhardus and Jakob Runge. High-recall causal discovery for autocorrelated time series with latent confounders. In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin (eds.), Advances in
Newral Information Processing Systems, volume 33, pp. 12615 1225.. Curran Associates, Inc., 2020. URL http:://prooeeiiggg..uurr..ppepp www.:ng....nnrii..coooooomo. 28d0b-Paper.po pdf
Satyajit Ghosh, Kshitij Khare, and George Michailidis. High-dhmensional postcrior consistency in bayesian vector autoregressivve models. Journal of the American Statistical Association, 114:735 748, 2018, doi:
10.1080/0111145..2001.131333
C. W. J. Granger. Investigating causal relations by econometric models and cross-spectral methods. Econo- metrica, 37(3):424 438, 1969. doi: http://di..rr/1002377///999/
L. Grippo and M. Sciandrone. On the convergence of the block nonlinear gauss seidel method under convex constraints. Operations Research Letters, 26(3):127-136 2000. ISSN 0167-6377. doi: https://doi.orrg/l .1016/ (001//$0667-6377()0007777 URL http:///pppp/.aaacaaccccccccoooa
37799000747
James Douglas Hamilton. Time series analysis. Princeton University Press, Princeton, NJ, January 1994. ISBN 9780691042893.
Vân Anh Huynh- Thu, Alexandre Irrthum, Louis Wehenkel, and Pierre Geurts. Inferring regulatory networks from expression data using tree-based methods. PLOS ONE, 2010. URL https://journals.. .loos. org/
plosone/articleeiee-l 1371/.2Fourrnal pone. 0012776
Elvin Isufi, Andreas Loukas, Nathanaël Perraudin, and Geert Leus. Forecasting time series with varma recursions on graphs. IEEE Transactions on Signal Processing, 67(18):4870-4885 2019. doi: 10.1109/TS
P.2019.2929930.
Anders Kock and Laurent Callot, Oracle inequalities for high dimensional vector autoregressions. Journal of Econometrics, 322-3343553344, 2015. URL httpps//EcooPaperss repec. wwteeeegrg//eeeeeeeooo.oo :ee::PEeeee::::ceoooo.
186.y:20166:::::::325---11
Sandeep Kumar, Jiaxi Ying, JosA@ VinAcius de M. Cardoso, and Daniel P. Palomar. A unilied framework for structured graph learning via spectral constraints. Journal of Machine Learning Research, 21(22):1-60 60,
2020. URL httpp:/jmmr. http:p////all.ooggpppars///////s https/ppprs//v///122 121m11127..hmm
Jiahe Lin and George Michailidis. Regularized estimation and testing for high-dimensional multi-block vector-autooregressive models. Journal of Machine Learning Research, 18(117):1 49, 2017. URL http:
w/jilr.org/paperss//e ott//appes///1/11000555 html
Jiahe Lin, Huitian Lei, and George Michailidis. Structural discovery with partial ordering information for time-dependent data with convergence guarantees. Journal of Computational and Graphical Statistics, pp. 1-20, 01 2024. doi: 10.1080//0611600.220 0t0.:/88600.02228811099771 2301097.
Helmut Lütkepohl. New Introduction to Multiple Time Series Analysis. Springer, 2005. CRL https:
Benjamin Marlin, Mark Schmidt, and Kevin Murphy. Group Sparse Priors for Covariance Estimation. arXiv e-prints, art. Nr.Xiv:1205.2626 May 2012. doi: 10.48550/arXiv... 1ww................
Jakob W. Messner and Pierre Pinson. Online adaptive lasso estimation vector autoregressive models for high dimensional wind power forecasting. International Journal of Forecasting, 35(4):1485-1499) 2019. doi: 10.001//.iiforrcas..11...2......
26

---

Published in Transactions on Machine Learning Research (06/2024)
Silvia Miranda Agrippino and Giovanni Ricco. Bayesian vector autoregressions. working paper or preprint, May 2018. URL htttps///sciencespoo hal. sciemccehhal-00455277
William B. Nicholson, Ines Wilms, Jacob Bien, and David S. Matteson. High dimensional forecasting via interpretable vector autoregression. Journal of Machine Learning Research, 21(166):1 -2, 2020. URL
http://pa...oo//pppppps/////1111 .tmll
Antonio Orlega, Pascal Frossard, Jelena Kovačević, José M. F. Moura, and Pierre Vandergheynst. Graph signal processing: Overview, challenges, and applications. Proceedings of the IEEE, 106(5):808 828, 2018. doi: J0.1100//PROC
Bastien Pasdeloup, Vincent Gripon, Grégoire Mercier, Dominique Pastor, and Michael G. Rabbat. Char- acterization and inference of graph diffusion processes from observations of stationary signals. IEEE Transactions on Signal and Information Processing over Networks, 4:481-496, 2016. URL https:
Jakob Runge. Discovering contemporancous and lagged causal relations in autocorrelated nonlincar time series datasets. In Conference on Uncerlainly in Artificial Intelligence, pp. 1388 1397. Pmir, 2020.
Jakob Runge, Peer Nowack, Marlene Kretschmer, Seth Flaxman, and Dino Sejdinovic. Detecting and quan- tifying causal associations large nonlinear time series datasets. Science advances, 5(1 ;(11))9a,,4996, 2019.
Anil K. Seth, Adam B. Barrett, and Lionel Barnett. Granger causality analysis in neuroscience and neu- roimaging. Journal of Neuroscience, 3598):3293 3297, 2015. ISSN 0270-6474. doi: 10.1523/JNEUROSC. 4399-14.2015. URL
Ali Shojaie and George Michailidis. Discovering graphical Granger causality using the truncating lasso penalty. Bioinformatics, 26(11):1517- 09 2010. ISSN 1367-4803. doi: 10.1093/ bioinformatics/bttee bt7777 URL http::///aabll.g/.....0/0//aaoo.
Christopher A. Sims. A Nine- Voriable Probabilistic Macroecomomiic Forecasting Model, pp. 179 212. Uni- versity of Chicago Press, January 1993. URL http://www wtt.//wwwwwbb....eeo wtggchaattrs//7192
Peter Spirtes and Clark Glymour. An algorithm for fast recovery of sparse causal graphs. Social science computer review, 9(1):62-72, 1991.
Peter Spirtes, Clark N Glymour, and Richard Scheines. Causation, prediction, and search. MIT press, 2000.
J.H. Stock and W. Watson. Dynamic Factor Models, Factor- Augioo-ed Vector Autoregressions, and Structural Vector Autoregrossions in Macrocconomics. In J. B. Taylor and Harald Uhlig (eds.), Handbook of Macrorconomicss volume 2 of Handbook of Macroeconomics, chapter 0, pp. 415 525. Elsevier, 2016. doi: http//.hhsma..206.... URL httts:///dda:: repec. httpppppc..gog///e/////cccccc//p 41m.hmml
Robert Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society: 10000............... URL http://rrra.oomneieeeerrrrrraa/ http::/ww.wlllaaao////////////// 2tt.::/....11111111 2517 tb00000
Series B (Methodoolopical), 58(1):267 288, 1996. doi:
onl
P. Tseng. Convergence of a block coordinate descent method for nondifferentiable minimization. Journal of Optimization Theory and Applications,, 109(3):475 494, jun 2001. doi: 10.1028/0:10133017011006 URL
https://doi.. 1t.rg/11..000 http://00.0000//////011110000000
Tong Tong Wu and Kenneth Lange. Coordinate descent algorithms for lasso penalized regression. The
Annals of Applied Statistics, 2(1), March 2008. ISSN 1932-6157. doi: 10.1214/07-a0ss447 URL http:
wax.doo..rgg/0... 1214111-0007071144
Shun Yao, Shinjac Yoo, and Dantong Yu. Prior knowledge driven causality analysis in gene regulatory network discovery. In 2013 IEEE 13th International Conference on Data Mining Workshops, pp. 124-130, 2013. doi: J0.100//CDMM.2203.100.
27

---

Published in Transactions on Machine Learning Research (06/2024)
Yujie Zhao and Xiaoming Huo. A survey of mumerical algorithms that solve the lasso problems. arXiv, 2023. URL httpt::ttt.....ooo v www.aa.ww............a.........a 03576
Cunlu Zou, Christophe Ladroue, Shuixia Guo, and Jianfeng Feng. Identifying interactions in the time and frequency domains in local and global networks granger causality approach. BMC Bioinformaties, 11 (1), 2010. doi:
11-337.
Hui Zou. The adaptive lasso and its oracle properties. Journal of the American Statistioal Association, 101 (476))1418-1429, 2006. doi: 10.119/000011100000007733 URL https::/abb.oog///0.111///0/22 6000000735
Álvaro Méndez Civieta, M. Carmen Aguilera-Morillo and Rosa E. Lillo. Asgl: A python package for
penalized linear and quantile regression. arXiv, 2021. URL attps://doi. org/ 10. 48550/arXiv 2303. 03 576
28

---

Published in Transactions on Machine Learning Research (06/2024)
A Proofs
A.1 Statistical model
Proof of 13
The MAP of the model
4 is given by:
A, C arg max I(A, C (X1(1 [[[[1]]]]=1 A A
(17)
A,C
subject to 2
(18)
where L(-) is the likelihood function. Using Bayes formula, ones have
L(A", C* X A A = PAA", C* | A A
=PC(C" X A prior
I(A" C X I:t A prior)
=IP(C" | A prior
x IPAA"
C A prior)
P(X IC') X P(C
prior
X P(A" C", A",AArioo
P(X Aprior)
P(XL# C*)
x (P(C* Aprior)
x P(A" C", Aprior)
P(X1# Aprior)
Then:
IP(C" Aprior) x P(A" C", A Aprior)
=P(C" AA", A Arior) x IPAA A prior)
=P(C" |A*) x P(A" A Arrior)
P(A")
=P(C* A*) x P(A prinr A')
(A P(Aprior) prior)
Assuming that A follows the improper distribution (R), the MAP is finally given by:
max P(X1: A'CC*
max C*) IPCC | (A))
PrA A*)
N(X'; ox Id) Il
Laplace (7 0, A7.) x
II
2
Applying the log function the previous expression concludes the proof.
A.2 Proof of Theorem 12
Definition 14 (Regular function). We say that f is regular at domf if f'(z; 2 0, vd (d), dp), such that f2(2; (0, ..., dg,..., 0)) 2 = I....p. Remark 15. If f differentiable then ('(x; = d. So, if f' ..., dr .))...,))))))))) 2 0 4k. 1,.., P. we have that:
P (z1;; (di,..,, Vfxx) (di..., de...,, dp) Σ (0.,) 20 20
(19)
Thus a differentiable function is regular.
Definition 16 (Stationary points). We say that z is a stationary point of f if 2 E domf and f((2; Z
0, Vd.
We say that z is coordinatewise minimum point of f if domf and f(z (0, ..., du ... 0)) 2 f(z),
Yd, E R"* for all I 1,.,,, N.
Remark 17. If z is a coordinatewise minimum point of f, z is a stationary point of f whenever f regular
al z.
29

---

Published in Transactions on Machine Learning Research (06/2024)
Using this framework, the following theorem holds.
Theorem 18 (Theorem 4.1 in (Tseng 2001)). Assume that the level set Lx = fz f(x) S f(zo)} is compact (where To E IR and that f is conlinuous on Lx. Then, the sequence {x" {=")r=1..... generaled by the Block Coordinate Descent method is defined and bounded. Moreover, the following statements hold:
If [(x1,.., x())) has at most one minimum in Ik for l E (2..., - 1}, then every cluster point a of {x") Ir=N-1 mod is a coordinatewise minimum point of f. In addition, if f is regular 2, then is stationary point of f-
Actually the case 2 is very simple and we do not need any supplementary assumptions that the continuity
of f and the compactness of Lx to prove the convergence of the algorithm. Indeed, since f2...., N - 1} = for the case N = 2, the point (3) can be directly applied without satisfying any assumptions about the number of minimum so ones obtain directly the following corollary.
Corollary 19. Assume that the level set Lx is compact, that f is continuous Lx and that N = Then, the sequence {x"), ,==1,.... generated by the BCD method is defined and bounded. Moreover, every cluster point z of {x') [r")r=N-1 mod N is a coordinatewise minimum point of f. In addition, if f is regular at 2, then 2 is stationary point of f-
A simplified version of the proof provided in (Teeng 2001 for the case N = 2 can be found in Appendix A.2h
Now we need to prove that the function f of our model (the MAP) satisfies the conditions of the previous theorem.
Proposition 20. The function f: TRP² x TR. R defined by
1 N
N | Cx |xx"
Ca.i + 1Cril +1 H
f(C,A) =
log(2A.,)) + A (20)
Aij
is regular.
Proposition 21. The function I defined in 20 but constrained on IRP" x [e, xo0= with 0 satisfies
assumptions in 19
Finally, using 20 and 21 show that objective function in 7 verifies the assumptions in
19 and
applying the Theorem, we obtain the following result.
Proof 19 Let (x") 1,=1,,,. the sequence generated with the BCD algorithm. By definition of the algorithm ones have that f1 & f(x") for all r and x1+1 E Ex for all r. Since Cx is compact, {z JreR converges towards = == zt. In the same way, we can assume w.loo.g that {*r+1) JPer converges towards z² (taking a sub-sequence). First, [f(x"))EER is decreasing (and bounded) so it converges and ones have that
f( lim
f(x') f())= f(2²).
FER +
Now, we assume that for every r E R, r" argmin ), egg for all
f(x") & f(x" + (d), Adi f(x 1 < f(x 11 + (0, di)), Adz
of + where x" = (2)
30

---

Published in Transactions on Machine Learning Research (06/2024)
Since f is continuous on Lx, get
f(z!) S f(=' (di, hpA
ʃ(x²) s f(=2 + (0,d2), vd2
zî =1
Then, for all d2
f(=') = l(=
[((=), + + (0,d2))
[((=}, =5) + (0,d2))
= =[((x,,)) + (0, =)) + (0,d2)
f(='' + (0,d2))
Since 22 2, we proved that for all di, d2
f(z) S f(z (d1, 4di
[(z) < [(= + (0, d2)), Vdz
is componentwise minimum of f.
Finally, if f is regular, the previous incqualities become
(z1 2 0, 4d1
(0, 2 0, Vd2
and by definitions z is a stationary point of f-
Proof of 20 The only points where f is not differentiable are () Bi,j; Ct.j 0) because of the ab-
solute value. Let's write::
S(C,A) (CC,)) + h (C,A)+
+
where
(C, 1 N
Σ
[ [t] Cx 1]1 1
N n=1
hc,(C, |Ci,l |yjal
As,i
I(C, = A log(2Ai.)) + 7/1A A
9 and are differentiable so we have for all (C,A) E Ry x " for all (Dc, DA) E TRP² RR² such that
(C+ Do, A DA) E domf:
g' DA)) = [ .0)) + (O
31

---

Published in Transactions on Machine Learning Research (06/2024)
I((CC, (Dc,, E (D((CA))(()) 1.
( .0) (); D
where D(.)) is the matrix with zero values everywhere except the coefficient (i,j) which is equal to D) If C.. # #0 and C5 # 0, then has is differentiable in (C,A) so we have the same result. Otherwise, we need compute the lower directional derivative of his (C, W) with Ct,- = 0 or C₁ 0.
Ones can compute that:
(D. DA)) [| if Ci.j 0 and Cj.
A.,
A),((C, DA)) the same for the last case.
(ign(C,,)Do
if Ci. 0 and C.j & 0 and we can do D
A.,
An,
A
Thus we still have that
(C,(CC.) E 1
((C.A)) .0) (CC,
Finally, by definition of regular function, ones have that f is regular at all (C, E Rd x R.
Proof of 21 First it is clear that f is continuous on LX.
Let's s consider again the decomposition
() - A(C,A)= 9(C) +C,(()))+((A))
1.j
where
(CC) 1
N Σ | cx -1]1
Ci,ll 1Ctss
his (C,(C
Aij
((C,,,,
Σ i<j
log(2A,.)) + 7/ A 1
It is clear that lim ((C) = too and ||0|-++o
|| lim I(A) = +00.
+
Then, since bigj (C, 2 for all C, A for all i, we have that.
lim J(C, A) = +00. [[C.A)++O
(21)
We proved that f was coercive, it follows that Lx is bounded. Moreover, lim I(A) = +oc, for i,j E [[], Additionally, f is continuous so f-1( - 0o, ((x(0))) closed A and finally Lx is compact.
Proposition 22. The function defined in 20 is not conver.
Proof of 22 The function 20 is a function of (C, (A4)) so is sufficient to prove that it is not convex in jij for fixed and j and for fixed value of C. The second derivative wrt Ai, is
X
Cr.l
1C
+
27
(22)
A2
A3
jj
32

---

Published in Transactions on Machine Learning Research (06/2024)
it has the same sign than the degree 3 polynomial:
-AAss + 2(1ccsl++ + 2
(23)
Then, the minimum of Chis polynomial on [0, + reached in v by and take the value C3al) which can be negative for small values of |Ca,Il + [Caal- Thus the second derivative can reach negative value with certain value of Ci. and Cj.--
223/2 + 2([Ct,]
3/07
7
A.3 Time complexity
Lemma 23 (C update time complexity). Let Casso(p, the time cost for training Lasso estimator to fit VAR(1) parameters in dimension with N samples, then the time complexity of C update is in (p x N + Ctassolp,N).
Proof of 233 The C' update is done solving an Adaptive Lasso problem, so recalling that an Adaptive Lasso problem can be writlen as a Lasso problem in x N multiplications and that we solve Adaptive Lasso problems (cf 3.4.1 the time complexity of this step is in (p² x N CLasso (p, (p, N)).
Lemma 24 (A update time complexity). The time complexily of A updale is in O (p²).
Proof of 24 The A update is done by computing in O(1) the closed form given in 3.4.2 for each value Aij
so this update is in O(p²).
In order to completely express the time complexity, we need to compute CLasso .N.. Note that when using gradient descent based methods to solve an optimization problem, whereas a convergence rate analysis can be conducted (cf (Zhao & Huo 2023]), the time complexity depends on the stop criterion of the algorithm. Thus we conduct the time complexity analysis assuming that ADMM utilized with a fixed number of iterations NADMM.
Lemma 25 (ADMM complexity for Lasso). The time complerity of the ADMM algorithm (with fized number of steps) to solve one Lasso problem in dimension is O(p² + p²).
Proof of 25 The updated formula of the ADMM given in 3.4.1 require to multiply N mairix which is in O(N x p²) and to inverse p matrix which is in O
N matrix with
Theorem 26. Let CAALasso (P, the time cost for truining a AALasso estimator to fit VAR(1) parameters in dimension p with N samples, then:
CAALaSSo(,, N) P.N
OCCrassopp N)) = O x N +
(24)
P.N
Proof of [ Summing the time complexity of A step 24 and C step 23 gives complexity in O (p² + N (p. Since the matrix inversion need to be performed only time for the Lasso problems at each step, the complexity of (p, N) using 25 becomes O(p² + N x p²), finally
resulting in complexity in O (p² +p² N
+ = O (p² x N x = OCCasso N)).
B Additional experiments
B.1 Number of iterations
For synthetic experiments in Section f we motivated the choice of Nute = 10 by Figure S While this parameter allows good results in various scenarios, it can interesting to understand whether this parameter is related the datasot. Convergence analysis results are shown in Figures [14 and [5] trend seems appear the larger N, the faster the convergence. Moreover, while the noise impacts the performances of A ALasso, the convergence rate seems to remain unchanged for a 0.1 a 0.25.
33

---

Published in Transactions on Machine Learning Research (06/2024)
Lasso LS
0.20
Lasso
Reconstruction Error
LS + ALasso
AALasso
0.18
AALasso
F1 Score
0.50
0.12
6.45
0.10
0.0
5.0
10.0
12.5
15.0
17.5
0.0
2.5
5.0
7.5
10.0
12.
15.0
17.5
N iterations
iterations
N
40,
0.1
Lasso
0.08
LS
Reconstructio Error
0.75
0.07
F1 Score
Reconstruction
0.70
Lasso
+ALasoo
AALasso
0.60
0.0
5.0
15.0
17.5
7.5
12.5
N iterations
N iterations
(b) N = 40, 0 = 0.25
Lasso LS ALasso AALasso
0.50
Reconstruction Error
F1 Score
0.75
Lasso
LS
0.70
AALasso
0.9
5.0
10.0
12.5
15.0
17.5
0.6
2.5
5.0
7.5
10.0
125
15.9
17.5
N iterations
N iterations
(c) N 100, σ = 0.1
Figure 14: F1-score and Prediction error through AALasso iterations for N {40, 100, 140) and σ 0.1.
B.2
Runtime
B.3
Initialization impact
Let's recall the algorithm to train AALasso
The results presented in the previous sections were obtained initializing (0)) with Arroor However, it can be interesting to test the algorithm with other initializations. We therefore conducted experiments initializing ( with only 1 values (denoted AALasso-ones! A obtained by solving the Least Squares problem (denotes A ALasso-T.S) or with random values (denoted AA AATssso-random). Our method is proven to converge towards
34

---

Published in Transactions on Machine Learning Research (06/2024)
Lasso
0.20
LS ALasso
Reconstruction Error
AALasso
0.18
F1 Score
Lasso
LS5ALa
AALasso
0.40
0.0
25
5.0
7.5
10.0
12.5
15.0
17.5
0.0
2.5
5.0
7.5
10.0
12.5
15.0
17.5
iterations
N iterations
(a) N = 40, σ
0.1
Lasso LS=ALasso
Reconstruction Error
0.725
AALasso
F1 Score
0.700
0.07
Lasso
0.675
LS+ALasso
0.550
0.06
AALasso
0.05
5.0
10.0
15.0
17.5
10.0
17.5
N iterations
N iterations
(b) N 40, 0
0.25
Lasso
LS+ALasso AALasso
Reconstruction Error
6.050
F1 Score
0.725
6.045
0.700
0.040
ALssoo
AAlasso
6.035
0.850
0.625
0.030
0.025
2.5
5.0
7.5
10.0
12.5
15.0
17.5
2.5
5.0
10.0 N iiterations
7.5
12.5
15.0
17.5
N iterations
(c) N 100,
Figure 15: F1-score and Prediction error through AALasso iterations for N E (40, 100, 140) and σ 0.25.
stationary point, no matter the initialization of A However, since our objective function is not convex, better local minimum could be found starting from 1 vector random vector rather than the A prior The results for all the scenarios tested are presented in Figures 16 and 17 Globally, the results are similar for all initializations tested and it is not clear if one of the initializations is better than the other. However we remark than in settings with very few samples (N = 40), the random initialization surprisingly outperforms the other with a gain of 2.5 in the Fl-scores compared to AALasso for all noise levels.
35

---

Published in Transactions on Machine Learning Research (06/2024)
Algorithm 3: Fitting algorithm.
input :Nite:: X, 7,
output: C, W
W()) Subdiagonal values of A
for i <-1 to Nater do
C( +- fo(C, W(i=l)) where fc denotes the update in W + fw( (C where fw denotes the update in 3.4.2 return ((Nrw)) ((Nim)))
3.4.
0.8
0.8
F1 SCORE
F1 SCORE
0.7
0.7
AALasso-LS
AALasso-LS
AALasso-random
AALasso-random
AALasso-ones
AALasso-
0.5
0.5
AALasso
AALasso
LS+ALasso
LS+ALasso
100
150
250
100
200
250
Samples
Samples
0.02
0
0.1
0.9
0.8
F1 SCORE
0.7
F1 SCORE
0.7
0.8
AALasso-LS
AALasso-LS
AALasso-random
AALasso-random
AALasso-ones
AALasso-ones
AALasso
AALasso
LS+ALasoo
LS+ALasso
Lasso
Lasso
100
250
100
250
Samples
Samples
T 0.25
A 0.35
Figure 16: Impact of the initialization of A (0) on the Fi-scores with σ E (0.02, 0.1, 0.25, 0.35).
B.4 Comparison with random prior graph
In this last experiment, we compare our results with the AALasso method taking a random A prior to check that the prior structure is well leveraged and that random L2 penalization can not achicve the same performances. The values prior sampled independently from uniform distribution in [0.2, 1). The results are presented Figure 18 and we see that results using a random A prior are very poor regarding both Fi-score and Prediction error. This reinforces the thesis that AALasso judiciously leverages the information provided by A prioo
36

---

Published in Transactions on Machine Learning Research (06/2024)
0.80
0.65
0.75
F1 SCORE
F1 SCORE
0.55
0.70
0.50
AALasso-LS
AALasso-LS
AALasso-random
0.03
AALasso-random
AALasso-ones
AALasso-ones
AALasso
AALasso
6.40
L5+ALasso
0.60
LS+ALasso
Lasso
Lasso
0,15
0.20
0.35
0.05
0.10
0.15
0.20
0.25
0.35
Noise
Noise
40
N 100
0.90
0.85
w
F1 SCORE
0.80
AALasso-LS
SCORE
0.80
AALassorrandom
AALasso-ones
0.75
AALasso
AALasso-LS
LS+ALasso
F1
AALasso-random
0.65
AALasso-ones
0.70
AALasso
AL
Lasso
0.65
0.05
0.10
0.15
0.20
0.25
0.30
0.35
0.05
0.10
0.15
0.20
0.25
0.35
Noise
Noise
N 180
N
250
Figure 17: Impact of the initialization of A on the F1-scores with N E (40, 100, 180, 250).
AALasso
AALasso random prior
Prediction Error
ALasso
Lasso
F1 SCORE
0.20
0.15
AALasso
AALasso random prior
LS+ ALasso
Lasso
100
200
250
200
Samples
Samples
Figure 18:
Comparison of AALasso, Lasso and LS+ALasso with AALasso taking prior random, with
a 0.1.
B.5 Runtime
Here, several scenarios tested to complete the Figure and show that the behavior regarding the Fl score in function of the runtime remains similar for some parameters choices.
37

---

Published in Transactions on Machine Learning Research (06/2024)
0.65
0.80
0.60
0.75
F1 score
F1 score
AALasso (1 iterations)
AALasso (1 iterations)
0.50
(2 iterations)
AALasso iterations)
iterations)
AALasso iterations)
AALasso iterations)
AALasso (10 iterations)
AALasso (15 iterations)
AALasso (15 iterations)
0.46
Lasso
Lasso
LS + ALasso
LS+ ALasso
Runtime (in s)
Runtime (in s)
(a)
0.60
0.55
0.725
F1 score
F1 score
0:50
AALasso (1 iterations)
0.675
AALasso Iterations)
AALasso (2 iterations)
0.650
AALasso ittrations)
0.45
AALasso (5 iterations)
AALasso iterations)
AALasso (10 iterations)
(10 iterations)
AALasso (15 iterations)
0.600
AALasso (15 ittrationn)
Lasso
Lasso
LS
ALasso
Runtime (in s)
Runtime (in s)
(d)
Figure 19: Fl score in function of runtimes for Lasso, LS+ALass and A.ssss for Niter (0,5, 15.. (a) N 40 samples, σл 0.1. (b) 100 samples, σл 0.1. (c) N 40 samples, σ 0.25. (d) N 100 samples, σ 0.25.
38


Reviews:
[{"review_id": "VaQrv6fiBy", "reviewer": "Reviewer_Hjee", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n1) The paper provides an extensive mathematical study of the proposed two stage optimization method, and importantly convergence guarantees.\n\n2) The loss function of equation (7) of the main paper, is well-motivated and principled under the statistical model provided in Figure (1) of the main paper.\n\n3) The model is simple and intuitive.\n\n\nWeaknesses:\n\n1) In real-world scenarions, the problem setting where there exists even a noisy A_prior matrix is rather limiting. In most settings, such an assumption is difficult to be made. In this sense, the proposed approach is unable to work in the more general case where prior knowledge is not present.\n\n2) The baselines are rather limited. Comparison with more recent and competitive methods that aim on learning Granger causality graphs based on time series data should be provided. I am including just a few here:\n\ni) Amortized Causal Discovery: Learning to Infer Causal Graphs from Time-Series Data, Löwe et al, 2022\n\nii) Discovering Nonlinear Relations with Minimum Predictive Information Regularization, Wu et al, 2020", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper introduces a novel optimization problem that allows to estimate VAR models while incorporating the use of Granger causality graphs based on the assumption that there exists a noisy Prior causal graph of the specific multivariate times-series input. Importantly, the authors provide an efficient two-block coordinate descent algorithm while also providing guarantees for the convergence of such an algorithm. The results showcase the superiority of the method against the considered baselines. Lastly, they showcase the important of sparsity in the performance of the general model.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1) Please provide an extensive discussion on the problem setting, and on how common is the existence of a prior Granger causality information\n\n2) The VAR coefficient matrix C as stated in the main paper expresses the directed Granger causality graph. Given that, why is A and Aprior chosen to be undirected? Would it make more sense that also the prior is directed? Please elaborate more on this modeling decision.\n\n3) Please include more recent and complex baselines and compare how your method competes against them in discovering causality graphs.", "broader_impact_concerns": "No ethical concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "7TsRDrjbpc", "reviewer": "Reviewer_oNFM", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\nThe paper is well written, with clear arguments and proofs.\nThe model is interpretable which is very important for practitioners in the targeted fields.\nThe algorithm is sound and the convergence proof for a non-convex and non-differentiable function is a useful result.\nThe method is tested in both simulated and real data and it improves compared to other traditional methods.\nA visual explanation of the results facilitates the qualitative evaluation of the methods.\n\nWeaknesses:\nTo my understanding, the method aims to reconstruct the causal network. This is useful in cases where we search for new causal relationships while observing an already established causal network i.e. enriching the network. If this is the aim of the method, it is not perfectly clear from the title or the text and the experiments.\nSpecifically, the simulated data utilize a noisy representation as a prior and are evaluated based on the reconstruction error with the correct network. This is not very clear with the real data, as neither the noisy network nor the reconstruction error (accuracy in 4.5.1?) is clarified sufficiently.  Similar to section 4.1. an integral experiment is to remove edges for varying ratios in the real network, use it as a prior, and calculate how close is the retrieved network with the real one. \n\nAn even more important question is the role of A in eq. 7. If C is the retrieved network and A^prior the noisy observation, what does A represent? I imagine it can not be the actual causal network that we aim to identify with C, because in the real-world this is not observed.\n\nIn addition, the method is not tested sufficiently. Comparison with other causal network discovery models https://arxiv.org/pdf/2006.10833.pdf or algorithms that do not utilize priors or facilitate prediction, such as PCMCI https://www.science.org/doi/pdf/10.1126/sciadv.aau4996 and its variants  https://openreview.net/forum?id=dYeUvLUxBQ can be beneficial to argue about the usability of the model. In accordance, comparison with other forecasting methods https://arxiv.org/pdf/1704.04110.pdf and simple average window benchmarks for the accuracy need to be added.\n\nA final vague point pertains to the identification of confounders. The model seems to allow for causal relationships for example s1 -> s2, s1 -> s3 to create spurious edges s2 -> s3 in case s2  precedes s3 since s2 and s3 change behaviour according to s1.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The current paper proposes a method to build a vector autoregressive model for multi-variate forecasting of time series data while recovering an underlying causal network between the time series. \nThe problem is cast as a linear program to predict the current value of a time series by a linear combination of all the series’ history with a predefined depth. \nThe least squares cost function is added to a reformulation of the adaptive lasso such that the undirected edges that act as coefficients in the cost function are penalized based on a prior undirected network. The edges of the prior network are sampled uniformly from a distribution in [0.2,1]. In addition, an l2 penalty is set to adjust the prior graph according to the observation, similar to an l2 regularisation in linear regression which is known to be equivalent to setting a Gaussian prior in the learning parameters of the MLE problem. The authors develop a similar proof for the same argument on the proposed model.\nThe method is proved to be nonconvex and an ADMM algorithm is derived and is proved to converge in stationary points. The results indicate an improved performance in forecasting accuracy and network reconstruction error in several real-world datasets.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "“The closer the I and j are in the graph the lower the penalty”: My confusion here pertains to the nature of A. If A is a weighted adjacency the elements show strength, not proximity. If A indicates distance i.e. in a number of hops, it is A + A^2 + A^3 etc.. That however would mean that A is no longer sparse and the computational burden for real-world graphs would be a problem. This could be clarified further.\n\nSome simulations indicating how the model behaves under non-stationary setting would make the paper more complete. Moreover, it is not clear if equation 16 indicates a normalization of solely the input or the whole time series. If it is the latter doesn’t it cause a change of the results?\n\nAlthough explained in 3.5 and used in 4.5.2, it is not perfectly clear to me, does the model for d > 1 breaks down to d different models i.e. graphs learnt simultaneously?\n\nCausality is inherently a directed measure so the intuition behind setting an undirected graph as a prior is not fully clear.\n\nTo verify the intuition of section 3.4, wouldn’t sampling the prior edges from a gaussian instead of uniform distribution improve the final reconstruction?\n\nTypos: ant, wi,j, the generalize prior", "broader_impact_concerns": "There is no obvious ethical concern for the paper's broader impact.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "A4YHIY6m1t", "reviewer": "Reviewer_GGay", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength:\n\nS1: The work has a solid theoretical foundation and support.\n\nS2: The authors conducted extensive experiments to evaluate the performance of the proposed method.\n\nWeakness:\n\nW1: Although the work has a solid theoretical foundation, the work lacks high-level intuition to illustrate the idea, which establishes a barrier to a broad audience.\n\nW2: Although the authors provided extensive preliminaries to help the reviewer delve into the research. However, some details of the network Granger causality should be clarified.\n\nW3: It is unclear how to construct the prior matrix Aprior. The concept is suddenly introduced in problem formulation without a detailed explanation. \n\nW4: The assumption in the problem formulation (i.e. if two nodes (i, j) are far away in the graph, then there is no Granger causality between the two-time series) seems to be problematic and needs further discussion.\n\nW5: Compared to existing works, it is still unclear the advantages of the proposed one using the graph prior knowledge.\n\nW6: The authors claimed the design can deal with high dimensional variables with limited samples. This case often means n<<d, where n is the number of samples and d is the dimension. However, the experiments fall short of evaluating the case.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper introduces a model to learn network Granger causality by using prior knowledge in the form of an undirected graph. The study aims to mitigate the challenge of estimating Granger causality graphs in high-dimensional scenarios with limited samples. Experimental results from synthetic and real-world datasets demonstrate the effectiveness of the proposed method over existing alternatives. The work has a solid theoretical and technical foundation. However, some serious issues need to be clarified and resolved. I lean toward a major revision for the work.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "The authors can refer to the weakness to make changes. Specifically, the authors should make the following changes.\n\nC1: The authors should explain why non-zero matrix entry Cij can capture the Granger cause between Xi and Xj?\n\nC2: Please add high-level intuition to illustrate the idea.\n\nC3: Please provide a detailed explanation of how to construct the graph prior Aprior in the problem formulation.\n\nC4: Is there a long-distance dependency between two nodes of Aprior that are far away? If yes, then the assumption may be problematic. If not, please justify the point.\n\nC5: The authors should conduct experiments on extremely scarce data and higher dimensional data to validate performance\n\nC6: \"Equation (2) capture specific temporal dependencies between the p time series\". However, p is the dimension of variables. Is this a typo?\n\nC7: Compared to existing works, please illustrate the advantages of the proposed one and the reasons.", "broader_impact_concerns": "No ethical concerns", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The reviewers provided insightful comments and the authors have responded to those comments. The clarity of presentation is much improved. One reviewer expressed some concerns about the extent of comparisons to other state-of-the-art methods. These concerns are reasonable. However, the authors scope their claim to rest upon the evidence: "Thus the incorporation of prior information has proven instrumental in achieving superior accuracy and robustness when compared to classical algorithms in this domain." Overall, the reviewers were in favor of acceptance, and I agree.
Decision: Accept (poster)