# SG-PALM: a Fast Physically Interpretable Tensor Graphical Model

Yu Wang  Alfred Hero

###### Abstract

We propose a new graphical model inference procedure, called SG-PALM, for
learning conditional dependency structure of high-dimensional tensor-variate
data. Unlike most other tensor graphical models the proposed model is
interpretable and computationally scalable to high dimension. Physical
interpretability follows from the Sylvester generative (SG) model on which SG-
PALM is based: the model is exact for any observation process that is a
solution of a partial differential equation of Poisson type. Scalability
follows from the fast proximal alternating linearized minimization (PALM)
procedure that SG-PALM uses during training. We establish that SG-PALM
converges linearly (i.e., geometric convergence rate) to a global optimum of
its objective function. We demonstrate the scalability and accuracy of SG-PALM
for an important but challenging climate prediction problem: spatio-temporal
forecasting of solar flares from multimodal imaging data.

Machine Learning, ICML

##  1 Introduction

High-dimensional tensor-variate data arise in computer vision (video data
containing multiple frames of color images), neuroscience (EEG measurements
taken from different sensors over time under various experimental conditions),
and recommending system (user preferences over time). Due to the non-
homogeneous nature of these data, second-order information that encodes
(conditional) dependency structure within the data is of interest. Assuming
the data are drawn from a tensor normal distribution, a straightforward way to
estimate this structure is to vectorize the tensor and estimate the underlying
Gaussian graphical model associated with the vector. However, such an approach
ignores the tensor structure and requires estimating a rather high dimensional
precision matrix, often with insufficient sample size. For instance, in the
aforementioned EEG application the sample size is one if we aim to estimate
the dependency structure across different sensors, time and experimental
conditions for a single subject. To address such sample complexity challenges,
sparsity is often imposed on the covariance $\mathbf{\Sigma}$ or the inverse
covariance $\mathbf{\Omega}$, e.g., by using a sparse Kronecker product (KP)
or Kronecker sum (KS) decomposition of $\mathbf{\Sigma}$ or $\mathbf{\Omega}$.
The earliest and most popular form of sparse structured precision matrix
estimation approaches represent $\mathbf{\Omega}$, equivalently
$\mathbf{\Sigma}$, as the KP of smaller precision/covariance matrices (Allen &
Tibshirani, 2010; Leng & Tang, 2012; Yin & Li, 2012; Tsiligkaridis et al.,
2013; Zhou, 2014; Lyu et al., 2019). The KP structure induces a generative
representation for the tensor-variate data via a separable covariance/inverse
covariance model. Alternatively, Kalaitzis et al. (2013); Greenewald et al.
(2019) proposed to model inverse covariance matrices using a KS
representation. Rudelson & Zhou (2017); Park et al. (2017) proposed KS-
structured covariance model which corresponds to an errors-in-variables model.
The KS (inverse) covariance structure corresponds to the Cartesian product of
graphs (Kalaitzis et al., 2013; Greenewald et al., 2019), which leads to more
parsimonious representations of (conditional) dependency than the KP. However,
unlike the KP model, KS lacks an interpretable generative representation for
the data. Recently, Wang et al. (2020) proposed a new class of structured
graphical models, called the Sylvester graphical models, for tensor-variate
data. The resulting inverse covariance matrix has the KS structure in its
square-root factors. This square-root KS structure is hinted in the paper to
have a connection with certain physical processes, but no illustration is
provided.

A common challenge for structured tensor graphical models is the efficient
estimation of the underlying (conditional) dependency structures. KP-
structured models are generally estimated via extension of GLasso (Friedman et
al., 2008) that iteratively minimize the $\ell_{1}$-penalized negative
likelihood function for the matrix-normal data with KP covariance. This
procedure was shown to converge to some local optimum of the penalized
likelihood function (Yin & Li, 2012; Tsiligkaridis et al., 2013). Similarly,
Kalaitzis et al. (2013) further extended GLasso to the KS-structured case for
$2$-way tensor data. Greenewald et al. (2019) extended this to multiway
tensors, exploiting the linearity of the space of KS-structured matrices and
developing a projected proximal gradient algorithm for KS-structured inverse
covariance matrix estimation, which achieves linear convergence (i.e.,
geometric convergence rate) to the global optimum. In Wang et al. (2020), the
Sylvester-structured graphical model is estimated via a nodewise regression
approach inspired by algorithms for estimating a class of vector-variate
graphical models (Meinshausen et al., 2006; Khare et al., 2015). However, no
theoretical convergence result for the algorithm was established nor did they
study the computational efficiency of the algorithm.

In the modern era of big data, both computational and statistical learning
accuracy are required of algorithms. Furthermore, when the objective is to
learn representations for physical processes, interpretablility is crucial. In
this paper, we bridge this “Statistical-to-Computational-to-Interpretable gap”
for Sylvester graphical models. We develop a simple yet powerful first-order
optimization method, based on the Proximal Alternating Linearized Minimization
(PALM) algorithm, for recovering the conditional dependency structure of such
models. Moreover, we provide the link between the Sylvester graphical models
and physical processes obeying differential equations and illustrate the link
with a real-data example. The following are our principal contributions:

  1. 1.

A fast algorithm that efficiently recovers the generating factors of a
representation for high-dimensional multiway data, significantly improving on
Wang et al. (2020).

  2. 2.

A comprehensive convergence analysis showing linear convergence of the
objective function to its global optimum and providing insights for choices of
hyperparameters.

  3. 3.

A novel application of the algorithm to an important multi-modal solar flare
prediction problem from solar magnetic field sequences. For such problems, SG-
PALM is physically interpretable in terms of the partial differential
equations governing solar activities proposed by heliophysicists.

##  2 Background and Notation

###  2.1 Notations

In this paper, scalar, vector and matrix quantities are denoted by lowercase
letters, boldface lowercase letters and boldface capital letters,
respectively. For a matrix
$\mathbf{A}=(\mathbf{A}_{i,j})\in\mathbb{R}^{d\times d}$, we denote
$\|\mathbf{A}\|_{2},\|\mathbf{A}\|_{F}$ as its spectral and Frobenius norm,
respectively. We define $\|\mathbf{A}\|_{1,\text{off}}:=\sum_{i\neq
j}|\mathbf{A}_{i,j}|$ as its off-diagonal $\ell_{1}$ norm. For tensor algebra,
we adopt the notations used by Kolda & Bader (2009). A $K$-th order tensor is
denoted by boldface Euler script letters, e.g,
$\bm{\mathcal{X}}\in\mathbb{R}^{d_{1}\times\dots\times d_{K}}$. The
$(i_{1},\dots,i_{K})$-th element of $\bm{\mathcal{X}}$ is denoted by
$\bm{\mathcal{X}}_{i_{1},\dots,i_{K}}$, and the vectorization of
$\bm{\mathcal{X}}$ is the $d$-dimensional vector
$\operatorname{vec}(\bm{\mathcal{X}}):=(\bm{\mathcal{X}}_{1,1,\dots,1},\bm{\mathcal{X}}_{2,1,\dots,1},\dots,\bm{\mathcal{X}}_{d_{1},1,\dots,1},\dots,\bm{\mathcal{X}}_{d_{1},d_{2},\dots,d_{k}})^{T}$
with $d=\prod_{k=1}^{K}d_{k}$. A fiber is the higher order analogue of the row
and column of matrices. It is obtained by fixing all but one of the indices of
the tensor. Matricization, also known as unfolding, is the process of
transforming a tensor into a matrix. The mode-$k$ matricization of a tensor
$\bm{\mathcal{X}}$, denoted by $\bm{\mathcal{X}}_{(k)}$, arranges the mode-$k$
fibers to be the columns of the resulting matrix. The $k$-mode product of a
tensor $\bm{\mathcal{X}}\in\mathbb{R}^{d_{1}\times\dots\times d_{K}}$ and a
matrix $\mathbf{A}\in\mathbb{R}^{J\times d_{k}}$, denoted as
$\bm{\mathcal{X}}\times_{k}\mathbf{A}$, is of size $d_{1}\times\dots\times
d_{k-1}\times J\times d_{k+1}\times\dots d_{K}$. Its entry is defined as
$(\bm{\mathcal{X}}\times_{k}\mathbf{A})_{i_{1},\dots,i_{k-1},j,i_{k+1},\dots,i_{K}}:=\sum_{i_{k}=1}^{d_{k}}\bm{\mathcal{X}}_{i_{1},\dots,i_{K}}A_{j,i_{k}}$.
For a list of matrices $\\{\mathbf{A}_{k}\\}_{k=1}^{K}$ with
$\mathbf{A}_{k}\in\mathbb{R}^{d_{k}\times d_{k}}$, we define
$\bm{\mathcal{X}}\times\\{\mathbf{A}_{1},\dots,\mathbf{A}_{K}\\}:=\bm{\mathcal{X}}\times_{1}\mathbf{A}_{1}\times_{2}\dots\times_{K}\mathbf{A}_{K}$.
Lastly, we define the $K$-way Kronecker product as
$\bigotimes_{k=1}^{K}\mathbf{A}_{k}=\mathbf{A}_{1}\otimes\cdots\otimes\mathbf{A}_{K}$,
and the equivalent notation for the Kronecker sum as
$\bigoplus_{k=1}^{K}\mathbf{A}_{k}=\mathbf{A}_{1}\oplus\dots\oplus\mathbf{A}_{K}=\sum_{k=1}^{K}\mathbf{I}_{[d_{k+1:K}]}\otimes\mathbf{A}_{k}\otimes\mathbf{I}_{[d_{1:k-1}]}$,
where
$\mathbf{I}_{[d_{k:\ell}]}=\mathbf{I}_{d_{k}}\otimes\dots\otimes\mathbf{I}_{d_{\ell}}$.
For the case of $K=2$,
$\mathbf{A}_{1}\oplus\mathbf{A}_{2}=\mathbf{I}_{d_{2}}\otimes\mathbf{A}_{1}+\mathbf{A}_{2}\otimes\mathbf{I}_{d_{1}}$.

###  2.2 Tensor Graphical Models

A random tensor $\bm{\mathcal{X}}\in\mathbb{R}^{d_{1}\times\dots\times d_{K}}$
follows the tensor normal distribution with zero mean when
$\operatorname{vec}(\bm{\mathcal{X}})$ follows a normal distribution with mean
$\mathbf{0}\in\mathbb{R}^{d}$ and precision matrix
$\mathbf{\Omega}:=\mathbf{\Omega}(\mathbf{\Psi}_{1},\dots,\mathbf{\Psi}_{K})$,
where $d=\prod_{k=1}^{K}d_{k}$. Here,
$\mathbf{\Omega}(\mathbf{\Psi}_{1},\dots,\mathbf{\Psi}_{K})$ is parameterized
by $\mathbf{\Psi}_{k}\in\mathbb{R}^{d_{k}\times d_{k}}$ via either Kronecker
product, Kronecker sum, or the Sylvester structure, and the corresponding
negative log-likelihood function (assuming $N$ independent observations
$\bm{\mathcal{X}}^{i},i=1,\dots,N$)

$-\frac{N}{2}\log|\mathbf{\Omega}|+\frac{N}{2}\operatorname*{tr}(\mathbf{S}\mathbf{\Omega}),$
(1)

where $\mathbf{\Omega}=\bigotimes_{k=1}^{K}\mathbf{\Psi}_{k}$,
$\bigoplus_{k=1}^{K}\mathbf{\Psi}_{k}$, or
$\Big{(}\bigoplus_{k=1}^{K}\mathbf{\Psi}_{k}\Big{)}^{2}$ for KP, KS, and
Sylvester models, respectively; and
$\mathbf{S}=\frac{1}{N}\sum_{i=1}^{N}\operatorname{vec}(\bm{\mathcal{X}}^{i})\operatorname{vec}(\bm{\mathcal{X}}^{i})^{T}$.
To encourage sparsity, penalized negative log-likelihood function is proposed

$-\frac{N}{2}\log|\mathbf{\Omega}|+\frac{N}{2}\operatorname*{tr}(\mathbf{S}\mathbf{\Omega})+\sum_{k=1}^{K}P_{\lambda_{k}}(\mathbf{\Psi}_{k}),$

where $P_{\lambda_{k}}(\cdot)$ is a penalty function indexed by the tuning
parameter $\lambda_{k}$ and is applied elementwise to the off-diagonal
elements of $\mathbf{\Psi}_{k}$. Popular choices for $P_{\lambda_{k}}(\cdot)$
include the lasso penalty (Tibshirani, 1996), the adaptive lasso penalty (Zou,
2006), the SCAD penalty (Fan & Li, 2001), and the MCP penalty (Zhang et al.,
2010).

###  2.3 The Sylvester Generating Equation

Wang et al. (2020) proposed a Sylvester graphical model that uses the
Sylvester tensor equation to define a generative process for the underlying
multivariate tensor data. The Sylvester tensor equation has been studied in
the context of finite-difference discretization of high-dimensional elliptical
partial differential equations (Grasedyck, 2004; Kressner & Tobler, 2010). Any
solution $\bm{\mathcal{X}}$ to such a PDE must have the (discretized) form:

$\displaystyle\sum_{k=1}^{K}\bm{\mathcal{X}}\times_{k}\mathbf{\Psi}_{k}=\bm{\mathcal{T}}$
$\displaystyle\Longleftrightarrow\Big{(}\bigoplus_{k=1}^{K}\mathbf{\Psi}_{k}\Big{)}\operatorname{vec}(\bm{\mathcal{X}})=\operatorname{vec}(\bm{\mathcal{T}}).$
(2)

where $\bm{\mathcal{T}}$ is the driving source on the domain, and
$\bigoplus_{k=1}^{K}\mathbf{\Psi}_{k}$ is a Kronecker sum of
$\mathbf{\Psi}_{k}$’s representing the discretized differential operators for
the PDE, e.g., Laplacian, Euler-Lagrange operators, and associated
coefficients. These operators are often sparse and structured.

For example, consider a physical process characterized as a function $u$ that
satisfies:

$\mathcal{D}u=f\quad\text{in}\quad\Omega,\quad
u(\Gamma)=0,\quad\Gamma=\partial\Omega.$

where $f$ is a driving process, e.g., a Wiener process (white Gaussian noise);
$\mathcal{D}$ is a differential operator, e.g, Laplacian, Euler-Lagrange;
$\Omega$ is the domain; and $\Gamma$ is the boundary of $\Omega$. After
discretization, this is equivalent to (ignoring discretization error) the
matrix equation

$\mathbf{D}\mathbf{u}=\mathbf{f}.$

Here, $\mathbf{D}$ is a sparse matrix since $\mathcal{D}$ is an infinitesimal
operator. Additionally, $\mathbf{D}$ admits Kronecker structure as a mixture
of Kronecker sums and Kronecker products.

The matrix $\mathbf{D}$ reduces to a Kronecker sum when $\mathcal{D}$ involves
no mixed derivatives. For instance, consider the Poisson equation in 2D, where
$u(x,y)$ on $[0,1]^{2}$ satisfies the elliptical PDE

$\mathcal{D}u=(\partial^{2}_{x}+\partial^{2}_{y})u=f.$

The Poisson equation governs many physical processes, e.g., electromagnetic
induction, heat transfer, convection, etc. A simple Euler discretization
yields $\mathbf{U}=(u(i,j))_{i,j}$, where $u(i,j)$ satisfies the local
equation (up to a constant discretization scale factor)

$\displaystyle 4u(i,j)$ $\displaystyle=u(i+1,j)+u(i-1,j)+u(i,j+1)$
$\displaystyle\quad+u(i,j-1)-f(i,j).$

Defining $\mathbf{u}=\text{vec}(\mathbf{U})$ and $\mathbf{A}$ (a tridiagonal
matrix)

$\mathbf{A}=\begin{bmatrix}-1&2&-1&&&\\\ &\ddots&\ddots&\ddots\\\
&&\ddots&\ddots&\ddots\\\ &&&-1&2&-1\end{bmatrix},$

then $(\mathbf{A}\oplus\mathbf{A})\mathbf{u}=\mathbf{f}$, which is the
Sylvester equation ($K=2$).

For the Poisson example, if the source $\mathbf{f}$ is a white noise random
variable, i.e., its covariance matrix is proportional to the identity matrix,
then the inverse covariance matrix of $\mathbf{u}$ has sparse square-root
factors, since
$\text{Cov}^{-1}(\mathbf{u})=(\mathbf{A}\oplus\mathbf{A})(\mathbf{A}\oplus\mathbf{A})^{T}$.
Other physical processes that are generated from differential equations will
also have sparse inverse covariance matrices, as a result of the sparsity of
general discretized differential operators. Note that similar connections
between continuous state physical processes and sparse “discretized”
statistical models have been established by Lindgren et al. (2011), who
elucidated a link between Gaussian fields and Gaussian Markov Random Fields
via stochastic partial differential equations.

The Sylvester generative (SG) model (2) leads to a tensor-valued random
variable $\bm{\mathcal{X}}$ with a precision matrix
$\mathbf{\Omega}=\Big{(}\bigoplus_{k=1}^{K}\mathbf{\Psi}_{k}\Big{)}^{2}$,
given that $\bm{\mathcal{T}}$ is white Gaussian. The Sylvester generating
factors $\mathbf{\Psi}_{k}$’s can be obtained via minimization of the
penalized negative log-pseudolikelihood

$\displaystyle\mathcal{L}_{\mathbf{\lambda}}(\mathbf{\Psi})=$
$\displaystyle-\frac{N}{2}\log|(\bigoplus_{k=1}^{K}\mathop{\text{diag}}(\mathbf{\Psi}_{k}))^{2}|$
(3)
$\displaystyle+\frac{N}{2}\operatorname*{tr}(\mathbf{S}\cdot(\bigoplus_{k=1}^{K}\mathbf{\Psi}_{k})^{2})+\sum_{k=1}^{K}\lambda_{k}\|\mathbf{\Psi}_{k}\|_{1,\text{off}}.$

This differs from the true penalized Gaussian negative log-likelihood in the
exclusion of off-diagonals of $\mathbf{\Psi}_{k}$’s in the log-determinant
term. (3) is motivated and derived directly using the Sylvester equation
defined in (2), from the perspective of solving a sparse linear system. This
maximum pseudolikelihood estimation procedure has been applied to vector-
variate Gaussian graphical models (see Khare et al. (2015) and references
therein). Detailed derivations and further discussions are provided in
Appendix A.

##  3 The SG-PALM Algorithm

Estimation of the generating parameters $\mathbf{\Psi}_{k}$’s of the SG model
is challenging since the sparsity penalty applies to the square root factors
of the precision matrix, which leads to a highly coupled likelihood function.
Wang et al. (2020) proposed an estimation procedure called SyGlasso, that
recovers only the off-diagonal elements of each Sylvester factor. This is a
deficiency in many applications where the factor-wise variances are desired.
Moreover, the convergence rate of the cyclic coordinate-wise algorithm used in
SyGlasso is unknown and the computational complexity of the algorithm is
higher than other sparse Glasso-type procedures. To overcome these
deficiencies, we propose a proximal alternating linearized minimization method
that is more flexible and versatile, called SG-PALM, for finding the minimizer
of (3). SG-PALM is designed to exploit structures of the coupled objective
function and yields simultaneous estimates for both off-diagonal and diagonal
entries.

The PALM algorithm was originally proposed to solve nonconvex optimization
problems with separable structures, such as those arising in nonnegative
matrix factorization (Xu & Yin, 2013; Bolte et al., 2014). Its efficacy in
solving convex problems has also been established, for example, in regularized
linear regression problems (Shefi & Teboulle, 2016), it was proposed as an
attractive alternative to iterative soft-thresholding algorithms (ISTA). The
SG-PALM procedure is summarized in Algorithm 1.

For clarity of notation we write

$\mathcal{L}_{\mathbf{\lambda}}(\mathbf{\Psi}_{1},\dots,\mathbf{\Psi}_{K})=H(\mathbf{\Psi}_{1},\dots,\mathbf{\Psi}_{K})+\sum_{k=1}^{K}G_{k}(\mathbf{\Psi}_{k}),$
(4)

where $H:\mathbb{R}^{d_{1}\times
d_{1}}\times\cdots\times\mathbb{R}^{d_{K}\times d_{K}}\rightarrow\mathbb{R}$
represents the log-determinant plus trace terms in (3) and
$G_{k}:\mathbb{R}^{d_{k}\times d_{k}}\rightarrow(-\infty,+\infty]$ represents
the penalty term in (3) for each axis $k=1,\dots,K$. For notational simplicity
we use $\mathbf{\Psi}$ (i.e., omitting the subscript) to denote the set
$\\{\mathbf{\Psi}_{k}\\}_{k=1}^{K}$ or the $K$-tuple
$(\mathbf{\Psi}_{1},\dots,\mathbf{\Psi}_{K})$ whenever there is no risk of
confusion. The gradient of the smooth function $H$ with respect to
$\mathbf{\Psi}_{k}$, $\nabla_{k}H(\mathbf{\Psi})$, is given by

$\displaystyle\mathop{\text{diag}}\Big{(}\Big{\\{}\operatorname*{tr}[(\mathop{\text{diag}}((\mathbf{\Psi}_{k})_{ii})+\bigoplus_{j\neq
k}\mathop{\text{diag}}(\mathbf{\Psi}_{j}))^{-1}]\Big{\\}}_{i=1}^{d_{k}}\Big{)}$
(5)
$\displaystyle\quad+\mathbf{S}_{k}\mathbf{\Psi}_{k}+\mathbf{\Psi}_{k}\mathbf{S}_{k}+2\sum_{j\neq
k}\mathbf{S}_{j,k}.$

Here, the first “$\mathop{\text{diag}}$” maps a $d_{k}$-vector to a
$d_{k}\times d_{k}$ diagonal matrix, the second one maps a scalar (i.e.,
$(\mathbf{\Psi}_{k})_{ii}$) to a $(\prod_{j\neq k}d_{j})\times(\prod_{j\neq
k}d_{j})$ diagonal matrix with the same elements, and the third operator maps
a symmetric matrix to a matrix containing only its diagonal elements. In
addition, we define:

$\displaystyle\mathbf{S}_{k}=\frac{1}{N}\sum_{i=1}^{N}\bm{\mathcal{X}}_{(k)}^{i}(\bm{\mathcal{X}}_{(k)}^{i})^{T},$
(6)
$\displaystyle\mathbf{S}_{j,k}=\frac{1}{N}\sum_{i=1}^{N}\mathbf{V}_{j,k}^{i}(\mathbf{V}_{j,k}^{i})^{T},$
$\displaystyle\mathbf{V}_{j,k}^{i}=\bm{\mathcal{X}}_{(k)}^{i}\Big{(}\mathbf{I}_{d_{1:j-1}}\otimes\mathbf{\Psi}_{j}\otimes\mathbf{I}_{d_{j:K}}\Big{)}^{T},\quad
j\neq k.$

A key ingredient of the PALM algorithm is a proximal operator associated with
the non-smooth part of the objective, i.e., $G_{k}$’s. In general, the
proximal operator of a proper, lower semi-continuous convex function $f$ from
a Hilbert space $\mathcal{H}$ to the extended reals $(-\infty,+\infty]$ is
defined by (Parikh & Boyd, 2014)

$\vspace{-1pt}\text{prox}_{f}(v)=\operatorname*{argmin}_{x\in\mathcal{H}}f(x)+\frac{1}{2}\|x-v\|^{2}_{2}$

for any $v\in\mathcal{H}$. The proximal operator well-defined as the
expression on the right-hand side above has a unique minimizer for any
function in this class. For $\ell_{1}$-regularized cases, the proximal
operator for the function $G_{k}$ is given by

$\text{prox}_{G_{k}}^{\lambda_{k}}(\mathbf{\Psi}_{k})=\mathop{\text{diag}}(\mathbf{\Psi}_{k})+\text{soft}(\mathbf{\Psi}_{k}-\mathop{\text{diag}}(\mathbf{\Psi}_{k}),\lambda_{k}),$
(7)

where the soft-thresholding operator
$\text{soft}_{\lambda}(x)=\text{sign}(x)\max(|x|-\lambda,0)$ has been applied
element-wise. For popular choices of non-convex $G_{k}$, the proximal
operators are derived in Appendix D.

0: Data tensor $\bm{\mathcal{X}}$, mode-$k$ Gram matrix $\mathbf{S}_{k}$,
regularizing parameter $\lambda_{k}$, backtracking constant $c\in(0,1)$,
initial step size $\eta_{0}$, initial iterate $\mathbf{\Psi}_{k}$ for each
$k=1,\dots,K$.

while not converged do

for $k=1,\dots,K$ do

Line search: Let $\eta^{t}_{k}$ be the largest element of
$\\{c^{j}\eta_{k,0}^{t}\\}_{j=1,\dots}$ such that condition (8) is satisfied.

Update:
$\mathbf{\Psi}_{k}^{t+1}\leftarrow\text{prox}_{\eta^{t}_{k}\lambda_{k}}^{G_{k}}\Big{(}\mathbf{\Psi}_{k}^{t}-\eta^{t}_{k}\nabla_{k}H(\mathbf{\Psi}_{i<k}^{t+1},\mathbf{\Psi}_{i\geq
k}^{t})\Big{)}$.

end for

Update initial step size: Compute Barzilai-Borwein step size
$\eta_{0}^{t+1}=\min_{k}\eta^{t+1}_{k,0}$, where $\eta^{t+1}_{k,0}$ is
computed via (9).

end while

Final iterates $\\{\mathbf{\Psi}_{k}\\}_{k=1}^{K}$.

Algorithm 1 SG-PALM

###  3.1 Choice of Step Size

In the absence of a good estimate of the blockwise Lipchitz constant, the step
size of each iteration of SG-PALM is chosen using backtracking line search,
which, at iteration $t$, starts with an initial step size $\eta_{0}^{t}$ and
reduces the size with a constant factor $c\in(0,1)$ until the new iterate
satisfies the sufficient descent condition:

$H(\mathbf{\Psi}_{i\leq k}^{t+1},\mathbf{\Psi}_{i>k}^{t})\leq
Q_{\eta^{t}}(\mathbf{\Psi}_{i\leq
k}^{t+1},\mathbf{\Psi}_{i>k}^{t};\mathbf{\Psi}_{i<k}^{t+1},\mathbf{\Psi}_{i\geq
k}^{t}).$ (8)

Here,

$\displaystyle
Q_{\eta}(\mathbf{\Psi}_{i<k},\mathbf{\Psi}_{k},\mathbf{\Psi}_{i>k};\mathbf{\Psi}_{i<k},\mathbf{\Psi}_{k}^{\prime},\mathbf{\Psi}_{i>k})$
$\displaystyle=H(\mathbf{\Psi}_{i<k},\mathbf{\Psi}_{k},\mathbf{\Psi}_{i>k})$
$\displaystyle+\operatorname*{tr}\Big{(}(\mathbf{\Psi}_{k}^{\prime}-\mathbf{\Psi}_{k})^{T}\nabla_{k}H(\mathbf{\Psi}_{i<k},\mathbf{\Psi}_{k},\mathbf{\Psi}_{i>k})\Big{)}$
$\displaystyle+\frac{1}{2\eta}\|\mathbf{\Psi}_{k}^{\prime}-\mathbf{\Psi}_{k}\|_{F}^{2}.$

The sufficient descent condition is satisfied with any $\frac{1}{\eta}=M_{k}$
and $M_{k}\geq L_{k}$, for any function that has a block-wise Lipschitz
gradient with constant $L_{k}$ for $k=1,\dots,K$. In other words, so long as
the function $H$ has block-wise gradient that is Lipschitz continuous with
some block Lipschitz constant $L_{k}>0$ for each $k$, then at each iteration
$t$, we can always find an $\eta^{t}$ such that the inequality in (8) is
satisfied. Indeed, we proved in Lemma C.1 in the Appendix that $H$ has the
desired properties. Additionally, in the proof of Theorem 4.2 we also showed
that the step size found at each iteration $t$ satisfies
$\frac{1}{\eta_{k}^{0}}\leq L_{k}\leq\frac{1}{\eta_{k}^{t}}\leq cL_{k}$.

In terms of the initialization, a safe step size (i.e., very small
$\eta_{0}^{t}$) often leads to slower convergence. Thus, we use the more
aggressive Barzilai-Borwein (BB) step (Barzilai & Borwein, 1988) to set a
starting $\eta_{0}^{t}$ at each iteration (see Appendix B for justifications
of the BB method). In our case, for each $k$, the step size is given by

$\eta_{k,0}^{t}=\frac{\|\mathbf{\Psi}_{k}^{t+1}-\mathbf{\Psi}_{k}^{t}\|_{F}^{2}}{\operatorname*{tr}(\mathbf{A})},$
(9)

where

$\displaystyle\mathbf{A}$
$\displaystyle=(\mathbf{\Psi}_{k}^{t+1}-\mathbf{\Psi}_{k}^{t})^{T}\times$
$\displaystyle(\nabla_{k}H(\mathbf{\Psi}_{i\leq
k}^{t+1},\mathbf{\Psi}_{i>k}^{t})-\nabla_{k}H(\mathbf{\Psi}_{i<k}^{t+1},\mathbf{\Psi}_{i\geq
k}^{t})).$

###  3.2 Computational Complexity

After pre-computing $\mathbf{S}_{k}$, the most significant computation for
each iteration in the SG-PALM algorithm is the sparse matrix-matrix
multiplications $\mathbf{S}_{k}\mathbf{\Psi}_{k}$ and $\mathbf{S}_{j,k}$ in
the gradient calculation. In terms of computational complexity, if $s_{k}$ is
the number of non-zeros per column in $\mathbf{\Psi}_{k}$, then the former and
latter can be computed using $O(s_{k}d_{k}^{2})$ and $O(N\sum_{j\neq
k}s_{j}d_{j}^{2})$ operations, respectively. Thus, each iteration of SG-PALM
can be computed using $O\Big{(}\sum_{k=1}^{K}(s_{k}d_{k}^{2}+N\sum_{j\neq
k}s_{j}d_{j}^{2})\Big{)}$ floating point operations, which is significantly
lower than competing methods.

For instance, other popular algorithms for tensor-variate graphical models,
such as the TG-ISTA presented in Greenewald et al. (2019) and the Tlasso
proposed in Lyu et al. (2019) both require inversion of $d_{k}\times d_{k}$
matrices, which is non-parallelizable and requires $O(d_{k}^{3})$ operations
for each $k$. In particular, TeraLasso’s TG-ISTA algorithm requires
$O\Big{(}Kd+\sum_{k=1}^{K}d_{k}^{3}\Big{)}$ operations. The TG-ISTA algorithm
requires matrix inversions that cannot easily exploit the sparsity of
$\mathbf{\Psi}_{k}$’s. In the sample-starved ultra-sparse setting ($N\ll d$
and $s_{k}\ll d_{k}$), the $O(N\sum_{j\neq k}s_{j}d_{j}^{2})$ terms in SG-PALM
are comparable to $O(Kd)$ in TG-ISTA, making it more appealing. The cyclic
coordinate-wise method proposed in (Wang et al., 2020) does not allow for
parallelization since it requires cycling through entries of each
$\mathbf{\Psi}_{k}$ in specified order. In contrast, SG-PALM can be
implemented in parallel to distribute the sparse matrix-matrix multiplications
because at no step do the algorithms require storing all dense matrices on a
single machine.

##  4 Convergence Analysis

In this section, we present the main convergence theorems. Detailed proofs are
included in the supplement. Here, we study the statistical convergence
behavior for the Sylvester graphical model with an $\ell_{1}$ penalty
function. The convergence behavior of the SG-PALM iterates is presented for
convex cases but similar convergence rate can be established for non-convex
penalties (see Appendix D).

We first establish statistical convergence of a global minimizer
$\hat{\mathbf{\Psi}}$ of (3) to its true value, denoted as
$\bar{\mathbf{\Psi}}$, under the correct statistical model.

######  Theorem 4.1.

Let $\mathcal{A}_{k}:=\\{(i,j):(\bar{\mathbf{\Psi}}_{k})_{i,j}\neq 0,i\neq
j\\}$ and $q_{k}:=|\mathcal{A}_{k}|$ for $k=1,\dots,K$. If
$N>O(\max_{k}q_{k}d_{k}\log d)$ and $d:=d_{N}=O(N^{\kappa})$ for some
$\kappa\geq 0$, and further, if the penalty parameter satisfies
$\lambda_{k}:=\lambda_{N,k}=O(\sqrt{\frac{d_{k}\log d}{N}})$ for all
$k=1,\dots,K$, then under conditions (A1-A3) in Appendix C.1, there exists a
constant $C>0$ such that for any $\eta>0$ the following events hold with
probability at least $1-O(\exp(-\eta\log d))$:

$\displaystyle\sum_{k=1}^{K}\|\text{offdiag}(\hat{\mathbf{\Psi}}_{k})-\text{offdiag}(\bar{\mathbf{\Psi}}_{k})\|_{F}$
$\displaystyle\leq C\sqrt{K}\max_{k}\sqrt{q_{k}}\lambda_{k}.$

Here $\text{offdiag}(\mathbf{\Psi}_{k})$ contains only the off-diagonal
elements of $\mathbf{\Psi}_{k}$. If further
$\min_{(i,j)\in\mathcal{A}_{k}}|(\bar{\mathbf{\Psi}}_{k})_{i,j}|\geq
2C\max_{k}\sqrt{q_{k}}\lambda_{k}$ for each $k$, then
sign($\hat{\mathbf{\Psi}}_{k}$)=sign($\bar{\mathbf{\Psi}}_{k}$).

Theorem 4.1 means that under regularity conditions on the true generative
model, and with appropriately chosen penalty parameters $\lambda_{k}$’s guided
by the theorem, one is guaranteed to recover the true structures of the
underlying Sylvester generating parameters $\mathbf{\Psi}_{k}$ for
$k=1,\dots,K$ with probability one, as the sample size and dimension grow.

We next turn to convergence of the iterates $\\{\mathbf{\Psi}^{t}\\}$ from SG-
PALM to a global optimum of (3).

######  Theorem 4.2.

The Let $\\{\mathbf{\Psi}^{(t)}\\}_{t\geq 0}$ be generated by SG-PALM. Then,
SG-PALM converges in the sense that

$\displaystyle\frac{\mathcal{L}_{\mathbf{\lambda}}(\mathbf{\Psi}^{(t+1)})-\min\mathcal{L}_{\mathbf{\lambda}}}{\mathcal{L}_{\mathbf{\lambda}}(\mathbf{\Psi}^{(t)})-\min\mathcal{L}_{\mathbf{\lambda}}}$
$\displaystyle\leq\Bigg{(}\frac{\alpha^{2}L_{\min}}{4Kc^{2}(\sum_{j=1}^{K}L_{j})^{2}+4c^{2}L_{\max}}+1\Bigg{)}^{-1},$

where $\alpha$, $L_{k},k=1,\dots,K$ are positive constants,
$L_{\min}=\min_{j}L_{j}$, $L_{\max}=\max_{j}L_{j}$, and $c\in(0,1)$ is the
backtracking constant defined in Algorithm 1.

Note that the term on the right hand side of the inequality above is strictly
less than $1$. This means that the SG-PALM algorithm converges linearly, which
is a strong results for a non-strongly convex objective (i.e.,
$\mathcal{L}_{\mathbf{\lambda}}$). Although similar convergence behaviors of
the PALM-type algorithms have been studied for other problems (Xu & Yin, 2013;
Bolte et al., 2014), such as nonnegative matrix/tensor factorization, the
analysis of this paper works for non-strongly block multi-convex objectives,
leveraging more recent analyses of multi-block PALM and a class of functions
satisfying the the Kurdyka - Łojasiewicz (KL) property (defined in Section C
of the Appendix). To the best of our knowledge, for first-order optimization
methods, our rate is faster than any other Gaussian graphical models having
non-strongly convex objectives (see Khare et al. (2015); Oh et al. (2014) and
references therein) and comparable with those having strongly-convex
objectives (see, for example, Guillot et al. (2012); Dalal & Rajaratnam
(2017); Greenewald et al. (2019)).

##  5 Experiments

Experiments in this section were performed in a system with 8-core Intel Xeon
CPU E5-2687W v2 3.40GHz equipped with 64GB RAM. Both SG-PALM and SyGlasso were
implemented in Julia v1.5 (https://github.com/ywa136/sg-palm). For real data
analyses, we used the Tlasso package implementation in R (Sun et al., 2016)
and the TeraLasso implementation in MATLAB
(https://github.com/kgreenewald/teralasso).

###  5.1 Synthetic Data

We first validate the convergence theorems discussed in the previous section
via simulation studies. Synthetic datasets were generated from true sparse
Sylvester factors $\\{\mathbf{\Psi}_{k}\\}_{k=1}^{K}$ where $K=\\{2,3\\}$ and
$d_{k}=\\{16,32,64,128\\}$ for all $k$. Instances of the random matrices used
here have uniformly random sparsity patterns with edge densities (i.e., the
proportion of non-zero entries) ranging from $0.1\%-30\%$ on average over all
$\mathbf{\Psi}_{k}$’s. For each $d$ and edge density combination, random
samples of size $N=\\{10,100,1000\\}$ were tested. For comparison, the initial
iterates, convergence criteria were matched between SyGlasso and SG-PALM.
Highlights of the results in run times are summarized in Table 1.

Table 1: Run time comparisons (in seconds with N/As indicating those exceeding $24$ hour) between SyGlasso and SG-PALM on synthetic datasets with different dimensions, sample sizes, and densities of the generating Sylvester factors. Note that the proposed SG-PALM has average speed-up ratios ranging from $1.5$ to $10$ over SyGlasso. $d$ |  $N$ |  NZ% | SyGlasso | SG-PALM
---|---|---|---|---
iter sec |  iter sec
$128^{2}$ |  $10^{1}$ |  $1.20$ |  $17$ $138.5$ |  $46$ $5.8$
$24.0$ |  $20$ $169.3$ |  $48$ $6.2$
$10^{2}$ |  $1.30$ |  $21$ $211.3$ |  $50$ $12.6$
$27.0$ |  $30$ $303.6$ |  $47$ $21.9$
$10^{3}$ |  $1.30$ |  $21$ $2045.8$ |  $50$ $80.1$
$25.0$ |  $47$ $4782.7$ |  $51$ $373.1$
$16^{3}$ |  $10^{1}$ |  $0.11$ |  $9$ $4.6$ |  $11$ $4.5$
$4.10$ |  $9$ $5.1$ |  $32$ $5.1$
$10^{2}$ |  $0.21$ |  $8$ $8.8$ |  $11$ $5.4$
$2.60$ |  $8$ $10.8$ |  $35$ $7.2$
$10^{3}$ |  $0.26$ |  $8$ $82.4$ |  $12$ $14.3$
$3.40$ |  $10$ $99.2$ |  $37$ $33.5$
$32^{3}$ |  $10^{1}$ |  $0.13$ |  $10$ $191.2$ |  $19$ $7.3$
$7.50$ |  $17$ $304.8$ |  $42$ $10.2$
$10^{2}$ |  $0.46$ |  $9$ $222.4$ |  $24$ $28.9$
$7.00$ |  $17$ $395.2$ |  $41$ $48.5$
$10^{3}$ |  $0.10$ |  $9$ $1764.8$ |  $22$ $226.4$
$6.90$ |  $19$ $3789.4$ |  $41$ $473.9$
$64^{3}$ |  $10^{1}$ |  $0.65$ |  $10$ $583.7$ |  $42$ $91.3$
$14.5$ |  $22$ $952.2$ |  $47$ $119.0$
$10^{2}$ |  $0.62$ |  $9$ $6683.7$ |  $41$ $713.9$
$14.4$ |  $21$ $15607.2$ |  $48$ $1450.9$
$10^{3}$ |  $0.85$ | N/A |  $39$ $6984.4$
$14.0$ | N/A |  $48$ $12968.7$

Convergence behavior of SG-PALM is shown in Figure 1 (a) for the datasets with
$d_{k}=32$, $N=\\{10,100\\}$, and edge densities roughly around $5\%$ and
$20\%$, respectively. Geometric convergence rate of the function value gaps
under Theorem 4.2 can be verified from the plot. Note an acceleration in the
convergence rate (i.e., a steeper slope) near the optimum, which is suggested
by the “localness” of the KL property of the objective function close to its
global optimum. Further for the same datasets, in Figure 1 (b), SG-PALM graph
recovery performances is illustrated, where the Matthew’s Correlation
Coefficients (MCC) is plotted against run time. Here, MCC is defined by

$\text{MCC}=\frac{\text{TP}\times\text{TN}-\text{FP}\times\text{FN}}{\sqrt{(\text{TP}+\text{FP})(\text{TP}+\text{FN})(\text{TN}+\text{FP})(\text{TN}+\text{FN})}},$

where TP is the number of true positives, TN the number of true negatives, FP
the number of false positives, and FN the number of false negatives of the
estimated edges (i.e., non-zero elements of $\mathbf{\Psi}_{k}$’s). An MCC of
$1$ represents a perfect prediction, $0$ no better than random prediction and
$-1$ indicates total disagreement between prediction and observation. The
results validate the statistical accuracy under Theorem 4.1. It also shows
that SG-PALM outperforms SyGlasso (indicated by blue/red solid dots) within
the same time budget.

(a) Cost gap vs. Iteration

(b) MCC vs. Run time

Figure 1: Convergence of SG-PALM algorithm under datasets with varying sample
sizes (solid and dashed) generated via matrices with different sparsity (red
and blue). The function value gaps on log-scale (left) verifies the geometric
convergence rate in all cases and the MCC over time (right) demonstrates the
algorithm’s accuracy and efficiency. Note that the SG-PALM reached almost
perfect recoveries (i.e., MCC of $1$) within $20$ seconds in all cases. In
comparison, SyGlasso (big solid dots with line-range) was only able to achieve
at lower MCCs for lower sample-size cases within $30$ seconds.

###  5.2 Solar Flare Imaging Data

A solar flare occurs when magnetic energy that has built up in the solar
atmosphere is suddenly released. Such events strongly influence space weather
near the Earth. Therefore, reliable predictions of these flaring events are of
great interest. Recent work (Chen et al., 2019; Jiao et al., 2019; Sun et al.,
2019) has shown the promise of machine learning methods for early forecasting
of these events using imaging data from the solar atmosphere. In this work, we
illustrate the viability of the SG-PALM algorithm for solar flare prediction
using data acquired by multiple instruments: the Solar Dynamics Observatory
(SDO)/Helioseismic and Magnetic Imager (HMI) and SDO/Atmospheric Imaging
Assembly (AIA). It is evident that these data contain information about the
physical processes that govern solar activities (see Appendix E for detailed
data descriptions).

The data samples are summarized in $d_{1}\times d_{2}\times d_{3}\times d_{4}$
tensors with $q=d_{1}\cdot d_{2}\cdot d_{3}=50\cdot 100\cdot 7=35000$ and
$p=d_{4}=13$. The first two modes represent the images’ heights and widths,
the third mode represents the HMI/AIA components/channels, and the last mode
represents the length of the temporal window. Previous studies (Chen et al.,
2019; Jiao et al., 2019) found that the time series of solar images from the
SDO/HMI data provide useful information for distinguishing strong solar flares
of M/X class from weak flares of A/B class roughly 24 to 12 hours prior to the
flare event. Thus, in this study we use a $13$-hour temporal window recorded
with $1$-hour cadence, prior to the occurrence of a solar flare. The task is
to predict the $p$th frame using the frames in each of the $p-1$ previous
hours (i.e., one hour ahead prediction). Each observation is a video with full
dimension $d=pq$, and each $p$-dimensional observation vector is formed by
concatenating the $p$ time-consecutive $q$-dimensional vectors (vectorization
of the matrices representing pixels of the multichannel images) without
overlapping the time segments. The training set contains two types (B- and MX-
class flares) of active regions producing flares. Each is distinguished by the
flaring intensities, and there are a total of $186$ B flares and $48$ MX
flares. Forward linear predictors were constructed using estimated precision
matrices in a multi-output least squares regression setting. Specifically, we
constructed the linear predictor of a frame from the $p-1$ previous frames in
the same video:

$\hat{\mathbf{y}}_{t}=-\mathbf{\Omega}_{2,2}^{-1}\mathbf{\Omega}_{2,1}\mathbf{y}_{t-1:t-(p-1)},$
(10)

where $\mathbf{y}_{t-1:t-(p-1)}\in\mathbb{R}^{(p-1)q}$ is the stacked set of
pixel values from the previous $p-1$ time instances and
$\mathbf{\Omega}_{2,1}\in\mathbb{R}^{q\times(p-1)q}$ and
$\mathbf{\Omega}_{2,2}\in\mathbb{R}^{q\times q}$ are submatrices of the
$pq\times pq$ estimated precision matrix:

$\hat{\mathbf{\Omega}}=\begin{pmatrix}\mathbf{\Omega}_{1,1}&\mathbf{\Omega}_{1,2}\\\
\mathbf{\Omega}_{2,1}&\mathbf{\Omega}_{2,2}\end{pmatrix}.$

The predictors were tested on the data containing flares observed from
different active regions than those in training set, so that the predictor has
never “seen” the frames that it attempts to predict, corresponding to $117$
observations of which $93$ are B-class flares and $24$ are MX-class flares.
Figure 2 shows the root mean squared error normalized by the difference
between maximum and minimum pixels (NRMSE) over the testing samples, for the
forecasts based on the SG-PALM estimator, TeraLasso estimator (Greenewald et
al., 2019), Tlasso estimator (Lyu et al., 2019), and IndLasso estimator. Here,
the TeraLasso and the Tlasso are estimation algorithms for a KS and a KP
tensor precision matrix model, respectively; the IndLasso denotes an estimator
obtained by applying independent and separate $\ell_{1}$-penalized regressions
to each pixel in $\mathbf{y}_{t}$. The SG-PALM estimator was implemented using
a regularization parameter
$\lambda_{N}=C_{1}\sqrt{\frac{\min(d_{k})\log(d)}{N}}$ for all $k$ with the
constant $C_{1}$ chosen by optimizing the prediction NRMSE on the training set
over a range of $\lambda$ values parameterized by $C_{1}$. The TeraLasso
estimator and the Tlasso estimator were implemented using
$\lambda_{N,k}=C_{2}\sqrt{\frac{\log(d)}{N\prod_{i\neq k}d_{i}}}$ and
$\lambda_{N,k}=C_{3}\sqrt{\frac{\log(d_{k})}{Nd}}$ for $k=1,2,3$,
respectively, with $C_{2},C_{3}$ optimized in a similar manner. Each sparse
regression in the IndLasso estimator was implemented and tuned independently
with regularization parameters chosen from a grid via cross-validation.

We observe that SG-PALM outperforms all three other methods, indicated by
NRMSEs across pixels. Figure 3 depicts examples of predicted images, comparing
with the ground truth. The SG-PALM estimates produced most realistic image
predictions that capture the spatially varying structures and closely
approximate the pixel values (i.e., maintaining contrast ratios). The latter
is important as the flares are being classified into weak (B-class) and strong
(MX-class) categories based on the brightness of the images, and stronger
flares are more likely to lead to catastrophic events, such as those damaging
spacecrafts. Lastly, we compare run times of the SG-PALM algorithm for
estimating the precision matrix from the solar flare data with SyGlasso. Table
2 in Appendix E illustrates that the SG-PALM algorithm converges faster in
wallclock time. Note that in this real dataset, which is potentially non-
Gaussian, the convergence behavior of the algorithms is different compare to
synthetic examples. Nonetheless, SG-PALM enjoys an order of magnitude speed-up
over SyGlasso.

Avg. NRMSE = $0.0379$, $0.0386$, $0.0579$, $0.1628$ (from left to right)
---
AR B
Avg. NRMSE = $0.0620$, $0.0790$, $0.0913$, $0.1172$ (from left to right)
AR M/X
Figure 2: Comparison of the SG-PALM, Tlasso, TeraLasso, IndLasso performances
measured by NRMSE in predicting the last frame of $13$-frame video sequences
leading to B- and MX-class solar flares. The NRMSEs are computed by averaging
across testing samples and AIA channels for each pixel. 2D images of NRMSEs
are shown to indicate that certain areas on the images (usually associated
with the most abrupt changes of the magnetic field/solar atmosphere) are
harder to predict than the rest. SG-PALM achieves the best overall NRMSEs
across pixels. B flares are generally easier to predict due to both a larger
number of samples in the training set and smoother transitions from frame to
frame within a video (see the supplemental material for details). Predicted
examples - B vs. M/X
---
AR B
AR B
AR M/X
AR M/X
Figure 3: Examples of one-hour ahead prediction of the first two AIA channels
of last frames of $13$-frame videos, leading to B- (first two rows) and MX-
class (last two rows) flares, produced by the SG-PALM, Tlasso, TeraLasso,
IndLasso algorithms, comparing to the real image (far left column). Note that
in general linear forward predictors tend to underestimate the contrast ratio
of the images. The proposed SG-PALM produced the best-quality images in terms
of both the spatial structures and contrast ratios.See the supplemental
material for examples of predicted images from the HMI instrument.

###  5.3 Physical Interpretability

To explain the advantages of the proposed model over other similar models
(e.g., Tlasso, TeraLasso), we provide further discussions here on the
connection between the Sylvester generating model and PDEs. Consider the 2D
spatio-temporal process $u(\mathbf{x},t)$:

$\partial u/\partial t=\theta\sum_{i=1}^{2}\partial^{2}u/\partial
x_{i}^{2}+\epsilon\sum_{i=1}^{2}\partial u/\partial x_{i},$ (11)

where $\theta,\epsilon$ are positive real (unknown) coefficients. This is the
basic form of a class of parabolic and hyperbolic PDEs, the Convection-
Diffusion equation that generalizes the Poisson equation presented in Section
2 by incorporating temporal evolution. These equations are closely related to
the Navier-Stokes equation commonly used in stochastic modelling for weather
and climate prediction. Coupled with Maxwell’s equations, they can be used to
model and study magneto-hydrodynamics (Roberts, 2006), which characterize
solar activities including flares.

After finite-difference discretization, Equation (11) is equivalent to the
Sylvester matrix equation
$\mathbf{A}_{\theta,\epsilon}\mathbf{U}_{t}+\mathbf{U}_{t}\mathbf{A}_{\theta,\epsilon}=\mathbf{U}_{t-1}$,
where $\mathbf{U}_{t}=(u((i,j),t))_{ij}$ and $\mathbf{A}_{\theta,\epsilon}$ is
a tridiagonal matrix with values that depend on the coefficients
$\theta,\epsilon$ and discretization step sizes. Assuming a linear Gaussian
state-space model for some observed process $\mathbf{X}_{t}$ governed by the
Convection-Diffusion dynamics:

$\displaystyle\mathbf{A}_{\theta,\epsilon}\mathbf{U}_{t}+\mathbf{U}_{t}\mathbf{A}_{\theta,\epsilon}=\mathbf{U}_{t-1},$
$\displaystyle\mathbf{X}_{t}=\mathbf{U}_{t}+\mathbf{V}_{t},$

where $\mathbf{V}_{t}\sim\mathcal{N}(\mathbf{0},\sigma^{2}\mathbf{I})$ is some
time-invariant white noise. Then the precision matrix of the true process
$\mathbf{U}_{t}$ evolves as
$\mathbf{\Omega}_{t}=\mathbf{A}_{\theta,\epsilon}\mathbf{\Omega}_{t-1}\mathbf{A}_{\theta,\epsilon}^{T}+\sigma^{2}\mathbf{I}$.
Note that this is not necessarily sparse as assumed by the Sylvester graphical
model, but the steady-state precision matrix satisfies
$\mathbf{\Omega}_{\infty}=\mathbf{A}_{\theta,\epsilon}\mathbf{\Omega}_{\infty}\mathbf{A}_{\theta,\epsilon}^{T}+\sigma^{2}\mathbf{I}$,
which is indeed sparse because $\mathbf{A}_{\theta,\epsilon}$ is tridiagonal.
This strong connection between the Sylvester graphical model and the
underlying physical processes governing solar activities make the proposed
approach particularly suitable for the case study presented in the previous
section.

Additionally, the learned generating factors $\mathbf{A}_{\theta,\epsilon}$
could be further used to interpret physical processes that involve both
unknown structure and unknown parameters. Particularly, in Equation (11), the
coefficients $\theta$ (diffusion constant) and $\epsilon$ (convective
constant) affect the dynamics. Similarly, with the estimated Sylvester
generating factors ($\mathbf{\Psi}_{k}$’s), we are not only able to extract
the sparsity patterns of the discretized differential operators but also
estimate the coefficients of the underlying magneto-hydrodynamics equation for
solar flares. Therefore, the SG-PALM can be used as a data-driven method for
PDE parameter estimation from physical observations.

##  6 Conclusion

We proposed SG-PALM, a proximal alternating linearized minimization method for
solving a pseudo-likelihood based sparse tensor-variate Gaussian precision
matrix estimation problem. Geometric rate of convergence of the proposed
algorithm is established building upon recent advances in the theory of PALM-
type algorithms. We demonstrated that SG-PALM outperforms the coordinate-wise
minimization method in general, and in ultra-high dimensional settings SG-PALM
can be faster by at least an order of magnitude. A link between the Sylvester
generating equation underlying the graphical model and the Convection-
Diffusion type of PDEs governing certain physical processes was established.
This connection was illustrated on a novel astrophysics application, where
multi-instrument imaging datasets characterizing solar flare events were used.
The proposed methodology was able to robustly forward predict both the
patterns and intensities of the solar atmosphere, yielding potential insights
to the underlying physical processes that govern the flaring events.

#### Acknowledgements

The authors thank Zeyu Sun and Xiantong Wang for their help in pre-processing
the solar flare datasets. The research was partially supported by US Army
grant W911NF-15-1-0479 and NASA grant 80NSSC20K0600.

## References

  * Allen & Tibshirani (2010) Allen, G. I. and Tibshirani, R.  Transposable regularized covariance models with an application to missing data imputation.  _The Annals of Applied Statistics_ , 4(2):764, 2010.
  * Attouch & Bolte (2009) Attouch, H. and Bolte, J.  On the convergence of the proximal algorithm for nonsmooth functions involving analytic features.  _Mathematical Programming_ , 116(1-2):5–16, 2009\.
  * Barzilai & Borwein (1988) Barzilai, J. and Borwein, J. M.  Two-point step size gradient methods.  _IMA journal of numerical analysis_ , 8(1):141–148, 1988.
  * Besag (1977) Besag, J.  Efficiency of pseudolikelihood estimation for simple gaussian fields.  _Biometrika_ , pp. 616–618, 1977.
  * Bolte et al. (2014) Bolte, J., Sabach, S., and Teboulle, M.  Proximal alternating linearized minimization for nonconvex and nonsmooth problems.  _Mathematical Programming_ , 146(1-2):459–494, 2014.
  * Chen et al. (2019) Chen, Y., Manchester, W. B., Hero, A. O., Toth, G., DuFumier, B., Zhou, T., Wang, X., Zhu, H., Sun, Z., and Gombosi, T. I.  Identifying solar flare precursors using time series of sdo/hmi images and sharp parameters.  _Space Weather_ , 17(10):1404–1426, 2019.
  * Dai & Liao (2002) Dai, Y.-H. and Liao, L.-Z.  R-linear convergence of the Barzilai and Borwein gradient method.  _IMA Journal of Numerical Analysis_ , 22(1):1–10, 2002.
  * Dalal & Rajaratnam (2017) Dalal, O. and Rajaratnam, B.  Sparse gaussian graphical model estimation via alternating minimization.  _Biometrika_ , 104(2):379–395, 2017.
  * Fan & Li (2001) Fan, J. and Li, R.  Variable selection via nonconcave penalized likelihood and its oracle properties.  _Journal of the American statistical Association_ , 96(456):1348–1360, 2001.
  * Fletcher (2005) Fletcher, R.  On the Barzilai-Borwein method.  In _Optimization and control with applications_ , pp. 235–256. Springer, 2005.
  * Friedman et al. (2008) Friedman, J., Hastie, T., and Tibshirani, R.  Sparse inverse covariance estimation with the graphical lasso.  _Biostatistics_ , 9(3):432–441, 2008.
  * Galvez et al. (2019) Galvez, R., Fouhey, D. F., Jin, M., Szenicer, A., Muñoz-Jaramillo, A., Cheung, M. C., Wright, P. J., Bobra, M. G., Liu, Y., Mason, J., et al.  A machine-learning data set prepared from the nasa solar dynamics observatory mission.  _The Astrophysical Journal Supplement Series_ , 242(1):7, 2019.
  * Grasedyck (2004) Grasedyck, L.  Existence and computation of low kronecker-rank approximations for large linear systems of tensor product structure.  _Computing_ , 72(3-4):247–265, 2004.
  * Greenewald et al. (2019) Greenewald, K., Zhou, S., and Hero III, A.  Tensor graphical lasso (teralasso).  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 81(5):901–931, 2019.
  * Guillot et al. (2012) Guillot, D., Rajaratnam, B., Rolfs, B. T., Maleki, A., and Wong, I.  Iterative thresholding algorithm for sparse inverse covariance estimation.  _arXiv preprint arXiv:1211.2532_ , 2012.
  * Jiao et al. (2019) Jiao, Z., Sun, H., Wang, X., Manchester, W., Hero, A., Chen, Y., et al.  Solar flare intensity prediction with machine learning models.  _arXiv preprint arXiv:1912.06120_ , 2019.
  * Kalaitzis et al. (2013) Kalaitzis, A., Lafferty, J., Lawrence, N. D., and Zhou, S.  The bigraphical lasso.  In _International Conference on Machine Learning_ , pp. 1229–1237, 2013.
  * Karimi et al. (2016) Karimi, H., Nutini, J., and Schmidt, M.  Linear convergence of gradient and proximal-gradient methods under the polyak-łojasiewicz condition.  In _Joint European Conference on Machine Learning and Knowledge Discovery in Databases_ , pp. 795–811. Springer, 2016.
  * Khare et al. (2015) Khare, K., Oh, S.-Y., and Rajaratnam, B.  A convex pseudolikelihood framework for high dimensional partial correlation estimation with convergence guarantees.  _Journal of the Royal Statistical Society: Series B: Statistical Methodology_ , pp. 803–825, 2015.
  * Kolda & Bader (2009) Kolda, T. G. and Bader, B. W.  Tensor decompositions and applications.  _SIAM review_ , 51(3):455–500, 2009.
  * Kressner & Tobler (2010) Kressner, D. and Tobler, C.  Krylov subspace methods for linear systems with tensor product structure.  _SIAM journal on matrix analysis and applications_ , 31(4):1688–1714, 2010.
  * Leng & Tang (2012) Leng, C. and Tang, C. Y.  Sparse matrix graphical models.  _Journal of the American Statistical Association_ , 107(499):1187–1200, 2012.
  * Li & Pong (2018) Li, G. and Pong, T. K.  Calculus of the exponent of kurdyka–łojasiewicz inequality and its applications to linear convergence of first-order methods.  _Foundations of computational mathematics_ , 18(5):1199–1232, 2018.
  * Lindgren et al. (2011) Lindgren, F., Rue, H., and Lindström, J.  An explicit link between gaussian fields and gaussian markov random fields: the stochastic partial differential equation approach.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 73(4):423–498, 2011.
  * Lourenço & Takeda (2019) Lourenço, B. F. and Takeda, A.  Generalized subdifferentials of spectral functions over euclidean jordan algebras.  _arXiv preprint arXiv:1902.05270_ , 2019.
  * Lyu et al. (2019) Lyu, X., Sun, W. W., Wang, Z., Liu, H., Yang, J., and Cheng, G.  Tensor graphical model: Non-convex optimization and statistical inference.  _IEEE transactions on pattern analysis and machine intelligence_ , 2019.
  * Meinshausen et al. (2006) Meinshausen, N., Bühlmann, P., et al.  High-dimensional graphs and variable selection with the lasso.  _Annals of statistics_ , 34(3):1436–1462, 2006\.
  * Oh et al. (2014) Oh, S., Dalal, O., Khare, K., and Rajaratnam, B.  Optimization methods for sparse pseudo-likelihood graphical model selection.  In _Advances in Neural Information Processing Systems_ , pp. 667–675, 2014.
  * Parikh & Boyd (2014) Parikh, N. and Boyd, S.  Proximal algorithms.  _Foundations and Trends in optimization_ , 1(3):127–239, 2014.
  * Park et al. (2017) Park, S., Shedden, K., and Zhou, S.  Non-separable covariance models for spatio-temporal data, with applications to neural encoding analysis.  _arXiv preprint arXiv:1705.05265_ , 2017.
  * Raydan (1993) Raydan, M.  On the Barzilai and Borwein choice of steplength for the gradient method.  _IMA Journal of Numerical Analysis_ , 13(3):321–326, 1993.
  * Raydan (1997) Raydan, M.  The Barzilai and Borwein gradient method for the large scale unconstrained minimization problem.  _SIAM Journal on Optimization_ , 7(1):26–33, 1997\.
  * Roberts (2006) Roberts, B.  Slow magnetohydrodynamic waves in the solar atmosphere.  _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 364(1839):447–460, 2006.
  * Rudelson & Zhou (2017) Rudelson, M. and Zhou, S.  Errors-in-variables models with dependent measurements.  _Electronic Journal of Statistics_ , 11(1):1699–1797, 2017.
  * Shefi & Teboulle (2016) Shefi, R. and Teboulle, M.  On the rate of convergence of the proximal alternating linearized minimization algorithm for convex problems.  _EURO Journal on Computational Optimization_ , 4(1):27–46, 2016.
  * Sun et al. (2019) Sun, H., Manchester, W., Jiao, Z., Wang, X., and Chen, Y.  Interpreting lstm prediction on solar flare eruption with time-series clustering.  _arXiv preprint arXiv:1912.12360_ , 2019.
  * Sun et al. (2016) Sun, W. W., Wang, Z., Lyu, X., Liu, H., and Cheng, G.  _Tlasso: Non-Convex Optimization and Statistical Inference for Sparse Tensor Graphical Models_ , 2016.  URL https://CRAN.R-project.org/package=Tlasso.  R package version 1.0.1.
  * Tibshirani (1996) Tibshirani, R.  Regression shrinkage and selection via the lasso.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 58(1):267–288, 1996.
  * Tsiligkaridis et al. (2013) Tsiligkaridis, T., Hero III, A. O., and Zhou, S.  On convergence of kronecker graphical lasso algorithms.  _IEEE transactions on signal processing_ , 61(7):1743–1755, 2013.
  * Varin et al. (2011) Varin, C., Reid, N., and Firth, D.  An overview of composite likelihood methods.  _Statistica Sinica_ , pp. 5–42, 2011.
  * Wang & Ma (2007) Wang, Y. and Ma, S.  Projected Barzilai-Borwein method for large-scale nonnegative image restoration.  _Inverse Problems in Science and Engineering_ , 15(6):559–583, 2007.
  * Wang et al. (2020) Wang, Y., Jang, B., and Hero, A.  The sylvester graphical lasso (syglasso).  In _The 23rd International Conference on Artificial Intelligence and Statistics (AISTATS)_ , pp. 1943–1953, 2020.  URL http://proceedings.mlr.press/v108/wang20d.html.
  * Wen et al. (2010) Wen, Z., Yin, W., Goldfarb, D., and Zhang, Y.  A fast algorithm for sparse reconstruction based on shrinkage, subspace optimization, and continuation.  _SIAM Journal on Scientific Computing_ , 32(4):1832–1857, 2010.
  * Wright et al. (2009) Wright, S. J., Nowak, R. D., and Figueiredo, M. A.  Sparse reconstruction by separable approximation.  _IEEE Transactions on signal processing_ , 57(7):2479–2493, 2009.
  * Xu & Yin (2013) Xu, Y. and Yin, W.  A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion.  _SIAM Journal on imaging sciences_ , 6(3):1758–1789, 2013.
  * Yin & Li (2012) Yin, J. and Li, H.  Model selection and estimation in the matrix normal graphical model.  _Journal of multivariate analysis_ , 107:119–140, 2012\.
  * Zhang et al. (2010) Zhang, C.-H. et al.  Nearly unbiased variable selection under minimax concave penalty.  _The Annals of statistics_ , 38(2):894–942, 2010\.
  * Zhang (2020) Zhang, H.  New analysis of linear convergence of gradient-type methods via unifying error bound conditions.  _Mathematical Programming_ , 180(1):371–416, 2020\.
  * Zhou (2014) Zhou, S.  GEMINI: Graph estimation with matrix variate normal instances.  _The Annals of Statistics_ , 42(2):532–562, 2014\.
  * Zou (2006) Zou, H.  The adaptive lasso and its oracle properties.  _Journal of the American statistical association_ , 101(476):1418–1429, 2006.

## Appendix

  * Section A provides detailed derivation of the log-pseudolikelihood function.

  * Section B provides justifications for the Barzilai-Borwein step sizes implemented in Algorithm 1.

  * Section C provides detailed proofs of Theorems 4.1 and 4.2.

  * Section D discusses extensions of Algorithm 1 and its convergence properties to non-convex cases.

  * Section E provides additional details of the solar flare experiments.

##  Appendix A Derivation of the Log-Pseudolikelihood

By rewriting the Sylvester tensor equation defined in (2) element-wise, we
first observe that

$\displaystyle\left(\sum_{k=1}^{K}(\mathbf{\Psi}_{k})_{i_{k},i_{k}}\right)\bm{\mathcal{X}}_{i_{[1:K]}}$
(12) $\displaystyle=-\sum_{k=1}^{K}\sum_{j_{k}\neq
i_{k}}(\mathbf{\Psi}_{k})_{i_{k},j_{k}}\bm{\mathcal{X}}_{i_{[1:k]},j_{k},i_{[k+1:K]}}+\bm{\mathcal{T}}_{i_{[1:K]}}.$

Note that the left-hand side of (12) involves only the summation of the
diagonals of the $\mathbf{\Psi}_{k}$’s and the right-hand side is composed of
columns of $\mathbf{\Psi}_{k}$’s that exclude the diagonal terms. Equation
(12) can be interpreted as an autogregressive model relating the
$(i_{1},\dots,i_{K})$-th element of the data tensor (scaled by the sum of
diagonals) to other elements in the fibers of the data tensor. The columns of
$\mathbf{\Psi}_{k}$’s act as regression coefficients. The formulation in (12)
naturally leads to a pseudolikelihood-based estimation procedure (Besag, 1977)
for estimating $\mathbf{\Omega}$ (see also Khare et al. (2015) for how this
procedure applied to vector-variate Gaussian graphical model estimation). It
is known that inference using pseudo-likelihood is consistent and enjoys the
same $\sqrt{N}$ convergence rate as the MLE in general (Varin et al., 2011).
This procedure can also be more robust to model misspecification (e.g., non-
Gaussianity) in the sense that it assumes only that the sub-models/conditional
distributions (i.e., $\bm{\mathcal{X}}_{i}|\bm{\mathcal{X}}_{-i}$) are
Gaussian. Therefore, in practice, even if the data is not Gaussian, the
Maximum Pseudolikelihood Estimation procedure is able to perform reasonably
well. Wang et al. (2020) also studied a different model misspecification
scenario where the Kronecker product/sum and Sylvester structures are
mismatched for SyGlasso.

From (12) we can define the sparse least-squares estimators for
$\mathbf{\Psi}_{k}$’s as the solution of the following convex optimization
problem:

$\begin{aligned}
&\min_{\begin{subarray}{c}\mathbf{\Psi}_{k}\in\mathbb{R}^{d_{k}\times
d_{k}}\\\ k=1,\dots
K\end{subarray}}-N\sum_{i_{1},\dots,i_{K}}\log\bm{\mathcal{W}}_{i_{[1:K]}}\\\
&\qquad+\frac{1}{2}\sum_{i_{1},\dots,i_{K}}\|(I)+(II)\|_{2}^{2}+\sum_{k=1}^{K}P_{\lambda_{k}}(\mathbf{\Psi}_{k}).\end{aligned}\
$

where $P_{\lambda_{k}}(\cdot)$ is a penalty function indexed by the tuning
parameter $\lambda_{k}$ and

$\displaystyle(I)$
$\displaystyle=\bm{\mathcal{W}}_{i_{[1:K]}}\bm{\mathcal{X}}_{i_{[1:K]}}$
$\displaystyle(II)$ $\displaystyle=\sum_{k=1}^{K}\sum_{j_{k}\neq
i_{k}}(\mathbf{\Psi}_{k})_{i_{k},j_{k}}\bm{\mathcal{X}}_{i_{[1:k]},j_{k},i_{[k+1:K]}},$

with
$\bm{\mathcal{W}}_{i_{[1:K]}}:=\sum_{k=1}^{K}(\mathbf{\Psi}_{k})_{i_{k},i_{k}}$.

The optimization problem above can be put into the following matrix form:

$\displaystyle\min_{\begin{subarray}{c}\mathbf{\Psi}_{k}\in\mathbb{R}^{d_{k}\times
d_{k}}\\\ k=1,\dots K\end{subarray}}$
$\displaystyle-\frac{N}{2}\log|(\text{diag}(\mathbf{\Psi}_{1})\oplus\dots\oplus\text{diag}(\mathbf{\Psi}_{K}))^{2}|$
$\displaystyle+$
$\displaystyle\frac{N}{2}\operatorname*{tr}(\mathbf{S}(\mathbf{\Psi}_{1}\oplus\dots\oplus\mathbf{\Psi}_{K})^{2})+\sum_{k=1}^{K}P_{\lambda_{k}}(\mathbf{\Psi}_{k})$

where $\mathbf{S}\in\mathbb{R}^{d\times d}$ is the sample covariance matrix,
i.e.,
$\mathbf{S}=\frac{1}{N}\sum_{i=1}^{N}\operatorname{vec}(\bm{\mathcal{X}}^{i})\operatorname{vec}(\bm{\mathcal{X}}^{i})^{T}$.
Note that this is equivalent to the negative log-pseudolikelihood function
that approximates the $\ell_{1}$-penalized Gaussian negative log-likelihood in
the log-determinant term by including only the Kronecker sum of the diagonal
matrices instead of the Kronecker sum of the full matrices.

##  Appendix B The Barzilai-Borwein Step Size

The BB method has been proven to be very successful in solving nonlinear
optimization problems. In this section we outline the key ideas behind the BB
method, which is motivated by quasi-Newton methods. Suppose we want to solve
the unconstrained minimization problem

$\min_{x}f(x),$

where $f$ is differentiable. A typical iteration of quasi-Newton methods for
solving this problem is

$x_{t+1}=x_{t}-B_{t}^{-1}\nabla f(x_{t}),$

where $B_{t}$ is an approximation of the Hessian matrix of $f$ at the current
iterate $x_{t}$. Here, $B_{t}$ must satisfy the so-called secant equation:
$B_{t}s_{t}=y_{t}$, where $s_{t}=x_{t}-x_{t-1}$ and $y_{t}=\nabla
f(x_{t})-\nabla f(x_{t-1})$ for $t\geq 1$. It is noted that in to get
$B_{t}^{-1}$ one needs to solve a linear system, which may be computationally
expensive when $B_{t}$ is large and dense.

One way to alleviate this burden is to use the BB method, which replaces
$B_{t}$ by a scalar matrix $(1/\eta_{t})\mathbf{I}$. However, it is hard to
choose a scalar $\eta_{t}$ such that the secant equation holds with
$B_{t}=(1/\eta_{t})\mathbf{I}$. Instead, one can find $\eta_{t}$ such that the
residual of the secant equation, i.e., $\|(1/\eta_{t})s_{t}-y_{t}\|_{2}^{2}$,
is minimized, which leads to the following choice of $\eta_{t}$:

$\eta_{t}=\frac{\|s_{t}\|_{2}^{2}}{s_{t}^{T}y_{t}}.$

Therefore, a typical iteration of the BB method for solving the original
problem is

$x_{t+1}=x_{t}-\eta_{t}\nabla f(x_{t}),$

where $\eta_{t}$ is computed via the previous formula.

For convergence analysis, generalizations and variants of the BB method, we
refer the interested readers to Raydan (1993, 1997); Dai & Liao (2002);
Fletcher (2005) and references therein. BB method has been successfully
applied for solving problems arising from emerging applications, such as
compressed sensing (Wright et al., 2009), sparse reconstruction (Wen et al.,
2010) and image processing (Wang & Ma, 2007).

##  Appendix C Proofs of Theorems

###  C.1 Proof of Theorem 4.1

We first state the regularity conditions needed for establishing convergence
of the SG-PALM estimators $\\{\hat{\mathbf{\Psi}}_{k}\\}_{k=1}^{K}$ to their
true value $\\{\bar{\mathbf{\Psi}}_{k}\\}_{k=1}^{K}$.

(A1 - Subgaussianity) The data
$\bm{\mathcal{X}}^{1},\dots,\bm{\mathcal{X}}^{N}$ are i.i.d subgaussian random
tensors, that is, $\operatorname{vec}(\bm{\mathcal{X}}^{i})\sim\mathbf{x}$,
where $\mathbf{x}$ is a subgaussian random vector in $\mathbb{R}^{d}$, i.e.,
there exist a constant $c>0$, such that for every
$\mathbf{a}\in\mathbb{R}^{d}$, $\mathbb{E}e^{\mathbf{a}^{T}x}\leq
e^{c\mathbf{a}^{T}\bar{\mathbf{\Sigma}}\mathbf{a}}$, and there exist
$\rho_{j}>0$ such that $\mathbb{E}e^{tx_{j}^{2}}\leq+\infty$ whenever
$|t|<\rho_{j}$, for $1\leq j\leq d$.

(A2 - Bounded eigenvalues) There exist constants
$0<\Lambda_{\min}\leq\Lambda_{\max}<\infty$, such that the minimum and maximum
eigenvalues of $\mathbf{\Omega}$ are bounded with
$\lambda_{\min}(\bar{\mathbf{\Omega}})=(\sum_{k=1}^{K}\lambda_{\max}(\mathbf{\Psi}_{k}))^{-2}\geq\Lambda_{\min}$
and
$\lambda_{\max}(\bar{\mathbf{\Omega}})=(\sum_{k=1}^{K}\lambda_{\min}(\mathbf{\Psi}_{k}))^{-2}\leq\Lambda_{\max}$.

(A3 - Incoherence condition) There exists a constant $\delta<1$ such that for
$k=1,\dots,K$ and all $(i,j)\in\mathcal{A}_{k}$

$|\bar{\mathcal{L}}_{ij,\mathcal{A}_{k}}^{{}^{\prime\prime}}(\bar{\mathbf{\Psi}})[\bar{\mathcal{L}}_{\mathcal{A}_{k},\mathcal{A}_{k}}^{{}^{\prime\prime}}(\bar{\mathbf{\Psi}})]^{-1}\text{sign}(\bar{\mathbf{\Psi}}_{\mathcal{A}_{k},\mathcal{A}_{k}})|\leq\delta,$

where for each $k$ and $1\leq i<j\leq d_{k}$, $1\leq k<l\leq d_{k}$,

$\bar{\mathcal{L}}_{ij,kl}^{{}^{\prime\prime}}(\bar{\mathbf{\Psi}}):=E_{\bar{\mathbf{\Psi}}}\Bigg{(}\frac{\partial^{2}\mathcal{L}(\mathbf{\Psi})}{\partial(\mathbf{\Psi}_{k})_{i,j}\partial(\mathbf{\Psi}_{k})_{k,l}}|_{\mathbf{\Psi}=\bar{\mathbf{\Psi}}}\Bigg{)},$

and

$\mathcal{L}(\mathbf{\Psi})=-\frac{N}{2}\log|(\bigoplus_{k=1}^{K}\mathop{\text{diag}}(\mathbf{\Psi}_{k}))^{2}|+\frac{N}{2}\operatorname*{tr}(\mathbf{S}\cdot(\bigoplus_{k=1}^{K}\mathbf{\Psi}_{k})^{2}).$

Given assumptions (A1-A3), the theorem follows from Theorem 3.3 in Wang et al.
(2020).

###  C.2 Proof of Theorem 4.2

We next turn to convergence of the iterates $\\{\mathbf{\Psi}^{t}\\}$ from SG-
PALM to a global optimum of (3). The proof leverages recent results in the
convergence of alternating minimization algorithms for non-strongly convex
objective (Bolte et al., 2014; Karimi et al., 2016; Li & Pong, 2018; Zhang,
2020). We outline the proof strategy:

  1. 1.

We establish Lipschitz continuity of the blockwise gradient
$\nabla_{k}H(\mathbf{\Psi})$ for $k=1,\dots,K$.

  2. 2.

We show that the objective function $\mathcal{L}_{\mathbf{\lambda}}$ satisfies
the Kurdyka - Łojasiewicz (KL) property. Further, it has a KL exponent of
$\frac{1}{2}$ (defined later in the proofs).

  3. 3.

The KL property (with exponent $\frac{1}{2}$) is equivalent to a generalized
Error Bound (EB) condition, which enables us to establish linear iterative
convergence of the objective function (3) to its global optimum.

######  Definition C.1 (Subdifferentials).

Let $f:\mathbb{R}^{d}\rightarrow(-\infty,+\infty]$ be a proper and lower
semicontinuous function. Its domain is defined by

$\text{dom}f:=\\{x\in\mathbb{R}^{d}:f(x)<+\infty\\}.$

If we further assume that $f$ is convex, then the subdifferential of $f$ at
$x\in\text{dom}f$ can be defined by

$\partial f(x):=\\{v\in\mathbb{R}^{d}:f(z)\geq f(x)+<v,z-x>,\forall
z\in\mathbb{R}^{d}\\}.$

The elements of $\partial f(x)$ are called subgradients of $f$ at $x$.

Denote the domain of $\partial f$ by $\text{dom}\partial
f:=\\{x\in\mathbb{R}^{d}:\partial f(x)\neq\emptyset\\}$. Then, if $f$ is
proper, semicontinuous, convex, and $x\in\text{dom}f$, then $\partial f(x)$ is
a nonempty closed convex set. In this case, we denote by $\partial^{0}f(x)$
the unique least-norm element of $\partial f(x)$ for $x\in\text{dom}\partial
f$, along with $\|\partial^{0}f(x)\|=+\infty$ for $x\notin\text{dom}\partial
f$. Points whose subdifferential contains $0$ are critical points, denoted by
$\textbf{crit}f$. For convex $f$, $\textbf{crit}f=\operatorname*{argmin}f$.

######  Definition C.2 (KL property).

Let $\Gamma_{c_{2}}$ stands for the class of functions
$\phi:[0,c_{2}]\rightarrow\mathbb{R}_{+}$ for $c_{2}>0$ with the properties:

  1. i.

$\phi$ is continuous on $[0,c_{2}]$;

  2. ii.

$\phi$ is smooth concave on $(0,c_{2})$;

  3. iii.

$\phi(0)=0,\phi^{\prime}(s)>0,\forall s\in(0,c_{2})$.

Further, for $x\in\mathbb{R}^{d}$ and any nonempty $Q\subset\mathbb{R}^{d}$,
define the distance function $d(x,Q):=\inf_{y\in Q}\|x-y\|$. Then, a function
$f$ is said to have the Kurdyka - Łojasiewicz (KL) property at point $x_{0}$,
if there exist $c_{1}>0$, a neighborhood $B$ of $x_{0}$, and
$\phi\in\Gamma_{c_{2}}$ such that for all

$x\in B(x_{0},c_{1})\cap\\{x:f(x_{0})<f(x)<f(x_{0})+c_{2}\\},$

the following inequality holds

$\phi^{\prime}\Big{(}f(x)-f(x_{0})\Big{)}\text{dist}(0,\partial f(x))\geq 1.$

If $f$ satisfies the KL property at each point of $\text{dom}\partial f$ then
$f$ is called a KL function.

We first present two lemmas that characterize key properties of the loss
function.

######  Lemma C.1 (Blockwise Lipschitzness).

The function $H$ is convex and continuously differentiable on an open set
containing $\text{dom}G$ and its gradient, is block-wise Lipschitz continuous
with block Lipschitz constant $L_{k}>0$ for each $k$, namely for all
$k=1,\dots,K$ and all
$\mathbf{\Psi}_{k},\mathbf{\Psi}_{k}^{\prime}\in\mathbb{R}^{d_{k}\times
d_{k}}$

$\displaystyle\|\nabla_{k}H(\mathbf{\Psi}_{i<k},\mathbf{\Psi}_{k},\mathbf{\Psi}_{i>k})-\nabla_{k}H(\mathbf{\Psi}_{i<k},\mathbf{\Psi}_{k}^{\prime},\mathbf{\Psi}_{i>k})\|$
$\displaystyle\leq L_{k}\|\mathbf{\Psi}_{k}-\mathbf{\Psi}_{k}^{\prime}\|,$

where $\nabla_{k}H$ denotes the gradient of $H$ with respect to
$\mathbf{\Psi}_{k}$ with all remaining $\mathbf{\Psi}_{i}$, $i\neq k$ fixed.
Further, the function $G_{k}$ for each $k=1,\dots,K$ is a proper lower
semicontinuous (lsc) convex function.

###### Proof.

For simplicity of notation, in this and the following proofs we use
$\mathbf{\Psi}$ (i.e., omitting the subscript) to denote the set
$\\{\mathbf{\Psi}_{k}\\}_{k=1}^{K}$ or the $K$-tuple
$(\mathbf{\Psi}_{1},\dots,\mathbf{\Psi}_{K})$ whenever there is no confusion.
Recall the blockwise gradient of the smooth part of the objective function $H$
with respect to $\mathbf{\Psi}_{k}$, for each $k=1,\dots,K$, is given by

$\displaystyle\nabla_{k}H(\mathbf{\Psi})$
$\displaystyle=\mathop{\text{diag}}\Big{(}\Big{[}\operatorname*{tr}\\{(\mathop{\text{diag}}((\mathbf{\Psi}_{k}))_{ii}+\bigoplus_{j\neq
k}\mathop{\text{diag}}(\mathbf{\Psi}_{j}))^{-1}\\}\quad
i=1:d_{k}\Big{]}\Big{)}$
$\displaystyle\quad+\mathbf{S}_{k}\mathbf{\Psi}_{k}+\mathbf{\Psi}_{k}\mathbf{S}_{k}+2\sum_{j\neq
k}\mathbf{S}_{j,k}.$

Then for $\mathbf{\Psi}_{k},\mathbf{\Psi}^{\prime}_{k}$,

$\displaystyle\|\mathbf{S}_{k}\mathbf{\Psi}_{k}+\mathbf{\Psi}_{k}\mathbf{S}_{k}+2\sum_{j\neq
k}\mathbf{S}_{j,k}-(\mathbf{S}_{k}\mathbf{\Psi}^{\prime}_{k}+\mathbf{\Psi}^{\prime}_{k}\mathbf{S}_{k}+2\sum_{j\neq
k}\mathbf{S}_{j,k})\|$
$\displaystyle=\|\mathbf{S}_{k}\mathbf{\Psi}_{k}+\mathbf{\Psi}_{k}\mathbf{S}_{k}-\mathbf{S}_{k}\mathbf{\Psi}^{\prime}_{k}-\mathbf{\Psi}^{\prime}_{k}\mathbf{S}_{k}\|$
$\displaystyle\leq
2\|\mathbf{S}_{k}\|\|\mathbf{\Psi}_{k}-\mathbf{\Psi}^{\prime}_{k}\|.$

To prove Lipschitzness of the remaining parts, we consider the case of $K=2$
for simplicity of notations. The arguments easily carry over cases of $K>2$.
In this case, denote $\mathbf{A}=(a_{ij}):=\mathbf{\Psi}_{1}$ and
$\mathbf{B}=(b_{kl}):=\mathbf{\Psi}_{2}$. Let
$f(\mathbf{A}):=\frac{\partial}{\partial\mathbf{A}}\log|\mathop{\text{diag}}(\mathbf{A}\oplus\mathbf{B})|$,
then

$f(\mathbf{A})-f(\mathbf{A}^{\prime})=\mathop{\text{diag}}\Big{(}\Big{[}\sum_{i=1}^{m_{2}}(a_{jj}+b_{ii})^{-1}-\sum_{i=1}^{m_{2}}(a^{\prime}_{jj}+b_{ii})^{-1}\quad
j=1,\dots,m_{1}\Big{]}\Big{)}$

and

$\displaystyle\|f(\mathbf{A})-f(\mathbf{A}^{\prime})\|_{F}$
$\displaystyle=\Bigg{(}\sum_{j=1}^{m_{1}}\Big{(}\sum_{i=1}^{m_{2}}(a_{jj}+b_{ii})^{-1}-\sum_{i=1}^{m_{2}}(a^{\prime}_{jj}+b_{ii})^{-1}\Big{)}^{2}\Bigg{)}^{1/2}$
$\displaystyle\leq\Bigg{(}\sum_{j=1}^{m_{1}}m_{2}\sum_{i=1}^{m_{2}}\Big{(}(a_{jj}+b_{ii})^{-1}-(a^{\prime}_{jj}+b_{ii})^{-1}\Big{)}^{2}\Bigg{)}^{1/2}$
$\displaystyle=\Big{(}m_{2}\sum_{j=1}^{m_{1}}\sum_{i=1}^{m_{2}}(c_{ji}^{-1}-(c^{\prime}_{ji})^{-1})^{2}\Big{)}^{1/2}$
$\displaystyle=\Big{(}m_{2}\sum_{j=1}^{m_{1}}\sum_{i=1}^{m_{2}}(c^{\prime}_{ji})^{-2}(c^{\prime}_{ji}-c_{ji})^{2}c_{ji}^{-2}\Big{)}^{1/2}$
$\displaystyle=\Big{(}m_{2}\sum_{j=1}^{m_{1}}(a_{jj}-a^{\prime}_{jj})^{2}\sum_{i=1}^{m_{2}}(c^{\prime}_{ji}c_{ji})^{-2}\Big{)}^{1/2}$
$\displaystyle\leq\Big{(}Cm_{2}\sum_{j=1}^{m_{1}}\sum_{i=1}^{m_{2}}(c^{\prime}_{ji}c_{ji})^{-2}\Big{)}^{1/2}\|\mathbf{A}-\mathbf{A}^{\prime}\|_{F},$

where the first inequality is due to Cauchy-Schwartz inequality; the third
line is due to $c_{ji}:=a_{jj}+b_{ii}$; and in the last inequality we upper-
bound each $(a_{jj}-a^{\prime}_{jj})^{2}$ by its maximum over all $j$, which
is absorbed in a constant $C$. Note that the first term in the last line above
is finite as long as the summations of the diagonal elements of the factors
$\mathbf{A}$ and $\mathbf{B}$ are finite, which is implied if the precision
matrix $\mathbf{\Omega}$ defined by the Sylvester generating equation as
$(\mathbf{A}\oplus\mathbf{B})^{2}$ has finite diagonal elements. This follows
from Theorem 3.1 of Oh et al. (2014), who proved that if a symmetric matrix
$\mathbf{\Omega}$ satisfying $\mathbf{\Omega}\in\mathcal{C}_{0}$, where

$\mathcal{C}_{0}=\Big{\\{}\mathbf{\Omega}|\mathcal{L}_{\mathbf{\lambda}}(\mathbf{\Omega})\leq\mathcal{L}_{\mathbf{\lambda}}(\mathbf{\Omega}^{(0)})=M\Big{\\}},$

and $\mathbf{\Omega}^{(0)}$ is an arbitrary initial point with a finite
function value $\mathcal{L}_{\mathbf{\lambda}}(\mathbf{\Omega}^{(0)}):=M$, the
diagonal elements of $\mathbf{\Omega}$ are bounded above and below by
constants which depend only on $M$, the regularization parameter
$\mathbf{\lambda}$, and the sample covariance matrix $\mathbf{S}$. Therefore,
we have

$\|f(\mathbf{A})-f(\mathbf{A}^{\prime})\|_{F}\leq\tilde{C}\|\mathbf{A}-\mathbf{A}^{\prime}\|_{F}$

for some constant $\tilde{C}\in(0,+\infty)$. Similarly, we can establish such
an inequality for $\mathbf{B}$, proving that the first term in $\nabla_{k}H$
is Lipschitz continuous. ∎

As a consequence of Lemma C.1, the gradient of $H$, defined by $\nabla
H=(\nabla_{1}H,\dots,\nabla_{K}H)$ is Lipschitz continuous on bounded subsets
$\mathbb{B}_{1}\times\cdots\times\mathbb{B}_{K}$ of $\mathbb{R}^{d_{1}\times
d_{1}}\times\cdots\times\mathbb{R}^{d_{K}\times d_{K}}$ with some constant
$L>0$, such that for all
$(\mathbf{\Psi}_{k},\mathbf{\Psi}_{k}^{\prime})\in\mathbb{B}_{k}\times\mathbb{B}_{k}$,

$\displaystyle\|(\nabla_{1}H(\mathbf{\Psi}_{1},\mathbf{\Psi}_{i>1})-\nabla_{1}H(\mathbf{\Psi}_{1}^{\prime},\mathbf{\Psi}_{i>1}^{\prime}),\dots,$
$\displaystyle\nabla_{K}H(\mathbf{\Psi}_{i<K}^{\prime},\mathbf{\Psi}_{K}^{\prime})-\nabla_{K}H(\mathbf{\Psi}_{i<K}^{\prime},\mathbf{\Psi}_{K}^{\prime}))\|$
$\displaystyle\leq
L\|(\mathbf{\Psi}_{1}-\mathbf{\Psi}_{1}^{\prime},\dots,\mathbf{\Psi}_{K}-\mathbf{\Psi}_{K}^{\prime})\|,$

and we have $L\leq\sum_{k=1}^{K}L_{k}$.

######  Lemma C.2 (KL property of $\mathcal{L}_{\mathbf{\lambda}}$).

The objective function $\mathcal{L}_{\mathbf{\lambda}}(\mathbf{\Psi})$ defined
in (3) satisfies the KL property. Further, $\phi$ in this case can be chosen
to have the form $\phi(s)=\alpha s^{1/2}$, where $\alpha$ is some positive
real number. Functions satisfying the KL property with this particular choice
of $\phi$ is said to have a KL exponent of $\frac{1}{2}$.

###### Proof.

This can be established in a few steps:

  1. 1.

It can be shown that the function (of $\mathbf{X}$)
$\operatorname*{tr}(\mathbf{S}\mathbf{X}^{2})+\|\mathbf{X}\|_{1,\text{off}}$
satisfies the KL property with exponent $\frac{1}{2}$ (Karimi et al., 2016).
We then apply the calculus rules of the KL exponent (compositions and
separable summations) studied in Li & Pong (2018) to prove that
$\operatorname*{tr}(\mathbf{S}(\bigoplus_{j}\mathbf{\Psi}_{j})^{2})$ and
$\sum_{j}\|\mathbf{\Psi}_{j}\|_{1,\text{off}}$ are also KL functions with
exponent $\frac{1}{2}$.

  2. 2.

The
$-\log\det\Big{(}\bigoplus_{j}\mathop{\text{diag}}(\mathbf{\Psi}_{j})\Big{)}$
term can be shown to be KL with exponent $\frac{1}{2}$ using a transfer
principle studied in Lourenço & Takeda (2019).

  3. 3.

Finally, using the calculus rules of KL exponent one more time, we combine the
first two results and establish that $\mathcal{L}_{\mathbf{\lambda}}$ has KL
exponent of $\frac{1}{2}$.

Karimi et al. (2016) proved that the following function, parameterized by some
symmetric matrix $\mathbf{X}$, satisfies the KL property with KL exponent
$\frac{1}{2}$:

$\operatorname*{tr}(\mathbf{S}\mathbf{X}^{2})+\|\mathbf{X}\|_{1,\text{off}}=\|\mathbf{A}\mathbf{X}\|_{F}^{2}+\|\mathbf{X}\|_{1,\text{off}}$

for $\mathbf{S}=\mathbf{A}\mathbf{A}^{T}$, even when $\mathbf{A}$ is not of
full rank.

We apply the calculus rules of the KL exponent studied in Li & Pong (2018) to
prove that
$\operatorname*{tr}(\mathbf{S}(\bigoplus_{j}\mathbf{\Psi}_{j})^{2})$ and
$\sum_{j}\|\mathbf{\Psi}_{j}\|_{1,\text{off}}$ are KL functions with exponent
$\frac{1}{2}$. Particularly, we observe that
$\operatorname*{tr}\Big{(}\mathbf{S}(\bigoplus_{j}\mathbf{\Psi}_{j})^{2}\Big{)}$
is the composition of functions
$\mathbf{X}\rightarrow\operatorname*{tr}(\mathbf{S}\mathbf{X})$ and
$(\mathbf{X}_{1},\dots,\mathbf{X}_{K})\rightarrow\bigoplus_{j}\mathbf{X}_{j}$;
and $\sum_{j}\|\mathbf{\Psi}_{j}\|_{1,\text{off}}$ is a separable block
summation of functions
$\mathbf{X}_{j}\rightarrow\|\mathbf{X}_{j}\|_{1,\text{off}}$.

Thus, by Theorem 3.2. (exponent for composition of KL functions) in Li & Pong
(2018), since the Kronecker sum operation is linear and hence continuously
differentiable, the trace function is KL with exponent $\frac{1}{2}$, and the
mapping
$(\mathbf{X}_{1},\dots,\mathbf{X}_{K})\rightarrow\bigoplus_{j}\mathbf{X}_{j}$
is clearly one to one, the function
$\operatorname*{tr}(\mathbf{S}(\bigoplus_{j}\mathbf{\Psi}_{j})^{2})$ has the
KL exponent of $\frac{1}{2}$. By Theorem 3.3. (exponent for block separable
sums of KL functions) in Li & Pong (2018), since the function
$\|\cdot\|_{1,\text{off}}$ is proper, closed, continuous on its domain, and is
KL with exponent $\frac{1}{2}$, the function
$\|\mathbf{X}_{j}\|_{1,\text{off}}$ is KL with an exponent of $\frac{1}{2}$.

It remains to prove that
$-\log\det\Big{(}\bigoplus_{j}\mathop{\text{diag}}(\mathbf{\Psi}_{j})\Big{)}$
is also a KL function with an exponent of $\frac{1}{2}$. By Theorem 30 in
Lourenço & Takeda (2019), if we have $f:\mathbb{R}^{r}\rightarrow\mathbb{R}$ a
symmetric function and $F:\mathcal{E}\rightarrow\mathbb{R}$ the corresponding
spectral function, the followings hold

  1. (i).

$F$ satisfies the KL property at $\mathbf{X}$ iff $f$ satisfies the KL
property at $\lambda(\mathbf{X})$, i.e., the eigenvalues of $\mathbf{X}$.

  2. (ii).

$F$ satisfies the KL property with exponent $\alpha$ iff $f$ satisfies the KL
property with exponent $\alpha$ at $\lambda(\mathbf{X})$.

Here, take
$f(\lambda(\mathbf{X})):=-\sum_{i=1}^{r}\log(\lambda_{i}(\mathbf{X}))$, and
$F(\mathbf{X}):=-\log\det(\mathbf{X})$ the corresponding spectral function.
Then, the function $f$ is symmetric since its value is invariant to
permutations of its arguments, and it is a strictly convex function in its
domain, so it satisfies the KL property with an exponent of $\frac{1}{2}$.
Therefore, $F$ satisfies the KL property with the same KL exponent of
$\frac{1}{2}$. Now, we apply the calculus rules for KL functions again. As
both the Kronecker sum and the $\mathop{\text{diag}}$ operators are linear, we
conclude that
$-\log\det\Big{(}\bigoplus_{j}\mathop{\text{diag}}(\mathbf{\Psi}_{j})\Big{)}$
is a KL function with an exponent of $\frac{1}{2}$.

Overall, we have that the negative log-pseudolikelihood function
$\mathcal{L}(\mathbf{\Psi})$ satisfies the KL property with an exponent of
$\frac{1}{2}$. ∎

Now we are ready to prove Theorem 4.2. We follow Zhang (2020) and divide the
proof into three steps.

##### Step 1.

We obtain a sufficient decrease property for the loss function $\mathcal{L}$
in terms of the squared distance of two successive iterates:

$\mathcal{L}(\mathbf{\Psi}^{(t)})-\mathcal{L}(\mathbf{\Psi}^{(t+1)})\geq\frac{L_{\min}}{2}\|\mathbf{\Psi}^{(t)}-\mathbf{\Psi}^{(t+1)}\|^{2}.$
(13)

Here and below,
$\mathbf{\Psi}^{(t+1)}:=(\mathbf{\Psi}_{1}^{(t+1)},\dots,\mathbf{\Psi}_{K}^{(t+1)})$
and $L_{\min}:=\min_{k}L_{k}$. First note that at iteration $t$, the line
search condition is satisfied for step size $\frac{1}{\eta_{k}^{(t)}}\geq
L_{k}$, where $L_{k}$ is the Lipschitz constant for $\nabla_{k}H$. Further, it
follows that for SG-PALM with backtracking one has for every $t\geq 0$ and
each $k=1,\dots,K$,

$\frac{1}{\eta_{k}^{(0)}}\leq\frac{1}{\eta_{k}^{(t)}}\leq cL_{k},$

where $c>0$ is the backtracking constant.

Then by Lemma 3.1 in Shefi & Teboulle (2016), we get

$\displaystyle\mathcal{L}(\mathbf{\Psi}^{(t)})-\mathcal{L}(\mathbf{\Psi}^{(t+1)})$
$\displaystyle\geq\frac{1}{2\eta_{\min}^{(t+1)}}\|\mathbf{\Psi}^{(t)}-\mathbf{\Psi}^{(t+1)}\|^{2}$
$\displaystyle\geq\frac{L_{\min}}{2}\|\mathbf{\Psi}^{(t)}-\mathbf{\Psi}^{(t+1)}\|^{2}$

for $\eta_{\min}^{(t)}:=\min_{k}\eta_{k}^{(t)}$.

##### Step 2.

By Lemma C.2, $\mathcal{L}$ satisfies the KL property with an exponent of
$\frac{1}{2}$. Then from Definition C.2, this suggests that at
$x=\mathbf{\Psi}^{t+1}$ and $f(x_{0})=\min\mathcal{L}$

$\|\partial^{0}\mathcal{L}(\mathbf{\Psi}^{t+1})\|\geq\alpha\sqrt{\mathcal{L}(\mathbf{\Psi}^{t+1})-\min\mathcal{L}},$
(14)

where $\alpha>0$ is a fixed constant defined in Lemma C.2. This property is
equivalent to the error bound condition,
($\partial^{0}\mathcal{L},\alpha,\Omega$)-(res-obj-EB), defined in Definition
5 in Zhang (2020), for $\Omega\subset\text{dom}\partial\mathcal{L}$. This is
strictly weaker than strong convexity (see Section 4 in Zhang (2020)).

At iteration $t+1$, there exists $\xi_{k}^{(t+1)}\in\partial
G_{k}(\mathbf{\Psi}_{k}^{(t+1)})$ satisfying the optimality condition:

$\nabla_{k}H(\mathbf{\Psi}_{i<k}^{(t+1)},\mathbf{\Psi}_{i\geq
k}^{(t)})+\frac{1}{\eta_{k}^{(t+1)}}(\mathbf{\Psi}_{k}^{(t+1)}-\mathbf{\Psi}_{k}^{(t)})+\xi_{k}^{(t+1)}=0.$

Let $\xi^{(t+1)}:=(\xi_{1}^{(t+1)},\dots,\xi_{K}^{(t+1)})$. Then,

$\nabla
H(\mathbf{\Psi}^{(t+1)})+\xi^{(t+1)}\in\partial\mathcal{L}(\mathbf{\Psi}^{(t+1)})$

and hence the error bound condition becomes

$\mathcal{L}(\mathbf{\Psi}^{(t+1)})-\min\mathcal{L}\leq\frac{\|\partial^{0}\mathcal{L}(\mathbf{\Psi}^{(t+1)})\|^{2}}{\alpha^{2}}\leq\frac{\|\nabla
H(\mathbf{\Psi}^{(t+1)})+\xi^{(t+1)}\|^{2}}{\alpha^{2}}.$

It follows that

$\displaystyle\|\nabla H(\mathbf{\Psi}^{(t+1)})+\xi^{(t+1)}\|^{2}$
$\displaystyle=\sum_{k=1}^{K}\|\nabla_{k}H(\mathbf{\Psi}^{(t+1)})-\nabla_{k}H(\mathbf{\Psi}_{i<k}^{(t+1)},\mathbf{\Psi}_{i\geq
k}^{(t)})-\frac{1}{\eta_{k}^{(t+1)}}(\mathbf{\Psi}_{k}^{(t+1)}-\mathbf{\Psi}_{k}^{(t)})\|^{2}$
$\displaystyle\leq\sum_{k=1}^{K}2\|\nabla_{k}H(\mathbf{\Psi}^{(t+1)})-\nabla_{k}H(\mathbf{\Psi}_{i<k}^{(t+1)},\mathbf{\Psi}_{i\geq
k}^{(t)})\|^{2}+\sum_{k=1}^{K}\frac{2}{(\eta_{k}^{(t+1)})^{2}}\|\mathbf{\Psi}_{k}^{(t+1)}-\mathbf{\Psi}_{k}^{(t)}\|^{2}$
$\displaystyle\leq\sum_{k=1}^{K}2\|\nabla H(\mathbf{\Psi}^{(t+1)})-\nabla
H(\mathbf{\Psi}_{i<k}^{(t+1)},\mathbf{\Psi}_{i\geq
k}^{(t)})\|^{2}+\sum_{k=1}^{K}\frac{2}{(\eta_{k}^{(t+1)})^{2}}\|\mathbf{\Psi}_{k}^{(t+1)}-\mathbf{\Psi}_{k}^{(t)}\|^{2}$
$\displaystyle\leq\sum_{k=1}^{K}2\Big{(}\sum_{j=1}^{K}\frac{1}{\eta_{j}^{(t+1)}}\Big{)}^{2}\|\mathbf{\Psi}^{(t+1)}_{i\geq
k}-\mathbf{\Psi}^{(t)}_{i\geq
k}\|^{2}+\sum_{k=1}^{K}\frac{2}{(\eta_{k}^{(t+1)})^{2}}\|\mathbf{\Psi}_{k}^{(t+1)}-\mathbf{\Psi}_{k}^{(t)}\|^{2}$
$\displaystyle\leq\Bigg{(}2Kc^{2}\Big{(}\sum_{j=1}^{K}L_{j}\Big{)}^{2}+2c^{2}L_{\max}\Bigg{)}\|\mathbf{\Psi}^{(t+1)}-\mathbf{\Psi}^{(t)}\|^{2}.$

Therefore, we get

$\mathcal{L}(\mathbf{\Psi}^{(t+1)})-\min\mathcal{L}\leq\frac{\Bigg{(}2Kc^{2}\Big{(}\sum_{j=1}^{K}L_{j}\Big{)}^{2}+2c^{2}L_{\max}\Bigg{)}}{\alpha^{2}}\|\mathbf{\Psi}^{(t+1)}-\mathbf{\Psi}^{(t)}\|^{2}.$
(15)

##### Step 3.

Combining (13) and (15), we have

$\displaystyle\mathcal{L}(\mathbf{\Psi}^{(t)})-\min\mathcal{L}$
$\displaystyle=\Big{(}\mathcal{L}(\mathbf{\Psi}^{(t)})-\mathcal{L}(\mathbf{\Psi}^{(t+1)})\Big{)}+\Big{(}\mathcal{L}(\mathbf{\Psi}^{(t+1)})-\min\mathcal{L}\Big{)}$
$\displaystyle\geq\frac{L_{\min}}{2}\|\mathbf{\Psi}^{(t)}-\mathbf{\Psi}^{(t+1)}\|^{2}+\Big{(}\mathcal{L}(\mathbf{\Psi}^{(t+1)})-\min\mathcal{L}\Big{)}$
$\displaystyle\geq\Bigg{(}\frac{\alpha^{2}L_{\min}}{4Kc^{2}(\sum_{j=1}^{K}L_{j})^{2}+4c^{2}L_{\max}}+1\Bigg{)}\Big{(}\mathcal{L}(\mathbf{\Psi}^{(t+1)})-\min\mathcal{L}\Big{)}.$

This completes the proof.

##  Appendix D SG-PALM with Non-Convex Regularizers

The estimation algorithm for non-convex regularizer is largely the same as
Algorithm 1, except with an additional term added to the gradient term.
Specifically, the updates are of the form

$\mathbf{\Psi}_{k}^{(t+1)}=\text{prox}_{\eta^{t}_{k}\lambda_{k}}^{\|\cdot\|_{1,\text{off}}}\Big{(}\mathbf{\Psi}_{k}^{t}-\eta^{t}_{k}\nabla_{k}\bar{H}(\mathbf{\Psi}_{i<k}^{t+1},\mathbf{\Psi}_{i\geq
k}^{t})\Big{)},$

where

$\bar{H}(\mathbf{\Psi})=H(\mathbf{\Psi})+\sum_{k=1}^{K}\sum_{i\neq
j}\Big{(}g_{\lambda_{k}}([\mathbf{\Psi}_{k}]_{i,j})-\lambda_{k}|[\mathbf{\Psi}_{k}]_{i,j}|\Big{)}.$

Here, the formulation covers a range of non-convex regularizations.
Particularly, the SCAD penalty (Fan & Li, 2001) with parameter $a>2$ is given
by

$g_{\lambda}(t)=\begin{cases}\lambda|t|,\quad\text{if}\quad|t|<\lambda\\\
-\frac{t^{2}-2a\lambda|t|+\lambda^{2}}{2(a-1)},\quad\text{if}\quad\lambda<|t|<a\lambda\\\
\frac{(a+1)\lambda^{2}}{2},\quad\text{if}\quad a\lambda<|t|,\end{cases}$

which is linear for small $|t|$, constant for large $|t|$, and a transition
between the two regimes for moderate $|t|$.

The MCP penalty (Zhang et al., 2010) with parameter $a>0$ is given by

$g_{\lambda}(t)=\text{sign}(t)\lambda\int_{0}^{|t|}\Big{(}1-\frac{z}{a\lambda}\Big{)}_{+}dz,$

which gives a smoother transition between the approximately linear region and
the constant region ($t>a\lambda$) as defined in SCAD.

The updates can also be written as

$\mathbf{\Psi}_{k}^{(t+1)}=\text{prox}_{\eta^{t}_{k}\lambda_{k}}^{\|\cdot\|_{1,\text{off}}}\Bigg{(}\mathbf{\Psi}_{k}^{t}-\eta^{t}_{k}\nabla_{k}\Big{(}H(\mathbf{\Psi}_{i<k}^{t+1},\mathbf{\Psi}_{i\geq
k}^{t})+Q_{\lambda_{k}}^{\prime}(\mathbf{\Psi}_{k})\Big{)}\Bigg{)},$

where $q_{\lambda}^{\prime}(t):=\frac{d}{dt}(g_{\lambda}(t)-\lambda|t|)$ for
$t\neq 0$ and $q_{\lambda}^{\prime}(0)=0$ and $Q_{\lambda}^{\prime}$ denotes
$q_{\lambda}^{\prime}$ applied elementwise to a matrix argument. These updates
can be inserted into the framework of Algorithm 1. The details are summarized
in Algorithm 2.

0: Data tensor $\bm{\mathcal{X}}$, mode-$k$ Gram matrix $\mathbf{S}_{k}$,
regularizing parameter $\lambda_{k}$, backtracking constant $c\in(0,1)$,
initial step size $\eta_{0}$, initial iterate $\mathbf{\Psi}_{k}$ for each
$k=1,\dots,K$.

while not converged do

for $k=1,\dots,K$ do

Line search: Let $\eta^{t}_{k}$ be the largest element of
$\\{c^{j}\eta_{k,0}^{t}\\}_{j=1,\dots}$ such that condition (8) is satisfied
for
$\mathbf{\Psi}_{k}^{t+1}=\text{prox}_{\eta^{t}_{k}\lambda_{k}}^{\|\cdot\|_{1,\text{off}}}\Bigg{(}\mathbf{\Psi}_{k}^{t}-\eta^{t}_{k}\nabla_{k}\Big{(}H(\mathbf{\Psi}_{i<k}^{t+1},\mathbf{\Psi}_{i\geq
k}^{t})+Q_{\lambda_{k}}^{\prime}(\mathbf{\Psi}_{k})\Big{)}\Bigg{)}$.

Update:
$\mathbf{\Psi}_{k}^{t+1}\longleftarrow\text{prox}_{\eta^{t}_{k}\lambda_{k}}^{\|\cdot\|_{1,\text{off}}}\Bigg{(}\mathbf{\Psi}_{k}^{t}-\eta^{t}_{k}\nabla_{k}\Big{(}H(\mathbf{\Psi}_{i<k}^{t+1},\mathbf{\Psi}_{i\geq
k}^{t})+Q_{\lambda_{k}}^{\prime}(\mathbf{\Psi}_{k})\Big{)}\Bigg{)}$.

end for

Next initial step size: Compute Barzilai-Borwein step size
$\eta_{0}^{t+1}=\min_{k}\eta^{t+1}_{k,0}$, where $\eta^{t+1}_{k,0}$ is
computed via (9).

end while

Final iterates $\\{\mathbf{\Psi}_{k}\\}_{k=1}^{K}$.

Algorithm 2 SG-PALM with non-convex regularizer

###  D.1 Convergence Property

Consider a sequence of iterate $\\{\mathbf{x}^{t}\\}_{t\in\mathbb{N}}$
generated by a generic PALM algorithm for minimizing some objective function
$f$. Specifically, assume

  * $(\mathcal{H}_{1})$ $\inf f>-\infty$.

  * $(\mathcal{H}_{2})$ The restriction of the function to its domain is a continuous function.

  * $(\mathcal{H}_{3})$ The function satisfies the KL property.

Then, as in Theorem 2 of Attouch & Bolte (2009), if this objective function
satisfying $(\mathcal{H}_{1}),(\mathcal{H}_{2}),(\mathcal{H}_{3})$ in addition
satisfies the KL property with

$\phi(s)=\alpha s^{1-\theta},$

where $\alpha>0$ and $\theta\in(0,1]$. Then, for $\mathbf{x}^{\ast}$ some
critical point of $f$, the following estimations hold

  1. (i).

If $\theta=0$ then the sequence of iterates converges to $\mathbf{x}^{\ast}$
in a finite number of steps.

  2. (ii).

If $\theta\in(0,\frac{1}{2}]$ then there exist $\omega>0$ and $\tau\in[0,1)$
such that $\|\mathbf{x}^{t}-\mathbf{x}^{\ast}\|\leq\omega\tau^{t}$.

  3. (iii).

If $\theta\in(\frac{1}{2},1)$ then there exist $\omega>0$ such that
$\|\mathbf{x}^{t}-\mathbf{x}^{\ast}\|\leq\omega
t^{-\frac{1-\theta}{1\theta-1}}$.

In the case of SG-PALM with non-convex regularizations, so long as the non-
convex $\mathcal{L}$ satisfies the KL property with an exponent in
$(0,\frac{1}{2}]$, the algorithm remains linearly convergent (to a critical
point). We argue that this is true for SG-PALM with MCP or SCAD penalty. Li &
Pong (2018) showed that penalized least square problems with such penalty
functions satisfy the KL property with an exponent of $\frac{1}{2}$. The proof
strategy for the convex case can be easily adopted, incorporating the KL
results for MCP and SCAD in Li & Pong (2018), to show that the new
$\mathcal{L}$ still has KL exponent of $\frac{1}{2}$. Therefore, SG-PALM with
MCP or SCAD penalty converges linearly in the sense outlined above.

##  Appendix E Additional Details of the Solar Flare Experiments

###  E.1 HMI and AIA Data

The Solar Dynamics Observatory (SDO)/Helioseismic and Magnetic Imager (HMI)
data characterize solar variability including the Sun’s interior and the
various components of magnetic activity; the SDO/Atmospheric Imaging Assembly
(AIA) data contain a set of measurements of the solar atmosphere spectrum at
various wavelengths. In general, HMI produces data that is particularly useful
in determining the mechanisms of solar variability and how the physical
processes inside the Sun that are related to surface magnetic field and
activity. AIA contains structural information about solar flares, and the the
high AIA pixel values are correlated with the flaring intensities. We are
interested in examining if combination of multiple instruments enhances our
understanding of the solar flares, comparing to the case of single instrument.
Both HMI and AIA produce multi-band (or multi-channel) images, for this
experiment we use all three channels of the HMI images and
$9.4,13.1,17.1,19.3$ nm wavelength channels of the AIA images. For a detailed
descriptions of the instruments and all channels of the images, see
https://en.wikipedia.org/wiki/Solar_Dynamics_Observatory and the references
therein. Furthermore, for training and testing involved in this study, we used
the data described in (Galvez et al., 2019), which are further pre-processed
HMI and AIA imaging data for machine learning methods.

###  E.2 Classification of Solar Flares/Active Regions (AR)

The classification system for solar flares uses the letters A, B, C, M or X,
according to the peak flux in watts per square metre ($W/m^{2}$) of X-rays
with wavelengths $100$ to $800$ picometres ($1$ to $8$ angstroms), as measured
at the Earth by the GOES spacecraft
(https://en.wikipedia.org/wiki/Solar_flare#Classification). Here, A usually
refers to a “quite” region, which means that the peak flux of that region is
not high enough to be classified as a real flare; B usually refers a “weak”
region, where the flare is not strong enough to have impact on spacecrafts,
earth, etc; and M or X refers to a “strong” region that is the most
detrimental. Differentiating between a weak and a strong flare/region ahead of
time is a fundamental task in space physics and has recently attracted
attentions from the machine learning community (Chen et al., 2019; Jiao et
al., 2019; Sun et al., 2019). In our study, we also focus on B and M/X flares
and attempt to predict the videos that lead to either one of these two types
of flares.

###  E.3 Run Time Comparison

We compare run times of the SG-PALM algorithm for estimating the precision
matrix from the solar flare data with SyGlasso. Table 2 illustrates that the
SG-PALM algorithm converges faster in wallclock time. Note that in this real
dataset, which is potentially non-Gaussian, the convergence behavior of the
algorithms is different compare to synthetic examples. Nonetheless, SG-PALM
enjoys an order of magnitude speed-up over SyGlasso.

Table 2: Run time (in seconds) comparisons between SyGlasso and SG-PALM on solar flare data for different regularization parameters. Note that the SG-PALM is an order of magnitude faster that SyGlasso. $\lambda$ | SyGlasso | SG-PALM
---|---|---
iter sec |  iter sec
$0.28$ |  $47$ $5772.1$ |  $89$ $583.7$
$0.41$ |  $43$ $5589.0$ |  $86$ $583.4$
$0.54$ |  $45$ $5673.7$ |  $85$ $568.8$
$0.67$ |  $42$ $5433.0$ |  $77$ $522.6$
$0.79$ |  $39$ $4983.2$ |  $82$ $511.4$
$0.92$ |  $40$ $5031.9$ |  $72$ $498.0$
$1.05$ |  $39$ $4303.7$ |  $76$ $452.2$
$1.18$ |  $41$ $4234.7$ |  $64$ $437.6$
$1.30$ |  $40$ $4039.5$ |  $58$ $406.9$
$1.43$ |  $35$ $3830.7$ |  $64$ $364.9$

###  E.4 Examples of Predicted Magnetogram Images

Figure 4 depicts examples of the predicted HMI channels by SG-PALM. We observe
that the proposed method was able to reasonably capture various components of
the magnetic field and activity. Note that the spatial behaviors of the HMI
components are quite different from those of AIA channels, that is, the
structures tend to be less smooth and continuous (e.g., separated holes and
bright spots) in HMI.

Predicted HMI examples - B vs. M/X
---
AR B
AR B
AR B
AR M/X
AR M/X
AR M/X
Figure 4: Examples of one-hour ahead prediction of the first three channels
(HMI components) of ending frames of $13$-frame videos, leading to B- (first
three rows) and MX-class (last three rows) flares, produced by the SG-PALM,
comparing to the real image (left column). Similarly to AIA predictions,
linear forward predictors tend to underestimate the contrast ratio of the
images. Nonetheless, the SG-PALM algorithm was able to both capture the
spatial structures and intensities of the underlying magnetic fields. Note
that the HMI images tend to be harder to predict, as indicated by the
increased number and decreased degree of smoothness of features in the images,
signifying the underlying magnetic activity on the solar surface.

###  E.5 Multi-instrument vs. Single Instrument Prediction

To illustrate the advantages of multi-instrument analysis, we compare the
NRMSEs between an AIA-only (i.e., last four channels of the dataset) and an
HMI&AIA (i.e., all seven channels of the dataset) study in predicting the last
frames of $13$-frame AIA videos, for each flare class, respectively, using the
proposed SG-PALM. The results are depicted in Figure 5, where the average,
standard deviation, and range of the NRMSEs across pixels are also shown for
each error image. By leveraging the cross-instrument correlation structure,
there is a $0.5\%-1\%$ drop in the averaged error rates and a $2\%-4\%$ drop
in the range of the errors.

AIA Avg. NRMSE $=0.0379$ (with HMI channels) AIA Avg. NRMSE $=0.0479$ (without
HMI channels)
---
AR B
AIA Avg. NRMSE $=0.0620$ (with HMI channels) AIA Avg. NRMSE $=0.0674$ (without
HMI channels)
AR MX
Figure 5: Comparison of the SG-PALM performance measured by NRMSE in
predicting the AIA channels (i.e., last four channels) of the ending frame of
$13$-frame videos leading to B- and MX-class solar flares, by using all
HMI&AIA channels (left column) and AIA-only channels (right column). The
NRMSEs are computed by averaging across both testing samples and channels for
each pixel. Note that there are improvements in both the averaged errors rates
and the uncertainty in those errors (i.e., range of the errors) by including
multi-instrument image channels.

###  E.6 Illustration of the Difficulty of Predictions for Two Flares Classes

We demonstrate the difficulty of forward predictions of video frames. Figure 6
depicts two different channels of multiple frames from two videos leading to
MX-class solar flares. Note that the current frame is the $13$th frame in the
sequence that we are trying to predict. We observe that the prediction task is
particularly difficult if there is a sudden transition of either the
brightness or spatial structure of the frames near the end of the video. These
sudden transitions are more frequent for MX flares than for B flares. In
addition, as MX flares are generally considered as rare events (i.e., less
frequent than B flares), it is harder for SG-PALM or related methods to learn
a common correlation structures from training data.

On the other hand, typical image sequences leading to B flares exhibit much
smoother transitions from frame to frame. As shown in Figure 7, the SG-PALM
was able to produce remarkably good predictions of the current frames.

Predicted examples - M/X
---
Figure 6: Examples of frames at various timestamps of videos preceding the
predictions of the last frames (last column) that lead to MX flares. Here, the
first two rows correspond to the same video as the last two rows in Figure 3.
Note that the prediction tasks are difficult in these two extreme cases, where
there are dramatic changes from the $12$th to the current ($13$th) frames.
Predicted examples - B
---
Figure 7: Examples of frames at various timestamps of videos preceding the
predictions of the last frames (last column) that lead to B flares. Here, the
first two rows correspond to the same video as the first two rows in Figure 3.
Note that the prediction tasks are easier than those illustrated in Figure 6,
since the transitions near the end of the videos are much smoother.

###  E.7 Illustration of the Estimated Sylvester Generating Factors

Figure 8 illustrates the patterns of the estimated Sylvester generating
factors ($\mathbf{\Psi}_{k}$’s) for each flare class. Here, the videos from
both classes appear to form Markov Random Fields, that is, each pixel only
depends on its close neighbors in space and time given all other pixels. This
is demonstrated by observing that the temporal or each of the spatial
generating factor, which can be interpreted as conditional dependence graph
for the corresponding mode, has its energies concentrate around the diagonal
and decay as the nodes move far apart (in space or time).

The spatial patterns are similar for different flares. Although the exact
spatial patterns are different from one frame to another, they always have
their energies being concentrated at certain region (i.e., the brightest spot)
that is usually close to the center of the images. This is due to the way how
these images were curated and pre-processed before analysis. On the other
hand, the temporal structures are quite different. Specifically, B flares tend
to have longer range dependencies, as the frames leading to these types flares
are smooth, which is consistent with results from the previous section.

(a) Estimated precision matrices - B flares

(b) Estimated precision matrices - M/X flares

Figure 8: Estimated spatial and two (longitude and latitude) temporal
Sylvester generating factors for B and MX solar flares, along with their off-
diagonal sparsity patterns (second row in each subplot). Both classes exhibit
autoregressive dependence structures (across time or space). Note the
significant difference in the temporal components, where the B flares exhibit
longer range dependency. This is consistent with the smooth transition
property of the corresponding videos as illustrated previously.