Fairness-Aware Estimation of Graphical Models

Authors: zhuoping  zhou
davoud ataee tarzanagh
bojian  hou
qi  long
li  shen
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 05:40:52

Abstract: This paper examines the issue of fairness in the estimation of graphical models (GMs), particularly Gaussian, Covariance, and Ising models. These models play a vital role in understanding complex relationships in high-dimensional data. However, standard GMs can result in biased outcomes, especially when the underlying data involves sensitive characteristics or protected groups. To address this, we introduce a comprehensive framework designed to reduce bias in the estimation of GMs related to protected attributes. Our approach involves the integration of the pairwise graph disparity error and a tailored loss function into a nonsmooth multi-objective optimization problem, striving to achieve fairness across different sensitive groups while maintaining the effectiveness of the GMs. Experimental evaluations on synthetic and real-world datasets demonstrate that our framework effectively mitigates bias without undermining GMs' performance.

Main paper: 
 introduction: 
 Graphical models (GMs) are probabilistic models that use graphs to represent dependencies between random variables [34] . They are essential in domains such as gene expression [91] , social networks [17] , computer vision [33] , and recommendation systems [8] . The capacity of GMs to handle complex dependencies makes them crucial across various data-intensive disciplines. Therefore, as our society's reliance on machine learning grows, ensuring the fairness of these models becomes increasingly paramount; see Section 1.1 for further discussions. While significant research has addressed fairness in supervised learning [29] , the domain of unsupervised learning, particularly in the estimation of GMs, remains less explored.We address the fair estimation of sparse GMs where the number of variables P is much larger than the number of observations N [22, 16, 43] . We focus on three types of GMs: I. Gaussian Graphical Model: Rows X 1: , . . . , X N : in the data matrix X ∈ R N ×P are i.i.d. from a multivariate Gaussian distribution N (0, Σ). The conditional independence graph is determined by the sparsity of the inverse covariance matrix Σ -1 , where (Σ -1 ) jj ′ = 0 indicates conditional independence between the jth and j ′ th variables. II. Gaussian Covariance Graph Model: Rows X 1: , . . . , X N : are i.i.d. from N (0, Σ). The marginal independence graph is determined by the sparsity of the covariance matrix Σ, where Σ jj ′ = 0 indicates marginal independence between the jth and j ′ th variables. III. Binary Ising Graphical Model: Rows X 1: , . . . , X N : are binary vectors and i.i.d. with p(x; Θ) = (Z(Θ))-1 expEQUATIONHere, Θ is a symmetric matrix, and Z(Θ) normalizes the density. θ jj ′ = 0 indicates conditional independence between the jth and j ′ th variables. The sparsity pattern of Θ reflects the conditional independence graph. In a data matrix X ∈ R N ×P , each column corresponds to a node in a graph G = (V, E), where V = {1, 2, . . . , P } are vertices and E ⊆ V × V are edges. Column X :i (i ∈ {1, 2, . . . , P }) is a vector of length N , representing the observations for the i-th variable across all N samples. The graph G, represented by the symmetric matrix Θ, has nonzero entries indicating edges and reflects the graph's independence properties. To obtain a sparse and interpretable graph estimate, we consider minimize Θ L (Θ; X) + λ∥Θ∥ 1 subj. to Θ ∈ M.Here, L is a loss function; λ∥ • ∥ 1 is the ℓ 1 -norm regularization with parameter λ > 0; and M is a convex constraint subset of R P ×P . For example, in a Gaussian GM, L(Θ; X) = -log det(Θ) + trace(SΘ), where S = n -1 n i=1 X ⊤ i: X i: , and M is the set of P × P positive definite matrices.Our motivations for obtaining a fair GM estimation are summarized as follows. i) Equitable Representation: Standard group-specific GM models may improve accuracy for targeted groups but do not ensure fairness and can reinforce biases present in the data [48] . A unified approach considering the entire dataset is essential for mitigating biases and promoting fairness across all groups. ii) Legal and Ethical Compliance: Ethical and legal considerations [12] require explicit consent for processing sensitive attributes in model selection. Thus, constructing a fair estimation approach that adheres to fairness practices, uses data with consent, and excludes sensitive attribute information during deployment ensures privacy and legal compliance. iii) Generalization across Groups: A unified fair GM captures differences across groups without segregating the model, enhancing generalizability and preventing overfitting to a specific group [30] , a risk in training separate models for each group.For further discussion, we compare a GM with its proposed Fair variant, as illustrated in Figure 1 . Panel (a) shows the entire dataset, divided into Group Blue and Group Orange in panels (b) and (c). Panels (d) and (e) detail the GM for each group, highlighting variable relationships. Panel (f) demonstrates a conventional GM applied to the full dataset, revealing a bias towards Group Blue. Panel (g) introduces a Fair GM, including modifications (red dashed lines) to reduce bias and ensure balanced representation. These adjustments correct relationships within the model, promoting fairness by preventing disproportionate favor towards any group. This illustration highlights the bias challenge in GMs and the steps Fair GMs take to ensure fair and equal modeling outcomes across groups.Our contributions are summarized as follows:⋄ We propose a framework to mitigate bias in Gaussian, Covariance, and Ising models related to protected attributes. This is achieved by incorporating pairwise graph disparity error and a tailored loss function into a nonsmooth multi-objective optimization problem, striving to achieve fairness across different sensitive groups while preserving GMs performance. ⋄ We develop a proximal gradient method with non-asymptotic convergence guarantees for nonsmooth multi-objective optimization, applicable to Gaussian, Covariance, and Ising models (Theorems 6-8). To our knowledge, this is the first work providing a multi-objective proximal gradient method for GM estimation, in contrast to existing single-objective GM methods [3, 87, 13] . ⋄ We provide extensive experiments to validate the effectiveness of our GM framework in mitigating bias while maintaining model performance on synthetic data, the Credit Dataset, the Cancer Genome Atlas Dataset, Alzheimer's Disease Neuroimaging Initiative (ADNI), and the binary LFM-1b Dataset for recommender systems 3 . 
related work: 
 Estimation of Graphical Models. The estimation of network structures from high-dimensional data [89, 51, 1, 84, 19, 84] is a well-explored domain with significant applications in biomedical and social sciences [44, 59, 25] . Given the challenge of parameter estimation with limited samples, sparsity is imposed via regularization, commonly through an ℓ 1 penalty to encourage sparse network structures [22, 37, 25] . However, these approaches may overlook the complexity of real-world networks, which often have varying structures across scales, including densely connected subgraphs or communities [13, 26, 23] . Recent work extends beyond simple sparsity to estimate hidden communities within networks, reflecting homogeneity within and heterogeneity between communities [47] . This includes inferring connectivity and performing graph estimation when community information is known, as well as considering these tasks in the context of heterogeneous observations [42, 80, 24] .Fairness. Fairness research in machine learning has predominantly focused on supervised methods [11, 4, 15, 39, 92, 79, 28] . Our work broadens this scope to unsupervised learning, incorporating insights from [65, 77, 56, 9, 10] . Notably, [41] has developed algorithms for fair clustering using the Laplacian matrix. Our approach diverges by not presupposing any graph and Laplacian structures. The most relevant works to this study are [78, 93, 52, 53] . Specifically, [78] initiated the learning of fair GMs using an ℓ 1 -regularized pseudo-likelihood method for joint GMs estimation and fair community detection. [93, 94] proposed a fair spectral clustering model that integrates graph construction, fair spectral embedding, and discretization into a single objective function. Unlike these models, which assume community structures, our study formulates fair GMs without such assumption. Concurrently with this work, [52] proposed a regularization method for fair Gaussian GMs assuming the availability of node attributes. Their methodology significantly differs from ours, as we focus on developing three classes of fair GMs (Gaussian, Covariance, and Ising models) for imbalanced groups without node attributes, aiming to automatically ensure fairness through non-smooth multi-objective optimization.Σ 1 is symmetric, it is adjusted to (Σ 1 + Σ ⊤ 1 )/2.To ensure it is positive definite, we further adjust Σ 1 as:EQUATIONwhere λj represents max(λ j (Σ 1 ), 10 -5 ), and v j is the corresponding eigenvector. 3. For each subsequent group (k = 2, . . . , K), the covariance matrix Σ k is initially set equal to Σ k-1 . Then, two (one for sensitivity analysis) of its sub-matrices, which have not been altered yet, are reset to the identity matrix. 
fair estimation of graphical models: 
 Notation. R d denotes the d-dimensional real space, and R d + and R d ++ its positive and negative orthants. Vectors and matrices are in bold lowercase and uppercase letters (e.g., a, A), with elements a i and a ij . Rows and columns of A are A i: and A :j , respectively. For symmetric A, A ≻ 0 and A ⪰ 0 denote positive definiteness and semi-definiteness. Λ i (A) is the ith smallest eigenvalue of A. The matrix norms are defined as∥A∥ 1 = ij |a ij | and ∥A∥ F = ( ij |a ij | 2 ) 1/2. For any positive integer n, [n] := {1, . . . , n}. Any notation is defined upon its first use and summarized in Table 3 .To evaluate the effects of joint GMs learning on different groups, we compare models trained on group-specific data with those trained on a combined dataset. Let a dataset X be divided into K sensitive groups, with the data for group k ∈ [K] represented as X k ∈ R N k ×P , where N k is the sample size for group k, and N = n k=1 N k . The performance of a GM, denoted by Θ, for group k is measured by the loss function L(Θ; X k ). Our goal is to find a global model Θ * that minimizes performance discrepancies across groups. We define graph disparity error to quantify fairness: Definition 1 (Graph Disparity Error). Given a dataset X ∈ R N ×P with K sensitive groups, where X k represents the data for group k ∈EQUATIONThe graph disparity error for group k is then:EQUATIONThis measures the loss difference between a global graph matrix Θ and the optimal local graph matrix Θ * k for each group's data X k . A fair GM, under Definition 1, seeks to balance E k across all groups.Algorithm 1 Fair Estimation of GMs (Fair GMs) Require: Data Matrix X = X 1 ∪ X 2 ∪ • • • ∪ X K ; Parameters λ > 0, ϵ > 0, T > 0,EQUATIONTo address the imbalance in graph disparity error among all groups, we introduce the idea of pairwise graph disparity error, which quantifies the variation in graph disparity between different groups. Definition 3 (Pairwise Graph Disparity Error). Let ϕ : R → R + be a penalty function such as ϕ(x) = exp(x) or ϕ(x) = 1 2 x 2 . The pairwise graph disparity error for the group k is defined asEQUATIONThe motivation for Definition 3 follows from the work of [35, 65, 95] in PCA and CCA. In our convergence analysis, we focus on smooth functions ϕ, such as squared or exponential functions, while nonsmooth choices, such as ϕ(x) = |x|, can be explored in the experimental evaluations.This section introduces a framework designed to mitigate bias in GMs (including Gaussian, Covariance, and Ising) related to protected attributes by incorporating pairwise graph disparity error into a nonsmooth multi-objective optimization problem. Smooth multi-objective optimization tackles fairness challenges in unsupervised learning [35, 95] , proving particularly useful when decision-making involves multiple conflicting objectives.We use non-smooth multi-objective optimization to balance two key factors: the loss in GMs and the pairwise graph disparity errors. To achieve this, letEQUATIONwhere g (Θ) := λ∥Θ∥ 1 for some λ > 0.Consequently, we propose the following multi-objective optimization problem for Fair GMs:EQUATIONHere, M is a convex constraint subset of R P ×P andF : Ω → R M is a multi-objective function. Assumption A. For some L > 0, all Θ, Φ ∈ M, k ∈ [M ], ∥∇f k (Φ)-∇f k (Θ) ∥ F ≤ L∥Φ-Θ∥ F .Note that Assumption A holds for smooth ϕ functions such as squared or exponential, as specified in Definition 6, and when L is a smooth loss function. We demonstrate in Appendix C that this assumption holds for the Gaussian, Covariance, and Ising models studied in this work. To proceed, we provide the following definitions; see [20, 75, 73] for more details. Definition 4 (Pareto Optimality). In Problem (8), a solution Θ * ∈ M is Pareto optimal if there is noΘ ∈ M such that F(Θ) ⪯ F(Θ * ) and F(Θ) ̸ = F(Θ * ). It is weakly Pareto optimal if there is no Θ ∈ M such that F(Θ) ≺ F(Θ * ).Definition 5 (Pareto Stationary). We define a point Θ ∈ R P ×P as Pareto stationary (or critical) if it satisfies the following condition:max k∈[M ] F ′ k ( Θ; D) := lim α→0 F k ( Θ + αD) -F k ( Θ) α ≥ 0 for all D ∈ R P ×P .To solve Problem (8), we use the proximal gradient method and establish its convergence to a Pareto stationary point for the nonsmooth Problem (8) . The procedure for our fairness-aware GMs (Fair GMs) is detailed in Algorithm 1. Given local graph estimates {Θ * k } K k=1 obtained in S1., and ℓ > L, where L is a Lipschitz constant defined in Assumption A, the update of the global fair graph estimate Θ is produced in S2. by solving:EQUATIONEQUATIONNote that the convexity of g(Θ) = λ∥Θ∥ 1 ensures a unique solution for Problem (9) . We provide a simple yet efficient approach to solve Subproblem (9) through its dual in Appendix B. In addition, Proposition 11 in Appendix B characterizes the weak Pareto optimality for Problem (8).We apply Algorithm 1 to Gaussian, Covariance, and Ising models and provide theoretical guarantees.Fair Graphical Lasso (Fair GLasso). Consider X 1: , . . . , X N : as i.i.d. samples from N (0, Σ). In the GLasso method [22] , the loss is defined as L G (Θ; X) := -log det(Θ) + trace(SΘ), where Θ is constrained to the set M = {Θ : Θ ≻ 0, Θ = Θ ⊤ } and S = n -1 n i=1 X i: X ⊤ i: ∈ R N ×N is the empirical covariance matrix of X. Extending this to fair GLasso and following (8) , the multi-objective optimization problem is formulated as:minimize Θ F(Θ) = [L G (Θ; X) + λ∥Θ∥ 1 , F 2 (Θ), • • • , F M (Θ)] subj. to Θ ∈ M = {Θ : Θ ≻ 0, Θ = Θ ⊤ }.(Fair GLasso)Assumption B. Let N * be the set of weakly Pareto optimal points for (8), and Ω F (α) := {Θ ∈ S | F(Θ) ⪯ α} denote the the level set of F for α ∈ R M . For all Θ ∈ Ω F (F(Θ (0) )), there exists Θ * ∈ N * such that F(Θ * ) ⪯ F(Θ) and R := supF * ∈F(N * ∩Ω F (F(Θ 0 ))) inf Θ∈F -1 ({F * }) ∥Θ -Θ (0) ∥ 2 F < ∞.This assumption is satisfied when Ω F (F(Θ (0) )) is bounded [75, 73] . When M = 1, it holds if the problem has at least one optimal solution. If Ω F (F(Θ (0) )) is bounded, Assumption B also holds, such as when F k is strongly convex for some k ∈ [M ].Theorem 6. Suppose Assumptions A and B hold. Let {Θ (t) } t≥1 be the sequence generated by Algorithm 1 for solving (Fair GLasso). Then,sup Θ∈M min k∈[M ] F k Θ (t) -F k (Θ) ≤ ℓR 2t .Fair Covariance Graph (Fair CovGraph). For the Fair CovGraph, we assume X 1: , . . . , X N : are i.i.d. samples from N (0, Σ). We use a sparse estimator for the covariance matrix, ensuring it remains positive definite and specifies the marginal independence graph. Following [62] , we define the estimator's loss function as L C (Σ, X) := 1 2 ∥Σ -S∥ 2 F -τ log det(Σ) with τ > 0. Building on this and using (7) and (8) , for some nonegative constants γ C and λ, we introduce the Fair CovGraph optimization problem, formulated as follows:minimize Σ F (Σ) = [F 1 (Σ), F 2 (Σ), . . . , F M (Σ)] subj. to Σ ∈ M = {Σ : Σ ≻ 0, Σ = Σ ⊤ }. (Fair CovGraph)Here, following (8), we havef 1 (Σ) = L C (Σ; X) and f k (Σ) = ∆ k-1 (Σ) for 2 ≤ k ≤ K + 1. Also, F 1 (Σ) = f 1 (Σ) + λ∥Σ∥ 1 and F k (Σ) = f k (Σ) + λ∥Σ∥ 1 + γ C ∥Σ∥ 2 F for 2 ≤ k ≤ M = K + 1.The parameter γ C is used to convexify (Fair CovGraph) and is crucial for ensuring the convergence of Algorithm 1. The following theorem establishes the convergence of Algorithm 1 for (Fair CovGraph). Theorem 7. Under conditions similar to Theorem 6, by replacing Θ with Σ and L G (Θ; X) with L C (Σ, X), for the sequence {Σ (t) } t≥1 generated by Algorithm 1 applied to (Fair CovGraph), and for γ C ≥ max{0, -Λ min (∇ 2 f k (Σ))} for all k ∈ [K], we have:sup Σ∈M min k∈[M ] F k Σ (t) -F k (Σ) ≤ ℓR 2t .Fair Binary Ising Network (Fair BinNet). In this section, we focus on the binary Ising Markov random field as described by [58] . The model considers binary-valued, i.i.d. samples with probability density function defined in (1) . Following [31] , we consider the following loss function:EQUATION)Given some nonegative constants γ I and λ, the Fair BinNet objective is defined as:minimize Θ F (Θ) = [F 1 (Θ) , F 2 (Θ) , . . . , F M (Θ)] subj. to Θ ∈ M = {Θ : Θ = Θ ⊤ }. (Fair BinNet)Here, following (8), we have f 1 (Θ) = L I (Θ; X), andf k (Θ) = ∆ k-1 (Θ) for 2 ≤ k ≤ K +1. Also, F 1 (Θ) = f 1 (Θ) + λ∥Θ∥ 1 , and F k (Θ) = f k (Θ) + λ∥Θ∥ 1 + γ I ∥Θ∥ 2 F for 2 ≤ k ≤ M = K + 1.The parameter γ I convexifies Problem (Fair BinNet) and ensures Algorithm 1 converges. The following theorem establishes the convergence of Algorithm 1 for Problem (Fair BinNet). Theorem 8. Suppose Assumptions A and B hold, and thatγ I ≥ max{0, -Λ min (∇ 2 f k (Θ))} for all k ∈ [K].Then, the sequence {Θ (t) } t≥1 generated by Algorithm 1 for (Fair BinNet) satisfiessup Θ∈M min k∈[M ] F k Θ (t) -F k (Θ) ≤ ℓR 2t .Remark 9 (Iteration Complexity of Algorithm 1). Theorems 6, 7, and 8 establish the global convergence rates of O(1/t) for Algorithm 1 for Gaussian, Covariance, and Ising models, respectively. In contrast to Theorem 6, Theorems 7 and 8 necessitate the inclusion of an additional convex regularization term with parameters γ C and γ I , respectively, to achieve Pareto optimality. Remark 10 (Computational Complexity of Algorithm 1). Given the iteration complexity to achieve ϵ-accuracy is O(ϵ -1 ), the overall time complexity of our optimization procedure becomes O ϵ -1 max(N P 2 , P 3 ) . Assuming a small number of groups (K << N, P, 1/ϵ), the complexity aligns with that of standard proximal gradient methods used for covariance and inverse covariance estimation, making it feasible for large N and P . To further support the theoretical analysis, sensitivity analysis experiments are conducted to investigate the impact of varying P , N , K, and group imbalance on the performance of the proposed methods. Note that the complexity of Algorithm 1 applied to (Fair BinNet) depends on the choice of subproblem solver (e.g., first or second order) due to the nonlinearity of (10) . Further experiments and discussions are detailed in Appendices D.6-D.9.PCEE := j,j ′ ∈[P ] 1{ Θjj ′ ≥ λ and |Θ jj ′ | ≥ λ} / j,j ′ ∈[P ] 1{|Θ jj ′ | ≥ λ} ,where 1{•} is the indicator function, Θ and Θ are the groundtruth and estimated graph. 
experimental setup: 
 Baseline. The Iterative Shrinkage-Thresholding Algorithm (ISTA) is widely used for sparse inverse covariance estimation [60] due to its simplicity and efficiency. We adapt ISTA for the Covariance and Ising models and use them as a baseline to compare with our proposed Fair GMs. Note that our Fair GMs reduce to ISTA for Gaussian, Covariance, and Ising models if M = 1 in (8 Evaluation Criteria. In our experiments, we introduce three metrics to evaluate the performance of our methods and the baseline methods:1. Value of the objective function of GM: F 1 := L(Θ; X) + λ∥Θ∥ 1 .2. Summation of pairwise graph disparity error for fairness:∆ := K k=1 ∆ k .In the simulation, we construct two 100 × 100 block diagonal covariance matrices, Σ 1 and Σ 2 (Figures 2c and 2d ). These matrices correspond to two sensitive groups and are created following the rigorous process in Appendix D.2. Each graph has three consistent diagonal blocks, with Group 1 also featuring two distinct blocks indicating bias. For each group, we derive the ground truth graphs by 2a and 2b ). Datasets are generated from normal distributions: 1000 samples from N (0, Σ 1 ) for the first group, and 1000 samples from N (0, Σ 2 ) for the second.Θ 1 = Σ -1 1 and Θ 2 = Σ -1 2 (FiguresResults. Figure 2g shows the global graph derived using Standard GLasso on the entire dataset, where the two top-left blocks are not distinctly marked, suggesting bias towards Θ 2 . In contrast, Figure 2h shows a graph from our method that enhances block visibility, reducing bias. This improvement is supported by the results in Table 1 , where the PCEE difference of Fair GLasso is smaller than that of Standard GLasso. Comparable efficacy in bias reduction for CovGraph is shown in Figures 2c, 2d, 2i , 2j, and Table 1 , demonstrating our methods' effectiveness in achieving fairness.We provided two simulation networks: Θ 1 for Group 1 with P = 50 nodes and three hub nodes, and Θ 2 for Group 2 with one hub node (see Appendix D.3 for details). Adjacency matrices are shown in Figures 2e and 2f . We generate N 1 = 500 and N 2 = 1000 observations via Gibbs sampling, updating each variable x (t+1) j at iteration t + 1 using the Bernoulli distribution:x (t+1) j ∼ Bernoulli(z θ /(1 + z θ )), where z θ = exp(θ jj + j ′ ̸ =j θ jj ′ x (t) j ′ ).The first 10,000 iterations are designated as the burn-in period to ensure statistical independence among observations. Finally, observations are collected at every 100th iteration.Results. Figure 2k illustrates the global graph from Standard BinNet, which is predominantly biased towards Θ 2 by identifying only one hub node. In contrast, Figure 2l , derived from Fair BinNet, presents a more balanced structure. While this improvement might not be visually evident, the quantitative results in Table 1 and Table 2 confirm it. Table 1 reveals that PCEE for Group 1 improved significantly, increasing from 0.4444 to 0.7485. Conversely, PCEE for Group 2 exhibited a decrease from 0.9481 to 0.7662. This convergence in performance metrics between the two groups indicates a more balanced distribution of predictive errors, thus enhancing the overall fairness of the model.We apply GLasso to analyze RNA-Seq data from TCGA, focusing on lung adenocarcinoma. The data includes expression levels of 60,660 genes from 539 patients. From these, 147 KEGG pathway genes [36] are selected to construct a gene regulatory network. GLasso reveals conditional dependencies, aiding in understanding cancer genetics and identifying therapeutic targets. However, initially, this method, without accounting for sex-based differences, risks overlooking critical biological disparities, potentially skewing drug discovery and health outcomes across genders. Therefore, we divide the patient cohort into two groups based on sex: 248 males and 291 females. This stratification enables the use of Fair GLasso, which creates a more equitable gene regulatory network by accounting for these differences. The parameter λ is set to 0.03 for this experiment. Additionally, each variable is normalized to achieve a zero mean and a unit variance.Results. The gene networks identified by GLasso and Fair GLasso are presented in Figures 3a-3b . GLasso identified several hub nodes, including NCOA1, BRCA1, FGF8, AKT1, NOTCH4, and CSNK1A1L. In contrast, Fair GLasso uniquely detected PIK3CD, suggesting its potential relevance in capturing sex-specific differences in lung adenocarcinoma. Although direct evidence linking PIK3CD exclusively to sex-specific traits in cancer is limited, this finding aligns with recent insights into sex-specific regulatory mechanisms in cancer [63, 45] . PIK3CD is a key component of the PI3K/Akt signaling pathway, which is involved in cell regulation and frequently implicated in various malignancies. The identification of PIK3CD by Fair GLasso demonstrates its potential to uncover biologically relevant genes that may be overlooked in conventional analyses, enhancing our understanding of lung adenocarcinoma and facilitating the development of personalized therapies.The performance of GLasso and Fair GLasso is evaluated using AV45 and AV1451 PET imaging data from the Alzheimer's Disease Neuroimaging Initiative (ADNI) [85, 86] , focusing on amyloid-β and tau deposition in the brain. The dataset includes standardized uptake value ratios of AV45 and AV1451 tracers in 68 brain regions, as defined by the Desikan-Killiany atlas [14] , collected from 1,143 participants. An amyloid (or tau) accumulation network [68] is constructed to investigate the pattern of amyloid (or tau) accumulation. GLasso and Fair GLasso are used to uncover conditional dependencies between brain regions, providing insights into Alzheimer's disease progression and identifying potential biomarkers for early diagnosis and treatment response monitoring. To examine the influence of sensitive attributes on the network structure, marital status, and race are incorporated as exemplary sensitive attributes due to their reported association with dementia risk [71, 49] . Comprehensive details regarding the experiments, results, and analysis are provided in Appendix 4.5.The performance of Fair CovGraph is evaluated using the Credit Datasets [90] from the UCI Machine Learning Repository [2] . These datasets have been previously used in research on Fair PCA [55, 83] , which shows potential for improvement through sparse covariance estimation. The dataset composition is detailed in Table 5 in Appendix D.5, with categorizations based on gender, marital status, and education level. For this experiment, the parameters τ and λ are set to 0.01 and 0.1, respectively. Each variable in the dataset is standardized to have a mean of zero and a variance of one. As shown in Table 2 , our Fair CovGraph achieves a 53.75% increase in fairness with only a 0.42% decrease in the graph objective, demonstrating the strong ability of our method to attain fairness. Table 2 : Outcomes in terms of the value of the objective function (F 1 ), the summation of the pairwise graph disparity error (∆), and the average computation time in seconds (± standard deviation) from 10 repeated experiments. "↓" means the smaller, the better, and the best value is in bold. Note that both F 1 and ∆ are deterministic. Fair BinNet provides a more diversified recommendation network.Dataset F 1 ↓ %F 1 ↑ ∆ ↓ %∆ ↑ Runtime ↓ GM Fair GM GMLFM-1b Dataset4 contains over one billion listening events intended for use in recommendation systems [66] . In this experiment, we use the user-artist play counts dataset to construct a recommendation network of artists. Our analysis focuses on 80 artists intersecting the 2016 Billboard Artist 100 and 1,807 randomly selected users who listened to at least 400 songs. We transform the play counts into binary datasets for BinNet models, setting play counts above 0 to 1 and all others to 0.This experiment examines male and female categories, stratifying the dataset into two groups with 1,341 and 466 samples, respectively. We set the BinNet models' parameter, λ, to 1e -5.Results. Figures 3c-3d show the recommendation networks of the 2016 Billboard Top 10 popular music artists based on BinNet's and Fair BinNet's outputs. The comparative analysis reveals that Fair BinNet provides a more diversified recommendation network, particularly for the artist The Weeknd. Enhancing fairness fosters cross-group musical preference exchange, breaks the echo chamber effect, and broadens users' exposure to potentially intriguing music, enhancing the user-friendliness of the music recommendation system.In fairness studies, the trade-off between fairness and model accuracy presents a fundamental challenge. An effective fair method should achieve equitable outcomes while maintaining strong accuracy performance. We evaluate this balance by analyzing the percentage changes in both accuracy and fairness metrics. Specifically, we define these changes as:%F 1 = -F1 of Fair GM -F1 of GM F1 of GM × 100%, and %∆ = -∆ of Fair GM -∆ of GM ∆ of GM × 100%.Our empirical results (Tables 2, 4 , and Figure 4 ) demonstrate that Fair GMs substantially reduce disparity error, thereby improving fairness, while incurring only minimal degradation in the objective function's value. This favorable trade-off validates the effectiveness of our approach. However, Fair GMs do face computational challenges, primarily stemming from two sources: local graph computation and multi-objective optimization. 2 and 4 ). %F 1 is slight, while %∆ changes are substantial, signifying fairness improvement without serious accuracy sacrifice.To address these limitations, we propose several solutions. The local graph learning phase can be accelerated using advanced graphical model algorithms such as QUIC [32] , SQUIC [7] , PISTA [69] , GISTA [60] , OBN [57] , or ALM [67] . Moreover, to mitigate the increased complexity from multiple objectives, we introduce a stochastic objective selection strategy, randomly sampling a subset of objectives in each iteration. This approach effectively reduces computational overhead while maintaining model fairness and performance. To validate these computational considerations, we conducted additional experiments using GLasso, with detailed results presented in the Appendix D.10. 
conclusion: 
 In this paper, we tackle fairness in graphical models (GMs) such as Gaussian, Covariance, and Ising models, which interpret complex relationships in high-dimensional data. Standard GMs exhibit bias with sensitive group data. We propose a framework incorporating a pairwise graph disparity error term and a custom loss function into a nonsmooth multi-objective optimization. This approach enhances fairness without compromising performance, validated by experiments on synthetic and real-world datasets. However, it increases computational complexity and may be sensitive to the choice of loss function and balancing multiple objectives. Future research can include:F1. Integrating our Fair GMs approach with supervised methods for downstream tasks, including spectral clustering [82] , graph regularized dimension reduction [76] .F2. Developing novel group fairness notions based on sensitive attributes within our nonsmooth multi-objective optimization framework.F3. Extending fairness to ordinal data models, which are crucial for socioeconomic and healthrelated applications [27] , neighborhood selection [50] , and partial correlation estimation [40] .Despite some limitations of Fair GMs for larger group sizes, this work demonstrates the potential of nonsmooth multi-objective optimization as a powerful tool for mitigating biases and promoting fairness in high-dimensional graph-based machine learning, contributing to the development of more equitable and responsible AI systems across a wide range of domains. Frobenius norm: 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments? Answer: [Yes] Justification: Yes, we repeated all the experiments 10 times. We provide a full description of this statistical analysis in our implementation details provided in Section 4. Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified.• For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 10. Broader Impacts Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed? Answer: [Yes] Justification: Yes, our work is about fairness in graphical models. Improving fairness in machine learning algorithms will definitively have positive societal impacts, which are discussed in the introduction and experiment sections. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: ( ij |a ij | 2 ) 1/2 [M ]The set {1, 2, . . . , M } Λi (A) ith eigenvalue of A Λ min (A)The smallest eigenvalue of A M αh (x)Moreau envelope: min y h(y)+ 1 2α ∥x -y∥ 2 prox αh (x)Proximal operator: arg min y h(y) + K k=1 ∆ k φ ℓ max k=1,...,M ⟨∇f k (Θ), Φ -Θ⟩ + g(Φ) -g(Θ) + ℓ 2 ∥Φ -Θ∥ 2 FF Objective function of the multi-objective optimization problem F k kth objective in the multi-objective optimization problem P ℓ solutions of the min-max problem for each ℓ:arg min Φ∈M φ ℓ (Φ; Θ) R sup F * ∈F(N * ∩Ω F (F(Θ 0 ))) inf Θ∈F -1 ({F * }) Θ -Θ (0) 2 F M Convex constraint subset of R P ×P C Standard simplex: ρ ∈ R M : M k=1 ρ k = 1, ρ k ∈ [0, 1] ∀k ∈ [M ]B Addendum to Section 3 B.1 Dual Reformulation and Computation of Subproblem (9) In this section, we provide a dual method for solving the Subproblem (9) defined as:EQUATIONfor all ℓ > L where L is defined in Assumption A.For simplicity, letEQUATIONBy considering the standard simplex,EQUATIONwe reformulate (11) asEQUATIONBy leveraging the convexity of M, the compactness and convexity of C, and the convexity-concavity property of M k=1 ρ k ψ k,ℓ (Φ; Θ) with respect to Φ and ρ, respectively, we can invoke Sion's minimax theorem to reformulate (14) as follows:EQUATIONExpanding on the definition of ψ k,ℓ , we arrive at the following expression:EQUATIONwhere Moreau envelope M αh (x) and the proximal operator are defined asEQUATION)EQUATION)Problem ( 16) is equivalent to the dual problem:EQUATIONwhereEQUATIONUpon solving the dual Problem (18) , the optimal solution Φ * of ( 11) is obtained through:EQUATIONIn the implementation, for the given g (Θ) = λ∥Θ∥ 1 , prox 1 ℓ g is computed using soft thresholding η 1 ℓ λ , as shown below:EQUATIONTo provide a clear and logical summary of the iterative update process in Algorithm 1, we proceed as follows: At each iteration t, the update for Θ (t+1) is performed by inputting Θ (t) and solving the Subproblem (11) . This is achieved by utilizing the scipy.optimize.minimize function with the method="trust-constr" option to solve the dual problem. Specifically, for given constants ℓ > L and λ > 0, and the calculated ρ (t) ∈ C at the tth iteration, the update rule for Θ (t+1) is given by: As demonstrated in the analysis of Subproblem (9) in the beginning of this section, the proposition implies that the descent direction is the minimum norm matrix within the convex hull of the gradients of all objectives. Furthermore, this direction is non-increasing with respect to each individual objective function. This property ensures that the chosen descent direction simultaneously minimizes the overall norm while guaranteeing non-increasing behavior for each objective, thereby facilitating the optimization process in a multi-objective setting.EQUATIONTo update Θ (t) in Algorithm 1 applied to (Fair GLasso), the iterative update formula in Equation ( 21) is used at each iteration. The gradients for the functions f 1 and {f k+1 } M -1 k=1 are computed as follows: The gradient of f 1 with respect to Θ is given by:EQUATION)where S is the sample covariance matrix and Θ -1 is the inverse of the precision matrix Θ.The gradient of f k+1 for k = 1, . . . , M -1 with respect to Θ is given by:EQUATIONwhere K is the number of groups, S k and S s are the sample covariance matrices for groups k and s, respectively, and Θ * k and Θ * s are the optimal precision matrices for groups k and s obtained by solving the group-specific GLasso problems.To refine the iterative update rule for estimating the fair covariance matrix Σ in (Fair GLasso) using Algorithm 1, the gradients for f 1 and the set {f k+1 } M -1 k=1 are computed as follows: The gradient of f 1 with respect to Σ is given by:EQUATION)where S is the pooled sample covariance matrix, τ is the regularization parameter, and Σ -1 is the inverse of the covariance matrix Σ.The gradient of f k+1 for k = 1, . . . , M -1 with respect to Σ is given by:EQUATIONwhere K is the number of groups, S k and S s are the sample covariance matrices for groups k and s, respectively, Σ * k and Σ * s are the optimal covariance matrices for groups k and s obtained by solving the group-specific CovGraph problems.To simplify, denotez θ = exp   θ jj + j ′ ̸ =j θ jj ′ x ij ′   , and z ϕ = exp   ϕ jj + j ′ ̸ =j ϕ jj ′ x ij ′   .The gradients of the objectives of Fair BinNet that are utilized in the iterative update formula (21) are computed as follows:EQUATIONHere, L (Θ, X) = f 1 (Θ) is the negative log-likelihood of the Ising model, where Θ ∈ R P ×P is the interaction matrix, X ∈ {0, 1} N ×P is the binary data matrix with N samples and P variables, and θ jj ′ denotes the (j, j ′ )-th element of Θ.Lemma 12. Let {Θ (t) } be generated by Algorithm 1. Then for all k = 1, . . . , M , we haveEQUATIONProof. Let φ ℓ (Φ, Θ) be defined as in (9) . Following the proof of [73, Lemma 4.1], we haveEQUATIONIf ℓ > L, using the descent lemma [5, Proposition A.24], for all k = 1, . . . , M , we obtainEQUATIONSince the right-hand side of the above inequality is less than or equal to zero, it implies thatEQUATION)Lemma 13. Suppose Assumption A holds. Let f k and g have convexity parameters µ k ∈ R + and ν ∈ R + , respectively, and define µ := min k∈[M ] µ k . Then, for all Θ ∈ M, we haveEQUATIONk satisfies the following conditions:1. There exists η (t) ∈ ∂g Θ (t+1) such that- M k=1 ρ (t) k ∇f i Θ (t) + η (t) = ℓ Θ (t+1) -Θ (t) .2. ρ (t) ∈ C where C is defined in (13) .Proof. Assumption A yieldsEQUATIONFrom the convexity of f k and g, we haveEQUATIONCondition 1 and Condition 2 yieldEQUATIONC.2 Proof of Theorem 6 for Fair GLasso First, we present the convexity analysis and gradient Lipschitz continuity. Proposition 14 (Convexity of Fair GLasso). Each f k for k = 1, . . . , M and g defined in (Fair GLasso) of Fair GLasso are convex. Further, f 1 is strongly convex.Proof. First, consider f 1 in the first objective function of Fair GLasso:EQUATIONwhere Θ is a positive definite matrix.The gradient and Hessian of f 1 are, respectively:EQUATIONThe positive definiteness of Θ implies that Θ -1 is also positive definite. Therefore, the Hessian H f1 , being the Kronecker product of Θ -1 with itself, is positive definite. This establishes that the objective function f 1 is strongly convex.Next, the functions f k+1 for k = 1, . . . , M -1 are defined as:EQUATIONwhich are convex due to the linearity of the trace operator in the loss function difference L(Θ; X k ) -L(Θ; X s ) = trace((S k -S s )Θ), leading to a strong convexity parameter of 0.In addition, g(Θ) = λ∥Θ∥ 1 is identified as a closed, proper, and convex function.Proposition 15 (Gradient Lipschitz Continuity of Fair GLasso). The gradients of f k for k = 1, . . . , M defined in (Fair GLasso) are Lipschitz continuous.Proof. First, we present the gradient and Hessian of functions f 1 and {f k+1 } M -1 k=1 as follows:f 1 (Θ) = -log det(Θ) + trace(SΘ), ∇f 1 (Θ) = S -Θ -1 , H f1 = Θ -1 ⊗ Θ -1 ; f k+1 (Θ) = s∈[K],s̸ =k 1 2 (trace(S k Θ) -trace(S s Θ) + log det(Θ * k ) -trace(S k Θ * k ) -log det(Θ * s ) + trace(S s Θ * s )) 2 , ∇f k+1 (Θ) = s∈[K],s̸ =k (trace(S k Θ) -trace(S s Θ) + log det(Θ * k ) -trace(S k Θ * k ) -log det(Θ * s ) + trace(S s Θ * s )) (S k -S s ) , H f k+1 (Θ) = s∈[K],s̸ =k (S k -S s ) ⊗ (S k -S s ) . (38) Define L 1 = Λ max (H f1 ) and L k+1 = Λ max (H f k+1 ) for k = 1, . . . , M -1. Given that {f k } M k=1are convex (as proven in Proposition 14) and twice differentiable, their gradients satisfy Lipschitz continuity with Lipschitz constants {L k } M k=1 .Next, we present the proof for Theorem 6.proof for Theorem 6. From Proposition 14 and Proposition 15, convexity and gradient Lipschitz continuity of objective functions {f k } M k=1 are verified. Hence, Assumption A holds. From Lemma 13 and the convexity of f i and g, for all Θ ∈ M, we obtainEQUATIONAdding up the above inequality (39) from t = 0 to t = t, we haveEQUATIONLemma 12 implies that F k Θ ( t+1) ≤ F k Θ (t+1) for all t ≤ t andEQUATIONLetEQUATIONSince ρt k ≥ 0 andEQUATIONNow, following the proof of [75, Theorem 5.1] and using Assumption B, we obtainEQUATION)EQUATION)EQUATION)The inequalityF k Θ (t) ≤ F k Θ (0) from Lemma 12 implies that sup Θ∈Ω F( F(Θ (0) )) min k∈[M ] F k Θ ( t+1) -F k (Θ) = sup Θ∈Ω F( F(Θ t+1 )) min k∈[M ] F k Θ ( t+1) -F k (Θ) = sup Θ∈M min k∈[M ] F k Θ ( t+1) -F k (Θ) . (47)Moreover, from Assumption B that for all Θ ∈ Ω F F Θ (0) , there exists Θ * ∈ N * such that F (Θ * ) ⪯ F (Θ), it follows:EQUATION)Therefore, from ( 44) and ( 48), we can conclude thatEQUATIONFirst, we present the convexity analysis and gradient Lipschitz continuity. Proposition 16 (Convexity of Fair CovGraph). Through incorporating a convex regularization termγ C ∥Θ∥ 2 F for some γ C ≥ max{0, -Λ min (∇ 2 f k (Θ))} into each f k for k = 2, .. . , M , each f k for k = 1, . . . , M and g defined in the multi-objective optimization problem (Fair CovGraph) of Fair CovGraph are guaranteed to be convex. In particular, f 1 is strongly convex.Proof. In Fair CovGraph (Fair CovGraph), the function f 1 , its gradient, and Hessian are defined as follows:f 1 (Σ) = 1 2 ∥Σ -S∥ 2 F -τ log det (Σ) , ∇f 1 (Σ) = Σ -S -τ Σ -1 , H f1 = I P 2 + τ Σ -1 ⊗ Σ -1 .The positive definiteness of the covariance matrix Σ guarantees that the Hessian matrix H f1 is also positive definite, establishing the strong convexity of the function f 1 . Next, consider:EQUATION)This difference is necessarily convex, as it is the difference between two convex functions.To ensure the convexity of the functions f k for k = 2, . . . , M , a convexity regularization term γ C ∥Θ∥ 2 F is added to f k , denoted by fk , where γ C is chosen to be γ C ≥ max{0, -Λ min (∇ 2 f k (Θ))} such that Λ min (H fk ) ≥ 0 for k = 2, . . . , M . This regularization term guarantees that the minimum eigenvalue of the Hessian matrix H fk is non-negative, thereby ensuring the convexity of f k . Furthermore, the function g(Σ) is a closed, proper, and convex function. Proof. We detail the gradient and Hessian of functions f 1 and {f k+1 } M -1 k=1 as follows:EQUATION)Then given that f 1 and ∂f1 ∂Σ are Lipschitz continuous and bounded on the set {Σ ∈ M|∥Σ∥ 1 < ∞}, the function sequence {f k+1 } M -1 k=1 is also Lipschitz continuous.Proof of Theorem 7. Note that forEQUATIONthe problem (Fair CovGraph) is convex. Now, from Proposition 16 and Proposition 17, convexity and gradient Lipschitz continuity of objective functions {f k } M k=1 are verified. Then, the proof of Theorem 7 is a slightly modified version of the proof of Theorem 6.First, we present the convexity analysis and gradient Lipschitz continuity. Proposition 18 (Convexity of Fair BinNet). In the multi-objective optimization Problem (Fair BinNet), the functions f 1 and g are convex. Furthermore, by incorporating a convex regularization term γ I ∥Θ∥ 2 F for someγ I ≥ | min{ 1 2 Λ min (∇ 2 f k ), 0}| into each f k for k = 2, . . . , M , the set of functions {f k } Mk=2 are ensured to be convex as well.Proof. The function f 1 for Fair BinNet is defined as:EQUATION)To demonstrate the convexity of f 1 , observe that h(x) = log(1+exp(x)) is convex and nondecreasing. Since convexity is preserved under linear combination and summation, f 1 is convex by construction. Also, g(Σ) is a closed, proper, and convex function.Consider fk+1 (Θ) = f k+1 (Θ) + γ I ∥Θ∥ 2 F for k = 1, . . . , M -1, its Hessian matrix is given by:EQUATION)If γ I is chosen to be | min{ 1 2 Λ min (∇ 2 f k ), 0}| such that γ I I P 2 dominates any negative curvature in H f k+1 (Θ), then H fk+1 (Θ) will be positive semidefinite, leading the convexity of { fk+1 } M k=1 .Proposition 19 (Gradient Lipschitz Continuity of Fair BinNet). The gradients of f k for k = 1, . . . , M defined in the multi-objective optimization Problem (Fair BinNet) are Lipschitz continuous.Proof. For notational simplicity, we introduce the following substitutions: utilize L (Θ, X) in placeof f 1 (Θ), denote z θ = exp θ jj + j ′ ̸ =j θ jj ′ x ij ′ , z ϕ = exp ϕ jj + j ′ ̸ =j ϕ jj ′ x ij ′ .Then, we proceed to evaluate the gradient of the function f 1 in the context of Fair BinNet as follows:EQUATION)Given that h 1 (x) = exp(x)/ (1 + exp(x)) is Lipschitz continuous with Lipschitz constant 0.25, for any Θ, Φ ∈ M,EQUATION)It follows that there existsL 1 = N × P 2 ∈ R such that ∥∇f 1 (Θ) -∇f 1 (Φ)∥ F ≤ L 1 ∥Θ -Φ∥ F .Subsequently, the gradients of functions {f k+1 } M -1 k=1 in Fair BinNet are evaluated as:EQUATION)Then given that L and ∂L ∂Θ are Lipschitz continuous and bounded on the set {Θ ∈ M|∥Θ∥ 1 < ∞}, the function sequence {f k+1 } M -1 k=1 is also Lipschitz continuous.Proof of Theorem 8 for Fair BinNet. Building on Proposition 18 and Proposition 19, proof of Theorem 8 can be viewed as a nuanced adaptation of the proof presented in Theorem 6.The computational complexity of the fair GLasso and fair CovGraph algorithm depends on both the number of variables P and the number of observations N . We aim to demonstrate that our algorithm has a complexity of O max(N P 2 ,P 3 ) ϵ , which is similar to standard graph learning methods when K << N, P . This is applicable to our experimental results, where K = 2, 8, 2, 000 ≤ N ≤ 15, 000, and 5 ≤ P ≤ 120. The computational complexity is primarily influenced by the following factors:1. Number of Variables (P ): The complexity scales as O(P 3 ) due to matrix inversion for computing the gradient at each step. 2. Number of Observations (N ): Computing the empirical covariance matrix from the data has a complexity of O(N P 2 ). 3. Global Fair GMs Complexity: Considering factors 1 and 2, the complexity of each proximal gradient step applied to global fair GM is O(max(N P 2 , P 3 )). 4. Local GMs Complexity: Applying factors 1 and 2 to group-specific data, the complexity of each local GM is max(N k P 2 , P 3 ) for all k = 1, . . . , K. The total complexity of the local GMs is K k=1 max(N k P 2 , P 3 ). As established in Theorem 6 and Theorem 8 for fair inverse covariance and covariance estimation, the iteration complexity of our algorithm to achieve ϵ-accuracy is O 1 ϵ . Combining this result with the per-iteration complexity of the algorithm, the total time complexity of our optimization procedure is O max(N P 2 ,P 3 ) ϵ . Including the Local GMs computation, the total time complexity of fair GMsis O K k=1 max(N k P 2 , P 3 ) + max(N P 2 ,P 3 ) ϵ .Under the assumption that the number of groups is small (i.e., K << N , K << P , and K << 1/ϵ), the complexity reduces to O max(N P 2 ,P 3 ) ϵ . This complexity is of the same order as the complexity of running the proximal gradient method applied to covariance estimation and inverse covariance estimation. Therefore, for large N and P and a small number of groups, the time complexity of our algorithm is comparable to the standard method. In addition to theoretical analysis, we also provide sensitivity analysis experiments on P , N , K, and group imbalance in Appendix D.6-D.9.D.1 Iterative Soft-Thresholding Algorithm (ISTA) ISTA for sparse inverse covariance estimation is initially introduced by [60] and demonstrates a closed-form linear convergence rate. We adapt this approach and extend it to other GMs, utilizing it in the generation of both baseline and local graphs. Specifically, for a GM characterized by the loss function L (Θ; X) + λ∥Θ∥ 1 , we employ the following detailed algorithm:Algorithm 2 ISTA for GMs Input: Sample matrix X, initial iterate Θ (0) , maximum iteration T , step size ζ, regularization parameter λ, tolerance ϵ. Set t = 0. for t = 0, 1, . . . , T -1 do Gradient Step: Θ (t+1) ← Θ (t) -ζ∇L Θ (t) ; X Soft-Thresholding Step: Θ (t+1) ← η ζρ (Θ (t+1) ), (η ζρ (Θ)) jj ′ = sign(θ jj ′ ) max(|θ jj ′ | - ζρ, 0) if ∥∇L Θ (t+1) , X ∥ 1 ≤ ϵ then Break end if end for Output: Θ (t+1)As a supplement to Section 4.2, we detail the process of generating K block diagonal covariance matrices of dimensions P × P , denoted as {Σ k } K k=1 , each corresponding to distinct sensitive groups. The procedure is as follows:1. Firstly, we assume that each Σ k contains Q blocks and P is divisible by Q. For the first group, the covariance matrix Σ 1 is constructed as a block diagonal matrix:EQUATIONHere, each block B q is a sub-matrix filled with values drawn from a normal distribution N (0.7, 0.2).We specify the process of generating synthetic data for the simulation study in Section 4.3. This process adapts a hub node-based network as proposed by [72] , aiming to generate a sequence of networks {Θ} k l=1 . The process comprises the following steps: 1. Initialize a P × P matrix A, setting A jj ′ = 1 with a probability of 0.01 for all j < j ′ and A jj ′ = 0 otherwise. Ensure the matrix is symmetric by assigning A j ′ j = A jj ′ . From the set of nodes, randomly select H hub nodes and modify their corresponding rows and columns in A to 1 with a 99% probability or to 0 otherwise. 2. Construct another P ×P matrix E, where each element In the experiments of applying GLasso to the ADNI dataset, we investigate the influence of sensitive attributes on brain networks associated with Alzheimer's disease (AD) pathology. Specifically, we focus on the amyloid accumulation network using AV45 (florbetapir) positron emission tomography (PET) data [88] and the tau accumulation network using AV1451 (flortaucipir) PET data [46] . For the amyloid network, we consider the sensitive attribute of marital status, as previous studies suggest that marriage may affect the progression of dementia due to factors such as social support, cognitive stimulation, and lifestyle habits [21, 61] . The dataset is divided into two groups based on marital status: a single group with 52 samples and a married group with 1,018 samples, creating an imbalanced and high-dimensional setting that poses challenges for network estimation. In the tau accumulation network, we explore the impact of the sensitive attribute race, which separates the dataset into two groups: the white group with 755 samples and the non-white group with 118 samples. This division allows us to investigate potential disparities in tau pathology across racial groups. Throughout the experiments, the regularization parameter λ, which controls the sparsity of the estimated networks, is fixed at 0.3 for the AV45 data and 0.2 for the AV1451 data based on empirical observations. Besides, the dataset is normalized such that it has a mean of zero and a standard deviation of one.E jj ′ is i.i.d.. Set E jj ′ = 0 if A jj ′ = 0. Otherwise, draw E jj ′ fromResults. The numerical results in Table 4 demonstrate that Fair GLasso effectively reduces disparity error compared to standard GLasso, enhancing fairness while maintaining a good objective value. Figure 5 reveals notable differences in the learned network structures for AV45 results. The presence of edges between the left caudal middle frontal gyrus and right medial orbitofrontal cortex and between the left superior frontal gyrus and left superior parietal lobule in the GLasso graph suggests an increased influence of emotional factors on executive function and higher-order cognitive processes on amyloid accumulation, respectively [64, 6, 54] . Conversely, the absence of an edge between the left pars opercularis and left supramarginal gyrus in the Fair GLasso graph indicates a weaker association between language deficits and sensorimotor impairments in amyloid accumulation [18, 81, 38] .In contrast, the AV1451 results show primarily numerical differences between the two graphs, with edges remaining largely unchanged, suggesting that the tau accumulation network is robust to the sensitive attribute of race. These findings highlight the importance of considering fairness in brain imaging data analysis and the potential of Fair GLasso to uncover more equitable and unbiased patterns of amyloid and tau accumulation in Alzheimer's disease. Further research is needed to validate these findings in larger and more diverse cohorts and explore the biological mechanisms and clinical implications of the observed differences between standard GLasso and Fair GLasso graphs.Table 5 summarizes the details of credit data utilized on Fair CovGraph mentioned in Section 4.6. To improve the clarity of the visualizations, weak edges have been removed, and edges that show significant differences in values between the two methods are highlighted. It is important to note that even though some edges may appear unchanged in the visual comparison, their actual values will differ between the standard GLasso and Fair GLasso methods. 7 ) with respect to sample size N . %F 1 is slight, while %∆ changes are substantial, signifying fairness improvement without significant accuracy sacrifice. (b) Runtime (mean ± std) (results from Table 7 ) with respect to sample size N .us to systematically assess the effects of increasing N on the performance metrics and computational efficiency of the algorithm.The specific results are presented in Table 7 and visualized in Figure 7 . From these, it is evident that the sample size does not significantly impact the objective value, pairwise graph disparity error, or runtime. Our proposed method consistently maintains its effectiveness across different sample sizes. This stability highlights the robustness of our approach under varying data quantities.D.8 Sensitivity Analysis to Sample Size Ratio N 2 /N 1 Group Size (K) 6 ) with respect to group size K. %F 1 is slight, while %∆ changes are substantial, signifying fairness improvement without significant accuracy sacrifice. (b) Runtime (mean ± std) (results from Table 6 ) with respect to group size K. value and pairwise graph disparity error, performance maintains a good balance, with a significant enhancement in fairness.This conclusion aligns with our theoretical analysis of algorithmic complexity. Notably, as the group size increases, the pairwise graph disparity error also significantly rises, as shown in Table 9 . Consequently, our proposed method effectively enhances fairness, albeit at the cost of sacrificing a greater portion of the objective value. This trade-off is a critical aspect of our approach, balancing computational performance with the desired ethical outcomes in machine learning applications.To address the computational complexity of Fair GMs, we explore a range of optimization methods tailored to GLasso and multi-objective optimization (MOO):• GLasso Optimization Methods -Preconditioned Iterative Soft Thresholding Algorithm (PISTA): Efficiently handles large-scale sparse matrix operations [69] . -Graphical Iterative Shrinkage Thresholding Algorithm (GISTA): Employs an iterative framework for sparsity-inducing penalty functions in high-dimensional settings [60] . -Orthant-Based Newton Method (OBN): Uses second-order information for faster convergence in structured sparsity constraints [57] .• MOO Optimization Methods -Fast Iterative Shrinkage-Thresholding Algorithm (FISTA): Provides globally optimal convergence rates for MOO objectives [74] . -Stochastic Objective Selection Approach (SOSA): Introduces a randomized selection technique for optimizing multi-objective functions under varying conditions [70] .We validate these methods through comprehensive experiments on synthetic datasets. Our first evaluation uses data with 100 variables across two subgroups, each containing 1000 observations, generated following the procedure in Appendix D.2. This experiment demonstrates that faster optimization methods improve time complexity for both GLasso and MOO while maintaining performance. All GLasso methods achieve optimal loss, while Fair GLasso variants successfully reduce pairwise graph disparity error without significant performance degradation.To assess scalability, we extend our analysis to a larger dataset with 200 variables, maintaining the same experimental setup. Furthermore, we evaluate the efficiency of our approach with increased group complexity using synthetic data containing 100 variables across ten subgroups, each with 1000 observations. In this setting, SOSA reduces training time by approximately 36% compared to the original approach while preserving model fairness and robustness.The numerical results presented in Justification: Yes, our work is about fairness in unsupervised graphical models, we describe the detailed sensitive attribute and the resulting group splitting. The optimization method is also well described. Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material.Code is available at https://github.com/PennShenLab/Fair_GMsAvailable at http://www.cp.jku.at/datasets/LFM-1b/. 


Decision: The reviewers were generally positive on this paper, mostly because the authors did an excellent job in addressing their concerns during the response period. I hope that they will update their paper in such a way that the main body of the paper is improved similarly.

Reviews: review: paper_summary: paper_summary: The paper investigates the issue of bias  in 3 particular graphical models: Gaussian, Gaussian Covariance, and Binary Ising. In this regard, the authors propose a framework to enhance fairness in the estimation of graph models. They incorporate the difference of loss between the protected groups and the accuracy of graph models into a multi-objective optimization problem and solve it. The experiments show that the approach reduce the bias.
strength_weakness: strength_weakness: 1. The authors only consider one criterion of fairness, which is achieving equal loss among subgroups. How do the authors justify focusing solely on this criterion? Additionally, how might incorporating other fairness criteria impact the results?

2. Many recent works utilize reinforcement learning to discover fair structures. For instance, the paper "Rethinking Bias Mitigation: Fairer Architectures Make for Fairer Face Recognition" by S. Dooley et al. (NeurIPS 2023) explores such an approach. How do the authors justify their method compared to these more comprehensive approaches?

3. The authors employ multi-objective optimization, suggesting that the loss in GMs and pairwise graph disparity are at odds. However, many papers indicate that this is not always the case. For example, the paper "Is There a Trade-Off Between Fairness and Accuracy? A Perspective Using Mismatched Hypothesis Testing" by S. Dutta et al. (ICML 2020) discusses this. How do the authors justify their approach in light of such findings?

4. The learning algorithm appears to be very time-consuming. Do the authors have any comments on this aspect of their method?

5. Why do the authors not use more widely accepted metrics for fairness comparison and instead refer only to the difference in loss between subgroups?
questions: questions: Please refer to the weaknesses as well as these questions:

1. Practical graphs and practical fairness problems have high K and P values. This will raise the computational complexity to the extent that it may question the applicability of the method. Do the authors have comment on this?

2. What if the number of protected groups (like in the case of age) is high? How do the authors comment on these cases with regards to complexity analysis?
limitations: limitations: 1. The applicability for real problems
2. It does not scale well
3. The study is limited to Gaussian, Gaussian Covariance, and Binary Ising models.
paper_summary: paper_summary: This paper introduces a framework to address fairness in the estimation of graphical models, particularly focusing on Gaussian, Covariance, and Ising models. The motivation stems from the potential bias in standard GMs when handling data involving sensitive characteristics or protected groups. The proposed framework integrates a pairwise graph disparity error and a tailored loss function into a nonsmooth multi-objective optimization problem to reduce bias while maintaining the effectiveness of GMs.
strength_weakness: strength_weakness: Fairness metric: The choice and justification of fairness metrics used in the evaluation could be more thoroughly discussed. This would provide a clearer understanding of how fairness is quantified and the implications of these choices.

Unconvincing experiment results: While the authors claim that the proposed method "approach enhances fairness without compromising performance" in line 306, Table 1 indicates that 6 out of 7 datasets experience lower F1 scores, suggesting worse performance. 

Limited baseline: The paper could benefit from more extensive comparisons with existing fairness-aware methods in graphical models other than ISTA (the only baseline used in Table 2). This would help to contextualize the performance and advantages of the proposed framework.
questions: questions: 1) Could you elaborate more on the choice of fairness metrics, disparity error, used in your evaluation? What might be the implications for using such metric compared to other metrics, like those stemming from counterfactual evaluation?

2) While the authors claim that the proposed method improves fairness without compromising performance, Table 1 indicates that 6 out of 7 datasets experience lower F1 scores, suggesting worse performance. Could you clarify this discrepancy and provide a more detailed interpretation of Table 1?

3) It's mentioned that ISTA is used as baseline. What are other possible baselines? Why in particular choose ISTA as baseline?
limitations: limitations: The paper discuss the complexity limitation of the work but did not solve it.
paper_summary: paper_summary: This paper explores the issue of fairness in the estimation of graphical models (GMs), specifically Gaussian, Covariance, and Ising models. The authors introduce a comprehensive framework aimed at mitigating bias in GM estimation concerning protected attributes. This framework integrates pairwise graph disparity error and a customized loss function into a non-smooth multi-objective optimization problem. The pairwise graph disparity error assesses loss discrepancies across all groups, while the tailored loss function evaluates GM performance. Through this approach, the authors seek to achieve fairness across diverse sensitive groups while upholding GM performance, supported by theoretical proofs and experimental validation on various real-world datasets.
strength_weakness: strength_weakness: In the introduction and related works section, the paper does not provide complete information about GMs. Additionally, the paper does not clearly explain why it focuses on only three types of GMs instead of providing a general framework. The author needs to clearly highlight any difficulties encountered when applying this framework to each type of model.

In the second paragraph of Section 2 (Related Work on Fairness), the author does not discuss related works concerning fairness in GMs estimation. By focusing solely on fairness in unsupervised learning, the paper may limit readers' ability to gain a comprehensive understanding of the research field and the distinctiveness of the proposed method.

Moreover, the paper lacks clarity regarding the statement: “Our paper addresses the challenge of learning GMs without any predefined assumptions on the graph’s structure.” In the referenced paper, "Fair Community Detection and Structure Learning in Heterogeneous Graphical Models," the authors assume the presence of community structures and aim to learn a sparse undirected graph where demographic groups are fairly represented within these communities. This concept could be relevant to the current paper when the number of clusters is set to one. Therefore, the paper should clarify the differences and advantages of learning GMs without any predefined assumptions in the theoretical section.

In Section 4.8, trade-off analysis, the author mentions the shortcomings in runtime without clear discussion. Although there is a significant difference in runtime compared to the standard models, no reasons or future related works are mentioned in the subsequent sections.
questions: questions: Please explain the vagueness and confusion in the “Weaknesses” section.
limitations: limitations: The paper introduces a new framework that addresses fairness in the estimation of GMs, but it has some limitations. Firstly, the introduction lacks a comprehensive overview of current GMs and does not explain why the focus is on only three types of GMs rather than GMs in general. Secondly, the discussion of related work on fairness in GM estimation is insufficient and needs clearer explanations. Lastly, the implementation of different GMs requires distinct theorems with specific assumptions, which could limit the development and extension of the framework.
paper_summary: paper_summary: This paper proposes a novel method for estimating graphical models, particularly Gaussian, Covariance, and Ising models, from data, taking fairness into consideration. Fairness is defined on the graph disparity error, i.e., the difference between the loss of the model for a particular group and the minimum loss for that group. The goal is to make the graph disparity errors among differnt groups equal or close. The paper proposes to use non-smooth multi-objective optimization to solve the problem. The paper shows that the method achieves the weak Pareto optimality, and the convergence rates analyzed. Experiments show the efficacy of the method.
strength_weakness: strength_weakness: -	The optimality analysis requires the convexity of the loss function, and it remains unclear how the performance would be affected if the loss function were non-convex, which is common in machine learning. For example, if we use an encoding network to convert the data X into a representation space and compute the loss in the representation space, then the loss function would a non-convex.

-	There is just one baseline used in the experiments, which was published in 2012. It could strengthen the experiment section if modern GM estimation methods are used as baselines.
questions: questions: - Can the proposed method be applied to the estimation of directed graphs, like causal graphs?

- In Algorithm 1, the first step is to initialize local graph estimates for all groups. When the number of groups is large, the sample sizes of certain groups could become small, which may affect the accuracy of the local graph estimates. How will the accuracy of the local graph estimates affect the global graph estimate?
limitations: limitations: N/A

Comments: title: Response to Reviewer 1oDX
Response to Reviewer nXfy
Response to Reviewer tjnm
Response to Reviewer fakh
Follow up on Our Rebuttal
Rebuttal by Authors (Part 4)
Rebuttal by Authors (Part 3)
Rebuttal by Authors (Part 2)
Rebuttal by Authors (Part 2)
Rebuttal by Authors (Part 2)
Rebuttal by Authors (Part 2)
comment: Thank you for your review and valuable suggestions.
> The study is limited to Gaussian, Gaussian Covariance, and Binary Ising models.


> I have noticed that the reviewer "1oDX" also asked about why the authors focus only on Gaussian case. This needs to be discussed in the paper. 

We sincerely appreciate your feedback and the time you’ve taken to review our work. We would like to clarify that our work encompasses more than just the Gaussian case. Specifically, we address three canonical graphical models in statistical learning:

* Inverse Covariance Estimation,
* Covariance Estimation, 
and
* Ising Model Estimation.


It is crucial to emphasize that the Ising model is not Gaussian; it represents a distinct type of graphical model, particularly used for binary data.

Moreover, many other types of graphical models or statistical graph learning approaches can be viewed as variants derived from these canonical models.

We further elaborate on this below.

Firstly, we review two canonical types of graphical models (we have eliminated the discussion on Covariance Estimation as it falls under the category of Gaussian models for continuous data):

**Ising Model for Binary Data:**  This model is defined by a probability distribution over binary variables.  The distribution is determined by the so-called *Boltzmann distribution*:

$$p(\mathbf{x}; \Theta) = \left(Z(\Theta)\right)^{-1} \exp ( \sum_{j=1}^{P} \theta_{jj}x_j + \sum_{1 \leq j < j' \leq P} \theta_{jj'}x_jx_{j'} ). 
$$
 Here, $\Theta$ (graph matrix) is a symmetric matrix, and $Z(\Theta)$ is the partition function that normalizes the density. 

These models are primarily designed for binary data [[HT2009](https://www.jmlr.org/papers/volume10/hoefling09a/hoefling09a.pdf), [RWL2010](https://projecteuclid.org/journals/annals-of-statistics/volume-38/issue-3/----Custom-HTML----High/10.1214/09-AOS691.short)].


**Gaussian Graphical Model for Continuous Data:**  This model is suitable for continuous variables. The Gaussian graphical model is defined by the following distribution:
$$
p(\mathbf{x}; \Theta) = \frac{1}{(2\pi)^{P/2} |\Sigma|^{1/2}} \exp\left(-\frac{1}{2} \mathbf{x}^T \Sigma^{-1} \mathbf{x}\right)
$$
where $\mathbf{x} \sim \mathcal{N}(0, \Sigma)$, $\Sigma$ is the covariance matrix, and $\Theta=\Sigma^{-1}$ ( graph matrix) represents the precision matrix.

This model is widely used for continuous data [[MB2006](https://projecteuclid.org/journals/annals-of-statistics/volume-34/issue-3/High-dimensional-graphs-and-variable-selection-with-the-Lasso/10.1214/009053606000000281.short),
[YL2007](https://academic.oup.com/jrsssb/article-abstract/68/1/49/7110631),
[ABG2008](https://epubs.siam.org/doi/abs/10.1137/060670985),
[FHT2008](https://academic.oup.com/biostatistics/article-abstract/9/3/432/224260),
[RZY2008](https://arxiv.org/abs/0807.3734)].

Secondly, as we mentioned earlier, the Ising model is not designed for Gaussian cases but specifically for binary data rather than continuous data. This distinction makes the Ising model particularly suitable for applications involving binary data, such as modeling interactions between binary states in statistical physics or capturing binary decisions in social networks.

Finally, many other types of graphical models or statistical graph learning approaches can be viewed as variants derived from these two canonical models. For example,

* *Graphical models for ordinal variables* [[GLM2015](https://doi.org/10.1080/10618600.2014.889023)] assumes that the (categorical) data are generated by discretizing the marginal distributions of a latent multivariate Gaussian distribution.

* *Mixed Graphical Models* [[CLL2017](https://doi.org/10.1080/10618600.2016.1237362)] combines Ising and Gaussian models to handle mixed data, including both continuous and discrete variables.

* *Structured Graphical Models* [[KYC2020](https://www.jmlr.org/papers/volume21/19-276/19-276.pdf)] combines Gaussian graphical models with spectral graph theory for joint clustering and graph learning, as well as bipartite graph learning. 


Our current framework effectively addresses these canonical forms of graphical models—Gaussian for continuous data and Ising for binary data—within a single paper. 
Additionally, it can be extended to other models derived from these canonical forms, some of which are listed above. 

We believe that our non-smooth multi-objective framework is novel for graphical model analysis and addressing fairness issues. It provides tools for a rich class of graphical model estimation, and corrects the biases present in traditional approaches to graphical model estimation.
Thank you for the response. I raised my rating.
Thank you for the thorough review and positive feedback.
We appreciate your time and valuable suggestions.
Thank you for your thorough responses. I don't have further questions.
I thank the authors for providing additional results. I think the experiment sections looks stronger now and is willing to raise my score by 1.
I have noticed that the reviewer "1oDX" also asked about why the authors focus only on Gaussian case. This needs to be discussed in the paper.
I thank the authors for addressing my concerns. I have read the rebuttal.

Here are my post-rebuttal comments:
Weakness 1 and 5:
The loss disparity between K sub-groups seems rational to me. Just please bold in the paper that you need to have access to the label of the sensitive attribute.
Weakness2:
What I meant was to using reinforcement learning instead of multi-objective optimization. Please discuss it in the related work section.
weakness 3: 
Nothing new is added by the authors. My question remains.
weakness 4:
I agree with the authors. Using the strategies, the computational overhead can indeed be reduced.
Limitation 1:
Thanks for bringing it into my attention.
Limitation 3:
Nothing new is added and the question remains. The authors can have a brief discussion about this in the paper.

Given all the comments, I stick to my rating.
Dear Reviewers,


We greatly appreciate the time and effort you have invested in reviewing our work. 


We have carefully addressed all the concerns you raised, including the Fairness Metrics and Computational Complexity for large (P) and (N) highlighted by Reviewer nXfy, and the Limited Baselines identified by Reviewers tjnm and fakh. Additionally, we have provided detailed comparisons to related work on Fairness as suggested by Reviewer 1oDX and have responded to the concerns about the Non-convex Objective raised by Reviewer fakh.


As the discussion period draws to a close, we want to ensure that our responses have comprehensively addressed all of your questions and concerns. Please let us know if there are any aspects of our work that still require clarification or further discussion.


Best,

Authors
> **Limitation 1:** The applicability for real problems

**Response:** Thank you for your comment. We have provided the real-world applications of our proposed method in Section 4 ( Experiment). Specifically, in Sections 4.4 and 4.5, we analyzed the application of the fair GLasso on the Gene Regulatory Network and the Brain Amyloid/Tau Accumulation Network. Section 4.6 explored the application of Fair CovGraph in Credit Dataset, which promotes a more equitable credit system. In addition, we demonstrated the application of the Fair binary Ising model on Music Recommendation Systems in Section 4.7. These real-world examples illustrate the practical applicability of our proposed method.

> **Limitation 3:** The study is limited to Gaussian, Gaussian Covariance, and Binary Ising models.

**Response:** Thank you for your comment. Our paper focuses on three types of graphical models: Gaussian Graphical Models, Gaussian Covariance Graph Models, and Binary Ising Graphical Models. These models were chosen because they represent a diverse set of widely used GMs with distinct characteristics and applications. Indeed, each of these models is of significant interest to the research community and has been extensively studied. For example:

* GLasso: Widely used for sparse inverse covariance estimation, as demonstrated by Friedman et al. (2008) [[FHT2008](https://academic.oup.com/biostatistics/article/9/3/432/224260)] and Witten et al. (2011)[[WFS2011](https://www.jstor.org/stable/23248939)].

* Sparse Covariance Estimation: Important for large covariance matrix estimation, with notable contributions from Bickel and Levina (2008)[[BL2008](https://projecteuclid.org/journals/annals-of-statistics/volume-36/issue-1/Regularized-estimation-of-large-covariance-matrices/10.1214/009053607000000758.full)] and Cai et al. (2011)[[CLL2011](https://www.tandfonline.com/doi/abs/10.1198/jasa.2011.tm10155)].

* Binary Ising Model: Crucial for high-dimensional discrete data, highlighted by the work of Ravikumar et al. (2010)[[RWL2010](https://projecteuclid.org/journals/annals-of-statistics/volume-38/issue-3/High-dimensional-Ising-model-selection-using-ℓ1-regularized-logistic-regression/10.1214/09-AOS691.full)].

While our algorithm is a general framework applicable to various types of graphs, we focused on GLasso, sparse covariance estimation, and binary Ising models to provide detailed parameters and convergence guarantees specific to each model. This focus ensures a comprehensive analysis, including specific parameter tuning and convergence properties, thus providing a stronger and more practical contribution to the field. Generalizing too broadly might dilute these specific details and overlook important model-specific insights.

In summary, our focus on these three specific GMs balances the depth and breadth of our contributions, ensuring detailed and practically relevant results. We believe this graphical model estimation approach is beneficial for the community, offering clear and actionable insights for these widely studied graphical models.
> **Weakness 4:** The learning algorithm appears to be very time-consuming. Do the authors have any comments on this aspect of their method?

> **Question 1:** Practical graphs and practical fairness problems have high K and P values. This will raise the computational complexity to the extent that it may question the applicability of the method. Do the authors have comment on this? 

> **Question 2:** What if the number of protected groups (like in the case of age) is high? How do the authors comment on these cases with regards to complexity analysis?

> **Limitation 2:** It does not scale well

**Response:** We appreciate the reviewer’s observation regarding the time-consuming nature of our learning algorithm. Given the multi-objective optimization framework we employ, it is indeed expected to be more time-intensive. However, our primary goal is to provide an automatic (without any additional regularization or related hyperparameter tuning) and effective framework for fair graph learning. We acknowledge the importance of improving the computational efficiency of our method and outline potential future directions to address this concern.

Firstly, the time complexity partly arises from the local graph learning phase. This can be significantly accelerated using existing fast graphical model algorithms such as QUIC[[HSD2014](https://jmlr.org/papers/volume15/hsieh14a/hsieh14a.pdf)], SQUIC[[BS2016](http://www.icm.tu-bs.de/~bolle/Publicat/squic.pdf)], PISTA[[STY2023](https://arxiv.org/pdf/2205.10027)], GISTA[[GRR2012](https://arxiv.org/pdf/1211.2532)], OBN[[OON2012](https://papers.nips.cc/paper/2012/file/b3967a0e938dc2a6340e258630febd5a-Paper.pdf)], and ALM[[SG2010](https://proceedings.neurips.cc/paper/2010/file/2723d092b63885e0d7c260cc007e8b9d-Paper.pdf)]. Incorporating faster algorithms into the multi-objective optimization process can also enhance the overall efficiency of our learning algorithm.

Secondly, the increase in time complexity can also be attributed to the growing number of objectives in the multi-objective optimization solver. To mitigate this, we propose selecting a subset of objectives randomly in each iteration. This approach can reduce the computational load without significantly compromising the model's fairness and performance.

To validate our discussion, we conducted additional experiments using GLasso. In the first experiment, we generated synthetic data following the procedure detailed in the appendix of our paper, which includes two subgroups, each having 1000 observations and 100 variables. We applied various optimization algorithms to both the local graph learning and the multi-objective optimization processes. The detailed numerical results, presented in **Table 2 of the attached one-page PDF** file, demonstrate that all tested optimization algorithms for GLasso achieved optimal loss while maintaining performance improvements in fairness. Notably, GISTA and OBN for GLasso, along with FISTA for multi-objective optimization, significantly improved learning efficiency. We also repeated this experiment for high P values (P = 200), and the detailed numerical results presented in **Table 2 of the attached one-page PDF** further validate the improved efficiency of faster graphical model algorithms.

In the second experiment, we generated a synthetic dataset with ten subgroups to validate our approach to reducing time complexity. In each iteration of the multi-objective optimization, we randomly selected three objectives. The numerical results in **Table 2 of the attached one-page PDF file** indicate that this strategy effectively reduces runtime without substantially sacrificing model performance.

These results highlight the potential of our proposed methods to improve the computational efficiency of the learning algorithm while maintaining fairness and performance. We appreciate the reviewer’s feedback and believe these enhancements will further strengthen our framework.
> **Weakness 3:** The authors employ multi-objective optimization, suggesting that the loss in GMs and pairwise graph disparity are at odds. However, many papers indicate that this is not always the case. For example, the paper "Is There a Trade-Off Between Fairness and Accuracy? A Perspective Using Mismatched Hypothesis Testing" by S. Dutta et al. (ICML 2020) discusses this. How do the authors justify their approach in light of such findings?

**Response:** We appreciate the reviewer's comment on the trade-off between fairness and accuracy. The paper by S. Dutta et al. (ICML 2020) argues that the observed trade-off may stem from the real-world application of biased datasets rather than an intrinsic limitation of the algorithms themselves. They emphasize the importance of considering ideal distributions when evaluating fairness and accuracy and highlight the potential of active data collection to improve fairness without compromising accuracy.
In our work, we focus on the algorithm's fairness. Classical graphical models aim to achieve optimal loss for the given data distribution but may propagate and amplify existing biases. Our fairness-aware optimization framework aims to mitigate these biases systematically at the algorithmic level, independent of the data distribution. 

While classical graphical models achieve optimal loss, they do not inherently address fairness concerns. Our multi-objective optimization framework explicitly incorporates a fairness criterion, balancing the trade-off between minimizing loss and reducing pairwise graph disparity. Multi-objective optimization is a promising approach for integrating fairness considerations into machine learning models. Prior works, such as [[MBS2020](https://arxiv.org/pdf/2011.01821.pdf)], have formulated group fairness as a multi-objective optimization problem with separate objectives for each sensitive group’s risk. This strategy aligns closely with works like [[PAG2021](https://arxiv.org/pdf/2009.04441.pdf)] and [[MSG2021](https://arxiv.org/pdf/2110.01951.pdf)], which enhance fairness in classification. Further exploration by [[PBD2022](https://arxiv.org/pdf/2006.06137.pdf)], and [[KHF2018](https://arxiv.org/pdf/1911.04931.pdf)]  demonstrates the effectiveness of multi-objective optimization in improving fairness in unsupervised PCA algorithms. 

These efforts underscore the growing consensus on the potential of multi-objective optimization to address fairness issues in machine learning. On the other hand, existing solutions may not adequately address fairness in nonsmooth or high-dimensional settings, where our framework offers novel insights and solutions. To our knowledge, ours is the first nonsmooth multi-criteria method for fairness, applicable to other problems involving nonsmooth regularizations such as sparse fair PCA or fair supervised learning based on ERM with sparsity/norm-one regularization.
> **Question 2:** In Algorithm 1, the first step is to initialize local graph estimates for all groups. When the number of groups is large, the sample sizes of certain groups could become small, which may affect the accuracy of the local graph estimates. How will the accuracy of the local graph estimates affect the global graph estimate?

**Response:** Thank you for your insightful question. In Algorithm 1, we initialize local graph estimates for all groups, and we recognize that small sample sizes in certain groups could impact the accuracy of these local estimates. The accuracy of local graph estimates is indeed crucial, as these estimates form the foundation for the global graph estimate. When local estimates are based on small sample sizes, they may be less reliable, potentially introducing bias or variance that can propagate to the global graph estimate.

One possibility for handling the challenges of limited samples in certain subgroups is the bilevel optimization approach [[SXC2022](https://proceedings.neurips.cc/paper_files/paper/2022/file/dc96134e169de5aea1ba1fc34dfb8419-Paper-Conference.pdf)]. In this framework, the lower-level optimization learns subgroup-specific local graph matrices using the limited data available, guided by a fair model informed by the overall dataset. This approach is theoretically sound and empirically validated to improve fairness without sacrificing accuracy, even when some subgroups have limited samples. However, the  bilevel formulation for graph learning can be computationally expensive as it requires hyper-gradient computations. Further, this approach requires novel estimation and convergence analysis. Therefore, we will add this as a future direction.
To our knowledge, There are indeed a few works on fair graphical models that significantly differ from this work. The most relevant works to this study are [[68](https://arxiv.org/abs/2112.05128), [ZW2023](https://arxiv.org/abs/2311.13766), [NRB2024](https://arxiv.org/abs/2403.15591), [ZW2024](https://www.sciencedirect.com/science/article/pii/S0925231224009810)]. Specifically,  [[68](https://arxiv.org/abs/2112.05128)] initiated the learning of fair GMs using an $\ell_1$-regularized pseudo-likelihood method for joint GMs estimation and fair community detection. [[ZW2023](https://arxiv.org/abs/2311.13766), [ZW2024](https://www.sciencedirect.com/science/article/pii/S0925231224009810)] proposed a fair spectral clustering model that integrates graph construction, fair spectral embedding, and discretization into a single objective function. Unlike these models, which assume community structures, our study formulates fair GMs without such assumptions. Concurrently with this work, [[NRB2024](https://arxiv.org/abs/2403.15591)] proposed a regularization method for mitigating subpopulation bias for fair network topology inference. To our knowledge, this methodology significantly differs from ours, as we focus on developing three classes of fair GMs (Gaussian, Covariance, and Ising models) for imbalanced groups *without node attributes*, aiming to *automatically* ensure fairness through non-smooth multi-objective optimization.



> **Weakness 3:** Moreover, the paper  ... predefined assumptions in the theoretical section.


**Response:** Thank you for your comment.  

Firstly, the method in the referenced paper, ``Fair Community Detection and Structure Learning in Heterogeneous Graphical Models,'' is specifically designed for scenarios where there are inherent community or cluster structures in the data.  If the clustering were based on the number of observations (N) instead of the number of variables (P), setting the number of clusters to one would essentially imply that there is only one community. This would not address the issue of fairness as intended by the method. In such a scenario, the concept of fair community detection loses its relevance since there would be no differentiation between groups to balance fairness. Hence, the method is not applicable when the number of clusters is equal to one.

Secondly, the method assumes that the nodes are clusterable. However, in our case, we do not make such assumptions, as many applications involve graphs with nested structures rather than distinct communities. For example, consider a sparse graph with hub nodes [[TML2014](https://www.jmlr.org/papers/volume15/tan14b/tan14b.pdf)] or star graphs. In this setting, the nodes may not form distinct communities, challenging the method's assumptions.

Finally, we do not assume the existence of attributes on nodes, whereas the referenced method does. This can restrict applications such as gene network analysis (Section 5.1), where nodes could be genes that do not have attributes. In brain amyloid/tau accumulation network analysis (Section 5.2), nodes are regions of interest (ROIs) that do not have sensitive attributes. The requirement for nodes to have attributes limits the method's applicability in these contexts.

We will clarify these distinctions in the related work section. 

> **Weakness 4:** In Section 4.8, ... in the subsequent sections.  

**Response:** We appreciate the reviewer’s feedback regarding the trade-off analysis in Section 4.8. We acknowledge that the discussion on the shortcomings in runtime and potential future work was insufficiently detailed. We have expanded our discussion to include specific strategies to address these issues and outline potential future work.

One promising direction to accelerate the graphical model part of our algorithm is the adoption of faster optimization techniques. Specifically, the Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) has shown significant improvements in optimization efficiency in various contexts. By integrating FISTA into our framework, we can achieve considerable runtime reductions without compromising the accuracy and convergence properties of our models, as presented in **Table 2 in the attached one-page PDF file**.

Another factor impacting runtime is the number of groups considered. Reducing the number of groups or selecting representative sample groups can effectively manage computational complexity. One potential approach is the method proposed by [[SXC2022](https://openreview.net/pdf?id=YsRH6uVcx2l)], which suggests selecting a sample of groups randomly in each iteration of the optimization process to balance computational efficiency and model accuracy. By implementing a strategic sampling of groups, we can maintain the robustness and fairness of our models while significantly reducing the runtime. This approach will be explored in future work to optimize the balance between computational demands and model performance.
> **Weakness 3:** Limited baseline: The paper could benefit from more extensive comparisons with existing fairness-aware methods in graphical models other than ISTA (the only baseline used in Table 2). This would help to contextualize the performance and advantages of the proposed framework.

> **Question 3:** It's mentioned that ISTA is used as baseline. What are other possible baselines? Why in particular choose ISTA as baseline?

**Response:** Thank you for your valuable feedback. In response to your concern, we have conducted additional experiments to apply our multi-objective optimization framework to several state-of-the-art GLasso algorithms. We specifically focused on synthetic data to showcase the robustness and effectiveness of our framework beyond the ISTA algorithm.

We initially chose ISTA as a baseline due to its simplicity and widespread use in sparse inverse covariance estimation, making it an ideal reference point. However, understanding the necessity for a more comprehensive evaluation, we have now included comparisons with advanced baselines such as PISTA[[STY2023](https://arxiv.org/pdf/2205.10027)], GISTA[[GRR2012](https://arxiv.org/pdf/1211.2532)], and OBN[[OON2012](https://papers.nips.cc/paper/2012/file/b3967a0e938dc2a6340e258630febd5a-Paper.pdf)].

Our additional experiments, detailed in the revised version of our paper, demonstrate that our framework consistently outperforms these advanced baselines in enhancing fairness while maintaining competitive performance. Specifically, all new baseline methods reached the optimal loss for GLasso. When applying our multi-objective optimization framework, Fair GLasso achieved comparable results and significantly reduced pairwise graph disparity error. This comprehensive evaluation underscores the advantages of our approach in various scenarios, reinforcing its robustness and effectiveness.


> **Limitation:** The paper discuss the complexity limitation of the work but did not solve it.

**Response:** We appreciate the reviewer’s observation regarding the time-consuming nature of our learning algorithm. Given our use of a multi-objective optimization framework, it is expected that the process would be more time-intensive. Nonetheless, our primary goal remains to offer a simple yet effective framework for fair graph learning. We recognize the need to enhance the computational efficiency of our method and outline potential directions to achieve this.

Firstly, the time complexity is partly due to the local graph learning phase. This can be significantly accelerated by utilizing fast graphical model algorithms such as QUIC[[HSD2014](https://jmlr.org/papers/volume15/hsieh14a/hsieh14a.pdf)], SQUIC[[BS2016](http://www.icm.tu-bs.de/~bolle/Publicat/squic.pdf)], PISTA[[STY2023](https://arxiv.org/pdf/2205.10027)], GISTA[[GRR2012](https://arxiv.org/pdf/1211.2532)], OBN[[OON2012](https://papers.nips.cc/paper/2012/file/b3967a0e938dc2a6340e258630febd5a-Paper.pdf)], and ALM[[SG2010](https://proceedings.neurips.cc/paper/2010/file/2723d092b63885e0d7c260cc007e8b9d-Paper.pdf)]. Integrating these faster algorithms into the multi-objective optimization process can substantially improve the overall efficiency of our learning algorithm.

Secondly, the growing number of objectives in the multi-objective optimization solver can also increase time complexity. To address this, we propose randomly selecting a subset of objectives in each iteration. This approach can reduce the computational load while preserving the model's fairness and performance.

To support our discussion, we conducted additional experiments using GLasso. In the first experiment, we generated synthetic data, including two subgroups, as described in the appendix of our paper. We applied various optimization algorithms to both the local graph learning and multi-objective optimization processes. The detailed numerical results, presented in **Table 2 of the attached one-page PDF file**, show that all tested optimization algorithms for GLasso achieved optimal loss while maintaining fairness improvements. Notably, GISTA and OBN for GLasso, along with FISTA for multi-objective optimization, significantly enhanced learning efficiency.

In the second experiment, we generated a synthetic dataset with ten subgroups to validate our strategy for reducing time complexity. In each iteration of the multi-objective optimization, we randomly selected three objectives. The numerical results in **Table 2 of the attached one-page PDF file** indicate that this strategy effectively reduces runtime without substantially compromising model performance.

These findings underscore the potential of our proposed methods to improve the computational efficiency of the learning algorithm while maintaining fairness and performance. We appreciate the reviewer’s feedback and believe these enhancements will further strengthen our framework.