Learning Product Graphs Underlying
Smooth Graph Signals

The existing works in graph signal processing can mainly be divided into four chronological categories. The first set of works in GSP introduced the idea of information processing over graphs [4, 2, 5]. These works highlight the advantages and superior performance of graph-based signal processing approach (with known underlying graph) over classical signal processing. The second wave of research in this area built upon the first one to exploit knowledge of the underlying graph for graph signal recovery from samples obtained over a subset of graph nodes or from noisy observations of all nodes [15, 16]. Through these works, the idea of bandlimitidness was extended to graph signals and the concept of smooth graph signals was introduced. Since the underlying graph is not always available beforehand, the third wave of GSP analyzed the problem of recovering the underlying graph through observations over the graph [6, 8, 7, 9, 10, 11, 12]. Finally, the fourth wave of research in GSP has focused on joint identification of the underlying graph and graph signals from samples/noisy observations using interrelated properties of graph signals and graphs [17, 18, 19, 12].

Within the third set of papers in GSP, our work falls in the category of combinatorial graph Laplacian estimation [9, 10, 11, 12] from graph signals. Combinatorial graph Laplacian refers to the unnormalized Laplacian of an unstructured graph with no self loops [11]. The earliest of works in this category [12] aims to jointly denoise noisy graph signals observed at the nodes of the graph and also learn the graph from these denoised graph signals. The authors pose this problem as a multiconvex problem in the graph Laplacian and the denoised graph signals, and then solve it via an alternating minimization approach that solves a convex problem in each unknown. The authors in [8] examine the problem of graph Laplacian learning when the eigenvectors of the graph Laplacian are known beforehand. They achieve this by formulating a convex program to learn a valid graph Laplacian (from a feasible set) that is diagonalized by the noiseless and noisy Laplacian eigenvectors.

The work in [9] takes a slightly different route and learns a sparse unweighted graph Laplacian matrix from noisy graph signals through an alternating minimization approach that restricts the number of edges in the graph. In contrast to earlier work, [6] focuses on learning a graph diffusion process from observations of stationary signals on graphs through convex formulations. In this regard, the authors also consider different criteria in addition to searching for a valid Laplacian and devise specialized algorithms to infer the diffusion processes under these criteria. In [13, 10] the graph learning problem is adressed by posing it in terms of learning a sparse weighted adjacency matrix from the observed graph signals. Finally, the authors in [11] provide a comprehensive unifying framework for inferring several types of graph Laplacians from graph signals. They also make connection with the state-of-the-art and describe where the past works fit in light of their proposed framework.

It should be mentioned here that since Laplacian matrices (and thus adjacency matrices) are related to precision matrices, defined as the (pseudo-)inverses of covariance matrices, imposing a structure on the graph adjacency matrix amounts to imposing a structure on the covariance of data. Earlier works in the field have already made comparisons of Laplacian learning approaches with those for learning precision matrices from data [12, 11], and established the superior graph recovery performance of Laplacian-based learning. Some recent works have also investigated learning structured covariance and precision matrices, and their usefulness in efficiently representing real-world datasets [20, 21, 22]. While these models work well in practice, we will demonstrate through our experiments that there are scenarios where graph-based learning outperforms structured covariance-based learning (see Sec. VI for details).

Our first contribution in this work is a novel formulation of the graph learning problem as a linear program. We show, both theoretically and empirically, that graph adjacency matrices can be learned through a simple and fast linear program. We then shift our attention towards learning structured graphs. Most of the prior works regarding graph learning have considered arbitrary graphs with either some connectivity constraints [10], or no constraints at all [7, 12]. In all cases, the complexity of the graph learning procedure and the number of free parameters scale quadratically with the number of nodes in the graph, which can be prohibitively large in real-world scenarios. In contrast, our work focuses on inferring the underlying graph from graph signals in the context of structured graphs. Specifically, we investigate graphs that can be represented as Kronecker, Cartesian, and strong products of several smaller graphs. We first show how, for these product graphs, the graph adjacency matrix, the graph Laplacian, the graph Fourier transform, and the graph smoothness measure can be represented with far fewer parameters than required for arbitrary graphs. This reduction in number of parameters to be learned results in reduced sample complexity and helps avoid overfitting. Afterwards, we outline an algorithm to learn these product graphs from the data and provide convergence guarantees for the proposed algorithm in terms of the estimation error of factor graphs. We validate the performance of our algorithm with numerical experiments on both synthetic and real data.

The following notational convention is used throughout the rest of this paper. We use bold lower-case and bold-upper case letters to represent vectors and matrices, respectively. Calligraphic letters are used to represent tensors, which are arrays of three or more dimensions. For a tensor $\mathcal{T}$, $\mathcal{T}_{(i)}$ represents its matricization (flattening) in the $i$-th mode and $\mathrm{vec}(\mathcal{T})$ represents its vectorization along the first mode [23]. Also, “$:$” represents the scalar product or double dot product between two tensors, which results in a scalar [23]. The Hadamard product (elemntwise product) of two vectors or matrices is denoted by “$\circ$”. For a matrix $\mathbf{A}$, $\|\mathbf{A}\|_{F}$ represents its Frobenius norm, $\|\mathbf{A}\|$ represents its spectral norm, and $\mathbf{A}^{\dagger}$ represents its Moore-Penrose inverse. Moreover, $\|\mathbf{A}\|_{1}$ represents the $\ell_{1}$-norm of the entries of $\mathbf{A}$, while $\|\mathbf{A}\|_{1,\text{off}}$ represents the $\ell_{1}$-norm of the off-diagonal entries of $\mathbf{A}$. The Kronecker and Cartesian products of two matrices $\mathbf{A}$ and $\mathbf{B}$ are denoted by $\mathbf{A}\otimes\mathbf{B}$ and $\mathbf{A}\oplus\mathbf{B}$, respectively [24]. The strong product of two matrices is denoted by $\mathbf{A}\boxtimes\mathbf{B}=\mathbf{A}\otimes\mathbf{B}~{}+~{}\mathbf{A}\oplus\mathbf{B}$, which is the sum of Cartesian and Kronecker products of the respective matrices. Furthermore, $\underset{\mathbf{s}}{\otimes}$, $\underset{\mathbf{s}}{\oplus}$, and $\underset{\mathbf{s}}{\boxtimes}$, respectively denote the Kronecker, Cartesian, and strong products taken over the indices provided by the entries of the vector $\mathbf{s}$. Finally, $\times_{i}$ represents matrix multiplication in the $i$-th mode of a tensor and $\underset{\mathbf{s}}{\times}$ represents matrix multiplications in the modes of a tensor specified in the entries of the vector $\mathbf{s}$.

We use $\mathrm{diag}(\mathbf{x})$ to denote a diagonal matrix with diagonal entries given by the entries of the vector $\mathbf{x}$, $\mathbf{1}$ to denote a vector of all ones with appropriate length, and $\mathds{1}$ to denote a tensor of all ones of appropriate size. We denote the set with elements $\{1,2,\dots,K\}$ as $[K]$, and $[K]\setminus k$ represents the same set without the element $k$. The sets of valid Laplacian and weighted adjacency matrices are represented by $\mathcal{L}$ and $\mathcal{W}$, respectively. The set of weighted adjacency matrices with any product structure is denoted by $\mathcal{W}_{p}$.

The rest of this paper is organized as follows. In Sec. II we give a probabilistic formulation of the graph learning problem in line with existing literature. Then, we propose our novel formulation of the graph learning problem as a liner program in Sec. III. Sec. IV describes the motivation for product graphs and formulates the graph learning problem in the context of product graphs. In Sec. V we propose an algorithm for learning product graphs from data and derive error bounds on estimated factor graphs. We present our numerical experiments with synthetic and real datasets in Sec. VI, and the paper is concluded in Sec. VII.

We now present an algorithm, named Graph learning with Linear Programming (GLP), for solving the linear graph learning problem (III). To proceed, note that the objective term in the graph learning problem can be reformulated as:

where the matrix $\widetilde{\mathbf{S}}=(\underset{m=1}{\overset{M_{0}}{\sum}}\mathbf{1}\bar{\mathbf{x}}_{m}^{T}-\underset{m=1}{\overset{M_{0}}{\sum}}\mathbf{x}_{m}\mathbf{x}_{m}^{T})/M_{0}$, the vector $\bar{\mathbf{x}}_{m}^{T}=\mathbf{x}_{m}^{T}\mathrm{diag}(\mathbf{x}_{m})=\mathbf{x}_{m}^{T}\circ\mathbf{x}_{m}^{T}$, $\mathrm{vec}(\mathbf{W})=\mathbf{M}\mathbf{w}$, $\mathbf{w}$ is a vector of distinct elements of upper triangular part of the symmetric matrix $\mathbf{W}$, and $\mathbf{M}$ is the duplication matrix that duplicates the elements from $\mathbf{w}$ to generate a vectorized version of $\mathbf{W}$. With this rearrangement of the objective and $\mathbf{W}$, our graph learning problem can be posed as:

where $\mathbf{A}$ is a matrix that represents the equality constraints from (III) in terms of equality constraints on $\mathbf{w}$, $\mathbf{b}=[\mathbf{0}^{T},n]^{T}$, and $F$ is the set containing the indices of the off-diagonal elements in $\mathbf{w}$. Once the solution $\widehat{\mathbf{w}}$ is obtained, it can be converted to the symmetric adjacency matrix $\widehat{\mathbf{W}}$, which can then be used to get $\widehat{\mathbf{L}}$.

The standard way of solving a linear program with mixed (equality and inequality) constraints is through interior point methods whose complexity scales quadratically with the problem dimension [25]. A better alternative is to deploy a first-order method whose per-iteration complexity is linear in the number of nonzero entries of $\mathbf{A}$. However, a first-order method would exhibit slow convergence for a linear program because of the lack of smoothness and strong convexity in linear programs [26]. To overcome these issues, linear programs have been solved through the Alternating Direction Method of Multiplier (ADMM) [27, 26]. To solve our proposed linear formulation of graph learning, we follow a recent algorithm proposed in [26]. This ADMM-based algorithm for linear programs proposed a new variable splitting scheme that achieves a convergence rate of $\mathcal{O}(\|\mathbf{A}\|^{2}\log(1/\epsilon))$, where $\epsilon$ is the desired accuracy of the solution. To this end, we start by modifying the original graph learning problem with the introduction of an additional variable $\mathbf{y}$ as follows:

where $\mathbf{c}=\alpha\mathbf{M}^{T}\widetilde{\mathbf{s}}$. The corresponding augmented Lagrangian can then be expressed as follows:

where $h(\mathbf{y})$ denotes the non-negativity constraint on the entries of $\mathbf{y}$ indexed by $F$, i.e., $\forall i\in F$, $h(\mathbf{y})=0$ when $y_{i}\geq 0$, and $h(\mathbf{y})=\infty$ when $y_{i}<0$. Moreover, $\mathbf{z}=[\mathbf{z}_{\mathbf{w}}^{T},\mathbf{z}_{\mathbf{y}}^{T}]^{T}$, $\mathbf{z}_{\mathbf{w}}$ and $\mathbf{z}_{\mathbf{y}}$ are the Lagrange multipliers, $\mathbf{A}_{\mathbf{w}}=[\mathbf{A}^{T},\mathbf{I}]^{T}$, $\mathbf{A}_{\mathbf{y}}=[\mathbf{0},-\mathbf{I}]^{T}$, and finally $\widetilde{\mathbf{b}}=[\mathbf{b}^{T},\mathbf{0}^{T}]^{T}$. One can then use ADMM to go through the steps outlined in Algorithm 1 until convergence to obtain $\widehat{\mathbf{w}}$.

In Algorithm 1, $[\cdot]_{\underset{F}{\geqslant}\mathbf{0}}$ is entrywise thresholding that projects the entries with indices in $F$ to the nonnegative orthant. As we can see from the algorithm, all updates have closed-form solutions and the most computationally expensive step is the $\mathbf{w}^{(t+1)}$ update that involves matrix inversion. This matrix inversion, however, can be computed efficiently using the identity $(\mathbf{I}+\mathbf{A}^{T}\mathbf{A})^{-1}=\mathbf{I}-\mathbf{A}^{T}(\mathbf{I}+\mathbf{A}\mathbf{A}^{T})^{-1}\mathbf{A}$. Since matrix $\mathbf{A}$ is a fat matrix, $\mathbf{A}\mathbf{A}^{T}$ has smaller dimensions than $\mathbf{A}^{T}\mathbf{A}$. Moreover, one can easily see that $\mathbf{A}\mathbf{A}^{T}$ is a matrix of dimensions $(n+1)\times(n+1)$, and

where $c_{n}=2n^{2}-n$. In addition, the inverse only needs to be computed once at the start of the algorithm since this matrix is deterministic and depends only on the size of the adjacency matrix being estimated.

The number of parameters that one needs to learn a graph adjacency matrix is $\frac{n(n+1)}{2}$. This implies that the number of unknown parameters scales quadratically with the number of nodes in the graph. Additionally, the per-iteration computational complexity of the proposed method also scales quadratically with the number of nodes [26]. The same computational and memory complexities also hold for the existing state-of-the-art graph learning algorithms [12, 10, 13, 9, 11]. However, while these complexities are manageable for small graphs, for real-world datasets with even hundreds of nodes current methods become prohibitive. To overcome these issues, we next examine the problem of learning product graphs.

For the Kronecker product of graphs $G_{k}$, for $k=1,\dots,K_{0}$, with adjacency matrices $\mathbf{W}_{k}$, the Kronecker product graph can be expressed as $G=\underset{[K_{0}]}{\otimes}G_{k}=G_{K_{0}}\otimes G_{K_{0}-1}\otimes\dots\otimes G_{1}$. The respective Kronecker-structured adjaceny matrix of the resultant graph can written in terms of component/factor adjacency matices as $\mathbf{W}=\underset{[K_{0}]}{\otimes}\mathbf{W}_{k}$. Additionally, if the factor adjacency matrix $\mathbf{W}_{k}$ can be expressed via eigenvalue decomposition (EVD) as $\mathbf{W}_{k}=\mathbf{U}_{k}\boldsymbol{\Lambda}_{k}\mathbf{U}_{k}^{T}$, then the Kronecker adjaceny matrix can be written as (using properties in [3]):

One can see that both the eigenmatrix and the eigenvalue matrix of the Kronecker adjacency matrix have a Kronecker structure in terms of the component eigenmatrices and component eigenvalue matrices, respectively. Given the number of edges in the component graphs are $|E_{k}|$, the number of edges in the Kronecker graph are $|E|=2^{K_{0}-1}\overset{K_{0}}{\underset{k=1}{\prod}}|E_{i}|$.

An example of Kronecker product graph is the bipartite graph of a recommendation system like Netflix [28] where the graph between users and movies can be seen as a Kronecker product of two smaller factor graphs. In fact, the adjacency matrix of any bipartite graph can be represented in terms of a Kronecker product of appropriate factor matrices [29]. As adjacency matrices are also closely related to precision matrices (inverse covariance matrices), and inverse of Kronecker product is Kronecker product of inverses [3], imposing Kronecker structure on the adjacency matrix also amounts to imposing a Kronecker structure on the covariance matrix of the data.

The optimization problem in (III) can be specialized to the case of learning Kronecker graphs by explicitly imposing the Kronecker product structure on the adjacency matrix and posing the problems in terms of the individual factor adjacency matrices, rather than the bigger adjacency matrix produced after the product. This leads us to the following nonconvex problem for learning Kronecker graphs:

The Cartesian product (also called Kronecker sum product) of graphs $G_{k}$ is represented as $G=\underset{[K_{0}]}{\oplus}G_{k}=G_{K_{0}}\oplus G_{K_{0}-1}\oplus\dots\oplus G_{1}$. The correspoding cartesian adjacency matrix can be written in terms of the component adjacency matrices as $\mathbf{W}=\underset{[K_{0}]}{\oplus}\mathbf{W}_{k}$. Furthermore, with the EVD of the component adjacency matrices, the Cartesian adjacency matrix can be decomposed as [3]:

This means that the eigenmatrix and the eigenvalue matrix of the Cartesian adjacency matrix are represented, respectively, as Kronecker and Cartesian products of component eigenmatrices and eigenvalue matrices. The number of edges in the Cartesian graph can be found as $|E|=\sum_{k=1}^{K_{0}}(\underset{[K_{0}]\setminus k}{\otimes}n_{i})|E_{k}|$, where $|E_{k}|$ represents the number of edges and $n_{k}$ represents the number of vertices in the $k$-th component graph.

A typical exmaple of a Cartesian product graph is images. Images reside on two-dimensional rectangular grids that can be represented as the Cartesian product between two line graphs pertaining to the rows and columns of the image [3]. A social network can also be approximated as a Cartesian product of an inter-community graph with an intra-community graph [3].

Similar to the previous discussion, the optimization problem in (III) can be specialized to learning Cartesian graphs by explicitly imposing the Cartesian structure and posing the problem in terms of the factor adjacency matrices as follows:

The strong product of graphs $G_{k}$ can be represented as $G=\underset{[K_{0}]}{\boxtimes}G_{k}=G_{K_{0}}\boxtimes G_{K_{0}-1}\boxtimes\dots\boxtimes G_{1}$. The respective strong adjacency matrix of the resultant strong graph is given in terms of the component adjacency matrices as $\mathbf{W}=\underset{[K_{0}]}{\boxtimes}\mathbf{W}_{k}$, and can be further expressed as:

in terms of EVD of the component adjacency matrices.

The strong product graphs can be seen as the sum of Kronecker and Cartesian products of the factor adjacency matrices. An example of data conforming to the strong product graph is a spatiotemporal sensor network graph, which consists of a strong product of a spatial graph and a temporal graph (representing the temporal dependencies of the sensors). The spatial graph has as many nodes as the number of sensors in the sensor network and represents the spatial distribution of sensors. On the other hand, the temporal graph has as many nodes as the number of temporal observations of the whole sensor network and represents the overall temporal dynamics (changes in connectivity over time) of the network [3].

By making the strong product structure explicit in terms of the factor adjacency matrices, the optimization problem for learning strong graphs can be expressed as the following nonconvex problem:

One can see from (IV-A), (IV-B), and (IV-C) that the graph Fourier transform of a product graph (which is the eigenmatrix of the product adjacency matrix), has a Kronecker structure in terms of the eigenmatrices of the component graph adjacency matrices: $\mathbf{U}=\underset{[K_{0}]}{\otimes}\mathbf{U}_{k}$. In terms of the implementation of the graph Fourier transform, this structure provides an efficient implementation of the graph Fourier as (using the properties of Kronecker product and tensors [23]):

where $\mathbf{x}\in\mathbb{R}^{n_{1}n_{2}\dots n_{K_{0}}}$ is an arbitrary graph signal on the product graph, and $\mathcal{X}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{K_{0}}}$ represents appropriately tensorized version of the signal $\mathbf{x}$. Because of this, one does not need to form the huge Fourier matrix $\mathbf{U}$ and can avoid costly matrix multiplications by just applying the component graph Fourier matrices to each respective mode of the tensorized observation $\mathcal{X}$ and then vectorizing the result.

Smoothness of a graph signal is one of the core concepts in graph signal processing [1, 2, 3, 4, 5] and product graph Laplacians provide an efficient representation for the notion of smoothness. The smoothness of a graph signal can be measured through the Dirchlet energy defined as $\mathbf{x}^{T}\mathbf{L}\mathbf{x}$. The Dirichlet energy can be reexpressed as: $\mathbf{x}^{T}\mathbf{L}\mathbf{x}=\mathbf{x}^{T}(\mathbf{D}-\mathbf{W})\mathbf{x}=\mathbf{x}^{T}\mathbf{D}\mathbf{x}-\mathbf{x}^{T}\mathbf{W}\mathbf{x}.$ Let us now focus on each term separately in the context of product graphs. For the term involving $\mathbf{W}$ we have:

Similarly, the term involving $\mathbf{D}$ can be re-expressed as:

which can be computed efficiently along the lines of (IV-E). With this reformulation, one circumvents the need to explicitly form the prohibitively large eigenmatrix $\mathbf{U}$ and can evaluate the Dirichlet energy much more efficiently with just mode-wise products with the smaller component eigenmatrices.

Let us consider an unknown graph $G$ with the number of nodes $|V|=n=\prod_{k=1}^{K_{0}}n_{k}$, where $n_{k}$ represents the number of nodes in each component graph and $K_{0}$ is total number of component graphs. If one were to learn this graph by means of an arbitrary adjacency matrix, the number of parameters that needs to be estimated would be $\frac{n(n+1)}{2}$ (since the graph adjacency matrix is a symmetric matrix). On the other hand, for the same graph, by utilizing the product model of the graph adjacency matrix, one would need to estimate only $\sum_{k=1}^{K_{0}}\frac{n_{k}(n_{k}+1)}{2}$ parameters. Considering the special case of $n_{1}=n_{2}=\cdots=n_{K_{0}}=\bar{n}$, and imposing the product structure on graph adjacency matrix reduces the number of parameters needed to be learned by $\bar{n}^{K_{0}-1}/K_{0}$!

Since Algorithm 2 utilizes factor-wise minimization, we can characterize the error for product graph learning in terms of the factor-wise errors of each factor adjacency matrix. The error bounds for learning Kronecker graphs are provided in the following theorem with the proof in Appendix A.

The following theorem characterizes the error of the factor-wise minimization of the Cartesian graph learning problem.

The proof of this theorem is given in Appendix B. The theorem states, remarkably, that the objective function for learning Cartesian product graphs is convex and separable in each factor adjacency matrix, i.e., the objective function can be represented as a sum of linear terms, each of which is dependent on only one factor adjacency matrix. As a consequence of this fact, for learning Cartesian graphs, one can optimize over all factor adjacency matrices in parallel!

The following theorem, with its proof in Appendix C, characterizes the behavior of factor-wise minimization for learning strong product graphs.

As a byproduct of Theorem 1, we can also obtain an error bound for arbitrary graph learning problem. To this end, see that the objective in (III) for learning arbitrary graphs can be expressed as $\alpha/M_{0}\sum_{m=1}^{M_{0}}\mathbf{x}_{m}^{T}\mathbf{L}\mathbf{x}_{m}=\alpha\mathrm{trace}(\mathbf{D}\bar{\mathbf{S}})-\alpha\mathrm{trace}(\mathbf{W}\mathbf{S})$, where $\mathbf{S}$ and $\bar{\mathbf{S}}$ are similiar to the definitions in Appendix A. Following along the lines of Theorem 1, by solving (III) via Algorithm 1 one is guaranteed to converge to the original generating adjacency matrix of the unstructured graph with the error $\|\widehat{\mathbf{W}}-\mathbf{W}^{*}\|_{F}=\mathcal{O}_{P}\Big{(}\sqrt{\log(n)/M_{0}}\Big{)}$, with a high probability as $M\rightarrow\infty$.

The computational complexity of solving each factor-wise problem scales quadratically with the number of nodes in the graph. This implies that when the product structure is imposed, one only has to solve $K_{0}$ smaller problems each with computational complexity of $\mathrm{O}(\bar{n}^{2})$, assuming the special case of $n_{1}=n_{2}=\dots=n_{K_{0}}=\bar{n}$. In contrast, for learning unstructured graphs the computational complexity would scale as $\mathrm{O}(n^{2})=\mathrm{O}(\prod_{k=1}^{K_{0}}n_{k}^{2})\sim\mathrm{O}(\bar{n}^{2K_{0}})$. Thus, the computational gains are huge in comparison to the original problem for learning unstructured graphs!

This computational gain is only made possible due to the way that we pose the original graph learning problem as a linear program. This way of posing the original problem made the objective amenable to factor-wise minimization, which would not be possible if the original graph learning problem was posed in the traditional way from the existing literature.

The overall convergence of the algorithm to a stationary point of the problem can be established through the following theorem, whose proof is provided in Appendix D.

We now provide an insight into our theorems statements with regards to the required number of observation. Theorem 1, 2, and 3 claim that the estimated factor lies within a ball of radius $\frac{\log(n_{k})}{(n/n_{k})M_{0}}$ around the true factor. The accuracy of the estimate increases with the number of available observations $M_{0}$, and the product of the dimensions of the other factors $n/n_{k}=\prod_{j\neq k}n_{j}$. Moreover, the accuracy decreases with the increasing dimensions of the factor being estimated.

Taking a closer look at the error bounds for learning arbitrary and structured graphs reveals an important point. The denominator in the error bound is the number of observations available to estimate the graph. For $M_{0}$ observed graph signals, the number of observations available to arbitrary graph learning are (obviously) $M_{0}$; however, for estimating the $k$-th factor adjacency when learning product graphs the effective number of observations are $\prod_{j\neq k}n_{j}\times M_{0}$. This means that imposing the product structure results in an increased number of effective observations to estimate each factor adjacency matrix. This combined with the reduced number of parameters required to learn these graphs makes product graphs very attractive for real world applications.

To showcase the performance of our new formulation for graph learning, we run synthetic experiments on a graph with $n=64$ nodes. We generate various different graph types such as: (i) a sparse random graph with Gaussain weights, (ii) an Erdos-Renyi random graph with edge probability 0.7, (iii) a scale-free graph with preferential attachment (6 edges at each step), (iv) a random regular graph where each node is connected to $0.7n$ other nodes, (v) a uniform grid graph, (vi) a spiral graph , (vii) a community graph , and (viii) a low stretch tree graph on a grid of points. The related details of how the graphs are simulated can be found in [30, 12]. For each kind of graph, we generate 20 different realizations, and for each realization we generate observations using a degenerate multivariate Gaussian distribution with the graph Laplacian as the precision matrix [12, 11, 13].

We compare the performance of our proposed method with two other state-of-the-art methods for arbitrary graph learning: (i) combinatorial graph learning from [11] (which we refer to as CGL), and (ii) graph learning method from [10] (which we refer to as LOG), which also aims to learn a combinatorial graph Laplacian through a slightly different optimization problem than [11]. We choose $\alpha$ for our algorithm in the range $0.75^{\{0:14\}}\times\sqrt{\frac{\log(n)}{M_{0}}}$, as dictated by the error bounds for learning graphs in App. A and by the existing literature [12, 11, 13]. Furthermore, we choose $\rho=0.75/\log(M_{0})$ as the value that works best for most cases. For each algorithm, in the prescribed range of the optimization parameters, we choose the parameters that produce the best results.

The results of our experiments are shown as F-measure values in Fig. 1. F-measure is the harmonic mean of precision and recall, and signifies the overall accuracy of the algorithm [11]. Precision here denotes the fraction of true graph edges recovered among all the recovered edges, and recall signifies the fraction of edges recovered from the true graph edges. One can see that our algorithm (except for community graphs) performs just as well or better than the existing state-of-the-art algorithms. Moreover, the average performance over all graphs in Fig. 2 shows that on average we outperform the existing algorithms. A runtime comparison of all algorithms in Fig. 2 also reveals competitive run time for our proposed scheme. The LOG algorithm, with the smallest runtime, has a huge computational overhead incurred by the construction of a matrix of pairwise distances of all rows of data matrix $\mathbf{X}$. In contrast, other algorithms work with the graph signal observations directly and do not require preprocessing steps.

We now present the results of our numerical experiments involving synthetic data for product graphs. We run experiments for random Erdos-Renyi factor graphs with $n=n_{1}n_{2}n_{3}=12\times 12\times 12$ nodes, and having either Cartesian, Kronecker or strong structure. We then use our proposed algorithms to learn the generated graphs with varying number of observations and compare the performance with LOG as its performance was the second best in Fig. 1. The results for all three types of product graphs are shown in Fig. 3 (top). For a fixed number of observations, Cartesian product graphs can be learned with the highest F-measure score followed by strong and then Kronecker graphs. The figure also shows that for each graph, imposing product structure on the learned graph drastically improves the performance of the learning algorithm. Fig. 3 (bottom) also shows the runtimes comparison of our approach BPGL with the algorithm in [10]. Even for a graph of this size, with total number of nodes $n=1728$, we can see a considerable reduction in runtimes. Thus, our learning algorithm that explicitly incorporates the product structure of the graph enjoys superior performance, reduced computational complexity and faster runtimes.

The first real data we use for experimentation is NCEP wind speed data. The NCEP wind speed data [31] represents wind conditions in the lower troposphere and contains daily averages of U (east-west) and V (north-south) wind components over the years $1948-2012$. Similar to the experiments in [22] with preprocessed data, we use a grid of $n_{1}n_{2}=10\times 10=100$ stations to extract the data available for the United States region. From this extracted data, we choose the years $2003-2007$ for training and keep the years $2008-2012$ for testing purposes. Using a non-overlapping window of length $n_{3}=8$, which amounts to dividing the data into chunks of 8 days, we obtain $M_{0}=228$ samples, each of length $n=n_{1}n_{2}n_{3}$. Therefore, for graph learning we have 228 samples, where each sample contains spatiotemporal data for 100 stations over 8 days. Next, the testing procedure consists of introducing missing values in each sample of the test data by omitting the data for the 8th day, and then using a linear minimum-mean-square-error (LMMSE) estimator [32, Chapter 4] to predict the missing values. As for learning, 228 samples are obtained for testing through the same procedure. Our proposed method estimates the (structured) adjacency matrix of the graph (which is related to the precision matrix of the data), and we use $\mathbf{W}+\mathbf{I}$ in place of the data covariance for the LMMSE estimator.

We make a comparison of the following approaches: (1) sample covariance matrix (SCM) learning, (2) the permuted rank-penalized least squares (PRLS) approach [22] with $r=6$ Kronecker components, (3) PRLS with $r=2$ Kronecker components, (4) time-varying graph learning (TVGL) approach from [33] (which was shown to outperform the approach in [34]), (5) spatiotemporal strong graph with BPGL with a spatial component of size $n_{1}n_{2}$ and a temporal component of size $n_{3}$, and (6) spatiotemporal Cartesian graph with BPGL of the same dimensions. The parameters for PRLS were chosen for optimal performance as given in [22]. The optimization parameters for TVGL and BPGL were manually tuned for best performance. It should be noted here that SCM and PRLS aim to learn a covariance matrix and a structured covariance matrix from the data, respectively.

The SCM aims to estimate $\frac{n(n+1)}{2}$ parameters, while the number of parameters that PRLS needs to estimate is $r(\frac{n_{1}n_{2}(n_{1}n_{2}+1)}{2}+\frac{n_{3}(n_{3}+1)}{2})$. On the other hand, TVGL estimates $n_{3}(\frac{n_{1}n_{2}(n_{1}n_{2}+1)}{2})$ parameters, while BPGL needs to learn only $\frac{n_{1}n_{2}(n_{1}n_{2}+1)}{2}+\frac{n_{3}(n_{3}+1)}{2}$parameters for both strong and Cartesian graphs. The mean prediction root-mean-squared errors (RMSE) for all methods are shown in Table I. One can see that our proposed method outperforms PRLS and TVGL while estimating far fewer parameters than both. The table also shows that learning a strong graph for this data results in a higher RMSE reduction over the baseline (SCM), and is thus better suited for this data than the Cartesian product graph.

The second real data that we use as an application for our proposed algorithm is a part of the ABIDE fMRI dataset [35, 36]. Our aim is to learn the graphs over the fMRI data of control and autistic subjects and to use the learned graphs to higlight the differences in the control and autistic brains. The data we obtain is already preprocessed to remove various fMRI artifacts and for controlization of the obtained scans [37]. The final preprocessed data consists of measurements from $n_{1}=111$ brain regions scanned over $116$ time instances for each subject. The data contains scans for control and autistic subjects. To avoid class imbalance we randomly choose $47$ subjects for each class. Out of the $47$ subjects for each class, we then randomly choose $30$ subjects for training and keep the remaining $17$ for testing purposes. We use a non-overlapping window length of $n_{2}=29$ which results into $M_{0}=120$ samples of length $n=n_{1}\times n_{2}$.

As before, we compare the performance of our proposed approach with SCM and PRLS. Table II shows the results of our experiments. One can see that our approach performs very similar to PRLS for both Cartesian and strong product graphs, all the while using much fewer parameters (five times fewer). We also see that strong product graphs are more suited to model brain activity. The work in [37] suggests that autistic brains exhibit hypoconnectivity in different regions of the brain as compared to control subjects. The results from our graph learning procedure go a step further and bring more insight into the spatiotemporal dynamics of the brain. Firstly, as already suggested in [37], we see clear evidence of spatial hypoconnectivity (see Fig. 4). More importantly, our learned graphs in Fig. 5 reveal that, in addition to spatial hypoconnectivity, autistic brains also suffer from temporal hypoconnectivity.

The final dataset that we experiment on is the estrogen receptor data [38, 39], which consists of 157 samples of 693 probe sets related to estrogen receptor pathway. We aim to learn a Kronecker structured graph on this data using $120$ randomly selected samples for training and the remaining $37$ for testing. We choose Kronecker structured graph for this data because transposable models, i.e., models that learn a Kronecker structured data covariance, have been shown to work well for genomic data in the existing literature [28]. And as pointed out in Sec. IV-A, Kronecker structured adjacency matrix corresponds to a Kronecker structured data covariance. For testing purposes, we follow a procedure similar to the previous subsections.

We compare our graph learning approach with SCM, PRLS and sparse covariance estimation (SEC) from [40]. Optimization parameters are manually tuned for best results for each method. We learn a graph through our method (and covariance through PRLS), with a Kronecker structure composed of two factor matrices of dimensions $n_{1}=21$ and $n_{2}=33$. We then use LMMSE estimator to predict $33$ probe set measurements removed from the test data. PRLS, SEC and BPGL result in an improvement of $\mathbf{0.91347}$ dB, $\mathbf{0.93598}$ dB, and $\mathbf{1.0242}$ dB over SCM, respectively. This demonstrates that our method outperforms the state-of-the-art unstructured and structured sparse covariance estimation techniques, and provides a better model for real datasets.