Rebuttal for fwTq:

We greatly appreciate your positive feedback along with your valuable questions and suggestions. We hope that our responses below maintain your positive assessment.

Answer to Question 1.

Thank you for your question. For graph-based works, the predominant choice of bias metric is demographic parity (DP). Thus, we approach the nascent task of graphical model estimation with a familiar bias metric to verify our approach with established measurements. However, our formulation is suited to other bias metrics such as equalized odds (EO), defined in the global response. While both DP and EO are popular fairness definitions, we cannot compute EO for the true precision matrix since it is conditioned on the ground truth connections. As the crux of our theoretical results relies on measuring the bias in the true precision matrix, we choose DP as our bias metric.

Answer to Question 2.

We present a brief analysis for bias metric $H\_\\mathrm{node}(\\mathbf{\\Theta})$, but a similar analysis holds for other penalties such as $H(\\mathbf{\\Theta})$. The first step of Algorithm 1 is a proximal gradient step. Computing the gradient requires an inverse $( \\mathbf{\\Theta} + \\epsilon \\mathbf{I} )^{-1}$ and product $\\mathbf{A} \\mathbf{\\Theta}\_{\\bar{\\mathcal{D}}}$, both incurring $\\mathcal{O}(p^3)$ operations. The gradient step and soft-thresholding enjoy entry-wise computations with complexities $\\mathcal{O}(p^2)$. The projection step onto the set of positive semidefinite matrices involves an eigendecomposition of $\\dot{\\mathbf{\\Theta}}^{(k)}$ with complexity $\\mathcal{O}(p^3)$. Finally, the step size update $t^{(k)}$ only requires scalar operations, and the accelerated update of $\\check{\\mathbf{\\Theta}}^{(k+1)}$ involves $\\mathcal{O}(p^2)$ operations, so they can be neglected.

If accepted, we will include these additional details in the revised version of the manuscript to strengthen its clarity.

Answer to Question 3.

In Table 2 of our manuscript, we estimate real-world networks from real data, such as the discrete graph signals of the School and Friendship social networks. Thus, we empirically observe satisfactory performance even when observations are non-Gaussian, as queried by the reviewer. However, you also ask an important question: can Fair GLASSO be extended to non-Gaussian graphical models? Gaussianity specifies the loss in the objective of Fair GLASSO, but we may consider other distributions in the optimization problem, such as the Ising negative log-likelihood. Such a substitution requires altering our theoretical results, and in future work we will explore the fairness-accuracy relationship for non-Gaussian distributions.

Answer to Question 4.

We thank you for this question. If accepted, we will provide a detailed discussion on the assumptions for Fair GLASSO. Here, we share a brief summary.

First, AS1 merely defines the cardinality of the support of the true precision matrix.

For the spectrum of the true covariance matrix $\\mathbf{\\Sigma}\_0$, eigenvalues with finite magnitudes as in AS3 is reasonable for any practical setting. AS2 may be violated for a rank-deficient $\\mathbf{\\Sigma}\_0$. However, letting $\\epsilon>0$ addresses both theoretical and implementation concerns, where $\\mathbf{\\Sigma}\_0$ need not be positive definite for convergence or the error bound. As would be expected, the error bound becomes perturbed based on the magnitude of $\\epsilon>0$.

The final assumption can be relaxed to require only that no group vanishes as $p\\rightarrow\\infty$. If a group vanishes, then edges cannot achieve perfect balance across all pairs of groups, which will result in a lower bound on the second term of the error upper bound.

Moreover, we present additional simulations in the attached document illustrating how the assumptions affect performance.

Answer to Question 5.

Our metrics for error and bias are described more thoroughly in Appendix G. We employ the normalized squared Frobenius error, a standard metric for network inference works. We follow similar intuition for bias, where we normalize measurements to compare biases across networks without being affected by changes in graph size or edge weights across graphs.

Answer to Question 6.

In real-world social network analysis, common network characteristics such as homophily can lead to negative outcomes across gender in both social and academic settings [A4]. Hence, the School, Friendship, and Co-authorship networks require scrutiny with respect to fairness, where gender is a critical consideration. For demonstrative purposes, we use publication type as the sensitive attribute for Co-authorship as an example of data with more than two groups. Additionally, biases in recommendation systems can reproduce and even exacerbate existing harmful stereotypes [A5]. The MovieLens dataset, a common benchmark for fair graph machine learning, exemplifies our ability to form unbiased models from networks used for recommendation systems. If our work is accepted, we will include this relevant discussion in the final version.

Answer to Question 7.

Thank you for asking about the reasons for our choice. We followed a common approach of numbering equations to which we refer, which serves to provide lighter notation. However, we agree with the reviewer that numbering all equations may better facilitate the understanding of our results. If the reviewer deems it necessary, we will gladly make this change upon acceptance of this work.

We also appreciate your valuable suggestion regarding further discussion of ethical implications, and we will augment our broader impact discussion if accepted.

Rebuttal for V8H3:

We thank you for your thorough review and detailed questions. We are grateful for your kind words regarding the strengths of our work.

Answer to Weakness 1.

Group-wise modularity is computed as in reference [19] of the paper, that is,

Q(\\mathbf{\\Theta}) = \\sum\_{a=1}^g \\frac{\\mathbf{z}\_a^\\top \\mathbf{\\Theta}\_{\\bar{\\mathcal{D}}} \\mathbf{z} }{2s} - \\sum\_{a=1}^g \\left(\\frac{\\mathbf{z}\_a^\\top \\mathbf{\\Theta}\_{\\bar{\\mathcal{D}}} \\mathbf{z} }{2s}\\right)^2,

where $s$ denotes the number of nonzero entries in $\\mathbf{\\Theta}$. To estimate partial correlation, we apply graphical lasso without a bias penalty, and entries of the resultant precision matrix denote estimates of partial correlation for every pair of variables.

We thank you for your comment. We agree that these definitions are necessary for completeness. We will add more detailed versions of these definitions in the revised paper.

Answer to Weakness 2.

Your point is well taken, and we clarify the meaning of the theoretical result. The result in Theorem 1 holds with high probability as either $n$ or $p$ increase; it is the probability of the error bound holding that converges to 1 as $p$ increases to infinity. We understand that this presentation appears vague. Thanks to your question, we will state this fact more clearly in the final version if accepted.

Answer to Weakness 3.

Thank you for this comment. For lack of space, we outline the effects of violating our assumptions here, but we will provide a comprehensive discussion in the updated paper if accepted. We also provide empirical demonstrations of the assumptions on estimation performance in the attached document, described in the global response.

Observe that since AS1 has no restrictions on permissible values of $s$, AS1 merely defines the number of edges in the true graphical model.

AS2 and AS3 bound the eigenvalues of the true covariance matrix $\\mathbf{\\Sigma}_0$. In realistic settings, these eigenvalues will have finite magnitudes as in AS3, but the covariance may indeed be rank deficient, violating AS2.
Theoretically, this affects the log-determinant bound in equation (18) of our paper, as we cannot finitely bound the difference above with the reciprocal of the smallest eigenvalue. However, in such a case, we may instead consider Theorem 1 with $\\epsilon > 0$. We then obtain a similar error bound, albeit perturbed by the value of $\\epsilon$, that is, the magnitude of $\\epsilon$ increases with the number of zero-valued eigenvalues of $\\mathbf{\\Sigma}_0$.

Finally, AS4 can be relaxed to require only that group sizes be asymptotically similar, so no group vanishes as $p\\rightarrow\\infty$. If a group vanishes, then the error bound will contain persistent terms corresponding to average edge magnitudes of the remaining groups, which cannot be balanced in connection across all groups.

Answer to Weakness 4.

As in our response to your Question 2, we note that our theoretical results do not require that $n$ be fixed. Moreover, you are correct that with insufficient samples, the empirical covariance matrix may not be invertible. Indeed, a major advantage of graphical lasso is that it is suitable in a low sample regime. The sparsity penalty implements prior assumptions of parsimonious entries in the precision matrix to supplement inadequate information in the sample covariance matrix.

Answer to Weakness 5.

We appreciate your valuable suggestion. We propose the following approach to show conditions on parameters $\\mu_1$ and $\\mu_2$ to guarantee Pareto optimality of Fair GLASSO. Say that we assume that our Fair GLASSO estimate $\\mathbf{\\Theta}^*$ has error $\\|\\mathbf{\\Theta}^* - \\mathbf{\\Theta}_0\\|_F^2 = \\delta$ for some $\\delta > 0$. For any feasible precision matrix $\\mathbf{\\Theta} \\in \\mathcal{M}$ such that $\\|\\mathbf{\\Theta} - \\mathbf{\\Theta}_0\\|_F^2 \\leq \\delta$, we wish to determine the selection of $\\mu_1$ and $\\mu_2$ that guarantees $H(\\mathbf{\\Theta}^*) < H(\\mathbf{\\Theta})$, that is, that we have a Pareto optimal estimate $\\mathbf{\\Theta}^*$. Our proposed approach is to exploit the optimality of $\\mathbf{\\Theta}^*$, which yields

$$\\mu\_2 \\left( H(\\mathbf{\\Theta}) - H(\\mathbf{\\Theta}^*) \\right) \\geq \\mathrm{tr}(\\hat{\\mathbf{\\Sigma}}(\\mathbf{\\Theta}^*-\\mathbf{\\Theta})) + \\log \\det \\mathbf{\\Theta} - \\log \\det \\mathbf{\\Theta}^* + \\mu\_1 \\left( \\| \\mathbf{\\Theta}^*\_{\\bar{\\mathcal{D}}} \\|\_1 - \\| \\mathbf{\\Theta}\_{\\bar{\\mathcal{D}}} \\|\_1 \\right),$$

from which we derive bounds on $\\mu_1$ and $\\mu_2$ to ensure that the right-hand side is positive. If our paper is accepted, we will provide the full version of the above Pareto optimality result, as requested by the reviewer.

We appreciate your high standards regarding our paper. We will indeed elaborate on the limitations of our proposed work, the intuition of which your questions have helped us to develop.

Rebuttal for QDSi:

We sincerely thank you for your review and your kind words. We appreciate your insightful comments and how you link our approach to existing works.

Answer to Question 1.

Thank you for your interesting question. Your observation is key; similar to group lasso or joint graph inference methods, we optimize group-wise submatrices of the precision matrix, but our formulation differs subtly with important consequences for fairness.

Joint graph inference methods such as [R1], [R4], and [R5] assume each graph lies on the same node set, where we may treat each submatrix of $\\mathbf{\\Theta}$ for every group pair as a different graph. Our topological fairness only requires that the support and signs of edges are balanced on average across all group pairs. However, joint graph inference yields similar sparsity patterns for all submatrices, which is far more restrictive than balancing edges in expectation. Moreover, the group lasso penalties in [R3], [R4], and [R5] can achieve fairness in support, but they ignore the signs of the edges, thus they cannot balance correlation biases. Thus, existing group lasso or joint graph learning approaches are either too restrictive or do not consider both kinds of graphical model bias. In contrast, our fairness metrics can balance structure in both sparsity patterns and signed edges.

Answer to Question 2.

The main differences between our approach and other graphical model FISTA algorithms lie in (i) the fairness penalty in the objective function; (ii) the constraints for positive semidefinite precision matrices; and (iii) the type of convergence analysis. We first note that our Fair GLASSO algorithm and the work in [R2] have different objective functions and thus different gradient updates. Moreover, while [R2] solves an unconstrained optimization problem, we require projection onto the feasible set after soft-thresholding. On top of this, we guarantee convergence of the optimization variable, which is stronger than guaranteeing convergence of the objective function, as is done for [R2] and classical FISTA algorithms.

To show that the Hessian of $f$ is positive semidefinite (PSD), here we supplement the proof of Theorem 2 (Appendix F). Note that $f$ can be split as $f = f\_1(\\mathbf{\\Theta}) + R\_H(\\mathbf{\\Theta})$, and recall that $\\mathbf{\\Theta}$ is PSD and $\\epsilon > 0$. Then, from (39), it follows that the Hessian of $f\_1(\\mathbf{\\Theta})$ is PSD since it is given by the Kronecker product of two PSD matrices. Next, note that the terms $R\_H(\\mathbf{\\Theta})$ promoting fairness are convex, so their Hessian is also PSD, rendering the Hessian of $f$ PSD.

Finally, the convergence rate in Theorem 2 depends on the Lipschitz constant $L$ and the strong convexity constant $\\alpha$, respectively associated with the largest and smallest eigenvalue of the Hessian of $f$. Thus, it follows that increasing the smallest eigenvalues or decreasing the largest eigenvalues of the Hessian will improve the rate of convergence, and the converse is also true.

Based on the feedback provided, we will highlight these relevant distinctions in the revised manuscript.

Answer to Question 3.

You are correct that the proof is similar to that of the original result in [32]. Since our proof is not limited by space, we aim to provide a self-contained result with enumerated steps for clarity, which can increase accessibility for audiences outside of statistics. We are also careful to identify the effects of bias mitigation. For example, the trace and sparsity differences do adhere to the proof in [32]. However, the log-determinant difference in (18) of our paper deviates slightly from the original due to the additional bias term in the right-hand side of equation (13) in our manuscript versus equation (8) in [32]. Moreover, our goal differs from [R1], [R3], [R4], and [R5], which show the effects on error when graphical lasso is modified to improve performance, while we formalize the tradeoff between accurate estimation versus unbiased solutions.

Answer to Question 4.

For Karate club, School, Friendship, and Co-authorship, nodes represent individuals, while MovieLens nodes represent movies. Regarding samples, the Karate club dataset is the only one without real graph signals, so we generate synthetic Gaussian data given the ground truth network. MovieLens graph signals represent movie ratings. Edges in the Co-authorship network denote collaborations, and we use keyword frequencies as signals. Finally, for both the School and Friendship datasets, the edges represent pairwise student interactions over all time, and the graph signals are sums of interactions over windows of time, that is, the graph signals represent time-varying node degrees as interactions vary.

Detailed descriptions of the real datasets can be found in the Appendix, and we will augment Table 3 to include the interpretation of nodes and signals.

Answer to Question 5.

Your concern is warranted since real data simulations validate our approach in practical scenarios. Social network analysis is a common application of graphical models, but often no ground truth graph exists. However, we must confirm the effectiveness of Fair GLASSO. The Friendship and School datasets thus demonstrate our method for realistic applications, which we can verify with ground truth graphs.

Critically, we emphasize that the Karate club dataset is the only real network for which we generate synthetic graph data. All other real-world datasets are each accompanied by a set of real graph signals. For example, the social networks possess real discrete graph signals, yet we still observe satisfactory performance in Table 2 of our paper. We hope this assuages your concern about the realism of our simulations.

We sincerely thank you for your feedback, and we will refer to your thoughtful comments to clarify our method upon updating the manuscript.

I appreciate the authors' response. Below are my replies to your response.

1. In my opinion, the fair Glasso presented in this paper is indeed a variant of a graphical model with a group penalty, and further investigation with different choices of penalty formulations is needed. I also respectfully disagree with the statement that "group lasso penalties can achieve fairness in support, but they ignore the signs of the edges, thus they cannot balance correlation biases" without providing experimental or theoretical results. Indeed, there are variants of group penalties that show improvement in support recovery for sign-coherent groups. For example, (https://arxiv.org/pdf/1103.2697).

2. I believe an extensive discussion on the novelty of the method in the introduction, regarding both the penalty and the optimization method (FISTA for GLasso), is needed. For example, the phrases "For this purpose, we propose ..." and "we also propose a stronger alternative metric" in Lines 127 and 134 should be rephrased to explicitly mention that these penalties are from previous works such as [15, 18]. Further, the use of FISTA for smooth + $\ell_1$ objectives, as well as graphical lasso, is standard, and the convergence analysis in Theorem 2 follows from existing literature.

3. In my opinion, having no space limits does not allow us to repeat others' proofs. I recommend that the authors remove the identical parts (e.g., trace difference, log-determinant difference, and sparsity penalties) and frequently reference the modified parts of the estimation proof (e.g., from Theorem 1 of [32]).

4. I still did not receive a response to my question, "What is the data matrix $X$ for the Karate Club, School, and Friendship datasets?" Is this data matrix generated using the ground truth graph, or is it available as real data for graph estimation? I don't think the statement "For both the School and Friendship datasets, the edges represent pairwise student interactions ..." answers my question. If I understand correctly, you are using this edge information to construct the ground truth graph, but it is not clear to me what $X$ in Eq. (4) represents for these datasets and how you constructed it. Also, I am uncertain whether these datasets (Karate Club, School, and Friendship datasets) are suitable for Gaussian graphical models, or whether another form of graph learning method should be applied. Do we have good literature on applying GLasso to these datasets, or should we consider using another type of graph learning?

Rebuttal for pQPd:

We are grateful for your positive review and clear questions. Indeed, your review helps us clarify the utility of our approach beyond conceptual discussions. We are glad to hear that you find our work novel, thorough, and comprehensive.

Answer to Question 1.

Thank you for your question. Indeed, as with many graphical lasso modifications, we require selection of appropriate weights for additional penalties. Figures 1 and 2 of the attached document present Fair GLASSO performance as hyperparameters $\\mu_1$ and $\\mu_2$ vary when estimating a fair or unfair precision matrix, respectively. Observe that when the true precision matrix is unfair, increasing $\\mu_2$ encourages a fairer estimate and thus increases the error. Moreover, smaller values of $\\mu_2$ yield greater bias for the unfair setting in Figure 2 than the fair setting in Figure 1. While a larger $\\mu_2$ decreases the bias in both settings, the effect is greater in Figure 2 for a true precision matrix that is unfair. As hyperparameter tuning is a highly practical consideration, we will add these results to the final version if accepted.

Answer to Question 2.

Your question allows us to clarify the flexibility of our bias metrics for graphical model estimation. Indeed, Gaussianity is widely used for its ubiquity and appealing statistical properties, yielding copious theoretical guarantees that we exploit for our results. However, note that Gaussianity only dictates the loss to be minimized, thus any optimization-based graph learning approach is amenable to our fairness penalties. Moreover, our FISTA algorithm is still applicable under other distributions as long as the associated loss is convex and differentiable, such as the negative log-likelihood of the Ising model that you mentioned.

Answer to Question 3.

You raise an excellent point; notions of fairness may differ by application, thus consideration of other definitions is critical. Since fair graphical model estimation is a novel task, we follow the custom of existing graph fairness works that employ demographic parity (DP). However, the penalty for Fair GLASSO is flexible and allows for other bias metrics. Moreover, as long as the chosen fairness penalty is convex and differentiable, our convergent FISTA algorithm remains applicable. Consider for instance the following adaptation of equalized odds (EO),

\\mathbb{P}[ \\Theta^*\_{ij} \~|\~ [\\mathbf{\\Theta}\_0]\_{ij},Z\_{ia}=1, Z\_{ja}=1 ] = \\mathbb{P}[ \\Theta^*\_{ij} \~|\~ [\\mathbf{\\Theta}\_0]\_{ij},Z\_{ia}=1, Z\_{jb}=1 ] ~\\forall a,b\\in[g],

where $\\mathbf{\\Theta}^*$ denotes an estimate of the true precision matrix $\\mathbf{\\Theta}\_0$, and $\\mathbf{Z}$ is the group membership indicator matrix. Note that graphical EO is conditioned on the true precision matrix. Thus, we can measure biases in estimated precision matrices but not the true one. For this reason, we emphasize DP for group fairness since a measure of bias in the true precision matrix is critical to our theoretical interpretation of the fairness-accuracy tradeoff.

Answer to Question 4.

While we show estimation of graphs on the order of 1000 nodes, the complexity indeed scales more than linearly with the number of optimization variables, thus our approach may not be viable for millions of nodes. However, existing works can efficiently estimate very large graphical models, potentially under additional assumptions [A1,A2,A3]. We can combine these approaches with our proposed fairness metrics, although we may lose our guarantee of convergence.

Answer to Weakness 3.

Finally, we wish to clarify a concern you expressed in Weakness 3. We stress that Karate club is the only real network for which we generate synthetic graph signals; all remaining real-world datasets possess data that are not synthetically generated. Indeed, as you noted, we wish to demonstrate the viability of Fair GLASSO in realistic settings, such as social networks with non-Gaussian data.

We thank the reviewer for your questions. You pinpointed vital practical concerns that are crucial to address in order to continue exploring fair graph learning in realistic scenarios.