figurec

# Fair Structure Learning in Heterogeneous Graphical Models

Davoud Ataee Tarzanagh, Laura Balzano, and Alfred O. Hero
Department of Electrical Engineering and Computer Science, University of
Michigan, Ann Arbor, MI 48109, USA (E-mail<EMAIL_ADDRESS>hero@eecs.umich.edu).

###### Abstract

Inference of community structure in probabilistic graphical models may not be
consistent with fairness constraints when nodes have demographic attributes.
Certain demographics may be over-represented in some detected communities and
under-represented in others. This paper defines a novel $\ell_{1}$-regularized
pseudo-likelihood approach for fair graphical model selection. In particular,
we assume there is some community or clustering structure in the true
underlying graph, and we seek to learn a sparse undirected graph and its
communities from the data such that demographic groups are fairly represented
within the communities. Our optimization approach uses the demographic parity
definition of fairness, but the framework is easily extended to other
definitions of fairness. We establish statistical consistency of the proposed
method for both a Gaussian graphical model and an Ising model for,
respectively, continuous and binary data, proving that our method can recover
the graphs and their fair communities with high probability.

Key words: Graphical Models; Fairness; Community Detection; Generalized
Pseudo-Likelihood

##  1 Introduction

Probabilistic graphical models have been applied in a wide range of machine
learning problems to infer dependency relationships among random variables.
Examples include gene expression (Peng et al., 2009; Wang et al., 2009),
social interaction networks (Tan et al., 2014; Tarzanagh and Michailidis,
2018), computer vision (Hassner and Sklansky, 1981; Laferté et al., 2000;
Manning and Schutze, 1999), and recommender systems (Kouki et al., 2015; Wang
et al., 2015). Since in most applications the number of model parameters to be
estimated far exceeds the available sample size, it is necessary to impose
structure, such as sparsity or community structure, on the estimated
parameters to make the problem well-posed. With the increasing application of
structured graphical models and community detection algorithms in human-
centric contexts (Tan et al., 2013; Song et al., 2011; Glassman et al., 2014;
Burke et al., 2011; Pham et al., 2011; Das et al., 2014), there is a growing
concern that, if left unchecked, they can lead to discriminatory outcomes for
protected groups. For instance, the proportion of a minority group assigned to
some community can be far from its underlying proportion, even if detection
algorithms do not take the minority sensitive attribute into account in
decision making (Chierichetti et al., 2017). Such an outcome may, in turn,
lead to unfair treatment of minority groups. For example, in precision
medicine, patient-patient similarity networks over a biomarker feature space
can be used to cluster a cohort of patients and support treatment decisions on
particular clusters (Parimbelli et al., 2018; Lafit et al., 2019). If the
clusters learned by the algorithm are demographically imbalanced, this
treatment assignment may unfairly exclude under-represented groups from
effective treatments.

To the best of our knowledge, the estimation of fair structured graphical
models has not previously been addressed. However, there is a vast body of
literature on learning structured probabilistic graphical models. Typical
approaches to impose structure in graphical models, such as
$\ell_{1}$-regularization, encourage sparsity structure that is uniform
throughout the network and may therefore not be the most suitable choice for
many real world applications where data have clusters or communities, i.e.,
groups of graph nodes with similar connectivity patterns or stronger
connections within the group than to the rest of the network. Graphical models
with these properties are called heterogeneous.

It is known that if the goal is structured heterogeneous graph learning,
structure or community inference and graph weight estimation should be done
jointly. In fact, performing structure inference before weight estimation
results in a sub-optimal procedure (Marlin and Murphy, 2009). To overcome this
issue, some of the initial work focused on either inferring connectivity
information or performing graph estimation in case the connectivity or
community information is known a priori (Danaher et al., 2014; Guo et al.,
2011b; Gan et al., 2019; Ma and Michailidis, 2016; Lee and Liu, 2015). Recent
developments consider the two tasks jointly and estimate the structured
graphical models arising from heterogeneous observations (Kumar et al., 2020;
Hosseini and Lee, 2016; Hao et al., 2018; Tarzanagh and Michailidis, 2018;
Kumar et al., 2019; Gheche and Frossard, 2020; Pircalabelu and Claeskens,
2020; Cardoso et al., 2020; Eisenach et al., 2020).

In this paper, we develop a provably convergent penalized pseudo-likelihood
method to induce fairness into clustered probabilistic graphical models. More
specifically, we

  * •

Formulate a novel version of probabilistic graphical modeling that takes
fairness/bias into consideration. In particular, we assume there is some
community structure in our graph, and we seek to learn an undirected graph
from the data such that demographic groups are fairly represented within the
communities of the graph.

  * •

Provide a rigorous analysis of our algorithms showing that they can recover
fair communities with high probability. Furthermore, it is shown that the
estimators are asymptotically consistent in high dimensional settings for both
a Gaussian graphical model and an Ising model under standard regularity
assumptions.

  * •

Conclude by giving experimental results on synthetic and real-world datasets
where proportional clustering can be a desirable goal, comparing the
proportionality and objective value of standard graphical models to our
methods. Our experiments confirm that our algorithms tend to better estimate
graphs and their fairer communities compared to standard graphical models.

The remainder of the paper is organized as follows: Section 2 gives a general
framework for fair structure learning in graphs. Section 3 gives a detailed
statement of the proposed fair graphical models for continuous and binary
datasets. In Sections 4 and 5, we illustrate the proposed framework on a
number of synthetic and real data sets, respectively. Section 6 provides some
concluding remarks.

##### Notation.

For a set $\mathcal{S}$, $|\mathcal{S}|$ is the cardinality of $\mathcal{S}$,
and $\mathcal{S}^{c}$ is its complement. The reals and nonnegative reals
fields are denoted as ${\mathbb{R}}$ and ${\mathbb{R}}_{+}$, respectively. We
use lower-case and upper-case bold letters such as ${\bf{x}}$ and ${\bf{X}}$
to represent vectors and matrices, respectively, with $x_{i}$ and $x_{ij}$
denoting their elements. If all coordinates of a vector ${\bf{x}}$ are
nonnegative, we write ${\bf{x}}\geq 0$. The notation ${\bf{x}}>0$, as well as
${\bf{X}}\geq 0$ and ${\bf{X}}>0$ for matrices, are defined similarly. For a
symmetric matrix ${\bf{X}}\in{\mathbb{R}}^{n\times n}$, we write
${\bf{X}}\succ 0$ if ${\bf{X}}$ is positive definite, and ${\bf{X}}\succeq 0$
if it is positive semidefinite. ${\bf{I}}_{p}$, ${\bf{J}}_{p}$, and
${\bf{0}}_{p}$ denote the $p\times p$ identity matrix, matrix of all ones, and
matrix of all zeros, respectively. We use $\Lambda_{i}({\bf{X}})$,
$\Lambda_{\max}({\bf{X}})$ and $\Lambda_{\min}({\bf{X}})$ to denote the
$i$-th, maximum, and minimum singular values of ${\bf{X}}$, respectively. For
any matrix ${\bf{X}}$, we define $\|{\bf{X}}\|_{\infty}:=\max_{ij}|x_{ij}|$,
$\|{\bf{X}}\|_{1}:=\sum_{ij}|x_{ij}|$,
$\|{\bf{X}}\|:=\Lambda_{\max}({\bf{X}})$, and
$\|{\bf{X}}\|_{\text{F}}:=\sqrt{\sum_{ij}|x_{ij}|^{2}}$.

##  2 Fair Structure Learning in Graphical Models

We introduce a fair graph learning method that simultaneously accounts for
fair community detection and estimation of heterogeneous graphical models.

Let ${\bf{Y}}$ be an $n\times p$ matrix, with columns
${\bf{y}}_{1},\ldots,{\bf{y}}_{p}$. We associate to each column in ${\bf{Y}}$
a node in a graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$, where
$\mathcal{V}=\\{1,2,\ldots,p\\}$ is the vertex set and
$\mathcal{E}\in\mathcal{V}\times\mathcal{V}$ is the edge set. We consider a
simple undirected graph, without self-loops, and whose edge set contains only
distinct pairs. Graphs are conveniently represented by a $p\times p$ matrix,
denoted by $\Theta$, whose nonzero entries correspond to edges in the graph.
The precise definition of this usually depends on modeling assumptions,
properties of the desired graph, and application domain.

In order to obtain a sparse and interpretable graph estimate, many authors
have considered the problem

$\displaystyle~{}\begin{array}[]{ll}\underset{\mbox{\boldmath$\Theta$}}{\text{minimize}}&\begin{array}[]{c}\hskip
5.69046ptL(\mbox{\boldmath$\Theta$};{\bf{Y}})+\rho_{1}\|\mbox{\boldmath$\Theta$}\|_{1,\textnormal{off}}\end{array}\\\
\text{subj. to}&\begin{array}[t]{l}\hskip
5.69046pt\mbox{\boldmath$\Theta$}\in\mathcal{M}.\end{array}\end{array}$ (5)

Here, $L$ is a loss function;
$\rho_{1}\|\mbox{\boldmath$\Theta$}\|_{1,\textnormal{off}}$ is the
$\ell_{1}$-norm regularization applied to off-diagonal elements of $\Theta$
with parameter $\rho_{1}>0$; and $\mathcal{M}$ is a convex constraint subset
of ${\mathbb{R}}^{p\times p}$. For instance, in the case of a Gaussian
graphical model, we could take
$L(\mbox{\boldmath$\Theta$};{\bf{Y}})=-\log\det(\mbox{\boldmath$\Theta$})+\operatorname{trace}({\bf{S}}\mbox{\boldmath$\Theta$})$,
where ${\bf{S}}=n^{-1}\sum_{i=1}^{n}{\bf{y}}_{i}{\bf{y}}_{i}^{\top}$ and
$\mathcal{M}$ is the set of $p\times p$ positive definite matrices. The
solution to (5) can then be interpreted as a sparse estimate of the inverse
covariance matrix (Banerjee et al., 2008; Friedman et al., 2008). Throughout,
we assume that $L(\mbox{\boldmath$\Theta$};{\bf{Y}})$ and $\mathcal{M}$ are
convex function and set, respectively.

###  2.1 Model Framework

We build our fair graph learning framework using (5) as a starting point. Let
$\mathcal{V}$ denote the set of nodes. Suppose there exist $K$ disjoint
communities of nodes denoted by
$\mathcal{V}=\mathcal{C}_{1}\cup\cdots\cup\mathcal{C}_{K}$ where
$\mathcal{C}_{k}$ is the subset of nodes from $\mathcal{G}$ that belong to the
$k$–th community. For each candidate partition of $n$ nodes into $K$
communities, we associate it with a _partition matrix_
${\bf{Q}}\in[0,1]^{p\times p}$, such that $q_{ij}=1/|C_{k}|$ if and only if
nodes $i$ and $j$ are assigned to the $k$th community. Let $\mathcal{Q}_{pK}$
be the set of all such partition matrices, and $\bar{{\bf{Q}}}$ the true
partition matrix associated with the ground-truth clusters
$\\{\bar{\mathcal{C}}_{k}\\}_{k=1}^{K}$.

Assume the set of nodes contains $H$ demographic groups such that
$\mathcal{V}=\mathcal{D}_{1}\cup\cdots\cup\mathcal{D}_{H}$, potentially with
overlap between the groups. Chierichetti et al. (2017) proposed a model for
fair clustering requiring the representation in each cluster to preserve the
global fraction of each demographic group $\mathcal{D}_{h}$, i.e.,

$\frac{|\mathcal{D}_{h}\cap\mathcal{C}_{k}|}{|\mathcal{C}_{k}|}=\frac{|\mathcal{D}_{h}|}{p}~{}~{}\textnormal{for
all}~{}~{}k\in[K].$ (6)

Let ${\bf{R}}\in\\{0,1\\}^{p\times p}$ be such that $r_{ij}=1$ if and only if
nodes $i$ and $j$ are assigned to the same demographic group, with the
convention that $r_{ii}=1,\forall i$. One will notice that (6) is equivalent
to ${\bf{R}}({\bf{I}}-{\bf{1}}{\bf{1}}^{\top}/p){\bf{Q}}=0$. Let
${\bf{A}}_{1}:={\bf{R}}({\bf{I}}-{\bf{1}}{\bf{1}}^{\top}/p)$ and
${\bf{B}}_{1}:=\text{diag}(\mbox{\boldmath$\epsilon$}){\bf{J}}_{p}$ for some
$\mbox{\boldmath$\epsilon$}>0$ that controls how close we are to exact
demographic parity. Under this setting, we introduce a general optimization
framework for fair structured graph learning via a trace regularization and a
fairness constraint on the partition matrix ${\bf{Q}}$ as follows:

$\displaystyle\begin{array}[]{ll}\underset{\mbox{\boldmath$\Theta$},~{}\mbox{\boldmath$Q$}}{\text{minimize}}&\begin{array}[]{c}\hskip
14.22636ptL(\mbox{\boldmath$\Theta$};{\bf{Y}})+\rho_{1}\|\mbox{\boldmath$\Theta$}\|_{1,\textnormal{off}}+\rho_{2}\operatorname{trace}\big{(}({\bf{S}}+{\bf{Q}})G(\mbox{\boldmath$\Theta$})\big{)}\end{array}\\\
\text{subj. to}&\begin{array}[t]{l}\hskip
5.69046pt\mbox{\boldmath$\Theta$}\in\mathcal{M},~{}~{}{\bf{A}}_{1}{\bf{Q}}\leq{\bf{B}}_{1},~{}~{}\text{and}~{}~{}{\bf{Q}}\in\cup_{K}\mathcal{Q}_{pK}.\end{array}\end{array}$
(11)

Here, $G(\mbox{\boldmath$\Theta$}):\mathcal{M}\rightarrow\mathcal{M}$ is a
function of $\Theta$ (introduced in Sections 3.1 and 3.2).

We clarify the purpose of each component of the minimization (11). The term
$\rho_{1}\|\mbox{\boldmath$\Theta$}\|_{1,\textnormal{off}}$ shrinks small
entries of $\Theta$ to $0$, thus enforcing sparsity in $\Theta$ and
consequently in $\mathcal{G}$. This term controls the presence of edges
between any two nodes irrespective of the community they belong to, with
higher values for $\rho_{1}$ forcing sparser estimators. The polyhedral
constraint is the fairness constraint, enforcing that every community contains
the $\epsilon$-approximate proportion of elements from each demographic group
$\mathcal{D}_{h}$, $h\in[H]$, matching the overall proportion. The term
$\rho_{2}\operatorname{trace}\left(({\bf{S}}+{\bf{Q}})G(\mbox{\boldmath$\Theta$})\right)$
enforces community structure in a similarity graph,
$G(\mbox{\boldmath$\Theta$})$. A similar linear trace term, i.e.,
$\operatorname{trace}({\bf{Q}}\mbox{\boldmath$\Theta$})$ is used as an
objective function in (Cai et al., 2015; Amini et al., 2018; Hosseini and Lee,
2016; Pircalabelu and Claeskens, 2020; Eisenach et al., 2020) when estimating
communities of networks. However, the perturbation of the membership matrix
with either a sample covariance for which the population inverse covariance
satisfies Assumption (A2) or some positive definite matrix is necessary for
developing a consistent fair graphical model.

###  2.2 Relaxation

The problem (11) is in general NP-hard due to its constraint on ${\bf{Q}}$.
However, it can be relaxed to a computationally feasible problem. To do so, we
exploit algebraic properties of a community matrix ${\bf{Q}}$. By definition,
we see that ${\bf{Q}}$ must have the form
${\bf{Q}}=\mbox{\boldmath$\Psi$}\mbox{\boldmath$\Gamma$}\mbox{\boldmath$\Psi$}^{\intercal}$,
where ${\bf{\Gamma}}$ is block diagonal with size ${p_{k}}\times{p_{k}}$
blocks on the diagonal with all entries equal to $1$ associated with $k$-th
community, $\Psi$ is some permutation matrix, and the number of communities
$K$ is unknown. The set of all matrices ${{\bf{Q}}}$ of this form is non-
convex. The key observation is that any such ${\bf{Q}}$ satisfies several
convex constraints such as (i) all entries of ${\bf{Q}}$ are nonnegative, (ii)
all diagonal entries of ${\bf{Q}}$ are $1$, and (iii) ${\bf{Q}}$ is positive
semi-definite (Cai et al., 2015; Amini et al., 2018; Li et al., 2021). Without
loss of generality, we assume that the permutation matrix corresponding to the
ground truth communities is the identity, i.e.,
$\mbox{\boldmath$\Psi$}={\bf{I}}$. Now, let

${\bf{A}}:=[{\bf{A}}_{1};{\bf{I}}_{p}]~{}~{}\text{and}~{}~{}{\bf{B}}:=[{\bf{B}}_{1};{\bf{J}}_{p}].$

Thus, we propose the following relaxation:

$\displaystyle~{}\begin{array}[]{ll}\underset{\mbox{\boldmath$\Theta$},~{}{\bf{\mbox{\boldmath$Q$}}}}{\text{minimize}}&\begin{array}[]{c}\hskip
14.22636ptL(\mbox{\boldmath$\Theta$};{\bf{Y}})+\rho_{1}\|\mbox{\boldmath$\Theta$}\|_{1,\textnormal{off}}+\rho_{2}\operatorname{trace}\big{(}({\bf{S}}+{\bf{Q}})G(\mbox{\boldmath$\Theta$})\big{)}\end{array}\\\
\text{subj. to}&\begin{array}[t]{l}\hskip
14.22636pt\mbox{\boldmath$\Theta$}\in\mathcal{M},~{}~{}\textnormal{and}~{}~{}{\bf{Q}}\in\mathcal{N},\end{array}\end{array}$
(12e) where $\displaystyle\mathcal{N}$
$\displaystyle=\left\\{{\bf{Q}}\in{\mathbb{R}}^{p\times
p}:~{}{\bf{Q}}\succeq{\bf{0}},~{}{\bf{0}}\leq{\bf{A}}{\bf{Q}}\leq{\bf{B}},~{}q_{ii}=1~{}\textnormal{for}~{}1\leq
i\leq p\right\\}.$ (12f)

The solution of (12) jointly learns the fair community matrix ${\bf{Q}}$ and
the network estimate $\Theta$. We highlight the following attractive
properties of the formulation (12): (i) the communities are allowed to have
significantly different sizes; (ii) the number of communities $K$ may grow as
$p$ increases; (iii) the knowledge of $K$ is not required for fair community
detection, and (iv) the objective function (12e) is convex in $\Theta$ given
${\bf{Q}}$ and conversely.

###  2.3 Algorithm

In order to solve (12), we use an _alternating direction method of
multipliers_ (ADMM) algorithm (Boyd et al., 2011). ADMM is an attractive
algorithm for this problem, as it allows us to decouple some of the terms in
(12) that are difficult to optimize jointly. In order to develop an ADMM
algorithm for (12) with guaranteed convergence, we reformulate it as follows:

$\displaystyle\operatorname*{minimize}_{\mbox{\boldmath$Q$},\mbox{\boldmath$\Theta$},\mbox{\boldmath$\Omega$}}~{}~{}~{}~{}$
$\displaystyle
L(\mbox{\boldmath$\Theta$};{\bf{Y}})+\rho_{2}\operatorname{trace}\big{(}{\bf{S}}G(\mbox{\boldmath$\Theta$})\big{)}+\rho_{1}\|\mbox{\boldmath$\Omega$}\|_{1,\textnormal{off}}+\rho_{2}\operatorname{trace}\big{(}{\bf{Q}}G(\mbox{\boldmath$\Omega$})\big{)}$
subj. to
$\displaystyle\mbox{\boldmath$\Theta$}=\mbox{\boldmath$\Omega$},~{}~{}\mbox{\boldmath$\Theta$}\in\mathcal{M},~{}~{}\textnormal{and}~{}~{}{\bf{Q}}\in\mathcal{N}.$
(13)

The scaled augmented Lagrangian function for (2.3) takes the form

$\begin{split}\Upsilon_{\gamma}(\mbox{\boldmath$\Theta$},\mbox{\boldmath$\Omega$},{\bf{Q}},{\bf{W}})&:=L(\mbox{\boldmath$\Theta$};{\bf{Y}})+\rho_{2}\operatorname{trace}\big{(}{\bf{S}}G(\mbox{\boldmath$\Theta$})\big{)}+\rho_{1}\|\mbox{\boldmath$\Omega$}\|_{1,\textnormal{off}}\\\
&+\rho_{2}\operatorname{trace}\big{(}{\bf{Q}}G(\mbox{\boldmath$\Omega$})\big{)}+\frac{\gamma}{2}\|\mbox{\boldmath$\Theta$}-\mbox{\boldmath$\Omega$}+{\bf{W}}\|_{\textnormal{F}}^{2},\\\
\end{split}$ (14)

where $\mbox{\boldmath$\Theta$}\in\mathcal{M}$, $\Omega$, and
${\bf{Q}}\in\mathcal{N}$ are the primal variables; ${\bf W}$ is the dual
variable; and $\gamma>0$ is a dual parameter. We note that the scaled
augmented Lagrangian can be derived from the usual Lagrangian by adding a
quadratic term and completing the square (Boyd et al., 2011, Section 3.1.1).

The proposed ADMM algorithm requires the following updates:

$\displaystyle{\bf{Q}}^{(t+1)}$ $\displaystyle\leftarrow$
$\displaystyle\underset{{\bf{Q}}\in\mathcal{N}}{\text{argmin}}~{}\Upsilon_{\gamma}\left({\bf{Q}},\mbox{\boldmath$\Omega$}^{(t)},\mbox{\boldmath$\Theta$}^{(t)},{\bf{W}}^{(t)}\right),$
(15a) $\displaystyle\mbox{\boldmath$\Omega$}^{(t+1)}$
$\displaystyle\leftarrow$
$\displaystyle\underset{\mbox{\boldmath$\Omega$}}{\text{argmin}}~{}\Upsilon_{\gamma}\left({\bf{Q}}^{(t+1)},\mbox{\boldmath$\Omega$},\mbox{\boldmath$\Theta$}^{(t)},{\bf{W}}^{(t)}\right),$
(15b) $\displaystyle\mbox{\boldmath$\Theta$}^{(t+1)}$
$\displaystyle\leftarrow$
$\displaystyle\underset{\mbox{\boldmath$\Theta$}\in\mathcal{M}}{\text{argmin
}}\Upsilon_{\gamma}\left({\bf{Q}}^{(t+1)},\mbox{\boldmath$\Omega$}^{(t+1)},\mbox{\boldmath$\Theta$},{\bf{W}}^{(t)}\right),$
(15c) $\displaystyle{\bf{W}}^{(t+1)}$ $\displaystyle\leftarrow$
$\displaystyle{\bf{W}}^{(t)}+\mbox{\boldmath$\Theta$}^{(t+1)}-\mbox{\boldmath$\Omega$}^{(t+1)}.$
(15d)

A general algorithm for solving (12) is provided in Algorithm 1. Note that
update for $\Theta$, ${\bf{Q}}$, and $\Omega$ depends on the form of the
functions $L$ and $G$, and is addressed in Sections 3.1 and 3.2. We also note
that ${\bf{Q}}$ sub-problem in S1. can be solved via a variety of convex
optimization methods such as CVX (Grant and Boyd, 2014) and ADMM (Cai et al.,
2015; Amini et al., 2018). In the following sections, we consider special
cases of (12) that lead to estimation of Gaussian graphical model and an Ising
model for, respectively, continuous and binary data.

Algorithm 1 Fair Graph Learning via Alternating Direction Method of
Multipliers

Initialize primal variables $\mbox{\boldmath$\Theta$}^{(0)},{\bf{Q}}^{(0)}$,
and $\mbox{\boldmath$\Omega$}^{(0)}$; dual variable ${\bf{W}}^{(0)}$; and
positive constants $\gamma,\nu$.
Iterate until the stopping criterion
$\max\left\\{\frac{\|{\mathbf{\Theta}}^{(t+1)}-{\mathbf{\Theta}}^{(t)}\|_{F}^{2}}{\|{\mathbf{\Theta}}^{(t)}\|_{\textnormal{F}}^{2}},\frac{\|{\mathbf{Q}}^{(t+1)}-{\mathbf{Q}}^{(t)}\|_{\textnormal{F}}^{2}}{\|{\mathbf{Q}}^{(t)}\|_{\textnormal{F}}^{2}}\right\\}\leq\nu$
is met, where ${\bf\Theta}^{(t+1)}$ and ${\bf Q}^{(t+1)}$ are the value of
$\bf\Theta$ and $\bf Q$, respectively, obtained at the $t$-th iteration:

  1. S1.

${\bf{Q}}^{(t+1)}\leftarrow\underset{{\bf{Q}}\in\mathcal{N}}{\text{argmin
}}~{}\operatorname{trace}\big{(}{\bf{Q}}G(\mbox{\boldmath$\Omega$}^{(t)})\big{)}$.

  2. S2.

$\mbox{\boldmath$\Omega$}^{(t+1)}\leftarrow\underset{\mbox{\boldmath$\Omega$}}{\text{argmin
}}\rho_{2}\operatorname{trace}\big{(}{\bf{Q}}^{(t+1)}G(\mbox{\boldmath$\Omega$})\big{)}+\rho_{1}\|\mbox{\boldmath$\Omega$}\|_{1}+\frac{\gamma}{2}\big{\|}\mbox{\boldmath$\Theta$}^{(t)}-\mbox{\boldmath$\Omega$}+{\bf{W}}^{(t)}\big{\|}_{\textnormal{F}}^{2}.$

  3. S3.

$\mbox{\boldmath$\Theta$}^{(t+1)}\leftarrow\underset{\mbox{\boldmath$\Theta$}\in\mathcal{M}}{\arg\min}~{}L(\mbox{\boldmath$\Theta$};{\bf{Y}})+\rho_{2}\operatorname{trace}\big{(}{\bf{S}}G(\mbox{\boldmath$\Theta$})\big{)}+\frac{\gamma}{2}\|\mbox{\boldmath$\Theta$}-\mbox{\boldmath$\Omega$}^{(t+1)}+{\bf{W}}^{(t)}\|_{\textnormal{F}}^{2}$.

  4. S4.

${\bf{W}}^{(t+1)}\leftarrow{\bf{W}}^{(t)}+\mbox{\boldmath$\Theta$}^{(t+1)}-\mbox{\boldmath$\Omega$}^{(t+1)}$.

We have the following global convergence result for Algorithm 1.

######  Theorem 1.

Algorithm 1 converges globally for any sufficiently large $\gamma$111The lower
bound is given in (Wang et al., 2019, Lemma 9)., i.e., starting from any
$(\mbox{\boldmath$\Theta$}^{(0)},\mbox{\boldmath$\Omega$}^{(0)},{\bf{Q}}^{(0)},{\bf{W}}^{(0)})$,
it generates a sequence
$(\mbox{\boldmath$\Theta$}^{(t)},\mbox{\boldmath$\Omega$}^{(t)},{\bf{Q}}^{(t)},{\bf{W}}^{(t)})$
that converges to a stationary point of (14).

In Algorithm 1, S1. dominates the computational complexity in each iteration
of ADMM (Cai et al., 2015; Amini et al., 2018). In fact, an exact
implementation of this subproblem of optimization requires a full SVD, whose
computational complexity is $O(p^{3})$. When $p$ is as large as hundreds of
thousands, the full SVD is computationally impractical. An open question is
how to facilitate the implementation, or whether there exists a surrogate that
is computationally inexpensive. A possible remedy is to apply an iterative
approximation method in each iteration of ADMM, the full SVD is replaced by a
partial SVD where only the leading eigenvalues and eigenvectors are computed.
Although this type of method may get stuck in local minimizers, given the fact
that SDP implementation can be viewed as a preprocessing before K-means
clustering, such a low-rank iterative method might be helpful. It is worth
mentioning that when the number of communities $K$ is known, the computational
complexity of ADMM is much smaller than $O(p^{3})$; see Remark 4 for further
discussion.

###  2.4 Related Work

To the best of our knowledge, the fair graphical model proposed here is the
first model that can jointly learn fair communities simultaneously with the
structure of a conditional dependence network. Related work falls into two
categories: graphical model estimation and fairness.

##### Estimation of graphical models.

There is a substantial body of literature on methods for estimating network
structures from high-dimensional data, motivated by important biomedical and
social science applications (Liljeros et al., 2001; Robins et al., 2007; Guo
et al., 2011c; Danaher et al., 2014; Friedman et al., 2008; Tan et al., 2014;
Guo et al., 2015; Tarzanagh and Michailidis, 2018). Since in most applications
the number of model parameters to be estimated far exceeds the available
sample size, the assumption of sparsity is made and imposed through
regularization of the learned graph. An $\ell_{1}$ penalty on the parameters
encoding the network edges is the most common choice (Friedman et al., 2008;
Meinshausen et al., 2006; Karoui, 2008; Cai and Liu, 2011; Xue et al., 2012;
Khare et al., 2015; Peng et al., 2009). This approach encourages sparse
uniform network structures that may not be the most suitable choice for real
world applications, that may not be uniformly sparse. As argued in (Danaher et
al., 2014; Guo et al., 2011b; Tarzanagh and Michailidis, 2018) many networks
exhibit different structures at different scales. An example includes a
densely connected community in the social networks literature. Such structures
in social interaction networks may correspond to groups of people sharing
common interests or being co-located (Tarzanagh and Michailidis, 2018), while
in biological systems to groups of proteins responsible for regulating or
synthesizing chemical products and in precision medicine the communities may
be patients with common disease suceptibilities. An important part of the
literature therefore deals with the estimation of hidden communities of nodes,
by which it is meant that certain nodes are linked more often to other nodes
which are similar, rather than to dissimilar nodes. This way the nodes from
communities that are more homogeneous within the community than between
communities where there is a larger degree of heterogeneity. Some of the
initial work focused on either inferring connectivity information (Marlin and
Murphy, 2009) or performing graph estimation in case the connectivity or
community information is known a priori (Danaher et al., 2014; Guo et al.,
2011b; Gan et al., 2019; Ma and Michailidis, 2016; Lee and Liu, 2015), but not
both tasks simultaneously. Recent developments consider the two tasks jointly
and estimate the structured graphical models arising from heterogeneous
observations (Kumar et al., 2020; Hosseini and Lee, 2016; Hao et al., 2018;
Tarzanagh and Michailidis, 2018; Kumar et al., 2019; Gheche and Frossard,
2020; Pircalabelu and Claeskens, 2020; Cardoso et al., 2020; Eisenach et al.,
2020).

##### Fairness.

There is a growing body of work on fairness in machine learning. Much of the
research is on fair supervised methods; see, Chouldechova and Roth (2018);
Barocas et al. (2019); Donini et al. (2018) and references therein. Our paper
adds to the literature on fair methods for unsupervised learning tasks
(Chierichetti et al., 2017; Celis et al., 2017; Samadi et al., 2018;
Tantipongpipat et al., 2019; Oneto and Chiappa, 2020; Caton and Haas, 2020;
Kleindessner et al., 2019). We discuss the work on fairness most closely
related to our paper. Chierichetti et al. (2017) proposed the notion of
fairness for clustering underlying our paper: namely, that each cluster has
proportional representation from different demographic groups (Feldman et al.,
2015; Zafar et al., 2017). Chierichetti et al. (2017) provides approximation
algorithms that incorporate this fairness notion into $K$-center as well as
$K$-median clustering. Kleindessner et al. (2019) extend this to K-means and
provide a provable fair spectral clustering method; they implement $K$-means
on the subspace spanned by the smallest fair eigenvectors of Laplacian matrix.
Unlike these works, which assume that the graph structure and/or the number of
communities is given in advance, an appealing feature of our method is to
learn fair community structure while estimating heterogeneous graphical
models.

##  3 The Fair Graphical Models

In the following subsections, we consider two special cases of (12) that lead
to estimation of graphical models for continuous and binary data.

###  3.1 Fair Pseudo-Likelihood Graphical Model

Suppose ${\bf{y}}_{i}=(y^{1}_{i},\ldots,y^{p}_{i})$ are i.i.d. observations
from $N(\mbox{\boldmath$0$},\mbox{\boldmath$\Sigma$})$, for $i=1,\ldots,n$.
Denote the sample of the $i$th variable as
${\bf{y}}^{i}=(y^{i}_{1},\ldots,y^{i}_{n})$. Let
$\omega_{ij}=-\theta_{ij}/\theta_{ii}$, for all $j\neq i$. We note that the
set of nonzero coefficients of $\omega_{ij}$ is the same as the set of nonzero
entries in the row vector of $\theta_{ij}~{}(i\neq j)$, which defines the set
of neighbors of node $\theta_{ij}$. Using an $\ell_{1}$-penalized regression,
Meinshausen et al. (2006) estimates the zeros in $\Theta$ by fitting separate
Lasso regressions for each variable ${\bf{y}}^{i}$ given the other variables
as follows

$F_{i}(\mbox{\boldmath$\Theta$};{\bf{Y}})=\|{\bf y}^{i}-\sum_{j\neq
i}\omega_{ij}{\bf y}^{j}\|^{2}+\rho_{1}\sum_{1\leq i<j\leq
p}\left|\omega_{ij}\right|,~{}~{}\textnormal{where}~{}~{}\omega_{ij}=-\theta_{ij}/\theta_{ii}.$
(16)

These individual Lasso fits give neighborhoods that link each variable to
others. Peng et al. (2009) improve this neighborhood selection method by
taking the natural symmetry in the problem into account (i.e.,
$\theta_{ij}=\theta_{ij}$ ), and propose the following joint objective
function (called SPACE):

$\displaystyle F(\mbox{\boldmath$\Theta$};{\bf{Y}})$
$\displaystyle=\frac{1}{2}\sum_{i=1}^{p}\big{(}-n\log\theta_{ii}+w_{i}\|{\bf
y}^{i}-\sum_{j\neq
i}\dot{\omega}_{ij}\sqrt{\frac{\theta_{jj}}{\theta_{ii}}}{\bf
y}^{j}\|^{2}\big{)}+\rho\sum_{1\leq i<j\leq p}\left|\dot{\omega}_{ij}\right|$
$\displaystyle=\frac{1}{2}\sum_{i=1}^{p}\big{(}-n\log\theta_{ii}+w_{i}\|{\bf
y}^{i}+\sum_{j\neq i}\frac{\theta_{ij}}{\theta_{ii}}{\bf
y}^{j}\|^{2}\big{)}+\rho_{1}\sum_{1\leq i<j\leq
p}\left|\dot{\omega}_{ij}\right|,$ (17)

where $\\{w_{i}\\}_{i=1}^{p}$ are nonnegative weights and
$\dot{\omega}_{ij}=-\frac{\theta_{ij}}{\sqrt{\theta_{ii}\theta_{jj}}}$ denotes
the partial correlation between the $i$th and $j$th variables for $1\leq i\neq
j\leq p$. Note that $\dot{\omega}_{ij}=\dot{\omega}_{ji}$ for $i\neq j$.

It is shown in (Khare et al., 2015) that the above expression is not convex.
Setting $w_{i}=\theta_{ii}^{2}$ and putting the $\ell_{1}$-penalty term on the
partial covariances $\theta_{ij}$ instead of on the partial correlations
$\dot{\omega}_{ij}$, they obtain a convex pseudo-likelihood approach with good
model selection properties called CONCORD. Their objective takes the form

$\displaystyle F(\mbox{\boldmath$\Theta$};{\bf{Y}})$
$\displaystyle=\sum_{i=1}^{p}\big{(}-n\log\theta_{ii}+\frac{1}{2}\|\theta_{ii}{\bf
y}^{i}+\sum_{j\neq i}\theta_{ij}{\bf
y}^{j}\|_{2}^{2}\big{)}+\rho_{1}\sum_{1\leq i<j\leq p}|\theta_{ij}|.$ (18)

Note that the penalized matrix version of the CONCORD objective can be
obtained by setting
$L(\mbox{\boldmath$\Theta$};{\bf{Y}})=n/2[-\log|\text{diag}(\mbox{\boldmath$\Theta$})^{2}|+\text{trace}({\bf{S}}\mbox{\boldmath$\Theta$}^{2})]$
in (5).

Our proposed fair graphical model formulation (called FCONCORD) is a fair
version of CONCORD from (18). In particular, letting
$G(\mbox{\boldmath$\Theta$})=\mbox{\boldmath$\Theta$}^{2}$ and

$\mathcal{M}=\left\\{\mbox{\boldmath$\Theta$}\in{\mathbb{R}}^{p\times
p}:~{}\theta_{ij}=\theta_{ji},\text{ and }\theta_{ii}>0,\mbox{ for every
}1\leq i,j\leq p\right\\},$

in (12), our problem takes the form

$\displaystyle\operatorname*{minimize}_{\mbox{\boldmath$\Theta$},~{}\mbox{\boldmath$Q$}}~{}~{}~{}~{}F(\mbox{\boldmath$\Theta$},{\bf{Q}};{\bf{Y}})$
$\displaystyle:=\frac{n}{2}\big{[}-\log|\text{diag}(\mbox{\boldmath$\Theta$})^{2}|$
$\displaystyle+\operatorname{trace}\big{(}((1+\rho_{2}){\bf{S}}+\rho_{2}{\bf{Q}})\mbox{\boldmath$\Theta$}^{2}\big{)}\big{]}+\rho_{1}\|\mbox{\boldmath$\Theta$}\|_{1,\textnormal{off}}$
subj. to
$\displaystyle\mbox{\boldmath$\Theta$}\in\mathcal{M}~{}~{}\text{and}~{}~{}{\bf{Q}}\in\mathcal{N}.$
(19)

Here, $\mathcal{M}$ and $\mathcal{N}$ are the graph adjacency and fairness
constraints, respectively.

######  Remark 2.

When $\rho_{2}=0$, i.e., without a fairness constant and the second trace
term, the objective in (3.1) reduces to the objective of the CONCORD
estimator, and is similar to those of SPACE (Peng et al., 2009), SPLICE (Rocha
et al., 2008), and SYMLASSO (Friedman et al., 2010). Our framework is a
generalization of these methods to fair graph learning and community
detection, when the demographic group representation holds.

Problem (3.1) can be solved using Algorithm 1. The update for $\Omega$ and
$\Theta$ in S2. and S3. can be derived by minimizing

$\displaystyle\Upsilon_{1,\gamma}(\mbox{\boldmath$\Omega$})$
$\displaystyle:=\frac{n\rho_{2}}{2}\operatorname{trace}\big{(}{\bf{Q}}\mbox{\boldmath$\Omega$}^{2}\big{)}+\rho_{1}\|\mbox{\boldmath$\Omega$}\|_{1,\textnormal{off}}+\frac{\gamma}{2}\left\|\mbox{\boldmath$\Theta$}-\mbox{\boldmath$\Omega$}+{\bf{W}}\right\|_{\textnormal{F}}^{2},$
(20a) $\displaystyle\Upsilon_{2,\gamma}(\mbox{\boldmath$\Theta$})$
$\displaystyle:=\frac{n}{2}\big{[}-\log|\text{diag}(\mbox{\boldmath$\Theta$})^{2}|+\operatorname{trace}\big{(}(1+\rho_{2}){\bf{S}}\mbox{\boldmath$\Theta$}^{2}\big{)}\big{]}+\frac{\gamma}{2}\left\|\mbox{\boldmath$\Theta$}-\mbox{\boldmath$\Omega$}+{\bf{W}}\right\|_{\textnormal{F}}^{2},$
(20b)

with respect to $\Omega$ and $\Theta$, respectively.

For $1\leq i\leq j\leq p$, define the matrix function
${\bf{T}}_{ij}:\mathcal{M}\rightarrow\mathcal{M}$ by

$\displaystyle{\bf{T}}_{ij}(\mbox{\boldmath$\Omega$})$
$\displaystyle\leftarrow\operatorname*{argmin}_{{\tilde{\mbox{\boldmath$\Omega$}}}}\big{\\{}\Upsilon_{1,\gamma}({\tilde{\mbox{\boldmath$\Omega$}}}):~{}{\tilde{\omega}}_{kl}=\omega_{kl}\;\forall(k,l)\neq(i,j)\big{\\}},$
(21a) $\displaystyle{\bf{T}}_{ij}(\mbox{\boldmath$\Theta$})$
$\displaystyle\leftarrow\operatorname*{argmin}_{{\tilde{\mbox{\boldmath$\Theta$}}}}\big{\\{}\Upsilon_{2,\gamma}({\tilde{\mbox{\boldmath$\Theta$}}}):~{}{\tilde{\theta}}_{kl}=\theta_{kl}\;\forall(k,l)\neq(i,j)\big{\\}}.$
(21b)

For each $(i,j)$, ${\bf{T}}_{ij}(\mbox{\boldmath$\Omega$})$ and
${\bf{T}}_{ij}(\mbox{\boldmath$\Theta$})$ updates the $(i,j)$-th entry with
the minimizer of (20a) and (20b) with respect to $\omega_{ij}$ and
$\theta_{ij}$, respectively, holding all other variables constant. Given
${\bf{T}}_{ij}(\mbox{\boldmath$\Omega$})$ and
${\bf{T}}_{ij}(\mbox{\boldmath$\Theta$})$, the update for $\Omega$ and
$\Theta$ in S2. and S3. can be obtained by a similar coordinate-wise descent
algorithm proposed in Peng et al. (2009); Khare et al. (2015). Closed form
updates for ${\bf{T}}_{ij}(\mbox{\boldmath$\Omega$})$ and
${\bf{T}}_{ij}(\mbox{\boldmath$\Theta$})$ are provided in Lemma 3.

######  Lemma 3.

Let $\gamma_{n}:=\gamma n$. For $1\leq i\leq p$, define

$\displaystyle a_{i}:=(1+\rho_{2})s_{ii}+\gamma_{n}\qquad
b_{i}:=(1+\rho_{2})\sum_{j\neq
i}\theta_{ij}s_{ij}+\gamma_{n}(\omega_{ii}-w_{ii}),$ $\displaystyle
c_{i}:=\rho_{2}q_{ii}+\gamma_{n},\qquad\qquad d_{i}:=\rho_{2}\sum_{j\neq
i}q_{ij}\omega_{ij}+\gamma_{n}(w_{ii}+\theta_{ii})$

Further, for $1\leq i<j\leq p$, let

$\displaystyle a_{ij}:=(1+\rho_{2})(s_{ii}+s_{jj})+\gamma_{n}\quad
b_{ij}:=(1+\rho_{2})\big{(}\sum_{j^{\prime}\neq
j}\theta_{ij^{\prime}}s_{jj^{\prime}}+\sum_{i^{\prime}\neq
i}\theta_{i^{\prime}j}s_{ii^{\prime}}\big{)}+\gamma_{n}(w_{ij}-\omega_{ij})$
$\displaystyle c_{ij}:=\rho_{2}(q_{ii}+q_{jj})+\gamma_{n},\qquad
d_{ij}:=\rho_{2}\big{(}\sum_{j^{\prime}\neq
j}\omega_{ij^{\prime}}q_{jj^{\prime}}+\sum_{i^{\prime}\neq
i}\omega_{i^{\prime}j}q_{ii^{\prime}}\big{)}+\gamma_{n}(w_{ij}+\theta_{ij})$

Then, we have

$\displaystyle\big{(}{\bf{T}}_{ii}(\mbox{\boldmath$\Omega$})\big{)}_{ii}=-\frac{d_{i}}{c_{i}},\qquad\qquad\qquad\big{(}{\bf{T}}_{ii}(\mbox{\boldmath$\Theta$})\big{)}_{ii}=\frac{-b_{i}+\sqrt{b_{i}^{2}+4a_{i}}}{2a_{i}},~{}~{}\textnormal{for
all}~{}~{}1\leq i\leq p,$ (22a)
$\displaystyle\big{(}{\bf{T}}_{ij}(\mbox{\boldmath$\Omega$})\big{)}_{ij}=S\big{(}-\frac{d_{ij}}{c_{ij}},\frac{\rho_{1}}{\gamma_{n}}\big{)},\quad\big{(}{\bf{T}}_{ij}(\mbox{\boldmath$\Theta$})\big{)}_{ij}=-\frac{b_{ij}}{a_{ij}},~{}~{}\textnormal{for
all}~{}~{}1\leq i<j\leq p,$ (22b)

where $S(\alpha,\beta):=\textnormal{sign}(\alpha)\max(|\alpha|-\beta,0)$.

######  Remark 4.

In the case when $K$ and $H$ are known, the complexity of Algorithm 1 is of
the same order as the CONCORD estimator, SPACE, SPLICE, and SYMLASSO. In fact,
computing the fair clustering matrix ${\bf{Q}}$ requires $O((p-H+1)^{2}K)$
operations. On the other hand, it follows from (Khare et al., 2015, Lemma 5)
that $\Theta$ and $\Omega$ updates can be performed with complexity
$\min\left(O(np^{2}),O(p^{3})\right)$. This shows that when the number of
communities is known, the computational cost of each iteration of FCONCORD is
$\max\left(\min\left(O(np^{2}),O(p^{3})\right),(p-H+1)^{2}K\right)$.

####  3.1.1 Large Sample Properties of FCONCORD

We show that under suitable conditions, the FCONCORD estimator achieves both
model selection consistency and estimation consistency.

As in other studies (Khare et al., 2015; Peng et al., 2009), for the
convergence analysis we assume that the diagonal of the graph matrix $\Theta$
and partition matrix ${\bf{Q}}$ are known. Let
$\mbox{\boldmath$\theta$}^{o}=(\theta_{ij})_{1\leq i<j\leq p}$ and
${\bf{q}}^{o}=(q_{ij})_{1\leq i<j\leq p}$ denote the vector of off-diagonal
entries of $\Theta$ and ${\bf{Q}}$, respectively. Let
$\mbox{\boldmath$\theta$}^{d}$ and ${\bf{q}}^{d}$ denote the vector of
diagonal entries of $\Theta$ and ${\bf{Q}}$, respectively. Let
$\bar{\mbox{\boldmath$\theta$}}^{o},\bar{\mbox{\boldmath$\theta$}}^{d},\bar{{\bf{q}}}^{o}$,
and $\bar{{\bf{q}}}^{d}$ denote the true value of
${\mbox{\boldmath$\theta$}}^{o},{\mbox{\boldmath$\theta$}}^{d},{\bf{q}}^{o}$,
and ${\bf{q}}^{d}$, respectively. Let $\mathcal{B}$ denote the set of non-zero
entries in the vector $\bar{\mbox{\boldmath$\theta$}}^{o}$ and define

$\displaystyle q$
$\displaystyle:=|\mathcal{B}|,~{}~{}~{}~{}\Psi(p,H,K):=(p-H+1)((p-H+1)/K-1).$
(23)

In our consistency analysis, we let the regularization parameters
$\rho_{1}=\rho_{1n}$ and $\rho_{2}=\rho_{2n}$ vary with $n$.

The following standard assumptions are required:
Assumption A

  1. (A1)

The random vectors ${\bf{y}}_{1},\dots,{\bf{y}}_{n}$ are _i.i.d._ sub-Gaussian
for every $n\geq 1$, i.e., there exists $M>0$ such that
$\|{\bf{u}}^{\top}{\bf{y}}_{i}\|_{\psi_{2}}\leq
M\sqrt{{\mathbb{E}}({\bf{u}}^{\top}{\bf{y}}_{i})^{2}}$,
$\forall{\bf{u}}\in{\mathbb{R}}^{p}$. Here,
$\|{\bf{y}}\|_{\psi_{2}}=\sup_{t\geq
1}({\mathbb{E}}|{\bf{y}}|^{t})^{\frac{1}{t}}/\sqrt{t}$.

  2. (A2)

There exist constants $\tau_{1},\tau_{2}\in(0,\infty)$ such that

$\tau_{1}<\Lambda_{\min}(\bar{\mbox{\boldmath$\Theta$}})\leq\Lambda_{\max}(\bar{\mbox{\boldmath$\Theta$}})<\tau_{2}.$

  3. (A3)

There exists a constant $\tau_{3}\in(0,\infty)$ such that

$0\leq\Lambda_{\min}(\bar{{\bf{Q}}})\leq\Lambda_{\max}(\bar{{\bf{Q}}})<\tau_{3}.$

  4. (A4)

For any $K,H\in[p]$, we have $K\leq p-H+1$.

  5. (A5)

There exists a constant $\delta\in(0,1]$ such that, for any
$(i,j)\in\mathcal{B}^{c}$

$\displaystyle\left|\bar{{\bf{H}}}_{ij,\mathcal{B}}\bar{{\bf{H}}}_{\mathcal{B},\mathcal{B}}^{-1}\operatorname{sign}(\bar{\mbox{\boldmath$\theta$}}^{o}_{\mathcal{B}})\right|\leq(1-\delta),$

where for $1\leq i,j,t,s\leq p$ satisfying $i<j$ and $t<s$,

$\displaystyle\bar{{\bf{H}}}_{ij,ts}:={\mathbb{E}}_{\bar{\mbox{\boldmath$\theta$}}^{o}}\Big{(}\frac{\partial^{2}L(\bar{\mbox{\boldmath$\theta$}}^{d},,{\mbox{\boldmath$\theta$}}^{o};{\bf{Y}})}{\partial\theta_{ij}\partial\theta_{ts}}\Big{|}_{\mbox{\boldmath$\theta$}^{o}=\bar{\mbox{\boldmath$\theta$}}^{o}}\Big{)}.$
(24)

Assumptions (A2)–(A3) guarantee that the eigenvalues of the true graph matrix
$\bar{\mbox{\boldmath$\Theta$}}$ and those of the true membership matrix
$\bar{{\bf{Q}}}$ are well-behaving. Assumption (A4) links how $H,K$ and $p$
can grow with $n$. Note that $K$ is limited in order for fairness constraints
to be meaningful; if $K>p-H+1$ then there can be no community with $H$ nodes
among which we enforce fairness. Assumption (A5) corresponds to the
incoherence condition in Meinshausen et al. (2006), which plays an important
role in proving model selection consistency of $\ell_{1}$ penalization
problems. Zhao and Yu (2006) show that such a condition is almost necessary
and sufficient for model selection consistency in lasso regression, and they
provide some examples when this condition is satisfied. Note that Assumptions
(A1), (A2), and (A5) are identical to Assumptions (C0)–(C2) in Peng et al.
(2009). Further, it follows from Peng et al. (2009) that under Assumption (A5)
for any $(i,j)\in\mathcal{B}^{c}$,

$\|\bar{{\bf{H}}}_{ij,\mathcal{B}}\bar{{\bf{H}}}_{\mathcal{B},\mathcal{B}}^{-1}\|\leq
M(\bar{\mbox{\boldmath$\theta$}}^{o})$ (25)

for some finite $M(\bar{\mbox{\boldmath$\theta$}}^{o})$.

Next, inspired by Peng et al. (2009); Khare et al. (2015), we prove estimation
consistency for the nodewise FCONCORD.

######  Theorem 5.

Suppose Assumptions (A1)–(A5) are satisfied. Assume further that
$p=O(n^{\kappa})$ for some $\kappa>0$, $\rho_{1n}=O(\sqrt{\log p/n})$,
$n>O(q\log(p))$ as $n\rightarrow\infty$, $\rho_{2n}=O(\sqrt{\log(p-H+1)/n})$,
$\rho_{2n}\leq\delta\rho_{1n}/((1+M(\bar{\mbox{\boldmath$\theta$}}^{o}))\tau_{2}\tau_{3})$,
and $\mbox{\boldmath$\epsilon$}=\mbox{\boldmath$0$}$. Then, there exist finite
constants $C(\bar{\mbox{\boldmath$\theta$}}^{o})$ and $D(\bar{{\bf{q}}}^{o})$,
such that for any $\eta>0$, the following events hold with probability at
least $1-O(\exp(-\eta\log p))$:

  * •

There exists a minimizer
$({\widehat{\mbox{\boldmath$\theta$}}}^{o},{\widehat{{\bf{q}}}}^{o})$ of (3.1)
such that

$\max\left(\|\widehat{\mbox{\boldmath$\theta$}}^{o}-\bar{\mbox{\boldmath$\theta$}}^{o}\|,\|{\widehat{{\bf{q}}}}^{o}-\bar{{\bf{q}}}^{0}\|\right)\\\
\leq\max\left(C(\bar{\mbox{\boldmath$\theta$}}^{o})\rho_{1n}\sqrt{q}/(1+\rho_{2n}),D(\bar{{\bf{q}}}^{o})\rho_{2n}\sqrt{\Psi(p,H,K)}\right),$

where $q$ and $\Psi(p,H,K)$ are defined in (23).

  * •

If $\min_{(i,j)\in\mathcal{B}}\bar{\theta}_{ij}\geq
2C(\bar{\mbox{\boldmath$\theta$}}^{o})\rho_{1n}\sqrt{q}/(1+\rho_{2n})$, then
${\widehat{\mbox{\boldmath$\theta$}}}^{o}_{\mathcal{B}^{c}}=0$.

Theorem 6 provides sufficient conditions on the quadruple $(n,p,H,K)$ and the
model parameters for the FCONCORD to succeed in consistently estimating the
neighborhood of every node in the graph and communities simultaneously. Notice
that if $H=1$ (no fairness) and $K=p$ (no clustering) we recover the results
of Khare et al. (2015); Peng et al. (2009).

###  3.2 Fair Ising Graphical Model

In the previous section, we studied the fair estimation of graphical models
for continuous data. Next, we focus on estimating an Ising Markov random field
(Ising, 1925), suitable for binary or categorical data. Let
${\bf{y}}=(y_{1},\ldots,y_{p})\in\\{0,1\\}^{p}$ denote a binary random vector.
The Ising model specifies the probability mass function

$p({\bf{y}})=\frac{1}{\mathcal{W}(\mbox{\boldmath$\Theta$})}\exp\Big{(}\sum_{j=1}^{p}\theta_{jj}y_{j}+\sum_{1\leq
j<j^{\prime}\leq p}\theta_{jj^{\prime}}y_{j}y_{j^{\prime}}\Big{)}.$ (26)

Here, $\mathcal{W}(\mbox{\boldmath$\Theta$})$ is the partition function, which
ensures that the density function in (26) integrates to one; $\Theta$ is a
$p\times p$ symmetric matrix that specifies the graph structure:
$\theta_{jj^{\prime}}=0$ implies that the $j$th and $j^{\prime}$th variables
are conditionally independent given the remaining ones.

Several sparse estimation procedures for this model have been proposed. Lee et
al. (2007) considered maximizing an $\ell_{1}$-penalized log-likelihood for
this model. Due to the difficulty in computing the log-likelihood with the
expensive partition function, alternative approaches have been considered. For
instance, Ravikumar et al. (2010) proposed a neighborhood selection approach
which involves solving $p$ logistic regressions separately (one for each node
in the network), which leads to an estimated parameter matrix that is in
general not symmetric. In contrast, others have considered maximizing an
$\ell_{1}$-penalized pseudo-likelihood with a symmetric constraint on $\Theta$
(Guo et al., 2011c, a; Tan et al., 2014; Tarzanagh and Michailidis, 2018).
Under the probability model above, the negative $\log$-pseudo-likelihood for
$n$ observations takes the form

$L(\mbox{\boldmath$\Theta$};{\bf{Y}})=-\sum_{j=1}^{p}\sum_{j^{\prime}=1}^{p}\theta_{jj^{\prime}}s_{jj^{\prime}}+\frac{1}{n}\sum_{i=1}^{n}\sum_{j=1}^{p}\log\big{(}1+\exp(\theta_{jj}+\sum_{j^{\prime}\neq
j}\theta_{jj^{\prime}}y_{ij^{\prime}})\big{)}.$ (27)

We propose to additionally impose the fairness constraints on $\Theta$ in (27)
in order to obtain a sparse binary network with fair communities. This leads
to the criterion

$\displaystyle\operatorname*{minimize}_{\mbox{\boldmath$\Theta$},\mbox{\boldmath$Q$}}~{}~{}~{}~{}~{}~{}$
$\displaystyle
F(\mbox{\boldmath$\Theta$},{\bf{Q}};{\bf{Y}})+\rho_{1}\|\mbox{\boldmath$\Theta$}\|_{1,\textnormal{off}}:=-n\sum_{j=1}^{p}\sum_{j^{\prime}=1}^{p}\theta_{jj^{\prime}}(s_{jj^{\prime}}+\rho_{2}q_{jj^{\prime}})$
$\displaystyle+\sum_{i=1}^{n}\sum_{j=1}^{p}\log\big{(}1+\exp(\theta_{jj}+\sum_{j^{\prime}\neq
j}\theta_{jj^{\prime}}y_{ij^{\prime}})\big{)}+\rho_{1}\sum_{1\leq i<j\leq
p}|\theta_{ij}|,$ subj. to
$\displaystyle\mbox{\boldmath$\Theta$}\in\mathcal{M}~{}~{}\text{and}~{}~{}{\bf{Q}}\in\mathcal{N}.$
(28)

Here, $\mathcal{M}$ and $\mathcal{N}$ are the graph and fairness constraints,
respectively.

We refer to the solution to (3.2) as the _Fair Binary Network_ (FBN). An
interesting connection can be drawn between our technique and a fair variant
of Ising block model discussed in Berthet et al. (2016), which is a
perturbation of the mean field approximation of the Ising model known as the
Curie-Weiss model: the sites are partitioned into two blocks of equal size and
the interaction between those within the same block is stronger than across
blocks, to account for more order within each block. One can easily seen that
the Ising block model is a special case of (3.2).

An ADMM algorithm for solving (3.2) is given in Algorithm 1. The update for
$\Omega$ in S2. can be obtained from (20a) by replacing
$\mbox{\boldmath$\Omega$}^{2}$ with $\Omega$. We solve the update for $\Theta$
in S3. using a relaxed variant of Barzilai-Borwein method (Barzilai and
Borwein, 1988). The details are given in (Tarzanagh and Michailidis, 2018,
Algorithm 2).

####  3.2.1 Large Sample Properties of FBN

In this section, we present the model selection consistency property for the
separate regularized logistic regression. The spirit of the proof is similar
to Ravikumar et al. (2010), but since their model does not include membership
matrix ${\bf{Q}}$ and fairness constraints that are significant differences.

Similar to Section 3.1.1, let
$\mbox{\boldmath$\theta$}^{o}=(\theta_{ij})_{1\leq i<j\leq p}$ and
${\bf{q}}^{o}=(q_{ij})_{1\leq i<j\leq p}$ denote the vector of off-diagonal
entries of $\Theta$ and ${\bf{Q}}$, respectively. Let
$\mbox{\boldmath$\theta$}^{d}$ and ${\bf{q}}^{d}$ denote the vector of
diagonal entries of $\Theta$ and ${\bf{Q}}$, respectively. Let
$\bar{\mbox{\boldmath$\theta$}}^{o},\bar{\mbox{\boldmath$\theta$}}^{d},\bar{{\bf{q}}}^{o},\bar{{\bf{q}}}^{d}$,
$\bar{\mbox{\boldmath$\Theta$}}$ and $\bar{{\bf{Q}}}$ denote the true value of
${\mbox{\boldmath$\theta$}}^{o},{\mbox{\boldmath$\theta$}}^{d},{\bf{q}}^{o},{\bf{q}}^{d}$,
$\Theta$ and ${\bf{Q}}$, respectively. Let $\mathcal{B}$ denote the set of
non-zero entries in the vector $\bar{\theta}^{o}$, and let $q=|\mathcal{B}|$.

Denote the log-likelihood for the $i$-th observation by

$L_{i}(\mbox{\boldmath$\theta$}^{d},{\mbox{\boldmath$\theta$}}^{o};{\bf{Y}})=-\sum_{j=1}^{p}y_{ij}\big{(}\theta_{jj}+\sum_{j\neq
j^{\prime}}\theta_{jj^{\prime}}y_{ij^{\prime}}\big{)}+\log\big{(}1+\exp(\mbox{\boldmath$\theta$}_{jj}\sum_{j\neq
j^{\prime}}\theta_{jj^{\prime}}y_{ij^{\prime}})\big{)}.$ (29)

The population Fisher information matrix of $L$ at
$(\bar{\mbox{\boldmath$\theta$}}^{d},\bar{\mbox{\boldmath$\theta$}}^{o})$ can
be expressed as
$\bar{{\bf{H}}}={\mathbb{E}}(\nabla^{2}L_{i}(\bar{\mbox{\boldmath$\theta$}}^{d},\bar{\mbox{\boldmath$\theta$}}^{o};{\bf{Y}}))$
and its sample counterpart
$\bar{{\bf{H}}}^{n}=1/n\sum_{i=1}^{n}\nabla^{2}L_{i}(\bar{\mbox{\boldmath$\theta$}}^{d},\bar{\mbox{\boldmath$\theta$}}^{o};{\bf{Y}})$.
Let

$\displaystyle
v_{ij}=\dot{v}_{ij}(1-\dot{v}_{ij}),~{}~{}\textnormal{where}~{}~{}\dot{v}_{ij}=\frac{\exp(\theta_{jj}+\sum_{j^{\prime}\neq
j}\theta_{jj^{\prime}}y_{ij^{\prime}})}{1+\exp(\theta_{jj}+\sum_{j^{\prime}\neq
j}\theta_{jj^{\prime}}y_{ij^{\prime}})}.$

Let
$\widetilde{{\bf{y}}}_{j}=(\sqrt{v_{1j}}y_{1j}-\dot{{\bf{y}}}_{j},\ldots,\sqrt{v_{nj}}y_{nj}-\dot{{\bf{y}}}_{j})^{\top}$,
where $\dot{{\bf{y}}}_{j}=1/n\sum_{i=1}^{n}\sqrt{v_{ij}}y_{ij}$. We use
$\tilde{{\bf{X}}}=(\tilde{{\bf{X}}}_{(1,2)},\cdots,\tilde{{\bf{X}}}_{(p-1,p)})$
to denote an $np$ by $p(p-1)/2$ matrix, with

$\tilde{{\bf{X}}}_{(j,j^{{}^{\prime}})}=({\bf{0}}_{n},...,{\bf{0}}_{n},\underbrace{\widetilde{{\bf{y}}}_{j}}_{j-\textnormal{th
block}},{\bf{0}}_{n},...,{\bf{0}}_{n},\underbrace{\widetilde{{\bf{y}}}_{j^{{}^{\prime}}}}_{j^{{}^{\prime}}-\textnormal{th
block}},{\bf{0}}_{n},\ldots,{\bf{0}}_{n})^{\top},$

where ${\bf{0}}_{n}$ is an $n$-dimensional column vector of zeros. Let
$\tilde{{\bf{X}}}^{(i,j)}$ be the $[(j-1)n+i]$-th row of $\tilde{{\bf{X}}}$
and
$\tilde{{\bf{X}}}^{(i)}=(\tilde{{\bf{X}}}^{(i,1)},\ldots,\tilde{{\bf{X}}}^{(i,p)})$.
Let
${\bf{T}}={\mathbb{E}}(\tilde{{\bf{X}}}^{(i)}(\tilde{{\bf{X}}}^{(i)})^{\top})$
and
${\bf{T}}^{n}=1/n\sum_{i=1}^{n}\tilde{{\bf{X}}}^{(i)}(\tilde{{\bf{X}}}^{(i)})^{\top}$
as its sample counterpart.

Our results rely on Assumptions (A3)–(A4) and the following regularity
conditions:
Assumption B

  1. (B1)

There exist constants $\tau_{4},\tau_{5}\in(0,\infty)$ such that

$\Lambda_{\min}(\bar{{\bf{H}}}_{\mathcal{B}\mathcal{B}})\geq\tau_{4}~{}~{}~{}\text{and}~{}~{}~{}\Lambda_{\max}({\bf{T}})\leq\tau_{5}.$

  2. (B2)

There exists a constant $\delta\in(0,1]$, such that

$\|\bar{{\bf{H}}}_{\mathcal{B}^{c}\mathcal{B}}\left(\bar{{\bf{H}}}_{\mathcal{B}\mathcal{B}}\right)^{-1}\|_{\infty}\leq(1-\delta).$
(30)

Under these assumptions and (A3)–(A4), we have the following result:

######  Theorem 6.

Suppose Assumptions (A3)–(A4) and (B1)–(B2) are satisfied. Assume further that
$\rho_{1n}=O(\sqrt{\log p/n})$, $n>O(q^{3}\log p)$ as $n\rightarrow\infty$,
$\rho_{2n}=O(\sqrt{\log(p-H+1)/n})$,
$\rho_{2n}\leq\delta\rho_{1n}/(4(2-\delta)\tau_{3}\|\bar{\mbox{\boldmath$\theta$}}^{o}\|_{\infty})$,
and $\mbox{\boldmath$\epsilon$}=\mbox{\boldmath$0$}$. Then, there exist finite
constants $\check{C}(\bar{\mbox{\boldmath$\theta$}}^{o})$,
$\check{D}(\bar{{\bf{q}}}^{o})$, and $\eta$ such that the following events
hold with probability at least $1-O(\exp(-\eta\rho_{1n}^{2}n))$:

  * •

There exists a local minimizer
$({\widehat{\mbox{\boldmath$\theta$}}}^{o},{\widehat{{\bf{q}}}}^{o})$ of (3.2)
such that

$\max\left(\|\widehat{\mbox{\boldmath$\theta$}}^{o}-\bar{\mbox{\boldmath$\theta$}}^{o}\|,\|{\widehat{{\bf{q}}}}^{o}-\bar{{\bf{q}}}^{0}\|\right)\\\
\leq\max\left(\check{C}(\bar{\mbox{\boldmath$\theta$}}^{o})\rho_{1n}\sqrt{q},\check{D}(\bar{{\bf{q}}}^{o})\rho_{2n}\sqrt{\Psi(p,H,K)}\right),$
(31)

where $q$ and $\Psi(p,H,K)$ are defined in (23).

  * •

If $\min_{(i,j)\in\mathcal{B}}\bar{\theta}_{ij}^{o}\geq
2\check{C}(\bar{\mbox{\boldmath$\theta$}}^{o})\sqrt{q}\rho_{1n}$, then
${\widehat{\mbox{\boldmath$\theta$}}}^{o}_{\mathcal{B}^{c}}=0$.

Theorem 6 provides sufficient conditions on the quadruple $(n,p,H,K)$ and the
model parameters for the FBN to succeed in consistently estimating the
neighborhood of every node in the graph and communities simultaneously. We
note that if $H=1$ (no fairness) and $K=p$ (no clustering), we recover the
results of Ravikumar et al. (2010).

###  3.3 Consistency of Fair Community Labeling in Graphical Models

In this section, we aim to show that our algorithms recover the fair ground-
truth community structure in the graph. Let ${\widehat{{\bf{V}}}}$ and
$\bar{{\bf{V}}}$ contain the orthonormal eigenvectors corresponding to the $K$
leading eigenvalues of ${\widehat{{\bf{Q}}}}$ and $\bar{{\bf{Q}}}$ as columns,
respectively. It follows from (Lei et al., 2015, Lemma 2.1) that if any rows
of the matrix $\bar{{\bf{V}}}$ are same, then the corresponding nodes belong
to the same cluster. Consequently, we want to show that up to some orthogonal
transformation, the rows of ${\widehat{{\bf{V}}}}$ are close to the rows of
$\bar{{\bf{V}}}$ so that we can simply apply K-means clustering to the rows of
the matrix ${\widehat{{\bf{V}}}}$. In particular, we consider the K-means
approach Lei et al. (2015) defined as

$({\widehat{{\bf{U}}}},{\widehat{{\bf{O}}}})=\operatorname*{argmin}_{{\bf{U}},{\bf{O}}}\|{\bf{U}}{\bf{O}}-{\widehat{{\bf{V}}}}\|_{\textnormal{F}}^{2},~{}~{}\text{subject
to}~{}~{}{\bf{U}}\in\mathbb{M}_{p,K},~{}~{}{\bf{O}}\in{\mathbb{R}}^{K\times
K},$ (32)

where $\mathbb{M}_{p,K}$ is the set of $p\times K$ matrices that have on each
row a single $1$, indicating the fair community to which the node belongs, and
all other values on the row set at $0$, since a node belongs to only one
community. Finding a global minimizer for the problem (32) is NP-hard (Aloise
et al., 2009). However, there are polynomial time approaches (Kumar et al.,
2004) that find an approximate solution
$({\widehat{{\bf{U}}}},{\widehat{{\bf{O}}}})\in\mathbb{M}_{p,K}\times{\mathbb{R}}^{K\times
K}$ such that

$\|{\widehat{{\bf{U}}}}{\widehat{{\bf{O}}}}-{\widehat{{\bf{V}}}}\|_{\textnormal{F}}^{2}\leq(1+\xi)\operatorname*{argmin}_{({\bf{U}},{\bf{O}})\in\mathbb{M}_{p,K}\times{\mathbb{R}}^{K\times
K}}\|{\bf{U}}{\bf{O}}-{\widehat{{\bf{V}}}}\|_{\textnormal{F}}^{2}.$ (33)

Next, similar to (Lei et al., 2015, Theorem 3.1), we quantify the errors when
performing $(1+\xi)$-approximate K-means clustering on the rows of
$\widehat{{\bf{V}}}$ to estimate the community membership matrix. To do so,
let $\mathcal{E}_{k}$ denote the set of misclassified nodes from the $k$-th
community. By
$\bar{\mathcal{C}}=\cup_{k\in[K]}({\mathcal{C}_{k}}\backslash{\mathcal{E}_{k}})$
we denote the set of all nodes correctly classified across all communities and
by $\bar{{\bf{V}}}_{\bar{\mathcal{C}}}$ we denote the submatrix of
$\bar{{\bf{V}}}$ formed by retaining only the rows indexed by the set
$\bar{\mathcal{C}}$ of correctly classified nodes and all columns. Theorem 7
relates the sizes of the sets of misclassified nodes for each fair community
and specify conditions on the interplay between $n$, $p$, $H$, and $K$.

######  Theorem 7.

Let ${\widehat{{\bf{U}}}}$ be the output of $(1+\xi)$-approximate K-means
given in (33). If

$(2+\xi)\phi(p,H,K)\sqrt{\frac{K}{n}}<\pi$

for some constant $\pi>0$, then there exist subsets
$\mathcal{E}_{k}\subset\mathcal{C}_{k}$ for $k=1,\ldots,K$, and a permutation
matrix $\Phi$ such that
${\widehat{{\bf{V}}}}_{\bar{\mathcal{C}}}\mbox{\boldmath$\Phi$}=\bar{{\bf{V}}}_{\bar{\mathcal{C}}}$,
and

$\sum_{k=1}^{K}|\mathcal{E}_{k}|/|\mathcal{C}_{k}|\leq\pi^{-1}(2+\xi)\phi(p,H,K)\sqrt{\frac{K}{n}}$

with probability tending to 1.

##  4 Simulation Study

###  4.1 Tuning Parameter Selection

We consider a _Bayesian information criterion_ (BIC)-type quantity for tuning
parameter selection in (11). Recall from Section 2.1 that objective function
(11) decomposes the parameter of interest into
$(\mbox{\boldmath$\Theta$},{\bf{Q}})$ and places $\ell_{1}$ and trace
penalties on $\Theta$ and ${\bf{Q}}$, respectively. Specifically, for the
graphical Lasso, i.e., problem in (11) with $\rho_{2}=0$, Yuan and Lin (2006)
proposed to select the tuning parameter $\rho_{1}$ such that
$\hat{\mbox{\boldmath$\Theta$}}$ minimizes the following quantity:

$n\left(-\log\det(\hat{\mbox{\boldmath$\Theta$}})+\cdot\text{trace}(\mbox{\boldmath$S$}\hat{\mbox{\boldmath$\Theta$}})\right)+\log(n)\cdot|\hat{\mbox{\boldmath$\Theta$}}|.$

Here, $|\hat{\mbox{\boldmath$\Theta$}}|$ is the cardinality of
$\hat{\mbox{\boldmath$\Theta$}}$, i.e., the number of unique non-zeros in
$\hat{\mbox{\boldmath$\Theta$}}$. Using a similar idea, we consider minimizing
the following BIC-type criteria for selecting the set of tuning parameters
$(\rho_{1},\rho_{2})$ for (11):

$\mathrm{BIC}(\hat{\mbox{\boldmath$\Theta$}},\hat{{\bf{Q}}}):=\sum_{k=1}^{K}n_{k}\left(-\log|\hat{\mathbf{\Theta}}_{k}^{2}|+\text{trace}\left(\left((1+c)\mathbf{S}_{k}+\hat{{\bf{Q}}}_{k}\right)\hat{\mathbf{\Theta}}_{k}^{2}\right)\right)+\log(n_{k})\cdot|\hat{\mbox{\boldmath$\Theta$}}_{k}|),$

where $\hat{\mbox{\boldmath$\Theta$}}_{k}$ is the $k$-th estimated inverse
covariance matrix.

Note that when the constant $c$ is large,
$\mathrm{BIC}(\hat{\mbox{\boldmath$\Theta$}},\hat{{\bf{Q}}})$ will favor more
graph estimation in $\hat{\mbox{\boldmath$\Theta$}}$. Throughout, we take
$c=0.25$.

###  4.2 Notation and Measures of Performance

We define several measures of performance that will be used to numerically
compare the various methods. To assess the clustering performance, we compute
the following clustering error (CE) which calculates the distance between an
estimated community assignment $\hat{z}_{i}$ and the true assignment $z_{i}$
of the $i$th node:

$\textnormal{CE}:=\frac{1}{p}\left|\left\\{(i,j):{\bf{1}}(\hat{z}_{i}=\hat{z}_{j})\neq{\bf{1}}(z_{i}=z_{j}),i<j\right\\}\right|.$

To measure the estimation quality, we calculate the proportion of correctly
estimated edges (PCEE):

$\textnormal{PCEE}:=\frac{\sum_{j^{\prime}<j}\left(1_{\\{|\hat{\theta}_{jj^{\prime}}|>10^{-5}\text{
and }|\theta_{jj^{\prime}}|\neq
0\\}}\right)}{\sum_{j^{\prime}<j}\left(1_{\\{|\theta_{jj^{\prime}}|\neq
0\\}}\right)}.$

Finally, we use balance as a fairness metric to reflect the distribution of
the fair clustering (Chierichetti et al., 2017). Let
$\mathcal{N}_{i}=\\{j:r_{ij}=1\\}$ be the set of neighbors of node $i$ in
${\bf{R}}$. For a set of communities $\\{\mathcal{C}_{k}\\}_{k=1}^{K}$, the
balance coefficient is defined as

$\textnormal{Balance}:=\frac{1}{p}\sum_{i=1}^{p}\min_{k,\ell\in[K]}\frac{|\mathcal{C}_{k}\cap\mathcal{N}_{i}|}{|\mathcal{C}_{\ell}\cap\mathcal{N}_{i}|}.$

The balance coefficient, called simply the balance, is used to quantify how
well the selected edges can eliminate discrimination–the selected edges are
considered fairer if they can lead to a balanced community structure that
preserve proportions of protected attributes.

###  4.3 Data Generation

In order to demonstrate the performance of the proposed algorithms, we create
several synthetic datasets based a special random graph with community and
group structures. Then the baseline and proposed algorithms are used to
recover graphs (i.e., graph-based models) from the artificially generated
data. To create a dataset, we first construct a graph, then its associated
matrix, $\Theta$, is used to generate independent data samples from the
distribution $N(0,\mbox{\boldmath$\Theta$}^{-1})$. A graph (i.e., $\Theta$) is
constructed in two steps. In the first step, we determine the graph structure
(i.e., connectivity) based on the random modular graph also known as
stochastic block model (SBM) (Holland et al., 1983; Kleindessner et al.,
2019).

The SBM takes, as input, a function $\pi_{c}:[p]\rightarrow[K]$ that assigns
each vertex $i\in\mathcal{V}$ to one of the $K$ clusters. Then, independently,
for all node pairs $(i,j)$ such that $i>j$,
${\mathbb{P}}(a_{ij}=1)=b_{\pi_{c}(i)\pi_{c}(j)}$, where
${\bf{B}}\in[0,1]^{K\times K}$ is a symmetric matrix. Each $b_{k\ell}$
specifies the probability of a connection between two nodes that belong to
clusters $\mathcal{C}_{k}$ and $\mathcal{C}_{\ell}$, respectively. A commonly
used variant of SBM assumes $b_{kk}=\xi_{1}$ and $b_{k\ell}=\xi_{2}$ for all
$k,\ell\in[K]$ such that $k\neq\ell$. Let $\pi_{d}:[p]\rightarrow[H]$ be
another function that assigns each vertex $i\in\mathcal{V}$ to one of the $H$
protected groups. We consider a variant of SBM with the following
probabilities:

${\mathbb{P}}(a_{ij}=1)=\begin{cases}\zeta_{4}&\text{if
}\pi_{c}(i)=\pi_{c}(j)~{}\text{ and }~{}\pi_{d}(i)=\pi_{d}(j),\\\
\zeta_{3}&\text{if }\pi_{c}(i)\neq\pi_{c}(j)~{}\text{ and
}~{}\pi_{d}(i)=\pi_{d}(j),\\\ \zeta_{2}&\text{if
}\pi_{c}(i)=\pi_{c}(j)~{}\text{ and }~{}\pi_{d}(i)\neq\pi_{d}(j),\\\
\zeta_{1}&\text{if }\pi_{c}(i)\neq\pi_{c}(j)~{}\text{ and
}~{}\pi_{d}(i)\neq\pi_{d}(j).\end{cases}$ (34)

Here, $1\geq\zeta_{i+1}\geq\zeta_{i}\geq 0$ are probabilities used for
sampling edges. In our implementation, we set $\zeta_{i}=0.1i$ for all
$i=1,\ldots,4$. We note that when vertices $i$ and $j$ belong to the same
community, they have a higher probability of connection between them for a
fixed value of $\pi_{d}$; see, (Kleindessner et al., 2019) for further
discussions.

In the second step, the graph weights (i.e., node and edge weights) are
randomly selected based on a uniform distribution from the interval $[0.1,3]$
and the associated (Laplacian) matrix $\Theta$ is constructed. Finally, given
the graph matrix $\Theta$, we generate the data matrix ${\bf{Y}}$ according to
${\bf{y}}_{1},\ldots,{\bf{y}}_{n}\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}N(\mbox{\boldmath$0$},\mbox{\boldmath$\Theta$}^{-1})$.

###  4.4 Comparison to Graphical Lasso and Neighbourhood Selection Methods

We consider four setups for comparing our methods with community-based
graphical models (GM):

  1. I.

Three-stage approach for which we (i) use a GM to estimate precision matrix
$\hat{\mbox{\boldmath$\Theta$}}$, (ii) apply a community detection approach
(Cai et al., 2015) to compute partition matrix $\hat{{\bf{Q}}}$, and (iii)
employ a K-means clustering to obtain clusters.

  2. II.

Two-stage approach for which we (i) use (12) without fairness constraint to
simultaneously estimate precision and partition matrices and (ii) employ a
K-means clustering to obtain clusters.

  1. FI.

Three-stage approach for which we (i) use a GM to estimate precision matrix
$\hat{\mbox{\boldmath$\Theta$}}$, (ii) apply a community detection approach
(Cai et al., 2015) to compute partition matrix $\hat{{\bf{Q}}}$, and (iii)
employ a fair K-means clustering (Chierichetti et al., 2017) to obtain
clusters.

  2. FII.

Two-stage approach for which we (i) use (12) to simultaneously estimate
precision and partition matrices and (ii) employ the K-means clustering to
obtain clusters.

The main goal of Setups I. and II. is to compare the community detection
errors without fairness constraints under different settings of $L$ and $G$
functions.

We consider three type of GMs in Setups I.–FII.:

  1. A.

A graphical Lasso-type method (Friedman et al., 2008) implemented using
$L(\mbox{\boldmath$\Theta$};{\bf{Y}})=n/2[-\log\det(\mbox{\boldmath$\Theta$})+\operatorname{trace}({\bf{S}}\mbox{\boldmath$\Theta$})]$
and $G(\mbox{\boldmath$\Theta$})=\mbox{\boldmath$\Theta$}$.

  2. B.

A neighborhood selection-type method (Khare et al., 2015) implemented using
$L(\mbox{\boldmath$\Theta$};{\bf{Y}})=n/2[-\log|\text{diag}(\mbox{\boldmath$\Theta$}^{2})|+\text{trace}({\bf{S}}\mbox{\boldmath$\Theta$}^{2})]$
and $G(\mbox{\boldmath$\Theta$})=\mbox{\boldmath$\Theta$}^{2}$.

  3. C.

A neighborhood selection-type method (Ravikumar et al., 2010) implemented
using
$L(\mbox{\boldmath$\Theta$};{\bf{Y}})=-\sum_{j=1}^{p}\sum_{j^{\prime}=1}^{p}\theta_{jj^{\prime}}s_{jj^{\prime}}+1/n\sum_{i=1}^{n}\sum_{j=1}^{p}\log\big{(}1+\exp(\theta_{jj}+\sum_{j^{\prime}\neq
j}\theta_{jj^{\prime}}y_{ij^{\prime}})\big{)}$ and
$G(\mbox{\boldmath$\Theta$})=\mbox{\boldmath$\Theta$}$.

In our results, we name each one “GM-Type-Setup” to refer to the GM type and
the setup above. For example, GM-A.II. refers to a graphical Lasso used in the
first step of the two stage approach in II. It is worth mentioning that
GM-A.II. can be seen as a variant of the recently developed cluster-based
GLASSO (Pircalabelu and Claeskens, 2020; Kumar et al., 2020; Hosseini and Lee,
2016).

Sample Size | Method | CE | PCEE | Balance
---|---|---|---|---
$n=300$ | GM-B.I. | 0.464(0.015) | 0.773(0.010) | 0.143(0.005)
| GM-B.II. | 0.426(0.015) | 0.791(0.010) | 0.152(0.005)
| GM-B.FI. | 0.197(0.005) | 0.773(0.010) | 0.301(0.005)
| GM-B.FII.(FCONCORD) | 0.060(0.005) | 0.835(0.010) | 0.516(0.010)
$n=450$ | GM-B.I. | 0.423(0.015) | 0.813(0.010) | 0.187(0.009)
| GM-B.II. | 0.418(0.011) | 0.863(0.010) | 0.207(0.009)
| GM-B.FI. | 0.163(0.015) | 0.813(0.010) | 0.335(0.010)
| GM-B.FII.(FCONCORD) | 0.011(0.005) | 0.891(0.012) | 0.520(0.010)

Table 1: Simulation results of neighborhood selection-type GMs on SBM network.
The results are for $p=600$, $H=3$, and $K=2$.

We repeated the procedure 10 times and reported the averaged clustering
errors, proportion of correctly estimated edges, and balance. Tables 1 and 2
is for SBM with $p=300$. As shown in Tables 1 and 2, GM-A.I., GM-A.II.,
GM-B.I., and GM-B.I. have the largest clustering error and the proportion of
correctly estimated edges. GM-A.II. and GM-B.II. improve the performance of
GM-A.I. and GM-B.I. in the precision matrix estimation. However, they still
incur a relatively large clustering error and balance since they ignore the
similarity across different community matrices and employ a standard K-means
clustering. In contrast, our FCONCORD and FGLASSO algorithms achieve the best
clustering accuracy, estimation accuracy, and balance for both scenarios. This
is due to our simultaneous clustering and estimation strategy as well as the
consideration of the fairness of precision matrices across clusters. This
experiment shows that a satisfactory fair community detection algorithm is
critical to achieve accurate estimations of heterogeneous and fair GMs, and
alternatively good estimation of GMs can also improve the fair community
detection performance. This explains the success of our simultaneous method in
terms of both fair clustering and GM estimation.

Sample Size | Method | CE | PCEE | Balance
---|---|---|---|---
$n=300$ | GM-A.I. | 0.431(0.015) | 0.731(0.011) | 0.139(0.005)
| GM-A.II. | 0.416(0.015) | 0.761(0.011) | 0.167(0.005)
| GM-A.FI. | 0.182(0.006) | 0.731(0.011) | 0.322(0.010)
| GM-A.FII. (FGLASSO) | 0.097(0.006) | 0.815(0.050) | 0.504(0.010)
$n=450$ | GM-A.I. | 0.470(0.011) | 0.792(0.010) | 0.154(0.005)
| GM-A.II. | 0.406(0.011) | 0.803(0.010) | 0.187(0.015)
| GM-A.FI. | 0.163(0.005) | 0.792(0.010) | 0.366(0.005)
| GM-A.FII.(FGLASSO) | 0.050(0.005) | 0.863(0.012) | 0.513(0.005)

Table 2: Simulation results of graphical Lasso-type algorithms on SBM network.
The results are for $p=600$, $H=3$, and $K=2$.

Next, we consider a natural composition of SBM and the Ising model called
Stochastic Ising Block Model (SIBM); see, e.g., Berthet et al. (2016); Ye
(2020) for more details. In SIBM, we take SBM similar to (34), where $p$
vertices are divided into clusters and subgroups and the edges are connected
independently with probability $\\{\xi_{i}\\}_{i=1}^{4}$. Then, we use the
graph $\mathcal{G}$ generated by the SBM as the underlying graph of the Ising
model and draw $n$ i.i.d. samples from it. The objective is to exactly recover
the fair clusters in SIBM from the samples generated by the Ising model,
without observing the graph $\mathcal{G}$.

Tables 3 reports the averaged clustering errors, the proportion of correctly
estimated edges, and balance for SIBM with $p=100$. The standard GM-C.I.
method has the largest clustering error and balance due to its ignorance of
the network structure in the precision matrices. GM-C.FI. improves the
clustering performance of the GM-C.I. by using the method of (Ravikumar et
al., 2010) in the precision matrix estimation and the robust community
detection approach (Cai et al., 2015) for computing partition matrix
$\hat{{\bf{Q}}}$. GM-C.FII. is able to achieve the best balance and clustering
performance due to the procedure of simultaneous fair clustering and
heterogeneous GM estimation.

Sample Size | Method | CE | PCEE | Balance
---|---|---|---|---
$n=200$ | GM-C.I. | 0.489(0.011) | 0.601(0.010) | 0.104(0.005)
| GM-C.II. | 0.461(0.005) | 0.651(0.010) | 0.107(0.005)
| GM-C.FI. | 0.219(0.005) | 0.601(0.010) | 0.212(0.010)
| GM-C.FII.(FBN) | 0.117(0.009) | 0.738(0.005) | 0.385(0.010)
$n=400$ | GM-C.I. | 0.434(0.010) | 0.637(0.010) | 0.113(0.005)
| GM-C.II. | 0.461(0.010) | 0.681(0.010) | 0.107(0.005)
| GM-C.FI. | 0.219(0.005) | 0.637(0.010) | 0.382(0.010)
| GM-C.FII.(FBN) | 0.104(0.005) | 0.796(0.005) | 0.405(0.010)

Table 3: Simulation results of binary neighborhood selection-type GMs on SIBM
network. The results are for $p=100$, $H=2$, and $K=2$.

##  5 Real Data Application

Recommender systems (RS) model user-item interactions to provide personalized
item recommendations that will suit the user’s taste. Broadly speaking, two
types of methods are used in such systems—content based and collaborative
filtering. Content based approaches model interactions through user and item
covariates. Collaborative filtering (CF), on the other hand, refers to a set
of techniques that model user-item interactions based on user’s past response.

A popular class of methods in RS is based on clustering users and/or items
(Ungar and Foster, 1998; O’Connor and Herlocker, 1999; Sarwar et al., 2001;
Schafer et al., 2007). Indeed, it is more natural to model the users and the
items using clusters (communities), where each cluster includes a set of like-
minded users or the subset of items that they are interested in. The overall
procedure of this method, called cluster CF (CCF), contains two main steps.
First, it finds clusters of users and/or items, where each cluster includes a
group of like-minded users or a set of items that these users are particularly
interested in. Second, in each cluster, it applies traditional CF methods to
learn users’ preferences over the items within this cluster. Despite
efficiency and scalability of these methods, in many human-centric
applications, using CCF in its original form can result in unfavorable and
even harmful clustering and prediction outcomes towards some demographic
groups in the data.

It is shown in Schafer et al. (2007); Mnih and Salakhutdinov (2008) that using
item-item similarities based on “who-rated-what” information is strongly
correlated with how users explicitly rate items. Hence, using this information
as user covariates helps in improving predictions for explicit ratings.
Further, one can derive an item graph where edge weights represent movie
similarities that are based on global “who-rated-what” matrix (Kouki et al.,
2015; Wang et al., 2015; Agarwal et al., 2011; Mazumder and Agarwal, 2011).
Imposing sparsity on such a graph and finding its fair communities is
attractive since it is intuitive that an item is generally related to only a
few other items. This can be achieved through our fair GMs. Such a graph gives
a fair neighborhood structure that can also help better predict explicit
ratings. In addition to providing useful information to predict ratings, we
note that using who-rated-what information also provides information to study
the fair relationships among items based on user ratings.

The goal of our analysis is to understand the balance and prediction accuracy
of fair GMs on RS datasets as well as the relationships among the items in
these datasets. We compare the performance of our fair GMs implemented in the
framework of standard CCF and its fair K-means variant. In particular, we
consider the following algorithms:

  * •

FGLASSO (FCONCORD)+CF: A two-stage approach for which we first use FGLASSO
(FCONCORD) to obtain the fair clusters and then apply traditional CF to learn
users’ preferences over the items within each cluster. We set $\rho_{1}=1$,
$\rho_{2}=0.05$, $\gamma=0.01$, and $\epsilon=1e-3$ in our implementations.

  * •

CCF (Fair CCF): A two-stage approach for which we first use K-means (fair
K-means Chierichetti et al. (2017)) clustering to obtain the clusters and then
apply CF to learn users’ preferences within each cluster (Ungar and Foster,
1998).

###  5.1 MovieLens Data

We use the MovieLens 10K dataset222http://www.grouplens.org. Following
previous works (Koren, 2009; Kamishima et al., 2012; Chen et al., 2020), we
use year of the movie as a sensitive attribute and consider movies before 1991
as old movies. Those more recent are considered new movies. Koren (2009)
showed that the older movies have a tendency to be rated higher, perhaps
because only masterpieces have survived. When adopting year as a sensitive
attribute, we show that our fair graph-based RS enhances the neutrality from
this masterpiece bias. The clustering balance and root mean squared error
(RMSE) have been used to evaluate different modeling methods on this dataset.
Since reducing RMSE is the goal, statistical models assume the response
(ratings) to be Gaussian for this data (Kouki et al., 2015; Wang et al., 2015;
Agarwal et al., 2011).

Figure 1: RMSE (left) and Balance (right) of standard CCF, Fair CCF,
FGLASSO+CF, and FCONCORD+CF on MovieLens 10K data.

Experimental results are shown in Figure 1. As expected, the baseline with no
notion of fairness–CCF–results in the best overall RMSEs, with our two
approaches (FGLASSO+CF and FCONCORD+CF) providing performance fairly close to
CCF. Figure 1 (right) demonstrates that compared to fair CCF, FGLASSO+CF and
FCONCORD+CF significantly improve the clustering balance. Hence, our fair
graph-based RSs successfully enhanced the neutrality without seriously
sacrificing the prediction accuracy.

GM-B.II. The pair of movies Partial correlation The Godfather (1972) The
Godfather: Part II (1974) 0.592 Grumpy Old Men (1993) Grumpier Old Men (1995)
0.514 Patriot Games (1992) Clear and Present Danger (1994) 0.484 The Wrong
Trousers (1993) A Close Shave (1995) 0.448 Toy Story (1995) Toy Story 2 (1999)
0.431 Star Wars: Episode IV–A New Hope (1977) Star Wars: Episode V–The Empire
Strikes Back (1980) 0.415 GM-B.FII. (FCONCORD) The pair of movies Partial
correlation The Godfather (1972) The Godfather: Part II (1974) 0.534 Grumpy
Old Men (1993) Grumpier Old Men (1995) 0.520 Austin Powers: International Man
of Mystery (1997) Austin Powers: The Spy Who Shagged Me (1999) 0.491 Toy Story
(1995) Toy Story 2 (1999) 0.475 Patriot Games (1992) Clear and Present Danger
(1994) 0.472 The Wrong Trousers (1993) A Close Shave (1995) 0.453

Table 4: Pairs of movies with top 5 absolute values of partial correlations in
the precision matrix from GM-B.II. and GM-B.FII.(FCONCORD).

Fair GMs also provide information to study the relationships among items based
on user ratings. To illustrate this, the top-5 movie pairs with the highest
absolute values of partial correlations are shown in Table 4. If we look for
the highly related movies to a specific movie in the precision matrix, we find
that FCONCORD enhances the balance by assigning higher correlations to more
recent movies such as “The Wrong Trousers” (1993) and “A Close Shave” (1995).

Figure 2: Subgraphs of the precision matrices estimated by GM-B.II.(left) and
GM-B.FII.(right). Nodes represent the movies labeled by their titles. Circle
markers denote movies within a single community in each subgraph, and square
markers denote isolated movies. Blue nodes are new movies and red nodes old
movies within each community. The width of a link is proportional to the
magnitude of the corresponding partial correlations. GM-B.FII.(FCONCORD)
enhances the neutrality from year bias by replacing the old movies within each
community by new ones.

In addition, the estimated communities for two sub-graphs of movies are also
shown in Figure 2. From both networks, we can see that the estimated
communities mainly consist of mass marketed commercial movies, dominated by
action films. Note that these movies are usually characterized by high
production budgets, state of the art visual effects, and famous directors and
actors. Examples in this communities include “The Godfather” (1972),
“Terminator 2” (1991), and “Return of the Jedi” (1983), “Raiders of Lost Ark”
(1981), etc. As expected, movies within the same series are most strongly
associated. Figure 2 (right) shows that FCONCORD enhances the neutrality from
the old movies bias by replacing them with new ones such as “Jurassic Park
(1993),” “The Wrong Trousers” (1993) and “A Close Shave” (1995).

###  5.2 Music Data

Music RSs are designed to give personalized recommendations of songs,
playlists, or artists to a user, thereby reflecting and further complementing
individual users’ specific music preferences. Although accuracy metrics have
been widely applied to evaluate recommendations in music RS literature,
evaluating a user’s music utility from other impact oriented perspectives,
including their potential for discrimination, is still a novel evaluation
practice in the music RS literature (Epps-Darling et al., 2020; Chen et al.,
2020; Shakespeare et al., 2020). Next, we center our attention on artists’
gender bias for which we want to estimate if standard music RSs may exacerbate
its impact.

Figure 3: Input Preference Ratio (PR) distributions of LFM-360K dataset.

To illustrate the impact of artists gender bias in RSs, we use the freely
available binary LFM-360K music dataset333http://www.last.fm. The LFM-360K
consists of approximately 360,000 users listening histories from Last.fm
collected during Fall 2008. We generate recommendations for a sample of all
users for which gender can be identified. We limit the size of this sample to
be 10% randomly chosen of all male and female users in the whole dataset due
to computational constraints. Let $\mathcal{U}$ be the set of $n$ users,
$\mathcal{I}$ be the set of $p$ items and ${\bf{Y}}$ be the $n\times p$ input
matrix, where $y_{ui}=1$ if user $u$ has selected item $i$, and zero
otherwise. Given the matrix ${\bf{Y}}$, the input preference ratio (PR) for
user group $\mathcal{D}$ on item category $\mathcal{C}$ is the fraction of
liked items by group $\mathcal{D}$ in category $\mathcal{C}$, defined as the
following (Shakespeare et al., 2020):

$\textnormal{PR}(\mathcal{D},\mathcal{C}):=\frac{\sum_{u\in\mathcal{D}}\sum_{i\in\mathcal{C}}y_{ui}}{\sum_{u\in\mathcal{D}}\sum_{i\in\mathcal{I}}y_{ui}}.$
(35)

Figure 3 represents the distributions of users’ input PR towards male and
female artist groups. It shows that only around 20% of users have a PR towards
male artists lower than 0.8. On the contrary, 80% of users have a PR lower
than 0.2 towards female artists. This shows that commonly deployed state of
the art CF algorithms may act to further increase or decrease artist gender
bias in user-artist RS.

Figure 4: RMSE (left) and Balance (right) of standard CCF, Fair CCF, and
FBN+CF on LFM-360K music data.

Next, we study the balance and prediction accuracy of fair GMs on music RSs.
Figure 4 indicates that the proposed FBN+CF has the best performance in terms
of RMSE and balance. As expected, the baseline with no notion of
fairness–CCF–results in the best overall precision. Of the two fairness-aware
approaches, the fair K-means based approach–Fair CCF–performs considerably
below FBN+CF. This suggests that recommendation quality can be preserved, but
leaves open the question of whether we can improve fairness. Hence, we turn to
the impact on fairness of the three approaches. Figure 4(right) presents the
balance. We can see that fairness-aware approaches–Fair CCF and FBN+CF–have a
strong impact on the balance in comparison with standard CCF. And for RMSE, we
see that FBN+CF achieves much better ratings difference in comparison with
Fair CCF, indicating that we can induce aggregate statistics that are fair
between the two sides of the sensitive attribute (male vs. female).

##  6 Conclusion

In this work we developed a novel approach to learning fair GMs with community
structure. Our goal is to motivate a new line of work for fair community
learning in GMs that can begin to alleviate fairness concerns in this
important subtopic within unsupervised learning. Our optimization approach
used the demographic parity definition of fairness, but the framework is
easily extended to other definitions of fairness. We established statistical
consistency of the proposed method for both a Gaussian GM and an Ising model
proving that our method can recover the graphs and their fair communities with
high probability.

We applied the proposed framework to the tasks of estimating a Gaussian
graphical model and a binary network. The proposed framework can also be
applied to other types of graphical models, such as the Poisson graphical
model (Allen and Liu, 2012) or the exponential family graphical model (Yang et
al., 2012).

## Acknowledgment

D. Ataee Tarzanagh and L. Balzano were supported by NSF BIGDATA award
#1838179, ARO YIP award W911NF1910027, and NSF CAREER award CCF-1845076. A.O.
Hero was supported by US Army Resesarch Office grants #W911NF-15-1-0479 and
#W911NF-19-1-0269.

## References

  * Agarwal et al. (2011) Agarwal, D., L. Zhang, and R. Mazumder (2011).  Modeling item-item similarities for personalized recommendations on yahoo! front page.  The Annals of applied statistics, 1839–1875.
  * Allen and Liu (2012) Allen, G. I. and Z. Liu (2012).  A log-linear graphical model for inferring genetic networks from high-throughput sequencing data.  In 2012 IEEE International Conference on Bioinformatics and Biomedicine, pp. 1–6. IEEE.
  * Aloise et al. (2009) Aloise, D., A. Deshpande, P. Hansen, and P. Popat (2009).  Np-hardness of euclidean sum-of-squares clustering.  Machine learning 75(2), 245–248.
  * Amini et al. (2018) Amini, A. A., E. Levina, et al. (2018).  On semidefinite relaxations for the block model.  The Annals of Statistics 46(1), 149–179.
  * Banerjee et al. (2008) Banerjee, O., L. El Ghaoui, and A. d’Aspremont (2008).  Model selection through sparse maximum likelihood estimation for multivariate gaussian or binary data.  The Journal of Machine Learning Research 9, 485–516.
  * Barocas et al. (2019) Barocas, S., M. Hardt, and A. Narayanan (2019).  Fairness and Machine Learning.  fairmlbook.org.  http://www.fairmlbook.org.
  * Barzilai and Borwein (1988) Barzilai, J. and J. M. Borwein (1988).  Two-point step size gradient methods.  IMA journal of numerical analysis 8(1), 141–148.
  * Berthet et al. (2016) Berthet, Q., P. Rigollet, and P. Srivastava (2016).  Exact recovery in the Ising block model.  arXiv:1612.03880.
  * Boyd et al. (2011) Boyd, S., N. Parikh, E. Chu, B. Peleato, and J. Eckstein (2011).  Distributed optimization and statistical learning via the alternating direction method of multipliers.  Foundations and Trends® in Machine Learning 3(1), 1–122.
  * Burke et al. (2011) Burke, R., A. Felfernig, and M. H. Göker (2011).  Recommender systems: An overview.  Ai Magazine 32(3), 13–18.
  * Cai and Liu (2011) Cai, T. and W. Liu (2011).  Adaptive thresholding for sparse covariance matrix estimation.  Journal of the American Statistical Association 106(494), 672–684.
  * Cai et al. (2015) Cai, T. T., X. Li, et al. (2015).  Robust and computationally feasible community detection in the presence of arbitrary outlier nodes.  Annals of Statistics 43(3), 1027–1059.
  * Cardoso et al. (2020) Cardoso, J. V. d. M., J. Ying, and D. P. Palomar (2020).  Algorithms for learning graphs in financial markets.  arXiv preprint arXiv:2012.15410.
  * Caton and Haas (2020) Caton, S. and C. Haas (2020).  Fairness in machine learning: A survey.  arXiv preprint arXiv:2010.04053.
  * Celis et al. (2017) Celis, L. E., D. Straszak, and N. K. Vishnoi (2017).  Ranking with fairness constraints.  arXiv preprint arXiv:1704.06840.
  * Chen et al. (2020) Chen, J., H. Dong, X. Wang, F. Feng, M. Wang, and X. He (2020).  Bias and debias in recommender system: A survey and future directions.  arXiv preprint arXiv:2010.03240.
  * Chierichetti et al. (2017) Chierichetti, F., R. Kumar, S. Lattanzi, and S. Vassilvitskii (2017).  Fair clustering through fairlets.  In Advances in Neural Information Processing Systems, pp. 5029–5037.
  * Chouldechova and Roth (2018) Chouldechova, A. and A. Roth (2018).  The frontiers of fairness in machine learning.  arXiv preprint arXiv:1810.08810.
  * Danaher et al. (2014) Danaher, P., P. Wang, and D. M. Witten (2014).  The joint graphical lasso for inverse covariance estimation across multiple classes.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 76(2), 373–397.
  * Das et al. (2014) Das, J., P. Mukherjee, S. Majumder, and P. Gupta (2014).  Clustering-based recommender system using principles of voting theory.  In 2014 International conference on contemporary computing and informatics (IC3I), pp. 230–235. IEEE.
  * Donini et al. (2018) Donini, M., L. Oneto, S. Ben-David, J. Shawe-Taylor, and M. Pontil (2018).  Empirical risk minimization under fairness constraints.  arXiv preprint arXiv:1802.08626.
  * Eisenach et al. (2020) Eisenach, C., F. Bunea, Y. Ning, and C. Dinicu (2020).  High-dimensional inference for cluster-based graphical models.  Journal of machine learning research 21(53).
  * Epps-Darling et al. (2020) Epps-Darling, A., R. T. Bouyer, and H. Cramer (2020).  Artist gender representation in music streaming.  In Proceedings of the 21st International Society for Music Information Retrieval Conference (Montréal, Canada)(ISMIR 2020). ISMIR, pp. 248–254.
  * Feldman et al. (2015) Feldman, M., S. A. Friedler, J. Moeller, C. Scheidegger, and S. Venkatasubramanian (2015).  Certifying and removing disparate impact.  In proceedings of the 21th ACM SIGKDD international conference on knowledge discovery and data mining, pp. 259–268.
  * Friedman et al. (2008) Friedman, J., T. Hastie, and R. Tibshirani (2008).  Sparse inverse covariance estimation with the graphical lasso.  Biostatistics 9(3), 432–441.
  * Friedman et al. (2010) Friedman, J., T. Hastie, and R. Tibshirani (2010).  Applications of the lasso and grouped lasso to the estimation of sparse graphical models.  Technical report.
  * Gan et al. (2019) Gan, L., X. Yang, N. N. Nariestty, and F. Liang (2019).  Bayesian joint estimation of multiple graphical models.  Advances in Neural Information Processing Systems 32.
  * Gheche and Frossard (2020) Gheche, M. E. and P. Frossard (2020).  Multilayer clustered graph learning.  arXiv preprint arXiv:2010.15456.
  * Glassman et al. (2014) Glassman, E. L., R. Singh, and R. C. Miller (2014).  Feature engineering for clustering student solutions.  In Proceedings of the first ACM conference on Learning@ scale conference, pp. 171–172.
  * Grant and Boyd (2014) Grant, M. and S. Boyd (2014).  Cvx: Matlab software for disciplined convex programming, version 2.1.
  * Guo et al. (2015) Guo, J., J. Cheng, E. Levina, G. Michailidis, and J. Zhu (2015).  Estimating heterogeneous graphical models for discrete data with an application to roll call voting.  The Annals of Applied Statistics 9(2), 821.
  * Guo et al. (2010) Guo, J., E. Levina, G. Michailidis, and J. Zhu (2010).  Joint structure estimation for categorical markov networks.  Unpublished manuscript 3(5.2), 6.
  * Guo et al. (2011a) Guo, J., E. Levina, G. Michailidis, and J. Zhu (2011a).  Asymptotic properties of the joint neighborhood selection method for estimating categorical markov networks.  arXiv preprint math.PR/0000000.
  * Guo et al. (2011b) Guo, J., E. Levina, G. Michailidis, and J. Zhu (2011b).  Joint estimation of multiple graphical models.  Biometrika 98(1), 1–15.
  * Guo et al. (2011c) Guo, J., E. Levina, G. Michailidis, and J. Zhu (2011c).  Joint estimation of multiple graphical models.  Biometrika 98(1), 1–15.
  * Hao et al. (2018) Hao, B., W. W. Sun, Y. Liu, and G. Cheng (2018).  Simultaneous clustering and estimation of heterogeneous graphical models.  Journal of Machine Learning Research.
  * Hassner and Sklansky (1981) Hassner, M. and J. Sklansky (1981).  The use of markov random fields as models of texture.  In Image Modeling, pp. 185–198. Elsevier.
  * Holland et al. (1983) Holland, P. W., K. B. Laskey, and S. Leinhardt (1983).  Stochastic blockmodels: First steps.  Social networks 5(2), 109–137.
  * Hosseini and Lee (2016) Hosseini, M. J. and S.-I. Lee (2016).  Learning sparse gaussian graphical models with overlapping blocks.  In Advances in Neural Information Processing Systems, Volume 30, pp. 3801–3809.
  * Ising (1925) Ising, E. (1925).  Beitrag zur theorie des ferromagnetismus.  Zeitschrift für Physik A Hadrons and Nuclei 31(1), 253–258.
  * Kamishima et al. (2012) Kamishima, T., S. Akaho, H. Asoh, and J. Sakuma (2012).  Enhancement of the neutrality in recommendation.  Citeseer.
  * Karoui (2008) Karoui, N. E. (2008).  Operator norm consistent estimation of large-dimensional sparse covariance matrices.  The Annals of Statistics, 2717–2756.
  * Khare et al. (2015) Khare, K., S.-Y. Oh, and B. Rajaratnam (2015).  A convex pseudolikelihood framework for high dimensional partial correlation estimation with convergence guarantees.  Journal of the Royal Statistical Society: Series B: Statistical Methodology, 803–825.
  * Kleindessner et al. (2019) Kleindessner, M., S. Samadi, P. Awasthi, and J. Morgenstern (2019).  Guarantees for spectral clustering with fairness constraints.  In International Conference on Machine Learning, pp. 3458–3467. PMLR.
  * Koren (2009) Koren, Y. (2009).  Collaborative filtering with temporal dynamics.  In Proceedings of the 15th ACM SIGKDD international conference on Knowledge discovery and data mining, pp. 447–456.
  * Kouki et al. (2015) Kouki, P., S. Fakhraei, J. Foulds, M. Eirinaki, and L. Getoor (2015).  Hyper: A flexible and extensible probabilistic framework for hybrid recommender systems.  In Proceedings of the 9th ACM Conference on Recommender Systems, pp. 99–106.
  * Kumar et al. (2004) Kumar, A., Y. Sabharwal, and S. Sen (2004).  A simple linear time $(1+/splepsiv/)$-approximation algorithm for k-means clustering in any dimensions.  In 45th Annual IEEE Symposium on Foundations of Computer Science, pp. 454–462. IEEE.
  * Kumar et al. (2019) Kumar, S., J. Ying, J. V. Cardoso, de Miranda Cardoso, and D. Palomar (2019).  Structured graph learning via laplacian spectral constraints.  In Advances in Neural Information Processing Systems, pp. 11651–11663.
  * Kumar et al. (2020) Kumar, S., J. Ying, J. V. de Miranda Cardoso, and D. P. Palomar (2020).  A unified framework for structured graph learning via spectral constraints.  Journal of Machine Learning Research 21(22), 1–60.
  * Laferté et al. (2000) Laferté, J.-M., P. Pérez, and F. Heitz (2000).  Discrete markov image modeling and inference on the quadtree.  IEEE Transactions on image processing 9(3), 390–404.
  * Lafit et al. (2019) Lafit, G., F. Tuerlinckx, I. Myin-Germeys, and E. Ceulemans (2019).  A partial correlation screening approach for controlling the false positive rate in sparse gaussian graphical models.  Scientific reports 9(1), 1–24.
  * Lee et al. (2007) Lee, S.-I., V. Ganapathi, and D. Koller (2007).  Efficient structure learning of Markov networks using $\ell_{1}$-regularization.  In Advances in neural Information processing systems, pp. 817–824.
  * Lee and Liu (2015) Lee, W. and Y. Liu (2015).  Joint estimation of multiple precision matrices with common structures.  The Journal of Machine Learning Research 16(1), 1035–1062.
  * Lei et al. (2015) Lei, J., A. Rinaldo, et al. (2015).  Consistency of spectral clustering in stochastic block models.  Annals of Statistics 43(1), 215–237.
  * Li et al. (2021) Li, X., Y. Chen, and J. Xu (2021).  Convex relaxation methods for community detection.  Statistical Science 36(1), 2–15.
  * Liljeros et al. (2001) Liljeros, F., C. R. Edling, L. A. N. Amaral, H. E. Stanley, and Y. Åberg (2001).  The web of human sexual contacts.  Nature 411(6840), 907–908.
  * Ma and Michailidis (2016) Ma, J. and G. Michailidis (2016).  Joint structural estimation of multiple graphical models.  The Journal of Machine Learning Research 17(1), 5777–5824.
  * Manning and Schutze (1999) Manning, C. and H. Schutze (1999).  Foundations of statistical natural language processing.  MIT press.
  * Marlin and Murphy (2009) Marlin, B. M. and K. P. Murphy (2009).  Sparse gaussian graphical models with unknown block structure.  In Proceedings of the 26th Annual International Conference on Machine Learning, pp. 705–712.
  * Mazumder and Agarwal (2011) Mazumder, R. and D. K. Agarwal (2011).  A flexible, scalable and efficient algorithmic framework for primal graphical lasso.  arXiv preprint arXiv:1110.5508.
  * Meinshausen et al. (2006) Meinshausen, N., P. Bühlmann, et al. (2006).  High-dimensional graphs and variable selection with the lasso.  Annals of statistics 34(3), 1436–1462.
  * Mnih and Salakhutdinov (2008) Mnih, A. and R. R. Salakhutdinov (2008).  Probabilistic matrix factorization.  In Advances in neural information processing systems, pp. 1257–1264.
  * Oneto and Chiappa (2020) Oneto, L. and S. Chiappa (2020).  Fairness in machine learning.  In Recent Trends in Learning From Data, pp. 155–196. Springer.
  * O’Connor and Herlocker (1999) O’Connor, M. and J. Herlocker (1999).  Clustering items for collaborative filtering.  In Proceedings of the ACM SIGIR workshop on recommender systems, Volume 128. Citeseer.
  * Parimbelli et al. (2018) Parimbelli, E., S. Marini, L. Sacchi, and R. Bellazzi (2018).  Patient similarity for precision medicine: A systematic review.  Journal of biomedical informatics 83, 87–96.
  * Peng et al. (2009) Peng, J., P. Wang, N. Zhou, and J. Zhu (2009).  Partial correlation estimation by joint sparse regression models.  Journal of the American Statistical Association 104(486), 735–746.
  * Pham et al. (2011) Pham, M. C., Y. Cao, R. Klamma, and M. Jarke (2011).  A clustering approach for collaborative filtering recommendation using social network analysis.  J. UCS 17(4), 583–604.
  * Pircalabelu and Claeskens (2020) Pircalabelu, E. and G. Claeskens (2020).  Community-based group graphical lasso.  Journal of Machine Learning Research 21(64), 1–32.
  * Ravikumar et al. (2010) Ravikumar, P., M. J. Wainwright, and J. D. Lafferty (2010).  High-dimensional Ising model selection using $\ell_{1}$-regularized logistic regression.  The Annals of Statistics 38(3), 1287–1319.
  * Ravikumar et al. (2010) Ravikumar, P., M. J. Wainwright, J. D. Lafferty, et al. (2010).  High-dimensional ising model selection using $\ell_{1}$-regularized logistic regression.  The Annals of Statistics 38(3), 1287–1319.
  * Robins et al. (2007) Robins, G., P. Pattison, Y. Kalish, and D. Lusher (2007).  An introduction to exponential random graph $p*$ models for social networks.  Social networks 29(2), 173–191.
  * Rocha et al. (2008) Rocha, G. V., P. Zhao, and B. Yu (2008).  A path following algorithm for sparse pseudo-likelihood inverse covariance estimation (splice).  arXiv preprint arXiv:0807.3734.
  * Samadi et al. (2018) Samadi, S., U. Tantipongpipat, J. Morgenstern, M. Singh, and S. Vempala (2018).  The price of fair pca: One extra dimension.  arXiv preprint arXiv:1811.00103.
  * Sarwar et al. (2001) Sarwar, B., G. Karypis, J. Konstan, and J. Riedl (2001).  Item-based collaborative filtering recommendation algorithms.  In Proceedings of the 10th international conference on World Wide Web, pp. 285–295.
  * Schafer et al. (2007) Schafer, J. B., D. Frankowski, J. Herlocker, and S. Sen (2007).  Collaborative filtering recommender systems.  In The adaptive web, pp. 291–324. Springer.
  * Shakespeare et al. (2020) Shakespeare, D., L. Porcaro, E. Gómez, and C. Castillo (2020).  Exploring artist gender bias in music recommendation.  arXiv preprint arXiv:2009.01715.
  * Song et al. (2011) Song, Q., J. Ni, and G. Wang (2011).  A fast clustering-based feature subset selection algorithm for high-dimensional data.  IEEE transactions on knowledge and data engineering 25(1), 1–14.
  * Tan et al. (2014) Tan, K. M., P. London, K. Mohan, S.-I. Lee, M. Fazel, and D. M. Witten (2014).  Learning graphical models with hubs.  Journal of Machine Learning Research 15(1), 3297–3331.
  * Tan et al. (2013) Tan, P.-N., M. Steinbach, and V. Kumar (2013).  Data mining cluster analysis: basic concepts and algorithms.  Introduction to data mining, 487–533.
  * Tantipongpipat et al. (2019) Tantipongpipat, U., S. Samadi, M. Singh, J. H. Morgenstern, and S. S. Vempala (2019).  Multi-criteria dimensionality reduction with applications to fairness.  Advances in neural information processing systems (32).
  * Tarzanagh and Michailidis (2018) Tarzanagh, D. A. and G. Michailidis (2018).  Estimation of graphical models through structured norm minimization.  Journal of Machine Learning Research 18(209), 1–48.
  * Ungar and Foster (1998) Ungar, L. H. and D. P. Foster (1998).  Clustering methods for collaborative filtering.  In AAAI workshop on recommendation systems, Volume 1, pp. 114–129. Menlo Park, CA.
  * Wang et al. (2015) Wang, H., N. Wang, and D.-Y. Yeung (2015).  Collaborative deep learning for recommender systems.  In Proceedings of the 21th ACM SIGKDD international conference on knowledge discovery and data mining, pp. 1235–1244.
  * Wang et al. (2009) Wang, P., D. L. Chao, and L. Hsu (2009).  Learning networks from high dimensional binary data: An application to genomic instability data.  arXiv preprint arXiv:0908.3882.
  * Wang et al. (2019) Wang, Y., W. Yin, and J. Zeng (2019).  Global convergence of admm in nonconvex nonsmooth optimization.  Journal of Scientific Computing 78(1), 29–63.
  * Xue et al. (2012) Xue, L., S. Ma, and H. Zou (2012).  Positive-definite $\ell_{1}$-penalized estimation of large covariance matrices.  Journal of the American Statistical Association 107(500), 1480–1491.
  * Yang et al. (2012) Yang, E., P. Ravikumar, G. I. Allen, and Z. Liu (2012).  Graphical models via generalized linear models.  In NIPS, Volume 25, pp. 1367–1375.
  * Ye (2020) Ye, M. (2020).  Exact recovery and sharp thresholds of stochastic ising block model.  arXiv preprint arXiv:2004.05944.
  * Yu et al. (2015) Yu, Y., T. Wang, and R. J. Samworth (2015).  A useful variant of the davis–kahan theorem for statisticians.  Biometrika 102(2), 315–323.
  * Yuan and Lin (2006) Yuan, M. and Y. Lin (2006).  Model selection and estimation in regression with grouped variables.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 68(1), 49–67.
  * Zafar et al. (2017) Zafar, M. B., I. Valera, M. G. Rogriguez, and K. P. Gummadi (2017).  Fairness constraints: Mechanisms for fair classification.  In Artificial Intelligence and Statistics, pp. 962–970. PMLR.
  * Zhao and Yu (2006) Zhao, P. and B. Yu (2006).  On model selection consistency of lasso.  The Journal of Machine Learning Research 7, 2541–2563.

Supplementary Material to “Fair Structure Learning in Heterogeneous Graphical
Models”

Davoud Ataee Tarzanagh, Laura Balzano, and Alfred O. Hero

  * •

Appendix A.1 provides some preliminaries used in the proof of main theorems.

  * •

Appendix A.2 provides large sample properties of FCONCORD, i.e., the proof of
Theorem 5.

  * •

Appendix A.3 provides large sample properties of FBN, i.e., the proof of
Theorem 6.

  * •

Appendix A.4 gives the consistency of fair community labeling in graphical
models.

  * •

Appendix A.5 provides the detailed derivation of the updates for Algorithm 1.

###  A.1 Preliminaries

Let the notation
${F}_{n}(\mbox{\boldmath$\theta$}^{d},\mbox{\boldmath$\theta$}^{o},{\bf{q}}^{d},{\bf{q}}^{o};{\bf{Y}})$
stands for $\frac{{F}}{n}$ in (3.1). We introduce a restricted version of
criterion (3.1) as below:

$\operatorname*{minimize}_{\mbox{\boldmath$\theta$}^{o},\mbox{\boldmath$q$}^{o}}~{}~{}F_{n}(\bar{\mbox{\boldmath$\theta$}}^{d},\mbox{\boldmath$\theta$}^{o},\bar{{\bf{q}}}^{d},{\bf{q}}^{o};{\bf{Y}})+\rho_{1n}\|\mbox{\boldmath$\theta$}^{o}\|_{1},~{}~{}~{}\text{subj.
to}~{}~{}\mbox{\boldmath$\theta$}^{o}_{\mathcal{B}^{c}}=0.$ (A.1)

We define a linear operator
$\mathcal{A}:{\mathbb{R}}^{p(p-1)/2}\rightarrow{\mathbb{R}}^{p\times p}$,
${\bf{w}}\rightarrow\mathcal{A}{\bf{w}}$, satisfying

$[\mathcal{A}{\bf{w}}]_{ij}=\begin{cases}-w_{i+d_{j}}\;&\;\;i>j,\\\
[\mathcal{A}{\bf{w}}]_{ji}\;&\;\;i<j,\\\ \sum_{i\neq
j}[\mathcal{A}{\bf{w}}]_{ij}\;&\;\;i=j,\end{cases}$ (A.2)

where $d_{j}=-j+\frac{j-1}{2}(2p-j).$ An example for $\mathcal{A}{\bf{w}}$ on
a vector ${\bf{w}}=[w_{1},w_{2},\cdots,w_{6}]^{\top}$ is given below

$\mathcal{A}{\bf{w}}=\left[\begin{array}[]{cccc}0&w_{1}&w_{2}&w_{3}\\\
w_{1}&0&w_{4}&w_{5}\\\ w_{2}&w_{4}&0&w_{6}\\\
w_{3}&w_{5}&w_{6}&0\end{array}\right].$

We derive the adjoint operator $\mathcal{A}^{*}$ of $\mathcal{A}$ by making
$\mathcal{A}^{*}$ satisfy
$\langle\mathcal{A}{\bf{w}},{\bf{z}}\rangle=\langle{\bf{w}},\mathcal{A}^{*}{\bf{z}}\rangle$;
see, Kumar et al. (2020) for more details.

Let ${\bf{Q}}^{o}=\text{diag}(\bar{{\bf{q}}}^{d})+\mathcal{A}{\bf{q}}^{o}$ and
$\mbox{\boldmath$\Theta$}^{o}=\text{diag}(\bar{\mbox{\boldmath$\theta$}}^{d})+\mathcal{A}\mbox{\boldmath$\theta$}^{o}$.
Since by our assumption $\mbox{\boldmath$\epsilon$}={\bf{0}}$, we obtain

$\displaystyle~{}\begin{array}[]{ll}\underset{\mbox{\boldmath$\theta$}^{o},\mbox{\boldmath$q$}^{o}}{\text{minimize}}&\begin{array}[]{c}\hskip
14.22636ptL_{n}(\bar{\mbox{\boldmath$\theta$}}^{d},\mbox{\boldmath$\theta$}^{o};{\bf{Y}})+\rho_{1n}\|\mbox{\boldmath$\theta$}^{o}\|_{1}+\rho_{2n}\operatorname{trace}(\mbox{\boldmath$\Theta$}^{o}{\bf{Q}}^{o}\mbox{\boldmath$\Theta$}^{o}),\end{array}\\\
\text{subj. to}&\begin{array}[t]{l}\hskip
14.22636pt\mbox{\boldmath$\theta$}^{o}_{\mathcal{B}^{c}}=0,~{}~{}{\bf{Q}}^{o}\succeq{\bf{0}},~{}~{}{\bf{A}}_{1}{\bf{Q}}^{o}={\bf{0}},~{}~{}{\bf{0}}\leq{\bf{Q}}^{o}\leq{\bf{J}}_{p}.\end{array}\end{array}$
(A.7)

It is easy to see that $\textnormal{rank}({\bf{A}}_{1})=H-1$. Let
${\bf{N}}\in{\mathbb{R}}^{p\times(p-H+1)}$ be a matrix whose rows form an
orthonormal basis of the nullspace of ${\bf{A}}_{1}$. We can substitute
${\bf{Q}}^{o}={\bf{N}}{\bf{R}}^{o}{\bf{N}}^{\top}$ for
${\bf{R}}^{o}\in{\mathbb{R}}^{(p-H+1)\times(p-H+1)}$, and then, using that
${\bf{N}}^{\top}{\bf{N}}={\bf{I}}_{(p-H+1)}$, Problem (A.7) becomes

$\displaystyle~{}\begin{array}[]{ll}\underset{\mbox{\boldmath$\theta$}^{o},\mbox{\boldmath$r$}^{o}}{\text{minimize}}&\begin{array}[]{c}\hskip
14.22636ptF_{n}(\bar{\mbox{\boldmath$\theta$}}^{d},\mbox{\boldmath$\theta$}^{o},\bar{{\bf{r}}}^{d},{\bf{r}}^{o};{\bf{Y}})+\rho_{1n}\|\mbox{\boldmath$\theta$}^{o}\|_{1}\end{array}\\\
\text{subj. to}&\begin{array}[t]{l}\hskip
14.22636pt\mbox{\boldmath$\theta$}^{o}_{\mathcal{B}^{c}}=0,~{}~{}{\bf{Q}}^{o}\succeq{\bf{0}},~{}~{}{\bf{0}}\leq{\bf{Q}}^{o}\leq{\bf{J}}_{p},\end{array}\end{array}$
(A.12)

where
$F_{n}(\bar{\mbox{\boldmath$\theta$}}^{d},\mbox{\boldmath$\theta$}^{o},\bar{{\bf{r}}}^{d},{\bf{r}}^{o};{\bf{Y}})=L_{n}(\bar{\mbox{\boldmath$\theta$}}^{d},\mbox{\boldmath$\theta$}^{o};{\bf{Y}})+\rho_{2n}\operatorname{trace}(\mbox{\boldmath$\Theta$}^{o}{\bf{N}}{\bf{R}}^{o}{\bf{N}}^{\top}\mbox{\boldmath$\Theta$}^{o})$.

Throughout, we use $\bar{{\bf{g}}}^{n}$ and $\bar{{\bf{H}}}^{n}$ to denote the
gradient and the Hessian of
$L_{n}(\bar{\mbox{\boldmath$\theta$}}^{d},\mbox{\boldmath$\theta$}^{o};{\bf{Y}})$.
We also define the population gradient and Hessian as follows: For $1\leq
i<j\leq p$

$\bar{{\bf{g}}}_{ij}:={\mathbb{E}}_{\bar{\mbox{\boldmath$\theta$}}^{o}}\left(\frac{\partial
L(\bar{\mbox{\boldmath$\theta$}}^{d},\mbox{\boldmath$\theta$}^{o},{\bf{Y}})}{\partial\theta^{o}_{ij}}\Bigl{|}_{\mbox{\boldmath$\theta$}^{o}=\bar{\mbox{\boldmath$\theta$}}^{o}}\right),$

and for $1\leq i<j\leq p$ and $1\leq t<s\leq p$,

$\bar{{\bf{H}}}_{ij,ts}:={\mathbb{E}}_{\bar{\mbox{\boldmath$\theta$}}^{o}}\Big{(}\frac{\partial^{2}L({\bar{\mbox{\boldmath$\theta$}}}^{d},{\mbox{\boldmath$\theta$}}^{o};{\bf{Y}})}{\partial\theta_{ij}^{o}\partial\theta_{ts}^{o}}\Big{|}_{\mbox{\boldmath$\theta$}^{o}=\bar{\mbox{\boldmath$\theta$}}^{o}}\Big{)}.$

###  A.2 Large Sample Properties of FCONCORD

We list some properties of the loss function.

######  Lemma 8.

(Peng et al., 2009) The following is true for the loss function:

  1. (L1)

There exist constants $0<M_{1}\leq M_{2}<\infty$ such that

$M_{1}(\bar{\mbox{\boldmath$\theta$}}^{o})\leq\Lambda_{\min}(\bar{{\bf{H}}})\leq\Lambda_{\max}(\bar{{\bf{H}}})\leq
M_{2}(\bar{\mbox{\boldmath$\theta$}}^{o}).$

  2. (L2)

There exists a constant $M_{3}(\bar{\mbox{\boldmath$\theta$}}^{o})<\infty$
such that for all $1\leq i<j\leq p$, $\bar{{\bf{H}}}_{ij,ij}\leq
M_{3}(\bar{\mbox{\boldmath$\theta$}}^{o})$.

  3. (L3)

There exist constants $M_{4}(\bar{\mbox{\boldmath$\theta$}}^{o})$ and
$M_{5}(\bar{\mbox{\boldmath$\theta$}}^{o})<\infty$, such that for any $1\leq
i<j\leq p$

$\operatorname{Var}_{\bar{\mbox{\boldmath$\theta$}}^{o}}(\bar{{\bf{g}}}^{n}_{ij})\leq
M_{4}(\bar{\mbox{\boldmath$\theta$}}^{o}),~{}~{}~{}~{}~{}\operatorname{Var}_{\bar{\mbox{\boldmath$\theta$}}^{o}}(\bar{{\bf{H}}}^{n}_{ij,ij})\leq
M_{5}(\bar{\mbox{\boldmath$\theta$}}^{o}).$

  4. (L4)

There exists a constant $0<M_{6}(\bar{\mbox{\boldmath$\theta$}}^{o})<\infty$,
such that for all $(i,j)\in\mathcal{B}$

$\bar{{\bf{H}}}_{ij,ij}-\bar{{\bf{H}}}_{ij,\mathcal{B}_{ij}}\bar{{\bf{H}}}_{\mathcal{B}_{ij},\mathcal{B}_{ij}}^{-1}\bar{{\bf{H}}}_{\mathcal{B}_{ij},ij}\geq
M_{4}(\bar{\mbox{\boldmath$\theta$}}^{o}),\quad\textnormal{where}\quad\mathcal{B}_{ij}:=\mathcal{B}/\\{(i,j)\\}.$

  5. (L5)

There exists a constant $M_{7}(\bar{\mbox{\boldmath$\theta$}}^{o})<\infty$,
such that for any $(i,j)\in\mathcal{B}^{c}$

$\|\bar{{\bf{H}}}_{ij,\mathcal{B}}\bar{{\bf{H}}}_{\mathcal{B},\mathcal{B}}^{-1}\|\leq
M_{7}(\bar{\mbox{\boldmath$\theta$}}^{o}).$ (A.13)

######  Lemma 9.

(Peng et al., 2009) Suppose Assumptions (A1)–(A2) hold, then for any $\eta>0$,
there exist constants $c_{0}$–$c_{3}$, such that for any
${\bf{v}}\in{\mathbb{R}}^{q}$ the following events hold with probability at
least $1-O(\exp(-\eta\log p))$ for sufficiently large $n$:

  1. (L1)

$\|\bar{{\bf{g}}}^{n}_{\mathcal{B}}\|\leq c_{0}\sqrt{q\frac{\log p}{n}}$.

  2. (L2)

$|{\bf{v}}^{\top}\bar{{\bf{g}}}^{n}_{\mathcal{B}}|\leq
c_{1}\|{\bf{v}}\|\sqrt{q\frac{\log p}{n}}$.

  3. (L3)

$|{\bf{v}}^{\top}(\bar{{\bf{H}}}^{n}_{\mathcal{B},\mathcal{B}}-\bar{{\bf{H}}}_{\mathcal{B},\mathcal{B}}){\bf{v}}|\leq
c_{2}\|{\bf{v}}\|^{2}q\sqrt{\frac{\log p}{n}}$.

  4. (L4)

$\|(\bar{{\bf{H}}}^{n}_{\mathcal{B},\mathcal{B}}-\bar{{\bf{H}}}_{\mathcal{B},\mathcal{B}}){\bf{v}}\|\leq
c_{3}\|{\bf{v}}\|q\sqrt{\frac{\log p}{n}}$.

######  Lemma 10.

Suppose Assumptions (A1)–(A4) are satisfied. Assume further that
$\rho_{1n}=O(\sqrt{\log p/n})$, $n=O(q\log(p))$,
$\rho_{2n}=O(\sqrt{\log(p-H+1)/n})$, and
$\mbox{\boldmath$\epsilon$}=\mbox{\boldmath$0$}$. Then, there exist finite
constants $C_{1}(\bar{\mbox{\boldmath$\theta$}}^{o})$ and
$D_{1}(\bar{{\bf{q}}}^{o})$, such that for any $\eta>0$, there exists a
(local) minimizer of the restricted problem (A.1) within the disc:

$\displaystyle\Big{\\{}(\widehat{\mbox{\boldmath$\theta$}}^{o},\widehat{{\bf{q}}}^{o}):\max\left(\|\widehat{\mbox{\boldmath$\theta$}}^{o}_{\mathcal{B}}-\bar{\mbox{\boldmath$\theta$}}^{o}_{\mathcal{B}}\|,\|{\widehat{{\bf{q}}}}^{o}-\bar{{\bf{q}}}^{o}\|\right)\leq\max\left(C_{1}(\bar{\mbox{\boldmath$\theta$}}^{o})\rho_{1n}\sqrt{q},D_{1}(\bar{{\bf{q}}}^{o})\rho_{2n}\sqrt{\Psi(p,H,K)}\right)\Big{\\}}$
(A.14)

with probability at least $1-O(\exp(-\eta\log p))$ for sufficiently large $n$.

###### Proof.

Let $\mu_{1n}=\rho_{1n}\sqrt{q}$ with $q=|\mathcal{B}|$ and
$\mu_{2n}=\rho_{2n}\sqrt{\Psi(p,H,K)}$. Let $C_{1}>0$ and
${\bf{w}}\in\mathbb{R}^{p(p-1)/2}$ such that ${\bf{w}}_{\mathcal{B}^{c}}=0$,
$\|{\bf{w}}\|_{2}=C_{1}$. Further, assume
${\bf{z}}\in\mathbb{R}^{(p-H+1)(p-H)/2}$ be an arbitrary vector with finite
entries and $\|{\bf{z}}\|=D_{1}$. Let

$F_{n}(\bar{\mbox{\boldmath$\theta$}}^{d},\bar{\mbox{\boldmath$\theta$}}^{o}+\mu_{1n}{\bf{w}},\bar{{\bf{r}}}^{d},\bar{{\bf{r}}}^{o}+\mu_{1n}{\bf{z}};{\bf{Y}})-F_{n}(\bar{\mbox{\boldmath$\theta$}}^{d},\bar{\mbox{\boldmath$\theta$}}^{o},\bar{{\bf{r}}}^{d},\bar{{\bf{r}}}^{o};{\bf{Y}})=I_{1}+I_{2}+I_{3},$

where

$\displaystyle I_{1}$
$\displaystyle:=L_{n}(\bar{\mbox{\boldmath$\theta$}}^{d},\bar{\mbox{\boldmath$\theta$}}^{o}+\mu_{1n}{\bf{w}};{\bf{Y}})-L_{n}(\bar{\mbox{\boldmath$\theta$}}^{d},\bar{\mbox{\boldmath$\theta$}}^{o};{\bf{Y}}),$
$\displaystyle I_{2}$
$\displaystyle:=\rho_{1n}(\|\bar{\mbox{\boldmath$\theta$}}^{o}+\mu_{1n}{\bf{w}}\|_{1}-\|\bar{\mbox{\boldmath$\theta$}}^{o}\|_{1}),~{}~{}\textnormal{and}$
$\displaystyle I_{3}$
$\displaystyle:=\rho_{2n}\operatorname{trace}\left((\bar{\mbox{\boldmath$\Theta$}}+\mu_{1n}\mathcal{A}{\bf{w}})^{2}({\bf{S}}+{\bf{N}}(\bar{{\bf{R}}}+\mu_{2n}\mathcal{A}{\bf{z}}){\bf{N}}^{\top})-\bar{\mbox{\boldmath$\Theta$}}^{2}\bar{{\bf{Q}}}\right).$

Following Peng et al. (2009), we first provide lower bounds for $I_{1}$ and
$I_{2}$. For the term $I_{1}$, it follows from Lemma 9 that

$\displaystyle{\bf{w}}_{\mathcal{B}}^{\top}\bar{{\bf{H}}}_{\mathcal{B},\mathcal{B}}{\bf{w}}_{\mathcal{B}}\geq\Lambda_{\min}(\bar{{\bf{H}}}_{\mathcal{B},\mathcal{B}})\|{\bf{w}}_{\mathcal{B}}\|_{2}^{2}\geq
M_{1}C^{2}_{1},$ (A.15)

which together with Lemma 9 gives

$\displaystyle I_{1}$
$\displaystyle=\mu_{1n}{\bf{w}}^{\top}_{\mathcal{B}}\bar{{\bf{g}}}^{n}_{\mathcal{B}}+\frac{1}{2}\mu_{1n}^{2}{\bf{w}}^{\top}_{\mathcal{B}}\bar{{\bf{H}}}^{n}_{\mathcal{B},\mathcal{B}}{\bf{w}}_{\mathcal{B}}$
$\displaystyle=\mu_{1n}({\bf{w}}_{\mathcal{B}})^{\top}\bar{{\bf{g}}}^{n}_{\mathcal{B}}+\frac{1}{2}\mu_{1n}^{2}{\bf{w}}_{\mathcal{B}}^{\top}\bar{{\bf{H}}}^{n}_{\mathcal{B},\mathcal{B}}{\bf{w}}_{\mathcal{B}}$
$\displaystyle=\mu_{1n}{\bf{w}}_{\mathcal{B}}^{\top}\bar{{\bf{g}}}^{n}_{\mathcal{B}}+\frac{1}{2}\mu_{1n}^{2}{\bf{w}}_{\mathcal{B}}^{\top}(\bar{{\bf{H}}}^{n}_{\mathcal{B},\mathcal{B}}-\bar{{\bf{H}}}_{\mathcal{B},\mathcal{B}}){\bf{w}}_{\mathcal{B}}+\frac{1}{2}\mu_{1n}^{2}{\bf{w}}_{\mathcal{B}}^{\top}\bar{{\bf{H}}}_{\mathcal{B},\mathcal{B}}{\bf{w}}_{\mathcal{B}}$
$\displaystyle\geq\frac{1}{2}\mu_{1n}^{2}M_{1}C^{2}_{1}-\mu_{1n}c_{1}\|{\bf{w}}_{\mathcal{B}}\|_{2}^{2}\sqrt{q\frac{\log
p}{n}}-\frac{1}{2}\mu_{1n}^{2}c_{2}\|{\bf{w}}_{\mathcal{B}}\|_{2}^{2}q\sqrt{\frac{\log
p}{n}}.$

For sufficiently large $n$, by assumption that $\rho_{1n}\sqrt{n/\log
p}\rightarrow\infty$ if $p\rightarrow\infty$ and $\sqrt{\log p/n}=o(1)$, the
second term in the last line above is
$O(\mu_{1n}\sqrt{q}\rho_{n})=o(\mu_{1n}^{2})$; the last term is
$o(\mu_{1n}^{2})$. Thus, for sufficiently large $n$, we have

$I_{1}\geq\frac{1}{2}\mu_{1n}^{2}M_{1}C^{2}_{1}.$ (A.16)

For the term $I_{2}$, by Cauchy-Schwartz and triangle inequality, we have

$|I_{2}|=\rho_{1n}|\|\bar{\mbox{\boldmath$\theta$}}^{o}+\mu_{1n}{\bf{w}}\|_{1}-\mu_{1n}\|{\bf{w}}\|_{1}|\leq
C_{1}\rho_{1n}\mu_{1n}\sqrt{q}=C_{1}\mu_{1n}^{2}.$ (A.17)

Next, we provide an upper bound for $I_{3}$. Note that

$\displaystyle I_{3}$
$\displaystyle:=\rho_{2n}\operatorname{trace}\left((\bar{\mbox{\boldmath$\Theta$}}+\mu_{1n}\mathcal{A}{\bf{w}})^{2}{\bf{S}}\right)$
$\displaystyle+\rho_{2n}\operatorname{trace}\left((\bar{\mbox{\boldmath$\Theta$}}+\mu_{1n}\mathcal{A}{\bf{w}})^{2}({\bf{N}}(\bar{\mbox{\boldmath$R$}}+\mu_{2n}\mathcal{A}{\bf{z}}){\bf{N}}^{\top})-\bar{\mbox{\boldmath$\Theta$}}^{2}\bar{{\bf{Q}}}\right)$
$\displaystyle=\rho_{2n}\operatorname{trace}\left(\bar{\mbox{\boldmath$\Theta$}}^{2}\bar{{\bf{S}}}\right)+\rho_{2n}\operatorname{trace}\left((\bar{\mbox{\boldmath$\Theta$}}+\mu_{1n}\mathcal{A}{\bf{w}})^{2}({\bf{S}}-\bar{{\bf{S}}})\right)$
$\displaystyle+\rho_{2n}\operatorname{trace}\left((\bar{\mbox{\boldmath$\Theta$}}+\mu_{1n}\mathcal{A}{\bf{w}})^{2}({\bf{N}}(\bar{\mbox{\boldmath$R$}}+\mu_{2n}\mathcal{A}{\bf{z}}){\bf{N}}^{\top})-\bar{\mbox{\boldmath$\Theta$}}^{2}\bar{{\bf{Q}}}\right)$

Now, we have

$\displaystyle\left|\operatorname{trace}\left((\bar{\mbox{\boldmath$\Theta$}}+\mu_{1n}\mathcal{A}{\bf{w}})^{2}({\bf{S}}-\bar{{\bf{S}}})\right)\right|$
$\displaystyle+\left|\operatorname{trace}\left((\bar{\mbox{\boldmath$\Theta$}}+\mu_{1n}\mathcal{A}{\bf{w}})^{2}({\bf{N}}(\bar{\mbox{\boldmath$R$}}+\mu_{2n}\mathcal{A}{\bf{z}}){\bf{N}}^{\top})-\bar{\mbox{\boldmath$\Theta$}}^{2}\bar{{\bf{Q}}}\right)\right|$
$\displaystyle=\left|\operatorname{trace}\left((\bar{\mbox{\boldmath$\Theta$}}+\mu_{1n}\mathcal{A}{\bf{w}})^{2}({\bf{N}}(\bar{\mbox{\boldmath$R$}}+\mu_{2n}\mathcal{A}{\bf{z}}){\bf{N}}^{\top})-\bar{\mbox{\boldmath$\Theta$}}^{2}\bar{{\bf{Q}}}\right)\right|+I_{3,0}$
$\displaystyle=\left|\operatorname{trace}\left((\bar{\mbox{\boldmath$\Theta$}}^{2}+2\mu_{1n}\bar{\mbox{\boldmath$\Theta$}}\mathcal{A}{\bf{w}})({\bf{N}}(\bar{\mbox{\boldmath$R$}}+\mu_{2n}\mathcal{A}{\bf{z}}){\bf{N}}^{\top})-\bar{\mbox{\boldmath$\Theta$}}^{2}\bar{{\bf{Q}}}\right)\right|+I_{3,1}+I_{3,0}$
$\displaystyle=\left|\operatorname{trace}\left(\bar{\mbox{\boldmath$\Theta$}}^{2}({\bf{N}}(\bar{\mbox{\boldmath$R$}}+\mu_{2n}\mathcal{A}{\bf{z}}){\bf{N}}^{\top}-\bar{{\bf{Q}}})\right)\right|+I_{3,2}+I_{3,1}+I_{3,0}$
$\displaystyle=I_{3,3}+I_{3,2}+I_{3,1}+I_{3,0}.$

From Assumption (A2), we have

$\displaystyle I_{3,0}$ $\displaystyle\leq O(\|{\bf{S}}-\bar{{\bf{S}}}\|),$
$\displaystyle I_{3,1}$
$\displaystyle=\mu_{1n}^{2}\left|\operatorname{trace}\left((\mathcal{A}{\bf{w}})^{2}{\bf{N}}\bar{\mbox{\boldmath$R$}}{\bf{N}}^{\top}\right)\right|+\mu_{1n}^{2}\mu_{2n}\left|\operatorname{trace}\left((\mathcal{A}{\bf{w}})^{2}{\bf{N}}\mathcal{A}{\bf{z}}{\bf{N}}^{\top}\right)\right|,$
$\displaystyle\leq\tau_{3}^{2}C_{1}^{2}\mu_{1n}^{2}+C_{1}^{2}D_{1}\mu_{1n}^{2}\mu_{2n},$
$\displaystyle I_{3,2}$
$\displaystyle=2\mu_{1n}\left|\operatorname{trace}\left(\bar{\mbox{\boldmath$\Theta$}}\mathcal{A}{\bf{w}}{\bf{N}}\bar{\mbox{\boldmath$R$}}{\bf{N}}^{\top}\right)\right|+2\mu_{1n}\mu_{2n}\left|\operatorname{trace}\left(\bar{\mbox{\boldmath$\Theta$}}\mathcal{A}{\bf{w}}{\bf{N}}\mathcal{A}{\bf{z}}{\bf{N}}^{\top}\right)\right|$
$\displaystyle\leq
2\tau_{2}\tau_{3}C_{1}\mu_{1n}+2\mu_{1n}\mu_{2n}\tau_{2}C_{1}D_{1}$
$\displaystyle I_{3,3}$
$\displaystyle:=\left|\operatorname{trace}\left(\bar{\mbox{\boldmath$\Theta$}}^{2}({\bf{N}}(\bar{\mbox{\boldmath$R$}}+\mu_{2n}\mathcal{A}{\bf{z}}){\bf{N}}^{\top}-\bar{{\bf{Q}}})\right)\right|\leq\tau_{2}^{2}C_{1}^{2}\mu_{2n}.$

Hence, for sufficiently large $n$, we have

$I_{3}\geq\tau_{1}\rho_{2n}-o(\rho_{2n}).$ (A.18)

Now, by combining (A.16)-(A.18), for sufficiently large $n$, we obtain

$\displaystyle
F_{n}(\bar{\mbox{\boldmath$\theta$}}^{d},\bar{\mbox{\boldmath$\theta$}}^{o}+\mu_{1n}{\bf{w}},\bar{{\bf{r}}}^{d},\bar{{\bf{r}}}^{o}+\mu_{1n}{\bf{z}};{\bf{Y}})-F_{n}(\bar{\mbox{\boldmath$\theta$}}^{d},\bar{\mbox{\boldmath$\theta$}}^{o},\bar{{\bf{r}}}^{d},\bar{{\bf{r}}}^{o};{\bf{Y}})$
$\displaystyle\geq\frac{1}{2}M_{1}\mu_{1n}^{2}C^{2}_{1}-\mu_{1n}^{2}C_{1}+\tau_{1}\rho_{2n}\geq
0.$

Here, the first inequality follows $\rho_{1n}=\sqrt{\log p/n}$ and the last
inequality follows by setting $C_{1}>2/M_{1}$.

Now, let
$\mathcal{S}_{{\bf{w}},{\bf{z}}}=\left\\{({\bf{w}},{\bf{z}}):{\bf{w}}_{\mathcal{B}^{c}}=0,\|{\bf{w}}\|=C_{1},\|{\bf{z}}\|=D_{1}\right\\}$.
Then, for $n$ sufficiently large, the following holds

$\inf_{({\bf{w}},{\bf{z}})\in\mathcal{S}_{{\bf{w}},{\bf{z}}}}~{}F_{n}(\bar{\mbox{\boldmath$\theta$}}^{d},\bar{\mbox{\boldmath$\theta$}}^{o}+\mu_{1n}{\bf{w}},\bar{{\bf{r}}}^{d},\bar{{\bf{r}}}^{o}+\mu_{2n}{\bf{z}};{\bf{Y}})>F_{n}(\bar{\mbox{\boldmath$\theta$}}^{d},\bar{\mbox{\boldmath$\theta$}}^{o},\bar{{\bf{r}}}^{d},\bar{{\bf{r}}}^{o};{\bf{Y}}),$

with probability at least $1-O(\exp(-\eta\log p))$.

Thus, any solution to the problem defined in (A.1) is within the disc (A.14)
with probability at least $1-O(\exp(-\eta\log p))$. Finally, since
$\bar{{\bf{Q}}}={\bf{N}}\bar{{\bf{R}}}{\bf{N}}^{\top}$ and
${\bf{N}}^{\top}{\bf{N}}={\bf{I}}_{(p-H+1)}$, we have that
$\|{\widehat{{\bf{q}}}}^{o}-\bar{{\bf{q}}}^{0}\|=\|{\widehat{{\bf{r}}}}^{o}-\bar{{\bf{r}}}^{0}\|$.
This completes the proof. ∎

######  Lemma 11.

Assume conditions of Lemma 12 hold and
$\rho_{2n}<\delta\rho_{1n}/(\tau_{2}\tau_{3})$. Then, there exists a constant
$C_{2}>0$, such that for any $\eta>0$, for sufficiently large $n$, the
following event holds with probability at least $1-O(\exp(-\eta\log p))$: for
any $\mbox{\boldmath$\theta$}^{o}$ satisfying

$\|\mbox{\boldmath$\theta$}^{o}-\bar{\mbox{\boldmath$\theta$}}^{o}\|\geq
C_{2}\sqrt{q}\rho_{1n},~{}~{}~{}\mbox{\boldmath$\theta$}^{o}_{\mathcal{B}^{c}}=0,$
(A.19)

we have
$\|\nabla_{\mbox{\boldmath$\theta$}^{o}}F_{n}(\bar{\mbox{\boldmath$\theta$}}^{d},\hat{\mbox{\boldmath$\theta$}}^{o}_{\mathcal{B}},\bar{{\bf{r}}}^{d},\hat{{\bf{r}}}^{o};{\bf{Y}})\|>\sqrt{q}\rho_{1n}$.

###### Proof.

The proof follows the idea of (Peng et al., 2009, Lemma S-4). For
$\mbox{\boldmath$\theta$}^{o}=\hat{\mbox{\boldmath$\theta$}}^{o}$ satisfying
(A.19), we have
$\hat{\mbox{\boldmath$\theta$}}^{o}=\bar{\mbox{\boldmath$\theta$}}^{o}+\mu_{1n}{\bf{w}}$,
with ${\bf{w}}_{\mathcal{B}^{c}}=0$ and $\|{\bf{w}}\|\geq C_{2}$. We have

$\displaystyle\nabla_{\mbox{\boldmath$\theta$}^{o}}F_{n}(\bar{\mbox{\boldmath$\theta$}}^{d},\hat{\mbox{\boldmath$\theta$}}^{o}_{\mathcal{B}},\bar{{\bf{r}}}^{d},\hat{{\bf{r}}}^{o};{\bf{Y}})$
$\displaystyle=(1+\rho_{2n})\hat{{\bf{g}}}^{n}_{\mathcal{B}}+\rho_{2n}({\bf{N}}\hat{{\bf{R}}}^{o}{\bf{N}}^{\top})\mathcal{A}^{*}\mathcal{A}\hat{\mbox{\boldmath$\theta$}}^{o}_{\mathcal{B}}$
$\displaystyle=(1+\rho_{2n})\bar{{\bf{g}}}^{n}_{\mathcal{B}}+(1+\rho_{2n})\mu_{1n}\bar{{\bf{H}}}_{\mathcal{B},\mathcal{B}}^{n}{\bf{w}}_{\mathcal{B}}+\rho_{2n}({\bf{N}}\hat{{\bf{R}}}^{o}{\bf{N}}^{\top})\mathcal{A}^{*}\mathcal{A}\bar{\mbox{\boldmath$\theta$}}^{o}_{\mathcal{B}}$
$\displaystyle\geq(1+\rho_{2n})\bar{{\bf{g}}}^{n}_{\mathcal{B}}+(1+\rho_{2n})\mu_{1n}\bar{{\bf{H}}}_{\mathcal{B},\mathcal{B}}^{n}{\bf{w}}_{\mathcal{B}}$
$\displaystyle+\rho_{2n}({\bf{N}}\hat{{\bf{R}}}^{o}{\bf{N}}^{\top})\mathcal{A}^{*}\mathcal{A}\bar{\mbox{\boldmath$\theta$}}^{o}_{\mathcal{B}}+\rho_{2n}\mu_{1n}({\bf{N}}\hat{{\bf{R}}}^{o}{\bf{N}}^{\top})\mathcal{A}^{*}\mathcal{A}{\bf{w}}_{\mathcal{B}},$

where the second inequality follows from Taylor expansion of
$\nabla_{\mbox{\boldmath$\theta$}^{o}}L_{n}(\bar{\mbox{\boldmath$\theta$}}^{d},\hat{\mbox{\boldmath$\theta$}}^{o}_{\mathcal{B}};{\bf{Y}})$.

Let

$\hat{{\bf{A}}}:={\bf{N}}\hat{{\bf{R}}}^{o}{\bf{N}}^{\top}\mathcal{A}^{*}\mathcal{A}^{\top}\bar{\mbox{\boldmath$\theta$}}^{o}~{}~{}\textnormal{and}~{}~{}\bar{{\bf{A}}}:={\bf{N}}\bar{{\bf{R}}}^{o}{\bf{N}}^{\top}\mathcal{A}^{*}\mathcal{A}^{\top}\bar{\mbox{\boldmath$\theta$}}^{o}.$
(A.20)

Then, we have

$\displaystyle\|\hat{{\bf{A}}}_{\mathcal{B}}\|$
$\displaystyle\leq\rho_{2n}\|{\bf{N}}\bar{{\bf{R}}}^{o}{\bf{N}}^{\top}\|\|\mathcal{A}^{*}\mathcal{A}\bar{\mbox{\boldmath$\theta$}}^{o}\|+\rho_{2n}\|\bar{{\bf{A}}}-\hat{{\bf{A}}}\|\leq
2\tau_{2}\tau_{3}\sqrt{q}\rho_{2n},$ (A.21)

where the last inequality follows since

$\displaystyle\|\bar{{\bf{A}}}\|\leq\Lambda_{\max}({\bf{N}}\bar{{\bf{R}}}^{o}{\bf{N}}^{\top})\|\mathcal{A}^{*}\mathcal{A}\bar{\mbox{\boldmath$\theta$}}^{o}\|\leq\tau_{2}\tau_{3}\sqrt{q},$
$\displaystyle\|\bar{{\bf{A}}}-\hat{{\bf{A}}}\|=o(\rho_{2n}).$

Now, let $\mu_{1n}=\sqrt{q}\rho_{1n}$. By triangle inequality and similar
proof strategies as in Lemma 10, for sufficiently large $n$, we obtain

$\displaystyle\|\nabla_{\mbox{\boldmath$\theta$}^{o}}F_{n}(\bar{\mbox{\boldmath$\theta$}}^{d},\hat{\mbox{\boldmath$\theta$}}^{o}_{\mathcal{B}},\bar{{\bf{r}}}^{d},\hat{{\bf{r}}}^{o};{\bf{Y}})\|$
$\displaystyle\geq(1+\rho_{2n})\mu_{1n}\|\bar{{\bf{H}}}_{\mathcal{B},\mathcal{B}}{\bf{w}}_{\mathcal{B}}\|-c_{0}(1+\rho_{2n})(q^{1/2}n^{-1/2}\sqrt{\log
p})$ $\displaystyle-
c_{3}(1+\rho_{2n})\|{\bf{w}}_{\mathcal{B}}\|_{2}(\mu_{1n}qn^{-1/2}\sqrt{\log
p})-2\tau_{2}\tau_{3}\sqrt{q}\rho_{2n}-o(\rho_{2n})$
$\displaystyle\geq(1+\rho_{2n})M_{1}C_{2}\sqrt{q}\rho_{1n}-2\tau_{2}\tau_{3}\sqrt{q}\rho_{2n},$

with probability at least $1-O(\exp(-\eta\log p))$. Here, the first inequality
uses Lemma 9 and the last inequality follows from Lemma 8 where we have
$\|\bar{{\bf{H}}}_{\mathcal{B},\mathcal{B}}{\bf{w}}_{\mathcal{B}}\|\geq
M_{1}\|{\bf{w}}_{\mathcal{B}}\|$. Now, taking

$C_{2}=\frac{1+2\delta}{(1+\rho_{2n})M_{1}+\epsilon}$ (A.22)

for some $\epsilon>0$, completes the proof. ∎

Next, inspired by Peng et al. (2009); Khare et al. (2015), we prove estimation
consistency for the nodewise FCONCORD, restricted to the true support, i.e.,
$\mbox{\boldmath$\theta$}^{o}_{\mathcal{B}^{c}}=0$.

######  Lemma 12.

Suppose Assumptions (A1)–(A4) are satisfied. Assume further that
$\rho_{1n}=O(\sqrt{\log p/n})$, $n>O(q\log(p))$ as $n\rightarrow\infty$,
$\rho_{2n}=O(\sqrt{\log(p-H+1)/n})$,
$\rho_{2n}\leq\delta\rho_{1n}/((1+M_{7}(\bar{\mbox{\boldmath$\theta$}}^{o}))\tau_{2}\tau_{3})$,
and $\mbox{\boldmath$\epsilon$}=\mbox{\boldmath$0$}$. Then, there exist finite
constants $C(\bar{\mbox{\boldmath$\theta$}}^{o})$ and $D(\bar{{\bf{q}}}^{o})$,
such that for any $\eta>0$, the following events hold with probability at
least $1-O(\exp(-\eta\log p))$:

  * •

There exists a local minimizer
$({\widehat{\mbox{\boldmath$\theta$}}}^{o}_{\mathcal{B}},{\widehat{{\bf{q}}}}^{o})$
of (A.1) such that

$\max\left(\|\widehat{\mbox{\boldmath$\theta$}}^{o}_{\mathcal{B}}-\bar{\mbox{\boldmath$\theta$}}^{o}_{\mathcal{B}}\|,\|{\widehat{{\bf{q}}}}^{o}-\bar{{\bf{q}}}^{o}\|\right)\\\
\leq\max\left(C(\bar{\mbox{\boldmath$\theta$}}^{o})\rho_{1n}\sqrt{q}/(1+\rho_{2n}),D(\bar{{\bf{q}}}^{o})\rho_{2n}\sqrt{\Psi(p,H,K)}\right),$

where $q$ and $\Psi(p,H,K)$ are defined in (23).

  * •

If $\min_{(i,j)\in\mathcal{B}}\bar{\theta}_{ij}\geq
2C(\bar{\mbox{\boldmath$\theta$}}^{o})\rho_{1n}\sqrt{q}/(1+\rho_{2n})$, then
${\widehat{\mbox{\boldmath$\theta$}}}^{o}_{\mathcal{B}^{c}}=0$.

###### Proof.

By the KKT condition, for any solution
$(\hat{\mbox{\boldmath$\theta$}}^{o},\hat{{\bf{r}}}^{o})$ of (A.1), it
satisfies

$\displaystyle\|\nabla_{\mbox{\boldmath$\theta$}^{o}}F_{n}(\bar{\mbox{\boldmath$\theta$}}^{d},\hat{\mbox{\boldmath$\theta$}}^{o}_{\mathcal{B}},\bar{{\bf{r}}}^{d},\hat{{\bf{r}}}^{o};{\bf{Y}})\|_{\infty}$
$\displaystyle\leq\ \rho_{1n}.$

Thus, for $n$ sufficiently large, we have

$\displaystyle\|\nabla_{\mbox{\boldmath$\theta$}^{o}}F_{n}(\bar{\mbox{\boldmath$\theta$}}^{d},\hat{\mbox{\boldmath$\theta$}}^{o}_{\mathcal{B}},\bar{{\bf{r}}}^{d},\hat{{\bf{r}}}^{o};{\bf{Y}})\|\leq\sqrt{q}\rho_{1n},$

Let $C(\bar{\mbox{\boldmath$\theta$}}^{o}):=C_{2}$. Using (A.22) and Lemma 11,
we obtain

$\displaystyle\|\hat{\mbox{\boldmath$\theta$}}^{o}-\bar{\mbox{\boldmath$\theta$}}^{o}\|_{2}\leq
C(\bar{\mbox{\boldmath$\theta$}}^{o})\sqrt{q}\rho_{1n}/(1+\rho_{2n}),~{}~{}~{}\mbox{\boldmath$\theta$}^{o}_{\mathcal{B}^{c}}=0$

with probability at least $1-O(\exp(-\eta\log p))$.

Now, if $\min_{(i,j)\in\mathcal{B}}\bar{\mbox{\boldmath$\theta$}}^{o}_{ij}\geq
2C(\bar{\mbox{\boldmath$\theta$}}^{o})\sqrt{q}\rho_{1n}/(1+\rho_{2n})$, then

$\displaystyle 1-O(\exp(-\eta\log p))$ $\displaystyle\leq
P_{\bar{\mbox{\boldmath$\theta$}}^{o}}\left(\|{\widehat{\mbox{\boldmath$\theta$}}}^{o}_{\mathcal{B}}-\bar{\mbox{\boldmath$\theta$}}^{o}_{\mathcal{B}}\|\leq
C(\bar{\mbox{\boldmath$\theta$}}^{o})\sqrt{q}\rho_{1n}/(1+\rho_{2n}),~{}~{}\min_{(i,j)\in\mathcal{B}}\bar{\mbox{\boldmath$\theta$}}^{o}_{ij}\geq
2C(\bar{\mbox{\boldmath$\theta$}}^{o})\sqrt{q}\rho_{1n}/(1+\rho_{2n})\right)$
$\displaystyle\leq
P_{\bar{\mbox{\boldmath$\theta$}}^{o}}\left(\text{sign}(\hat{\theta}^{o}_{ij})=\text{sign}(\bar{{\theta}}^{o}_{ij}),~{}~{}\forall(i,j)\in\mathcal{B}\right).$

∎

The following Lemma 13 shows that no wrong edge is selected with probability
tending to one.

######  Lemma 13.

Suppose that the conditions of Lemma 12 and Assumption (A5) are satisfied.
Suppose further that $p=O(n^{\kappa})$ for some $\kappa>0$. Then for $\eta>0$,
and for $n$ sufficiently large, the solution of (A.1) satisfies

$P\left(\|\nabla_{\mbox{\boldmath$\theta$}^{o}}F_{n}(\bar{\mbox{\boldmath$\theta$}}^{d},\hat{\mbox{\boldmath$\theta$}}^{o}_{{\mathcal{B}^{c}}},\bar{{\bf{r}}}^{d},\hat{{\bf{r}}}^{o};{\bf{Y}})\|_{\infty}<1\right)\geq
1-O(\exp(-\eta\log p)).$ (A.23)

###### Proof.

Let
$\mathcal{E}_{n,k}=\\{\text{sign}({\widehat{\mbox{\boldmath$\theta$}}}^{o}_{ij,\mathcal{B}})=\text{sign}(\bar{\mbox{\boldmath$\theta$}}^{o}_{ij,\mathcal{B}})\\}$.
Then by Lemma 12, $P_{\bar{\mbox{\boldmath$\theta$}}^{o}}(\mathcal{E}_{n})\geq
1-O(\exp(-\eta\log p))$ for large $n$. Define the sign vector $\hat{{\bf{t}}}$
for $\hat{\mbox{\boldmath$\theta$}}^{o}$ to satisfy the following properties,

$\begin{cases}\hat{t}_{ij}=\textnormal{sign}(\hat{\theta}^{o}_{ij}),&\textnormal{if
}\theta_{ij}^{o}\neq 0\ ,\\\ \rvert\hat{t}_{ij}\rvert\leq 1,&\textnormal{if
}\hat{\theta}_{ij}^{o}=0\ .\end{cases}$ (A.24)

for all $1\leq i<j\leq p$.

On $\mathcal{E}_{n,k}$, by the KKT condition and the expansion of $F_{n}$ at
$(\bar{\mbox{\boldmath$\theta$}}^{d},\hat{\mbox{\boldmath$\theta$}}^{o},\bar{{\bf{r}}}^{o},\hat{{\bf{r}}}^{o})$,
we have

$\displaystyle(1+\rho_{2n})\hat{{\bf{g}}}^{n}_{\mathcal{B}}+\rho_{1n}\hat{{\bf{t}}}_{\mathcal{B}}+\rho_{2n}({\bf{N}}\hat{{\bf{R}}}^{o}{\bf{N}}^{\top})\mathcal{A}^{*}\mathcal{A}\hat{\mbox{\boldmath$\theta$}}^{o}_{\mathcal{B}}=0.$
(A.25)

where
$\hat{{\bf{g}}}^{n}=\nabla_{\mbox{\boldmath$\theta$}^{o}}L_{n}(\bar{\mbox{\boldmath$\theta$}}^{d},\hat{\mbox{\boldmath$\theta$}}^{o};{\bf{Y}})$.
Then, we can write

$(1+\rho_{2n})(\hat{{\bf{g}}}^{n}_{\mathcal{B}}-\bar{{\bf{g}}}^{n}_{\mathcal{B}})=-\rho_{1n}\hat{{\bf{t}}}_{\mathcal{B}}-\rho_{2n}{\bf{N}}\hat{{\bf{R}}}^{o}{\bf{N}}^{\top}\mathcal{A}^{*}\mathcal{A}^{\top}\hat{\mbox{\boldmath$\theta$}}^{o}_{\mathcal{B}}-(1+\rho_{2n})\bar{{\bf{g}}}^{n}_{\mathcal{B}}.$

Let $\tilde{\mbox{\boldmath$\theta$}}^{o}$ denote a point in the line segment
connecting $\hat{\mbox{\boldmath$\theta$}}^{o}$ and
$\bar{\mbox{\boldmath$\theta$}}^{o}$. Applying the Taylor expansion, we obtain

$(1+\rho_{2n})\bar{{\bf{H}}}\left(\bar{\mbox{\boldmath$\theta$}}^{o}_{\mathcal{B}}-\hat{\mbox{\boldmath$\theta$}}^{o}_{\mathcal{B}}\right)=-(1+\rho_{2n})\bar{{\bf{g}}}^{n}_{\mathcal{B}}-\rho_{1n}\hat{{\bf{t}}}_{\mathcal{B}}+{\bf{L}}^{n}_{\mathcal{B}}-\rho_{2n}\hat{{\bf{A}}}_{\mathcal{B}}.$
(A.26)

where
${\bf{L}}^{n}:=\left(\bar{{\bf{H}}}^{n}-\bar{{\bf{H}}}\right)\left(\hat{\mbox{\boldmath$\theta$}}^{o}-\bar{\mbox{\boldmath$\theta$}}^{o}\right)$
and
$\hat{{\bf{A}}}:={\bf{N}}\hat{{\bf{R}}}^{o}{\bf{N}}^{\top}\mathcal{A}^{*}\mathcal{A}^{\top}\hat{\mbox{\boldmath$\theta$}}^{o}$.

Now, by utilizing the fact that
$\hat{\mbox{\boldmath$\theta$}}_{\mathcal{B}^{c}}^{o}=\bar{\mbox{\boldmath$\theta$}}_{\mathcal{B}^{c}}^{o}=0$,
we have

$\displaystyle(1+\rho_{2n})\bar{{\bf{H}}}_{\mathcal{B}^{c}\mathcal{B}}(\bar{\mbox{\boldmath$\theta$}}_{\mathcal{B}}^{o}-\hat{\mbox{\boldmath$\theta$}}_{\mathcal{B}}^{o})$
$\displaystyle=$
$\displaystyle-(1+\rho_{2n})\bar{{\bf{g}}}^{n}_{\mathcal{B}^{c}}-\rho_{1n}\hat{{\bf{t}}}_{\mathcal{B}^{c}}+{\bf{L}}^{n}_{\mathcal{B}^{c}}-\rho_{2n}\hat{{\bf{A}}}_{\mathcal{B}^{c}},$
(A.27)
$\displaystyle(1+\rho_{2n})\bar{{\bf{H}}}_{\mathcal{B}\mathcal{B}}(\bar{\mbox{\boldmath$\theta$}}_{\mathcal{B}}^{o}-\hat{\mbox{\boldmath$\theta$}}_{\mathcal{B}}^{o})$
$\displaystyle=$
$\displaystyle-(1+\rho_{2n})\bar{{\bf{g}}}^{n}_{\mathcal{B}}-\rho_{1n}\hat{{\bf{t}}}_{\mathcal{B}}+{\bf{L}}^{n}_{\mathcal{B}}-\rho_{2n}\hat{{\bf{A}}}_{\mathcal{B}}.$
(A.28)

Since $\bar{{\bf{H}}}^{n}_{\mathcal{B}\mathcal{B}}$ is invertible by
assumption, we get Now, using results from Lemmas 14 and 16, we obtain

$\displaystyle\rho_{1n}\|\hat{{\bf{t}}}_{\mathcal{B}^{c}}\|_{\infty}$
$\displaystyle=$
$\displaystyle\big{\|}\bar{{\bf{H}}}_{\mathcal{B}^{c}\mathcal{B}}(\bar{{\bf{H}}}_{\mathcal{B}\mathcal{B}})^{-1}(-(1+\rho_{2n})\bar{{\bf{g}}}^{n}_{\mathcal{B}}-\rho_{1n}\hat{{\bf{t}}}_{\mathcal{B}}$
(A.29) $\displaystyle+$