{'fraction_non_alphanumeric': 0.04874720044792833, 'fraction_numerical': 0.017882138857782753, 'mean_word_length': 4.518926226342217, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Learning Graphical Factor Models with Riemannian Optimization

Graphical models and factor analysis are well-established tools in multivariate statistics. While these models can be both linked to structures exhibited by covariance and precision matrices, they are generally not jointly leveraged within graph learning processes. This paper therefore addresses this issue by proposing a flexible algorithmic framework for graph learning under low-rank structural constraints on the covariance matrix. The problem is expressed as penalized maximum likelihood estimation of an elliptical distribution (a generalization of Gaussian graphical models to possibly heavy-tailed distributions), where the covariance matrix is optionally constrained to be structured as low-rank plus diagonal (low-rank factor model). The resolution of this class of problems is then tackled with Riemannian optimization, where we leverage geometries of positive definite matrices and positive semi-definite matrices of fixed rank that are well suited to elliptical models. Numerical experiments on real-world data sets illustrate the effectiveness of the proposed approach.

Within Gaussian graphical models (GGMs), graph learning boils down to the problem of estimating the precision (inverse covariance) matrix of a Gaussian Markov random field (Friedman et al., 2008;Lake and Tenenbaum, 2010). In practice, achieving an accurate covariance matrix estimation is often a difficult task due to low sample support. Thus, it is common to introduce prior assumptions on the structure of this matrix that guarantee a correct estimation with fewer samples. A popular approach is related to lowrank factorizations, which relies on the assumption that the data is driven by an underlying low-dimensional linear model, corrupted by an independent perturbation. The resulting covariance matrix decomposition then involves a core that is a low-rank positive semi-definite matrix. Such model is ubiquitous in statistics, and for example, at the heart of probabilistic principal component analysis (Tipping and Bishop, 1999), low-rank factor analysis (Robertson and Symons, 2007;Khamaru and Mazumder, 2019;Rubin and Thayer, 1982), and their many generalizations.
Since GGMs and low-rank factorizations share a common root in structured covariance (or precision) matrix estimation, it appears desirable to leverage both approaches in a unified graph learning formulation. On one hand, linear dimension reduction approaches rely on particular spectral structures that can be beneficial for graph learning (Kumar et al., 2020). On the other hand, it also opens the way to graph-oriented view of sparse principal component analysis (Yoshida and West, 2010;Meng et al., 2014). Though theoretically appealing, such unification is challenging because it formulates optimization problems with objective functions and constraints that apply both on the covariance matrix and its inverse. Thus, deriving single-step learning algorithms for these models has only recently been addressed (Chandra et al., 2021).
In this paper, we propose a new family of methods for graph learning with low-rank constraints on the covariance matrix, hereafter referred to as graphical factor models arXiv:2210.11950v1 [stat.ML] 21 Oct 2022 (GFM). First, we reformulate graph learning as a problem that encompasses both elliptical distributions and low-rank factor models. The main interest of generalizing Gaussian graphical models to elliptical ones is to ensure robustness to underlying heavy-tailed distributions (Vogel and Fried, 2011;Finegold and Drton, 2014;Zhang et al., 2013;de Miranda Cardoso et al., 2021). Moreover, additionally considering low-rank factor models allows for an effective dimensionality reduction. The main novelty of our approach is to tackle the resulting class of constrained and penalized maximum likelihood estimation in a unified way with Riemannian optimization (Absil et al., 2008;Boumal, 2020). To do so, we leverage geometries of of both the positive definite matrices (Bhatia, 2009), and positive semi-definite matrices of fixed rank (Bonnabel and Sepulchre, 2009;Bouchard et al., 2021) that are well suited to the considered models. The corresponding tools allows us to develop optimization methods that ensure the desired structures for the covariance matrix, thus providing a flexible framework for learning graphical factor models.
Finally, experiments conducted on real-world data sets demonstrate the interest of considering both elliptical distributions and factor model structures in a graph learning process. We observe that the proposed algorithms offer better interpretability of the obtained graphs compared to unstructured models. Notably, the factor model approaches compare well with the current state-of-the-art on Laplacian-constrained graph learning method that can set the number of components as additional prior information (Kumar et al., 2020;de Miranda Cardoso et al., 2021).
2 Background and proposed framework 2.1 Gaussian graphical and related models Gaussian graphical models assume that the random vector For the corresponding Gaussian Markov random field, an undirected graph is matched to the variables as follows: each variable corresponds to a vertex, and an edge is present between two vertices if the corresponding random variables are conditionally dependent given the others (Dempster, 1972;Lauritzen, 1996). The support of the precision matrix Θ = Σ −1 directly accounts for this conditional dependency, since Hence, a non-zero entry Θ q implies a conditional dependency between the variables x q and x , underlined by an edge between vertices q and of the graph. Within Gaussian graphical models, graph learning is therefore tied to the problem of estimating the precision matrix Θ from a set of observations {x i } n i=1 ∈ (R p ) n . In order to exhibit such correlation structure, a standard approach is to resort to regularized maximum likelihood estimation, i.e., solving the problem: where S = 1 n n i=1 x i x i is the sample covariance matrix, h is a regularization penalty, and λ ∈ R + is a regularization parameter. The 1 -norm if often used as penalty in order to promote a sparse structure in Θ, which is at the foundation of the well-known GLasso algorithm (Friedman et al., 2008;Lake and Tenenbaum, 2010;Mazumder and Hastie, 2012;Fattahi and Sojoudi, 2019). Many other convex on non-convex penalties have been considered in this stetting (Lam and Fan, 2009;Shen et al., 2012;Benfenati et al., 2020). Depending on the assumptions, it can also be beneficial to consider penalties that promote certain types of structured sparsity patterns (Heinävaara et al., 2016;Tarzanagh and Michailidis, 2018). Another main example of structure within the conditional correlations is total positivity, also known as attractive Gaussian graphical models, that assumes Θ q ≤ 0, ∀q = (Fallat et al., 2017;Lauritzen et al., 2019). In attractive Gaussian graphical models, the identifiability of the precision matrix to the graph Laplacian (Chung, 1997), has also attracted recent interest from the point of view of graph learning (Egilmez et al., 2017;Ying et al., 2020;Kumar et al., 2020).

Elliptical distributions
A first shortcoming of graph learning as formulated in (2) is its lack of robustness to outliers, or heavy-tailed distributed samples. This is consequence of the underlying Gaussian assumption, that cannot efficiently describe such data. A possible remedy is to consider a larger family of multivariate distributions. In this context, elliptical distributions (Anderson and Fang, 1990;Kai-Tai and Yao-Ting, 1990) offer a good alternative, that are well-known to provide robust estimators of the covariance matrix (Maronna, 1976;Tyler, 1987;Wald et al., 2019). In our present context, this framework has been successfully used to extend graphical models (Vogel and Fried, 2011;Finegold and Drton, 2014;Zhang et al., 2013;de Miranda Cardoso et al., 2021), as well as low-rank structured covariance matrices models (Zhao and Jiang, 2006;Zhou et al., 2019;Bouchard et al., 2021) for corrupted or heavy-tailed data.
A vector is said to follow a centered elliptically symmetric distribution of scatter matrix Σ and density generator g, denoted x ∼ ES(0, Σ, g), if its density has the form which yields the negative log-likelihood . Notice that using g(t) = exp(−t/2) allow us to recover the Gaussian case. However, the density generator g enables more flexibility, and notably to encompass many heavy-tailed multivariate distributions. Among popular choices, elliptical distributions include the multivariate tdistribution with degree of freedom ν > 2, which is obtained with g(t) = (1 + t/ν) − ν+p 2 . For this distribution, the parameter ν measures the rate of decay of the tails. Also notice that the Gaussian case corresponds to the limit case ν → ∞.

Low-rank factor models
A second limitation of (2) is that it does not account for other potential structure exhibited by the covariance or precision matrices. This can problematic when the sample support is low (n p, or n < p) as the sample covariance matrix is not a reliable estimate in these setups (Ledoit and Wolf, 2004;Smith, 2005;Vershynin, 2012). In this context, a popular approach consists in imposing low-rank structures on the covariance matrix. These structures arise from the assumption that the data is driven by an underlying low-dimensional linear model, i.e., x = W s + , where W is a rank-k factor loading matrix, and s ∈ R k are ∈ R p are independent random variables. Thus the resulting covariance matrix is of the form Σ ∆ = H + Ψ, where H = W E[ss ]W belongs to the set of positive semi-definite matrices of rank k, denoted S + p,k = {Σ ∈ S p , Σ 0, rank(Σ) = k}. This model is well known in statistics, and for example, at the core of probabilistic principal component analysis, that assumes Ψ ∝ I p (Tipping and Bishop, 1999). Also notice that in Laplacian-constrained models, a rank-k precision matrix implies a (p−k)-component graph (Chung, 1997). Hence it is also interesting to leverage such spectral structures from the graph learning perspective (Kumar et al., 2020).
In this paper, we will focus on the low-rank factor analysis models (Robertson and Symons, 2007;Rubin and Thayer, 1982;Khamaru and Mazumder, 2019), that consider the more general case Ψ ∈ D ++ p , where D ++ denotes the space of positive definite diagonal matrices. Given this model, the covariance matrix belongs to the space of rank-k plus diagonal matrices, denoted as Notice that this parameterization reduces the dimension of the estimation problem reduces from p(p + 1)/2 to p(k + 1) − k(k − 3)/2, which is why it is often used in regimes with few samples, or high dimensional settings.

Learning elliptical graphical factor models
We cast the problem of graph learning for elliptical graphical factor models as where L is the negative-log likelihood in (4), and M p,k is the space of rank-k plus diagonal positive definite matrices in (5). The penalty h is a smooth function that promotes a sparse structure on the graph, and λ ∈ R + is a regularization parameter. Although this formulation take into account many options, we focus on the usual element-wise penalty applied to the precision matrix Θ = Σ −1 , defined as: It is important to notice that the considered optimization framework will require h to be smooth. This is why we presently use a surrogate of the 1 -norm defined as with ε > 0. The 1 -norm can be recovered with lim ε→0 φ(t). The problem in (6) thus accounts for a lowrank structure in a robust graph learning formulation, and is also expected to be more robust to data following an underlying heavy-tailed distributions. We then introduce four main cases, and their corresponding acronyms: GGM/GGFM: Gaussian graphical factor models (GGFM) are obtained with the Gaussian log-likelihood, i.e., setting g(t) = exp(−t/2) in (4). The standard Gaussian graphical models (GGM) as in (2) are recovered when dropping the constraint Σ ∈ M p,k .
EGM/EGFM: elliptical graphical factor models (EGFM) are obtained for the more general case where L defines a negative log-likelihood of an elliptical distribution. Dropping the constraint Σ ∈ M p,k also yields a relaxed problem that we refer to as elliptical graphical models (EGM).

Learning graphs with Riemmanian optimization
The optimization problem in (6) is non-convex and difficult to address. Indeed, it involve objective functions and constraints that apply on both the covariance matrix and its inverse. An approach to handle such multiple intricate constraints is to resort to variable splitting and alternating direction method of multipliers (Kovnatsky et al., 2016) which has, for example, been considered for Laplacianconstrained models in (Zhao et al., 2019;de Miranda Cardoso et al., 2021). In this work, we propose a new and more direct approach by harnessing Riemannian optimization (Absil et al., 2008;Boumal, 2020).
For GGM and EGM, the resolution of (6) requires to consider optimization on S ++ p . This will be presented in Section 3.1, which will serve as both a short introduction to optimization on smooth Riemannian manifold, and as a building block for the resolution of (6) for GGFM and EGFM. For these models, we will derive Riemannian optimization algorithms on Σ ∈ M p,k in Section 3.2, which will be done by leveraging a suitable geometry for this space.
3.1 Learning GGM/EGM: optimization on S ++ p When relaxing the constraint in (6), the problem still requires to be solved on S ++ can simply be identified as S p , for any Σ ∈ S ++ p . To be able to perform optimization, the first step is to define a Riemannian metric, i.e., an inner product on this tangent space. In the case of S ++ p , the most natural choice is the affine-invariant metric (Skovgaard, 1984;Bhatia, 2009), that corresponds to the Fisher information metric of the Gaussian distribution and features very interesting properties from a geometrical point of view. It is defined for all Σ ∈ S ++ p , ξ, η ∈ S p as ξ, η The Riemannian gradient (Absil et al., 2008) at Σ ∈ S ++ p of an objective function f : S ++ p → R is the only tangent vector such that for all ξ ∈ S p , where D f (Σ)[ξ] denotes the directional derivative of f at Σ in the direction ξ. The Riemannian gradient ∇ S ++ p f (Σ) can also be obtained from the Euclidean gradient ∇ E f (Σ) through the formula where sym(·) returns the symmetrical part of its argument. The Riemannian gradient is sufficient in order to define a descent direction of f , i.e., a tangent vector inducing a decrease of f , hence to define a proper Riemannian steepest descent algorithm. However, if one wants to resort to a more sophisticated optimization algorithm such as conjugate gradient or BFGS, vector transport (Absil et al., 2008) is required to transport a tangent vector from one tangent space to another. In the case of S ++ p , we can employ the most natural one, i.e., the one corresponding to the parallel transport associated with the affine-invariant metric (Jeuris et al., 2012). The transport of the tangent vector ξ of Σ onto the tangent space of Σ is given by Now that we have all the tools needed to obtain a proper descent direction of some objective function f , it remains to be able to get from the tangent space back onto the manifold S ++ p . This is achieved by means of a retraction map. The best solution on S ++ p in order to ensure numerical stability while taking into account the chosen geometry is This retraction corresponds to a second order approximation of the geodesics of S ++ p (Jeuris et al., 2012), which generalize the concept of straight lines for a manifold. We now have all the necessary elements to perform the optimization of f on S ++ p . For instance, the sequence of iterates {Σ t } and descent directions {ξ t } generated by a Riemannian conjugate gradient algorithm is where α t is a stepsize that can be computed through a linesearch (Absil et al., 2008) and β t can be computed using the rule in (Hestenes and Stiefel, 1952) for example.
From there, to obtain a specific algorithm that solves (6) on S ++ p , it only remains to provide the Riemannian gradients of the negative log-likelihood L and penalty h.
where u(t) = −2g (t)/g(t). For the Gaussian distribution, we have u(t) = 1. For the t-distribution with degree of freedom ν, we have u(t) = (ν + p)/(ν + t). Concerning h, the q element of its Riemannian gradient at Σ ∈ S ++ where φ (t) = tanh(t/ε). In the next section, the Euclidean gradients of L and h in S ++ p are used. They are easily obtained from (11) and by noticing that the Euclidean gradient is already symmetric.

Learning GGFM/EGFM: optimization on M p,k
When considering the factor model, we aim at finding a structured matrix Σ in M p,k , which is a smooth submanifold of S ++ p . The Riemannian geometry of M p,k is not straightforward. In fact, defining an adequate geometry for low-rank matrices is a rather difficult task and, while there were many attempts, see e.g., (Bonnabel and Sepulchre, 2009;Meyer et al., 2011;Vandereycken et al., 2012;Massart and Absil, 2018;Neuman et al., 2021;Bouchard et al., 2021), no perfect solution has been found yet. Here, we choose the parametrization considered in (Bonnabel and Sepulchre, 2009;Meyer et al., 2011;Bouchard et al., 2021).
From a geometrical perspective, the novelty here is to adopt the particular structure of M p,k , i.e., by combining a lowrank positive semi-definite matrix and a diagonal positive definite matrix. Each of the two manifolds of this product are well-known and the proofs of the following can easily be deduced from the proofs in (Bouchard et al., 2021). All Σ ∈ M p,k can be written as where It follows that M p,k = ϕ(N p,k ). Therefore, to solve (6) on M p,k , one can exploit ϕ and solve it on N p,k . However, N p,k contains invariance classes with respect to ϕ, that is for any O ∈ O k (orthogonal group with k × k matrices), As a consequence, the space that best corresponds to M p,k is the quotient manifold N p,k /O k induced by equivalence classes on N p,k where . To efficiently solve an optimization problem on M p,k with the chosen parametrization, we thus need to consider the quotient N p,k /O k . In practice, rather than dealing with the quite abstract manifold N p,k /O k directly, we will manipulate objects in N p,k , implying that the different points of an equivalence class (20) are, in fact, one and only point.
Our first task is to provide the tangent space of θ = (V , Λ, Ψ) ∈ N p,k . It is obtained by aggregating tangent spaces T V St p,k , T Λ S ++ k and T Ψ D ++ p , i.e., Then, we define a Riemannian metric on N p,k . It needs to be invariant along equivalence classes, i.e., for all θ ∈ N p,k , ξ, η ∈ T θ N p,k and O ∈ O k ξ, η In this work, the Riemannian metric N p,k is chosen as the sum of metrics on St p,k , S ++ k and D ++ p . The metric on St p,k is the so-called canonical metric (Edelman et al., 1998), which yields the simplest geometry of St p,k . The metrics on S ++ k and D ++ p are the affine-invariant ones. Hence, our chosen metric on N p,k is At θ ∈ N p,k , T θ N p,k contains tangent vectors inducing a move along the equivalence class π(θ). In our setting, these directions are to be eliminated. They are contained in the so-called vertical space V θ , which is the tangent space T θ π(θ) to the equivalence class. In our case, it is where S ⊥ k denotes the vector space of skew-symmetric matrices. Now that the unwanted vectors have been identified, we can deduce the ones of interest: they are contained in the orthogonal complement of V θ according to metric (23). This space is called the horizontal space H θ and is equal to (25) The Riemannian gradient of an objective function f on N p,k is defined through the chosen metric as for S ++ p in (10). Again, it is possible to obtain the Riemannian gradient ∇ N p,k f (θ) on N p,k from the Euclidean one ∇ E f (θ). This is achieved by (26) where ∇ E f (θ) = (G V , G Λ , G Ψ ) and ddiag(·) cancels the off-diagonal elements of its argument. Furthermore, since the objective functions we are interested in are invariant along equivalence classes (20), their Riemannian gradient is naturally on the horizontal space H θ and we can use it directly. In this work, our focus is objective functions f : N p,k → R such that f = f • ϕ, where f : S ++ p → R; see (6). Hence, it is of great interest to be able to get the gradient of f directly from the one of f . The Euclidean gradient ∇ E f (θ) as a function of the Euclidean gradient While it is not necessary for the gradient, the vector transport requires to be able to project elements from the ambient space R p×k × R k×k × R p×p onto the horizontal space. The first step is to provide the orthogonal projection from the ambient space onto T θ N p,k . For all θ ∈ N p,k and From there we can obtain the orthogonal projection from T θ N p,k onto H θ . Given ξ ∈ T θ N p,k , it is where Ω ∈ S ⊥ k is the unique solution to We can now define an adequate vector transport operator. Given θ, θ ∈ N p,k and ξ ∈ H θ , it is simply The last object that we need is a retraction on N p,k . To address the invariance requirement, the retraction needs to be invariant along equivalence classes, i.e., for all θ ∈ N p,k , In this paper, we choose a second order approximation of geodesics on N p,k given by where uf(·) returns the orthogonal factor of the polar decomposition.
We finally have all the tools needed to solve problem (6) on N p,k and thus to ensure the structure of M p,k . With all the objects given in this section, many Riemannian optimization algorithms can be employed to achieve the minimization (6), e.g., steepest descent, conjugate gradient or BFGS. In this work, we use in practice a conjugate gradient algorithm.

Experiments
In this section, the performances of the proposed algorithms (GGM/GGFM and EGM/EGFM) are evaluated on three real-world data sets consisting of (i) animal species, (ii) GNSS time series and (iii) concepts. The purpose of the following experiments is to verify whether the factor model, which accounts for potential structures in the covariance matrix, provides an improved version of the learned graph topology. The latter is evaluated using both graph modularity m, for which high values measure high separability of a graph into sub-components (Newman, 2006), and visual inspection 1 .
Benchmarks: The proposed algorithms are compared with state-of-the-art approaches for (i) connected graphs: GLasso (Friedman et al., 2008), which uses 1 -norm as a sparse-promoting penalty, and NGL , which use a Gaussian Laplacian-constrained model with a concave penalty regularization; (ii) multi-component graphs: SGL (Kumar et al., 2020), which constraints the number of graph components in an NGL-type formulation, and StGL (de Miranda Cardoso et al., 2021), which generalizes the above to t-distributed data.
The proposed methods are tuned as follows: for GGM/GGFM, the initial guess is I p ; for EGM/EGFM, it is fixed to V I k V +I p where V ∈ St p,k is an orthogonal matrix obtained from the EVD of 1 n i x i x i . For these methods, we then consider an edge to be active when the conditional correlation in (1) is positive and over a tolerance threshold of 10 −2 . For a fair comparison, the parameters of all tested algorithms (including competing methods) are tuned to display the best results. We also notice that all these methods are not very sensitive to a slight change of the parameters, e.g., the rank k, or the order of magnitude of regularization parameter λ.
Animals data: In the animals data set (Osherson et al., 1991;Lake and Tenenbaum, 2010) each node represents an animal, and is associated with binary features (categorical non-Gaussian data) consisting of answers to questions such as "has teeth?", "is poisonous?", etc. In total, there are p = 33 unique animals and n = 102 questions. Figure 1 displays the learned graphs of the animals data set. Following the recommendations of (Kumar et al., 2020;Egilmez et al., 2017), the SGL algorithm is applied with β = 0.5 and α = 0, GLasso with α = 0.05 and the input for SGL and NGL is set to 1 n i x i x i + 1 3 I p . While EGM yields results similar to GLasso, GGM provides a clearer structure with (trout, salmon) and (bee, butterfly) clustered together. Interestingly, with no assumptions regarding the number of components, GGFM and EGFM reach similar structure and modularity compared to SGL and StGL, that set a priori the number of components in the graph.

GNSS data:
We present an application on Earth surface displacement data collected by a network of Global Navigation Satellite System (GNSS) receivers from the volcanological observatory of Piton de la Fournaise (OVPF-IPGP), located at La Réunion Island. The presented network is composed of p = 22 receivers recording daily vertical displacements from January 2014 to March 2017 (Smittarello et al., 2019), with a total of n = 1106 observations. During this time, vertical displacements induced by volcano eruptions have been recorded by the receivers network, sometimes leading to abrupt motion changes over time. Depending on their spatial position, some receivers might move in a particular direction (upwards or downwards), thus indicating (thin) spatial correlations.
Results of the learned graphs are presented in Figure 2. A general observation is that all graphs are mainly clustered into three groups: two located West (receivers 'GITG', 'FREG', etc.) and East ('BORG', 'FJAG', etc.) of the summit crater, and one extending from lower altitudes to the seashore ('CASG', 'TRCG', etc.). We call these groups west, east and low components, respectively. As described in (Peltier et al., 2017), the four 2015 eruptions (February, May, July and August) are characterized by asymmetric displacement patterns with respect to the North-South eruptive fissures which extend slightly westward from the summit crater. Interestingly, this corresponds to the separation between west and east graph components, which is best evidenced by factor model-based algorithms, especially EGFM, since it can be more robust to abrupt changes (heavy tails in the data distribution). The low graph component corresponds to receivers with small or no displacement. Note that the 'PRAG' receiver is also included in this group, probably because it did not record a significant motion during this period. Finally, GGM, EGM and StGL lead to similar results but with spurious connections from the west side of the crater to seashore receivers (e.g., from 'BORG' to 'CASG' for StGL, and 'PRAG' connected to some west and/or east receivers).
Concepts data: The concepts data set (Lake and Tenenbaum, 2010) collected by Intel Labs, includes p = 1000 nodes 2 and n = 218 semantic features. As stated in (Cai et al., 2022), "each node denotes a concept such as 'house', 'coat', and 'whale', and each semantic feature is a question such as 'Can it fly?', 'Is it alive?', and 'Can you use it?' The answers are on a five-point scale from 'definitely no' to 'definitely yes', conducted on Amazon Mechanical Turk". This dataset will only be used to assess the suitability of our factor model approaches, which are expected to be more efficient for such high-dimensional setting. Since GGFM and EGFM lead to very similar results, only the former is shown. In comparison, GGM and EGM did not provide as clean interpretable results. Because the "true" number of components is unclear in this data, methods that require prior setting of this variable were also found hard to exploit. The graph learned by the GGFM algorithm is presented in Figure 3a. This graph is mostly composed of connected sub-components, each falling into interpretable categories of similar concepts ("tools" in grey, "body parts" in blue, etc.). A closer look on the graph in Figure 3b shows interesting nodes that can be interpreted as links between the clusters of concepts (e.g., the node 'body' linking the cluster of "body parts" to the one of "animals"). Another interpretation is that concepts with ambivalent meaning also act as nodes between these clusters ('gas' can be associated to "drinks" as well as to "natural elements").
In summary, we have proposed a family of graph learning algorithms in a unified formulation involving elliptical distributions, low-rank factor models, and Riemannian optimization. Our experiments demonstrated that the proposed approaches can evidence interpretable correlation structures in exploratory data analysis for various applications.