Constructing networks by filtering correlation matrices: A null model approach

TL;DR

This paper introduces a null model-based regularization method for constructing networks from correlation matrices, effectively filtering out spurious correlations and improving predictive accuracy over traditional thresholding methods.

Contribution

It proposes a flexible algorithm that integrates null models with regularization for network construction, allowing model selection among candidates and demonstrating improved results on economic data.

Findings

Configuration null model often preferred for economic data

Method outperforms standard correlation thresholding in predicting exports

Flexible null model integration enhances network filtering

Abstract

Network analysis has been applied to various correlation matrix data. Thresholding on the value of the pairwise correlation is probably the most straightforward and common method to create a network from a correlation matrix. However, there have been criticisms on this thresholding approach such as an inability to filter out spurious correlations, which have led to proposals of alternative methods to overcome some of the problems. We propose a method to create networks from correlation matrices based on optimisation with regularization, where we lay an edge between each pair of nodes if and only if the edge is unexpected from a null model. The proposed algorithm is advantageous in that it can be combined with different types of null models. Moreover, the algorithm can select the most plausible null model from a set of candidate null models using a model selection criterion. For three economic data sets, we find that the configuration model for correlation matrices is often preferred to standard null models. For country-level product export data, the present method better predicts main products exported from countries than sample correlation matrices do.

Figures

Tables

Table 1. Table 1: Normalised EBIC values for correlation matrices and precision matrices combined with the different null models. We divide the EBIC value for each null model by that for the sample correlation matrix. Matrices 𝑪 WN superscript 𝑪 WN {\bm{C}}^{\text{WN}} , 𝑪 HQS superscript 𝑪 HQS {\bm{C}}^{\text{HQS}} and 𝑪 con superscript 𝑪 con {\bm{C}}^{\text{con}} indicate the white noise model, HQS model and configuration model, respectively, as null models. For null precision matrices, we adopt the inverse of the three null correlation matrices and construct networks using a variant of the Scola (Appendix  C ).

		Null model
Data		Correlation matrix				Precision matrix
		${𝑪}^{\text{WN}}$	${𝑪}^{\text{HQS}}$	${𝑪}^{\text{con}}$		${𝑪}^{\text{WN}}$	${𝑪}^{\text{HQS}}$	${𝑪}^{\text{con}}$
Product space
‘72–‘91		0.226	0.202	0.187		0.232	0.190	0.208
‘92–‘11		0.300	0.301	0.226		0.304	0.238	0.256
S&P 500
‘00–‘07		0.600	0.562	0.553		0.640	0.557	0.587
‘08–‘15		0.664	0.558	0.524		0.596	0.539	0.607
Nikkei		0.991	0.874	0.861		1.001	0.837	0.882

Equations

$$\mathbf{C}={\mathbf{C}}^{\text{null}}+\mathbf{W}.$$

$$P\left(\mathbf{X}|\mathbf{C}\right)$$

$$\ln P\left(\mathbf{X}|\mathbf{C}\right)$$

$$\mathcal{L}\left(\mathbf{W}\right)$$

$${\mathbf{W}}^{\text{MLE}}\equiv {\mathbf{C}}^{\text{sample}}-{\mathbf{C}}^{\text{null}}.$$

$$\hat{\mathcal{L}}\left(\mathbf{W}|\mathbf{\lambda }\right)\equiv \mathcal{L}\left(\mathbf{W}\right)-\frac{L}{2}\sum _{i=1}^N\sum _{j=1}^{i-1}{\lambda }_{ij}|W_{ij}|,$$

$${\lambda }_{ij}=\overline{\lambda }|W_{ij}^{\text{MLE}}{|}^{-\gamma },$$

$$-\int P(\tilde{\mathbf{X}})\ln P(\tilde{\mathbf{X}})\mathrm{d}\tilde{\mathbf{X}},$$

$$\int {\tilde{C}}_{ii}^{\text{sample}}P(\tilde{\mathbf{X}})\mathrm{d}\tilde{\mathbf{X}}=C_{ii}^{\text{sample}},\text{for 1≤i≤N}.$$

$$\int \sum _{j=1}^N{\tilde{C}}_{ij}^{\text{sample}}P(\tilde{\mathbf{X}})\mathrm{d}\tilde{\mathbf{X}}=\sum _{j=1}^N{C}_{ij}^{\text{sample}},\text{for 1≤i≤N}.$$

$$\underset{{\mathbf{C}}^{\text{con}}}{\max }\ln \det \left({\mathbf{C}}^{\text{con}}\right)$$

$$C_{ii}^{\text{con}}=C_{ii}^{\text{sample}}\text{and}\sum _{j=1}^N{C}_{ij}^{\text{con}}=\sum _{j=1}^N{C}_{ij}^{\text{sample}},\text{for 1≤i≤N}.$$

$${\text{EBIC}}_{\beta }=-2\ln P\left(\mathbf{X}|{\hat{\mathbf{C}}}^{(\overline{\lambda })}\right)+\left(M+K\right)\ln L+4\beta \left(M+K\right)\ln N,$$

$$F_{\overline{\mathbf{W}}}\left(\mathbf{W}\right)\equiv -\frac{L}{2}\ln \det \left({\mathbf{C}}^{\text{null}}+\overline{\mathbf{W}}\right)-\frac{L}{2}\text{tr}\left[{\left({\mathbf{C}}^{\text{null}}+\overline{\mathbf{W}}\right)}^{-1}\left({\mathbf{C}}^{\text{null}}+\mathbf{W}\right)\right]$$

$$-\frac{L}{2}\text{tr}\left[{\mathbf{C}}^{\text{sample}}{\left({\mathbf{C}}^{\text{null}}+\mathbf{W}\right)}^{-1}\right]-\frac{NL}{2}\ln (2\pi )-\frac{L}{2}\sum _{i=1}^N\sum _{j=1}^{i-1}{\lambda }_{ij}|W_{ij}|+\frac{NL}{2}.$$

$${\tilde{\mathbf{W}}}^{(k+1)}$$

$$\displaystyle\left({\cal S}_{\epsilon{\bm{\lambda}}}(Z)\right)_{ij}=\left\{\begin{array}[]{ll}Z_{ij}-\epsilon\lambda_{ij}&(\epsilon\lambda_{ij}<Z_{ij}),\\ 0&(-\epsilon\lambda_{ij}\leq Z_{ij}\leq\epsilon\lambda_{ij}),\\ Z_{ij}+\epsilon\lambda_{ij}&(Z_{ij}<-\epsilon\lambda_{ij}).\end{array}\right.$$

$$R_{cp}^{(t)}\equiv \frac{V_{cp}^{(t)}/{\sum }_{p^{\prime }=1}^{N_p}{V}_{c{p}^{\prime }}^{(t)}}{{\sum }_{c^{\prime }=1}^{N_c}{V}_{c^{\prime }p}^{(t)}/{\sum }_{c^{\prime }=1}^{N_c}{\sum }_{p^{\prime }=1}^{N_p}{V}_{c^{\prime }{p}^{\prime }}^{(t)}},$$

$${\mathbf{x}}_{(c,t)}\equiv \left[{\overline{R}}_{c1}^{(t)},\dots ,{\overline{R}}_{c{N}_p}^{(t)}\right],$$

$${\mathbf{C}}^{\text{sample}}$$

$$\displaystyle=\left[\begin{array}[]{cc}{\bm{C}}^{\text{sample}}_{0,0}&{\bm{C}}^{\text{sample}}_{0,+10}\\ ({\bm{C}}^{\text{sample}}_{0,+10})^{\top}&{\bm{C}}^{\text{sample}}_{+10,+10}\end{array}\right].$$

$$\mathbb{E}[{\mathbf{x}}_{(c,t+10)}|{\mathbf{x}}_{(c,t)}]$$

$$={\mathbf{x}}_{(c,t)}{\mathbf{C}}_{0,0}^{-1}{\mathbf{C}}_{0,+10},$$

$${\overline{\lambda }}_{\text{upper}}^{(1)}=\underset{\begin{array}{c}i,j\\ i>j\end{array}}{\max }\left|{\left[{\left({\mathbf{C}}^{\text{null}}\right)}^{-1}-{\left({\mathbf{C}}^{\text{null}}\right)}^{-1}{\mathbf{C}}^{\text{sample}}{\left({\mathbf{C}}^{\text{null}}\right)}^{-1}\right]}_{ij}\times {\left(W_{ij}^{\text{MLE}}\right)}^2\right|.$$

$$\mathbf{\Theta }\equiv {\mathbf{\Theta }}^{\text{null}}+\mathbf{W},$$

$$\hat{\mathcal{L}}\left(\mathbf{W}|\mathbf{\lambda }\right)$$

$$-\frac{L}{2}\sum _{i=1}^N\sum _{j=1}^{i-1}{\lambda }_{ij}|W_{ij}|.$$

$${\tilde{\mathbf{W}}}^{(k+1)}$$

Full text

I Introduction

Many networks have been constructed from correlation matrices. For instance, asset graphs are networks in which a node represents a stock of a company and an edge between a pair of nodes indicates strong correlations between two stock prices Onnela2004 ; Tse2010 . A variety of tools for network analysis, such as centralities, network motifs and community structure, can be used for studying properties of correlation networks Newman2018 . Network analysis may provide information that is not revealed by other analysis methods for correlation matrices such as principal component analysis and factor analysis. Network analysis on correlation data is commonly accepted across various domains Onnela2004 ; Rubinov2010 ; VanWijk2010 ; Zanin2012 ; DeVicoFallani2014 ; Kose2001 .

The present paper proposes a new method for constructing networks from correlation matrices. There exist various methods for generating sparse networks from correlation matrices such as those based on the minimum spanning tree Mantegna1999 and its variant Tumminello2005 . Perhaps the most widely used technique is thresholding, i.e., retaining pairs of variables as edges if and only if the correlation or its absolute value is larger than a threshold value Kose2001 ; Onnela2004 ; Rubinov2010 ; VanWijk2010 ; Zanin2012 ; DeVicoFallani2014 . Because the structure of the generated networks may be sensitive to the threshold value, a variety of criteria for choosing the threshold value have been proposed Onnela2004 ; Rubinov2010 ; VanWijk2010 ; Zanin2012 ; DeVicoFallani2017 . An alternative is to analyse a collection of networks generated with different threshold values Zanin2012 .

The thresholding method is problematic when large correlations do not imply dyadic relationships between nodes Jackson1991 ; Vigen2015 . For example, in a time series of stock prices, global economic trends (e.g., recession and inflation) simultaneously affect different stock prices, which can lead to a large correlation between various pairs of stocks Bouchaud2003 ; MacMahon2015 . Additionally, the correlation between two nodes may be accounted for by other nodes. A major instance of this phenomenon is that, if node $v_{1}$ is strongly correlated with nodes $v_{2}$ and $v_{3}$, then $v_{2}$ and $v_{3}$ would be correlated even if there is no direct relationship between them Marrelec2006 ; Zalesky2012 ; Masuda2018b . For example, the murder rate (corresponding to node $v_{2}$) is positively correlated with the ice cream sales (node $v_{3}$), which is accounted for by the fluctuations in the temperature (node $v_{1}$), i.e., people are more likely to interact with others and purchase ice creams when it is hot Vigen2015 .

The sparse covariance estimation Bien2011 and graphical Lasso Banerjee2008 ; Friedman2008 ; Fan2009 ; Foygel2010 ; VanBorkulo2014 were previously proposed for filtering out spurious correlations. These methods consider the correlation between nodes to be spurious if the correlation is accounted for by random fluctuations under a white noise null model, which assumes that observations at all nodes are independent of each other. However, nodes tend to be correlated with each other in empirical data owing to trivial factors (e.g., global trend), which calls for different null models that emulate different types of spurious correlations Laloux1999 ; Plerou1999 ; Hirschberger2007 ; MacMahon2015 ; Masuda2018 .

To incorporate such null models into network inference, we present an algorithm named the Scola, standing for Sparse network construction from COrrelational data with LAsso. The Scola places edges between node pairs if the correlations are not accounted for by a null model of choice. A main advantage of the Scola is that it leaves the choice of null models to users, enabling them to filter out different types of trivial relationships between nodes. Furthermore, the Scola can select the most plausible null model for the given data among a set of null models using a model selection framework. A Python code for the Scola is available at scola-code

II Methods

II.1 Construction of networks from correlation matrices

Consider $N$ variables, which we refer to as nodes. We aim to construct a network on $N$ nodes, in which edges indicate the correlations that are not attributed to some trivial properties of the system. Let ${\bm{C}}=({C}_{ij})$ be the $N\times N$ correlation matrix, where ${C}_{ij}$ is the correlation between nodes $i$ and $j$, i.e., $-1\leq{C}_{ij}\leq 1$ and ${C}_{ii}=1$ for $1\leq i\leq N$. We write ${\bm{C}}$ as

[TABLE]

Matrix ${\bm{C}}^{\text{null}}=\left({C_{ij}^{\text{null}}}\right)$ is an $N\times N$ correlation matrix, where ${C_{ij}^{\text{null}}}$ is the mean value of the correlation between nodes $i$ and $j$ under a null model. For example, if every node is independent of each other under the null model, then one sets ${\bm{C}}^{\text{null}}={\bm{I}}$, where ${\bm{I}}$ is the $N\times N$ identity matrix. We introduce three null models in Section II.2. Matrix ${{\bm{W}}}=\left({W_{ij}}\right)$ is an $N\times N$ matrix representing the deviation from the null model. We place an edge between nodes $i$ and $j$ (i.e., $W_{ij}\neq 0$) if and only if the correlation $C_{ij}$ is sufficiently different from that for the null model. We note that edges are undirected and may have positive or negative weights. The network is assumed not to have self-loops (i.e., ${W_{ii}}=0$) because the diagonal elements of ${\bm{C}}^{\text{null}}$ and ${\bm{C}}^{\text{sample}}$ are always equal to one.

We estimate ${{\bm{W}}}$ from data as follows. Assume that we have $L$ samples of data observed at the $N$ nodes, based on which the correlation matrix is calculated. Denote by ${\bm{X}}=(x_{\ell i})$ the $L\times N$ matrix, in which $x_{\ell i}$ is the value observed at node $i$ in the $\ell$th sample. Let ${\bm{x}}_{\ell}=\left[x_{\ell 1},x_{\ell 2},\ldots,x_{\ell N}\right]$ be the $\ell$th sample. We assume that each sample ${\bm{x}}_{\ell}$ is independently and identically distributed according to an $N$-dimensional multivariate Gaussian distribution with mean zero. For mathematical convenience, we assume that ${\bm{X}}$ is preprocessed such that the average and the variance of $x_{\ell i}$ over the $L$ samples are zero and one respectively, i.e., $\sum_{\ell=1}^{L}x_{\ell i}/L=0$ and $\sum_{\ell=1}^{L}x_{\ell i}^{2}/L=1$ for $1\leq i\leq N$.

Our goal is to find ${{\bm{W}}}$ that maximises the likelihood of ${\bm{X}}$, i.e.,

[TABLE]

where $\top$ is the transposition. It should be noted that ${\bm{C}}$ may not be equal to the sample Pearson correlation matrix, denoted by ${\bm{C}}^{\text{sample}}\equiv{\bm{X}}^{\top}{\bm{X}}/L$ Kollo2006 , if $L$ is finite. The log likelihood is given by

[TABLE]

Using $\sum_{\ell=1}^{L}{\bm{x}}_{\ell}{\bm{C}}^{-1}{\bm{x}}_{\ell}^{\top}=\text{tr}\left({\bm{X}}^{\top}{\bm{X}}{\bm{C}}^{-1}\right)$, one obtains

[TABLE]

Substitution of Eq. (1) into Eq. (4) leads to the log likelihood of ${{\bm{W}}}$ as follows.

[TABLE]

The log likelihood ${\cal L}$ is a concave function with respect to ${{\bm{W}}}$. Therefore, one obtains the maximiser of ${\cal L}$, denoted by ${{\bm{W}}}^{\text{MLE}}$, by solving $\partial{\cal L}/\partial{{\bm{W}}}=0$, i.e.,

[TABLE]

In practice, ${{\bm{W}}}^{\text{MLE}}$ overfits the given data, leading to a network with many spurious edges. This is because the number of samples, $L$, is often smaller than the number of elements in ${{\bm{W}}}^{\text{MLE}}$, i.e., $L<N(N-1)/2$, as is the case for the estimation of correlation matrices and precision matrices Banerjee2008 ; Bien2011 . To prevent overfitting, we impose penalties on the number of non-zero elements in ${{\bm{W}}}$ using the Lasso Tibshirani1996 . The Lasso is commonly used for regression analysis to obtain a model with a small number of non-zero regression coefficients. The Lasso is also used for estimating sparse covariance matrices (i.e., ${\bm{C}}$) Bien2011 and sparse precision matrices (i.e., ${\bm{C}}^{-1}$) Banerjee2008 with a small number of samples. (We discuss these methods in Section IV.) Here we apply the Lasso to obtain a sparse ${{\bm{W}}}$. Specifically, we maximise penalized likelihood function

[TABLE]

where $\lambda_{ij}\geq 0$ is the Lasso penalty for ${W_{ij}}$. Large values of $\lambda_{ij}$ yield sparse ${{\bm{W}}}$. Because $\hat{\cal L}$ is not concave with respect to ${{\bm{W}}}$, one cannot analytically find the maximum of $\hat{\cal L}$. Therefore, we numerically maximise $\hat{\cal L}$ using an extension of a previous algorithm Bien2011 , which is described in Section II.4.

The penalized likelihood $\hat{\cal L}$ contains $N(N-1)/2$ Lasso penalty parameters, $\lambda_{ij}$. A simple choice is to use the same value for all $\lambda_{ij}$. However, this method is problematic Fan2001 ; Zou2006 . If one imposes the same penalty to all node pairs, one tends to obtain either a sparse network with small edge weights or a dense network with large edge weights. However, sparse networks with large edge weights or dense networks with small edge weights may yield a larger likelihood. A remedy for this problem is the adaptive Lasso Zou2006 , which sets

[TABLE]

where $\overline{\lambda}\geq 0$ and $\gamma>0$ are hyperparameters. With the adaptive Lasso, a small penalty is imposed on a pair of nodes $i$ and $j$ if $W^{\text{MLE}}_{ij}$ is far from zero, allowing the edge to have a large weight. If one has sufficiently many samples, the adaptive Lasso correctly identifies zero and non-zero regression coefficients (i.e., ${W_{ij}}$ in our case) given an appropriate $\overline{\lambda}$ value and any positive $\gamma$ value Zou2006 . The estimated ${W_{ij}}$ values did not much depend on $\gamma$ in our numerical experiments. Therefore, we set $\gamma=2$, which is a typical value Zou2006 ; Fan2009 ; Hui2015 . Hyperparameter $\overline{\lambda}$ controls the number of edges in the network (i.e., the number of nonzero elements in ${{\bm{W}}}$). A large $\overline{\lambda}$ yields sparse networks. We describe how to choose $\overline{\lambda}$ in Section II.3.

II.2 Null models for correlation matrices

The Scola accepts various null correlation matrices, i.e., ${\bm{C}}^{\text{null}}$. Nevertheless, we mainly focus on the configuration model for correlation matrices Masuda2018 . Although we also examined two other null models in numerical experiments, the configuration model was mostly chosen in the model selection (Section III.2).

Our configuration model is based on the principle of maximum entropy. With this method, one determines the probability distribution of data, denoted by $P(\tilde{\bm{X}})$, by maximising the following Shannon entropy given by

[TABLE]

with respect to $P(\tilde{\bm{X}})$ under constraints, where $\tilde{\bm{X}}=(\tilde{x}_{\ell i})$ is an $L\times N$ matrix such that $\tilde{x}_{\ell i}$ is the value at node $i$ in the $\ell$th sample. Let $\tilde{\bm{C}}^{\text{sample}}\equiv\tilde{\bm{X}}^{\top}\tilde{\bm{X}}/L$ be the sample covariance matrix for $\tilde{\bm{X}}$. In the configuration model, we impose that each node has the same expected variance as that for the original data, i.e.,

[TABLE]

We also impose that the row sum (equivalently, the column sum) of ${\bm{C}}^{\text{sample}}$ is preserved, i.e.,

[TABLE]

Equation (11) is analogous to the case of the configuration model for networks, which by definition preserves the row sum of the adjacency matrix of a network, or equivalently the degree of each node. We note that the $i$th row sum of ${\bm{C}}^{\text{sample}}$ is proportional to the correlation between the observation at the $i$th node and the average of the observations over all nodes MacMahon2015 .

Denote by ${\bm{C}}^{\text{con}}$ the expectation of $\tilde{\bm{C}}^{\text{sample}}$. We compute ${\bm{C}}^{\text{con}}$, which we use as the null correlation matrix (i.e., ${\bm{C}}^{\text{null}}$), as follows. In Ref. Masuda2018 , we showed that $P(\tilde{\bm{X}})$ is a multivariate Gaussian distribution with mean zero. Under $P(\tilde{\bm{X}})$, ${\bm{C}}^{\text{con}}$ is equal to the variance parameter of $P(\tilde{\bm{X}})$ Gupta2000 ; Kollo2006 ; Masuda2018 . By substituting Eq. (2) into Eqs. (9), (10) and (11), we rewrite the maximisation problem as

[TABLE]

subject to

[TABLE]

Equation (12) is concave with respect to ${\bm{C}}^{\text{con}}$. Moreover, the feasible region defined by Eq. (13) is a convex set. Therefore, the maximisation problem is a convex problem such that one can efficiently find the global optimum. We compute ${\bm{C}}^{\text{con}}$ using an in-house Python program available at dmccalg .

In addition to the configuration model, we consider two other null models. The first model is the white noise model, in which the signal at each node is independent of each other and has the same variance with that in the original data. The null correlation matrix for the white noise model, denoted by ${\bm{C}}^{\text{WN}}$, is given by ${C}^{\text{WN}}_{ij}=0$ for $i\neq j$ and ${C}^{\text{WN}}_{ii}=1$ for $1\leq i\leq N$. The white noise model is often used in the analysis of correlation networks Fisher1915 ; MacMahon2015 .

Another null model is the Hirschberger-Qi-Steuer (HQS) model Hirschberger2007 , in which each node pair is assumed to be correlated to the same extent as expectation. As is the case for the configuration model for correlation matrices, the original HQS model provides a probability distribution of covariance matrices. The HQS model preserves the variance of the signal at each node averaged over all the ndoes as expectation. Moreover, the HQS model preserves the average and variance of the correlation values over different pairs of nodes in the original correlation matrix as expectation. The HQS model is analogous to the Erdős-Rényi random graph for networks, in which each pair of nodes is adjacent with the same probability Masuda2018 . We use the expectation of the correlation matrix generated by the HQS model, denoted by ${\bm{C}}^{\text{HQS}}$. We obtain ${C}_{ij}^{\text{HQS}}=\left[N(N-1)/2\right]^{-1}\sum_{i^{\prime}=1}^{N}\sum_{j^{\prime}=1}^{i^{\prime}-1}C_{i^{\prime}j^{\prime}}^{\text{sample}}$ for $i\neq j$ and ${C}_{ii}^{\text{HQS}}=1$ Hirschberger2007 .

II.3 Model selection

We determine the value of $\overline{\lambda}$ based on a model selection criterion. Commonly used criteria, Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC), favour an excessively rich model if the model has many parameters relative to the number of samples Chen2008 . As discussed in Section II.1, this is often the case for the estimation of correlation matrices Banerjee2008 ; Foygel2010 ; Bien2011 .

We use the extended BIC (EBIC) to circumvent this problem Chen2008 ; Foygel2010 . Let $\hat{\bm{C}}^{(\overline{\lambda})}$ be the estimated correlation matrix at $\overline{\lambda}$, i.e., $\hat{\bm{C}}^{(\overline{\lambda})}={\bm{C}}^{\text{null}}+{{\bm{W}}}^{(\overline{\lambda})}$, where ${{\bm{W}}}^{(\overline{\lambda})}$ is the network one estimates by maximising $\hat{\cal L}({{\bm{W}}}|\overline{\lambda})$. We adopt $\hat{\bm{C}}^{(\overline{\lambda})}$ that minimises

[TABLE]

with respect to $\overline{\lambda}$. In Eq. (14), $M$ is the number of edges in the network (i.e., the half of the number of nonzero elements in ${{\bm{W}}}$), $K$ is the number of parameters of the null model and $\beta$ is a parameter. The white noise model introduced in Section II.2 does not have parameters. Therefore, $K=0$. The HQS model has $K=1$ parameter, i.e., the average of the off-diagonal elements. The configuration model yields $K=N$, i.e., the row sum of each node. It should be noted that we compute the number of parameters by exploiting the fact that ${\bm{C}}^{\text{sample}}$ is a correlation matrix, i.e., the diagonal entries of ${\bm{C}}^{\text{sample}}$ are always equal to one.

Parameter $\beta\in[0,1]$ determines the prior distribution for a Bayesian inference. The prior distribution affects the sparsity of networks; a large $\beta$ value would yield sparse networks. Nevertheless, the effect of the prior distribution on the EBIC value diminishes when the number of samples (i.e., $L$) increases. We adopt a typical value, i.e., $\beta=0.5$, which provided reasonable results for linear regressions and the estimation of precision matrices Chen2008 ; Foygel2010 ; VanBorkulo2014 . We adopt the golden-section search method to find the $\overline{\lambda}$ value that yields the minimum EBIC value Press2007 (Appendix A).

In addition to selecting the $\overline{\lambda}$ value, the EBIC can be used for selecting a null model among different types of null models. In this case, we compute ${{\bm{W}}}$ for each null model with the $\overline{\lambda}$ value determined by the golden-section search method. Then, we select the pair of ${\bm{C}}^{\text{null}}$ and ${{\bm{W}}}$ that minimises the EBIC value.

II.4 Maximising the penalized likelihood

To maximise $\hat{\cal L}$ in terms of ${{\bm{W}}}$, we use the minorise-maximise (MM) algorithm Bien2011 . Although the MM algorithm may not find the global maximum, it converges to a local optimum. The MM algorithm starts from initial guess ${{\bm{W}}}={\bm{0}}$ and iterates rounds of the following minorisation step and the maximisation step.

In the minorisation step, we approximate $\hat{\cal L}$ around the current estimate $\overline{{\bm{W}}}$ by Bien2011

[TABLE]

It should be noted that $F_{\overline{{\bm{W}}}}$ is a minoriser of $\hat{\cal L}$ satisfying $F_{\overline{{\bm{W}}}}(\overline{{\bm{W}}})=\hat{\cal L}(\overline{{\bm{W}}})$ and $F_{\overline{{\bm{W}}}}({{\bm{W}}})\leq\hat{\cal L}({{\bm{W}}})$ for all ${{\bm{W}}}$ provided that ${\bm{C}}^{\text{null}}+{{\bm{W}}}$ is positive definite Bien2011 . This property of $F_{\overline{{\bm{W}}}}$ ensures that maximising $F_{\overline{{\bm{W}}}}$ yields a $\hat{\cal L}$ value larger than or equal to $\hat{\cal L}(\overline{{\bm{W}}})$.

In the maximisation step, we seek the maximiser of $F_{\overline{{\bm{W}}}}$. Function $F_{\overline{{\bm{W}}}}$ is a concave function, which allows us to find the global maximum with a standard gradient descent rule for the Lasso. Specifically, starting from $\widetilde{{\bm{W}}}^{(0)}=\overline{{\bm{W}}}$, we iterate the following update rule until convergence:

[TABLE]

where $\widetilde{{\bm{W}}}^{(k)}$ is the tentative solution at the $k$th iteration, $\epsilon$ is the learning rate, and ${\cal S}$ is the element-wise soft threshold operator given by

[TABLE]

If $\hat{\cal L}(\widetilde{{\bm{W}}}^{(k_{\text{final}})})>\hat{\cal L}(\overline{{\bm{W}}})$, where $k_{\text{final}}$ is the iteration number after a sufficient convergence, we set $\overline{{\bm{W}}}=\widetilde{{\bm{W}}}^{(k_{\text{final}})}$. Then, we perform another round of a minorisation step and a maximisation step. Otherwise, the algorithm finishes.

The learning rate $\epsilon$ mainly affects the speed of the convergence of the iterations in the maximisation step. We set $\epsilon$ using the ADAptive Moment estimation (ADAM), a gradient descent algorithm used in various machine learning algorithms Kingma2014 . ADAM adapts the learning rate at each update (Eq. (16)) based on the current and past gradients.

The MM algorithm requires computational time of ${\cal O}(N^{3})$ because Eq. (16) involves the inversion of $N\times N$ matrices. In our numerical experiments (Section III), the entirely of the Scola consumed approximately three hours on average for a network with $N=$1,930 nodes with 16 threads running on the Intel 2.6 GHz Sandy Bridge processors with 4GB memory.

III Numerical results

III.1 Prediction of country-level product exports

The product space is a network of products, where a pair of products is defined to be adjacent if both of them have a large share in the export volume of the same country Hidalgo2007 ; Hartmann2017 . Here we construct the product space with the Scola and use it for predicting product exports from countries.

We use the data set provided by the Observatory of Economic Complexity product-space , which contains annual export volumes of $N_{p}=988$ products from $N_{c}=249$ countries between 1962 and 2014. The data set also contains the Standard International Trade Classification (SITC) code for each product, which indicates the product type. We quantify the level of sophistication of the product types using the PRODY index Hausmann2007 . We compute the PRODY index in 1991 for each product, where a product with a large PRODY index is considered to be sophisticated. Then, we average the PRODY index over the products of the same type. The product types in descending order of the PRODY index are as follows: “Machinery & transport equipment”, “Chemicals”, “Miscellaneous manufactured articles”, “Manufactured goods by material”, “Miscellaneous”, “Mineral fuels, lubricants & related materials”, “Beverages & tobacco”, “Animal & vegetable oils, fats & waxes”, “Food & live animals” and “Raw materials”.

In previous studies, the product space was constructed based on the so-called revealed comparative advantage (RCA) defined by

[TABLE]

where $V_{cp}^{(t)}$ is the total export volume of product $p$ from country $c$ in year $t$. Then, it was assumed that two products $p$ and $p^{\prime}$ were adjacent in the product space if and only if $R_{cp}^{(t)},R_{cp^{\prime}}^{(t)}>1$ for at least one country $c$ Hidalgo2007 .

In contrast to this approach, we construct the product space as follows. First, countries may export different products in different years. To mitigate the effect of temporal changes, we split the data into two halves, i.e., those between $T_{\text{s}}=1972$ and $T_{\text{f}}=1991$, and those between $T_{\text{s}}=1992$ and $T_{\text{f}}=2011$.

Second, for each time window, we apply the Box-Cox transformation Box1964 to each $R_{cp}^{(t)}$ to make the distribution closer to a standard normal distribution (Appendix B). This preprocessing is crucial because the Scola assumes that the given data, ${\bm{X}}$, is distributed according to a multivariate Gaussian distribution. We note that the sample average and variance of the transformed RCA values based on the $N_{p}$ products are equal to zero and one, respectively.

Third, we define a sample ${\bm{x}}_{(c,t)}$ for each combination of country $c$ and year $t$ by

[TABLE]

where $\overline{R}_{cp}^{(t)}$ is the transformed value of $R_{cp}^{(t)}$.

Fourth, we compute a sample Pearson correlation matrix for concatenated samples $\left[{\bm{x}}_{(c,t)},{\bm{x}}_{(c,t+10)}\right]$. Specifically, for each of the two time windows, we compute ${\bm{C}}^{\text{sample}}$ by

[TABLE]

In Eq. (25), ${\bm{C}}^{\text{sample}}_{0,0}$ and ${\bm{C}}^{\text{sample}}_{+10,+10}$ are the correlation matrices for the products within the same year. The off-diagonal block ${\bm{C}}^{\text{sample}}_{0,+10}$ contains the correlations between the products with a time lag of ten years.

Fifth, we generate networks by applying either the thresholding method or the Scola to ${\bm{C}}^{\text{sample}}$. For the Scola, we adopt the configuration model as the null model. For the thresholding method, we set the threshold value such that the number of edges in the network is equal to that in the network generated by the Scola. We set the weight of each edge to one for the network generated by the thresholding method.

The sample correlation matrices are shown in Figs. 1(a) and 1(b). The thresholding method places a majority of edges in the on-diagonal blocks for both time windows (Figs. 1(c) and (d)), suggesting strong correlations between the products exported within the same year as compared to those between different years. This is also true for the networks generated by the Scola (Figs. 1(e) and (f)).

The two methods place few edges (less than 5%) within the off-diagonal blocks for ‘72–‘91. For ‘92–‘11, more than 22% of edges are placed within the off-diagonal blocks by both methods. Although both methods place a similar number of edges within the off-diagonal blocks, the distribution of edges is considerably different. To see this, we compute the fraction of edges between product types for ‘92–‘11 (Fig. 2). We do not show the result for ‘72–‘91 owing to a small number of edges within the off-diagonal blocks. We find that, within the off-diagonal blocks, the thresholding method places relatively many edges between two nodes that correspond to sophisticated products in terms of the PRODY index (Fig. 2(a)). Examples include “Machinery & transport equipment”, “Chemicals” and “Manufactured goods by material”. This result suggests that sophisticated products are strongly correlated with the same or other sophisticated products ten years apart. In contrast, the Scola finds many edges between nodes that correspond to less sophisticated products such as “Raw materials”, “Foods & live animals” and “Animal & vegetable oils, fats & waxes” (Fig. 2 (b)). We find a similar result for the on-diagonal blocks for ‘92–‘11 (Fig. 1(f)).

Highly correlated nodes may not be adjacent in the network generated by the Scola. To examine this issue, we plot the weight of edges estimated by the Scola against the correlation value of the corresponding node pair in Fig. 3. For time window ‘72–‘91, there are 3,021 node pairs with a correlation above the threshold in magnitude, of which 1,700 (56%) node pairs are not adjacent in the network generated by the Scola. We find qualitatively the same result for time window ‘92–‘11; there are 6,834 node pairs with a correlation above the threshold in magnitude, of which 6,343 node pairs ($93\%$) are not adjacent in the network generated by the Scola. The weights of edges are correlated strongly with the original correlation values for time window ‘72–‘91 but weakly for ‘92–‘11 (Spearman correlation coefficients are 0.87 and 0.21, respectively).

We further demonstrate the use of the generated networks for predictions. We aim to predict RCA values in year $t+10$ given those in year $t$. We make predictions using vector autoregression Kilian2006 , which consists in computing conditional probability distribution $P({\bm{x}}_{(c,t+10)}|{\bm{x}}_{(c,t)})=P({\bm{x}}_{(c,t)},{\bm{x}}_{(c,t+10)})/P({\bm{x}}_{(c,t)})$. Note that the joint probability distribution $P({\bm{x}}_{(c,t)},{\bm{x}}_{(c,t+10)})$ is given by Eq. (2). The joint probability distribution $P({\bm{x}}_{(c,t)},{\bm{x}}_{(c,t+10)})$ has covariance matrix ${\bm{C}}$ as parameter. We set ${\bm{C}}$ to either the sample correlation matrix (${\bm{C}}={\bm{C}}^{\text{sample}}$) or that provided by the Scola (i.e., ${\bm{C}}={\bm{C}}^{\text{null}}+{{\bm{W}}}$). We do not use the thresholding method in this prediction task because the thresholding method does not provide correlation matrices. We make a prediction by the conditional expected value of ${\bm{x}}_{(c,t+10)}$ under $P({\bm{x}}_{(c,t+10)}|{\bm{x}}_{(c,t)})$, which is given by Bishop2006

[TABLE]

where ${\bm{C}}_{0,0}$ and ${\bm{C}}_{0,+10}$ are the blocks of ${\bm{C}}$ defined analogously to Eq. (25).

We carry out five-fold cross-validation, where we split sample indices $\{1,\ldots,L\}$ into five subsets of almost equal sizes. We estimate ${\bm{C}}$ using the training set, which is the union of four out of the five subsets. Then, we perform predictions for the test set, which is the remaining subset. We carry out this procedure five times such that each of the subsets is used once as the test set.

The joint distribution of the actual and predicted $\overline{R}_{cp}^{(t+10)}$ values, where each tuple $(c,p,t)$ is regarded as an individual data point, is shown in Fig. 4. (The joint distributions for other null models are shown in the Supplementary Materials.) The perfect prediction would yield all points on the diagonal line. Between ‘72–‘91, the Scola realises better predictions than the sample correlation matrix does (Figs. 4(a) and 4(c)). In fact, the mean squared error (MSE) is approximately three times smaller for the Scola than for the sample correlation matrix. Between ‘92–‘11, the MSE for the Scola is approximately 1.5 times smaller than that for the sample correlation matrix (Figs. 4(b) and 4(d)). A probable reason for the poor prediction performance for the sample correlation matrices is overfitting. Matrix ${\bm{C}}^{\text{sample}}$ contains more than $1.9\times 10^{6}$ elements, whereas there are only $L=$1,400 samples on average. The Scola represents the correlation matrix with a considerably smaller number of parameters (i.e., $M+K=$8,185 on average), which mitigates overfitting.

III.2 Model selection

What are appropriate null models for correlation networks? To address this question, we carry out model selection based on the EBIC to compare the Scola with different null models. We examine the three null models, i.e., the white noise model, the HQS model and the configuration model. We also compare the performance of the Scola with estimators of sparse precision matrices with different null models.

To construct a network from a precision matrix, we use a variant of the Scola (Appendix C). We adopt the inverse of the correlation matrices for the white noise model, the HQS model and the configuration model as the null models for precision matrices. Because the white noise model is the identity matrix, its inverse is also the identity matrix, providing the white noise model for precision matrices. It should be noted that the variant is equivalent to the graphical Lasso Banerjee2008 ; Friedman2008 when one uses the white noise model as the null model.

In addition to the country-level export data used in Section III.1, we use the time series of stock prices in the US and Japanese markets. The US data comprise the time series of the daily closing prices of $N=$1,023 companies in the list of the Standard & Poor’s 500 (S&P 500) index on 4,174 days between 01/01/2000 and 29/12/2015. We compute the logarithmic return of the daily closing price by $x_{ti}\equiv\log\left[z_{i}(t+1)/z_{i}(t)\right]$, where $z_{i}(t)$ is the closing price of stock $i$ on day $t$. We split the time series of $x_{ti}$ into two halves of eight years. Then, for each half, we exclude the stocks that have at least one missing value. (The stock prices are missing at the time points at which companies are not included in S&P 500.) We compute the correlation matrix for the logarithmic returns of the remaining stocks for each half.

The Japanese data comprise the time series of the daily closing prices of 264 stocks in the first section of the Tokyo Stock Exchange. We retrieve the stock prices for 5,324 days between 12/03/1996 and 29/02/2016 using Nikkei NEED NIKKEI , where we exclude the stocks that do not have transactions on at least one day during the period. As is the case for the US data, we compute the correlation matrix for the logarithmic returns. We used the same correlation matrix for the Japanese data in our previous study Masuda2018 . We refer to the US and Japanese stock data as S&P 500 and Nikkei, respectively.

The EBIC values for the networks generated by the Scola and its variant for precision matrices combined with the three null models are shown in Table 1. For the product space and S&P 500, the EBIC value for the configuration model is the smallest. For Nikkei, the EBIC value for the HQS model for precision matrices is the smallest. This result suggests that the Scola does not always outperform its variant for precision matrices. It should be noted that the graphical Lasso, which is equivalent to the variant of the Scola for precision matrices combined with the white noise model, is among the poorest across the different data sets.

IV Discussion

We presented the Scola to construct networks from correlation matrices. We defined two nodes to be adjacent if the correlation between them is significantly different from that expected for a null model. The Scola yielded insights that were not revealed by the thresholding method such as a positive correlation between less sophisticated products across a decade. The generated networks also better predicted country-level product exports after ten years than the mere sample correlation matrices.

Null models that have to be fed to the Scola are not limited to the three models introduced in Section II.2. Another major family of null models for correlation matrices is those based on random matrix theory, which preserves a part of spectral properties of given correlation matrices Vasiliki2002 ; Utsugi2003 ; MacMahon2015 . The Scola cannot employ this family of null models because it requires the null correlation matrix to be invertible. Null matrices based on random matrix theory are often not invertible because they leave out some of the eigenmodes. A remedy to this problem is to use a pseudo inverse.

The Scola is equivalent to the previous algorithm Bien2011 if one adopts the white noise model as the null model, i.e., ${\bm{C}}^{\text{null}}={\bm{C}}^{\text{WN}}$. Another method closely related to the Scola is the graphical Lasso Banerjee2008 ; Friedman2008 ; Fan2009 ; Foygel2010 ; VanBorkulo2014 , which provides sparse precision matrices (i.e., ${\bm{C}}^{-1}$). In contrast to correlation matrices, precision matrices indicate correlations between nodes with the effect of other nodes being removed. We focused on correlation matrices because null models for correlation matrices are relatively well studied Vasiliki2002 ; Utsugi2003 ; Hirschberger2007 ; MacMahon2015 ; Masuda2018 , while studies on null models for precision matrices are still absent to the best of our knowledge.

The inverse of the correlation matrices provided by the HQS model and the configuration model may be reasonable null models for precision matrices. To explore this direction, we developed a variant of the Scola for the case of precision matrices, which is equivalent to the graphical Lasso if the white noise model is the null model. The variant of the Scola generates networks better than the original Scola for the Nikkei data in terms of the EBIC (Section III.2). We do not claim that the Scola is generally a strong performer. More comprehensive comparisons of the Scola and competitors warrant future work.

Although the thresholding method has been widely employed Kose2001 ; Onnela2004 ; Rubinov2010 ; VanWijk2010 ; Zanin2012 ; DeVicoFallani2014 , the overfitting problem inherent in this method has received much less attention than it deserves. In many cases, the number of observations based on which one computes the correlation matrix is of the same order of the number of nodes, which is much smaller than the number of the entries in the correlation matrix Banerjee2008 ; Bien2011 . In this overfitting situation, if one removes or adds a small number of observations, one may obtain a substantially different correlation matrix and the resulting network.

We have assumed that the input data obey a multivariate Gaussian distribution. However, this assumption may not hold true for empirical data. A remedy commonly used in machine learning is to transform data using an exponential function, which is referred to as power transformations. The Box-Cox transformation, which we used in the analysis of the product space, is a standard power transformation. Other transformation techniques include the Fisher transformation Fisher1915 and the inverse hyperbolic transformation Borbidge1988 . Alternatively, one may assume other probability distributions for input data, as with a graphical Lasso for binary data VanBorkulo2014 .

Although we illustrated the Scola on economic data, the method can be applied to correlation data in various fields including neuroimaging Rubinov2010 ; VanWijk2010 ; Zalesky2012 ; Zanin2012 ; DeVicoFallani2014 , psychology VanBorkulo2014 , climate Tsonis2006 , metabolomics Kose2001 and genomics Delafuente2002 . For example, in neuroscience, the correlation data are often used to construct functional brain networks, analysis of which is expected to provide insights into how brains operate and cognition occurs. Applications of the Scola with the aim of finding insights that have not been obtained by thresholding methods, which has conventionally been used for these data, warrant future work.

Acknowledgments

The Standard & Poor’s 500 data were provided by CheckRisk LLP in the UK. We thank Yukie Sano for providing the Japanese stock data used in the present paper. N. M. acknowledges the support provided through JST, CREST, Grant Number JPMJCR1304.

Competing interests

The authors declare no competing interests.

Appendix A Golden-section search

We adopt the golden-section search to find the $\overline{\lambda}$ value that minimises the EBIC Press2007 . Let $\overline{\lambda}^{*}$ be the $\overline{\lambda}$ value yielding the minimum EBIC value. Let $\text{EBIC}(\overline{\lambda})$ be the EBIC value at $\overline{\lambda}$. In most cases, when one increases $\overline{\lambda}$ from 0, the EBIC monotonically decreases for $\overline{\lambda}\leq\overline{\lambda}^{*}$, reaches the minimum at $\overline{\lambda}=\overline{\lambda}^{*}$ and monotonically increases for $\overline{\lambda}>\overline{\lambda}^{*}$. The golden-section search method exploits this property to find $\overline{\lambda}^{*}$. Specifically, suppose that one knows a lower bound and an upper bound of $\overline{\lambda}^{*}$, i.e., $\overline{\lambda}_{\text{lower}}\leq\overline{\lambda}^{*}\leq\overline{\lambda}_{\text{upper}}$. Then, consider $\overline{\lambda}_{1}$ and $\overline{\lambda}_{2}$, where $\overline{\lambda}_{\text{lower}}<\overline{\lambda}_{1}<\overline{\lambda}_{2}<\overline{\lambda}_{\text{upper}}$. If $\text{EBIC}(\overline{\lambda}_{1})<\text{EBIC}(\overline{\lambda}_{2})$, it indicates $\overline{\lambda}^{*}\leq\overline{\lambda}_{2}$, yielding a tighter bound $\overline{\lambda}^{*}\in[\overline{\lambda}_{\text{lower}},\overline{\lambda}_{2}]$. In contrast, $\text{EBIC}(\overline{\lambda}_{1})>\text{EBIC}(\overline{\lambda}_{2})$ indicates $\overline{\lambda}^{*}\geq\overline{\lambda}_{1}$, yielding $\overline{\lambda}^{*}\in[\overline{\lambda}_{1},\overline{\lambda}_{\text{upper}}]$.

Based on this observation, the golden-section search method iterates the following rounds. In round $k=0$, we set the initial lower and upper bounds following a previous study on regression analysis with Lasso Friedman2010 . Specifically, we set $\overline{\lambda}^{(1)}_{\text{lower}}=0$ and $\overline{\lambda}_{\text{upper}}^{(1)}$ to the minimum value of $\overline{\lambda}$ that satisfies ${\cal S}_{\bm{\lambda}}\left(\partial F_{\overline{{\bm{W}}}}(\overline{{\bm{W}}})/\partial\overline{{\bm{W}}}\right)={\bm{0}}$ at $\overline{{\bm{W}}}={\bm{0}}$, i.e.,

[TABLE]

In round $k\geq 1$, one sets $\overline{\lambda}_{1}^{(k)}=\overline{\lambda}_{\text{lower}}^{(k)}+h^{(k)}/\phi^{2}$ and $\overline{\lambda}_{2}^{(k)}=\overline{\lambda}_{\text{lower}}^{(k)}+h^{(k)}/\phi$, where $h^{(k)}\equiv\overline{\lambda}^{(k)}_{\text{upper}}-\overline{\lambda}^{(k)}_{\text{lower}}$ and $\phi=(1+\sqrt{5})/2$ is the golden ratio. Then, one updates the bound, i.e., $(\overline{\lambda}^{(k+1)}_{\text{lower}},\overline{\lambda}^{(k+1)}_{\text{upper}})=(\overline{\lambda}_{\text{lower}}^{(k)},\overline{\lambda}^{(k)}_{2})$ if $\text{EBIC}(\overline{\lambda}_{1}^{(k)})<\text{EBIC}(\overline{\lambda}_{2}^{(k)})$ and $(\overline{\lambda}^{(k+1)}_{\text{lower}},\overline{\lambda}^{(k+1)}_{\text{upper}})=(\overline{\lambda}_{1}^{(k)},\overline{\lambda}^{(k)}_{\text{upper}})$ if $\text{EBIC}(\overline{\lambda}_{1}^{(k)})>\text{EBIC}(\overline{\lambda}_{2}^{(k)})$. If $\text{EBIC}(\overline{\lambda}_{1}^{(k)})=\text{EBIC}(\overline{\lambda}_{2}^{(k)})$, we adopt either bound pair with an equal probability. If $h^{(k+1)}<h^{(1)}/100$, we stop the rounds of iteration and take the smaller of $\text{EBIC}(\overline{\lambda}_{\text{lower}}^{(k+1)})$ or $\text{EBIC}(\overline{\lambda}^{(k+1)}_{\text{upper}})$ as the output. Otherwise, we carry out the $(k+1)$th round.

Appendix B Preprocessing RCA values

The distribution of RCA may be considerably skewed. To make the distribution closer to a normal distribution, we perform a Box-Cox transformation Box1964 , i.e., $\log\left(R^{(t)}_{cp}+10^{-8}\right)$. Then, for each product, we compute the $z$-score for the transformed values, generating normalised samples $\overline{R}^{(t)}_{cp}$ that have average zero and variance one for each product.

For the prediction task, we transform the RCA values separately for each time window as follows. First, we perform the Box-Cox transformation. Then, we split the transformed RCA values into a training set and test set. We compute the $z$-score for the training samples, generating normalised samples. For the test samples, we compute $z$-score for each product using the average and variance for the training samples because we should assume that statistics of the test samples are unknown when predicting the values of the test samples.

Appendix C Constructing networks from precision matrices

The precision matrices are the inverse of the correlation matrix and contain the correlation between nodes with the effect of other nodes being removed. The precision matrix is sensitive to noise in data, which calls for robust estimators such as the graphical Lasso Banerjee2008 ; Friedman2008 ; Fan2009 ; Foygel2010 ; VanBorkulo2014 . The graphical Lasso implicitly assumes a null model, where every node is conditionally independent of each other, as is the case for the white noise model for correlation matrices. Other null models for precision matrices have not been proposed to the best of our knowledge. Nevertheless, one may use the inverse of the null models for correlation matrices as the null models for precision matrices.

Therefore, we develop a variant of the Scola for precision matrices as follows. Denote by ${{\bm{\Theta}}}$ the precision matrix (i.e., ${{\bm{\Theta}}}={\bm{C}}^{-1}$). We write the precision matrix as

[TABLE]

where ${{\bm{\Theta}}^{\text{null}}}$ is the null precision matrix. Note that we have redefined ${{\bm{W}}}$ by the deviation of the precision matrix from the null precision matrix ${{\bm{\Theta}}^{\text{null}}}$. The penalized log likelihood function (Eq. (7)) is rewritten as

[TABLE]

We note that $\hat{\cal L}$ is a concave function with respect to ${{\bm{W}}}$, which is different from the case of correlation matrices. By exploiting the concavity, we maximise $\hat{\cal L}$ using a gradient descent algorithm instead of the MM algorithm. Specifically, starting from an initial solution $\widetilde{{\bm{W}}}^{(0)}={\bm{0}}$, we update tentative solution $\widetilde{{\bm{W}}}^{(k)}$ until convergence. The update equation is given by

[TABLE]

Each update requires the inversion of an $N\times N$ matrix, resulting in time complexity of the entire algorithm of ${\cal O}(N^{3})$. One may be able to use conventional efficient optimisation algorithms for the graphical Lasso to save time Friedman2008 . However, we do not explore this direction.