Efficient Network Automatic Relevance Determination

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.

Dataset diversity

We 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.

For 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.

For 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.

Through these 2 experimental datasets, we have demonstrated the effectiveness of our method in both environmental and financial domains.

Since 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.

Table R4: Associations of A-shares stocks. ($m=3032, d=60640, N=1696$)

Interpretability in biological applications

Sparsity-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.

Question about Table 2

We 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.

Linear assumptions

We 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.

Numerical stability

We 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.

Discussion about Hybrid NARD

We have included further analysis on this issue in our rebuttal to Reviewer AaWq, which provides a more detailed discussion.

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.

Baseline comparison

As 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.

In 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.

Table R2: Performance Comparison of Various Methods.

We also include the results of experiments on aging phenotype data as shown in Table R3.

Table R3: Associations under different algorithms.

[1] Covariate-adjusted precision matrix estimation with an application in genetical genomics, Biometrika 2013.

[2] A generalized likelihood based Bayesian approach for scalable joint regression and covariance selection in high dimensions, Statistics and computing 2022.

Non-linearity

In 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.

Additionally, our method can be naturally extended to address nonlinearity through kernel method. Specifically, we consider the model

where $\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.

To 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.

As 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.

Numerical stability in Surrogate NARD

We 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.

Thank you for your thoughtful and constructive feedback on our submission.

Edge case of Theorem 3.1

Recall 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.

Tuning the $\lambda_{\text{glasso}}$ Parameter in the ARD Framework

To 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:

Here $\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. The 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.

Discussion about Hybrid NARD

In 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.

We 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.

We appreciate your suggestion and will explore these scenarios further in future work.

Table R1: Performance Comparison of Our Methods on Small datasets.