theoremTheorem lemmaLemma corollaryCorollary

# Making Multi-Axis Gaussian Graphical Models Scalable to Millions of Samples
and Features

Bailey Andrew<EMAIL_ADDRESS>
School of Computing
University of Leeds
Leeds, LS2 9JT, UK David R. Westhead<EMAIL_ADDRESS>
Faculty of Biological Sciences
University of Leeds
Leeds, LS2 9JT, UK Luisa Cutillo<EMAIL_ADDRESS>
School of Mathematics
University of Leeds
Leeds, LS2 9JT, UK

###### Abstract

Gaussian graphical models can be used to extract conditional dependencies
between the features of the dataset. This is often done by making an
independence assumption about the samples, but this assumption is rarely
satisfied in reality. However, state-of-the-art approaches that avoid this
assumption are not scalable, with $O(n^{3})$ runtime and $O(n^{2})$ space
complexity. In this paper, we introduce a method that has $O(n^{2})$ runtime
and $O(n)$ space complexity, without assuming independence.

We validate our model on both synthetic and real-world datasets, showing that
our method’s accuracy is comparable to that of prior work. We demonstrate that
our approach can be used on unprecedentedly large datasets, such as a real-
world 1,000,000-cell scRNA-seq dataset; this was impossible with previous
approaches. Our method maintains the flexibility of prior work, such as the
ability to handle multi-modal tensor-variate datasets and the ability to work
with data of arbitrary marginal distributions. An additional advantage of our
method is that, unlike prior work, our hyperparameters are easily
interpretable.

Keywords: gaussian graphical models, multi-axis models, transcriptomics,
multi-omics, scalability

##  1 Introduction

Figure 1: Dependencies between (gene, cell) pairs. In the same cell, two genes
are connected only if the genes are related to each other. An analogous
statement holds for cell-cell relations. Cross terms, such as (A, X) $\sim$
(B, Y), may exist - but they are indirect effects. Conditioning out the
mediators, in this case (B, X) and (A, Y), will remove the relation.

It is often the case that we want to find networks (‘graphs’) that describe
our data, whether it be for gene regulatory networks, tumor microenvironment
structure, social networks, or as a preprocessing step for a clustering
algorithm. There are several types of graphs, such as nearest-neighbors or
correlation relationships, but conditional dependency graphs are especially
attractive due to their sparsity, interpretability, and relationship to
causality. Most methods to find such graphs make independence assumptions, and
those that don’t are not scalable to millions of samples and features.

This paper focuses on the challenge of avoiding independence assumptions while
preserving scalability. To demonstrate the need for such methods, we will
focus on the case of single-cell transcriptomics (scRNA-seq). scRNA-seq can be
represented with a matrix whose rows correspond to individual cells and whose
columns correspond to genes. The number of cells in such datasets is
increasing exponentially as sequencing techniques improve: there is a clear
need for scalability. Likewise, neither genes nor cells could be reasonably
considered to be independent; genes regulate one another, and cells interact
with each other.

Conditional dependence can be defined as a measure of statistical dependence
once one conditions over other features in the dataset. Formally, if
$\mathbf{D}$ is the set of features in a dataset, and $\mathbf{F}$ is any
subset that does not contain $x$ and $y$, we say that $x$ is conditionally
independent of $y$ (conditioning over $\mathbf{F}$) if $\mathbb{P}\left(x\cap
y|\mathbf{F}\right)=\mathbb{P}\left(x|\mathbf{F}\right)\mathbb{P}\left(y|\mathbf{F}\right)$.

The graph of conditional dependencies connects $x$ and $y$ if and only if they
are not conditionally independent. Intuitively, this idea captures the ‘direct
effect’ of one variable on another, as opposed to ‘indirect effects’ that pass
through other variables. Due to this, the graph of conditional dependencies
will be much sparser than a graph of correlations.

In this paper, we will focus on the case where
$\mathbf{F}=\mathbf{D}_{\backslash x,\backslash y}$ ($\mathbf{D}$ with $x$ and
$y$ removed). For normally distributed data, there is a nice relationship
between the conditional dependencies and the inverse covariance matrix, also
known as the precision matrix $\mathbf{\Psi}$:
$\mathbf{\Sigma}^{-1}=\mathbf{\Psi}$. $x$ and $y$ are conditionally
independent, given the rest of the dataset, if and only if the $(x,y)$th entry
of $\mathbf{\Psi}$ is 0.

$\mathbf{\Psi}_{xy}=0\iff\mathbb{P}\left(x\cap y\text{ }|\text{
}\mathbf{D}_{\backslash x\backslash
y}\right)=\mathbb{P}\left(x|\mathbf{D}_{\backslash x\backslash
y}\right)\mathbb{P}\left(y|\mathbf{D}_{\backslash x\backslash y}\right)$

If we denote our dataset as $\mathbf{D}_{n\times m}$, we could model it with a
normal distribution:
$\mathbf{D}_{i}\sim_{i.i.d.}\mathcal{N}\left(\boldsymbol{\mu},\mathbf{\Psi}^{-1}_{m\times
m}\right)$. It is well known how to solve for $\mathbf{\Psi}$ in this case.

However, in practice the $i.i.d.$ assumption may not hold. We could vectorize
it such that $\mathrm{vec}\left[\mathbf{D}_{n\times
m}\right]\sim\mathcal{N}\left(\boldsymbol{\mu},\mathbf{\Psi}^{-1}_{nm\times
nm}\right)$. This formulation makes absolutely no independence assumptions;
any element in $\mathbf{D}$ can be connected to any other element. This leads
to a precision matrix of $O(n^{2}m^{2})$ parameters. Not only is this
computationally intractible to estimate (a $100\times 100$ double-precision
dataset would require almost a gigabyte of space), but statistically
intractable as well.

However, not all of those element-to-element connections are needed. Does the
expression of gene $A$ in cell $X$ directly affect the expression of gene $B$
in cell $Y$? It is more plausible that it is an indirect affect, mediated by
the expression of gene $B$ in cell $X$. This situation is shown in Figure 1.

Based on this observation, the dependencies can be explained in terms of the
dependencies between genes and the dependencies between cells, rather than the
complicated network of dependencies between (gene, cell) pairs. Put another
way, we are assuming that there are two factor matrices $\mathbf{\Psi}_{gene}$
and $\mathbf{\Psi}_{cell}$ and some method $\zeta$ of combining them to get
the full $O(n^{2}m^{2})$ precision matrix:
$\mathbf{\Psi}=\zeta\left(\mathbf{\Psi}_{cell},\mathbf{\Psi}_{gene}\right)$.
For tensor-variate datasets, we can consider more than two factors, one for
each axis of the data:
$\mathbf{\Psi}=\zeta\left(\mathbf{\Psi}_{1},...,\mathbf{\Psi}_{K}\right)$. The
challenge is then to estimate the factor matrices
$\left\\{\mathbf{\Psi}_{\ell}\right\\}_{\ell}$ without constructing the full
matrix. An overview of choices of $\zeta$ is given in Section 1.1; this paper
considers the Cartesian product. In this paper, we will restrict the factor
matrices $\mathbf{\Psi}_{\ell}$ to be positive semidefinite; this allows
$\mathbf{\Psi}_{\ell}$ to maintain their interpretation as axis-wise precision
matrices alongside their interpretation as axis-wise graphs.

###  1.1 Related Work

A popular method of graph inference with the independence assumption is the
Graphical Lasso (Friedman et al., 2008). It addresses the following
optimization problem, where $\mathbf{S}$ is the sample covariance matrix:

$\displaystyle\hat{\mathbf{\Psi}}$
$\displaystyle=\mathrm{argmin}_{\mathbf{\Psi}\succ
0}\mathrm{tr}\left[\mathbf{\Psi}\mathbf{S}\right]-\mathrm{logdet}\mathbf{\Psi}-\rho\lVert\mathbf{\Psi}\rVert_{1,od}$

Here, $\lVert\cdot\rVert_{1,od}$ represents an L1 penalty on the off-diagonal
terms. The Graphical Lasso penalty corresponds to L1-penalized maximum
likelihood estimation for $\mathbf{\Psi}$ under the normality assumption.

Decomposition method | Definition | Graph Product |  Properties
---|---|---|---
Kronecker Product | $\bigotimes_{\ell}\mathbf{\Psi}_{\ell}$ | Tensor Product |  Non-convex
Kronecker Sum | $\bigoplus_{\ell}\mathbf{\Psi}_{\ell}$ | Cartesian Product |  Convex MLE, maximum entropy distribution
Squared Kronecker Sum | $\left(\bigoplus_{\ell}\mathbf{\Psi}_{\ell}\right)^{2}$ | Unnamed111The motivation for this decomposition was based on generative models rather than graph relations. When considered as a product of graphs, it uses second-order edge relations (that is, whether or not two vertices in the output are connected in the input depends on their second-order neighbors in the factor graphs). |  Convex pseudo-likelihood, generative model interpretation
Table 1: A comparison of the three graph decomposition methods used in
independence-free graphical models.

In order to remove the independence assumption in a computationally tractable
way, we need to pick a graph factorization method $\zeta$. Methods that do
this are also called ‘multi-axis’ methods, since they estimate the graphs for
all axes simultaneously. To our knowledge, only three choices have appeared in
the literature; the Kronecker product, the Kronecker sum (Kalaitzis et al.,
2013), and the Squared Kronecker sum (Wang et al., 2020). We summarize these
in Table 1. Among the options, the Kronecker sum stands out as it is
corresponds the maximum entropy distribution (Kalaitzis et al., 2013), and
thus in a sense represents the least strict assumption one can make about how
axes interact. It is also the only one that abides by the example given in
Figure 1, making it particularly attractive in the context of conditional
dependence.

The original method to use the Kronecker sum, BiGLasso, had a non-optimal
space complexity of $O(d_{1}^{2}d_{2}^{2})$, leading it to be unusable on
anything but the smallest problems. Later work, such as scBiGLasso (Li et al.,
2022) and EiGLasso (Yoon and Kim, 2020) focused on improving the runtime and
achieving the optimal $O(d_{1}^{2}+d_{2}^{2})$ space complexity. In
particular, TeraLasso (Greenewald et al., 2019) extended the problem to
tensor-variate datasets (those with arbirtrarily many axes).

Each of the aforementioned methods was iterative, and used an
eigendecomposition every iteration222BiGLasso did not use an
eigendecomposition, but was much slower than the others due to its space
complexity. Eigendecompositions help reduce this.. As an eigendecomposition is
a very slow operation, it is not ideal for scalability. Our previous work,
GmGM (Andrew et al., 2024), focused on removing this requirement, by making
only a single eigendecomposition per axis overall. GmGM also generalized the
problem to work with multi-modality datasets: those with multiple tensors
sharing some axes. However, we were still unable to consider datasets larger
than around ten thousand samples, due to memory constraints.

Other methods that work with multi-modal datasets, but make the independence
assumption, are often based on a regularization penalty, such as the group
graphical lasso and fused graphical lasso (Danaher et al., 2014; Cai et al.,
2016). Another method is that of using Gaussian chain graphical models
(McCarter and Kim, 2014).

We are primarily interested in the application of these methods to omics data,
which is increasingly frequently becoming multi-modal. For an overview of
Gaussian graphical models in this context, see the review by Altenbuchinger et
al. (2020). As this type of data is non-Gaussian, the nonparanormal skeptic
(Liu et al., 2012) has been used to circumvent the normality assumption,
replacing it with the weaker ‘Gaussian copula’ assumption (Li et al., 2022).

###  1.2 Our Contributions

We improve the GmGM method to have quadratic runtime and linear memory usage,
subject to additional (modest) assumptions (Section 2.1). This is an
improvement over prior work’s cubic runtime and quadratic memory usage. In
particular, our work enables the use of these types of methods on datasets
with millions of cells and/or features (Section 3.2.3): any dataset that can
fit in our personal computer’s 8GB RAM (even sparse datasets) can run in less
than an hour. This was previously impossible.

We make these improvements without sacrificing convenient properties of GmGM,
namely its compatibility with the nonparanormal skeptic, its extensibility to
multi-modality and tensor-variate datasets, and its ability to handle sparse
data without densifying (Section 2.6). We validate these claims on both
synthetic and real data (Section 3), and prove existence and uniqueness-up-to-
identifiability of the resultant graphs (Section 2.4). Our hyperparameters are
interpretable, with simple a priori choices for both of them (Section 2.7).
Finally, we give an expression for the Fisher information matrix, allowing
statistical hypothesis tests on the edges of our graph (Section 2.5).

###  1.3 Notation

Our notation primarily follows Kolda and Bader (2009) and that of prior multi-
axis work, with some simplifications made for sanity and consistency. Scalars
will be represented as lower-case letters $a$, vectors as lower-case bold
$\mathbf{v}$, matrices as upper-case bold $\mathbf{M}$, and tensors as upper-
case calligraphic $\mathcal{T}$. Axes will always be denoted with the letter
$\ell$ in subscripts; for example, the precision matrix for axis $\ell$ will
be represented $\mathbf{\Psi}_{\ell}$. (Potentially tensor-variate) datasets
will be represented as $\mathcal{D}$. A dataset may contain multiple tensors,
which we will call modalities. We represent modalities as $\gamma$ in the
superscript, so for example $\mathcal{D}^{\gamma}$ is the tensor corresponding
to modality $\gamma$. The number of axes in a modality is represented as
$L^{\gamma}$, and the number of modalities is $\Gamma$. The overall number of
axes is $L$.

As each modality will have its own set of axes, we will use the terminology
$\ell\in\gamma$ to represent axis $\ell$ being present in
$\mathcal{D}^{\gamma}$. These axes will follow some ordering in $\gamma$,
which may be different in different $\gamma^{\prime}$. Thus it does not make
sense to speak of $\ell<\ell^{\prime}$, but it does relative to a modality;
$\ell<^{\gamma}\ell^{\prime}$.

We will use $d_{\ell}$ to denote the size of the $\ell$th axis of a tensor;
typically, it will be used to represent the size of the axes of our input data
$\mathcal{D}^{\gamma}$. In this paper, we will also use partial
eigendecompositions; $k_{\ell}$ will denote the number of eigenvalues kept for
axis $\ell$. It will be useful to define the following related quantities:

$\displaystyle d^{\gamma}_{<\ell}$
$\displaystyle=\prod_{\ell^{\prime}\in\gamma|\ell^{\prime}<^{\gamma}\ell}d_{\ell}$
$\displaystyle d^{\gamma}_{>\ell}$
$\displaystyle=\prod_{\ell^{\prime}\in\gamma|\ell^{\prime}>^{\gamma}\ell}d_{\ell}$
$\displaystyle d^{\gamma}_{\backslash\ell}$
$\displaystyle=d^{\gamma}_{<\ell}d^{\gamma}_{>\ell}$ $\displaystyle
d^{\gamma}_{\forall}$ $\displaystyle=d^{\gamma}_{\backslash\ell}d_{\ell}$
$\displaystyle d_{\forall}$ $\displaystyle=\sum_{\gamma}d^{\gamma}_{\forall}$
$\displaystyle d_{\forall|\ell}$
$\displaystyle=\sum_{\gamma|\ell\in\gamma}d^{\gamma}_{\forall}$ $\displaystyle
d_{\forall\backslash\ell}$
$\displaystyle=\sum_{\gamma|\ell\in\gamma}d^{\gamma}_{\backslash\ell}$

Analogous quantities can be defined for $k_{\ell}$. Note that
$d_{\forall}^{\gamma}$ represents the total number of entries in tensor
$\gamma$, and $d_{\forall}$ represents the total number of entries in the
entire dataset. $d^{\gamma}_{\forall\backslash\ell}$ represents the effective
number of samples we have when we consider the axis $\ell$. We will use
$n_{\ell}$ to represent the number of edges kept per axis.

The $\ell$-matricization of a tensor is denoted as
$\mathrm{mat}_{\ell}\left[\mathcal{D}^{\gamma}\right]$. It is a
$d_{\ell}\times d_{\backslash\ell}^{\gamma}$ matrix whose main purpose in our
paper is as a definition for the empirical Gram matrices;
$\mathbf{S}_{\ell}^{\gamma}=\mathrm{mat}_{\ell}\left[\mathcal{D}^{\gamma}\right]\mathrm{mat}_{\ell}\left[\mathcal{D}^{\gamma}\right]^{T}$.
We then define the ‘effective Gram matrices’ to be
$\mathbf{S}_{\ell}=\sum_{\gamma|\ell\in\gamma}\mathbf{S}_{\ell}^{\gamma}$;
these are the sufficient statistics of the multi-modality Kronecker-sum-
structured normal distribution; their rank will be at most
$d_{\forall|\backslash\ell}$, with that bound being reached unless there are
linear dependencies in the samples. We will let its top $k_{\ell}$
eigenvectors be denoted $\mathbf{V}^{(k_{\ell})}_{\ell}$ and corresponding
eigenvalues be denoted $\mathbf{E}^{(k_{\ell})}_{\ell}$ (a diagonal matrix).
$\mathbf{V}^{(k_{\ell})}_{\ell}$ will turn out to be eigenvectors of
$\mathbf{\Psi}_{\ell}$ (Theorem 4), so we will not introduce separate
terminology for them, but the corresponding eigenvalues of
$\mathbf{\Psi}_{\ell}$ will be denoted $\mathbf{\Lambda}_{\ell}^{(k_{\ell})}$.

We will denote $\mathbf{I}_{a}$ to be an identity matrix of size $a\times a$.
The Kronecker sum $\bigoplus$ of a set of matrices $\mathbf{\Psi}_{\ell}$ is
defined as
$\bigoplus_{\ell}\mathbf{\Psi}_{\ell}=\sum_{\ell}\mathbf{I}_{d_{<\ell}}\otimes\mathbf{\Psi}_{\ell}\otimes\mathbf{I}_{d_{>\ell}}$.
A matrix of only ones will be denoted $\mathbf{J}$. A matrix of zeros except
for a 1 at position $(i,j)$ will be denoted $\mathbf{J}^{ij}$, and analogously
a vector of zeros except at $i$ will be $\mathbf{j}^{i}$. This is used to
define the stridewise-blockwise trace, which shows up often in derivatives
involving the inverse of a Kronecker sum:

$\displaystyle\mathrm{tr}\left[\mathbf{M}\left(\mathbf{I}_{a}\otimes\mathbf{J}^{ij}\otimes\mathbf{I}_{b}\right)\right]$
$\displaystyle=\mathrm{tr}^{a}_{b}\left[\mathbf{M}\right]_{ij}$

$\mathrm{tr}^{a}_{b}\left[\mathbf{M}\right]$ is then the matrix whose
$(i,j)$th element is $\mathrm{tr}^{a}_{b}\left[\mathbf{M}\right]_{ij}$.

##  2 Methods

In this section, we will demonstrate the need for additional assumptions to
derive the MLE efficiently (Section 2.1), and demonstrate how to use these
assumptions to get a performant algorithm (Section 2.2). We will then prove
additional properties of the algorithm, such as the existence and uniqueness
of the estimated graph (Section 2.3 \- 2.5). We next describe how to implement
the algorithm practically, emphasizing how to avoid creating dense
intermediate products that worsen the memory complexity (Section 2.6).
Finally, we give recommendations on how to set the hyperparameters (Section
2.7).

###  2.1 Assumptions and Their Justification

Our method inherits the nonparanormality assumption from prior work; i.e., the
data has a Gaussian copula. Intuitively, this corresponds to a multivariate
distribution with arbitrary marginals but whose variables interact with each
other ‘like Gaussian variables’. See Liu et al. (2012) for more details. This
assumption is made primarily for computational convenience. To our knowledge,
there has been no multi-axis work that does not make this assumption. If there
were, it would likely lack crucial properties of the Gaussian copula that make
it efficient to scale (namely the relationship between $\mathbf{\Psi}$ and the
graph of conditional dependencies) - although it is certainly a worthwhile
avenue for future work to explore.

The main bottleneck to scalability in the prior version of GmGM was that its
memory complexity was $O(\sum_{\ell}d_{\ell}^{2})$, due to its use of a full
eigendecomposition during its calculation. This is in fact the optimal memory
complexity, without formally assuming sparsity, so any method that improves on
this must both make a formal sparsity assumption and avoid doing a full
eigendecomposition. The ideal method would be eigendecomposition-free, but an
efficient eigendecomposition-free method has proved elusive. Of all Kronecker-
sum-based methods, only the original (Kalaitzis et al., 2013) is
eigendecomposition-free - and it is not scalable due to other factors.

Thus, we propose the following two assumptions:

  1. 1.

The datasets $\mathcal{D}^{\gamma}$ and true graphs $\mathbf{\Psi}_{\ell}$ are
well-approximated by low-rank matrices.

  2. 2.

The graphs are ‘linearly sparse’; graphs with $d_{\ell}$ vertices have
$O(d_{\ell})$ edges rather than $O(d_{\ell}^{2})$. This corresponds to
assumption A4 of Greenewald et al. (2019).

The first gives us a theoretical basis to avoid full eigendecompositions, and
the second formally allows us to avoid an $O(\sum_{\ell}d_{\ell}^{2})$ memory
requirement. It is worth noting that all previous algorithms make an informal
sparsity assumption - this is the rational for their use of the L1 penalty -
we are just formalizing a precise level of sparsity that we expect. Likewise,
low-rank approximations are common; every time one picks the top $K$ principal
components for downstream analysis, they are making a low-rank approximation.

The choice of at-most-linear sparsity follows from the assumption that “there
is a constant $A$ such that no vertex can have more than $A$ edges, regardless
of the size of the graph”. This is often a reasonable assumption; for example,
in a network of friendships, there is a finite amount of time a given person
has to maintain their friendships, putting a bound on their vertex’s degree
regardless of how many people live on the planet. Similar arguments can be
made for most networks of interest. From a computational complexity
perspective, there is no point in assuming sub-linear sparsity, as partial
eigendecompositions will become the memory bottleneck. This justifies our
choice of linear sparsity.

The assumption that has the most effect on our methodology is that of the low-
rank approximation. By assuming that the underlying graphs are low-rank, we
must work with a ‘singular’ probability distribution. To exemplify the issue,
let us look at the probability density function of the normal distribution:

$\displaystyle\mathrm{pdf}_{\text{normal}}(\mathbf{x})$
$\displaystyle=\frac{\sqrt{\mathrm{det}\left(\mathbf{\Psi}\right)}}{\left(2\pi\right)^{\frac{d}{2}}}e^{-\frac{1}{2}\mathbf{x}^{T}\mathbf{\Psi}\mathbf{x}}$

As $\mathbf{\Psi}$ is assumed to be low-rank, its determinant is zero. Thus,
the pdf is zero everywhere, and does not actually integrate to 1 - it is not
technically a pdf. The problem gets worse if one defines the distribution in
terms of its covariance instead. However, this is somewhat of an artificial
problem. $\mathbf{\Psi}$ being singular means that all the probability mass of
the distribution is concentrated in a lower-dimensional subspace. If we
restrict ourselves to this subspace, it is possible to define a density. In
this subspace, it is only the nonzero eigenvalues of $\mathbf{\Psi}$ that
matter, so we replace the determinant with the ‘pseudodeterminant’
$\mathrm{det}^{\dagger}$ (the product of nonzero eigenvalues). The normalizing
constant also changes to $\left(2\pi\right)^{\frac{k}{2}}$, where $k$ is the
rank of the lower-dimensional subspace - although this does not affect the MLE
of the distribution. This modified distribution is known as the ‘singular
normal distribution’.

$\displaystyle\mathrm{pdf}_{\text{normal}}(\mathbf{x})$
$\displaystyle=\frac{\sqrt{\mathrm{det}^{\dagger}\left(\mathbf{\Psi}\right)}}{\left(2\pi\right)^{\frac{k}{2}}}e^{-\frac{1}{2}\mathbf{x}^{T}\mathbf{\Psi}\mathbf{x}}$

It is then straightforward to present a Kronecker-sum-structured version of
this distribution, the ‘singular Kronecker-sum-structured normal
distribution’. We let $k_{\ell}$ be the rank chosen for each axis, and
$k_{\forall}^{\gamma}$ be the product of ranks of all axes occuring in
modality $\gamma$.

$\displaystyle\mathrm{pdf}_{\text{normal}}(\mathbf{\mathcal{D}}^{\gamma})$
$\displaystyle=\frac{\sqrt{\mathrm{det}^{\dagger}\left(\bigoplus_{\ell\in\gamma}\mathbf{\Psi}_{\ell}\right)}}{\left(2\pi\right)^{\frac{k_{\forall}^{\gamma}}{2}}}e^{-\frac{1}{2}\mathrm{vec}\left[\mathcal{D}^{\gamma}\right]^{T}\left(\bigoplus_{\ell\in\gamma}\mathbf{\Psi}_{\ell}\right)\mathrm{vec}\left[\mathcal{D}^{\gamma}\right]}$
$\displaystyle=\frac{\sqrt{\mathrm{det}^{\dagger}\left(\bigoplus_{\ell\in\gamma}\mathbf{\Psi}_{\ell}\right)}}{\left(2\pi\right)^{\frac{k_{\forall}^{\gamma}}{2}}}e^{-\frac{1}{2}\sum_{\ell\in\gamma}\mathrm{tr}\left[\mathbf{S}_{\ell}^{\gamma}\mathbf{\Psi}_{\ell}\right]}$

When given a multi-modal dateset, we have:

$\displaystyle\mathrm{pdf}_{\text{normal}}\left(\left\\{\mathbf{\mathcal{D}}^{\gamma}\right\\}\right)$
$\displaystyle=\prod_{\gamma}\frac{\sqrt{\mathrm{det}^{\dagger}\left(\bigoplus_{\ell\in\gamma}\mathbf{\Psi}_{\ell}\right)}}{\left(2\pi\right)^{\frac{k_{\forall}^{\gamma}}{2}}}e^{-\frac{1}{2}\sum_{\ell\in\gamma}\mathrm{tr}\left[\mathbf{S}_{\ell}^{\gamma}\mathbf{\Psi}_{\ell}\right]}$

The only assumption we will make about the structure of the modalities is that
they are tensors without ‘repeated axes’. An example of matrix with a repeated
axis would be an adjacency matrix; both the rows and the columns represent the
same concept (the vertices of the graph). Shared axes (an axis appearing in
two different modalities) are allowed; taking advantage of shared axes is the
whole point of of the multi-modal setup, after all. The only thing that is
prohibited is two or more instances of the same axis in the same modality.

This is a reasonable assumption, as tensors with repeated axes can always be
interpreted as a graph. This method is meant to find graphs, which is less
important when one already has a graph describing the data. Nevertheless, if
it is absolutely essential to work with repeated axes, there is a way. As our
method takes as input the eigendecompositions of the sufficient statistics
$\left\\{\mathbf{S}_{\ell}\right\\}$, methods such as that of Peshekhonov et
al. (2024) can be used to produce eigendecompositions even in the case of
repeated axes.

###  2.2 Derivation of the Method

It will first be useful to establish a few helpful lemmas about the
stridewise-blockwise trace. The first two are known lemmas, Lemma 1 of Dahl et
al. (2013) and the Cyclic Property of the stridewise-blockwise trace from
Andrew et al. (2024).

######  Lemma (Dahl et al. (2013) Lemma 1).

$\displaystyle\mathrm{tr}\left[\left(\mathbf{I}_{a}\otimes\mathbf{X}\otimes\mathbf{I}_{b}\right)\mathbf{M}\right]$
$\displaystyle=\mathrm{tr}\left[\mathbf{X}\mathrm{tr}^{a}_{b}\left[\mathbf{M}\right]\right]$

######  Lemma (Cyclic Property (Andrew et al., 2024)).

$\displaystyle\mathrm{tr}^{a}_{b}\left[\left(\mathbf{A}_{a}\otimes
I\otimes\mathbf{B}_{b}\right)\mathbf{M}\right]$
$\displaystyle=\mathrm{tr}^{a}_{b}\left[\mathbf{M}\left(\mathbf{A}_{a}\otimes
I\otimes\mathbf{B}_{b}\right)\right]$

The next lemma is a generalized version of lemmas from prior work; an
‘extraction’ property (a generalized Lemma 2 from Dahl et al. (2013)).

######  Lemma 1 (Extraction Property)

$\displaystyle\mathrm{tr}^{a}_{b}\left[\left(\mathbf{I}_{a}\otimes\mathbf{X}\otimes\mathbf{I}_{b}\right)\mathbf{M}\left(\mathbf{I}_{a}\otimes\mathbf{Y}^{T}\otimes\mathbf{I}_{b}\right)\right]$
$\displaystyle=\mathbf{X}\mathrm{tr}^{a}_{b}\left[\mathbf{M}\right]\mathbf{Y}^{T}$

###### Proof.

This proof follows the original by Dahl et al. (2013) closely; the only
differences are that it is proven in the tensor-variate case and that
$\mathbf{X}$ and $\mathbf{Y}$ are not constrained to be the same. The
differences are not significant, and the re-proof here is just for formality.

$\displaystyle\mathrm{tr}^{a}_{b}\left[\left(\mathbf{I}_{a}\otimes\mathbf{X}\otimes\mathbf{I}_{b}\right)\mathbf{M}\left(\mathbf{I}_{a}\otimes\mathbf{Y}^{T}\otimes\mathbf{I}_{b}\right)\right]$
$\displaystyle=\mathrm{tr}\left[\left(\mathbf{I}_{a}\otimes\mathbf{X}\otimes\mathbf{I}_{b}\right)\mathbf{M}\left(\mathbf{I}_{a}\otimes\mathbf{Y}^{T}\otimes\mathbf{I}_{b}\right)\left(\mathbf{I}_{a}\otimes\mathbf{J}^{ij}\otimes\mathbf{I}_{b}\right)\right]_{ij}$
$\displaystyle=\mathrm{tr}\left[\left(\mathbf{I}_{a}\otimes\mathbf{j}^{j}\mathbf{X}\otimes\mathbf{I}_{b}\right)\mathbf{M}\left(\mathbf{I}_{a}\otimes\mathbf{Y}^{T}\mathbf{j}^{i,T}\otimes\mathbf{I}_{b}\right)\right]_{ij}$
$\displaystyle=\mathrm{tr}\left[\left(\mathbf{I}_{a}\otimes\mathbf{X}_{j}\otimes\mathbf{I}_{b}\right)\mathbf{M}\left(\mathbf{I}_{a}\otimes\mathbf{Y}_{i}^{T}\otimes\mathbf{I}_{b}\right)\right]_{ij}$
$\displaystyle=\mathrm{tr}\left[\left(\mathbf{I}_{a}\otimes\otimes\mathbf{Y}_{i}^{T}\mathbf{X}_{j}\otimes\mathbf{I}_{b}\right)\mathbf{M}\right]_{ij}$
$\displaystyle=\mathrm{tr}\left[\mathbf{Y}_{i}^{T}\mathbf{X}_{j}\mathrm{tr}^{a}_{b}\left[\mathbf{M}\right]\right]_{ij}$
$\displaystyle=\mathrm{tr}\left[\mathbf{X}_{j}\mathrm{tr}^{a}_{b}\left[\mathbf{M}\right]\mathbf{Y}_{i}^{T}\right]_{ij}$

Note now that the value inside the trace is a scalar value, and hence we can
drop the outer trace.

$\displaystyle\mathrm{tr}^{a}_{b}\left[\left(\mathbf{I}_{a}\otimes\mathbf{X}\otimes\mathbf{I}_{b}\right)\mathbf{M}\left(\mathbf{I}_{a}\otimes\mathbf{Y}^{T}\otimes\mathbf{I}_{b}\right)\right]$
$\displaystyle=\left[\mathbf{X}_{j}\mathrm{tr}^{a}_{b}\left[\mathbf{M}\right]\mathbf{Y}_{i}^{T}\right]_{ij}$

Which leads us to the final result.

$\displaystyle\mathrm{tr}^{a}_{b}\left[\left(\mathbf{I}_{a}\otimes\mathbf{X}\otimes\mathbf{I}_{b}\right)\mathbf{M}\left(\mathbf{I}_{a}\otimes\mathbf{Y}^{T}\otimes\mathbf{I}_{b}\right)\right]$
$\displaystyle=\mathbf{X}\mathrm{tr}^{a}_{b}\left[\mathbf{M}\right]\mathbf{Y}^{T}$

∎

Finally, we need a new property of the stridewise-blockwise trace, the
‘Downsampling Property’. This is essential to be able to work with rectangular
matrices in the stridewise-blockwise trace, such as those arising from partial
eigendecomposition.

######  Lemma 2 (Downsampling Property)

Suppose $\mathbf{U}_{a\times x}$ and $\mathbf{V}_{a\times x}$ are matrices
such that $\mathbf{V}^{T}\mathbf{U}=\mathbf{I}_{x\times x}$, and
$\mathbf{W}_{b\times y}$ and $\mathbf{X}_{b\times y}$ are matrices such that
$\mathbf{X}^{T}\mathbf{W}=\mathbf{I}_{y\times y}$. Then we have the following.

$\displaystyle\mathrm{tr}^{a}_{b}\left[\left(\mathbf{U}\otimes\mathbf{I}\otimes\mathbf{W}\right)\mathbf{M}\left(\mathbf{V}\otimes\mathbf{I}\otimes\mathbf{X}\right)^{T}\right]$
$\displaystyle=\mathrm{tr}^{x}_{y}\left[\mathbf{M}\right]$

###### Proof.

$\displaystyle\mathrm{tr}^{a}_{b}\left[\left(\mathbf{U}\otimes\mathbf{I}\otimes\mathbf{W}\right)\mathbf{M}\left(\mathbf{V}\otimes\mathbf{I}\otimes\mathbf{X}\right)^{T}\right]$
$\displaystyle=\mathrm{tr}\left[\left(\mathbf{U}\otimes\mathbf{I}\otimes\mathbf{W}\right)\mathbf{M}\left(\mathbf{V}\otimes\mathbf{I}\otimes\mathbf{X}\right)^{T}\left(\mathbf{I}\otimes\mathbf{J}^{ij}\otimes\mathbf{I}\right)\right]_{ij}$
$\displaystyle=\mathrm{tr}\left[\mathbf{M}\left(\mathbf{V}^{T}\mathbf{U}\otimes\mathbf{J}^{ij}\otimes\mathbf{X}^{T}\mathbf{W}\right)\right]_{ij}$
$\displaystyle=\mathrm{tr}\left[\mathbf{M}\left(\mathbf{I}_{x\times
x}\otimes\mathbf{J}^{ij}\otimes\mathbf{I}_{y\times y}\right)\right]_{ij}$
$\displaystyle=\mathrm{tr}^{x}_{y}\left[\mathbf{M}\right]$

∎

With these lemmas stated, we can now derive our method. It is easy to see that
the negative log-likelihood of the singular Kronecker-sum-structured normal
distribution is:

$\displaystyle\mathrm{NLL}$
$\displaystyle=\sum_{\gamma}\left(\frac{d_{\forall}}{2}\log\left(2\pi\right)+\frac{1}{2}\sum_{\ell^{\prime}}\mathrm{tr}\left[\mathbf{S}_{\ell^{\prime}}^{\gamma}\mathbf{\Psi}_{\ell^{\prime}}\right]-\frac{1}{2}\mathrm{logdet}^{\dagger}\bigoplus_{\ell^{\prime}\in\gamma}\mathbf{\Psi}_{\ell^{\prime}}\right)$
$\displaystyle=\frac{1}{2}\sum_{\ell^{\prime}}\mathrm{tr}\left[\mathbf{S}_{\ell^{\prime}}^{\gamma}\mathbf{\Psi}_{\ell^{\prime}}\right]+\sum_{\gamma}\left(\frac{d_{\forall}}{2}\log\left(2\pi\right)-\frac{1}{2}\mathrm{logdet}^{\dagger}\bigoplus_{\ell^{\prime}\in\gamma}\mathbf{\Psi}_{\ell^{\prime}}\right)$

To derive the MLE eigenvectors, we proceed in a manner similar to Andrew et
al. (2024). A key difference is that we must now differentiate a log-
pseudodeterminant, although it behaves analagously to the derivative of a full
determinant. This is given by Theorem 2.15 of Holbrook (2018):

$\displaystyle\frac{\partial}{\partial\mathbf{A}}\mathrm{logdet}^{\dagger}\left(\mathbf{A}\right)$
$\displaystyle=\mathrm{tr}\left[\mathbf{A}^{\dagger}\right]$

###### Lemma 3

At the maximum likelihood of the singular Kronecker-sum-structured normal
distribution, we have that
$\mathbf{S}_{\ell}=\sum_{\gamma|\ell\in\gamma}\mathrm{tr}^{d^{\gamma}_{<\ell}}_{d^{\gamma}_{>\ell}}\left[\left(\bigoplus_{\ell^{\prime}\in\gamma}\mathbf{\Psi}_{\ell^{\prime}}\right)^{\dagger}\right]$.

###### Proof.

Note that $\mathbf{\Psi}_{\ell}$ is symmetric, and hence $\frac{\partial
f\left(\mathbf{\Psi}\right)}{\partial\mathbf{\Psi}}=\mathrm{sym}\left[\mathbf{X}\right]=\frac{\mathbf{X}+\mathbf{X}^{T}}{2}$,
where $\mathbf{X}$ is the derivative found without taking into account
symmetry (Srinivasan and Panda, 2023). It does not affect the location of the
MLE (although we include it in the derivations for completeness), but it does
affect the value of the gradient, and thus will be important to account for
when the gradient itself is of interest.

$\displaystyle\frac{\partial}{\partial\mathbf{\Psi}_{\ell}}\mathrm{NLL}$
$\displaystyle=\mathrm{sym}\left[\frac{\partial}{\partial\mathbf{\Psi}_{\ell}}\sum_{\ell^{\prime}}\frac{1}{2}\mathrm{tr}\left[\mathbf{S}_{\ell^{\prime}}\mathbf{\Psi}_{\ell^{\prime}}\right]-\sum_{\gamma}\left(\frac{\partial}{\partial\mathbf{\Psi}_{\ell_{ij}}}\frac{1}{2}\mathrm{logdet}^{\dagger}\left[\bigoplus_{\ell^{\prime}\in\gamma}\mathbf{\Psi}_{\ell^{\prime}}\right]_{ij}\right)\right]$
$\displaystyle=\mathrm{sym}\left[\sum_{\gamma}\frac{1}{2}\mathbf{S}_{\ell}^{\gamma}-\frac{1}{2}\sum_{\gamma|\ell\in\gamma}\mathrm{tr}\left[\left(\bigoplus_{\ell^{\prime}\in\gamma}\mathbf{\Psi}_{\ell^{\prime}}\right)^{\dagger}\frac{\partial}{\partial\mathbf{\Psi}_{\ell_{ij}}}\bigoplus_{\ell^{\prime}\in\gamma}\mathbf{\Psi}_{\ell^{\prime}}\right]_{ij}\right]$

It is also not too hard to show that
$\frac{\partial}{\partial\mathbf{\Psi}_{\ell_{ij}}}\bigoplus_{\ell^{\prime}\in\gamma}\mathbf{\Psi}_{\ell^{\prime}}=\mathbf{I}_{d^{\gamma}_{<\ell}}\otimes\mathbf{J}^{ij}_{d_{\ell}\times
d_{\ell}}\otimes\mathbf{I}_{d^{\gamma}_{>\ell}}$ (not accounting for
symmetry).

$\displaystyle\frac{\partial}{\partial\mathbf{\Psi}_{\ell_{ij}}}\bigoplus_{\ell^{\prime}\in\gamma}\mathbf{\Psi}_{\ell^{\prime}}$
$\displaystyle=\frac{\partial}{\partial\mathbf{\Psi}_{\ell_{ij}}}\sum_{\ell^{\prime}\in\gamma}\mathbf{I}_{d^{\gamma}_{<\ell^{\prime}}}\otimes\mathbf{\Psi}_{\ell^{\prime}}\otimes\mathbf{I}_{d^{\gamma}_{>\ell^{\prime}}}$
$\displaystyle=\mathbf{I}_{d^{\gamma}_{<\ell}}\otimes\frac{\partial}{\partial\mathbf{\Psi}_{\ell_{ij}}}\mathbf{\Psi}_{\ell}\otimes\mathbf{I}_{d^{\gamma}_{>\ell}}$
$\displaystyle=\mathbf{I}_{d^{\gamma}_{<\ell}}\otimes\mathbf{J}^{ij}_{d_{\ell}\times
d_{\ell}}\otimes\mathbf{I}_{d^{\gamma}_{>\ell}}$

By the definition of the stridewise-blockwise trace, we have that:

$\displaystyle\frac{\partial}{\partial\mathbf{\Psi}_{\ell}}\mathrm{NLL}$
$\displaystyle=\mathrm{sym}\left[\sum_{\gamma}\frac{1}{2}\mathbf{S}_{\ell}^{\gamma}-\frac{1}{2}\sum_{\gamma|\ell\in\gamma}\mathrm{tr}^{d^{\gamma}_{<\ell}}_{d^{\gamma}_{>\ell}}\left[\left(\bigoplus_{\ell^{\prime}\in\gamma}\mathbf{\Psi}_{\ell^{\prime}}\right)^{\dagger}\right]\right]$

Now, note that the maximum likelihood corresponds to the point where
$\frac{\partial}{\partial\mathbf{\Psi}_{\ell}}\mathrm{NLL}$ is $0$.

$\displaystyle\mathbf{S}_{\ell}=\sum_{\gamma|\ell\in\gamma}\mathrm{tr}^{d^{\gamma}_{<\ell}}_{d^{\gamma}_{>\ell}}\left[\left(\bigoplus_{\ell^{\prime}\in\gamma}\mathbf{\Psi}_{\ell^{\prime}}\right)^{\dagger}\right]$

∎

###### Theorem 4

Let the rank-$k_{\ell}$ partial eigendecomposition of $\mathbf{S}_{\ell}$ be
$\mathbf{V}_{\ell}^{(k)}\mathbf{E}_{\ell}^{(k_{\ell})}\mathbf{V}_{\ell}^{(k_{\ell}),T}$,
with
$k_{\ell}\geq\mathrm{max}\left(\mathrm{rank}\left[\mathbf{S}_{\ell}\right],\mathrm{rank}\left[\mathbf{\Psi}_{\ell}\right]\right)$.
Then $\mathbf{V}_{\ell}^{(k_{\ell})}$ are also eigenvectors of
$\mathbf{\Psi}_{\ell}$.

###### Proof.

In this proof, we will want to leave open the possibility of using partial
eigendecompositions. Suppose we had a rank-$k_{\ell}$ partial
eigendecomposition yielding a $d_{\ell}\times k_{\ell}$ matrix of eigenvectors
$\mathbf{V}^{(k_{\ell})}_{\ell}$; in this case,
$\mathbf{V}^{(k_{\ell}),T}_{\ell}\mathbf{V}^{(k_{\ell})}_{\ell}$ is still the
identity matrix (albeit of a smaller dimension); this will allow us to use
Lemma 2 to ‘downsample’ the stridewise-blockwise trace.

We will use $\mathbf{V}^{\gamma}_{<\ell}$ as a shorthand for
$\bigotimes_{\ell^{\prime}\in\gamma|\ell^{\prime}<^{\gamma}\ell}\mathbf{V}^{(k_{\ell}^{\prime})}_{\ell^{\prime}}$
(with an analogous definition for $\mathbf{V}^{\gamma}_{>\ell}$).

With this in mind, note that by Lemma 3 we have that:

$\displaystyle\mathbf{S}_{\ell}$
$\displaystyle=\sum_{\gamma|\ell\in\gamma}\mathrm{tr}^{d^{\gamma}_{<\ell}}_{d^{\gamma}_{>\ell}}\left[\left(\bigoplus_{\ell^{\prime}\in\gamma}\mathbf{\Psi}_{\ell^{\prime}}\right)^{\dagger}\right]$
$\displaystyle=\sum_{\gamma|\ell\in\gamma}\mathrm{tr}^{d_{<\ell}}_{d^{\gamma}_{>\ell}}\left[\left(\mathbf{V}^{\gamma}_{<\ell}\otimes\mathbf{V}_{\ell}^{(k_{\ell})}\otimes\mathbf{V}^{\gamma}_{>\ell}\right)\left(\bigoplus_{\ell^{\prime}\in\gamma}\mathbf{\Lambda}_{\ell^{\prime}}^{(k_{\ell^{\prime}})}\right)^{\dagger}\left(\mathbf{V}^{\gamma}_{<\ell}\otimes\mathbf{V}_{\ell}^{(k_{\ell})}\otimes\mathbf{V}^{\gamma}_{>\ell}\right)^{T}\right]$
$\displaystyle=\mathbf{V}_{\ell}^{(k_{\ell})}\sum_{\gamma|\ell\in\gamma}\mathrm{tr}^{k^{\gamma}_{<\ell}}_{k^{\gamma}_{>\ell}}\left[\left(\bigoplus_{\ell^{\prime}\in\gamma}\mathbf{\Lambda}_{\ell^{\prime}}^{(k_{\ell^{\prime}})}\right)^{\dagger}\right]\mathbf{V}_{\ell}^{(k_{\ell}),T}$
(Extraction and Downsampling Properties.)

The $\mathbf{V}_{<\ell}^{\gamma}$ and $\mathbf{V}^{\gamma}_{>\ell}$ matrices
depend on $\gamma$ as the precise tensor they are in affects which $\ell$s are
present and in what order. Thankfully, they disappear due to the downsampling
property of sb-trace, allowing us to express our equation as the product of
orthogonal matrices
$\mathbf{V}_{\ell}^{(k_{\ell})},\mathbf{V}^{(k_{\ell}),T}_{\ell}$ and a
diagonal matrix
$\sum_{\gamma|\ell\in\gamma}\mathrm{tr}^{k^{\gamma}_{<\ell}}_{k^{\gamma}_{>\ell}}\left[\left(\bigoplus_{\ell^{\prime}\in\gamma}\mathbf{\Lambda}_{\ell^{\prime}}^{(k_{\ell^{\prime}})}\right)^{\dagger}\right]$.
Thus, this is clearly a partial eigendecomposition of $\mathbf{S}_{\ell}$;
this completes the proof.

∎

###### Corollary 5

Observe the correspondance between how we calculate the eigenvectors and how
PCA establishes the principal components. When trying to establish the value
of $k_{\ell}$, one can use the same techniques as PCA (such as using explained
variance or looking at scree plots).

###### Corollary 6

The proof only required that
$k_{\ell}\geq\mathrm{rank}\left[\mathbf{\Psi}_{\ell}\right]$. Suppose
$k_{\ell}<\mathrm{rank}\left[\mathbf{S}_{\ell}\right]$. The proof is still
valid, but note that the requirement that

$\displaystyle\mathbf{S}_{\ell}=\mathbf{V}_{\ell}^{(k_{\ell})}\sum_{\gamma|\ell\in\gamma}\mathrm{tr}^{k^{\gamma}_{<\ell}}_{k^{\gamma}_{>\ell}}\left[\left(\bigoplus_{\ell^{\prime}\in\gamma}\mathbf{\Lambda}_{\ell^{\prime}}^{(k_{\ell^{\prime}})}\right)^{\dagger}\right]\mathbf{V}_{\ell}^{(k_{\ell}),T}$

can never be met. The left hand side has rank
$\mathrm{rank}\left[\mathbf{S}_{\ell}\right]>k_{\ell}$ and the right hand side
has rank at most $k_{\ell}$! This means that the MLE does not exist. This is
where the assumption that $\mathbf{S}_{\ell}$ can be well-approximated by a
low rank matrix is useful; we can instead use this low-rank approximation as
our sufficient statistic. We analyze existence and uniqueness of the MLE in
Section 2.4.

Suppose $\mathbf{\Psi}_{\ell}$ is a rank $k_{\ell}<d_{\ell}$ matrix. If we do
a full eigendecomposition, then some eigenvectors will correspond to zero
eigenvalues. These eigenvectors aren’t useful; they play no role in the value
of $\mathbf{\Psi}_{\ell}$. Instead, it would be better to do a partial
eigendecomposition.

All equalities in this proof remain exact if one does a decomposition of rank
$\mathrm{max}\left(\mathrm{rank}\left[\mathbf{S}_{\ell}\right],\mathrm{rank}\left[\mathbf{\Psi}_{\ell}\right]\right)$.
If one does a decomposition of smaller rank, the equalities become
approximate, as we are throwing away some information that contributes to the
value of $\mathbf{\Psi}_{\ell}$ or $\mathbf{S}_{\ell}$. However, in practice a
matrix can often be well-approximated by one of smaller rank, with some of the
nonzero principal components corresponding to noise. As we will see in Section
2.4, it is actually important to use a decomposition of rank at most
$\mathrm{min}\left(\mathrm{rank}\left[\mathbf{S}_{\ell}\right],\mathrm{rank}\left[\mathbf{\Psi}_{\ell}\right]\right)$
for the existence of an MLE. In practice,
$\mathrm{rank}\left[\mathbf{S}_{\ell}\right]=\sum_{\gamma|\ell\in\gamma}d^{\gamma}_{\backslash\ell}$.

We now recast the negative log-likelihood from a function of
$\mathbf{\Psi}_{\ell}$ into a function of its (partial) eigendecomposition, by
the functional invariance of the MLE. Given our knowledge of the eigenvectors,
we only need to find the eigenvalues.

###### Theorem 7

The gradient of the $\mathrm{NLL}$ with respect to the eigenvalues is given by
$\frac{1}{2}\left(\mathbf{E}_{\ell}^{(k_{\ell})}-\sum_{\gamma|\ell\in\gamma}\mathrm{tr}^{k^{\gamma}_{<\ell}}_{k^{\gamma}_{>\ell}}\left[\left(\bigoplus_{\ell^{\prime}\in\gamma}\mathbf{\Lambda}_{\ell^{\prime}}^{(k_{\ell^{\prime}})}\right)^{\dagger}\right]\right)$.

###### Proof.

As mentioned earlier, we will use the functional invariance of the MLE to
change the problem from finding when
$\frac{\partial}{\partial\mathbf{\Psi}_{\ell}}\mathrm{NLL}=0$ to one of
finding when:

$\displaystyle\frac{\partial}{\partial\mathbf{V}_{\ell}^{(k_{\ell})}}\mathrm{NLL}$
$\displaystyle=0$
$\displaystyle\frac{\partial}{\partial\mathbf{\Lambda}_{\ell}^{(k_{\ell})}}\mathrm{NLL}$
$\displaystyle=0$

We already know the conditions on $\mathbf{V}_{\ell}^{(k_{\ell})}$ at the MLE,
thanks to Theorem 4, so we will focus on
$\mathbf{\Lambda}_{\ell}^{(k_{\ell})}$. We follow a similar derivation to
Lemma 3.

$\displaystyle\frac{\partial}{\partial\mathbf{\Lambda}_{\ell}}\mathrm{NLL}$
$\displaystyle=\frac{\partial}{\partial\mathbf{\Lambda}_{\ell}^{(k_{\ell})}}\sum_{\ell^{\prime}}\frac{1}{2}\mathrm{tr}\left[\mathbf{S}_{\ell^{\prime}}\mathbf{\Psi}_{\ell^{\prime}}\right]-\sum_{\gamma}\left(\frac{\partial}{\partial\mathbf{\Lambda}_{\ell_{ii}}^{(k_{\ell})}}\frac{1}{2}\mathrm{logdet}^{\dagger}\left[\bigoplus_{\ell^{\prime}\in\gamma}\mathbf{\Psi}_{\ell^{\prime}}\right]_{ii}\right)$
$\displaystyle=\frac{\partial}{\partial\mathbf{\Lambda}_{\ell}^{(k_{\ell})}}\sum_{\ell^{\prime}}\frac{1}{2}\mathrm{tr}\left[\mathbf{E}_{\ell^{\prime}}^{(k_{\ell^{\prime}})}\mathbf{\Lambda}_{\ell^{\prime}}^{(k_{\ell^{\prime}})}\right]-\sum_{\gamma}\left(\frac{\partial}{\partial\mathbf{\Lambda}_{\ell_{ii}}^{(k_{\ell})}}\frac{1}{2}\mathrm{logdet}^{\dagger}\left[\bigoplus_{\ell^{\prime}\in\gamma}\mathbf{\Lambda}_{\ell^{\prime}}^{(k_{\ell^{\prime}})}\right]_{ii}\right)$
$\displaystyle=\sum_{\gamma}\frac{1}{2}\mathbf{E}_{\ell}^{(k_{\ell}),\gamma}-\frac{1}{2}\sum_{\gamma|\ell\in\gamma}\mathrm{tr}\left[\left(\bigoplus_{\ell^{\prime}\in\gamma}\mathbf{\Lambda}_{\ell^{\prime}}^{(k_{\ell^{\prime}})}\right)^{\dagger}\frac{\partial}{\partial\mathbf{\Lambda}_{\ell_{ii}}^{(k_{\ell})}}\bigoplus_{\ell^{\prime}\in\gamma}\mathbf{\Lambda}_{\ell^{\prime}}^{(k_{\ell^{\prime}})}\right]_{ii}$
$\displaystyle=\frac{1}{2}\left(\mathbf{E}_{\ell}^{(k_{\ell})}-\sum_{\gamma|\ell\in\gamma}\mathrm{tr}^{k^{\gamma}_{<\ell}}_{k^{\gamma}_{>\ell}}\left[\left(\bigoplus_{\ell^{\prime}\in\gamma}\mathbf{\Lambda}_{\ell^{\prime}}^{(k_{\ell^{\prime}})}\right)^{\dagger}\right]\right)$

In the derivation, we made use of the fact that $\mathbf{S_{\ell}}$ and
$\mathbf{\Psi}_{\ell}$ had the same eigenvectors at the MLE, the cyclic
property of the trace, and the fact that pseudodeterminants only depend on the
nonzero eigenvalues. ∎

###### Corollary 8

The maximum likelihood estimate for $\mathbf{\Psi}_{\ell}$ can be found by
first using Theorem 4 to find $\mathbf{V}_{\ell}$ and then restricting our
problem to just finding the eigenvalues, using Theorem 7.
$\mathbf{\Psi}_{\ell}=\mathbf{V}_{\ell}^{(k_{\ell})}\mathbf{\Lambda}_{\ell}^{(k_{\ell})}\mathbf{V}_{\ell}^{(k_{\ell}),T}$.

###  2.3 Identifiability

The Kronecker sum decomposition produces a non-identifiable representation of
$\bigoplus_{\ell}\mathbf{\Psi}_{\ell}$. Let $\\{c_{\ell}\\}_{\ell\in\gamma}$
be a set of constants such that $\sum_{\ell\in\gamma}c_{\ell}=0$.

$\displaystyle\zeta\left(\\{\mathbf{\Psi}_{\ell}+c_{\ell}\mathbf{I}\\}_{\ell\in\gamma}\right)$
$\displaystyle=\bigoplus_{\ell}\mathbf{\Psi}_{\ell}+c_{\ell}\mathbf{I}$
$\displaystyle=\sum_{\ell\in\gamma}\mathbf{I}_{d^{\gamma}_{<\ell}}\otimes\mathbf{\Psi}_{\ell}\otimes\mathbf{I}_{d^{\gamma}_{>\ell}}+\sum_{\ell\in\gamma}\mathbf{I}_{d^{\gamma}_{<\ell}}\otimes
c_{\ell}\mathbf{I}_{d_{\ell}}\otimes\mathbf{I}_{d^{\gamma}_{>\ell}}$
$\displaystyle=\sum_{\ell\in\gamma}\mathbf{I}_{d^{\gamma}_{<\ell}}\otimes\mathbf{\Psi}_{\ell}\otimes\mathbf{I}_{d^{\gamma}_{>\ell}}+\left(\sum_{\ell\in\gamma}c_{\ell}\right)\mathbf{I}_{d^{\gamma}_{\forall}}$
$\displaystyle=\sum_{\ell\in\gamma}\mathbf{I}_{d^{\gamma}_{<\ell}}\otimes\mathbf{\Psi}_{\ell}\otimes\mathbf{I}_{d^{\gamma}_{>\ell}}$
$\displaystyle=\zeta\left(\\{\mathbf{\Psi}_{\ell}\\}_{\ell\in\gamma}\right)$

Thus, our choice of parameterization of the full precision matrix has forced
us to be unable to estimate the diagonals. We might wonder if there is any
other source of non-identifiablity (other than a constant shift applied to
diagonals). Lemma 9 assures us that this is not so.

###### Lemma 9

The only source of non-identifiability in the model is that of the diagonals
of $\mathbf{\Psi_{\ell}}$.

###### Proof.

Let $\mathbf{\Omega}^{\gamma}=\bigoplus_{\ell\in\gamma}\mathbf{\Psi}_{\ell}$.
Without loss of generality, let us consider only the last axis in $\gamma$,
$\mathbf{\Psi}_{\Gamma}$.

$\displaystyle\mathbf{\Omega}^{\gamma}$
$\displaystyle=\bigoplus_{\ell\in\gamma}\mathbf{\Psi}_{\ell}$
$\displaystyle=\mathbf{I}_{d^{\gamma}_{\backslash\Gamma}}\otimes\mathbf{\Psi}_{\Gamma}+\sum_{\ell\in\gamma\neq\Gamma}\mathbf{I}_{d^{\gamma}_{<\ell}}\otimes\mathbf{\Psi}_{\ell}\otimes\mathbf{I}_{d^{\gamma}_{>\ell}}$
$\displaystyle=\mathbf{I}_{d^{\gamma}_{\backslash\Gamma}}\otimes\mathbf{\Psi}_{\Gamma}+\sum_{\ell\in\gamma\neq\Gamma}\mathbf{I}_{d^{\gamma}_{<\ell}}\otimes\mathbf{\Psi}_{\ell}\otimes\mathbf{I}_{d^{\gamma}_{>\ell,<\Gamma}}\otimes\mathbf{I}_{d_{\Gamma}}$

Note that the first summand is a block-diagonal matrix, whose blocks are
either $\mathbf{0}$ or $\mathbf{\Psi}_{\Gamma}$; this follows from the
definition of the Kronecker product. All other summands are block matrices
whose blocks are multiples of $\mathbf{I}_{d^{\gamma}_{\Gamma}}$. Thus, each
block of $\mathbf{\Omega}_{\gamma}$ is some linear combination of the two,
$\mathbf{\Psi}_{\Gamma}+a\mathbf{I}_{d^{\gamma}_{\Gamma}}$. Clearly, the off-
diagonal elements are not affected by the addition of an identity matrix.
Thus, all off-diagonal elements of $\mathbf{\Psi}_{\Gamma}$ are directly
determined by the elements of $\mathbf{\Omega}^{\gamma}$. ∎

Thankfully, non-identifiability of the diagonals is not a severe issue.
Firstly, we are interested in the graph of conditional dependencies, not the
precision matrix itself. This only depends on the support of the precision
matrix. The diagonals do not affect this. Secondly, the gradient only depends
on the identifiable $\bigoplus_{\ell\in\gamma}\mathbf{\Psi}_{\ell}$ (or rather
$\bigoplus_{\ell\in\gamma}\mathbf{\Lambda}_{\ell}$), so non-identifiability
has no effect on the optimization procedure (we will make this statement more
formal in Section 2.4). Thirdly, as Greenewald et al. (2019) point out, if we
are interested in partial correlations then, rather than the typical equation
$\frac{-\psi_{\ell_{ij}}}{\sqrt{\psi_{\ell_{ii}}\psi_{\ell_{jj}}}}$, the
correct equation’s denominators depend on the diagonals of
$\bigoplus_{\ell\in\gamma}\mathbf{\Psi}_{\ell}$.

Nevertheless, it would be beneficial to have an identifiable representation.
In the single-modality case, Greenewald et al. (2019) derived the following
representation:

$\displaystyle\mathbf{\Omega}^{\gamma}$
$\displaystyle=\tau^{\gamma}\mathbf{I}_{d^{\gamma}_{\forall}}+\bigoplus_{\ell\in\gamma}\tilde{\mathbf{\Psi}}_{\ell}$
$\displaystyle\textit{where
}\mathrm{tr}\left[\tilde{\mathbf{\Psi}}_{\ell}\right]$ $\displaystyle=0$

We can map the standard parameterization to this one with ease. As a
shorthand, let $\tau_{\ell}=\mathrm{tr}\left[\mathbf{\Psi}_{\ell}\right]$

$\displaystyle\tau^{\gamma}$
$\displaystyle=\sum_{\ell\in\gamma}d^{\gamma}_{\backslash\ell}\tau_{\ell}$
$\displaystyle\tilde{\mathbf{\Psi}}_{\ell}$
$\displaystyle=\mathbf{\Psi}_{\ell}-\tau_{\ell}\mathbf{I}_{d_{\ell}}$

Unfortunately, the question of identifiability becomes more complex in the
multi-modality case. As an easy example, consider the case where
$\gamma_{1}=(\ell_{1},\ell_{2})$ and
$\gamma_{2}=(\ell_{1},\ell_{2},\ell_{3})$. As $\tau_{\gamma_{1}}$ and
$\tau_{\gamma_{2}}$ are identifiable, we can use the formula
$\tau_{\gamma_{2}}=d_{\ell_{3}}\tau_{\gamma_{1}}+d_{\ell_{1}}d_{\ell_{2}}\tau_{\ell_{3}}$
to work out the trace of the $\ell_{3}$ axis! In this case, such a result is
intuitive; clearly knowledge about $(\ell_{1},\ell_{2})$ and
$(\ell_{1},\ell_{2},\ell_{3})$ would allow us to deduce $(\ell_{3})$. However,
many datasets will have less obvious interactions - what about the dataset
$\gamma_{1}=(\ell_{1},\ell_{2}),\gamma_{2}=(\ell_{2},\ell_{3}),\gamma_{3}=(\ell_{3},\ell_{1})$?

Note that $\\{\tau^{\gamma}\\}$, the traces of
$\\{\mathbf{\Omega}^{\gamma}\\}$, are always a linear combination of
$\\{\tau_{\ell}\\}$. We can express this as:

$\displaystyle\begin{bmatrix}\tau^{\gamma_{1}}\\\ \vdots\\\
\tau^{\gamma_{\Gamma}}\end{bmatrix}$
$\displaystyle=\mathbf{M}_{\boldsymbol{\gamma}}\begin{bmatrix}\tau_{\ell_{1}}\\\
\vdots\\\ \tau_{\ell_{L}}\end{bmatrix}$
$\displaystyle\mathbf{M}_{\boldsymbol{\gamma}}$
$\displaystyle=\begin{bmatrix}d^{\gamma_{1}}_{\backslash\ell_{1}}&...&d^{\gamma_{1}}_{\backslash\ell_{L}}\\\
\vdots&\ddots&\vdots\\\
d^{\gamma_{\Gamma}}_{\backslash\ell_{1}}&...&d^{\gamma_{\Gamma}}_{\backslash\ell_{L}}\end{bmatrix}$

Where for notational convenience we let $d^{\gamma}_{\backslash\ell}=0$ if
$\ell\notin\gamma$. Depending on the modalities and the lengths of each axis,
$\mathbf{M}_{\boldsymbol{\gamma}}$ could be practically any $\Gamma\times L$
matrix with nonnegative integer coefficients, so it is hard to say much about
it in general. However, one thing can be said about
$\mathbf{M}_{\boldsymbol{\gamma}}$ for sure: the output space $\mathfrak{I}$
corresponds to identifiable quantities
($\mathrm{tr}\left[\mathbf{\Omega}^{\gamma}\right]$), and the input space
$\mathfrak{L}$ corresponds to the true latent factors in the model generating
such identifiable quantities.

In particular, $\mathbf{M}_{\boldsymbol{\gamma}}$ tells us how our
identifiable quantities are linked through the latent factors. Vectors with
the property that $\mathbf{M}_{\boldsymbol{\gamma}}\mathbf{x}=\mathbf{0}$
correspond to the non-identifiabilities, as
$\mathbf{M}_{\boldsymbol{\gamma}}\left(\mathbf{x}+\mathbf{y}\right)=\mathbf{M}_{\boldsymbol{\gamma}}\mathbf{y}$.
The nonzero right singular vectors correspond to a basis for the rowspace
$\mathfrak{R}$, and thus a basis for $\mathfrak{L}$ modulo the non-
identifiabilities. We can use these basis vectors of $\mathfrak{R}$ as the
parameters of our model, as we know there is a correspondence between vectors
in $\mathfrak{R}$ and vectors in $\mathfrak{I}$.

Dataset | $\mathbf{M}_{\boldsymbol{\gamma}}$ | Basis Vectors | Parameterization
---|---|---|---
$\begin{matrix}\\\ \gamma_{1}=(\ell_{1},\ell_{2})\end{matrix}$ | $\begin{bmatrix}d_{\ell_{2}}&d_{\ell_{1}}\end{bmatrix}$ | $\begin{bmatrix}\mathbf{e}_{1}\end{bmatrix}=\begin{bmatrix}d_{\ell_{2}}&d_{\ell_{1}}\end{bmatrix}$ | $\mathbf{\Omega}^{\gamma_{1}}=t_{1}\mathbf{I}_{d^{\gamma_{1}}_{\forall}}+\bigotimes_{\ell\in\gamma_{1}}\tilde{\mathbf{\Psi}}_{\ell}$
$\begin{matrix}\\\ \gamma_{1}=(\ell_{1},...,\ell_{L})\end{matrix}$ | $\begin{bmatrix}d_{\backslash\ell_{1}}&...&d_{\backslash\ell_{L}}\end{bmatrix}$ | $\begin{bmatrix}\mathbf{e}_{1}\end{bmatrix}=\begin{bmatrix}d_{\backslash\ell_{1}}&...&d_{\backslash\ell_{L}}\end{bmatrix}$ | $\mathbf{\Omega}^{\gamma_{1}}=t_{1}\mathbf{I}_{d^{\gamma_{1}}_{\forall}}+\bigotimes_{\ell\in\gamma_{1}}\tilde{\mathbf{\Psi}}_{\ell}$
$\left\\{\begin{matrix}\\\ \gamma_{i}=(\ell_{0},\ell_{i})\\\ \hskip 5.0pt\end{matrix}\right\\}_{1\leq i\leq\Gamma}$ | $\begin{bmatrix}d_{\ell_{1}}&d_{\ell_{0}}&0&...\\\ d_{\ell_{2}}&0&d_{\ell_{0}}&...\\\ d_{\ell_{3}}&0&0&...\\\ \vdots&\vdots&\vdots&\ddots\end{bmatrix}$ | $\begin{bmatrix}\mathbf{e}_{1}\\\ \vdots\end{bmatrix}=\begin{bmatrix}d_{\ell_{1}}&d_{\ell_{0}}&0&...\\\ d_{\ell_{2}}&0&d_{\ell_{0}}&...\\\ d_{\ell_{3}}&0&0&...\\\ \vdots&\vdots&\vdots&\ddots\end{bmatrix}$ | $\begin{bmatrix}\mathbf{\Omega}^{\gamma_{1}}\\\ \vdots\end{bmatrix}=\begin{bmatrix}t_{1}\mathbf{I}_{d^{\gamma_{1}}_{\forall}}+\bigotimes_{\ell\in\gamma_{1}}\tilde{\mathbf{\Psi}}_{\ell}\\\ \vdots\end{bmatrix}$
$\begin{matrix}\\\ \gamma_{1}=(\ell_{1},\ell_{2})\\\ \gamma_{2}=(\ell_{2},\ell_{3})\\\ \gamma_{3}=(\ell_{3},\ell_{1})\end{matrix}$ | $\begin{bmatrix}d_{\ell_{2}}&d_{\ell_{1}}&0\\\ 0&d_{\ell_{3}}&d_{\ell_{2}}\\\ d_{\ell_{3}}&0&d_{\ell_{1}}\end{bmatrix}$ | $\begin{bmatrix}\mathbf{e}_{1}\\\ \mathbf{e}_{2}\\\ \mathbf{e}_{3}\end{bmatrix}=\begin{bmatrix}d_{\ell_{2}}&d_{\ell_{1}}&0\\\ 0&d_{\ell_{3}}&d_{\ell_{2}}\\\ d_{\ell_{3}}&0&d_{\ell_{1}}\end{bmatrix}$ | $\begin{bmatrix}\mathbf{\Omega}^{\gamma_{1}}\\\ \mathbf{\Omega}^{\gamma_{2}}\\\ \mathbf{\Omega}^{\gamma_{3}}\end{bmatrix}=\begin{bmatrix}t_{1}\mathbf{I}_{d^{\gamma_{1}}_{\forall}}+\bigotimes_{\ell\in\gamma_{1}}\tilde{\mathbf{\Psi}}_{\ell}\\\ t_{2}\mathbf{I}_{d^{\gamma_{2}}_{\forall}}+\bigotimes_{\ell\in\gamma_{2}}\tilde{\mathbf{\Psi}}_{\ell}\\\ t_{3}\mathbf{I}_{d^{\gamma_{3}}_{\forall}}+\bigotimes_{\ell\in\gamma_{3}}\tilde{\mathbf{\Psi}}_{\ell}\end{bmatrix}$
$\begin{matrix}\\\ \gamma_{1}=(\ell_{1},\ell_{2})\\\ \gamma_{2}=(\ell_{1},\ell_{2},\ell_{3})\\\ \gamma_{3}=(\ell_{2},\ell_{3},\ell_{4})\\\ \hskip 5.0pt\end{matrix}$ | $\begin{bmatrix}d_{\ell_{2}}&d_{\ell_{1}}&0&0\\\ d_{\ell_{2}}d_{\ell_{3}}&d_{\ell_{1}}d_{\ell_{3}}&d_{\ell_{1}}d_{\ell_{2}}&0\\\ 0&d_{\ell_{3}}d_{\ell_{4}}&d_{\ell_{2}}d_{\ell_{4}}&d_{\ell_{2}}d_{\ell_{3}}\end{bmatrix}$ | $\begin{bmatrix}\mathbf{e}_{1}\\\ \mathbf{e}_{2}\\\ \mathbf{e}_{3}\end{bmatrix}=\begin{bmatrix}d_{\ell_{2}}&d_{\ell_{1}}&0&0\\\ d_{\ell_{2}}d_{\ell_{3}}&d_{\ell_{1}}d_{\ell_{3}}&d_{\ell_{1}}d_{\ell_{2}}&0\\\ 0&d_{\ell_{3}}d_{\ell_{4}}&d_{\ell_{2}}d_{\ell_{4}}&d_{\ell_{2}}d_{\ell_{3}}\end{bmatrix}$ | $\begin{bmatrix}\mathbf{\Omega}^{\gamma_{1}}\\\ \mathbf{\Omega}^{\gamma_{2}}\\\ \mathbf{\Omega}^{\gamma_{3}}\end{bmatrix}=\begin{bmatrix}t_{1}\mathbf{I}_{d^{\gamma_{1}}_{\forall}}+\bigotimes_{\ell\in\gamma_{1}}\tilde{\mathbf{\Psi}}_{\ell}\\\ t_{2}\mathbf{I}_{d^{\gamma_{2}}_{\forall}}+\bigotimes_{\ell\in\gamma_{2}}\tilde{\mathbf{\Psi}}_{\ell}\\\ t_{3}\mathbf{I}_{d^{\gamma_{3}}_{\forall}}+\bigotimes_{\ell\in\gamma_{3}}\tilde{\mathbf{\Psi}}_{\ell}\end{bmatrix}$
$\begin{matrix}\\\ \gamma_{1}=(\ell_{1},\ell_{2})\\\ \gamma_{2}=(\ell_{2},\ell_{1})\end{matrix}$ | $\begin{bmatrix}d_{\ell_{2}}&d_{\ell_{1}}\\\ d_{\ell_{2}}&d_{\ell_{1}}\end{bmatrix}$ | $\begin{bmatrix}\mathbf{e}_{1}\end{bmatrix}=\begin{bmatrix}d_{\ell_{2}}&d_{\ell_{1}}\\\ \end{bmatrix}$ | $\begin{bmatrix}\mathbf{\Omega}^{\gamma_{1}}\\\ \mathbf{\Omega}^{\gamma_{2}}\end{bmatrix}=\begin{bmatrix}t_{1}\mathbf{I}_{d^{\gamma_{1}}_{\forall}}+\bigotimes_{\ell\in\gamma_{1}}\tilde{\mathbf{\Psi}}_{\ell}\\\ t_{1}\mathbf{I}_{d^{\gamma_{2}}_{\forall}}+\bigotimes_{\ell\in\gamma_{2}}\tilde{\mathbf{\Psi}}_{\ell}\end{bmatrix}$
$\begin{matrix}\\\ \gamma_{1}=(\ell_{1},\ell_{2})\\\ \gamma_{2}=(\ell_{3},\ell_{4})\\\ \gamma_{3}=(\ell_{1},\ell_{2},\ell_{3},\ell_{4})\\\ \hskip 5.0pt\end{matrix}$ | $\begin{bmatrix}d_{\ell_{2}}&d_{\ell_{1}}&0&0\\\ 0&0&d_{\ell_{4}}&d_{\ell_{3}}\\\ d_{\backslash\ell_{1}}&d_{\backslash\ell_{2}}&d_{\backslash\ell_{3}}&d_{\backslash\ell_{4}}&\end{bmatrix}$ | $\begin{bmatrix}\mathbf{e}_{1}\\\ \mathbf{e}_{2}\\\ \end{bmatrix}=\begin{bmatrix}d_{\ell_{2}}&d_{\ell_{1}}&0&0\\\ 0&0&d_{\ell_{4}}&d_{\ell_{3}}\end{bmatrix}$ | $\begin{bmatrix}t_{1}\mathbf{I}_{d^{\gamma_{1}}_{\forall}}+\bigotimes_{\ell\in\gamma_{1}}\tilde{\mathbf{\Psi}}_{\ell}\\\ t_{2}\mathbf{I}_{d^{\gamma_{2}}_{\forall}}+\bigotimes_{\ell\in\gamma_{2}}\tilde{\mathbf{\Psi}}_{\ell}\\\ (d_{\ell_{3}}d_{\ell_{4}}t_{1}+d_{\ell_{1}}d_{\ell_{2}}t_{2})\mathbf{I}_{d^{\gamma_{2}}_{\forall}}+\bigotimes_{\ell\in\gamma_{2}}\tilde{\mathbf{\Psi}}_{\ell}\end{bmatrix}$

Table 2: Identifiable parameterizations for several datasets. Here, $t_{i}$ is
the coefficient corresponding to basis vector $\mathbf{e}_{i}$, i.e, $t_{i}$
is one of the parameters, along with $\tilde{\mathbf{\Psi}}_{\ell}$. Note that
many different bases are possible; our selection here is somewhat arbitrary.
As can be seen, typically $\mathbf{M}_{\gamma}$ itself can be used for the
basis vectors, but sometimes there can be linear dependencies in the data. The
first three rows in this matrix correspond to the cases most likely to arise
in practice; a matrix-variate dataset, a tensor-variate dataset, and a multi-
modality matrix-variate dataset with a consistently shared axis.

Due to the generality of $\mathbf{M}_{\boldsymbol{\gamma}}$, we are not able
to say more about the basis of $\mathfrak{R}$. For any given
$\mathbf{M}_{\boldsymbol{\gamma}}$, calculating the rowspace is a solved
problem. Example parameterizations for several datasets are given in Table 2.
For most datasets, $\mathbf{M}_{\boldsymbol{\gamma}}$ has maximal row rank, in
which case the approach of Greenewald et al. (2019) is unchanged, and we can
parameterize by $\\{\tau^{\gamma}\\}$. However, there are sometimes linear
dependencies in $\\{\tau^{\gamma}\\}$, causing us to choose a new, smaller
parameterization.

The probability distribution function with respect to the identifiable
representation is given below. $\mathbf{a}^{\gamma}$ corresponds to the
coefficients of the trace parameters $\mathbf{t}$. In the first row of Table
2, $\mathbf{a}^{1}=\left[1\right]$, and in the last row
$\begin{bmatrix}\mathbf{a}^{1}\\\ \mathbf{a}^{2}\\\
\mathbf{a}^{3}\end{bmatrix}=\begin{bmatrix}1&0\\\ 0&1\\\
d_{\ell_{3}}d_{\ell_{4}}&d_{\ell_{1}}d_{\ell_{2}}\end{bmatrix}$.

$\displaystyle\mathrm{pdf}\left(\left\\{\mathcal{T}^{\gamma}\right\\}_{\gamma}\right)$
$\displaystyle=\prod_{\gamma}\frac{\sqrt{\mathrm{det}^{\dagger}\left[\mathbf{t}^{T}\mathbf{a}^{\gamma}\mathbf{I}_{d^{\gamma}_{\forall}}+\bigoplus_{\ell\in\gamma}\mathbf{\Psi_{\ell}}\right]}}{\left(2\pi\right)^{\frac{k^{\gamma}_{\forall}}{2}}}e^{-\frac{1}{2}\mathrm{vec}\left[\mathcal{T}^{\gamma}\right]^{T}\left(\mathbf{t}^{T}\mathbf{a}^{\gamma}\mathbf{I}_{d^{\gamma}_{\forall}}+\bigoplus_{\ell\in\gamma}\mathbf{\Psi}_{\ell}\right)\mathrm{vec}\left[\mathcal{T}^{\gamma}\right]}$
$\displaystyle=\prod_{\gamma}\frac{\sqrt{\mathrm{det}^{\dagger}\left[\mathbf{t}^{T}\mathbf{a}^{\gamma}\mathbf{I}_{d^{\gamma}_{\forall}}+\bigoplus_{\ell\in\gamma}\mathbf{\Psi_{\ell}}\right]}}{\left(2\pi\right)^{\frac{k^{\gamma}_{\forall}}{2}}}e^{-\frac{1}{2}\left(\sum_{\ell\in\gamma}\mathrm{tr}\left[\mathbf{S}^{\gamma}_{\ell}\mathbf{\Psi}_{\ell}\right]+\frac{\mathbf{t}^{T}\mathbf{a}^{\gamma}}{L^{\gamma}}\sum_{\ell\in\gamma}\mathrm{tr}\left[\mathbf{S}^{\gamma}\right]\right)}$
$\displaystyle\mathrm{NLL}$
$\displaystyle=\frac{-1}{2}\sum_{\gamma}\left(\log\mathrm{det}^{\dagger}\left[\mathbf{t}^{T}\mathbf{a}^{\gamma}\mathbf{I}_{d_{\forall}^{\gamma}}+\bigoplus_{\ell\in\gamma}\mathbf{\Psi}_{\ell}\right]-\sum_{\ell\in\gamma}\mathrm{tr}\left[\mathbf{S}_{\ell}^{\gamma}\mathbf{\Psi}_{\ell}\right]-\sum_{\ell\in\gamma}\frac{\mathbf{t}^{T}\mathbf{a}^{\gamma}}{L^{\gamma}}\mathrm{tr}\left[\mathbf{S}^{\gamma}_{\ell}\right]\right)$

We removed the constant $k^{\gamma}_{\forall}\log 2\pi$ term from the negative
log likelihood, as it has no effect on the minimum nor the derivatives.

In practice, this choice of parameterization is not important. As mentioned
earlier, it has no bearing on the graph of conditional dependencies, nor (as
we will prove in Section 2.4) does it affect the algorithm that finds them. As
noted by Yoon and Kim (2020), one does not need to convert to an identifiable
representation until after finding the MLE (w.r.t. the standard, non-
identifiable representation). However, identifiable representations will be
useful in the derivation of a statistical hypothesis test in Section 2.5.

###  2.4 Convexity, Existence, and Uniqueness

Theorem 7 suggests an iterative approach to finding the eigenvalues, via
gradient descent. However, this approach would not work if the problem does
not have a minimum. In this section, we will show that the problem is convex
and does indeed have a unique minimum, up to non-identifiability.

###### Lemma 10

This iterative approach corresponds to a convex optimization problem over the
domain $\boldsymbol{\Lambda}_{\ell}^{(k_{\ell})}\in\mathbb{R}_{k_{\ell}}^{++}$

###### Proof.

Recall that in Theorem 7, our negative log likelihood was found to be

$\displaystyle\mathrm{NLL}$
$\displaystyle=\sum_{\ell^{\prime}}\frac{1}{2}\mathrm{tr}\left[\mathbf{E}_{\ell^{\prime}}^{(k_{\ell^{\prime}})}\mathbf{\Lambda}_{\ell^{\prime}}^{(k_{\ell^{\prime}})}\right]-\sum_{\gamma}\left(\frac{1}{2}\mathrm{logdet}^{\dagger}\left[\bigoplus_{\ell^{\prime}\in\gamma}\mathbf{\Lambda}_{\ell^{\prime}}^{(k_{\ell^{\prime}})}\right]\right)$

The relevant matrices are diagonal matrices, so to make the properties of this
function more obvious we can write it in terms of their diagonals.

$\displaystyle\mathrm{NLL}$
$\displaystyle=\sum_{\ell^{\prime}}\frac{1}{2}\mathbf{e}_{\ell^{\prime}}^{(k_{\ell^{\prime}}),T}\boldsymbol{\lambda}_{\ell^{\prime}}^{(k_{\ell^{\prime}})}-\sum_{\gamma}\left(\frac{1}{2}\mathrm{logdet}^{\dagger}\left[\bigoplus_{\ell^{\prime}\in\gamma}\mathrm{diag}\boldsymbol{\lambda}_{\ell^{\prime}}^{(k_{\ell^{\prime}})}\right]\right)$

Note that
$\mathbf{e}_{\ell}^{(k_{\ell^{\prime}}),T}\boldsymbol{\lambda}_{\ell}^{(k_{\ell^{\prime}})}$
is linear, and thus convex. Kronecker sum is a linear operation, and negative
log-determinants are convex. Since the second summand is a convex function of
a linear transform, it is convex. Thus, the negative log-likelihood is a
convex function over the (convex) domain of
$\boldsymbol{\lambda}_{\ell}^{(k_{\ell^{\prime}})}\in\mathbb{R}_{k_{\ell}}^{++}$.

This completes the proof. Due to the non-identifiabilities, the Kronecker sum
is not a full-rank linear transform, which is what prevents this function from
being strictly convex. ∎

Lemma 10 establishes the convexity of the objective, but it does not establish
the existence of an MLE.

###### Lemma 11

There exists a minimizer of the negative log likelihood in the domain
$\boldsymbol{\Lambda}_{\ell}^{(k_{\ell})}\geq\epsilon_{\ell}$, where
$\epsilon_{\ell}$ is some strictly positive constant (for each $\ell$), as
long as $\mathbf{E}_{\ell}^{(k_{\ell})}\in\mathbb{R}^{++}_{k_{\ell}}$.

###### Proof.

First, note that $\bigtimes_{\ell}\bigtimes_{1\leq i\leq
k_{\ell}}[\epsilon_{\ell},\infty]$ is a compact subset of the extended reals,
and thus a minimum to any continuous function must exist on this set. As it is
convex, it will be a global minimum. However, the minimum could exist at
$\infty$, i.e., only in the limit. To rule out this case, we will establish an
upper bound on the solution.

Recall that any minimum must satisfy
$\mathbf{E}_{\ell}^{(k_{\ell})}=\sum_{\gamma|\ell\in\gamma}\mathrm{tr}^{k^{\gamma}_{<\ell}}_{k^{\gamma}_{>\ell}}\left[\left(\oplus_{\ell^{\prime}\in\gamma}\mathbf{\Lambda}_{\ell^{\prime}}^{(k_{\ell^{\prime}})}\right)^{\dagger}\right]$
(for each $\ell$). It will be helpful to see what this formula equates to
elementwise:

$\displaystyle\mathbf{E}_{\ell_{ii}}^{(k_{\ell})}$
$\displaystyle=\sum_{\gamma|\ell\in\gamma}\mathrm{tr}^{d^{\gamma}_{<\ell}}_{d^{\gamma}_{>\ell}}\left[\left(\oplus_{\ell^{\prime}\in\gamma}\mathbf{\Lambda}_{\ell^{\prime}}^{(k_{\ell^{\prime}})}\right)^{\dagger}\right]_{ii}$
$\displaystyle=\sum_{\gamma|\ell\in\gamma}\mathrm{tr}\left[\left(\oplus_{\ell^{\prime}\in\gamma}\mathbf{\Lambda}_{\ell^{\prime}}^{(k_{\ell^{\prime}})}\right)^{\dagger}\left(\mathbf{I}_{k_{<\ell}^{\gamma}}\otimes\mathbf{J}^{ii}\otimes\mathbf{I}_{k_{>\ell}^{\gamma}}\right)\right]$
$\displaystyle=\sum_{\gamma|\ell\in\gamma}\sum_{\text{ }\mathbf{j}\text{
indexing all non-$\ell$
axes}}\frac{1}{\lambda_{\ell_{i}}+\sum_{\ell^{\prime}}\lambda_{\ell^{\prime}_{j_{\ell^{\prime}}}}}$

Thus, every element of $\mathbf{E}_{\ell}^{(k_{\ell})}$ has a corresponding
element of $\mathbf{\Lambda}_{\ell}^{(k_{\ell})}$ that appears in all summands
composing it, and vice versa. Suppose
$\lambda_{\ell_{i}}=\frac{\sum_{\gamma}d_{\backslash\ell}^{\gamma}}{\mathbf{E}_{\ell_{ii}}}$
(for a given $\ell$, and every $i$). Note that if for all other
$\ell^{\prime}$, $\lambda_{\ell^{\prime}_{j_{\ell^{\prime}}}}$ were 0, then
the sum of these $\lambda_{\ell}$ terms would sum to exactly
$\mathbf{E}_{\ell_{ii}}$. However, as
$\boldsymbol{\Lambda}_{\ell}^{(k_{\ell})}\geq\epsilon_{\ell}$, each reciprocal
is smaller than the case where they are 0; thus, they cannot add up to
$\mathbf{E}_{\ell_{ii}}$.

As $\lambda_{\ell_{i}}$ increases above this bound, the problem only becomes
more severe. This is easy to see by checking the gradient, which indicates
that the negative log-likelihood increases as one increases
$\lambda_{\ell_{i}}$.

$\displaystyle\frac{\partial}{\partial\mathbf{\Lambda}_{\ell_{ii}}^{(k_{\ell})}}\mathrm{NLL}$
$\displaystyle=\frac{1}{2}\left(\mathbf{E}_{\ell_{ii}}-\sum_{\gamma|\ell\in\gamma}\sum_{\text{
}\mathbf{j}\text{ indexing all non-$\ell$
axes}}\frac{1}{\lambda_{\ell_{i}}+\sum_{\ell^{\prime}}\lambda_{\ell^{\prime}_{j_{\ell^{\prime}}}}}\right)$
$\displaystyle>\frac{1}{2}\left(\mathbf{E}_{\ell_{ii}}-\mathbf{E}_{\ell_{ii}}\right)$
$\displaystyle=0$

This implies that the minimum must occur between
$\frac{\sum_{\gamma}d_{\backslash\ell}^{\gamma}}{\mathbf{E}_{\ell_{ii}}}\geq\lambda_{\ell_{ii}}\geq\epsilon_{\ell}$.
Clearly, this rules out the unsavory case of the minimum occurring in the
limit. This completes the proof.

∎

###### Corollary 12

The assumption that $\mathbf{E}_{\ell}^{(k_{\ell})}$ is strictly positive is
necessary; as it is the sum of reciprocals of strictly positive numbers, it
could only be zero in the limit as such numbers go to infinity. This is
reflected in the fact that that the upper bound contains a division by the
elements of $\mathbf{E}_{\ell}^{(k_{\ell})}$.

One might wonder if a different choice of factorization $\zeta$ would avoid
identifiability issues. However, both the Kronecker product and Sylvester sum
methods are non-identifiable at the diagonals as well.

Due to the non-identifiability, the problem has multiple minima. One might
wonder: are all alternative minima due to non-identifiability?

###### Lemma 13

All minima of
$\mathrm{NLL}=-\log\prod_{\gamma}\frac{\sqrt{\mathrm{det}^{\dagger}\left(\bigoplus_{\ell\in\gamma}\mathbf{\Psi}_{\ell}\right)}}{\left(2\pi\right)^{\frac{k_{\forall}^{\gamma}}{2}}}e^{-\frac{1}{2}\mathrm{vec}\left[\mathcal{D}^{\gamma}\right]^{T}\left(\bigoplus_{\ell\in\gamma}\mathbf{\Psi}_{\ell}\right)\mathrm{vec}\left[\mathcal{D}^{\gamma}\right]}$
are related by the non-identifiabilities induced by the choice of $\zeta$
(over the domain of rank-$k_{\ell}$ positive semidefinite matrices).

###### Proof.

Let $\mathbf{\Omega}^{\gamma}=\bigoplus_{\ell\in\gamma}\mathbf{\Psi}_{\ell}$,
and suppose we were to consider the negative log likelihood a function of
$\mathbf{\Omega}^{\gamma}$ instead. This becomes the case of estimating the
precision matrix of a (singular, multi-dataset) normal distribution:

$\displaystyle\mathrm{NLL}$
$\displaystyle=-\log\prod_{\gamma}\frac{\sqrt{\mathrm{det}^{\dagger}\left(\bigoplus_{\ell\in\gamma}\mathbf{\Psi}_{\ell}\right)}}{\left(2\pi\right)^{\frac{k_{\forall}^{\gamma}}{2}}}e^{-\frac{1}{2}\mathrm{vec}\left[\mathcal{D}^{\gamma}\right]^{T}\left(\bigoplus_{\ell\in\gamma}\mathbf{\Psi}_{\ell}\right)\mathrm{vec}\left[\mathcal{D}^{\gamma}\right]}$
$\displaystyle=\frac{1}{2}\sum_{\gamma}\left(\mathrm{tr}\left[\mathbf{\Omega}^{\gamma}\mathrm{vec}\left[\mathcal{D}^{\gamma}\right]\mathrm{vec}\left[\mathcal{D}^{\gamma}\right]^{T}\right]-\log\mathrm{det}^{\dagger}\mathbf{\Omega}^{\gamma}+k_{\forall}^{\gamma}\log
2\pi\right)$
$\displaystyle=\frac{1}{2}\sum_{\gamma}\left(\mathrm{tr}\left[\mathbf{\Lambda}^{\gamma,(k_{\ell})}\left(\mathbf{V}^{\gamma,(k_{\ell}),T}\mathbf{S}^{\gamma}\mathbf{V}^{\gamma,(k_{\ell})}\right)\right]-\log\mathrm{det}\mathbf{\Lambda}^{\gamma,(k_{\ell})}+k_{\forall}^{\gamma}\log
2\pi\right)$

Note that
$\mathrm{tr}\left[\mathbf{\Lambda}^{\gamma,(k_{\ell})}\left(\mathbf{V}^{\gamma,(k_{\ell}),T}\mathbf{S}^{\gamma}\mathbf{V}^{\gamma,(k_{\ell})}\right)\right]$
is convex as it is linear, and that we no longer need to use the
pseudodeterminant as we have already kept only nonzero eigenvalues. The
negative log determinant is strictly convex. Thus our problem is the sum of
convex and strictly convex functions; this is strictly convex over the domain
of all $\mathbf{\Omega}^{\gamma}$.

It remains to show that, restricted to the domain of Kronecker-sum-
decomposable $\mathbf{\Omega}^{\gamma}$, the function is strictly convex. As
the set of KS-decomposable matrices is a linear (convex) subspace of the whole
set, strict convexity is preserved. This completes the proof. ∎

Note the nuance at the end; the NLL is not strictly convex over the space of
$\\{\mathbf{\Psi}_{\ell}\\}$. However, it is strictly convex over the space of
all Kronecker-sum-decomposable matrices, a space which could be (non-
identifiably) parameterized by $\\{\mathbf{\Psi}_{\ell}\\}$, but also has
identifiable parameterizations.

Putting together all of the existence, uniqueness, and identifiability
results, we can now state Theorem 14.

###### Theorem 14

There exists a unique MLE (up to the non-identifiability of the diagonals of
the Kronecker sum) that estimates the precision matrix of the singular
Kronecker-sum-structured normal distribution, provided each
$\mathbf{E}_{\ell}^{(k_{\ell})}$ contains only strictly positive eigenvalues
for each $\ell$, and there exists a strictly positive constant
$\epsilon_{\ell}$ such that
$\mathbf{\Lambda}_{\ell}^{(k_{\ell})}\geq\epsilon_{\ell}$ for every $\ell$.

###### Proof.

This follows directly from Lemmas 10, 11, and 13. ∎

###### Corollary 15

If we pick the number of components to keep $k_{\ell}$ by the formula
$k_{\ell}=\min\left(\mathrm{rank}\left[\mathbf{S}_{\ell}\right],\mathrm{rank}\left[\mathbf{\Psi}_{\ell}\right]\right)$,
we are guaranteed to have a unique solution (up to identifiability).

###  2.5 Hypothesis Testing

We are often interested in understanding how sure we can be in our results. In
this section, we will derive a formula for the p-value of each edge of our
estimated graph, under the null hypothesis of an empty graph. The method will
make use of the asymptotic normality of the MLE, which states roughly that,
when one has a lot of data, the error of any MLE follows a normal distribution
- this is the Wald Test.

$\displaystyle\hat{\mathbf{\Psi}}_{\ell}-\mathbf{\Psi}_{\ell}$
$\displaystyle\mathrel{\dot{\sim}}\mathcal{N}\left(0,\mathbf{F}^{-1}\right)$
$\displaystyle\mathbf{F}$
$\displaystyle=\mathbb{E}_{\mathbf{\Psi}_{\ell}}\left[\frac{\partial}{\partial\mathrm{vec}\left[\mathbf{\Psi}_{\ell}\right]\partial\mathrm{vec}\left[\mathbf{\Psi}_{\ell}\right]}\mathrm{NLL}\right]$

$\mathbf{F}$ is the Fisher information matrix. For simplicity of notation we
focused on a single axis, $\mathbf{\Psi}_{\ell}$, in the above equations. In
reality, all of the parameters of our model are jointly related by a single
information matrix $\mathbf{F}$. There are terms in $\mathbf{F}$ that
correspond to cross-links between different axes of the data, and links
between the axes and the trace terms discussed in Section 2.3. The task of
creating our hypothesis test boils down to the calculation of
$\mathbf{F}^{-1}$.

If we are not careful, $\mathbf{F}$ will end up singular. For example,
$\mathbf{\Psi}_{\ell}$ is not a matrix of $d_{\ell}^{2}$ independent
parameters, and if we treated it as such then there would be repeat rows in
$\mathbf{F}$. Taking into account symmetry, there are
$\frac{d_{\ell}^{2}-d_{\ell}}{2}$ off-diagonal independent parameters.
Further, because we constrain $\mathrm{tr}\left[\mathbf{\Psi}_{\ell}\right]=0$
to arrive at our identifiable parameterization, there are overall
$\frac{d_{\ell}^{2}-d_{\ell}}{2}+d_{\ell}-1$ parameters per axis. We will let
$\psi_{\ell_{11}}$ be the diagonal element that is removed from consideration.
The number of trace parameters $t_{i}$ is highly variable, as can be seen in
Table 2.

In general, the total number of parameters is
$\left|\left\\{t_{i}\right\\}_{i}\right|+\sum_{\ell}\frac{d_{\ell}^{2}-d_{\ell}}{2}+d_{\ell}-1$.
In fact, because in this paper we often consider non-full-rank matrices, the
number of parameters is actually smaller. However, this would lead to a far
more complicated situation for little gain; every low-rank matrix is an
arbitrarily small nudge away from a full-rank one. Thus, we will not factor in
the reduced degrees of freedom from the low-rank hypothesis; the results given
here should well-approximate the low-rank ones.

Calculating $\mathbf{F}$ requires the consideration of intra-axis, inter-axis,
inter-trace, and trace-axis derivatives. It is a multi-step, tedious process,
so we have delayed the derivation to Appendix A. All derivatives will be with
respect to the identifiable representation, not the standard representation,
and thus will be slightly different from those used in the algorithm. For
convenience, we repeat the negative log-likelihood of the identifiable
representation.

$\displaystyle\mathrm{NLL}$
$\displaystyle=\frac{-1}{2}\sum_{\gamma}\left(\log\mathrm{det}\left[\mathbf{t}^{T}\mathbf{a}^{\gamma}\mathbf{I}_{d_{\forall}^{\gamma}}+\bigoplus_{\ell\in\gamma}\mathbf{\Psi}_{\ell}\right]-\sum_{\ell\in\gamma}\mathrm{tr}\left[\mathbf{S}_{\ell}^{\gamma}\mathbf{\Psi}_{\ell}\right]-\sum_{\ell\in\gamma}\frac{\mathbf{t}^{T}\mathbf{a}^{\gamma}}{L^{\gamma}}\mathrm{tr}\left[\mathbf{S}^{\gamma}_{\ell}\right]\right)$

For a statistical test, we need a null hypothesis. Since we are interested in
testing whether an edge is statistically significant, our null hypothesis
should include that
$\mathrm{offdiag}\left[\mathbf{\Psi}_{\ell}\right]=\mathbf{0}$. We also need
to pick a value for the variance. It is not always possible to have unit
variance for every modality; this can be seen by considering the case
$\gamma_{1}=(\ell_{1},\ell_{2}),\gamma_{2}=(\ell_{3},\ell_{4}),\gamma_{3}=(\ell_{1},\ell_{2},\ell_{3},\ell_{4})$
given in Table 2. We will thus have the hypothesis that, for each modality
$\gamma$, the variance is on average $\left(\sigma^{\gamma}\right)^{2}$ for
some constant $\sigma^{\gamma}$. This is equivalent to assuming that
$\mathbf{t}^{T}\mathbf{a}^{\gamma}=\frac{1}{\left(\sigma^{\gamma}\right)^{2}}$,
or alternatively that
$\mathbf{\Omega}^{\gamma}=\frac{1}{\left(\sigma^{\gamma}\right)^{2}}\mathbf{I}_{d^{\gamma}_{\forall}}$.

As mentioned, there may be dependencies between $\sigma^{\gamma}$, restricting
the set of valid hypotheses. $\mathbf{a}^{\gamma}$ are constants that depend
on the dataset; for most datasets, they are either 1 or 0, but in theory they
can take any value. See Table 2.

Under this hypothesis, we can calculate $\mathbf{F}$. The calculations are
tedious, so we have deferred the proof to Appendix A.

###### Lemma 16

Under the null hypothesis where
$\mathbf{\Omega}^{\gamma}=\frac{1}{\left(\sigma^{\gamma}\right)^{2}}\mathbf{I}_{d^{\gamma}_{\forall}}$,
we have the following:

$\displaystyle\mathbf{F}_{t_{i}t_{j}}$
$\displaystyle=\sum_{\gamma}\frac{d^{\gamma}_{\forall}\left(\sigma^{\gamma}\right)^{4}}{2}\left(\mathbf{a}^{\gamma}\mathbf{a}^{\gamma,T}\right)_{ij}$
$\displaystyle\mathbf{F}_{t_{i}\psi_{\ell_{jk}}}$
$\displaystyle=\sum_{\gamma|\ell\in\gamma}\frac{a^{\gamma}_{i}d^{\gamma}_{\forall}\left(\sigma^{\gamma}\right)^{4}}{2d_{\ell}}\delta_{jk}$
$\displaystyle\mathbf{F}_{\psi_{\ell^{1}_{i^{1}j^{1}}}\psi_{\ell^{2}_{i^{2}j^{2}}}}$
$\displaystyle=\frac{1}{2}\sum_{\gamma|\ell^{1},\ell^{2}\in\gamma}\left(\sigma^{\gamma}\right)^{4}\frac{d_{\forall}^{\gamma}}{d_{\ell^{1}}d_{\ell^{2}}}\delta_{i^{2}j^{2}}\delta_{i^{1}j^{1}}$
($\ell^{1}\neq\ell^{2}$)
$\displaystyle\mathbf{F}_{\psi_{\ell^{1}_{i^{1}j^{1}}}\psi_{\ell^{2}_{i^{2}j^{2}}}}$
$\displaystyle=\frac{1}{2}\sum_{\gamma|\ell^{1},\ell^{2}\in\gamma}\left(\sigma^{\gamma}\right)^{4}\frac{d^{\gamma}_{\forall}}{d_{\ell}}\left(\delta_{i^{1}j^{1}i^{2}j^{2}}+(1-\delta_{i^{1}j^{1}i^{2}j^{2}})\frac{\delta_{i^{1}i^{2}}\delta_{j^{1}j^{2}}}{2}\right)$
($\ell^{1}=\ell^{2}$)

| $\mathbf{t}$ | Diagonal $\psi_{ii}$ | Off-diagonal $\psi_{ij}$
---|---|---|---
$\mathbf{t}$ | Yes333This depends on $\left\\{\mathbf{a}^{\gamma}\right\\}$. For most datasets, the most natural choice of identifiable parameterization results in the $\mathbf{t}$-$\mathbf{t}$ interactions forming a diagonal matrix (none of the $\mathbf{t}$ interact). | Yes | No
Diagonal $\psi_{ii}$ | Yes | Only in different axes | No
Off-diagonal $\psi_{ij}$ | No | No | No
Table 3: Allowed interactions between parameters in $\mathbf{F}$.

Recall that the only independent parameters are $\mathbf{t}$ and
$\psi_{\ell_{ij}}$ for $i\leq j$ and $(i,j)\neq(1,1)$. There are three types
of parameters here; $\mathbf{t}$, diagonal $\psi_{\ell_{ii}}$, and off-
diagonal $\psi_{\ell_{ij}}$. Of particular interest to us are the interactions
between off-diagonal elements $\psi_{\ell_{ij}}$ (where $\delta_{ij}=0$). It
is clear they cannot interact with $\mathbf{t}$ terms due to the $\delta_{ij}$
term in its formula. They also cannot interact with any
$\psi_{\ell^{\prime}_{i^{\prime}j^{\prime}}}$ in different axes, since the
$\delta_{ij}\delta_{i^{\prime}j^{\prime}}$ term is only nonzero for the
diagonal elements. Finally, the $\delta$ terms governing interactions within
the same axis are only nonzero when $(i,j)=(i^{\prime},j^{\prime})$ (self-
interactions) or $(i,j)=(j^{\prime},i^{\prime})$, which for off-diagonal
elements can never be true due to the $i\leq j$ requirement.

For all possible interactions, consult Table 3 \- these can be derived by
investigating the $\delta$ terms in the formula for $\mathbf{F}$. The fact
that off-diagonal elements do not interact with anything in $\mathbf{F}$ is
quite significant, as seen by the following corollary to Lemma 16.

###### Corollary 17

Suppose we order $\mathbf{F}$ such that the first rows/columns represent
$\mathbf{t}$ terms, the next represent diagonal $\psi_{ii}$ terms, and the
last represent the off-diagonal terms. Then, we have the following:

$\displaystyle\mathbf{F}$
$\displaystyle=\begin{bmatrix}\mathbf{A}&\mathbf{0}\\\
\mathbf{0}&\mathbf{B}\end{bmatrix}$

$\mathbf{B}$ represents the off-diagonal terms, and is a diagonal matrix. This
allows a convenient expression for $\mathbf{F}^{-1}$, by inverting blockwise.

We will now give the example of a dataset
$\gamma_{1}=(\ell_{1},\ell_{2}),\gamma_{2}=(\ell_{1},\ell_{3})$ where
$d_{\ell_{1}}=3$ and $d_{\ell_{2}}=d_{\ell_{3}}=2$. Note that
$\mathbf{a}^{\gamma_{i}}=\begin{bmatrix}\delta_{1i}&\delta_{2i}\end{bmatrix}$
in this case (and most cases). This leads to
$\mathbf{F}_{\mathbf{t},\mathbf{t}}$ being proportional to the identity
matrix. If we let
$A^{\gamma}=\frac{d^{\gamma}_{\forall}\left(\sigma^{\gamma}\right)^{4}}{2}$,
then it becomes easy to express the matrix. Exact derivations are left as an
exercise to the reader.

$\tiny\sum_{\gamma}A^{\gamma}\hskip
5.0pt\blockarray{cccccccccccc}t_{1}&t_{2}\begin{matrix}\ell_{1}\\\
(2,2)\end{matrix}\begin{matrix}\ell_{1}\\\
(3,3)\end{matrix}\begin{matrix}\ell_{2}\\\
(2,2)\end{matrix}\begin{matrix}\ell_{3}\\\
(2,2)\end{matrix}\begin{matrix}\ell_{1}\\\
(1,2)\end{matrix}\begin{matrix}\ell_{1}\\\
(1,3)\end{matrix}\begin{matrix}\ell_{1}\\\
(2,3)\end{matrix}\begin{matrix}\ell_{2}\\\
(1,2)\end{matrix}\begin{matrix}\ell_{3}\\\ (1,2)\end{matrix}\\\
\block{(ccccccccccc)c}10\frac{1}{d_{\ell_{1}}}\frac{1}{d_{\ell_{1}}}\frac{1}{d_{\ell_{2}}}\frac{1}{d_{\ell_{3}}}00000t_{1}\\\
01\frac{1}{d_{\ell_{1}}}\frac{1}{d_{\ell_{1}}}\frac{1}{d_{\ell_{2}}}\frac{1}{d_{\ell_{3}}}00000t_{2}\\\
\frac{1}{d_{\ell_{1}}}\frac{1}{d_{\ell_{1}}}\frac{1}{d_{\ell_{1}}}0\frac{1}{d_{\ell_{1}\ell_{2}}}\frac{1}{d_{\ell_{1}\ell_{3}}}00000\begin{matrix}\ell_{1}\\\
(2,2)\end{matrix}\\\
\frac{1}{d_{\ell_{1}}}\frac{1}{d_{\ell_{1}}}0\frac{1}{d_{\ell_{1}}}\frac{1}{d_{\ell_{1}\ell_{2}}}\frac{1}{d_{\ell_{1}\ell_{3}}}00000\begin{matrix}\ell_{1}\\\
(3,3)\end{matrix}\\\
\frac{1}{d_{\ell_{2}}}\frac{1}{d_{\ell_{2}}}\frac{1}{d_{\ell_{1}}d_{\ell_{2}}}\frac{1}{d_{\ell_{1}}d_{\ell_{2}}}\frac{1}{d_{\ell_{2}}}\frac{1}{d_{\ell_{2}\ell_{3}}}00000\begin{matrix}\ell_{2}\\\
(2,2)\end{matrix}\\\
\frac{1}{d_{\ell_{3}}}\frac{1}{d_{\ell_{3}}}\frac{1}{d_{\ell_{1}}d_{\ell_{3}}}\frac{1}{d_{\ell_{1}}d_{\ell_{3}}}\frac{1}{d_{\ell_{2}}d_{\ell_{3}}}\frac{1}{d_{\ell_{3}}}00000\begin{matrix}\ell_{3}\\\
(2,2)\end{matrix}\\\
000000\frac{1}{2d_{\ell_{1}}}0000\begin{matrix}\ell_{1}\\\
(1,2)\end{matrix}\\\
0000000\frac{1}{2d_{\ell_{1}}}000\begin{matrix}\ell_{1}\\\
(1,3)\end{matrix}\\\
00000000\frac{1}{2d_{\ell_{1}}}00\begin{matrix}\ell_{1}\\\
(2,3)\end{matrix}\\\
000000000\frac{1}{2d_{\ell_{2}}}0\begin{matrix}\ell_{2}\\\
(1,2)\end{matrix}\\\
0000000000\frac{1}{2d_{\ell_{3}}}\begin{matrix}\ell_{3}\\\
(1,2)\end{matrix}\\\ $

With everything in place, we can now define the hypothesis test.

###### Theorem 18

Under the null hypothesis
$\mathbf{\Omega}^{\gamma}=\frac{1}{\left(\sigma^{\gamma}\right)^{2}}\mathbf{I}_{d^{\gamma}_{\forall}}$,
we have the following distribution (where each sample is independent):

$\displaystyle\sqrt{\frac{\sum_{\gamma|\ell\in\gamma}d^{\gamma}_{\forall}\left(\sigma^{\gamma}\right)^{4}}{d_{\ell}}}\frac{\psi_{\ell_{ij}}}{2}$
$\displaystyle\mathrel{\dot{\sim}}\mathcal{N}\left(0,1\right)$

###### Proof.

As noted in Corollary 17, the edge-edge interactions in $\mathbf{F}$ form a
diagonal matrix, and there are no cross-terms connecting edges to diagonal and
$\mathbf{t}$ parameters. Thus, inversion is simple.

By Lemma 16, each edge-edge precision term for axis $\ell$ in $\mathbf{F}$ can
be expressed as
$\sum_{\gamma|\ell\in\gamma}\frac{d^{\gamma}_{\forall}\left(\sigma^{\gamma}\right)^{4}}{4d_{\ell}}$.
In $\mathbf{F}^{-1}$, this becomes
$\frac{4d_{\ell}}{\sum_{\gamma|\ell\in\gamma}d^{\gamma}_{\forall}\left(\sigma^{\gamma}\right)^{4}}$.

Thus, under the null hypothesis, we have that:

$\displaystyle\psi_{\ell_{ij}}$
$\displaystyle\mathrel{\dot{\sim}}\mathcal{N}\left(0,\frac{4d_{\ell}}{\sum_{\gamma|\ell\in\gamma}d^{\gamma}_{\forall}\left(\sigma^{\gamma}\right)^{4}}\right)$
$\displaystyle\sqrt{\frac{\sum_{\gamma|\ell\in\gamma}d^{\gamma}_{\forall}\left(\sigma^{\gamma}\right)^{4}}{d_{\ell}}}\frac{\psi_{\ell_{ij}}}{2}$
$\displaystyle\mathrel{\dot{\sim}}\mathcal{N}\left(0,1\right)$

∎

###### Corollary 19

The null hypothesis distribution for an edge estimate from a dataset of the
form $\gamma_{1}=(\ell_{1},...,\ell_{L})$ is:

$\displaystyle\sqrt{d_{\backslash\ell}}\sigma^{2}\frac{\psi_{\ell_{ij}}}{2}$
$\displaystyle\mathrel{\dot{\sim}}\mathcal{N}\left(0,1\right)$

###### Corollary 20

The null hypothesis distribution for an edge estimate from a dataset of the
form $\left\\{\gamma_{i}=(\ell_{0},\ell_{i})\right\\}$ is:

$\displaystyle\sqrt{\sum_{k}d_{\ell_{k}}\left(\sigma^{\gamma_{k}}\right)^{4}}\frac{\psi_{\ell_{ij}^{0}}}{2}$
$\displaystyle\mathrel{\dot{\sim}}\mathcal{N}\left(0,1\right)$
$\displaystyle\sqrt{d_{0}}\left(\sigma^{\gamma_{k}}\right)^{2}\frac{\psi_{\ell_{ij}^{k}}}{2}$
$\displaystyle\mathrel{\dot{\sim}}\mathcal{N}\left(0,1\right)$

If our null hypothesis has unit variance, the coefficient is exactly the
square root of the total number of samples available to the axis. This is the
same as when the data only has one axis, and thus is exactly what we would
have expected and hoped for.

Note that the Nonparanormal Skeptic maps data to have unit variance along each
axis. If one uses this, or some other method that standardizes the data, then
a unit-variance null hypothesis per axis is quite reasonable. This means the
precision matrix has diagonals $L^{\gamma}$ (since the terms are a sum of
terms in the composite axes). In this case,
$\left(\sigma^{\gamma}\right)^{2}=\frac{1}{L^{\gamma}}$ would be the correct
null hypothesis. The $\left(\sigma^{\gamma}\right)^{2}=1$ null hypothesis is
correct when one applies a standardization method that maps the whole dataset
to have unit variance. In practice, we have found this test has quite low
power, so thresholding may be preferred, especially if the graph is merely
used as a preprocessing step (such as for clustering) rather than an end in
and of itself.

###  2.6 Practical Implementation

1:$\\{\mathcal{D}^{\gamma}\\}$, $\\{\epsilon_{\ell}\\}$, $\\{n_{\ell}\\}$,
$\\{t_{\ell}\\}$, $\\{k_{\ell}\\}$

2:$\\{\mathbf{\Psi}_{\ell}\\}$

3:for $\ 1\leq\ell\leq K$

4:
$\mathbf{D}_{\ell}\leftarrow\begin{bmatrix}\mathrm{mat}_{\ell}\left[\mathcal{D}^{\gamma_{1}}\right]&\cdots&\mathrm{mat}_{\ell}\left[\mathcal{D}^{\gamma_{\Gamma}}\right]\end{bmatrix}$

5: $\mathbf{V}_{\ell}^{(k_{\ell})}\leftarrow\text{top }k_{\ell}\text{ left
singular vectors of }\mathbf{D}_{\ell}$ $\triangleright$ Theorem 4

6: $\mathbf{E}_{\ell}^{(k_{\ell})}\leftarrow\left(\text{top }k_{\ell}\text{
singular values of }\mathbf{D}_{\ell}\right)^{2}$

7:end for

8:$\mathbf{\Lambda}_{\ell}^{(k_{\ell})}\leftarrow\begin{bmatrix}\frac{1}{E_{\ell_{11}}^{(k_{\ell})}}&...&\frac{1}{E_{\ell_{k_{\ell}k_{\ell}}}^{(k_{\ell})}}\end{bmatrix}^{T}$

9:$\mu\leftarrow 1$

10:while not converged

11: for $\ 1\leq\ell\leq K$ $\triangleright$ Theorem 7

12:
$\mathbf{\Lambda}^{\prime(k_{\ell})}_{\ell}\leftarrow\frac{1}{2}\left(\mathbf{E}_{\ell}^{(k_{\ell})}-\sum_{\gamma|\ell\in\gamma}\mathrm{tr}^{k^{\gamma}_{<\ell}}_{k^{\gamma}_{>\ell}}\left[\left(\bigoplus_{\ell^{\prime}\in\gamma}\mathbf{\Lambda}_{\ell^{\prime}}^{(k_{\ell^{\prime}})}\right)^{-1}\right]\right)$

13: end for

14: for $\ \gamma\in\mathrm{modalities}$ $\triangleright$ Prevent overshooting

15: while
$\sum_{\ell\in\gamma}\mathrm{min}\mathbf{\Lambda}^{\prime(k_{\ell})}_{\ell}<\epsilon_{\ell}$

16: decrease $\mu$

17:
$\mathbf{\Lambda}^{\prime(k_{\ell})}_{\ell}\leftarrow\frac{1}{2}\left(\mathbf{E}_{\ell}^{(k_{\ell})}-\sum_{\gamma|\ell\in\gamma}\mathrm{tr}^{k^{\gamma}_{<\ell}}_{k^{\gamma}_{>\ell}}\left[\left(\bigoplus_{\ell^{\prime}\in\gamma}\mathbf{\Lambda}_{\ell^{\prime}}^{(k_{\ell^{\prime}})}\right)^{-1}\right]\right)$

18: end while

19: end for

20: for $\ 1\leq\ell\leq K$ $\triangleright$ Update parameters

21:
$\mathbf{\Lambda}^{(k_{\ell})}_{\ell}\leftarrow\mathbf{\Lambda}^{\prime(k_{\ell})}_{\ell}$

22: end for

23:end while

24:for $\ 1\leq\ell\leq K$

25:
$\mathbf{\Psi}_{\ell}\leftarrow\mathrm{thresh}_{n_{\ell}}\left[\mathbf{V}_{\ell}^{(k_{\ell})}\mathbf{\Lambda}_{\ell}^{(k_{\ell})}\mathbf{V}_{\ell}^{T,(k_{\ell})}\right]$
$\triangleright$ Threshold simultaneously

26:end for

Algorithm 1 The improved GmGM algorithm

We implemented the algorithm in Python, using Numpy, Numba, SciPy, and Dask
(Harris et al., 2020; Lam et al., 2015; Virtanen et al., 2020; Rocklin, 2015).
Pseudocode is provided in Algorithm 1. Theorems 4 and 7 form the backbone of
the algorithm. A sketch would be:

  1. 1.

Calculate the effective Gram matrices $\mathbf{S}_{\ell}$.

  2. 2.

(Partially) eigendecompose each $\mathbf{S}_{\ell}$ to get
$\mathbf{V}_{\ell}^{(k_{\ell})}$.

  3. 3.

Use update formula given by Theorem 7 until convergence, estimating
$\mathbf{\Lambda}_{\ell}$.

  4. 4.

‘Eigen-recompose’ to find the solution;
$\mathbf{\Psi}_{\ell}=\mathbf{V}_{\ell}^{(k_{\ell})}\mathbf{\Lambda}_{\ell}^{(k_{\ell})}\mathbf{V}_{\ell}^{(k_{\ell}),T}$.

  5. 5.

Threshold $\mathbf{\Psi}_{\ell}$ to enforce sparsity.

There are, however, several complicating factors. If we wish to avoid an
$O(\sum_{\ell}d_{\ell})$ memory cost, we can never calculate the full
$\mathbf{S}_{\ell}$ and $\mathbf{\Psi}_{\ell}$. The latter solution is
conceptually easy: we eigen-recompose and threshold simultaneously, only ever
storing values that surpass this threshold. The implementation of this would
depend on the exact type of thresholding required, but assuming it either
involves keeping the top $n_{\ell}$ edges overall or the top $n_{\ell}$ edges
per vertex, the computational complexity would be
$O(\sum_{\ell}d_{\ell}^{2}\log n_{\ell})$, with space complexity $O(n_{\ell})$
or $O(n_{\ell}d_{\ell})$ depending on the method. Its runtime is not drastic,
but enough to make up a significant portion, at around 20% of the runtime for
large problems (see Figure 4).

Avoiding calculation of $\mathbf{S}_{\ell}$ requires a bit more care. Our
strategy will be to directly calculate $\mathbf{V}_{\ell}^{(k_{\ell})}$ from
the data, using a (partial) singular value decomposition. For convenience, we
will assume that our modalities $\gamma$ are indexed from $\gamma=1$ to
$\gamma=\Gamma$.

$\displaystyle\mathbf{S}_{\ell}$
$\displaystyle=\sum_{\gamma}\mathbf{S}_{\ell}^{\gamma}$
$\displaystyle=\sum_{\gamma}\mathrm{mat}_{\ell}\left[\mathcal{D}^{\gamma}\right]\mathrm{mat}_{\ell}\left[\mathcal{D}^{\gamma}\right]^{T}$
$\displaystyle=\begin{bmatrix}\mathrm{mat}_{\ell}\left[\mathcal{D}^{1}\right]&...&\mathrm{mat}_{\ell}\left[\mathcal{D}^{\Gamma}\right]\end{bmatrix}\begin{bmatrix}\mathrm{mat}_{\ell}\left[\mathcal{D}^{1}\right]^{T}\\\
...\\\ \mathrm{mat}_{\ell}\left[\mathcal{D}^{\Gamma}\right]^{T}\end{bmatrix}$

Thus, $\mathbf{S}_{\ell}$ can be thought of as the Gram matrix for a
$d_{\ell}\times\sum_{\gamma}d^{\gamma}_{\backslash\ell}$ data matrix
“$\mathbf{D}_{\ell}$”. The left singular vectors of $\mathbf{D}_{\ell}$
correspond to the eigenvectors of $\mathbf{S}_{\ell}$. We can calculate this
directly from the input data, without any intermediate products.

Directly calculating the eigenvectors does have a drawback. The Nonparanormal
Skeptic, which is used to weaken the normality assumption, is a way of
estimating $\mathbf{S}_{\ell}$ \- if we want to both avoid quadratic memory
usage and weaken the normality assumption, it requires a modification. To
calculate nonparanormal eigenvectors, we use COCA (Han and Liu, 2014). COCA
also makes the nonparanormal assumption (the data follows a Gaussian copula
with arbitrary marginals). This is the same as replacing $\mathcal{D}_{\ell}$
with its ranks, mapping them to a normal distribution, and then calculating
the left singular vectors of this new dataset.

The addition of this (optional) step requires
$O(\max_{\ell}d_{\ell}\sum_{\gamma}d^{\gamma}_{\backslash\ell})$ space, as we
have to construct the matrix of ranks, and
$O\left(\sum_{\ell}\left(\sum_{\gamma}d^{\gamma}_{\backslash\ell}\right)d_{\ell}\log
d_{\ell}\right)$ time due to the cost of ranking the data.

Another difficulty arises due to the desire to take advantage of input
sparsity. It is often the case that transcriptomics and other -omics datasets
are highly sparse; the majority of the input is full of zeros. Ideally, we
would never ‘densify’ the matrix (explicitly store the zeros in memory). There
are many SVD routines for sparse data, so calculating the singular vectors
without the nonparanormal skeptic is as simple as using one of those routines.
However, the nonparanormal skeptic requires the calculation of ranks, and then
the mapping of those ranks to the normal distribution. This will almost always
require densification; sparse count-variate data, such as scRNA-seq data, will
always have 0 be the minimum value. It will thus have the smallest rank,
leading to it being mapped to the left tail of the Gaussian; whatever value it
achieves will certainly not be zero. Even if 0 is not the minimum value in the
dataset, there is no guarantee it will be mapped to the exact center of the
distribution after this process.

We will demonstrate the problem with an example, where $\Phi$ represents the
cumulative density function of the normal distribution.

$\displaystyle\begin{bmatrix}3&0&-2&0&4\\\ 0&-1&-3&0&-2\end{bmatrix}$
$\displaystyle\rightarrow\begin{bmatrix}4&2&1&2&5\\\ 4&3&1&4&2\end{bmatrix}$
(rank data)
$\displaystyle\rightarrow\begin{bmatrix}\Phi^{-1}\left(\frac{4}{6}\right)&\Phi^{-1}\left(\frac{2}{6}\right)&\Phi^{-1}\left(\frac{1}{6}\right)&\Phi^{-1}\left(\frac{2}{6}\right)&\Phi^{-1}\left(\frac{5}{6}\right)\\\
\Phi^{-1}\left(\frac{4}{6}\right)&\Phi^{-1}\left(\frac{3}{6}\right)&\Phi^{-1}\left(\frac{1}{6}\right)&\Phi^{-1}\left(\frac{4}{6}\right)&\Phi^{-1}\left(\frac{2}{6}\right)\end{bmatrix}$
(map to normal)
$\displaystyle\approx\begin{bmatrix}0.43&-0.43&-0.97&-0.43&0.97\\\
0.43&0&-0.97&0.43&-0.43\end{bmatrix}$

The exact values depend on how you handle ties when ranking data, but the
general idea remains the same regardless. The trick to avoiding the
densification is to realize that, in each row, all zero values get mapped to
the same final value. We can rewrite our example as:

$\displaystyle\begin{bmatrix}0.43+0.43&0&-0.97+0.43&0&0.97+0.43\\\
0&-0.43&-0.97-0.43&0&-0.43-0.43\end{bmatrix}+\begin{bmatrix}-0.43\\\
0.43\end{bmatrix}\mathbf{1}_{1\times 5}$

Let $z_{i}$ be the normally-distributed value that the zeros in row $i$ get
mapped to. We can express our transformation more generally:

$\displaystyle\begin{bmatrix}\mathbf{x}_{1}\\\ ...\\\
\mathbf{x}_{a}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\mathbf{x}_{1}-z_{1}\\\ ...\\\
\mathbf{x}_{a}-z_{a}\end{bmatrix}+\mathbf{z}\mathbf{1}^{T}$

Crucially, we have expressed our transformed data as the sum of a matrix with
the same sparsity structure as our original data, and a rank-one matrix. If we
compute the SVD on the first matrix, we can then use a rank-one update
algorithm to calculate the SVD of the sum of the two matrices (Brand, 2006).
It is not difficult to work out $\mathbf{z}$, and hence we can perform this
procedure without ever densifying our input matrix. This is necessary for
scalability; as we will see in Section 3.2, our algorithm can handle very
large sparse datasets; if densified, they would be in excess of 30 GB of
memory, but when sparse they fit in RAM.

Another decision that needs to be made is which edges of
$\mathbf{\Psi}_{\ell}$ to keep. One could keep the $n_{\ell}$ largest edges
per vertex, or the $n_{\ell}$ largest edges overall. We found that these
methods often over-focused on a subset of the vertices, to the detriment of
others. This could be counteracted by down-weighting edges of vertices with
high degree (by dividing by the total outgoing edge weight of the vertex).

The choice of thresholding can depend on the problem. In a graph of cells
(from a transcriptomics dataset), a common downstream task is to cluster them
and then assign cell types to the clusters. For this, we don’t want there to
be many, or any, singleton cells; they should all end up in some cluster.
Thresholding, just as with the L1 penalty, tends to lead to many singletons on
this type of data, so it may be worth keeping a minimum number of edges for
each vertex. For applications in which we expect many contaminants or
outliers, singletonhood may serve as a good marker for such properties; in
such cases, it would not be beneficial to keep a minimum number of edges per
vertex.

Ultimately, we found that overall thresholding, downweighting high-degree
vertices, and keeping a small number of minimum edges performed best on scRNA-
seq datasets. Alternatively, we can choose the parameters of this process
based on statistical significance, though the test suffers from low power.

The final challenge to implementation lies in enforcing the constraints
required by Theorem 14 to guarantee existence of a solution. During our
algorithm’s iterations, we must ensure
$\mathbf{\Lambda}_{\ell}^{(k_{\ell})}\geq\epsilon_{\ell}$. We get to choose
$\epsilon_{\ell}$, so if one is consistently running into this barrier, it
indicates that $\mathrm{rank}\left[\mathbf{\Psi}_{\ell}\right]<k_{\ell}$, and
we need to pick a smaller $k_{\ell}$. This is unlikely to happen. The
iterative portion of our algorithm takes negligible time on large samples (see
Figure 4, in which it takes less than 1% of the time when the problem size is
more than 1000 vertices).

###  2.7 Picking the Hyperparameters

Ignoring convergence-related hyperparameters, there are three hyperparameters:
how much to threshold, how many components to keep, and whether or not to use
the nonparanormal skeptic. The nonparanormal skeptic performs well enough that
the answer on whether or not to use it is almost always yes; it allows the
method to be effectively used on non-normal data, a very common scenario. The
main circumstance in which one might avoid it is for very large problems, as
you can squeeze a bit more scalability out of the method without it due to its
need to create an intermediate sparse matrix of similar size to the input.

The next choice is that of how to threshold. A good a priori choice is to keep
only the edges that are statistically significant ($p=0.05$) after the
Bonferroni correction. In practice, we have found that this test has low
power; for smaller datasets, no edges may be significant. When that is the
case, one needs to pick a thresholding method. As discussed in Section 2.6,
there are many types of thresholding methods.

We found overall thresholding (down-weighting vertices of high degree) tended
to perform well, but it is reasonable to try a range of techniques. We would
not recommend down-weighting vertices by the strength of the diagonal, as the
diagonals are non-identifiable. In practice, we found that keeping on average
10 edges per vertex resulted in good performance when used for clustering on
scRNA-seq data. When the graph is only wanted as a preprocessing step, rather
than as a goal in and of itself, it may be reasonable to also keep a minimum
number of edges for each vertex, to ensure for example that every vertex can
be assigned to a cluster. When the graph itself is the goal, choosing edges by
their statistical significance when possible is recommended.

We also need to pick an amount $k_{\ell}$ of eigenvectors to keep. In order to
satisfy Theorem 14 (and hence guarantee a minimum exists), we must have that
$k_{\ell}\leq\mathrm{rank}\left[\mathbf{S}_{\ell}\right]$, a known quantity.
We also must have that
$k_{\ell}\leq\mathrm{rank}\left[\mathbf{\Psi}_{\ell}\right]$, but in practice
this will be a larger upper bound. We will see in Section 3.2 that choosing
small values can still result in good performance on real data. As mentioned
in Corollary 5, the user can use the same techniques as in PCA to determine
the number of eigenvalues to keep, such as by looking at the total explained
variance, looking at a scree plot, or other methods.

##  3 Results

The experiments were run on a 13-inch MacBook Pro with an M1 chip and 8GB RAM.
23GB of hard disk storage were available when running the experiments; the
million cell experiment (Section 3.2.3) could not fit in 8GB of RAM, and thus
much of the free hard disk space was used as swap memory. We used Apple’s
Accelerate library for its LAPACK and BLAS methods. The code will still work
on other architectures; to squeeze out the most scalability, we would
recommend using a BLAS library tuned for the user’s system. Code for all
experiments is available at https://github.com/BaileyAndrew/gmgm-jmlr.

###  3.1 Synthetic Data

There are several decisions to make when testing on synthetic data: the
distribution of the ground truth graph for each axis, the way to combine the
axes, the distribution of the dataset derived from the graphs, and the manner
in which we evaluate performance. We opted to test our model on data that had
a Kronecker sum structure and was normally distributed; these were the
assumptions made by our model and prior work. Later, in Section 3.2, we will
test performance on real data, which will not fit these assumptions

We generated ground truth graphs from the Barabasi-Albert distribution. To
sample from the distribution, we used python-igraph’s implementation (Csardi
and Nepusz, 2005) with default parameters (namely, the power and zero_appeal
parameters were set to 1, corresponding with the standard Barabasi-Albert
distribution). We chose this distribution because it is a scale-free
distribution; many real-world networks, such as social networks and gene
regulatory networks, are approximately scale-free, and thus performance on
scale-free networks is of significant interest. Unlike our assumption of real-
world data being fundamentally low-rank, the Barabasi-Albert model’s rank will
grow with the size of the graph we generate. In general, it is very hard to
generate synthetic data that is both low-rank and sparse. Thus, the synthetic
data will likely violate our low-rank assumption more than real world data
will (but it will violate the other assumptions, such as normality, less).

To test performance, we report precision-recall (PR) curves for TeraLasso,
GmGM, and our proposed modifications to GmGM; when an algorithm correctly
identifies an edge, it is considered a true positive. When comparing to prior
work, we used their default convergence tolerance and max iterations.

(a) Precision-recall curves for a small (250-node) graph.

(b) Precision-recall curves for a medium-sized (5000-node) graph. TeraLasso
was omitted because it could not run in reasonable time on a graph of this
size.

Figure 2: Precision-recall curves on graphs of various sizes. On the right, we
report how much variance the used eigenvalues account for. For GmGM and
TeraLasso, this will always be 1, as they do not make a low-rank assumption.
The shaded region around each PR curve corresponds to the maximum and minimum
values reported over 10 runs; the central curve is the mean of the maximum and
minimum.

We report our PR curves in Figure 2. We can see that, while our variant does
worse than GmGM and TeraLasso, this matters less for larger graphs; as the
graph size increases, so does our algorithm’s performance. As our algorithm
was designed for graphs with hundreds of thousands of nodes, not just the 5000
nodes reported in Figure 2(b), this bodes well. We see that this improvement
in performance is maintained even when the explained variance drops; GmGM with
125 principal components in a 5000-node graph performs similarly to GmGM with
125 principal components in a 250-node graph, despite the former representing
10% of the explained variance and the latter representing more than 90%. With
just 40% of the explained variance on a 5000-node graph, one can perform as
well as using 100% of the variance on a 250-node graph.

Figure 3: Runtimes of various algorithms as the number of nodes in the graph
increases. In the top-left, we have focused in on the sub-1000-node region to
be able to show TeraLasso and DNNLasso. GmGM-50pcs-minimal corresponds to the
case in which we assume the number of edges in the graph is the same as the
number of nodes; for the other models, we kept the full graph.

If PR curve performance was all that mattered, prior work would have an edge
over ours. However, our algorithm’s speed more than makes up for this loss of
accuracy in high dimensions. In Figure 3, one can see the substantial
improvement in runtime that comes as a result of this tradeoff. One can also
see the cost of generating data, which takes about a minute once the problem
size reaches 5000 nodes. The largest graph tested, 10000 nodes, took more than
seven minutes to generate. This was the reason we created PR curves in Figure
2(b) at 5000 nodes rather than million nodes we aim to work with on real data.

Figure 4: The GmGM algorithm can be split into three parts; calculating the
eigenvectors, calculating the eigenvalues, and then combining these two
(‘recomposition’) to produce the final precision matrix. This plot shows the
relative amount of time spent in each of the three parts.

Finally, in Figure 4, we show the relative cost of each component that goes
into the algorithm. While the final step of ‘eigen-recomposition’ takes 20% of
the runtime for large graphs, the rest of the runtime is spent entirely on the
single partial singular value decomposition that the algorithm makes to
calculate the eigenvectors. Thus, to handle even larger datasets it is
necessary to find an SVD-free approach.

###  3.2 Real Data

All experiments on real world data used the nonparanormal skeptic as a
preprocessing step. Statistical significance is always reported after applying
the Bonferroni correction.

####  3.2.1 COIL Dataset

Figure 5: A comparison of several methods to find the frame graph. Connections | Adjacent | One off | Two off | Three off | More | Total |  Method
---|---|---|---|---|---|---|---
GmGM | 72 | 72 | 25 | 4 | 0 | 173 |  p=0.05
GmGM (10pcs) | 72 | 64 | 13 | 0 | 51 | 200 |  p=0.05
GmGM (10pcs, one modality) | 66 | 58 | 44 | 13 | 19 | 200 |  Thresholding
TeraLasso (one modality) | 33 | 31 | 26 | 22 | 85 | 197444We were not able to get exactly 200 edges |  Thresholding
Table 4: The performance of several methods on reconstructing the frame graph of COIL-20 videos. The adjacent column denotes how many edges were between adjacent frames (the ‘correct’ edges), whereas the later columns denote how many edges exist between frames $i$ and $i\pm(n+1)$ for various $n$. The More column summarizes all connections where $n>3$. Connections | $\begin{matrix}\text{Precision/Recall}\\\ \text{(exact)}\end{matrix}$ | $\begin{matrix}\text{Precision/Recall}\\\ \text{(one off)}\end{matrix}$ | $\begin{matrix}\text{Precision/Recall}\\\ \text{(two off)}\end{matrix}$ | $\begin{matrix}\text{Precision/Recall}\\\ \text{(three off)}\end{matrix}$
---|---|---|---|---
GmGM |  41.6%/100% |  83.2%/100% |  97.6%/78.2% |  100%/60%
GmGM (10pcs) | 36%/100% | 60%/94.4% | 74.5%/69% | 74.5%/51.7%
GmGM (10pcs, one modality) | 33%/91.7% | 62%/86% | 84%/77.8% | 90.5%/62.8%
TeraLasso (one modality) | 16.8%/45.8% | 32.5%/44.4% | 45.7%/41.7% | 56.7%/38.9%
Table 5: The precisions and recalls for each algorithm, when we consider
correct edges to be those between edges $i$ and $i\pm(n+1)$, for $n=0$
(‘exact’), $n=1$ (‘one off’), $n=2$ (‘two off’), and $n=3$ (‘three off’). In
each column, the best performing model’s results were given in bold for each
metric.

It is difficult to test on real-world data, as ground truth graphs are often
unknown. This is perhaps why the COIL-20 video dataset (Nene et al., ) has
become a standard test for multi-axis models.

The dataset consists of 20 videos of 72 frames and 128x128 pixels. Each video
consists of a simple object rotating 360 degrees, and thus the frame graph can
reasonably expected to correspond to this circular structure. To test this, we
considered the dataset to have the form
$\left\\{\gamma_{i}=(\ell_{\mathrm{frame}},\ell_{\mathrm{row}(i)},\ell_{\mathrm{col}(i)})\right\\}_{i}$,
that is, that there is a singular frame graph describing the rotational
structure of the images, and then individual row and column graphs for each
object (since the structure of each object may be different). We compared the
results using all eigenvectors and using the top 10.

In the aforementioned experiments, we kept only the edges that were
statistically significant to the 5% level after applying the Bonferroni
correction. To compare with prior work (TeraLasso), we also ran our algorithm
on just a single object (the duck object, “obj1”). No edges were statistically
significant, so instead we kept the same amount of edges as were found to be
significant when considering the whole dataset. We then did the same with
TeraLasso.

The results can be compared in Figure 5. These show that GmGM on the whole
dataset does the best; TeraLasso cannot handle data of this form. On data that
TeraLasso can handle (a single modality) our method still arguably does better
than it, even when restricted to only a rank-10 approximation. Numerical
values are given in Tables 4 and 5; we can see that, when considering
statistically significant edges in the whole dataset, our method finds all the
correct edges. The false positives it finds are typically close-to-correct,
which is reasonable.

Figure 6: A comparison of the object graphs produced by GmGM and TeraLasso
(with same sparsity level), when the COIL-20 dataset is taken as a
$(\ell_{\mathrm{object}},\ell_{\mathrm{frame}},\ell_{\mathrm{row}},\ell_{\mathrm{col}})$
tensor. Connections to non-adjacent objects are dashed; this is only to aid
visualization, as the objects have no natural adjacencies.

As our method and TeraLasso can handle order-4 tensors, it is tempting to
consider the entire dataset a $20\times 72\times 128\times 128$
$(\ell_{\mathrm{object}},\ell_{\mathrm{frame}},\ell_{\mathrm{row}},\ell_{\mathrm{col}})$
tensor. We could then look at the object graph. It is hard to have an a priori
guess on what this graph should look like; perhaps objects with similar shapes
or color schemes should be connected? We report our results, for both GmGM and
TeraLasso, in Figure 6, and leave their interpretation up to the reader. No
connection was statistically significant; we kept on average 3 edges per
vertex, and chose a corresponding regularization parameter of TeraLasso to
achieve the same sparsity.

####  3.2.2 Head and Neck Single-Cell Transcriptomics Dataset

Figure 7: The performance of various algorithmic configurations on various
subsets of genes. The data came from a multi-patient head-and-neck cancer
dataset (Cillo et al., 2020); the last two examples were created on a single
patient (patient HD_6) within the full dataset.

To test how our algorithm performs on real world data, we used a 150,000 cell
multi-patient dataset characterizing the immune landscape of head and neck
cancer (Cillo et al., 2020). We followed the same pre-processing steps in the
paper - namely keeping only genes expressed in more than 1% of cells, and
cells that both expressed more than 200 unique genes while having a total gene
expression of less than 20,000. After preprocessing, we were left with a
2200-gene dataset of 155,970 cells.

We first ran our algorithm keeping only the top 1000 principal components and
statistically significant edges (p=0.05). This corresponded to keeping more
than 90% of the variance. This resulted in a gene graph with 1048 edges. This
was very sparse, so we also present results in which we kept on average 10
edges per cell by thresholding, leading to a gene graph with 10995 edges. It
is impossible to compare our results to prior work on the whole dataset, as
TeraLasso is not scalable enough to run on it. Thus, we selected the patient
with the least cells, patient HD_6 (with 1,466 cells). We ran our algorithm on
this subset, also keeping on average 10 edges per cell, and ran TeraLasso as
well, with a regularization parameter chosen to mach the sparsity level of the
other graphs. No edges were statistically significant when we restricted
ourselves to the smaller dataset. Some genes in the dataset had multiple
alleles; for simplicity, we combined all alleles of each gene.

To validate our graphs, we investigated the structure of specific sub-graphs;
small ribosomal protein-coding genes, as well as four immune-related GO terms.
Ribosomal proteins work together to create the ribosome, and are highly
expressed; thus, we should expect there to be many connections in this subset.
GO terms group genes together by functional/biological significance, and thus
we can often expect relations between the genes as well. The subgraphs are
presented in Figure 7.

Algorithm | Cluster | Biological Process | p-value
---|---|---|---
| $\begin{matrix}\text{1}\\\ \text{(159 genes)}\end{matrix}$ | $\begin{matrix}\text{regulation of RNA splicing}\\\ \text{regulation of mRNA processing}\\\ \text{response to endoplasmic reticulum stress}\end{matrix}$ | $\begin{matrix}2.1\times 10^{-6}\\\ 7.7\times 10^{-6}\\\ 1.6\times 10^{-5}\end{matrix}$
| $\begin{matrix}\text{2}\\\ \text{(48 genes)}\end{matrix}$ | $\begin{matrix}\text{cytoplasmic translation}\\\ \text{translation}\\\ \text{peptide biosynthetic process}\end{matrix}$ | $\begin{matrix}7.2\times 10^{-76}\\\ 3.7\times 10^{-48}\\\ 1.8\times 10^{-47}\end{matrix}$
GmGM (p=0.05) | $\begin{matrix}\text{3}\\\ \text{(44 genes)}\end{matrix}$ | $\begin{matrix}\text{\color[rgb]{0,0,1}immune response}\\\ \text{\color[rgb]{0,0,1}immune system response}\\\ \text{\color[rgb]{0,0,1}leukocyte activation}\end{matrix}$ | $\begin{matrix}6.8\times 10^{-17}\\\ 2.9\times 10^{-15}\\\ 1.2\times 10^{-11}\end{matrix}$
| $\begin{matrix}\text{4}\\\ \text{(35 genes)}\end{matrix}$ | $\begin{matrix}\text{cytoplasmic translation}\\\ \text{translation}\\\ \text{peptide biosynthetic process}\end{matrix}$ | $\begin{matrix}1.8\times 10^{-17}\\\ 5.3\times 10^{-11}\\\ 8.7\times 10^{-11}\end{matrix}$
| $\begin{matrix}\text{5}\\\ \text{(31 genes)}\end{matrix}$ | $\begin{matrix}\text{\color[rgb]{0,0,1}GO:0002504}\\\ \text{\color[rgb]{0,0,1}GO:0002399}\\\ \text{\color[rgb]{0,0,1}GO:0002503}\end{matrix}$ | $\begin{matrix}4.7\times 10^{-13}\\\ 2.2\times 10^{-12}\\\ 2.2\times 10^{-12}\end{matrix}$
| $\begin{matrix}\text{1}\\\ \text{(1008 genes)}\end{matrix}$ | $\begin{matrix}\text{organonitrogen compound metabolic process}\\\ \text{oxidative phosphorylation}\\\ \text{protein metabolic process}\end{matrix}$ | $\begin{matrix}5.5\times 10^{-31}\\\ 1.2\times 10^{-23}\\\ 3.1\times 10^{-22}\end{matrix}$
| $\begin{matrix}\text{2}\\\ \text{(214 genes)}\end{matrix}$ | $\begin{matrix}\text{\color[rgb]{0,0,1}regulation of immune system process}\\\ \text{\color[rgb]{0,0,1}regulation of immune response}\\\ \text{\color[rgb]{0,0,1}positive regulation of immune response}\end{matrix}$ | $\begin{matrix}7.3\times 10^{-9}\\\ 1.4\times 10^{-8}\\\ 1.1\times 10^{-7}\end{matrix}$
TeraLasso (HD_6) | $\begin{matrix}\text{3}\\\ \text{(197 genes)}\end{matrix}$ | $\begin{matrix}\text{organonitrogen compound metabolic process}\\\ \text{aerobic electron transport chain}\\\ \text{mitochondrial ATP synthesis coupled electron transport}\end{matrix}$ | $\begin{matrix}4.7\times 10^{-6}\\\ 7.9\times 10^{-6}\\\ 1.8\times 10^{-5}\end{matrix}$
| $\begin{matrix}\text{4}\\\ \text{(114 genes)}\end{matrix}$ | $\begin{matrix}\text{\color[rgb]{0,0,1}defense response}\\\ \text{\color[rgb]{0,0,1}cytokine-mediated signalling pathway}\\\ \text{\color[rgb]{0,0,1}positive regulation of immune system process}\end{matrix}$ | $\begin{matrix}5.0\times 10^{-9}\\\ 5.8\times 10^{-9}\\\ 3.8\times 10^{-8}\end{matrix}$
| $\begin{matrix}\text{5}\\\ \text{(103 genes)}\end{matrix}$ | $\begin{matrix}\text{synapse pruning}\\\ \text{\color[rgb]{.5,.5,.5}N/A}\\\ \text{\color[rgb]{.5,.5,.5}N/A}\end{matrix}$ | $\begin{matrix}3.0\times 10^{-2}\\\ \color[rgb]{.5,.5,.5}>0.05\\\ \color[rgb]{.5,.5,.5}>0.05\end{matrix}$

Table 6: GO terms with long names have been replaced with their GO ID; these
were all major-histocompatibility-complex-associated GO terms. Immune-related
GO terms have been colored blue. These were identified using the Python API of
the GProfiler (Kolberg et al., 2023) functional enrichment tool.

The use of GO terms help us validate the local structure of our learned graph,
but the choice of GO terms is somewhat arbitrary and it does not tell us much
about the validity of the global structure of the graph. We used the Leiden
clustering algorithm (Traag et al., 2019) to investigate to what extent the
gene networks could be partitioned into coherent regions (‘modules’). Some of
the modules are quite large, so we do not display them here (although plots of
them are available in the repository for this paper). Rather, we report the
most associated GO terms in the biological process GO category for the 5
largest clusters; these are given in Table 6. For space reasons, we only
report the analyses for GmGM with statistically significant edges on the full
dataset and for TeraLasso on patient HD_6; the results for other algorithmic
setups are also available in our repository.

We would expect that a good graph neatly partitions into GO terms, and that
said GO terms should often relate to immunological properties. Both seem to
identify modules that are reasonably associated with coherent communities of
genes. The advantage of GmGM comes from the fact that it can handle the whole
dataset, whereas TeraLasso had to be restricted to a subset with 1,466 cells
(1% of the total). Thus, it can be expected to have captured processes that
hold true across all patients.

####  3.2.3 Million-Cell PBMC Transcriptomics Dataset

Figure 8: The performance of our algorithm on a million-cell dataset (Yazar et
al., 2022), when keeping only statistically significant edges. The sparsity of
the overall network is 0.2%.

At the start of this paper, we claimed our algorithm could run on million-cell
datasets. To prove this, we used a 1,248,980 cell, 36,571 gene PBMC dataset by
Yazar et al. (2022). For preprocessing, we removed any cells and genes that
had 0 counts everywhere. Then, we kept only the top 2,200 highly variable
genes555The value of 2,200 is was chosen because it was the same amount of
genes kept in our previous experiment in Section 3.2.2.. No cells had
uniformly 0 counts, so the post-preprocessing dataset size contained 1,248,980
cells and 2200 genes.

As in previous cases, we used the nonparanormal skeptic and kept only
statistically significant connections. Due to the size of our dataset, we kept
only the top 50 principal components (which corresponded to 50% of the cell
variance and 56% of the gene variance).

To validate the graph, we used the same methods as in Section 3.2.2. For local
structure validation, we used the same GO terms. The grand majority of the
ribosomal genes did not make it into the top 2200 highly variable genes, so
the graph is not particularly informative. Overall, 5890 edges were
statistically significant out of the approximately 2.4 million possible edges;
in other words, 0.2% of all possible edges are actualized in the resultant
network. In a group of 50 randomly selected genes, we would expect on average
only one edge to exist. Our method finds many more than this on the
immunology-related GO terms, see Figure 8.

Algorithm | Cluster | Biological Process | p-value
---|---|---|---
| $\begin{matrix}\text{1}\\\ \text{(67 genes)}\end{matrix}$ | $\begin{matrix}\text{cell-cell adhesion}\\\ \text{apoptotic process}\\\ \text{positive regulation of cellular process}\end{matrix}$ | $\begin{matrix}1.0\times 10^{-9}\\\ 3.1\times 10^{-9}\\\ 5.1\times 10^{-9}\end{matrix}$
| $\begin{matrix}\text{2}\\\ \text{(65 genes)}\end{matrix}$ | $\begin{matrix}\text{protein maturation}\\\ \text{response to endoplasmic reticulum stress}\\\ \text{organonitrogen compound metabolic process}\end{matrix}$ | $\begin{matrix}5.6\times 10^{-7}\\\ 1.3\times 10^{-6}\\\ 9.1\times 10^{-6}\end{matrix}$
GmGM (p=0.05) | $\begin{matrix}\text{3}\\\ \text{(60 genes)}\end{matrix}$ | $\begin{matrix}\text{\color[rgb]{0,0,1}immune response}\\\ \text{\color[rgb]{0,0,1}leukocyte mediated immunity}\\\ \text{\color[rgb]{0,0,1}immune system process}\end{matrix}$ | $\begin{matrix}9.5\times 10^{-15}\\\ 3.8\times 10^{-14}\\\ 8.7\times 10^{-14}\end{matrix}$
| $\begin{matrix}\text{4}\\\ \text{(43 genes)}\end{matrix}$ | $\begin{matrix}\text{GO:0048523}\\\ \text{GO:0051172}\\\ \text{GO:0031324}\end{matrix}$ | $\begin{matrix}4.7\times 10^{-5}\\\ 6.8\times 10^{-5}\\\ 8.9\times 10^{-5}\end{matrix}$
| $\begin{matrix}\text{5}\\\ \text{(27 genes)}\end{matrix}$ | $\begin{matrix}\text{\color[rgb]{0,0,1}positive regulation of immune system process}\\\ \text{\color[rgb]{0,0,1}regulation of immune response}\\\ \text{\color[rgb]{0,0,1}GO:0002504}\end{matrix}$ | $\begin{matrix}8.3\times 10^{-15}\\\ 2.3\times 10^{-14}\\\ 2.4\times 10^{-14}\end{matrix}$
| $\begin{matrix}\text{6}\\\ \text{(14 genes)}\end{matrix}$ | $\begin{matrix}\text{intracellular sequestering of iron ion}\\\ \text{sequestering of iron ion}\\\ \text{negative regulation of actin filament polymerization}\end{matrix}$ | $\begin{matrix}8.1\times 10^{-3}\\\ 1.1\times 10^{-2}\\\ 1.3\times 10^{-2}\end{matrix}$

Table 7: GO terms with long names have been replaced with their GO ID; most
were of the form “negative regulation of cellular/nitrogen compound
(metabolic) process”. Immune-related GO terms have been colored blue. These
were identified using the Python API of the GProfiler (Kolberg et al., 2023)
functional enrichment tool.

Finally, we verify the global structure of our graph by exploring the
clustering. Again, we used Leiden clustering. Of our 2200 genes, 1924 were
singletons in our graph. Ignoring them, we had 6 clusters; we report each of
their top 3 associated biological process GO terms in Table 7.

###  3.3 Computational Complexity

Calculating the top $k_{\ell}$ singular vectors for each axis’s
$d_{\ell}\times d_{\forall\backslash\ell}$ matrix $\mathbf{D}_{\ell}$ takes
$O(\sum_{\ell}k_{\ell}d_{\forall|\ell})$ time and
$O(\sum_{\ell}k_{\ell}d_{\ell})$ space. Per iteration, the convex optimization
step requires computing
$\left(\bigoplus_{\ell\in\gamma}\mathbf{\Lambda}_{\ell}^{(k_{\ell})}\right)^{-1}$,
which takes $O(\sum_{\gamma}\prod_{\ell\in\gamma}k_{\ell})$ time and space.
The eigen-recomposition step takes $O(\sum_{\ell}k_{\ell}d_{\ell}^{2})$ time,
with thresholding talking $O(\sum d_{\ell}^{2})$ time666This is true
regardless of whether we threshold overall, threshold per row, weight before
thresholding, or threshold based on statistical significance. \- although our
implementation when not thresholding by significance takes
$O(\sum_{\ell}d_{\ell}^{2}\log n_{\ell})$ to minimize the number of times we
pass through the data. As we threshold and eigen-recompose simultaneously, it
requires $O(n_{\ell})$ space.

The nonparanormal skeptic requires ranking the data, taking
$O(\sum_{\ell}d_{\forall|\ell}\log d_{\ell})$ time and $O(sd_{\forall})$
space, where $s$ is the sparsity of the input dataset. This leads to an
algorithm with the following complexity:

$\displaystyle
O\left(\sum_{\gamma}\prod_{\ell\in\gamma}k_{\ell}+\sum_{\ell}\left(k_{\ell}d_{\forall|\ell}+k_{\ell}d_{\ell}^{2}{\color[rgb]{.5,.5,.5}\log
n_{\ell}}+{\color[rgb]{0,0,1}d_{\forall|\ell}\log d_{\ell}}\right)\right)$
(time) $\displaystyle
O\left(\sum_{\gamma}\prod_{\ell\in\gamma}k_{\ell}+\sum_{\ell}\left(k_{\ell}d_{\ell}+n_{\ell}\right)+{\color[rgb]{0,0,1}sd_{\forall}}\right)$
(space)

Blue terms arise only if using the nonparanormal skeptic, and the grey term is
a product of our implementation; it could be dropped, but would in practice
make the algorithm slower. As this formula is complicated, we will give an
example of a special case. In the case of matrix-variate data with
$n_{\ell}\in O(d_{\ell})$ and constant $k_{\ell}$ (the case for many real-
world datasets), as well as for simplicity assuming $d=d_{\ell_{1}}\approx
d_{\ell_{2}}$ and thresholding by statistical significance, these become:

$\displaystyle O\left(d^{2}+{\color[rgb]{0,0,1}d^{2}\log d}\right)$ (time)
$\displaystyle O\left(d+{\color[rgb]{0,0,1}sd^{2}}\right)$ (space)

##  4 Conclusion

We have created an algorithm to solve the problem of (conditional dependency)
graph inference without making an independence assumption. Our method is
orders of magnitude faster than prior work on matrix-variate data, and is
memory-scalable as well; it can handle datasets with more than a million
samples on a personal computer. We have made these improvements while still
being applicable to the multi-modality, tensor-variate case, whereas most
prior work is limited to considering uni-modal matrix-variate data. To
validate our results, we have shown our performance on synthetic data, as well
as three real-world datasets. We have accompanied experimental validation with
rigorous proofs as to the existence and uniqueness of the results. Our
hyperparameters are interpretable, as the low-rank approximation can be
analyzed in terms of explained variance and the thresholding parameter can be
picked through a statistical hypothesis test that we derived in this paper.

In future work, we aim to generalize the model further. For example, our
framework cannot handle the case of multi-modality datasets with only partial
overlap of axes (that is, when some samples have both modalities while others
have only one modality), nor can any existing multi-axis algorithm in the
literature. We also aim to relax some of the assumptions of the model, namely
the Gaussian copula assumption. The prospect of an efficient
eigendecomposition-free algorithm is tantalizing, but currently out of reach.
If achieved, it would likely allow further scalability - more importantly, it
would allow scalability in the presence of regularization.

Acknowledgments and Disclosure of Funding

We would like to acknowledge the many people and animals who have donated
their time and data for the creation of scRNA-seq, and other, datasets. We
have no funding or competing interests to declare.

##  Appendix A Derivation of the Fisher Information

This section contains the derivation of the Fisher Information, which is too
long and tedious to put in the main paper. For convenience, we repeat the
definition of the $\mathrm{NLL}$.

$\displaystyle\mathrm{NLL}$
$\displaystyle=\frac{-1}{2}\sum_{\gamma}\left(\log\mathrm{det}\left[\mathbf{t}^{T}\mathbf{a}^{\gamma}\mathbf{I}_{d_{\forall}^{\gamma}}+\bigoplus_{\ell\in\gamma}\mathbf{\Psi}_{\ell}\right]-\sum_{\ell\in\gamma}\mathrm{tr}\left[\mathbf{S}_{\ell}^{\gamma}\mathbf{\Psi}_{\ell}\right]-\sum_{\ell\in\gamma}\frac{\mathbf{t}^{T}\mathbf{a}^{\gamma}}{L^{\gamma}}\mathrm{tr}\left[\mathbf{S}^{\gamma}_{\ell}\right]\right)$

Lemmas 21, 22, and 23 contain the calculations that derive the Hessian of the
$\mathrm{NLL}$.

###### Lemma 21

$\frac{\partial}{\partial\psi_{\ell_{ij}}}\mathrm{NLL}=-\frac{1}{2}\mathrm{sym}\left[\sum_{\gamma|\ell\in\gamma}\mathrm{tr}^{d^{\gamma}_{<\ell}}_{d^{\gamma}_{>\ell}}\left[\left(\mathbf{t}^{T}\mathbf{a}^{\gamma}\mathbf{I_{d^{\gamma}_{\forall}}}+\bigoplus_{\ell\in\gamma}\mathbf{\Psi}_{\ell}\right)^{-1}\right]_{ij}-\mathbf{S}_{\ell_{ij}}\right]$

###### Proof.

The derivation is nearly identical to that of Lemma 3. ∎

###### Lemma 22

$\displaystyle\frac{\partial}{\partial t_{i}}\mathrm{NLL}$
$\displaystyle=\sum_{\gamma}-\frac{a^{\gamma}_{i}}{2}\mathrm{tr}\left[\left(\mathbf{t}^{T}\mathbf{a}^{\gamma}\mathbf{I}_{d^{\gamma}_{\forall}}+\bigoplus_{\ell\in\gamma}\mathbf{\Psi}_{\ell}\right)^{-1}\right]+\frac{a^{\gamma}_{i}}{2}\mathrm{tr}\left[\mathbf{S}^{\gamma}\right]$
$\displaystyle\frac{\partial}{\partial t_{i}\partial t_{j}}\mathrm{NLL}$
$\displaystyle=\sum_{\gamma}\frac{a^{\gamma}_{i}a^{\gamma}_{j}}{2}\mathrm{tr}\left[\left(\mathbf{t}^{T}\mathbf{a}^{\gamma}\mathbf{I}_{d^{\gamma}_{\forall}}+\bigoplus_{\ell\in\gamma}\mathbf{\Psi}_{\ell}\right)^{-2}\right]$
$\displaystyle\frac{\partial}{\partial
t_{i}\partial\psi_{\ell_{jk}}}\mathrm{NLL}$