# Leveraging Standardization in Graph Learning

Thu Ha Phi, Alexandre Hippert-Ferrer, Florent Bouchard, Arnaud Breloy

► To cite this version:

Thu Ha Phi, Alexandre Hippert-Ferrer, Florent Bouchard, Arnaud Breloy. Leveraging Standardization in Graph Learning. 32nd European Signal Processing Conference (EUSIPCO), Aug 2024, Lyon, France. pp.2492-2496, 10.23919/EUSIPCO63174.2024.10715401. hal-04952046

HAL Id: hal-04952046

https://hal.science/hal-04952046v1

Submitted on 17 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Leveraging Standardization in Graph Learning

Thu Ha Phi

LEME, Univ. Paris Nanterre

phi.thuha@parisnanterre.fr

Alexandre Hippert-Ferrer

LASTIG, Univ. Gustave Eiffel

alexandre.hippert-ferrer@ensg.eu

Florent Bouchard

L2S, CNRS, Univ. Paris Saclay

florent.bouchard@centralesupelec.fr

Arnaud Breloy

CEDRIC, CNAM

arnaud.breloy@cnam.fr

**Abstract**—Standardization of the variables is a staple scaling tool that can be beneficial in most statistical learning tasks. In the context of graph learning, we propose to take further advantage of this pre-processing, as it yields an additional implicit normalization structure for the covariance matrix. For Gaussian graphical models, we propose a new algorithm that aims at estimating a sparse precision matrix while its inverse is constrained to be a (possibly low-rank structured) correlation matrix. The corresponding optimization problem is addressed by using the framework of Riemannian optimization. Simulations over various graph structures illustrate the interest of the proposed approach.

**Index Terms**—graph learning, Gaussian graphical models, low-rank factorization, Riemannian optimization.

## I. INTRODUCTION

Graphs are ubiquitous structures in signal processing and machine learning. For example, the knowledge of the underlying graph topology behind the data is a key component of data processing methodologies in graph signal processing [1], optimal transport [2], and graph neural networks [3]. When it is unknown, graph learning aims at estimating the graph topology from the data itself. This recovered structure can then be used as a compact representation of relationship between entities in data exploratory analysis (visualization, clustering), or to apply one of the aforementioned techniques.

In a statistical formalization of graph learning, the graph nodes represent variables $x_q$, $q \in [1, p]$ of a multivariate observation $\mathbf{x} = [x_1, \dots, x_p] \in \mathbb{R}^p$, and an edge between two nodes $x_q$ and $x_\ell$ represents a conditional correlation between these two variables, i.e., the that the coefficient

$$ \operatorname{corr}[x_q x_\ell | \mathbf{x} \in [1, p] \setminus \{q, \ell\}] = \frac{\Theta_{q\ell}}{\sqrt{\Theta_{qq} \Theta_{\ell\ell}}} \quad (1) $$

is not null (here, $\Theta$ denotes the precision matrix, i.e., $\Theta = \Sigma^{-1}$ for the covariance matrix $\Sigma = \mathbb{E}[\mathbf{x}\mathbf{x}^\top]$). Under the assumption that the data set $\{\mathbf{x}_i\}_{i=0}^n$ follows a centered Gaussian distribution $\mathbf{x} \sim \mathcal{N}(0, \Sigma)$, also referred to as Gaussian Graphical Model (GGM), graph learning is usually performed by minimizing a penalized log-likelihood, i.e., solving:

$$ \underset{\Theta \in S_p^{+}}{\operatorname{minimize}} \quad \operatorname{tr}(\mathbf{S}\Theta) - \log |\Theta| + \lambda h(\Theta) \quad (2) $$

where $\mathbf{S} = \frac{1}{n} \sum_{i=1}^n \mathbf{x}_i \mathbf{x}_i^\top$ is the sample covariance matrix, $\lambda$ is a sparsity inducing penalty, and $\lambda$ is a regularization parameter. The support of the solution $\Theta^*$ of (2) is then used to infer the graph topology. For example, GLasso [4],

[5] is a prominent algorithm to tackle (2) when $h$ is an $\ell_1$-norm penalty. Many other works considered leveraging some additional structure constraints on $\Theta$ or $\Sigma$ in (2). Major examples include factor models, where $\Sigma$ is constrained to be low-rank plus diagonal [6], [7], and Laplacian constrained models, where $\Theta$ is constrained to be a Laplacian matrix [8], [9]. This work will follow a similar approach and proposes a new class of graph learning algorithms that constrains the resulting estimate of $\Sigma$ to be normalized as correlation matrix.

Indeed, most of data processing chains are actually applied on data that is implicitly standardized (i.e., each variable is scaled to unit variance). This normalization is beneficial in terms of graph edge detection performance, as it mitigates issues due to the variance discrepancies between all variables. This normalization also preserves the support of the conditional correlation matrix (i.e., the graph topology). Moreover, when $n$ is large enough the standardization does not affect the assumption of Gaussianity. Hence, graph learning algorithms inspired by (2) can directly be applied by replacing $\mathbf{S}$ with the sample correlation matrix. However, this type of plug-in approach does not modify the assumption behind the GGM, which should additionally assume that the standardized data follows a Gaussian distribution where the covariance matrix is a correlation matrix. Given these considerations, this paper presents the following contributions:

* To assess for standardization in a graph learning process, we propose a correlation constrained GGM (CC-GGM), i.e., we consider (2) with the additional constraint that $\Sigma = \Theta^{-1}$ is structured as a correlation matrix. For this problem we consider a generalized parameterization that can additionally account for low-rank structured correlation matrices [10], [11].

* We derive an algorithm to solve the corresponding graph learning problem by using Riemannian optimization [12], [13] on the oblique manifold [14].

* We finally demonstrate the interest of the proposed approach on simulated data for several random graph models. As expected, standardizing the data improves the performance of the graph learning process when the original data exhibits variance discrepancies. Furthermore, the proposed algorithm outperforms GGM (2) applied to standardized data, i.e., there is indeed an interest in additionally constraining the structure of the estimate to be in accordance with the standardization. The low-rank counterpart of the proposed algorithm is however not as efficient as the full rank one: it induces a performance versus computational complexity trade-off that requires careful tuning.

This work was supported by the MASSILIA project (ANR- 21-CE23- 0038-01) of the French National Research Agency (ANR).
## II. CORRELATION CONSTRAINED GGM (CC-GGM)

We handle a standardized data set denoted {$\tilde{x}_i$}$_{i=1}^n$, and the corresponding sample correlation matrix is denoted $\tilde{S} = \frac{1}{n} \sum_{i=1}^n \tilde{x}_i \tilde{x}_i^\top$. Given the assumption that the initial data follows a GGM $\mathbf{x} \sim N(0, \mathbf{\Sigma})$, the CC-GGM states that $\tilde{\mathbf{x}} \sim (0, \mathbf{C})$, where $\mathbf{C} = \text{diag}(\mathbf{\Sigma})^{-1/2}\mathbf{\Sigma}\text{diag}(\mathbf{\Sigma})^{-1/2}$ is the correlation matrix related to the initial covariance matrix $\mathbf{\Sigma}$. Transposing from (2), leaning the graph in a CC-GGM can be done by solving the problem

$$ \underset{\mathbf{C} \in \mathcal{C}_p^+}{\text{minimize}} \quad \text{tr}(\tilde{\mathbf{S}}\mathbf{C}^{-1}) + \log |\mathbf{C}| + \lambda h(\mathbf{C}^{-1}) \quad (3) $$

where $\mathcal{C}_p^+$ denotes the set of $p \times p$ correlation matrices (positive definite matrices with unit diagonal elements), and where $h$ is a penalty that promotes sparsity in $\mathbf{C}^{-1}$. To ensure the constraint $\mathbf{C} \in \mathcal{C}_p^+$ in problem (3), we consider the following parameterization

$$ \mathbf{C}_{\mathbf{W}}^{\beta} = (1 - \beta)\mathbf{W}\mathbf{W}^{\top} + \beta\mathbf{I}_p \quad (4) $$

with $\beta \in [0, 1]$ and where $\mathbf{W}$ is a matrix in the oblique manifold $\mathcal{OB}_{p,k}$ [14], defined as

$$ \mathcal{OB}_{p,k} = \{ \mathbf{W} \in \mathbb{R}^{p \times k} : \text{ddiag}(\mathbf{W}\mathbf{W}^{\top}) = \mathbf{I}_p \}. \quad (5) $$

The constraint $\mathbf{W} \in \mathcal{OB}_{p,k}$ directly implies that $\mathbf{C}_{\mathbf{W}}^{\beta} \in \mathcal{C}_p^+$. If we set $k = p$ and $\beta = 0$, $\mathbf{C}_{\mathbf{W}}^{\beta}$ yields a parameterization of a full rank correlation matrix that can be used to addressed the problem (3). In case we set $k < p$, the matrix $\mathbf{C}_{\mathbf{W}}^{\beta}$ is further constrained, as it exhibits a low-rank structure. This factorization will thus be advantageous in order to derive algorithms that scale to large dimension $p$: in essence the computational bottleneck will be of complexity $O(pk^2)$, as opposed to $O(p^3)$ in the full rank case [10], [11]. If $k < p$, we rely on a shrinkage to the identity ($\beta > 0$) to ensure that this matrix is always invertible. Notice that this shrinkage limits the maximal correlation value in $\mathbf{C}_{\mathbf{W}}^{\beta}$ to be at most $1 - \beta$, so we use a small value such as $\beta = 0.05$ in practice.

Merging (3) and (4), we cast the CC-GGM learning problem as follows:

$$ \underset{\mathbf{W} \in \mathcal{OB}_{p,k}}{\text{minimize}} \quad \text{tr}(\tilde{\mathbf{S}}(\mathbf{C}_{\mathbf{W}}^{\beta})^{-1}) + \log |\mathbf{C}_{\mathbf{W}}^{\beta}| + \lambda h((\mathbf{C}_{\mathbf{W}}^{\beta})^{-1}). \quad (6) $$

It remains to explicitly set a penalty $h$ that promotes sparsity in $\mathbf{C}^{-1}$. In this work, we use a smooth surrogate of $\ell_1$-norm which was shown to be a valid option in graph learning problems in [6]. We therefore consider

$$ h(\mathbf{C}^{-1}) = \sum_{q \neq \ell} \phi([\mathbf{C}^{-1}]_{q\ell}) \quad (7) $$

with

$$ \phi(t) = \epsilon \log(\cosh(t/\epsilon)) \quad (8) $$

and $\epsilon > 0$. We obtain the $\ell_1$-norm as $\lim_{\epsilon \to 0} \phi(t)$, and we use $\epsilon = 10^{-6}$ in practice. As in [6], experiments evidenced that the obtained results are not highly sensitive to this parameter once it is set to a sufficiently small value.

<figure><img src="image_2.png" /><figcaption>Fig. 1: Illustration of needed geometrical objects for the first order Riemannian optimization framework on $\mathcal{OB}_{p,k}$.</figcaption></figure>

## III. CC-GGM ALGORITHM

In this section, an algorithm is derived to efficiently solve the optimization problem (6). It relies on a Riemannian gradient descent [12] on the oblique manifold $\mathcal{OB}_{p,k}$ defined in (5). To achieve this, we need two things: (i) a proper Riemannian optimization framework on $\mathcal{OB}_{p,k}$; and (ii) the gradient of the cost function corresponding to CC-GGM. These are provided in Subsections III-A and III-B, respectively.

### A. Riemannian Optimization on $\mathcal{OB}_{p,k}$

This subsection provides a generic Riemannian optimization framework on the oblique manifold $\mathcal{OB}_{p,k}$ defined in (5). The first order Riemannian optimization workflow of a cost function $f$ is quite straightforward: given some iterate $\mathbf{W}$, a descent direction of $f$ - i.e., a tangent vector $\xi$ such that the directional derivative $\mathrm{D}f(\mathbf{W})[\xi] < 0$ - is obtained from the Riemannian gradient of $f$ at $\mathbf{W}$. Then, a new point is achieved on the manifold thanks to a retraction at $\mathbf{W}$, which is a mapping from the tangent space back onto the manifold. For full details on Riemannian optimization, the reader is referred to [12], [13]. Here, the focus is on $\mathcal{OB}_{p,k}$, for which we employ the tools proposed in [14]. The geometrical objects that are presented here are illustrated in Figure 1.

To be able to obtain a descent direction of a cost function in $\mathcal{OB}_{p,k}$, the first task is to define the tangent space $T_{\mathbf{W}}\mathcal{OB}_{p,k}$ at $\mathbf{W} \in \mathcal{OB}_{p,k}$. It is given by

$$ T_{\mathbf{W}}\mathcal{OB}_{p,k} = \{\xi \in \mathbb{R}^{p \times k} : \text{ddiag}(\mathbf{W}\xi^{\top}) = \mathbf{0}\}. \quad (9) $$

One then needs to choose a Riemannian metric - i.e., smoothly varying inner products on tangent spaces. With the oblique constraint, several metric have been considered; see e.g., [14], [15]. In this work, as in [14], we select the easiest solution: the Euclidean metric defined, for $\mathbf{W} \in \mathcal{OB}_{p,k}$, $\xi$ and $\eta \in T_{\mathbf{W}}\mathcal{OB}_{p,k}$, by

$$ \langle \xi, \eta \rangle_{\mathbf{W}} = \text{tr}(\xi^{\top} \eta). \quad (10) $$

From there, one can define the orthogonal projection $P_{\mathbf{W}}$ - according to (10) - from the ambient space $\mathbb{R}^{p \times k}$ onto the tangent space $T_{\mathbf{W}}\mathcal{OB}_{p,k}$. Given $\mathbf{Z} \in \mathbb{R}^{p \times k}$, it is

$$ P_{\mathbf{W}}(\mathbf{Z}) = \mathbf{Z} - \text{ddiag}(\mathbf{Z}\mathbf{W}^{\top})\mathbf{W}. \quad (11) $$

It follows that the Riemannian gradient $\nabla f(\mathbf{W})$ of some cost function $f: \mathcal{OB}_{p,k} \to \mathbb{R}$ at $\mathbf{W} \in \mathcal{OB}_{p,k}$ - which is defined as the only tangent vector such that, for all $\xi \in T_{\mathbf{W}}\mathcal{OB}_{p,k}$,
**Algorithm 1** – CC-GGM algorithm

**Input:** data $\mathbf{S} \in S_p^{++}$, initial guess $\mathbf{W}_0 \in \mathcal{O}B_{p,k}$
Set $t = 0$

**repeat**

Compute $\varphi_{\beta}(\mathbf{W}_t)$ with (18)

Compute $\nabla^{\mathcal{E}}\mathcal{L}(\varphi_{\beta}(\mathbf{W}_t))$ with (20)

Compute $\nabla^{\mathcal{E}}f(\mathbf{W}_t)$ with Proposition 1

Compute $\nabla f(\mathbf{W}_t)$ with projection (12)

$\mathbf{W}_{t+1} = R_{\mathbf{W}_t}(-\alpha_t \nabla f(\mathbf{W}_t))$ with retraction (13)

Set $t = t + 1$

**until** convergence

**Return:** estimated solution $\mathbf{W}^* = \mathbf{W}_t$

$\mathcal{D}f(\mathbf{W})[\xi] = \langle\nabla f(\mathbf{W}), \xi\rangle_{\mathbf{W}}$ – can be obtained from the Euclidean gradient $\nabla^{\mathcal{E}}f(\mathbf{W})$ of $f$ at $\mathbf{W}$ through

$$\nabla f(\mathbf{W}) = P_{\mathbf{W}}(\nabla^{\mathcal{E}} f(\mathbf{W})). \quad (12)$$

The Riemannian gradient $\nabla f(\mathbf{W})$ allows to define a descent direction $\xi$. For instance, for the gradient descent, the descent direction is simply $\xi = -\nabla f(\mathbf{W})$. To achieve a new point on the manifold from the descent direction $\xi$ at $\mathbf{W}$, we resort to a retraction on $\mathcal{O}B_{p,k}$. At $\mathbf{W} \in \mathcal{O}B_{p,k}$, it is a mapping $R_{\mathbf{W}}: T_{\mathbf{W}}\mathcal{O}B_{p,k} \to \mathcal{O}B_{p,k}$ such that $R_{\mathbf{W}}(\xi) = \mathbf{W} + \xi + o(\|\xi\|)$; see e.g., [12] for the full details. In our case, we can simply choose a projection based retraction [16], which is in our case

$$R_{\mathbf{W}}(\xi) = \mathcal{P}(\mathbf{W} + \xi), \quad (13)$$

where $\mathcal{P}: \mathbb{R}^{p \times k} \to \mathcal{O}B_{p,k}$ is the projection such that, for $\mathbf{Z} \in \mathbb{R}^{p \times k}$,

$$\mathcal{P}(\mathbf{Z}) = \text{ddiag}(\mathbf{Z}\mathbf{Z}^{\top})^{-1/2}\mathbf{Z}. \quad (14)$$

All necessary ingredients to perform Riemannian gradient descent on $\mathcal{O}B_{p,k}$ have now been presented. Indeed, to optimize $f: \mathcal{O}B_{p,k} \to \mathbb{R}$, given some initial guess $\mathbf{W}_0$, the sequence of iterates $\{\mathbf{W}_t\}$ is recursively obtained through

$$\mathbf{W}_{t+1} = R_{\mathbf{W}_t}(-\alpha_t \nabla f(\mathbf{W}_t)), \quad (15)$$

where $\alpha_t$ is the stepsize that can be computed through a linesearch.

## B. Euclidean gradient for CC-GGM

With the tools of the previous subsection, to be able to actually solve (6), it only remains to get the Euclidean gradient of the corresponding cost function. The Riemannian one is then obtained with (12). The cost function of interest in problem (6) is $f: \mathcal{O}B_{p,k} \to \mathbb{R}$ can be expressed as

$$f(\mathbf{W}) = \mathcal{L} \circ \varphi_{\beta}(\mathbf{W}), \quad (16)$$

where $\mathcal{L}: S_p^{++} \to \mathbb{R}$ is defined as

$$\mathcal{L}(\mathbf{C}) = \text{tr}(\mathbf{S}\mathbf{C}^{-1}) + \log|\mathbf{C}| + \lambda h(\mathbf{C}^{-1}), \quad (17)$$

and $\varphi_{\beta}: \mathcal{O}B_{p,k} \to S_p^{++}$ is such that

$$\varphi_{\beta}(\mathbf{W}) = \mathbf{C}_{\mathbf{W}}^{\beta} = (1 - \beta)\mathbf{W}\mathbf{W}^{\top} + \beta\mathbf{I}_{p}. \quad (18)$$

First, in Proposition 1, the Euclidean gradient of $f$ at $\mathbf{W}$ is provided as a function of the Euclidean gradient of $\mathcal{L}$ at $\varphi_{\beta}(\mathbf{W})$.

**Proposition 1.** The Euclidean gradient $\nabla^{\mathcal{E}}f(\mathbf{W})$ of $f: \mathcal{O}B_{p,k} \to \mathbb{R}$ defined in (16) at $\mathbf{W}$ is given by

$$\nabla^{\mathcal{E}}f(\mathbf{W}) = 2(1 - \beta)\nabla^{\mathcal{E}}\mathcal{L}(\varphi_{\beta}(\mathbf{W}))\mathbf{W}, \quad (19)$$

where $\nabla^{\mathcal{E}}\mathcal{L}(\varphi_{\beta}(\mathbf{W}))$ is the Euclidean gradient of $\mathcal{L}: S_p^{++} \to \mathbb{R}$ defined in (17) at $\varphi_{\beta}(\mathbf{W})$.

*Proof.* Deriving $f(\mathbf{W}) = \mathcal{L}(\varphi_{\beta}(\mathbf{W}))$ in direction $\xi$ yields

$$\mathcal{D}f(\mathbf{W})[\xi] = \mathcal{D}\mathcal{L}(\varphi_{\beta}(\mathbf{W}))\left[\mathcal{D}\varphi_{\beta}(\mathbf{W})[\xi]\right],$$

where $\mathcal{D}\varphi_{\beta}(\mathbf{W})[\xi] = (1 - \beta)(\mathbf{W}\xi^{\top} + \xi\mathbf{W}^{\top})$. By definition of the gradient, we have

$$\text{tr}(\nabla^{\mathcal{E}}f(\mathbf{W})^{\top}\xi) = \text{tr}(\nabla^{\mathcal{E}}\mathcal{L}(\varphi_{\beta}(\mathbf{W}))\mathcal{D}\varphi_{\beta}(\mathbf{W})[\xi]).$$

Injecting the formula of $\mathcal{D}\varphi_{\beta}(\mathbf{W})[\xi]$ and basic manipulation of the trace yield

$$\text{tr}(\nabla^{\mathcal{E}}f(\mathbf{W})^{\top}\xi) = 2(1 - \beta) \text{tr}((\nabla^{\mathcal{E}}\mathcal{L}(\varphi_{\beta}(\mathbf{W}))\mathbf{W})^{\top}\xi).$$

The result is then obtained by identification. $\square$

Finally, the Euclidean gradient $\nabla^{\mathcal{E}}\mathcal{L}(\mathbf{C})$ of $\mathcal{L}$ at $\mathbf{C} \in S_p^{++}$, which can be found in [6], is

$$\nabla^{\mathcal{E}}\mathcal{L}(\mathbf{C}) = \mathbf{C}^{-1}[\mathbf{C} - \mathbf{S} - \lambda\nabla h(\mathbf{C}^{-1})]\mathbf{C}^{-1}, \quad (20)$$

where $\nabla h(\mathbf{C}^{-1}) \in \mathbb{R}^{p \times p}$, whose elements are $[\nabla h(\mathbf{C}^{-1})]_{q\ell} = \phi'(\mathbf{C}_{q\ell}^{-1}) = \tanh(\mathbf{C}_{q\ell}^{-1}/\epsilon)$. Our proposed method to solve (3) is summarized in Algorithm 1.

Once the solution $\mathbf{W}^*$ is obtained from Algorithm 1, we have $\mathbf{C}^* = (1 - \beta)\mathbf{W}^*(\mathbf{W}^*)^{\top} + \beta\mathbf{I}_p$, and the support of the solution $(\mathbf{C}^*)^{-1}$ can then be used to infer the graph topology.

## IV. EXPERIMENTS

In this section, we illustrate the interest of leveraging the data standardization and its corresponding correlation structure in graph learning problems. Experiments are conducted on simulated data following a GGM with various underlying graph structures. The main point is to compare the performance of the GGM method (2) applied on raw data (whose solution is computed with algorithm from [6]), standardized data (referred to as GGM-st), and the proposed structure constrained method for problem (3), i.e., CC-GGM computed with Algorithm 1 using $k = p$ and $\beta = 0$. Additionally, we compare the performance of CC-GGM algorithm using various ranks $k < p$ and $\beta = 0.05$ to assess for the effect of rank reduction on the CC-GGM model.

### A. Simulation setup

**Graph generation:** The number of nodes for every graph (data dimension) is set to $p = 50$. The graph typologies are sampled from three random graph models: Barabási-Albert (BA), Watts-Strogatz (WS, with probability of edge creation is 0.1, and 5 nearest neighbor nodes are associated to each node), and Erdős–Rényi (ER, with probability of edge
<figure><img src="image_3.png" /><figcaption>Fig. 2: Examples of sampled graphs using (from left to right) the BA, WS, and ER models. The first row displays the graph and the second row is its respective adjacency matrix.</figcaption></figure>

creation is 0.1). For all of these graph models, weights of the edges are sampled from $U(2, 5)$ as in [6], [9]. Each sampling yields a weighted adjacency matrix $\mathbf{A}$ such that $[\mathbf{A}]_{ii} = 0$, $[\mathbf{A}]_{ij(i \neq j)} \neq 0$ if the edge exists between two variables and $[\mathbf{A}]_{ij(i \neq j)} = 0$ otherwise. Figure 2 displays visual examples of sampled graphs for all considered graph models.

**Data generation:** For each sampling of a Monte-Carlo experiment we generate a graph adjacency $\mathbf{A}$ from the aforementioned models. Then, the data set $\{\mathbf{x}_i\}_{i=1}^n$ is drawn from a GGM model $\mathbf{x} \sim N(0, \hat{\mathbf{\Sigma}})$. The covariance matrix $\hat{\mathbf{\Sigma}}$ is constructed in accordance with the graph structure and so that the variance discrepancies between the variables is controlled. To do so, we first construct an intermediary precision matrix $\mathbf{\Theta} = \mathbf{L} + \kappa \mathbf{I}_p$ ($\kappa$ is set to $1e^{-1}$ so that $\mathbf{\Theta}$ is non-singular), where $\mathbf{L}$ is the graph Laplacian matrix $\mathbf{L} = \mathbf{D} - \mathbf{A}$, and where $\mathbf{D}$ is the degree matrix, i.e., $[\mathbf{D}]_{ii} = \sum_{i \neq j} [\mathbf{A}]_{ij}$ and $[\mathbf{D}]_{ij(i \neq j)} = 0$. We then compute $\mathbf{\Sigma} = \mathbf{\Theta}^{-1}$ and $\mathbf{C} = \text{diag}(\mathbf{\Sigma})^{-1/2} \mathbf{\Sigma} \text{diag}(\mathbf{\Sigma})^{-1/2}$. Finally, the covariance matrix is set as $\hat{\mathbf{\Sigma}} = \text{Diag}(\sigma) \text{CDiag}(\sigma)^{\top}$, where $\sigma \in \mathbb{R}^p$ is a random vector of variance coefficients, in which each entry is sampled uniformly in the range $[1, \sigma_{\max}]$. In this setup, we set $\sigma_{\max} = 10$, i.e., we can observe up to 20dB of power difference between the variables.

## B. Algorithm comparison

**Evaluation metrics:** Each method yields either an estimate of the covariance or precision matrix. From these estimates, we can always compute the corresponding conditional correlation matrix from (1). The resulting estimate of the graph topology is then obtained by applying a detection threshold in $[0, 1]$ to this quantity. For each methods, we then investigate the receiver operating characteristic (ROC) curves, that displays the false edge discovery rate (false positive) versus the true edge detection rate (true positive) as a function of the this threshold. The area under curve (AUC) of the ROC is also

used as a global metric to compare the methods.

**Fair tuning of $\lambda$:** We select the regularization parameter $\lambda$ that yields the highest AUC using a grid search over 100 Monte-Carlo trials. The displayed result corresponds to this value that is set independently for each setup and method.

## C. Results

Figure 3 displays the ROC of the compared graph learning algorithms for various setups of underlying graphs BA, WS, ER) and for two scenarios corresponding to a small sample support ($n = 45$, so $n \approx p$ but $n < p$) and larger sample support ($n = 105$, $n \approx 2p$). We first recall that the variance discrepancies between the variables are piloted by the parameter $\sigma_{\max}$. Notice that if $\sigma_{\max} = 1$ the covariance matrix of the data is actually a correlation matrix, and the standardization has little effect to the data. In that case, algorithms GGM, GGM-st, and CC-GGM perform identically, as it can be expected. Therefore, this setup is omitted in the displays. In the considered case where $\sigma_{\max} = 10$, we observe that standardizing the data is beneficial to the graph learning process, as GGM-st outperforms GGM in terms of ROC and AUC in all setups. Furthermore, the proposed method CC-GGM (with $k = p$ and $\beta = 0$) slightly outperforms GGM-st, which means that there is indeed an interest in additionally constraining the structure of the estimate to be in accordance with the standardization, i.e., solving (3), rather than (2), applied to $\tilde{\mathbf{S}}$. The low-rank counterpart of the proposed algorithm is however not as efficient as the full rank one: it induces a performance loss compared to CC-GGM and GGM-st. This can probably be explained by the fact that the data is sampled from a covariance matrix that is not constructed to naturally exhibit a low-rank structure (meaning that the estimate is restricted to hold a structure that is not present in the data). Still, in some setups, this rank-reduced method can achieve reasonable AUCs (up to 0.8 ~ 0.92) even for fixed $k \ll p$. Hence it still provides a performance versus computational trade-off that can be convenient in the large $p$ regime, in which a computational complexity of $O(p^3)$ (the one of GGM, and CC-GGM) is an issue.

## V. CONCLUSION

This paper investigated the effect of standardization in a graph learning process. It was already well known that this pre-processing could be beneficial in terms of edge detection performance. However, its additional normalization structure induced on the covariance matrix of the data was not yet leveraged in GGMs. We proposed CC-GGM: a new algorithm to estimate a sparse precision matrix while its inverse is constrained to be correlation matrix. We also considered the additional option for this matrix to be low-rank structured. The resulting optimization problems were tackled with Riemannian optimization on the oblique manifold. Simulations illustrated the interest of the approach under various setups on multiple graph structures.
<figure><img src="image_4.png" /><figcaption>Fig. 3: ROC/AUC of the compared graph learning algorithms on BA, WS, and ER graph models (from left to right) for $n = 105$ (top) and $n = 45$ (bottom). For every experiment, $p = 50$ and number of trials = 100, and $\sigma_{max} = 10$.</figcaption></figure>

REFERENCES

[1] Antonio Ortega, Pascal Frossard, Jelena Kovačević, José MF Moura, and Pierre Vandergheynst, “Graph signal processing: Overview, challenges, and applications,” *Proceedings of the IEEE*, vol. 106, no. 5, pp. 808–828, 2018.

[2] Hermina Petric Maretic, Mireille El Gheche, Giovanni Chierchia, and Pascal Frossard, “Got: an optimal transport framework for graph comparison,” *Advances in Neural Information Processing Systems*, vol. 32, 2019.

[3] Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and S Yu Philip, “A comprehensive survey on graph neural networks,” *IEEE transactions on neural networks and learning systems*, vol. 32, no. 1, pp. 4–24, 2020.

[4] Jerome Friedman, Trevor Hastie, and Robert Tibshirani, “Sparse inverse covariance estimation with the graphical lasso,” *Biostatistics*, vol. 9, no. 3, pp. 432–441, 2008.

[5] Rahul Mazumder and Trevor Hastie, “The graphical lasso: New insights and alternatives,” *Electronic journal of statistics*, vol. 6, pp. 2125, 2012.

[6] Alexandre Hippert-Ferrer, Florent Bouchard, Ammar Mian, Titouan Vayer, and Arnaud Breloy, “Learning graphical factor models with riemannian optimization,” in *Joint European Conference on Machine Learning and Knowledge Discovery in Databases*. Springer, 2023, pp. 349–366.

[7] Noirrit Kiran Chandra, Peter Mueller, and Abhra Sarkar, “Bayesian scalable precision factor analysis for massive sparse gaussian graphical models,” *arXiv preprint arXiv:2107.11316*, 2021.

[8] Hilmi E Egilmez, Eduardo Pavez, and Antonio Ortega, “Graph learning

from data under laplacian and structural constraints,” *IEEE Journal of Selected Topics in Signal Processing*, vol. 11, no. 6, pp. 825–841, 2017.

[9] Sandeep Kumar, Jiaxi Ying, José Vinícius de M Cardoso, and Daniel P Palomar, “A unified framework for structured graph learning via spectral constraints,” *Journal of Machine Learning Research*, vol. 21, no. 22, pp. 1–60, 2020.

[10] Raoul Pietersz 4 and Patrick JF Groenen, “Rank reduction of correlation matrices by majorization,” *Quantitative Finance*, vol. 4, no. 6, pp. 649–662, 2004.

[11] Thu Ha Phi, Alexandre Hippert-Ferrer, Florent Bouchard, and Arnaud Breloy, “Robust low-rank correlation fitting,” in *ICASSP 2024 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*. IEEE, 2024.

[12] P-A Absil, Robert Mahony, and Rodolphe Sepulchre, *Optimization algorithms on matrix manifolds*, Princeton University Press, 2008.

[13] Nicolas Boumal, *An introduction to optimization on smooth manifolds*, Cambridge University Press, 2023.

[14] P-A Absil and Kyle A Gallivan, “Joint diagonalization on the oblique manifold for independent component analysis,” in *2006 IEEE international conference on acoustics, speech and signal processing proceedings*. IEEE, 2006, vol. 5, pp. V–V.

[15] Florent Bouchard, Bijan Afsari, Jérôme Malick, and Marco Congedo, “Approximate joint diagonalization with riemannian optimization on the general linear group,” *SIAM Journal on Matrix Analysis and Applications*, vol. 41, no. 1, pp. 152–170, 2020.

[16] P-A Absil and Jérôme Malick, “Projection-like retractions on matrix manifolds,” *SIAM Journal on Optimization*, vol. 22, no. 1, pp. 135–158, 2012.