Matrix Graphical Models

• Matrix graphical models are statistical frameworks that represent dependencies in matrix or tensor-form data using structured precision or covariance operators across different axes.
• They leverage formulations such as the matrix-normal distribution, Kronecker products and sums, and penalized likelihood methods to encode conditional independence and manage non-independent noise.
• These models have broad applications, including cancer mortality mapping, spatio-temporal data analysis, and multi-modal inference in complex networks.

Matrix graphical models are graphical-statistical models for data that naturally take matrix or tensor form, with dependence represented through structured precision or covariance operators along rows, columns, frequencies, or other axes. In a basic matrix-normal construction, an observed matrix $X\in\mathbb R^{p_R\times p_C}$ satisfies

$$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$

so that separate undirected graphs on rows and columns impose zero constraints in $\Omega_R$ and $\Omega_C$ (Dobra et al., 2010). Related formulations replace a single Kronecker product by a sum of Kronecker products, a Kronecker-sum precision, or frequency-dependent inverse spectral factors, thereby extending the notion of graphical structure to non-independent noise, multimodal tensors, and dependent matrix-valued time series (Dahl et al., 2013, Andrew et al., 2022, Tugnait, 2024). The resulting literature treats “matrix graphical models” as a family of models for conditional independence, sparse estimation, and structured covariance learning rather than as a single estimator.

1. Gaussian matrix-normal foundations

A canonical formulation assumes i.i.d. matrices $X^{(i)}\in\mathbb R^{p\times q}$ with

$$X \sim MN_{p,q}(0,U,V), \qquad \mathrm{vec}(X)\sim N_{pq}(0,V\otimes U).$$

The corresponding row- and column-precision matrices are $\Theta_U=U^{-1}$ and $\Theta_V=V^{-1}$. In the undirected graph on row nodes $\{1,\dots,p\}$, an edge $(a,b)$ is present iff $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$0; the column graph is defined analogously through $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$1 (Shin et al., 7 Jul 2025).

In the matrix-variate Gaussian graphical model of Dobra, Lenkoski and Rodríguez, the row graph $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$2 imposes

$$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$3

which encodes conditional independence of row-vectors $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$4 rest, while the column graph $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$5 imposes

$$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$6

encoding conditional independence of column-vectors $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$7 rest (Dobra et al., 2010). At the entry level, in the matrix-valued Gaussian setting of dependent-data analysis, $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$8 is conditionally independent of $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$9 given all other entries if and only if at least one of $\Omega_R$0 or $\Omega_R$1, so the pair of precision factors jointly determines an undirected graph on the $\Omega_R$2 entries (Tugnait, 2024).

The separable Gaussian model is not the only formulation. In a non-independent-noise model for phenotype data, observations are written as

$\Omega_R$3

with

$\Omega_R$4

so that the total covariance is

$\Omega_R$5

which is not itself a single Kronecker product (Dahl et al., 2013). A different line of work assumes a Kronecker-sum precision,

$\Omega_R$6

and generalizes this construction to multiple axes and multiple modalities that share axes (Andrew et al., 2022). These distinctions are substantive: Kronecker-product covariance, sum-of-Kronecker covariance, and Kronecker-sum precision encode different conditional-independence semantics and lead to different optimization and identifiability constraints.

2. Bayesian inference and model determination

A Bayesian treatment places G-Wishart priors on graph-constrained precision matrices. For a graph $\Omega_R$7 on $\Omega_R$8 vertices, the G-Wishart $\Omega_R$9 has density

$\Omega_C$0

where $\Omega_C$1 is the cone of positive-definite matrices with zeros in all off-diagonal entries $\Omega_C$2 for $\Omega_C$3, and $\Omega_C$4 is finite for $\Omega_C$5 (Dobra et al., 2010). In the matrix-variate model, the prior specification is

$\Omega_C$6

with the identifiability constraint $\Omega_C$7 and auxiliary scalar $\Omega_C$8. Typical default choices are $\Omega_C$9 and $X^{(i)}\in\mathbb R^{p\times q}$0, which place a weakly informative, one-sample equivalent prior.

Posterior computation interleaves five updates: row-graph moves by reversible-jump add/delete proposals, Metropolis–Hastings updates of $X^{(i)}\in\mathbb R^{p\times q}$1 on free Cholesky entries, analogous reversible-jump moves for $X^{(i)}\in\mathbb R^{p\times q}$2, Metropolis–Hastings updates of $X^{(i)}\in\mathbb R^{p\times q}$3 with $X^{(i)}\in\mathbb R^{p\times q}$4 fixed, and a direct gamma draw for the auxiliary scalar $X^{(i)}\in\mathbb R^{p\times q}$5. Data-augmentation appears only in the introduction of $X^{(i)}\in\mathbb R^{p\times q}$6 to break the nonidentifiability of the Kronecker factorization. The same framework is extended to a sparse multivariate CAR model with

$X^{(i)}\in\mathbb R^{p\times q}$7

a Poisson observation model

$X^{(i)}\in\mathbb R^{p\times q}$8

and a sparse unknown graph on $X^{(i)}\in\mathbb R^{p\times q}$9, so that both the spatial autocorrelation $$X \sim MN_{p,q}(0,U,V), \qquad \mathrm{vec}(X)\sim N_{pq}(0,V\otimes U).$$0 and the inter-outcome precision $$X \sim MN_{p,q}(0,U,V), \qquad \mathrm{vec}(X)\sim N_{pq}(0,V\otimes U).$$1 are learned from the data (Dobra et al., 2010).

This Bayesian program was illustrated on simulated and real examples. With $$X \sim MN_{p,q}(0,U,V), \qquad \mathrm{vec}(X)\sim N_{pq}(0,V\otimes U).$$2, $$X \sim MN_{p,q}(0,U,V), \qquad \mathrm{vec}(X)\sim N_{pq}(0,V\otimes U).$$3, true non-decomposable graphs $$X \sim MN_{p,q}(0,U,V), \qquad \mathrm{vec}(X)\sim N_{pq}(0,V\otimes U).$$4, and $$X \sim MN_{p,q}(0,U,V), \qquad \mathrm{vec}(X)\sim N_{pq}(0,V\otimes U).$$5 samples, the MCMC recovered both graphs with high posterior edge-probabilities and accurately estimated the nonzero entries of $$X \sim MN_{p,q}(0,U,V), \qquad \mathrm{vec}(X)\sim N_{pq}(0,V\otimes U).$$6. In cancer mortality mapping for $$X \sim MN_{p,q}(0,U,V), \qquad \mathrm{vec}(X)\sim N_{pq}(0,V\otimes U).$$7 U.S. states and $$X \sim MN_{p,q}(0,U,V), \qquad \mathrm{vec}(X)\sim N_{pq}(0,V\otimes U).$$8 cancer sites, posterior summaries revealed interpretable edge patterns, including a clique among lung, larynx and oral cancers, while the spatial precision departed substantially from the initial ICAR prior. In a bivariate CAR model for verbal and math SAT means across $$X \sim MN_{p,q}(0,U,V), \qquad \mathrm{vec}(X)\sim N_{pq}(0,V\otimes U).$$9 states, the learned inter-score precision had 80% posterior probability of an edge in $\Theta_U=U^{-1}$0. Convergence diagnostics were based on trace-plots and Gelman–Rubin $\Theta_U=U^{-1}$1 statistics on key scalars, and the algorithm mixed in a few thousand iterations even for non-decomposable graphs (Dobra et al., 2010).

3. Penalized likelihood and regression-based graph recovery

A major frequentist direction treats latent signal and noise explicitly. In the model

$\Theta_U=U^{-1}$2

$\Theta_U=U^{-1}$3 is known, $\Theta_U=U^{-1}$4 is the unknown column-precision of interest, and $\Theta_U=U^{-1}$5 is an unknown noise precision (Dahl et al., 2013). With $\Theta_U=U^{-1}$6 latent, the EM-expected objective reduces to minimizing over $\Theta_U=U^{-1}$7, $\Theta_U=U^{-1}$8

$\Theta_U=U^{-1}$9

A common choice is $\Theta_V=V^{-1}$0 and $\Theta_V=V^{-1}$1, so that sparsity is induced only in the signal precision. In the special case $\Theta_V=V^{-1}$2 and $\Theta_V=V^{-1}$3, the M-step updates are

$\Theta_V=V^{-1}$4

Because $\Theta_V=V^{-1}$5 is a sum of two Kronecker terms, one needs a known $\Theta_V=V^{-1}$6 or a constraint on $\Theta_V=V^{-1}$7 versus $\Theta_V=V^{-1}$8 to resolve scale ambiguities. A naïve inversion in the E-step costs $\Theta_V=V^{-1}$9, while eigendecomposing $\{1,\dots,p\}$0 once reduces each E-step to $\{1,\dots,p\}$1; the dominant cost often becomes the Glasso solve, $\{1,\dots,p\}$2 worst-case but very fast in practice for $\{1,\dots,p\}$3. In simulations with $\{1,\dots,p\}$4 individuals in $\{1,\dots,p\}$5 sibships of size $\{1,\dots,p\}$6, $\{1,\dots,p\}$7 traits, and 40 replicates, the full EM method (G$\{1,\dots,p\}$8M) substantially outperformed vanilla Glasso and KronGlasso when $\{1,\dots,p\}$9 was non-iid, while remaining competitive when $(a,b)$0 truly was iid (Dahl et al., 2013).

A separate approach uses neighborhood selection rather than likelihood maximization. For row index $(a,b)$1, the population least-squares coefficient $(a,b)$2 in the regression of row $(a,b)$3 on the remaining rows satisfies

$(a,b)$4

so nonzero entries of $(a,b)$5 exactly mark the nonzero off-diagonals of $(a,b)$6 in the $(a,b)$7-th column up to scaling by $(a,b)$8 (Shin et al., 7 Jul 2025). Treating each of the $(a,b)$9 column-vectors of the observed matrices as independent samples, the estimator solves, for each node,

$$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$00

and then symmetrizes the nodewise neighborhoods by an AND-rule or OR-rule. The primal–dual witness argument gives no false inclusions once a dual certificate satisfies $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$01, and a $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$02-min condition yields no false exclusions. Under the stated regularity assumptions, sufficient sample-size scaling takes the form $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$03 and $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$04, with constants depending on $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$05, eigen-bounds of $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$06, and incoherence $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$07. In simulations on hub, band, and Erdős–Rényi random graphs with $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$08 and $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$09, matrixNS uniformly outperformed GEMINI and Likelihood, especially on banded graphs, and “global” tuning of $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$10 and $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$11 matched the more costly “individual” node-wise tuning (Shin et al., 7 Jul 2025).

4. Dependent observations and multi-axis graphical models

The assumption of i.i.d. matrix observations is restrictive for image sequences, spatio-temporal arrays, and related data. For a stationary matrix-valued Gaussian time series $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$12, with $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$13, the covariance sequence is modeled as

$$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$14

and the power spectral density becomes

$$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$15

Dahlhaus’s criterion implies that an edge in the graph on the $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$16 entries is absent if

$$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$17

Writing $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$18 and $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$19, one has $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$20, so zeros in $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$21 and zeros in $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$22 for all $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$23 define two factor graphs whose Kronecker product is the full conditional-independence graph (Tugnait, 2024).

Estimation is based on a Whittle-approximate negative log-likelihood with sparse-group lasso penalties across frequencies. The optimization is bi-convex in $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$24 and is solved by ADMM within a flip-flop scheme: fixing $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$25 and updating $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$26, then fixing $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$27 and updating $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$28. Each ADMM subproblem consists of a log-det eigen-update, a proximal step, and a dual update. Theoretical guarantees include a population limit, high-dimensional consistency of the $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$29- and $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$30-subproblems in Frobenius norm, and normalized consistency of the flip-flop iterates. The rates are

$$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$31

and

$$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$32

under the stated assumptions on local stationarity, summable autocovariances, sparsity, and uniform eigen-bounds. In synthetic experiments the dependent-data method significantly outperformed the i.i.d. Kronecker-graphical-lasso, and in Beijing air-quality data it recovered a sparse pollution graph and a site graph that clustered rural versus urban stations (Tugnait, 2024).

A further generalization is the Gaussian multi-Graphical Model (GmGM), which handles multiple tensors that share axes. For modality $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$33, the model uses a Kronecker-sum precision $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$34, and the global estimator is

$$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$35

where $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$36 aggregates mode-$$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$37 Gram matrices over all modalities sharing axis $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$38 (Andrew et al., 2022). The key algorithmic result is that if

$$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$39

then the optimizer has the same eigenvectors,

$$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$40

so only a single eigendecomposition per axis is required. For two-axis $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$41 data, GmGM has time $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$42, compared with $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$43 for TeraLasso/EiGLasso. Empirically, GmGM matched or slightly outperformed existing solvers in AUC, completed $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$44 problems in approximately $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$45 minute, and was applied to COIL-20 video, EchoNet-Dynamic echocardiograms, scRNA-seq, LifeLines-DEEP multi-omics, and joint scRNA+scATAC data (Andrew et al., 2022).

5. Structured, interpretable, and algebraic variants

Some matrix graphical models are defined directly as linear spaces of symmetric matrices. In a coloured graphical model, a graph $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$46 has partitions of vertices and edges into colours, and each colour $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$47 defines a $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$48–$$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$49 matrix $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$50. The associated concentration-matrix space is

$$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$51

equivalently specified by entry-equalities $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$52 for equally coloured edges and $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$53 for equally coloured vertices (Davies et al., 2020). A coloured Gaussian graphical model consists of real Gaussians whose concentration lies in $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$54. The reciprocal variety

$$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$55

encodes polynomial relations among covariance entries, and graph symmetries generate linear forms $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$56 that vanish on $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$57. For the uniform-coloured $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$58-cycle, complete graph $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$59, balanced complete bipartite graph $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$60, and $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$61-hyperoctahedral graph $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$62, all linear forms in the ideal come from symmetries. By contrast, a $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$63-cycle with one distinguished vertex and unbalanced $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$64 require additional relations beyond symmetry-induced ones (Davies et al., 2020).

Another structured direction focuses on correlation matrices rather than precision matrices. For a decomposable graph $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$65, latent “parent” nodes are introduced so that the marginal covariance of the observed children is

$$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$66

where $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$67 is the sparse precision on the augmented node set (Sterrantino et al., 2023). Parent nodes have variance parameters $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$68, and in the toy graph of Section 2.2 the induced child correlations are explicit rational functions such as

$$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$69

Penalized-complexity priors shrink each $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$70 toward the simpler base model $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$71, with distance

$$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$72

Posterior inference can be performed by MCMC, MAP, or INLA. In a four-cancer disease-mapping example on $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$73 German districts, a single PC-prior rate $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$74 yielded posterior correlations concordant with those seen empirically, and DIC, WAIC, and CPO showed a fit as good as the fully unstructured MCAR but at half the hyperparameter cost (Sterrantino et al., 2023).

A separate but related formalism introduces a matrix algebra for graphical statistical models on directed mixed graphs. The basic objects are matrices of sets of walks, one for directed edges and one for bidirected edges, and the algebra is closed under set-union, walk-concatenation, and transpose (Zhao, 2024). In this framework, latent projection corresponds to marginalization, $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$75-separation implies conditional independence in every non-singular Gaussian system on the graph, and the covariance matrix of a Gaussian linear system satisfies the trek rule

$$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$76

The same algebra is used to represent confounder adjustment and the augmentation criterion. This work is not a matrix-variate Gaussian model in the separable-covariance sense, but it provides a matrix-based language for graphical operations that appear throughout the broader theory of graphical models (Zhao, 2024).

6. Extremal models, applications, and recurring issues

Matrix graphical structure also appears outside the Gaussian domain. In the Hüsler–Reiss family for multivariate extremes, dependence is parameterized by a variogram matrix $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$77 with zero diagonal and conditional negative definiteness on $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$78. The associated Hüsler–Reiss precision matrix $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$79 is positive semi-definite, has rank $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$80, satisfies $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$81, and obeys

$$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$82

Equivalent relations include

$$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$83

with $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$84 (Hentschel et al., 2022). For any connected graph $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$85, there exists a unique completion of a partially specified variogram matrix such that the resulting $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$86 has zeros on the non-edges of $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$87. If $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$88 is known, consistent edge-wise variogram estimates $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$89 are completed by $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$90; if $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$91 is unknown, the procedure can be combined with extreme-MST or a penalized surrogate-likelihood on $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$92 (EGlearn). Applications to U.S. flight delays and Danube river flows showed that data-driven sparse graphs outperformed physical-network graphs in held-out likelihood (Hentschel et al., 2022).

Across the literature, several methodological issues recur. Identifiability is handled by explicit constraints: $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$93 with an auxiliary scalar $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$94 in Bayesian Kronecker-factor models, a known $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$95 or an alternative normalization in sum-of-Kronecker covariance models, gradient projection for diagonal non-identifiability in Kronecker sums, and Frobenius-norm normalization for dependent spectral factors (Dobra et al., 2010, Dahl et al., 2013, Andrew et al., 2022, Tugnait, 2024). Exact support-recovery results are tied to stated regularity conditions such as irrepresentability, degree bounds, eigen-bounds, incoherence, sparsity, or $$\mathrm{vec}(X^T)\mid \Omega_R,\Omega_C \sim N_{p_Rp_C}\bigl(0,(\Omega_R\otimes\Omega_C)^{-1}\bigr),$$96-min assumptions (Shin et al., 7 Jul 2025), while computational tractability depends on exploiting sparsity, clique structure, eigendecompositions, or local Cholesky updates (Dobra et al., 2010, Hentschel et al., 2022, Andrew et al., 2022).

The application range is correspondingly broad. Reported case studies include cancer mortality surveillance and SAT score regression under sparse MCAR priors (Dobra et al., 2010), multi-trait genetic studies with non-independent noise (Dahl et al., 2013), EEG recordings with simultaneous spatial and temporal graph estimation (Shin et al., 7 Jul 2025), Beijing air-quality monitoring under dependent matrix-valued time series (Tugnait, 2024), and multimodal single-cell, imaging, and video data under shared-axis tensor models (Andrew et al., 2022). Taken together, these developments suggest that the central theme of matrix graphical modelling is the use of matrix-structured operators—precision matrices, variograms, reciprocal varieties, or walk matrices—to encode sparse dependence in ways that remain computationally manageable and scientifically interpretable.