Efficient Network Automatic Relevance Determination
Authors: Hongwei Zhang
Ziqi Ye
Xinyuan Wang
Xin Guo
Zenglin Xu
Yuan Cheng
Zixin Hu
Yuan Qi
Venue: ICML 2025 poster (2025)
Keywords: Automatic Relevance Determination
Graphical Lasso
Primary area: general_machine_learning->scalable_algorithms

Abstract: We propose Network Automatic Relevance Determination (NARD), an extension of ARD for linearly probabilistic models, to simultaneously model sparse relationships between inputs $X \in \mathbb R^{d \times N}$ and outputs $Y \in \mathbb R^{m \times N}$, while capturing the correlation structure among the $Y$. NARD employs a matrix normal prior which contains a sparsity-inducing parameter to identify and discard irrelevant features, thereby promoting sparsity in the model. Algorithmically, it iteratively updates both the precision matrix and the relationship between $Y$ and the refined inputs. To mitigate the computational inefficiencies of the $\mathcal O(m^3 + d^3)$ cost per iteration, we introduce Sequential NARD, which evaluates features sequentially, and a Surrogate Function Method, leveraging an efficient approximation of the marginal likelihood and simplifying the calculation of determinant and inverse of an intermediate matrix. Combining the Sequential update with the Surrogate Function method further reduces computational costs. The computational complexity per iteration for these three methods is reduced to $\mathcal O(m^3+p^3)$, $\mathcal O(m^3 + d^2)$, $\mathcal O(m^3+p^2)$ respectively, where $p \ll d$ is the final number of features in the model. Our methods demonstrate significant improvements in computational efficiency with comparable performance on both synthetic and real-world datasets.

Efficient Network Automatic Relevance Determination
Hongwei Zhang 123 Ziqi 43
Xinyuan Wang 53 Xin Guo 3 Zenglin Xu Yuan Cheng 13 Zixin Hu 13
Yuan Qi 163
Abstract
1. INTRODUCTION
Multiple-input Multiple-output Regression is a powerful modeling framework widely applied in quantitative disci- plines such as finance and genomics (Zellner, 1962; Breiman & Friedman, 1997; Hastie et al., 2009; Wang, 2010; He et al., 2016). In biological research, this approach is of- ten used to explore how molecular-level features (micro- phenotypes) influence broader phenotypic traits (macro- phenotypes) (Stephens, 2013; Akbani et al., 2014). A typical example involves examining how gene expression levels or protein concentrations, influence disease states develop- mental outcomes. However, biological data often involve ultra-hight strr-higghdiimend ditthiiiennional features, with thousands of genes or proteins contributing to a limited number of observable sam- ples (Moon et al., 2019). Such high-dimensional settings present significant computational challenges.
We propose Network Automatic Relevance De- termination (NARD), an extension ARD for linearly probabilistic models, to simultancously model sparse relationships between inputs X RdxN and outputs Y R while capturing the correlation structure among the Y. NARD employs matrix normal prior which contains a pearsity-inducin parameter to identify and dis- card irrelevant features, thereby promoting spar- sity in the model. Algorithmically, it iteratively updates both the precision matrix and the rela- tionship between Y and the refined inputs. To mitigate the computational inefficiencies of the Om3 + dd) cost per iteration, we introduce Se- quential NARD, which evaluates features sequen- tially, and a Surrogate Function Method, lever- aging an efficient approximation of the marginal likelihood and simplifying the calculation of de- terminant and inverse of an intermediate matrix. Combining the Sequential update with the Surro- gate Function method further reduces computa- tional costs. The computational complexity per iteration for these three methods is reduced to Om + +p"). Om" + ) Om 3 + p²), respec- tively, where p < d is the final number of features in the model. Our methods demonstrate signif- icant improvements in computational efficiency with comparable performance on both synthetic and real-world datasets.
Empirical studies suggest that only a subset of specific molecular features significantly impacts the observed phe- notypes. This highlights the need for sparse regression models that focus on the most relevant features. reducing over-fitting and enhancing interpretability. Furthermore, macro-piie phenotypes often exhibit sparse network structures, where only a few phenotypes are connected via interaction relationship. This sparse structure highlights the need for approaches that capture not only the individual effects of in- puts but also the inter- detenddepensscces among multiple outputs. Emerging methodologies for multi-omics integration and cross-modal network information processing capitalize on this principle, reflecting the fact that biological mechanisms arise from the complex interplay of numerous molecular events and their interactions (Cohen al., 2022; Kristensen et al., 2014).
Artificial Intelligence Innovation and Incubation Institute, Fu- dan University, Shanghai, China "School of Data Science, Fudan University, Shanghai, China 'Shanghai Academy of Artificial In- telligence for Science, Shanghai, China "Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Cambridge, United Kingdom 'Eberly College of Science, Penn- sylvania State University, PA, United State "Zhongshan Hospital, Fudan University, Shanghai, China. Correspondence to: Zixin Hu <huzixin @ fudan.edu.cn.. Yuan Qi <qiyuan @ wudan.edu.cn...
Another related example is the study of expression quan- titative trait loci (eQTL). Some genetic variants can affect the expression of multiple genes, acting as potential con- founders in gene networks. Gene expression data alone are unable to fully capture the gene activities. Ignoring these effects can result in spurious associations, leading to false positives or false negatives. Incorporating covariates such as genetic variants from eQTL studies improves the estimation of relationships among genes at the transcriptional level.
Proceedings of the ndnd International Conference Machine Learning, Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).
In this paper, we focus on linearly probabilistic models

---

Efficient Network Automatic Relevance Determination
(Minka, 2000). Drawing inspiration from Automatic Rele- vance Determination (ARD) framework (MacKay, 1992; MacKay et al., 1994), we introduce an extension named Nes- work ARD (NARD) for multiple output regression. Specifi- cally, we place an ARD prior on the regression coefficient matrix, enabling us to determine which input features are relevant for predicting the outputs. Simultaneously, we apply an L1 penalty on the precision matrix to encourage sparsity, thereby modeling the dependencies among the out- put. Although ARD prior is effective for feature selection, it faces computational challenges in high-dimensional set- tings. Standard ARD methods optimized via Expectation- Maximization (EM) or type-ll II maximum likelihood incur an O(d)) computational cost due to matrix inversion, where is the number of features.
is the maximum likelihood estimator (MLE). The negative
log-likelihood function of Y, up constant, is given by
((W = Ar((Y-WXX)) T - log [S2]. (3)
It is noticed that 7(W, .2) is not jointly convex in W and Ω but is bi-convex, i.e., it is convex in W for fixed P and in ft for fixed W. A common approach to address this bi-convexity is to employ an alternating minimization strategy, also known as the block coordinate descent (BCD) method (Tseng, 2001). This technique iteratively updates the parameters by fixing one set while optimizing the other until convergence reached.
1.2. Related work
Joint mean-covariance estimation. The problem of joint multivariate variable and covariance selection has garnered significant attention in recent years. Rothman et al. (2010) proposed MRCE, which incorporates an 4 penalty to pro- mote sparsity in both W and Ω. This formulation leads to BCD algorithm for solving the problem. Formally, the penalized ongllikelihod is
To address this issue, design several novel algorithms within the NARD framework. Specifically, inspired by Tip- ping's greedy approach (Tipping & Faul, 2003), we develop the sequential update method, which sequentially adds and removes features. This approach allows the model to start with few features, enabling decisions about each new fea- ture based on its contribution to the overall evidence. Addi- tionally, we introduce a surrogate function that approximates the lower bound of the log marginal likelihood, avoiding the need for matrix inversion. By integrating these two ap- proaches, we provide a more efficient implementation of NARD, making scalable for high- -ieh-dimensional datasets.
(W.0) argmin (W,0)
((W, + A Σ lwnel + A2 Ayple
|uyl
(4)
where whe the elements of 5t wy are the elements of vectorized W. .,1, A2 > O are tuning parameters. This for- mulation yields an alternative lasso (Tibshirani, 1996) and a graphical lasso problem. Similar approaches have been explored in previous studies, such as (Cai et al., 2013; Chen et al., 2016; Lin et al., 2016; Zhang & Schneider, 2010; Zhao et al., 2020). Another approach is to adopt the Bayesian framework, where hierarchical model is constructed us- ing appropriate priors, followed by MCMC sampling for parameter updates and probabilistic uncertainty quantifica- tion, as seen in (Bhadra & Mallick, 2013; Deshpande et al., 2019; Li et al., 2021; Ha et al., 2021; Samanta et al., 2022). The Bayesian methods provide a principled way to handle model uncertainty while promoting sparsity. Compared to penalized methods, MCMC-based techniques tend to be computationally expensive, making them less efficient for large-scale problems.
In summary, we present the NARD framework, which jointly estimates the sparse regression coefficients and the precision matrix. To improve computational efficiency, we propose three novel extensions: Sequential NARD, Surro- gate NARD and Hybrid NARD, which reduce complexity to Om3 + p²), O(33 + d²) and O(m² + y²), respec- tively. These approaches achieve significant improvements in computational efficiency while maintaining comparable predictive performance in synthetic and real-world datasets.
1.1. Problem formulation
Under the framework of linearly probabilistic models, given an input x Idd and an output y E R", we consider
y = WI WI+
(L)
where W E Rxd is the regression coefficient matrix and (0, V) represents the error term, assumed to follow normal distribution with mean zero and covariance matrix V. Given N sample pairs, the model extends to:
ARD and SBL. Automatic Relevance Determination (ARD), closely related to Sparse Bayesian Learning (SBL) (Tipping, 2001; Faul & Tipping, 2001; Wipf & Rao, 2004), is framework designed to identify and discard irrelevant features from digediimeniional data. ARD leverages pa- rameterized, data-driven prior to promote sparsity, mitigat- ing the ill-posed nature of problems.
YmxN= = WmxdXaxN + EmxN,
(2)
where the i-th column of E es
Subsequent studies have theoretically extended ARD or SBL by establishing connections with iterative reweighted 4 minimization (Wipf & Nagarajan, 2007; 2010), compressive
Our goal is to jointly estimate the regression coefficient W and the precision matrix & = V-1. A typical approach
2

---

Efficient Network Automatic Relevance Determination
sensing (Babacan et al., 2009), stepwise regression (Ament & Gomes, 2021), by developing new efficient iterative algorithms (Al-Shoukairi et al., 2017; Zhou et al., 2021; Wang et al., 2024).
function (MLF) in Eq. (7). This is equivalent to minimize the negative logarithm of the MLF.
Recall ((W, 52) is not jointly convex in W and Ω we em- ploy the BCD method to iteratively update the parameters. This involves alternatively solving the sparse covariance es- timation problem and performing Bayesian linear regression with the ARD prior, cycling through the parameters until convergence is reached.
Graphical Lasso. The Graphical Lasso (GLasso) aims to estimate a sparse precision matrix for a given dataset (Friedman et al., 2008; Banerjee et al., 2008; Yuan & Lin, 2007). The general form of GLasso is defined as:
② = arg min (- log [S] + Tr(V2) + (x(()), (5) 12
2.2. Evidence approximation
Referring to the notation in (Minka, 2000), we define
where f2 is the estimated precision matrix, V is the empiri- cal covariance matrix and r (1) is the regularization term controlling the sparsity of the solution with strength param- eter A. When (S2) = " [wil| Eq. (5) represents the original GLasso (Friedman et al., 2008). Other repre- sentative choices for (52) include the adaptive LASSO or SCAD penalty (Fan et al., 2009) and graphical horseshoe (Li et al., 2019) in the Bayesian framework.
S XX T +K,
Syy = YYT,
(8)
S = YXT,
Syle Syn Syx az Sya-
By multiplying the likelihood function p(Y [W, x) with the conjugate prior Eq. (6) of W. i.e., p(WIV, K), we have:
2. NETWORK AUTOMATIC RELEVANCE DETERMINATION
In PY, W|X, K)
xTr V (WS,,, - 25 W + Syy]]
2.1. Framework of NARD
Tr V (W - Syz Sex (W- Syx) T
(9)
In this paper, impose ARD prior the regression efficient matrix and Li penalty on the precision matrix to encourage sparsity. Actually, there are other potential penal- ties that can be applied to the precision matrix, offering flexibility in model specification. These alternatives can be considered as plug-in options.
- 1/
As a result, the posterior distribution for the coefficient matrix W remains a matrix normal distribution and the MLF for Y can then be given by integrating out W from the joint distribution given in Eq. (9). Formally, we have:
Specifically, we impose the matrix normal distribution as the prior of W. The probability density function is given by
P( Y, V. K) = MN (m. V, PV = MV
(10)
W~ (N, Vmxm; Kild)
m///2
1 Tr TT(V WKWT)
(6)
where = Spx and Σ Saz are the posterior mean and column covariance of W, respectively. C I T KTK- "X is the columns covariance of Y.
(2/)md/2|V//4/2 x exp
2
where V and K-1 are two positive definite matrices rep- resenting the covariance matrices for rows and columns of W respectively. In this paper, we restrict K diag(ar, a.) be a diagonal matrix.
The negative logarithm of the MLF takes the form:
f = - In |Y|,,,
m In |C| + N |n TI(TT--YC--) (11)
The basic idea of ARD is to give the W.j independent pa- rameterized data-dependent priors. The hyperparameters V and K are trained from the data by maximizing the evidence P, V, K), which can be done via type-II maximum likelihood or Aopectation-Maxiiiannn
In practice we use Woodbury identity to calculate C-1:
--
(12)
(EM) algorithm.
Note that the term YC- in Eq. (11) can be re-expressed by introducing the latent variable W as follows:
P(Y,V, = / / Kdw.
(7)
YC-lyT (Y ,X) ~x + +uKØ
Here p(V) and pKK) are hyperpriors imposed on V and K, which will be specified later. To estimate V. K, W, an essential procedure is to maximize the marginal likelihood
= m
- WX) - -XX) + WKWT,
(13)
3

---

Efficient Network Automatic Relevance Determination
i.e. In V, K) In K,
Algorithm NARD
Input: Input data X, Y. Output: Estimated a, V. V.V W
Then graphical lasso to update V, i...
V, V.v1+ glasso V, A)
(14)
1: Initialize elements, V VMLL..
for & + 1 to T' do
where glasso refers to GLasso procedure, which takes the empirical covariance as imput and returns updated covari- ance V and precision matrix V-1. The parameter controls the sparsity the V-1 via the L penalty.
3: Compute Σ <- (K XXT)-1
4: Compute μr + YX
5: Update V according to (15).
V, V. + glasso (V,).
Update according to (16), (18) depending on whether the Gamma prior is included.
2.3. Algorithm summary
8: if max | then
We can incorporate hyperpriors for the parameters V and K in our model. The hicrarchical Bayesian approach allows for greater flexibility and robustness by capturing uncer- tainties in the prior distributions. We have p(W V) [PWW V)p(K)dd. If each a, follows common Gamma distribution, then the prior of W; is multidimen- sional Student distribution.
9: Break.
10: end if
11: Update W + μ..
12: end for
affect the main conclusions and derivations of the paper, may complicate the expressions for some variables. Due to space limitations, we focus solely on the prior for K. However, the technical conclusions of this paper remain applicable to the prior for V. A detailed discussion of the choice of hyperprior is provided in Appendix D.
When a specific prior is adopted, a and V reach their min- ima at the points where the gradient of the logarithm of the MLF is zero. It can be shown that
(Y -X) -x) + tKKu
Vnew
(15)
N
2.4. Extension to the nonlinear setting
Here we consider two different priors for a.
NARD can be naturally extended to address nonlinearity through kernel method. Specifically, we consider
Flat prior. If follows the flat prior, i.e., p(a) 1 then
new
m
(16)
ml ((Sa)) Da + (n V 11
= ( + + ((.) Polynomial, RBF,
where (XX) represents a nonlinear feature mapping that transforms the input space into a higher-dimensional space, enabling more flexible modeling of complex relationships.
Gamma prior. If a follows a Gamma prior Gamma(a, b):
p(a) II b
a
bo. -ba
(17)
r(a)
3. SEQUENTIAL NARD
where and b are the shape and rate parameters of the Gamma distribution, respectively. Then
3.1. Sequential update
Rather than pruning redundant or irrelevant features as in NARD, we employ greedy approach that sequentially adds and removes features. The key difference is that the original NARD requires (dd) computations at the beginning of training, whereas the sequential update method begins with almost empty model, consisting of only few features, which significantly reduces the initial computational burden. We adopt a fast sequential optimization method to efficiently update the hyperparameters inspired by (Faul & Tipping, 2001; Tipping & Faul, 2003; Ament Gomes, 2021).
+ 2a -
a (S )a + (μ T V 14)a +
(18)
Note that if any a, = 0, then W,, 0 and the corre- sponding feature effectively pruned from the model. In each iteration, we calculate = --2 = (K + XX
Y and update a, V as the above formula. We re- peat the process above until the number of iterations reaches ]
max (0) ot is significantly small, where l. is the iteration number. Algo- rithm 1 summarizes the proposed NARD.
the maximum max A 1 +)1
a
3.2. Fast optimization of evidence
Prior for W. conjugate prior for V is the inverse Wishart distribution. This will lead to a matrix-7 distribution (Gupta & Nagar, 2018) for the MLF of Y. Although this will not
To perform a sequential update on and K, we separate out the contribution of a single prior parameter a, from the MLF P(Y|X, V, K).

---

Efficient Network Automatic Relevance Determination
Thus, by rewriting C' C C + af "poo"" and using determinant and inverse lemma of matrix, the logarithm of MLF can be explicitly decomposed into two parts, one part denoted by ()),,, that does not depend on ag and another that does, i.e.,
Algorithm 2 Sequential NARD
Input: Input data X, Y, E.
Output: Estimated a, V, V-1
1: Initialization: model with few features, a elements,
V + VMLE. AL(a) < 00, f 4 O.
2: while AL(a) > and < T do
Tr(4q, 1-
((a) = m [Ina In(a; + 8:)] +
3: t
a3 +
(19)
4: select E (1,2..... d) randomly.
+ +(ax))
5: Calculate temporary Qi q. S, 8,, 7ho Oz.,
if ag 00. i..., pp is in the model then if 73s > then ma?
where pr ERN denotes the i-th row of X. 9: = YC
and 8; -P. Pi:
8:
+
else m
Theorem 3.1. Denote 1i. := Tr V-1) ) ms, then the global maximum of ((a) with respect to ax is
9:
10:
0g + oo, Delete pi from X.
11: end if
ms
> 0:
12: else
04
=
M
(20)
13:
if 1/: > 0 then
7/ S 0.
ma + m
14:
Add Φi to X.
15:
else
The detailed proof of Theorem 3.1 can be found in Appendix C.I. Furthermore, to simplify the maintenance and update of s; and qt we can exploit the following relations
16:
Skip pr continue.
17:
end if
18:
end if
19:
Calculate KA.K KA. CA. C7'.
a
8 = O,S,
VML.E Y'
q: =
(21)
20:
- S
04 S
21:
V.V- + glasso (VMLL,,)
where we define
22:
Calculate L, (La((a)) Inew-
Q Qc=YCTTPP; Pi:
S. pT CA Ap..
(22)
23:
if AT.(a) new > then
24:
Update X, q, 8, n, K,
Here A is the active subset of features, with
25:
else
26:
Undo all changes in this iteration.
(I+
(23)
27: Continue
28:
end if
XA and KAA are the sub-matrices of X and K corresponding to A. We denote S:A := X A + KA. Using Woodbury identity, we can write
29:
end while
Q: = Yyi YX_ (SAA) -1 XAPi
mm
(24)
means that the corresponding a elements are not co. Next, we randomly select i {1, 2,... .,, d), calculate d) and on. and update V using graphical lasso. Based on the values of т)i and a; and whether the i-th feature is currently included in the model, we decide whether to add it (o:: + me?
S: o. 9: v._ (( ( -APP
It can be shown that the computation in Eq. (24) involves only those features in the active set A that correspond to finite hyperparameters ag Let p < denote the final num- ber of features in the model, the computational complexity per iteration is O(p²).
). re- estimate it (a; + maž ), or delete it (a: < oo). Afterward, we compute the new MIF,, L(x) new- If increases, we retain the change; otherwise, we undo it. This process is repeated until L(a) converges.
ni
It is easy to verify that when a; + xo, both the prior and posterior of W exhibit a very high density around W, = Ō.
This indicates that when a; - co, we can remove the i-th feature of X. See details in Appendix C.2.
4. SURROGATE NARD
4.1. Overview of surrogate function method
3.3. Algorithm summary
An alternative approach to reducing computational complex- ity is to introduce a surrogate function. More specifically, we approximate p( V, K, W) with a surrogate function
Algorithm 2 summarizes the proposed Sequential NARD. We begin with a model containing only a few features, which
5

---

Efficient Network Automatic Relevance Determination
P V, W, W. W satisfying PYY W) max K, (25)
Thus, we derive the desired inequality by applying
Lemma 4.2, shown below:
W
Trig(W)] - Trig(W")| + Tr/Vo(W")) T (w - w
which establishes a tight lower bound on the original like- lihood function. In other words, our objective becomes to maximize the lower bound of the MLF. In fact, this is a common technique in optimization and variational inference (Hunter& & Lange, 2004; Kingma& & Welling, 2014; Sun et al., 2016). In particular, the most computationally expensive step of the original NARD is matrix inversion. To mitigate this, a lower bound is strategically chosen to approximate the matrices involved in the inversion as diagonal matrices. This approximation eliminates the bottleneck caused by ma- trix inversion, resulting in a substantial reduction in overall computational complexity.
+ 2
- [ - W))
-
(31)
Revisiting the MLF in Eq. (25), we obtain the approximate posterior density of W by substituting p(Y)X, V, K, W) with its lower bound p(Y| V.K, W, W') via Bayesian rules as follows:
X V. W. W (WWV, K) P(WIV, ~ J | V. W, W W")((WVV K)dw (32)
4.2. Lower bound of the MLF
MN V.Sz2).
Lemma 4.1. (Boyd & Vandenberghe, 2004) Suppose fX) is a function f Raxx. - R. the first-order Taylor approxi- mation with trace inner product is
with
1 = PW - W'XX - +YXT, S.x = pT + K.
(33)
[(X+V] R2 f(X) TX(S(() V).
(26)
Lemma 4.2. Let f: Raxa + be a continuously differ- entiable function with Lipschitz continuous gradient and Lipschitz constant L.. Then, for any U, V R
By substituting Eq. (32) into the marginal likelihood expres- sion, we arrive at the following formulation:
Proposition 4.4. The new marginal likelihood function, up to a constant, is
-V)] & |IU - V1I2. (f(U)=((V) - IVV -
2 (27)
L V,,,,,, W. W')
x m In + V + Tr(V- ((W,)")) (34) ++(V-WKWT).
Denote
R(W.W") = - -XX)(X) W
2(W - W W')x('XX T)
(28)
This reformulation leads to a new objective function involv- ing 4 sets of variables W, W", V, K. Direct joint optimiza- tion of these variables is intractable due to their interdepen- dence. Accordingly, we adopt the BCD method to optimize the negative logarithm of MLF, i.c.,
+ +(W - (W")W")")
where p E R denotes the largest eigenvalue of XXT, we have the following lemma:
Lemma 4.3. Let g(W)= = =Y - WX)(X - T then Tr()(W)) < Tr(,WW
(W(k), Wi(k), ((kk, k(()) E arg min L(W, W K.
(35)
Proof. We begin by analyzing the Lipschitz constant for the function g( gWW). First, we compute the gradient of gW) with respect to W:
Specifically, the BCD method is utilized to alternatively solve Eq. (35) as follows:
(k) E arg C(W, ((k-)), v(k-1), k(l-1)), W
@g(W) 2(Y WXXX
(29)
aw
Wh())) E arg min ((W()), W", V(k-1), x((k-1)),
According the definition of Lipschitz continuous gradient, the following inequality holds:
W ((k) e arg ((W(K), wr((k) V, k((h-1)).
(36)
(W) Og (W") )
(()) e arg ((W(K), ((()), v(a) K).
I L LI - W'll
aw aw
+ ||2(W - X'XXX || S L/I - W'll
(30)
The optimal value of ((A), wi(k), v(a) (k) can be ob- tained by setting its gradient to zero respectively, leading to
+ L = 2 2p.
6

---

Efficient Network Automatic Relevance Determination
the following update schemes:
Hybrid NARD retains the core structure and principles of Sequential NARD, with key variables such Q q, S, 8, and 1 preserved, and updates are performed only when the MLF increases. However, the key difference is that Hybrid NARD requires the computation of p and the execution of steps (37), (38), (39), (40), (41) and (42) at each iteration, rather than relying on the updates specified Sequential NARD, Consequently, W is updated in every iteration of Hybrid NARD, whereas in Sequential NARD, W is calculated only after the final values of X, K and V have been determined.
5th) - 1((-11) + pl:
(37)
W(K) = [oW - W + +XX (S
(38)
Wi(h) - w(k),
(39)
V(k) = (ww(),,wgg + w(t (()-))) (W())
N
(40)
V(A) (v(k))-1 + (V(),, 2),
(41)
k(bb)
m
6. EXPERIMENTS
diag
+ mm(5 )-1)
In this section, we evaluate the performance of the proposed methods on synthetic data using True Positive Rate (TPR), defined as TPR TF-FN- and False Positive Rate (FPR),
(42)
TP
4.3. Algorithm summary
defined as FPR = TN-FP- where TP, FP, TN and FN repre- sent true positives, false positives, true negatives, and false negatives, respectively. We run each experiment 10 times with different random seeds and report the average. For Aging phenotype data, where true labels are unavailable, we use the Jaccard index to measure the overlap of biologi- cal associations identified by different algorithms. We also highlight the efficiency of the NARD and its variants. On the TCGA data, we focus on the effectiveness of NARD in identifying biological associations. In addition, we con- ducted experiments on financial and air quality datasets to further demonstrate the versatility of our method; the de- tailed experimental results are provided in Appendix F.2. All experiments were performed on 32 Intel(R) Xeon(R) Platinum CPUs.
TP
Algorithm 3 Surrogate NARD
Input: Input data X, e Output: Estimated a, V, W
1: Initialize a elements, V < VALLE W()) t YX (X
XXT)-
2: for k k= 1 to K do
3: Perform equations (37). (38). (39), (40), (41), (42).
4: if [| - W(K- I then
5: Break.
6: end if
7: end for
Algorithm 3 presents a comprchensive overview of Surro- gate NARD. Since S.z a diagonal matrix, the update in (37) involves only the inversion of diagonal matrix, which has computational complexity of (d).
To select the optimal X, we employ a 5-fold cross- validaail procedure. The dataset is partitioned into 5 disjoint subsets, and in each iteration, subset held out as the validation set while the remaining subsets are used for model estimation. The objective function for selecting Aglasso is defined as:
5. Hybrid NARD
This section presents a hybrid method that integrates Se- quential NARD and Surrogate NARD to further reduce computational costs. Specifically, we utilize a sequential update approach for the iterative steps of Surrogate NARD. This approach initiates by assessing the relevance of new feature to determine its inclusion, akin to the Sequential NARD framework.
Aeaass arg min Σ
Tr 1og Σ (ggl)
(43)
Here Vi2 is the empirical covariance estimator computed from the training data excluding the (-th fold, and Ω_ is the estimated precision matrix based on this subset. The log-likelihood is computed for each fold, and the À that maximizes the cross-validated log-likelihool is chosen. A grid search is performed over a range of candidate values for X and the value that yields the best performance across all folds is selected for the final model and evaluation.
Once a feature is deemed relevant and included, we employ Surrogate NARD for subsequent calculations of the matrices V, W. This strategy eliminates the need for the costly matrix inversion operations typically required in full NARD imple- mentations. By combining the feature selection process of Sequential NARD with the efficient computations of Surro- gate NARD, we achieve a significantly lower computational complexity. The resulting method not only streamlines the iterative process but also ensures that we maintain compara- ble performance while handling high-dimensional data.
6.1. Synthetic datasets
The synthetic data are generated as follows. The covariance matrix V is constructed based on a graph generated by an
7

---

Efficient Network Automatic Relevance Determination
Erdös-Rényi random graph (Erdös& & Rényi, 1959) with sparsity parameter p. The entries of the precision matrix are sampled from uniform distribution, after which the matrix is symmetrized and adjusted to ensure positive definiteness by controlling its minimum eigenvalue. The matrix W is generated with a different sparsity level, and its non-zero elements are drawn from a uniform distribution. Both the data matrix X and the error term E are sampled from their corresponding multivariate normal distributions, and the outcome Y is computed as Y W + E.
Table 2. Impact of data size on performance. (p 0.02)
METHOO
d
m
N
TPR
FPR
TIME (1 STEP)
MRCE
1000
1000
5000
0,9802
0.0313
3.4
NARD
1000
1000
5000
0.9825
0.0274
3.5
SURROGATE NARD
1000
1000
5000
0.9818
0.0281
2.8
MRCE
2000
2000
10000
0,9797
0.0293
8.7
NARD
2000
2000
10000 D. 0,9864
0.0205
8.7
SURROGATE
2000
2000
10000 0.9838 0.0223
MRCE
5000 2000
20000 0,9698 0.0411
19.2
NARD
5000
2000 20000
0.9732 0.0323
17.4
SURROGATE NARD 5000
2000 20000
0,9732 0.0379
13.0
We choose two representative categories of baseline meth- ods: MRCE' (Rothman et al., 2010) and CAPME² (Cai et al., 2013) as frequency-bes based approaches, and HS-GHS (Li et al., 2021) and JRNS (Samanta et al., 2022) as Bayesian sampling-bes algorithms.
NARD
Time cmira
40
Surrogate NARD
2500
5000
75.00
10000
12500
15000
17500
20000
Table 1. Performance comparison of various methods. (p 0.1)
Dimensions ot
METHOD
m
N
TPR
FPR
TIME (TOTAL.)
Figure 1. Time comparison for different methods.
MRCE
5000
1500
1500
0.9083
0.0072
53
CAPME
1500
1500
0.8972
6.2. Baselines
HS-GHS
5000
1500 1500 0.9463 0.0033
>3000
JRNS
5000 1500 1500 0.9485 0.0037
>3000
We further evaluate the scalability of NARD on larger datasets. As shown in Table 3, the experimental results align with the theoretical complexity, demonstrating that as the data size increases, the time per iteration for NARD grows in a manner consistent with its expected computa- tional behavior. The Surrogate NARD approach consistently outperforms NARD, showing substantial improvements in computational efficiency, especially at larger scales.
NARD
5000 1500 1500 0.9483 0.0062
49
SEQUENTIAL NARD
5000 1500 1500 0.9459 0.0067
SURRODATE NARD
5000 1500 1500
0,9462 0.0072
HYBRID NARD
5000
1500
1500
0.9471 0.0068
23
Table presents the results with N 1500, m 1500, and d' 5000, highlighting estimation performance and CPU times for all methods. NARD and its variants perform similarly to HS-GHS, slightly outperforming MRCE in esti- mation, while demonstrating better computational efficiency. Among the methods, HS-GHS is the most time-consuming, while NARD- based approaches maintain relatively high effi- ciency, with performance gains becoming more pronounced as the dataset size increases. Table 2 compares performance under different combinations of c, m, and N, with focus on time performance relative to MRCE.
Table 3. Average time per iteration. (N 20000, m 2000)
d
MRCE
NARD
SURROGATE NARD
500
2.4
2.1
1.56
1000
2.5
2.4
1.96
2000
3.6
3.4
2.2
5000
12.2
10.7
3.7
10000
33.4
30.8
8.9
15000
77.5
70.9
20.8
Figure illustrates the running times of NARD and Surrogate NARD for varying values (d E (1500, 5000, 10000, 15000, 20000} with N = 1500 and
20000
201.9
168.6
33.7
30000
421.3
376.7
64.7
m = 1500). The results show that while the computational time of NARD increases significantly as d grows, Surrogate NARD exhibits lower time overhead. These trends are consistent with the theoretical complexity of the methods.
6.3. Aging phenotype data
We use data from a cohort of 1022 healthy individuals to construct a phenotypic network closely related to aging. In total, we identify 5641 phenotypes that exhibited a signif- icant Pearson correlation with age, including 1522 macro phenotypes and 4119 molecular phenotypes. Figure 2 illustrates the phenotypic network. Notably, fea- tures within the range of 280-370 exhibited a pronounced block structure, which corresponds to phenotypic data de-
MRCH/index. html
hhttp://www-stt stat weelat.whaaree www.:eeddeneeeeeeeeeeeeee s
htte::://githbb.com///ifff9//
'https://githinn
Stepwise
8

---

Efficient Network Automatic Relevance Determination
Table 4. Associations under different algorithms.
MITTIOD
MRCE HS-GHS NARD
NARD VARIANTS
SEQUENTIAL
SURROGATE HVBRID
15330
14963
15101
15014
15062 0.989
15039
R62
JAUCARD ANDER
0.979
0.988
Figure 3. UpSet plots illustrating the relationships among 10 path- ways across 7 cancer types.
Figure Phenotype network aging.
rived from the same tissue type. This structural coherence underscores the biological relevance of the identified pheno- types and suggests that similar traits may arise from shared underlying biological mechanisms. Furthermore, various al- gorithms demonstrate consistent performance on the dataset, with the Jaccard index for phenotype associations surpass- ing 98.5%. In terms of computational efficiency, the NARD method is the slowest, taking approximately 24 minutes to complete, while the sequential NARD approach take ap- proximately 14 minutes.
Figure 4. Protein network of COAD. Different colors represent different pathways.
Sparsity-inducing priors like ARD enhance interpretability in biological applications, such TCGA cancer data, by identifying key features. In our analysis across 7 tumor types, ARD highlighted important genes and proteins linked to signaling pathways. In Figure 3, sparsity revealed consis- tent pathways across cancer types, exposing cancer-specine translational effects. In Figure 4, for COAD, sparsity high- lighted critical protein interactions within pathways and cross-talk between them, aiding biological interpretation. In COAD, the PI3K/AKT pathway was highlighted by the interaction between GSK3ALPHABETAPSS2I and AK- TPS473 (Li et al., 2024). This association indicates a key regulatory role in tumor growth and survival. The AKT signaling axis, activated by various upstream kinases like GSK3, has been implicated in colon cancer progression, making it valuable target for further investigation and ther- apeutic development (Zhang al., 2021; Yao et al., 2020).
6.4. TCGA cancer data
To evaluate the effectiveness of NARD, we analyze data from The Cancer Genome Atlas (TCGA) (Weinstein al., 2013) across seven tumor types: colon adenocareinoma (COAD), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), ovarian serous cystadenocarcinoma (OV), rectum adenocarcinoma (READ), skin cutaneous melanoma (SKCM), and uterine corpus endometrial car- cinoma (UCEC). Each cancer type dataset includes mRNA expression profiles and RPPA-based proteomic data, reflect- ing the biological relationship where mRNA is translated into proteins. We choose 10 key signaling pathways based on recent studies (Akbani et al., 2014; Cherniack et al., 2017; Li et al., 2017) of RPPA-basee proteomic profiling across various tumor types.
7. Conclusion
In this paper, we introduce the NARD framework and pro- pose three variants alleviate its computational burden, significantly reducing the cost while maintaining perfor- mance. Experimental results confirm the effectiveness and efficiency of our methods across diverse domains. While the original model assumes linear relationships, we have demonstrated that NARD can be readily extended to nonlin- car settings via kernel-based techniques, further broadening its applicability.
Figure 3 visualizes the UpSet plots (Lex al., 2014) for 7 cancers in 10 pathways. The translational effects exhibit heterogeneity across cancer types, reflecting the expected biological differences. Figure 4 presents the protein associa- tion network for COAD, highlighting interactions between proteins within the same pathway, as well as significant cross-talk between different pathways. The network plots for the other cancer types are included in the Appendix F.4.

---

Efficient Network Automatic Relevance Determination
Acknowledgements
Cai, T. T., Li, H., Liu, W., and Xie, J. Covariate-adjusted pre- cision matrix estimation with an application in genetical genomics. Biometrika, (00(1):139-156 2013.
Authors acknowledge the constructive feedback of review- ers and the work of ICML'25 program and area chairs. This work was supported by the National Natural Science Foun- dation of China (Grant No. 82394432 and 92249302), and the Shanghai Municipal Science and Technology Major Project (Grant No. 2023SHZDZX02), Technical support was provided by the Human Phenome Data Center of Fudan University and the CFFF platform of Fudan University. This research was conducted during an internship at the Shanghai Academy of Artificial Intelligence for Science.
Chen, M., Ren, Z., Zhao, H., and Zhou, H. Asymptoti- cally normal and efficient estimation of covariate-adjuster gaussian graphical model. Journal of the American Sta- tistical Association, 1 111(513):394 406, 2016.
Cherniack, A. D., Shen, H., Walter, V., Stewart, C.,, Murray, B. A., Bowlby, R.,, Hu, X., Ling, S., Soslow, R. A., Broad- dus, R. R., et al. Integrated molecular characterization of uterine carcinosarcoma. Cancer Cell, 31(3):411-422 2017.
Impact Statement
Cohen, A. A., Ferrucci, L., Fülöp, T., Gravel, D., Hao, N., Kriete, A., Levine, M. E.,, Lipsitz, L. A., Olde Rikkert, M. G., Rutenberg, A., et al. A complex systems approach to aging biology. Nature Aging, 2(7):580-591, 2022.
This paper presents work whose goal is to advance the field of Machine learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.
Deshpande, S. K., Ročková, V., and George, E. I. Simulta- neous variable and covariance selection with the multi- variate spike-and-slab lasso. Journal of Computational and Graphical Statistics, 28(4): 2019.
References
Akbani, R., Ng, P. K. S., Werner, H. M., Shahmoradgoli, M., Zhang, F., Ju, Z., Liu, W., Yang, J.-Y., Yoshihara, K., Li, J., et al. A pan-cancer proteomic perspective the cancer genome atlas. Nature Communications, 5(1):3887, 2014.
Erdós, P. and Rényi A. On random graphs i. Publicationes Mathematicae Debrecen, 6(290-297):18, 1959.
Fan, J., Feng. Y., and Wu, Y. Network exploration via the adaptive lasso and scad penalties. The Annals of Applied Atatistics, 3(2):521, 2009.
Al-Shoukairi, M., Schniter, P., and Rao, B. D. A gamp- based low complexity sparse bayesian learning algorithm. IEEE Transactions on Signal Processing, 66(2):294-308, 2017.
Faul, A. and Tipping, M. Analysis of sparse bayesian leam- ing. Advances in Neural Information Processing Systems, 14, 2001.
Ament, E. and Gomes, C. Sparse bayesian learning via stepwise regression. In International Conference on Machine Learning, pp. 264-2744 PMLR, 2021.
Friedman, J., Hastie, T., and Tibshirani, R. Sparse inverse covariance estimation with the graphical lasso. Biostatis- tics, 9(3):432-444 441, 2008.
Babacan, S. D., Molina, R., and Katsaggelos, A. K. Bayesian compressive sensing using laplace priors. IEEE Transactions on Image Processing, ((1):53-633, 2009.
Gupta, A. K. and Nagar, D. K. Matrix variate distributions. Chapman and Hall/CRC, 2018.
Banerjee, O., El Ghaoui, L., and d' Aspremont, A. Model selection through sparse maximum likelihood estimation for multivariate gaussian or binary data. Journal of Ma- chine Learning Research, 9:485-516, 2008
Ha, M. J., Stingo, F. C., and Baludandayuthapani, V. Bayesian structure learning in multilayered genomic net- works. Journal of the American Statistical Association, 1 116(534):605-618 2021.
Bhadra, A. and Mallick, B. K. Joint high- -igeedimensional bayesian variable and covariance selection with an ap- plication equl analysis. Biometrics, 69(2):447-457 2013.
Hastie, T., Tibshirani, R., Friedman, J. H., and Friedman, J. H. The elements of statistical learning: data mining, inference, and prediction, volume Springer, 2009.
Boyd, S. and Vandenberghe, L. Convex optimization. Cam- bridge University Press, 2004.
He, D.,, Kuhn, D., and Parida, L. Novel applications of mul- titask learning and multiple output regression to multiple genetic trait prediction. Bioinformaties, 32(12):i37-i43, 2016.
Breiman, L. and Friedman, J. H. Predicting multivariate re- sponses in multiple linear regression. Journal of the Royal Statistical Society Series B: Statistical Methodology. 59 (1):3-54, 1997.
Hunter, D. R. and Lange, K. A tutorial on mm algorithms.
The American Statistician, 58(1)330337, 2004.
10

---

Efficient Network Automatic Relevance Determination
Kingma, D. P. and Welling, M. Auto-encoding variational bayes. In International Conference Learning Repre- sentations, 2014.
Computational and Graphical Statistics, 1944)9447-962
2010.
Samanta, S., Khare, K., and Michailidis, G. generalized likelihood-baser bayesian approach for scalable joint re- gression and covariance selection in high dimensions.
Kristensen, V. N., Lingjærde, O, C., Russnes, H. G., Vollan, H. K. M., Frigessi, A., and Børresen-Dale A.-L. Prin- ciples and methods of integrative genomic analyses in cancer. Nature Reviews Cancer, 14(5):299-311 2014.
Statistics and Computing, 32(3):47, 2022.
Stephens, M. A unified framework for association analysis with multiple related phenotypes. PloS One, 8(7):e65245, 2013.
Lex, A., Gehlenborg, N., Strobelt, H., Vuillemot, R., and Pfister, H. Upset: visualization of intersecting sets. IEEE Transactions on Visualization and Computer Gruphics, 20(12):1933-1992, 2014.
Sun, Y., Babu, P. and Palomar, D. P. Majorization- minimization algorithms in signal processing, commu- nications, and machine learning. IEEE Transactions on Signal Processing, 65(3):794-816 2016.
Li, J., Zhao, W. Akbani, R., Liu, W., Ju, Z., Ling, S., Vellano, C. P., Roebuck, P., Yu, Q.. Eterovic, A. K., et al. Characterization of human cancer cell lines by reverse- phase protein arrays. Cancer Cell, 31 2017.
Tibshirani, R. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society Series B: Statistical Methodology, 58(1):267-288, 1996.
Li, J., Liu, W., Mojumdar, K., Kim, H., Zhou, Z., Ju, Z.,, Kumar, S. V., Ng, P. K.-S., Chen, H., Davies, M. A., et al. A protein expression atlas on tissue samples and cell lines from cancer patients provides insights into tumor heterogeneity and dependencies. Nature Cancer, 5(10): 1579-1595, 2024.
Tipping, M. E. Sparse bayesian learning and the relevance vector machine. Journal of Machine Learning Research, 1(Jun):211 244, 2001.
Tipping, M. E. and Faul, A. C. Fast marginal likelihood maximisation for sparse bayesian models. In Interna- tional workshop on Artificial Intelligence and Statistics, pp. 276-283. PMLR, 2003.
Li, Y., Craig, B. A., and Bhadra, A. The graphical horseshoe estimator for inverse covariance matrices. Journal of Computational and Graphical Statistics, 28(3):747-777, 2019.
Tseng, P. Convergence a block coordinate descent method for nondifferentiable minimization. Journal of Optimiza- tion Theory and Applications, 109:475 4944, 2001.
Li, Y., Datta, J., Craig, B. A., and Bhadra, A. Joint mean- covariance estimation via the horseshoe. Journal of Mul- fivariate Analysis, 183:104716, 2021.
Wang, H. Sparse seemingly unrelated regression modelling: Applications in finance and econometrics. Computational Statistics & Data Analysis, 54(11):2866-288) -27777277777 2010.
Lin, J., Basu, S., Banerjee, M., and Michailidis, G. Pe- nalized maximum likelihood estimation of multi-layered gaussian graphical models. Journal of Machine Learning Research, 17(146):1-5 2016.
Wang, Y., Li, Yue, Z., and Yuan, Y. An iterative min-min optimization method for sparse Bayesian learning. In Proceedings of the 41st International Conference on Ma- chine Learning, volume 235 of Proceedings of Machine Learning Research, pp. 2024.
MacKay, D. J. Bayesian interpolation. Neural Computation, 4(3):415-447, 1992.
PMLR, 21-27 Jul
MacKay, J.
J. et al. Bayesian nonlinear modeling for the prediction competition. ASHRAE Transactions, 100(2): 1053-1062, 1994.
Weinstein, J. N., Collisson, E. A., Mills, G. B., Shaw, K. R., Ozenberger, B. A., Ellrott, K., Shmulevich, I., Sander, C., and Stuart, J. M. The cancer genome atlas pan-cancer analysis project. Nature Genetics, 45(10):1113-11 2013.
Minka, T. Bayesian linear regression. Technical report, Citeseer, 2000.
Moon, K. R., Van Dijk, D., Wang, Z. Gigante, S., Burkhardt, D. B., Chen, S., Yim, K., Flzen, A. d., Hirn, M. J., Coifmann R., Visualizing structure and transitions in digeeniimenion biological data. Nature Biotechnol- 2019.
Wipf. D. and Nagarajan, S. A new view of automatic rel- evance determination. Advances in Neural Information Processing Systems, 20, 2007.
Wipf, D. and Nagarajan, S. Iterative reweighted 4 and l2 methods for finding sparse solutions. IEEE Journal of Selected Topics in Signal Processing, 4(2):317-329, 2010.
ogy,
Rothman, A. J.. Levina, E., and Zhu, J. Sparse multivari- ate regression with covariance estimation. Journal of

---

Efficient Network Automatic Relevance Determination
Wipf, D. P. and Rao, B. D. Sparse bayesian learning for basis selection. IEEE Transactions Signal Processing, 52(8):2153-2164, 2004.
Yao, M., Li, R., Yang, Z., Ding, Y., Zhang, W., Li, W., Liu, M., Zhao, C., Wang, Y., Tang, H., al. Pp9, steroidal saponin, induces g2/m arrest and apoptosis in human colorectal cancer cells by inhibiting the pi3Wakk/gsk33 pathway. Chemico Cremico-Biological Interactions, 331:109246, 2020.
Yuan, M. and Lin, Model selection and estimation in the gaussian graphical model. Biometrika, 94(1):19-35, 2007.
Zellner, An efficient method of estimating seemingly un- related regressions and tests for aggregation bias. Journal of the American Statistical Association, 57(298):348-366 1962.
Zhang. N., Tian, Y.-N., Zhou, L.-N., Li, M.-Z., Chen, H.-D., Song, S.-S., Huan, X.J.,, Bao, X.B.. Zhang, A., Miao, Z.-H., et al. Glycogen synthase kinase 33 inhibition synergizes with parp inhibitors through the induction of homologous recombination deficiency in colorectal cancer. Cell death & disease, 12(2):183, 2021.
Zhang. Y. and Schneider, J. Learning multiple tasks with sparse matrix-normal penalty. Advances in Neural Infor- mation Processing Systems, 23, 2010.
Zhao, H., Stretcu, O.,. Smola, J., and Gordon, G. Effi- cient multitask feature and relationship learning. In Un- certainty in Artificial Intelligence, pp. 777-787, PMLR, 2020.
Zhou, W. Zhang, H.-T., and Wang, J. An efficient sparse bayesian learning algorithm based gaussian-scale mix- tures. IEEE Transactions on Neural Networks and Learn- ing Systems, 33(7):3065-300 -078066--778,,, 2021.
12

---

Efficient Network Automatic Relevance Determination
A. Matrix Computation
A.I. Trace and determinant
Tr(ABC) Tr(CAB)= Tr(BA)), IN + AB
(44)
|/M |/M+ABB)
VA,BERRMM
A.2. WoodBury identity
(A +
(45)
A.3. Derivatives of matrix-variate function
T BXC+B
0x
(T((AX)) = BT,
(46)
OX
0l X
(X-1)T.
ax
A.4. Contribution of
Let C = 1====XXK-- Gb a, pipp then we have
|C| |c + aj |p; ppl. pali
O pipt Ci7l
(47)
c- = G7.
ai + si C pi
B. Derivations in Detail
B.1. Details of Eq. (9)
pY, W|X, V. K) is the joint distribution of Y and W. By multiplying the likelihood function with the conjugate prior and omitting the constant, we get the kernel of pYY, W [X, K)).
Recall that the distribution of Y given under the model
N
p(]]X; = II
(ull,, W)
2 Σ (vi Wx:) T V-1 (vi Wx:)
1 V N/2
exp
(48)
- Tr (V-1(y- WXX(Y WX)
2x V N/2
exp
Tr (V-1 [WXXTWT. - 2 + vYT)))
N/2 exp
and the conjugate prior of W ~ MN Vmxm: KTla) is
K m2/
Kmm/2
Tr(V WKWW))
(W:V
(2T) md v
exp
(49)
13

---

Efficient Network Automatic Relevance Determination
Define
S = - K, YY Syx = YX Sylz Syy Syx S22 S
By multiplying the likelihood function pY|,W, with the conjugate prior of W, i.e., p VV K), and note that
W (W Syo) Sax (W- Szx WS WT - 25 WT + +Syz 22 ST
(50)
We have
X, WX
x Tr V- (WS W 2S yr WT +yuu
(51)
= =Tr [ (WS.zW WT WT- -2S - W
Syz Syn + Syiz)]
(52)
= =r V-1 (W- - Syr Sex (W- Syx + V-
(53)
B.2. Details of Eq. (10)
In Bayesian analysis with conjugate priors, the normalization constants are often omitted during intermediate steps, and only the kernel of the distribution is considered. Afterward, if the form of the kernel matches a known distribution, the normalization factor can be reintroduced. In our case, since the matrix normal distribution conjugate to the Gaussian likelihood, the posterior distribution of W also follows a matrix normal distribution. Therefore, we only need to match the kernel of the posterior distribution with the known form of the matrix normal distribution derive the posterior.
As result, the first term in Eq. (9) can be decomposed into two parts: the first part corresponds the kernel of the matrix normal distribution for W, and the second part corresponds to the kernel for Y.
Formally, we have p( W XX Y, V, K) ~ MMN (H, V, )), where column covariance of W, respectively. I+ K
Sye and E Sez are the posterior mean and is the column covariance of Y.
For Y, Sylz Syy Syn ST = =Y - - identity. Thus, we have X V, K) VN).
-C--YTT, where the third equality follows from WoodBury
2x y2
B.3. Details of Eq. (11)
~ (OOXV,K)~MM(....) (2n)mma///C(mm//////) exp
Tr(C Tr
(54)
Hence,
1 + (CC--yTv--y)
C= In P(Y|X, V, K) [m. 2n mlm Nn
(55)
x mmn1 + N V VV
(56)
B.4. Details of Eq. (13)
The term YC-lyT can be re-expressed by introducing the latent variable W as follows:
YC-lyT- = -Y(I S-IX
(57)
am
YY YX S-1 T
(58)
YYT -
(59)
= (Y - -xX)(Y - -XX) + +YXT TtT - ITT
(60)
(Y XX -X) + +WSral -XXX TT.
(61)
(Y -xX)(Y UX)T + 1( XX
(62)
= (V--XX) -VX)X - -xX) + TKKHT
(63)
min(Y - WX)(Y W + WKWT, W
(64)
14

---

Efficient Network Automatic Relevance Determination
where μ ST The last equality represents the variational form. The value of W that minimizes the expression is
W = 14, which be derived by taking the derivative of the objective function (Y (X WX) T + W WKWT with respect to W and setting it equal to zero.
B.5. Details of Eq. (15)
Under flat prior, it can be shown that
OC
NV +Y - =NNV+ xX)( - =xx =Kp T.
NV
(65)
0V-1
Hence,
Vnew (X) TX)T + PKU
(66)
N
B.6. Details of Eq. (23)
Eq. (21), Eq. (22) and Eq. (23) the standard quantity in sparse learning. To some extent, 8; is called the sparsity and qi is known as the quality of o, The sparsity measures the extent which basis function overlaps with the other basis vectors in the model, and the quality represents measure of the alignment of the basis vector with the error between the training set values and the vector of predictions that would result from the model with the vector excluded.
Since features are sequentially added to the model in Sequential NARD, we only need to consider the features already presented in the model at each step. So in this case, we define A as the corresponding active subset of features. Therefore, we add a subscript to the quantities in the analysis of this algorithm.
In Eq. (22), we have defined
Qi Q= Yc=Yc pa. 5
We have also defined qi = Y Pi and 8; = 2. c7. pi. However, these two quantities involves C which leads to the high computational cost.
In the proof of Eq. (21) below. We omit the subscript A presented in CA for readability.
CP'pi0 Cid
S = 1o_ [77
(67)
a4 pp C73 4p:
(0] C pi)2
= oo. G7. P2
(68)
0r + 8i
s7
5;
(69)
0G +
0;8i
(70)
04 + 84
By rearranging it, we can easily get 8,
a.S,
a.-3.
Qi = YC- pi
(71)
T.p cTl
=X
(72)
ai + qp C.t 8i
(Y pi) (p: To pi)
=YC 4i
(73)
a 8 8;
9:st
0:
(74)
+
15

---

Efficient Network Automatic Relevance Determination
By rearranging it, we get
(as + 8,)Qi
9i
(75)
0
0j + a.S.
Qi
a, S.
(76)
Oi
QiQi
(77)
03 S
C. Details of Sequential NARD
C.1. Proof of Theorem 3.1
Theorem C.1. Denote nι := Tr(q:g] VV1) )
ms,, then the global maximum of L(a) with respect ag
ma? 00
> 0
02
(78)
T S 0
Proof. The stationary points of marginal likelihood with respect to az occur when
@l(a)
1 m
m
Tr(a:q (V)) (a: 8;)2
0.
da:
2
05
Then we have a mai When Ti 0, of
0, then this is reasonable as value.
LL
m
m
2Tr(g:g, T V-1)
0o?
=
2
0
(a; + s,)2
(a: + 8;)3
(79)
2aa
s ms
at (a: 8;)3
Case 1. 7); 0:
When a; = aç the nominator is 2 m="s =
1/ 3 m
ms? m at
ms. So, it is negative.
Then, 1). ms?
mss; (already satisfied since here 7) 0). Also, denominator is positive. 84 (a < O. Therefore, at is the maximum point. Hence, da
{44(a))
Case 2. ?/: 5:: Then of & Dom(L) since we need a; for all i, so at can not be the optimal point here. Note that
aL
ms? + aays, a, Tr(g.T]
ms? a.77
day
(o: + s;)3
af + 8;)3
We can observe that the nominator and the denominator is always positive. So increasing with respect to a;, Hence, 0g + Do will make L as large as possible.
Recall the Theorem 3. 8; is called the sparsity and gi is known as the quality of pi, The sparsity measures the extent which basis function overlaps with the other basis vectors in the model, and the quality represents a measure of the alignment of the basis vector with the error between the training set values and the vector of predictions that would result from the model with the vector excluded. The term π) Tr(q:q V-1) ms; actually measures the trade-off between the alignment quality of the basis vector and its sparsity in relation to the covariance structure. For L(os), when > 0, the function exhibits an initial increase followed by decrease, with the maximum value occurring at a stationary point. When n S 0, the process is monotonically increasing, and the maximum value is asymptotically approached as a; - oc, consistent with the proof of Theorem 1.1. Furthermore, as a; + ao, the part of L(a;) that depends on a; diminishes, and T(az)) represents the situation where the corresponding feature can be pruned from the model.
16

---

Efficient Network Automatic Relevance Determination
C.2. Analysis of prior Observation. When a; co, W. - Ô. Tr(KWTV V- W)
Proof. and V are positive definite, so Kー1 and V-1 are also positive definite. For all R", T W - V-I 2 Hence, WVV- is positive semidefinite. Note, the trace of the product of two positive semidefinite matrices is nonnegative. So Tr (IKWVV-WW 2 0.
Prior
K m2/2
Tr (V WWW
K(
(2T) (2m)md//3 md/
(md))m//2/v/ 2//2//2
exp
When 0; and fix one [K] + 0o then A -
So,
0o const
W, 0 W, +0
Tr(KW W
Hence,
cxp
Tr(KW
(VK +
v-l
0.
if W. #0 if W2 +
2
const,
For matrix Raxb following a matrix normal distribution MM(MM,, Q), the density
U //2 m//2
1 Tr (A-
(2 (2m) exp
(80)
where M _aab is the expectation, IRaxa and Q R8bbb are two positive definite matrices representing the covariance matrices for rows and columns of A respectively.
D. Hyperprior
D.1. Hyperprior for
Flat prior. can be shown that
aC
-V + (Y -x)(Y TX)T + uKpT.
(81)
0V-1
Hence,
V0Ow (Y - PX)T + pK₁T
(82)
N
Inverse Wishart prior. If V follows inverse Wishart distribution:
p(V) W- (v.(vv))
| |
(++++))
- --
(83)
exp
The posterior for V is still inverse Wishart:
MV XVXXY ~ W + V.N N + +v),
(84)
and Y is matrix-7 T distribution.
Y Y~ -
(85)
17

---

Efficient Network Automatic Relevance Determination
Hence,
Vrew (Y 1V-1x)(0 (x)) + +uK + v
(86)
N + V
D.2. Hyperprior for a
Here we consider different prior for a.
Flat prior. follows the flat prior, i.e., p(o) 1, then
dc
m |as (Stz) - 1 + as (V)+))
(87)
Inf(r)
we have
m
0
neW
(88)
m STa )sc + (1 VTVVV p):
Gamma prior. If a follows a Gamma prior Gamma(a, b):
d
p(a)
a; a-1 1 ba.
eba.
(89)
P(a)
where a and are the shape and rate parameters of the Gamma distribution, respectively. Then
Oc m Ul((a:)
[as (S 1] + 2 z (n T Vo_V-=+) + - 1 1 -bas,
(90)
2
therefore,
2a +
(91)
m jú + (u V-1 2b
E. Relation with Other Work
noticed that the model in the paper "An "An Iterative Min-Min Optimization Method for Sparse Bayesian Learning" (Wang al., 2024) is designed for univariate regression, where the model is expressed as y Xw + 6, with y. R" X Raxd and IR.. In contrast, approach focuses multivariate regression model. Therefore, direct comparison between the two models is not entirely appropriate unless restrict the analysis to datasets with single outcome variable.
To extend the Min-Min SBL methods from that paper to the multivariate case, several modifications are necessary. For example, we need to adjust the prior distribution for the regression coefficient matrix. Additionally, to ensure that the negative logarithm of the marginal likelihood function can be decomposed as shown in Eq 12 Min-Min SBL, and to preserve the desirable properties of the Concave-Convex Procedure (CCCP) as discussed Lemmas 2.1 and 2.2, we may need introduce specific constraints into the model.
While extending the model may seem intuitively straightforward, the actual derivations not trivial and require additional effort. For instance, when we introduce the covariance matrix V, all iterative formulas must be re-derived within the CCP framework. ARD, in this sense, is more of a framework, where the solution process includes the classical Mackay update, the EM update, and methods like Min-Min SBL. our paper, we focus on improving the algorithmic complexity of the Mackay update step, while the extension offered by Min-Min SBL ARD are orthogonal scope of our current work.
F. Additional Experiments
F.1. Extension to the nonlinear setting
Our method can be naturally extended to address nonlinearity through kernel method. Specifically, we consider the model
y w (-) E Polynomial, RBF,
18

---

Efficient Network Automatic Relevance Determination
Table 5. Associations under different algorithms with kernels.
METHOD
MRCE
CAPME
HS-GHS
JRNS
NARD
NARD(POLYNOMIAL)
NARD(RBF)
# OF ASSOCIATION
15330
15094
14983
15066
15101 0.988
15094
15072
JACCARD INDEX
0.979
0.977
0.988
0.990
0.989
where (X) represents nonlinear feature mapping that transforms the input space into higher-dimennional space, allowing for more flexible modeling of complex relationships.
To explore this extension, consider 2 different kernel functions: the polynomial kernel and the Gaussian (RBF) kernel. We evaluate the performance on real-world aging phenotype data.
As shown in Table 5, our approach with polynomial and RBF kernels demonstrates competitive performance, achieving high Jaccard index values. The results are consistent with our expectation that kernel- based extensions allow the model to capture more complex, nonlinear relationships in the data, further validating the robustness and flexibility of our method.
F.2 Dataset diversity
We also expanded our experiments to include 2 non-biological datasets: Kuggle's ' air quality dataset and A-shares stock dataset.
For the air quality dataset, we performed data imputation and timestamp alignment, then analyzed the relationships among 11 key indicators. The results show strong correlation between PM2.5 and humidity, supporting the environmental principle that higher humidity promotes the adhesion of fine particles, leading to increased PM2.5 levels. This aligns with previous studies on atmospheric dynamics.
For the A-shares dataset, we collect nearly years of daily trading data and use the previous 5 days' information to predict the next day's opening price. This results dataset of 3032 stocks.
As shown in Table 6, experiments show that Bayesian methods, such HS-GHS and JRNS, were unable complete calculations in 4 days, while approach demonstrates excellent scalability. Analysis of the precision matrix reveals significant block structures, indicating that stocks from the same sector or industry tend show similar trends in price movement. Through these 2 experimental datasets, we have demonstrated the effectiveness of our method in both environmental and financial domains. Since there no ground truth, we used MRCE as baseline algorithm for comparison. We reported the Jaccard index as a benchmark. Additionally, we presented the computational time to highlight the computational advantages of our method.
Table 6. Associations of A-shares stocks. = 3032, 60640, N 1696)
METHOD
MRCE
CАРME
HS-GHS
JRNS
NARD
SEQUENTIAL NARD
SURROGATE NARD
NARD
ASSOCIATION
98335
99671
100309
99105
99475
JACCARD INDEX
0.869
0.881
0.891
0.893
0.890
TIME PER ITERATION (SECOND)
1200
1300
1000
255
TIME
16.5
14.5
5.5
3
E.3. Pathway of TCGA cancer data
Table 7 shows the sample size N, the number of predictors d, and the number of response variables m for the 7 datasets corresponding to each cancer type.
Table 8 presents a list of all pathways considered in the analysis of the TCGA cancer data this study, along with their respective protein members.
19

---

Efficient Network Automatic Relevance Determination
Table Datasets different cancer types.
CANCER
N
READ
121
73
76
LUAD
356
73
76
COAD
338
LUSC
309
73
86
OV
227
73
77
SKCM
333
73
76
UCEC
393
Table 8, Pathways and protein names.
PATHRAY SHORT PHORT
PROTEIN
CASPASEZCLEAVEDDIIS BADPS112, BCL2, BCLXL,
BREAST-REACTIVE
CAVEOLINI,
BETACATENIN, GAPDH. RBMIS
CHI CYCLE
CDXI CYCLINEI. CYCLINEL, CYCLINE2, P27PT157, P27PT198, PCNA, FOXM]
CORE-REACTIVE
CAVEOLINI BETACATENIN, RBM15, ECADHERIN, CLAUDIN7 CHKIPSS45, CHK2PT6S,
XRCCI
FIBRONECTIN, NCADHERIN, COLLAGENVI, CLAUDINT,
BETACATININ,
PЕΣKAKT
AKTPS473, AKTPT308, GSKIALLHABETAPSSISS GSK3PS9, P27PT157, P27PT198, PRAS40PT26, SUBERINPTI46S INPP4B, PTEN
RAS/MAPK
ARAFPS209, CIUNPS73, CRAFPSTS, INKPTISSYISS, MAPKPTHOY MEKIPS2/7SS 1PS2/75221, FISPTI80YI82 PORSKKTS599536, YRI PS102
SRCPY527
TSCMMOR
4EBPIPS65, 4EBPIPT37T46, 4EBPIPT70, P70S6KPT389, MTORPS248,, 56P82358236, 36PS2405244, RBPS807S8
E4. Protein network for different cancer
Figure 5 and Figure 6 illustrate the protein networks for UCEC and LUAD, respectively. These networks highlight the interactions between key proteins associated with each cancer type. Figure 7 and Figure show the networks for LUSC and OV, revealing distinct protein association patterns that may contribute to the unique characteristics of these cancers. Similarly, Figure and Figure 10 present the protein networks for READ and SKCM, providing further insight into the molecular landscape of these cancer types.
F.5. Further discussions
The Surrogate NARD sometimes demonstrate instability in its computations. This related to numerical issues that arise during the iterative optimization process. Specifically, when working with large, high-dimensional real-world datasets, the precision matrix estimation can become unstable due to the ill-conditioning of the covariance matrix or the challenges associated with ensuring that the precision matrix remains positive definite throughout the iterations. This numerical instability is not related to gradient-basee methods but rather stems from the nature of the data and the underlying optimization procedure. A potential solutions is that we can use more robust initialization for the precision matrix.
'nttps://archive ies. usi wwwp.eee/daattte/e//ea5/55/ http:///cccceeeoe.ou.o.eeeeeee//
/bei http://aatteeeeeeeeeeiiimiiiiiia tyydaa
20

---

Efficient Network Automatic Relevance Determination
Figure 5. Protein network of UCEC.
Figure 6. Protein network of LUAD.
Figure Protein network of LUSC.
21

---

Efficient Network Automatic Relevance Determination
Figure 8. Protein network of OV.
ETITE
Figure 9. Protein network READ.
Figure 10. Protein network of SKCM.
22


Reviews:
[{"review_id": "So6fg0xEVn", "reviewer": "Reviewer_AaWq", "summary": "This paper introduces Network Automatic Relevance Determination (NARD), an extension of Automatic Relevance Determination (ARD) designed for multiple-output regression in high-dimensional settings. NARD integrates a matrix normal prior with a sparsity-inducing mechanism to simultaneously select relevant input features and capture output dependencies.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful and constructive feedback on our submission.\n>Edge case of Theorem 3.1\n\nRecall the Theorem 3.1, $s_i$ is called the sparsity and $q_i$ is known as the quality of $\\varphi_i$, The sparsity measures the extent to which basis function overlaps with the other basis vectors in the model, and the quality represents a measure of the alignment of the basis vector with the error between the training set values and the vector of predictions that would result from the model with the vector excluded. The term $\\eta_i = \\text{Tr}(q_i q_i^{\\top}V^{-1}) - m s_i$ actually measures the trade-off between the alignment quality of the basis vector and its sparsity in relation to the covariance structure. For $L(\\alpha_i)$, when $\\eta_i > 0$, the function exhibits an initial increase followed by a decrease, with the maximum value occurring at a stationary point. When $\\eta_i \\le 0$, the process is monotonically increasing, and the maximum value is asymptotically approached as $\\alpha_i \\rightarrow \\infty$, consistent with the proof of Theorem 3.1. Furthermore, as $\\alpha_i \\to \\infty$, the part of $L(\\alpha_i)$ that depends on $\\alpha_i$ diminishes, and $L(\\alpha_i) = 0$ represents the situation where the corresponding feature can be pruned from the model. \n\n>Tuning the $\\lambda_{\\text{glasso}}$ Parameter in the ARD Framework\n\nTo select the optimal $\\lambda$,  we employ a 5-fold cross-validation procedure. The dataset is partitioned into 5 disjoint subsets, and in each iteration, 1 subset is held out as the validation set while the remaining 4 subsets are used for model estimation. The objective function for selecting $\\lambda_{\\text{glasso}}$ is defined as:\n$$\n\\lambda_{\\text{glasso}} = \\arg \\min_{\\lambda} \\sum_{l=1}^5 \n\\bigg[\n    \\textbf{Tr}(\\tilde{V_l}\\Omega_{-l}) -\\log |\\Omega_{-l}| + \\lambda \\sum_{\\substack{i \\neq j}} |\\omega_{ij}|\n\\bigg].\n$$\nHere $\\tilde{V_l}$ is the empirical covariance estimator computed from the training data excluding the $l$-th fold, and $\\Omega_{-l}$ is the estimated precision matrix based on this subset.\nThe log-likelihood is computed for each fold, and the $\\lambda$ that maximizes the cross-validated log-likelihood is chosen. A grid search is performed over a range of candidate values for $\\lambda$ , and the value that yields the best performance across all folds is selected for the final model and evaluation.\n\n>Discussion about Hybrid NARD\n\nIn our experiments, we did not observe cases where the Hybrid NARD underperformed compared to either Sequential or Surrogate NARD individually. Despite incorporating approximations, the hybrid approach consistently provided a good balance between accuracy and efficiency. In the synthetic data shown in Table 1, we found that NARD and its variants performed similarly, with no significant differences. However, in terms of time efficiency, the Hybrid NARD was approximately twice as fast as the standard NARD. This aligns with our theoretical expectations, where combining sequential optimization with surrogate modeling helps leverage their respective advantages.\n\nWe acknowledge that performance trade-offs may become more pronounced in extreme cases. For smaller datasets, we hypothesize that the hybrid approach may not always outperform NARD. For instance, when tested with $d=80, m=50,N=100$ and $d=50, m=20,N=50$, the results (shown in Table R1) demonstrated that Hybrid NARD still performed well, highlighting the robustness of our method. Despite our efforts to identify challenging edge cases, the algorithm exhibited a notable degree of stability.  \n\nWe appreciate your suggestion and will explore these scenarios further in future work.  \n\nTable R1: Performance Comparison of Our Methods on Small datasets.\n| Method | d | m | N | TPR | FPR |\n| --- | --- | --- | --- | --- | --- |\n | NARD | 80| 50| 100 | 0.9695| 0.0031 |\n | Sequential NARD | 80| 50| 100| 0.9689| 0.0033 |\n|Surrogate NARD |80 |50 |100| 0.9693| 0.0029 | \n |Hybrid NARD | 80| 50 | 100| 0.9689| 0.0035 |\n | NARD | 50| 20|50 | 0.9542| 0.0051 |\n | Sequential NARD | 50| 20| 50| 0.9544| 0.0046 |\n|Surrogate NARD |50 |20 |50| 0.9540| 0.0050 |\n |Hybrid NARD | 50| 20 | 50| 0.9540| 0.0049 |", "claims_and_evidence": "Yes.", "methods_and_evaluation_criteria": "Yes.", "theoretical_claims": "I have checked the theoretical derivations in Sections 3 and 4, and I find them reasonable. For example, in Section 4, I verified Lemma 4.2, which establishes an upper bound on $ \\text{Tr}[g(W)] $ using a majorization argument. The proof correctly applies Lipschitz continuity, with $ L = 2\\|XX^\\top\\| = 2\\rho $, and a first-order Taylor approximation to derive the bound. The use of $ \\rho $ (the largest eigenvalue of $ XX^\\top $) ensures a valid approximation, confirming the correctness of the inequality.", "experimental_designs_or_analyses": "I have checked the soundness and validity of the experimental designs and analyses, and I find them reasonable. For example, in the synthetic data experiments, the covariance and precision matrices were generated using an Erdős-Rényi random graph, ensuring a structured yet realistic sparsity pattern. The metrics used (TPR and FPR) are appropriate for evaluating feature selection.", "relation_to_broader_scientific_literature": "The paper extends Automatic Relevance Determination (ARD) by incorporating a matrix normal prior to model both feature sparsity and output dependencies, addressing limitations in traditional ARD. NARD in this paper further improves computational efficiency by introducing Sequential NARD, which employs a greedy approach that sequentially adds and removes features, and Surrogate NARD, which introduces a surrogate function to approximate the marginal likelihood.", "essential_references_not_discussed": "I think all related works have already been cited.", "other_strengths_and_weaknesses": "The paper presents a novel extension of Automatic Relevance Determination (ARD) by incorporating a matrix normal prior. The introduction of Sequential NARD and Surrogate NARD significantly improves computational efficiency over traditional ARD methods, which is a notable strength.\n\nThe paper could provide more clarity on hyperparameter selection. Additionally, the paper could further discuss potential trade-offs in Hybrid NARD.", "other_comments_or_suggestions": "Theorem 3.1 could benefit from additional explanation for edge cases and it would be helpful to provide more details on how hyperparameters were tuned.", "questions_for_authors": "The Hybrid NARD approach combines Sequential and Surrogate NARD. Are there any cases where this hybrid method performs worse than using either method individually? How does it balance efficiency vs. accuracy in practice? Do you have more experiments for this point?", "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": "", "extra_scores": {}}, {"review_id": "jIKnhTElxN", "reviewer": "Reviewer_gDpz", "summary": "This paper introduces the Network Automatic Relevance Determination (NARD) framework for linearly probabilistic models. It proposes three novel algorithms, i.e. Sequential NARD, Surrogate NARD, and Hybrid NARD, which significantly reduce the computational complexity. These methods maintain comparable performance on synthetic and real-world datasets, effectively handling the sparse relationships between inputs and outputs while capturing output correlations.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful feedback. We appreciate your positive remarks and the concerns raised, which will help us refine both the theoretical and empirical aspects of our work.\n\n>Baseline comparison\n\nAs suggested, we add additional baseline methods: CAPME[1] and JRNS[2]. MRCE and CAPME are frequency-based representative methods, while HS-GHS and JRNS are Bayesian sampling-based algorithms. Our experiments show that NARD still outperforms these methods, which further strengthens the case for the proposed approach. The results are shown in Table R2.\n\nIn the original paper, we carefully selected MRCE and HS-GHS as the comparison methods because they represent two well-established approaches in the field: MRCE is a frequentist method, and HS-GHS is a Bayesian approach. Both of these methods have been compared with several other techniques in their respective papers, and in those comparisons, MRCE and HS-GHS consistently performed well. We believe these methods serve as strong baselines for our study, providing a comprehensive comparison between frequentist and Bayesian paradigms in sparse multivariate regression and graphical model estimation. \n\nTable R2: Performance Comparison of Various Methods.\n| Method | d | m | N | TPR | FPR | Time |\n| --- | --- | --- | --- | --- | --- |--- |\n |MRCE | 5000| 1500|1500 | 0.9083| 0.0072 | 53 | \n|**CAPME** | 5000| 1500|1500 | 0.8972| 0.0124 | 52 | \n|HS-GHS | 5000| 1500|1500| 0.9463| 0.0033 | >3000 |\n|**JRNS** | 5000| 1500|1500| 0.9485| 0.0037 | >3000 |\n|NARD | 5000| 1500|1500 | 0.9483| 0.0062 | 49|\n|Sequential NARD | 5000| 1500| 1500| 0.9459| 0.0067 | 35|\n|Surrogate NARD |5000 |1500 |1500  | 0.9462| 0.0072 | 31|\n |Hybrid NARD | 5000| 1500 | 1500| 0.9471| 0.0068 | 23| \n\nWe also include the results of experiments on aging phenotype data as shown in Table R3.\n\nTable R3:  Associations under different algorithms.\n| Method | MRCE |CAPME |HS-GHS |JRNS | NARD |NARD(Polynomial) | NARD(RBF) |\n| --- | --- | --- | --- | --- | --- |--- | --- |\n| # of association | 15330 |15094| 14983|15066| 15101|15094| 15072|\n| Jaccard index | 0.979 |0.977|-| 0.988 |0.988 | 0.990 |0.989 |\n\n[1] Covariate-adjusted precision matrix estimation with an application in genetical genomics, Biometrika 2013.\n\n[2] A generalized likelihood based Bayesian approach for scalable joint regression and covariance selection in high dimensions, Statistics and computing 2022. \n\n>Non-linearity\n\nIn our paper, Section 6.1 focuses on synthetic data, which follows a linear structure by design. However, in Sections 6.2 and 6.3, we evaluate our method on real-world biological datasets with more complex, nonlinear relationships. Despite these nonlinearities, our approach performs well, suggesting that our sparse linear approximation effectively captures dominant interaction patterns. We will clarify this fact more clearly in the final version.\n\nAdditionally, our method can be naturally extended to address nonlinearity through kernel method. Specifically, we consider the model  \n$$\nY = W\\Phi(X) + \\mathcal{E}, \\quad \\Phi(\\cdot) \\in { \\text{Polynomial, RBF, ...} }\n$$\nwhere $\\Phi(X)$ represents a nonlinear feature mapping that transforms the input space into a higher-dimensional space, allowing for more flexible modeling of complex relationships.\n\nTo explore this extension, we consider 2 different kernel functions: the polynomial kernel and the Gaussian (RBF) kernel. We evaluate the performance on a real-world aging phenotype data.\n\nAs shown in Table R3, our approach with the polynomial and RBF kernels demonstrates competitive performance, achieving high Jaccard index values. The results are consistent with our expectation that kernel-based extensions allow the model to capture more complex, nonlinear relationships in the data, further validating the robustness and flexibility of our method. \n\n>Numerical stability in Surrogate NARD\n\nWe appreciate your attention to this issue, which we discuss briefly in Appendix F.4. It arises from numerical challenges encountered during the iterative optimization process, particularly when estimating the precision matrix. For large datasets, the covariance matrix can be ill-conditioned, leading to instability in the precision matrix estimation. This stems from the inherent properties of the data. As we mentioned in the discussion, one potential solution is to apply a more robust initialization for the precision matrix, which may help mitigate this issue. We will explore this aspect more deeply in future work.", "claims_and_evidence": "I didn't find any obvious problem.", "methods_and_evaluation_criteria": "Generally speaking, the proposed methods and evaluation criteria in the paper make sense for the problem at hand. Some of my concerns might stem from whether the proposed methods are SOTA. Given that the benchmark methods for comparison are relatively limited, as far as I know, there should be numerous methods for handling sparse multivariate regression and learning sparse graph structures. It would be better if the authors could explain why they chose MRCE and HS - GHS as the comparison techniques.", "theoretical_claims": "Sorry, I didn't check the correctness of proofs.", "experimental_designs_or_analyses": "1. The paper only compares NARD and its variants with two baseline methods (MRCE and HS-GHS). Given the large number of methods available for sparse multivariate regression and graphical model estimation, this limited comparison may not fully establish the superiority of the proposed methods. \n2. The authors assume linear relationships in their models. While this is a common starting point, real-world data, especially in biology, may contain complex nonlinear relationships. The experimental designs do not explore how well the methods perform in the presence of nonlinearity. \n3. The Surrogate NARD method shows instability in computations, especially when dealing with large, high-dimensional real-world datasets.", "relation_to_broader_scientific_literature": "In the related work section, the authors establish connections between the content of this paper and the broader scientific literature. Additionally, in the experimental part, authors compares NARD and its variants with two baseline methods (MRCE and HS-GHS).", "essential_references_not_discussed": "N/A", "other_strengths_and_weaknesses": "N/A", "other_comments_or_suggestions": "N/A", "questions_for_authors": "N/A", "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": "", "extra_scores": {}}, {"review_id": "qW4Ka4hDWj", "reviewer": "Reviewer_YET1", "summary": "The paper introduces **Network Automatic Relevance Determination (NARD)**, an extension of Automatic Relevance Determination (ARD) designed for linearly probabilistic models. NARD aims to simultaneously model sparse relationships between input features X and output responses Y, while capturing correlations among the outputs Y. The method employs a **matrix normal prior** with a sparsity-inducing parameter to identify and discard irrelevant features, thereby promoting sparsity in the model.\n\nThe algorithm iteratively updates both the precision matrix and the relationship between Y and refined inputs. To address computational inefficiencies associated with the high cost per iteration, the authors propose two enhancements:\n\n1. **Sequential NARD**: Evaluates features sequentially to reduce computational overhead.\n  \n2. **Surrogate Function Method**: Uses an efficient approximation of the marginal likelihood, simplifying determinant and matrix inverse calculations.\n  \n\nBy combining these approaches, the computational complexity is further reduced. The paper demonstrates that these methods achieve significant improvements in computational efficiency while maintaining comparable predictive performance on both synthetic and real-world datasets.\n\n\n## Update after rebuttal\nThank you for your response. I will maintain my positive score.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your valuable feedback and constructive comments. We appreciate your positive remarks on the novelty of our methods and the clarity of the paper.\n\n>Dataset diversity\n\nWe have expanded our experiments to include 2 non-biological datasets: Kaggle’s air quality dataset (https://archive.ics.uci.edu/dataset/501/beijing+multi+site+air+quality+data) and A-shares stock dataset.\n\nFor the air quality dataset, we performed data imputation and timestamp alignment, then analyzed the relationships among 11 key indicators. The results show a strong correlation between PM2.5 and humidity, supporting the environmental principle that higher humidity promotes the adhesion of fine particles, leading to increased PM2.5 levels. This aligns with prior studies on atmospheric dynamics.\n\nFor the A-shares dataset, we collect nearly 7 years of daily trading data and use the previous 5 days' information to predict the next day's opening price. This results in a dataset of 3032 stocks. Our experiments show that Bayesian methods, such as HS-GHS and JRNS, could not complete calculations within 4 days, while our approach demonstrate excellent scalability. Analysis of the precision matrix reveals significant block structures, indicating that stocks from the same sector or industry tend to show similar trends in price movement.\n\nThrough these 2 experimental datasets, we have demonstrated the effectiveness of our method in both environmental and financial domains.\n\nSince there is no ground truth, we used MRCE as a baseline algorithm for comparison. We reported the Jaccard index as a benchmark. Additionally, we presented the computational time to highlight the computational advantages of our method.\n\nTable R4:  Associations of A-shares stocks. ($m=3032, d=60640, N=1696$)\n| Method | MRCE |CAPME |HS-GHS |JRNS | NARD |Sequential NARD | Suggorate NARD | Hybrid NARD |\n| --- | --- | --- | --- | --- | --- |--- | --- |--- |\n| # of association |97939|98335| -|-|99671|100309|99105|99475|\n| Jaccard index |-|0.869|-|-|0.881|0.891|0.893|0.890|\n| Time per iteration (second)|\\~1200|\\~1300| -|-|~1000|-|255|-|-|\n| Time all (h)|\\~17|\\~16.5|-|-|\\~14.5|\\~8|~5.5|\\~3|\n\n>Interpretability in biological applications\n\nSparsity-inducing priors like ARD enhance interpretability in biological applications, such as TCGA cancer data, by identifying key features. In our analysis across 7 tumor types, ARD highlighted important genes and proteins linked to signaling pathways. In Figure 3, sparsity revealed consistent pathways across cancer types, exposing cancer-specific translational effects. In Figure 4, for COAD, sparsity highlighted critical protein interactions within pathways and cross-talk between them, aiding biological interpretation. In COAD, the PI3K/AKT pathway was highlighted by the interaction between GSK3ALPHABETAPS21S9 and AKTPS473. This association indicates a key regulatory role in tumor growth and survival. The AKT signaling axis, activated by various upstream kinases like GSK3, has been implicated in colon cancer progression, making it a valuable target for further investigation and therapeutic development.\n\n>Question about Table 2\n\nWe understand your concern about the exclusion of Hybrid NARD and Sequential NARD in Table 2. Table 2 shows single update step times, while Table 1 presents total computation time. Since both methods involve iterative updates with varying step times, direct comparison in Table 2 is difficult. We apologize for any confusion and will clarify this in the final version.\n\n>Linear assumptions\n\nWe have considered extending the framework to non-linear models, which can be easily adapted to sparse kernel regression without significant changes to NARD. Relevant experiments and results are discussed in our rebuttal to Reviewer gDpz.\n\n>Numerical stability\n\nWe appreciate your attention to this issue, which we have discussed in detail in Appendix F.4. We have provided additional analysis in our rebuttal to Reviewer gDpz, which you can refer to for a more in-depth discussion.\n\n>Discussion about Hybrid NARD\n\nWe have included further analysis on this issue in our rebuttal to Reviewer AaWq, which provides a more detailed discussion.", "claims_and_evidence": "The claims made in the paper are generally well-supported by clear evidence, as outlined below:\n\n1. **Sparse Feature Selection via ARD Prior**  \n  The paper claims that placing an ARD prior on the regression coefficient matrix enables effective feature selection by identifying relevant input features for predicting outputs.\n  \n2. **Sparsity in Output Dependencies via L1 Penalty**  \n  The use of an L1 penalty on the precision matrix to model dependencies among outputs is a reasonable and theoretically sound approach. Sparse precision matrices are widely used in multi-output regression to capture conditional independence relationships between outputs.\n  \n3. **Computational Challenges in High Dimensions**  \n  The claim that standard ARD methods incur $O(d^3)$ computational costs due to matrix inversion is accurate and consistent with the literature.\n  \n4. **Efficiency of Proposed Algorithms**  \n  The paper introduces Sequential NARD, Surrogate NARD, and Hybrid NARD to reduce computational complexity. The stated reductions to $O(m^3+p^3)$, $O(m^3+d^2)$, and $O(m^3+p^2)$, respectively, are plausible given the described modifications (e.g., sequential updates and surrogate function approximations).", "methods_and_evaluation_criteria": "The evaluation criteria and datasets are generally appropriate for the problem. TPR and FPR on synthetic data effectively measure feature selection performance, while the Jaccard index is suitable for comparing biological associations in the absence of ground truth. The use of TCGA data validates NARD’s real-world applicability. However, the empirical evaluation has notable gaps:\n\n1. **Limited Dataset Diversity**: The evaluation focuses heavily on biological datasets (e.g., TCGA and aging phenotype data). Including experiments on non-biological datasets (finance?) would better demonstrate the generalizability of NARD across diverse domains.\n  \n2. **Limited Baselines**: Comparisons are limited to MRCE and HS-GHS. Incorporating additional state-of-the-art baselines would strengthen validation and contextualize NARD's performance more comprehensively.\n  \n\nAddressing these issues would significantly enhance the robustness and completeness of the evaluation.", "theoretical_claims": "The theoretical claims in the paper are well-supported by the proofs provided in the text. Key claims include:\n\n1. Matrix Normal Prior and ARD Framework\n  \n2. Complexity Reductions\n  \n3. Theorem 3.1 (Sequential NARD)\n  \n4. Surrogate Function Approximation", "experimental_designs_or_analyses": "The experimental designs and analyses in the paper are sound and appropriate for the stated objectives:\n\n1. **Synthetic Data**: The use of synthetic datasets with controlled sparsity is well-suited for validating the feature selection and dependency modeling capabilities of NARD and its variants.\n  \n2. **Real-World Applications**: The use of Aging Phenotype data (evaluated via Jaccard index) and TCGA cancer data (focusing on biological associations) aligns with the goals of demonstrating NARD's applicability to high-dimensional, multi-output regression problems.\n  \n3. **Baseline Comparisons**: Comparisons with established methods like MRCE and HS-GHS are appropriate for benchmarking both computational efficiency and estimation performance. The inclusion of time as a metric strengthens the analysis of computational complexity claims.\n\nIssues are discussed below.", "relation_to_broader_scientific_literature": "The paper builds on foundational concepts in Automatic Relevance Determination (ARD), introduced by MacKay (1992), extending it to multi-output regression with a matrix normal prior and sparsity-inducing penalties. It addresses computational challenges in high-dimensional settings, improving efficiency through Sequential and Surrogate methods inspired by Tipping's work (2003). The framework aligns with broader literature on Bayesian sparse modeling, such as MRCE (Rothman et al., 2010) and graphical lasso techniques (Friedman et al., 2008), while contributing novel algorithmic advancements for scalable multi-output regression.", "essential_references_not_discussed": "N/A", "other_strengths_and_weaknesses": "**Strengths:**\n\n1. **Originality**: The paper introduces novel extensions to ARD, including Sequential NARD, Surrogate NARD, and Hybrid NARD, addressing computational inefficiencies in high-dimensional regression tasks.\n2. **Significance**: The framework is highly relevant for sparse modeling in biological and genomic datasets, making it impactful for real-world applications.\n3. **Clarity**: The paper clearly explains the theoretical foundations, algorithms, and experimental results, ensuring accessibility for readers familiar with Bayesian modeling.\n\n**Weaknesses:**\n\n1. **Numerical Stability**: Surrogate NARD exhibits instability during precision matrix estimation in high-dimensional datasets, which could limit its reliability in practice.\n2. **Linear Assumptions**: The reliance on linear models may restrict applicability to scenarios involving complex non-linear relationships.\n3. **Limited Dataset Diversity**: While synthetic and biological datasets are used, additional experiments on diverse real-world datasets could enhance generalizability.\n4. Comparisons with baseline methods (MRCE and HS-GHS) are limited in scope. Including more state-of-the-art baselines would provide a clearer picture of NARD's relative performance.", "other_comments_or_suggestions": "- Include experiments on non-biological datasets to evaluate generalizability across domains.\n  \n- Include more detailed comparisons with alternative methods beyond MRCE and HS-GHS to strengthen the empirical validation.\n  \n- Address the scalability of Hybrid NARD explicitly in larger datasets.\n  \n- Provide a more detailed analysis of how sparsity-inducing priors affect interpretability in biological applications.\n  \n- Suggestion for Images:  It would be better to move some images, such as protein networks for COAD, to the supplementary materials. This would declutter the main text and allow the remaining figures to be enlarged for improved readability.", "questions_for_authors": "- How does Hybrid NARD balance computational efficiency and predictive accuracy compared to Sequential and Surrogate NARD? Could you provide additional insights into its practical advantages?\n- Have you considered extending the framework to non-linear models or incorporating kernel methods for capturing complex relationships? If not, what challenges do you foresee?\n- Could you elaborate on how the sparsity-inducing priors affect interpretability in biological applications like TCGA cancer data?\n- To improve the empirical evaluation, could you expand the experiments to include non-biological datasets to assess generalizability across diverse domains, and explicitly address the scalability of Hybrid NARD on larger datasets, including runtime and performance metrics?\n- Why are Hybrid NARD and Sequential NARD not included in Table 2, where the impact of data size on performance is analyzed? Including these methods would provide a more complete comparison of their scalability and efficiency relative to other approaches.", "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": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: The paper re-visits classical probabilistic linear models (multi-input, multi-output) which can naturally identify relevance of model parameters. Specifically, a method is devised that promotes sparsity in the linear coefficients while capturing correlations among output variables. 

The reviewers were consistently positive about the paper and recommended to accept it. The method seems sound and compares favorably to existing baselines.

The reviewers raised some concerns which have been discussed in the rebuttal phase and which shall be addressed by the authors in the revision of the paper.
Decision: Accept (poster)