\section{Introduction}


In many real-world problems for statistical analysis, complex dependencies are found among samples
as well as among features.
For example, in tumor gene-expression data collected from cancer patients,
gene-expression levels are correlated both among genes in the same pathways and 
among patients with the same or similar cancer subtypes~\citep{breastcancersubtype}.
Other examples include multivariate time-series data such as temporally correlated stock prices 
for related
companies~\citep{stocktogether} or a sequence of images in video data~\citep{videopixel}.
Gaussian graphical models with $L_1$-regularization have been widely used to learn a 
sparse inverse covariance matrix that corresponds to a graph over 
features~\citep{glasso,quic,bigquic}.
However, they were limited in that samples were assumed to be independent. As an alternative,
a multivariate Gaussian distribution has been generalized to a matrix-variate Gaussian distribution,
where the sample and feature dependencies were modeled with two separate 
graphs.



Kronecker-product or Kronecker-sum operators have been used
to combine two matrices, each representing a graph over samples and a graph over features,
to form an inverse covariance matrix for matrix-variate Gaussian distribution.
The Kronecker product of two graphs leads to a hard-to-interpret dense 
graph and non-convex log-likelihood~\citep{lengandtang2012,YIN2012,Tsiligkaridis13,zhou2014}.
However, its bi-convexity provided a fast flip-flop optimization method.
In contrast, the Kronecker sum has the advantage of producing a sparse graph as the Cartesian
product of the two graphs and having a convex log-likelihood, but the existing optimization methods
do not scale to large datasets~\citep{biglasso, TeraLasso}. BiGLasso 
used GLasso 
as a subroutine to estimate one graph while fixing the other graph in each iteration~\citep{biglasso}.
TeraLasso significantly improved the scalability of BiGLasso with the gradient
descent method~\citep{TeraLasso}, but still did not scale to problems with more than a few hundred
samples and features.


Another main challenge with the Kronecker sum comes from
the unidentifiable parameters in the diagonals of the matrices for feature and sample graphs. 
BiGLasso did not estimate these unidentifiable parameters. 
TeraLasso employed a reparameterization to make the
parameters identifiable and projected the estimates to the reparameterized space in each iteration.
In regression, the Kronecker-sum inverse covariance has been used to model errors
in covariates, but the trace of one of the two graphs was assumed to be known~\citep{RudelsonZhou2017,Park2017,Zhang2020}.  
With the Kronecker product, the parameters are unidentifiable as well, but 
they can be made identifiable with a simple method after the optimization is complete~\citep{YIN2012}.


In this paper, we introduce eigen graphical Lasso (EiGLasso) that addresses these limiations of the
existing methods for estimating the Kronecker-sum inverse covariance matrix.
Our contribution is two fold. First,
we develop an efficient optimization method based on Newton's method
that empirically leads two to three orders-of-magnitude speed-up compared to the existing methods.
Second, we describe a simple new approach to identifying and estimating the unidentifiable parameters
that generalizes the strategy used in TeraLasso~\citep{TeraLasso}.






For an efficient optimization, we adopt the framework used in QUIC~\citep{quic,bigquic},
the state-of-the-art method for estimating sparse Gaussian graphical models,
and extend it to address additional challenges
that arise when combining two graphs with the Kronecker-sum operator.
Since the gradient and Hessian matrices are far larger with an inflated structure than in QUIC, 
to reduce the computation time,
we leverage the eigendecomposition of the parameters. 
This eigendecomposition leads to a strategy
for approximating the Hessian to further
improve the scalability. Extending the theoretical results on the convergence of QUIC, we show
EiGLasso achieves quadratic convergence with the exact Hessian and linear convergence with 
the approximate Hessian.
In our experiments, we show that 
EiGLasso with the approximate Hessian does not require significantly more iterations
than EiGLasso with the exact Hessian and achieves
two to three orders-of-magnitude speed-up compared to the state-of-the-art method, 
TeraLasso~\citep{TeraLasso}.

In addition, we introduce a new approach to determining the unidentifiable parameters
in Kronecker-sum inverse covariance estimation.
Our work is the first to point out that the unidentifiable parameters
are uniquely determined given the ratio of the traces
of the two parameter matrices for graphs. Building upon this observation, we show that all the quantities involved in 
Newton's method, including the gradient, Hessian, the descent directions, and the step size,
are uniquely determined regardless of the trace ratio. This allows us to perform the optimization
in the original space of the parameters without being concerned about the unidentifiability,
rather than in the reparameterized space as in TeraLasso.
Furthermore, we show that the parameters can be identified once with the given trace ratio after EiGLasso converges,
analogous to the simple approach used in the Kronecker-product inverse covariance 
estimation~\citep{YIN2012},
rather than in every iteration as in TeraLasso. 

Our preliminary work on EiGLasso appeared in \citet{yoon2020}. While
in our earlier work, we presented a flip-flop optimization method as in the Kronecker-product
inverse covariance estimation, in this work, we present EiGLasso with the exact Hessian,
along with EiGLasso with the approximate Hessian directly derived from the exact Hessian.
We investigate the properties of these Hessian matrices and analyze the convergence behavior
of EiGLasso. We perform a more extensive experimental evaluation of our method.

The rest of the paper is organized as follows. In Section \ref{sec:related}, we review
the previous work on statistical methods with the Kronecker-sum operator. 
We introduce our EiGLasso optimization method in Section \ref{sec:eiglasso}, study
the convergence behavior in Section \ref{sec:convg}, present experimental results
in Section \ref{sec:exp}, and conlude with future work in Section \ref{sec:conclusion}.



\section{Related Work}\label{sec:related}

A Gaussian distribution of a matrix random variable 
$\bm{Y}\in\mathbb{R}^{q\times p}$ for $q$ samples and $p$ features
with a Kronecker-sum inverse covariance matrix~\citep{biglasso} is given as
\begin{equation}\label{eqn:ksmatrixnormal}
    \tvec(\bm{Y}) \sim \mathcal{N}\left(\tvec(\bm{M}), \bm{\Omega}^{-1}\right),
\end{equation}
where $\bm{M}\in\mathbb{R}^{q\times p}$ is the mean and 
$\tvec(\cdot)$ is an operator that stacks the columns of a matrix into a vector. 
The $pq \times pq$ inverse covariance matrix $\bm{\Omega}$ in Eq. \eqref{eqn:ksmatrixnormal} 
is defined as the Kronecker-sum of two graphs, each given as a $p\times p$ matrix $\bm{\Theta}$ 
and a $q \times q$ matrix $\bm{\Psi}$, modeling dependencies across features and across samples, respectively:
\begin{equation*}
    \bm{\Omega} = \bm{\Theta}\oplus\bm{\Psi} = \bm{\Theta}\otimes\bm{I}_q + \bm{I}_p\otimes\bm{\Psi},
\end{equation*}
where $\otimes$ is the Kronecker-product operator and $\bm{I}_a$ is an $a\times a$ identity matrix.
A non-zero value in the $(i,j)$th element $[\bm{\Theta}]_{ij}$, $[\bm{\Psi}]_{ij}$, and $[\bm{\Omega}]_{ij}$
implies an edge between the $i$th and $j$th nodes in the corresponding graph.
The two graphs $\bm{\Theta}$ and $\bm{\Psi}$ are constrained to form
a positive-definite Kronecker-sum space as follows:
\begin{equation}\label{def:ksspace}
    \mathbb{KS}^{p,q}=
	\{ (\bm{\Theta},\bm{\Psi})| \bm{\Theta}\in\mathbb{S}^p, \bm{\Psi}\in\mathbb{S}^q, \bm{\Omega} = \bm{\Theta}\oplus\bm{\Psi} \in \mathbb{S}^{pq}_{++}\}.
\end{equation}
Then, $\bm{\Omega}$ models a graph over $pq$ nodes, where each node is associated with an observation
for the given sample and feature pair.
For the rest of this paper, we assume zero mean, i.e., $\tvec(\bm{M})=\bm{0}$ and focus on the 
inverse covariance structure.






Compared to the Kronecker product $\bm{\Theta}\otimes \bm{\Psi}$, 
one main advantage of the Kronecker sum 
is that as the Cartesian product of the two graphs $\bm{\Theta}$ and $\bm{\Psi}$, 
it leads to a more intuitively appealing sparse structure in the graph $\bm{\Omega}$ 
(Figure \ref{fig:collapseW}(a)). 
In $\bm{\Omega}$ with the Kronecker sum,
the sample graph $\bm{\Psi}$ is used only within the same feature
(diagonal blocks in $\bm{\Omega}$ in Figure \ref{fig:collapseW}(a)) 
and the feature graph $\bm{\Theta}$ is used only within the same sample
(diagonals of off-diagonal blocks in $\bm{\Omega}$ in Figure \ref{fig:collapseW}(a)), 
whereas the Kronecker-product leads to a dense graph, 
where the sample graph $\bm{\Psi}$ and the feature graph $\bm{\Theta}$ are used 
even across features and across sample samples.
Because of the sparse structure, 
the Kronecker-sum operator has been adopted in various other statistical methods to represent a 
sparse association between two sets of variables, such as
dependencies between genes and mutations in colorectal cancer \citep{Haupt2020} 
and between microRNAs and diseases \citep{Li2018}. 
It has been embedded in neural networks to model both row and 
column dependencies in a matrix~\citep{Zhang2018tensor,Gao2020} and 
in inputs for more general functions~\citep{Benzi2017}. 





\begin{figure}[t!]
    \centering
    \begin{tabular}{cc}
        \begin{subfigure}{0.482\linewidth}
            \includegraphics[width=\linewidth]{images/omega2.pdf}
            \caption*{$\qquad\qquad\qquad$(a)}
        \end{subfigure} &
        \begin{subfigure}{0.518\linewidth}
            \includegraphics[width=\linewidth]{images/Wcollapse6.pdf}
            \caption*{(b)$\quad\qquad\qquad\qquad$}
        \end{subfigure}
    \end{tabular}
    \caption{Illustration of the Kronecker-sum inverse covariance matrix and gradient computation in EiGLasso.
    (a) $\bm{\Theta}$ and $\bm{\Psi}$ (left) and their Kronecker sum $\bm{\Omega}=\bm{\Theta}\oplus\bm{\Psi}$ (right).
    $\bm{\Theta}$ appears in the diagonals of off-diagonal blocks of $\bm{\Omega}$, and $\bm{\Psi}$ appears 
    in the diagonal blocks of $\bm{\Omega}$. 
    (b) $\bm{W}=\bm{\Omega}^{-1}$ (left) and $\bm{W}_{\bm{\Theta}}$ and $\bm{W}_{\bm{\Psi}}$ in the gradient (right). 
    Computing $\bm{W}_{\bm{\Theta}}$ and $\bm{W}_{\bm{\Psi}}$ amounts to
    collapsing $\bm{W}$ based on the Kronecker-sum structure in Panel (a).
    The $(i,j)$th element of $\bm{W}_{\bm{\Theta}}$ is obtained by summing over the diagonals of the
    $(i,j)$th block (blue arrows), and
    $\bm{W}_{\bm{\Psi}}$ is obtained by summing over the blocks on the diagonals of $\bm{W}$ (yellow arrows).
    }
    \label{fig:collapseW}
\end{figure}


Given data $\{\bm{Y}^1,\,\ldots,\,\bm{Y}^n\}$, where $\bm{Y}^i \in \mathbb{R}^{q\times p}$ 
for $i=1,\,\ldots,\,n$ for $n$ observations, 
BiGLasso~\citep{biglasso} and TeraLasso~\citep{TeraLasso} obtained
a sparse estimate of $\bm{\Theta}$ and $\bm{\Psi}$
by minimizing the $L_1$-regularized negative log-likelihood of data:
\begin{equation} \label{eq:obj}
    \argmin_{\bm{\Theta},\, \bm{\Psi}}  f(\bm{\Theta},\bm{\Psi}) \quad \textrm{subject to} \quad\bm{\Theta}\oplus\bm{\Psi}\succ 0.
\end{equation}  
The objective $f(\bm{\Theta},\bm{\Psi})$ is
\begin{equation}
    f(\bm{\Theta},\bm{\Psi})= g(\bm{\Theta},\bm{\Psi}) + h(\bm{\Theta},\bm{\Psi}), \nonumber
\end{equation}
with the smooth log-likelihood function 
$g(\bm{\Theta},\bm{\Psi}) = q\tr(\bm{S}\bm{\Theta}) + p\tr(\bm{T}\bm{\Psi}) - \log|\bm{\Theta}\oplus\bm{\Psi}|$ and
non-smooth penalty $h(\bm{\Theta},\bm{\Psi}) = q \gamma_{\bm{\Theta}} \| \bm{\Theta} \|_{1,\text{off}} + p \gamma_{\bm{\Psi}} \| \bm{\Psi} \|_{1,\text{off}}$,
given the sample covariances $\bm{S} = \frac{1}{nq}\sum_{i=1}^n\bm{Y}^{iT}\bm{Y}^i$ and 
$\bm{T} = \frac{1}{np}\sum_{i=1}^n\bm{Y}^i\bm{Y}^{iT}$, the regularization parameters 
$\gamma_{\bm{\Theta}}$ and $\gamma_{\bm{\Psi}}$, and 
$\|\cdot\|_{1,\text{off}}$ 
for the $L_1$-regularization of the off-diagonal elements of the matrix.



One of the main challenges of solving Eq. \eqref{eq:obj} 
arises from the unidentifiable diagonals of $\bm{\Theta}$ and $\bm{\Psi}$ given $\bm{\Omega}$.
Although \citet{TeraLasso} showed that the objective in 
Eq. \eqref{eq:obj} has a unique global optimum  with respect to $\bm{\Omega}=\bm{\Theta}\oplus\bm{\Psi}$, 
this does not imply the uniqueness of the pair $(\bm{\Theta},\bm{\Psi})$ given $\bm{\Omega}$,
since the set in Eq. \eqref{def:ksspace} forms equivalence classes 
\begin{equation}\label{eqn:noniddiag}
    S[\bm{\Omega}] = \{ (\bm{\Theta},\bm{\Psi})\in \mathbb{KS}^{p,q}| 
	\left(\bm{\Theta}-c\bm{I}_p\right)\oplus\left(\bm{\Psi}+c\bm{I}_q\right)=\bm{\Omega},
		c\in \mathbb{R}\}.
\end{equation}
Thus, $\bm{\Theta}$ and $\bm{\Psi}$ whose diagonals are modified by a constant $c$ and $-c$ 
lead to the same $\bm{\Omega}$.
While BiGLasso did not estimate the diagonals of $\bm{\Theta}$ and $\bm{\Psi}$, 
TeraLasso formed an identifiable re-parameterization
 $\bm{\Omega} = \bar{\bm{\Theta}} \oplus \bar{\bm{\Psi}} + \tau \bm{I}_{pq}$ with
$\tau = \frac{\tr(\bm{\Theta} \oplus \bm{\Psi})}{pq}$, where
$\bar{\bm{\Theta}}$ and $\bar{\bm{\Psi}}$ are forced to have zero traces as follows:
\begin{align}\label{eq:reparam}
\bar{\bm{\Theta}}=\bm{\Theta} -\frac{\tr(\bm{\Theta})}{p}\bm{I}_p\qquad\text{and}\qquad\bar{\bm{\Psi}}=\bm{\Psi} -\frac{\tr(\bm{\Psi})}{q}\bm{I}_q,
\end{align}
Then, in each iteration of a gradient descent method,
TeraLasso projected the gradient to this re-parameterized space and distributed 
$\tau$ to $\bm{\Theta}$ and $\bm{\Psi}$ as $\bar{\bm{\Theta}} \gets \bar{\bm{\Theta}} + \frac{1}{2}\tau\bm{I}$ and 
$\bar{\bm{\Psi}} \gets \bar{\bm{\Psi}} + \frac{1}{2}\tau\bm{I}$. 
For efficient computation and projection of the gradient, TeraLasso
employed an eigendecomposition of $\bm{\Theta}$ and $\bm{\Psi}$:
$\bm{\Theta}=\bm{Q}_{\bm{\Theta}}\bm{\Lambda}_{\bm{\Theta}}\bm{Q}_{\bm{\Theta}}^T$ with
the eigenvector matrix $\bm{Q}_{\bm{\Theta}}$ 
and the diagonal eigenvalue matrix $\bm{\Lambda}_{\bm{\Theta}}$, and similarly, 
$\bm{\Psi} = \bm{Q}_{\bm{\Psi}} \bm{\Lambda}_{\bm{\Psi}} \bm{Q}_{\bm{\Psi}}^T$.
TeraLasso exploited the property that the eigendecomposition of the Kronecker sum of 
$\bm{\Theta}$ and $\bm{\Psi}$ is 
\begin{equation}\label{eq:omegaeigdecomp}
    \bm{\Theta} \oplus \bm{\Psi} = (\bm{Q}_{\bm{\Theta}}\otimes\bm{Q}_{\bm{\Psi}})(\bm{\Lambda}_{\bm{\Theta}}\oplus\bm{\Lambda}_{\bm{\Psi}})(\bm{Q}_{\bm{\Theta}}\otimes\bm{Q}_{\bm{\Psi}})^T.
\end{equation}
Then, the inverse of $\bm{\Theta} \oplus \bm{\Psi}$ can be obtained efficiently 
by inverting the diagonal eigenvalue matrix $\bm{\Lambda}_{\bm{\Theta}} \oplus \bm{\Lambda}_{\bm{\Psi}}$:
\begin{equation}\label{eq:Weigdecomp}
    \bm{W} = (\bm{\Theta} \oplus \bm{\Psi})^{-1} = (\bm{Q}_{\bm{\Theta}}\otimes\bm{Q}_{\bm{\Psi}})(\bm{\Lambda}_{\bm{\Theta}}\oplus\bm{\Lambda}_{\bm{\Psi}})^{-1}(\bm{Q}_{\bm{\Theta}}\otimes\bm{Q}_{\bm{\Psi}})^T.
\end{equation}
While TeraLasso is significantly faster than BiGLasso, its scalability is limited to graphs with
a few hundred nodes. 




\section{EiGLasso}\label{sec:eiglasso}

We introduce EiGLasso, an efficient method for estimating a sparse Kronecker-sum inverse covariance matrix. 
We begin by describing a simple new scheme for identifying the unidentifiable parameters in Section \ref{sec:tr}. 
We introduce EiGLasso, Newton's method with the exact Hessian in Section \ref{sec:exactH} 
and its modification with approximate Hessian in Section \ref{sec:approxH}. 
In Section \ref{sec:identify}, we describe a simple strategy for identifying the 
diagonal parameters during EiGLasso estimation.
In Section \ref{sec:block}, we discuss additional strategies for improving the computational efficiency. 


\subsection{Identifying Parameters with Trace Ratio}\label{sec:tr}

Theorem \ref{thm:uniquepair} below provides a simple approach to identifying the diagonal elements 
of $\bm{\Theta}$ and $\bm{\Psi}$, given $\bm{\Omega}$ and the ratio of the traces of $\bm{\Theta}$ and $\bm{\Psi}$.


\begin{theorem}\label{thm:uniquepair}
    Assume that the trace ratio of $\bm{\Psi}$ to $\bm{\Theta}$ is given as $\rho = \frac{\tr(\bm{\Psi})}{\tr(\bm{\Theta})}$.        
    Then, a Kronecker-sum matrix $\bm{\Omega}=\bm{\Theta}\oplus\bm{\Psi}$ can be mapped to a unique pair of symmetric matrices $(\bm{\Theta},\bm{\Psi})$ with the diagonals of $\bm{\Theta}$ and $\bm{\Psi}$ identified as
    \begin{equation}\label{eq:diagmapping}
    \begin{aligned}
        \diag(\bm{\Theta}) &= \frac{1}{q}\left[\diag\left( \sum_{i=1}^q(\bm{I}_p\otimes\bm{e}_{q,i})^T\bm{\Omega}(\bm{I}_p\otimes\bm{e}_{q,i}) \right)  - \frac{\rho}{q + \rho p}\tr(\bm{\Omega})\bm{I}\right],\\
        \diag(\bm{\Psi}) &= \frac{1}{p}\left[\diag\left( \sum_{j=1}^p(\bm{e}_{p,j}\otimes\bm{I}_q)^T\bm{\Omega}(\bm{e}_{p,j}\otimes\bm{I}_q) \right)  - \frac{1}{q + \rho p}\tr(\bm{\Omega})\bm{I}\right],
    \end{aligned}
    \end{equation}
    where $\bm{e}_{q,i}$ is a $q\times 1$ one-hot vector with a single $1$ at the $i$th entry and $0$'s elsewhere.
\end{theorem}
\begin{proof}
    Given $\rho=\frac{\tr(\bm{\Psi})}{\tr(\bm{\Theta})}$, $\bm{\Theta}$ and $\bm{\Psi}$ are unique, since for any $c\in\mathbb{R}\setminus\{0\}$,
        $\frac{\tr(\bm{\Psi})}{\tr(\bm{\Theta})} \neq \frac{\tr(\bm{\Psi}+c\bm{I}_q)}{\tr(\bm{\Theta}-c\bm{I}_p)}=\frac{\tr(\bm{\Psi}) + cq}{\tr(\bm{\Theta}) - cp}$. 
    From the definition of Kronecker sum, we notice
    \begin{align}\label{eq:omegaii}
        \sum_{i=(j-1)q + 1}^{jq} [\bm{\Omega}]_{ii} &= q[\bm{\Theta}]_{jj} + \tr(\bm{\Psi}),
    \end{align}
and 
    \begin{equation}\label{eq:tromega}
           \tr(\bm{\Omega}) =  q\tr(\bm{\Theta}) + p\tr(\bm{\Psi}).
    \end{equation}
    From the trace ratio and Eq. \eqref{eq:tromega}, we have $\tr(\bm{\Psi}) = \frac{\rho}{q + \rho p}\tr(\bm{\Omega})$. 
    Plugging this into Eq. \eqref{eq:omegaii}, we identify $[\bm{\Theta}]_{jj}$ as 
    \begin{align*}
        [\bm{\Theta}]_{jj} &= \frac{1}{q}\left(\sum_{i=(j-1)q + 1}^{jq} [\bm{\Omega}]_{ii} - \frac{\rho}{q + \rho p}\tr(\bm{\Omega})\right).
    \end{align*}
    The case for $\bm{\Psi}$ can be shown in a similar way.
   
   
\end{proof}
    

Theorem \ref{thm:uniquepair} suggests that it is possible to identify the diagonals 
of $\bm{\Theta}$ and $\bm{\Psi}$ by solving Eq. \eqref{eq:obj} with the linear constraint 
$\frac{\tr(\bm{\Psi})}{\tr(\bm{\Theta})}=\rho$:
\begin{align}
    \argmin_{\bm{\Theta},\bm{\Psi}}  f(\bm{\Theta},\bm{\Psi})  
    \quad \textrm{subject to} \quad\bm{\Theta}\oplus\bm{\Psi}\succ0, \;
\bm{R}^T
    \begin{bmatrix}
        \tvec(\bm{\Theta})\\
        \tvec(\bm{\Psi})
    \end{bmatrix} = 0,    \label{eq:objrho}
\end{align} 
where the second constraint with 
$\bm{R} = \begin{bmatrix} 
-\rho\tvec(\bm{I}_p)\\
\tvec(\bm{I}_q)
\end{bmatrix}$
is equivalent to $\frac{\tr(\bm{\Psi})}{\tr(\bm{\Theta})}=\rho$.
To handle this additional constraint explicitly, the substitution method or Newton's method 
for equality-constrained optimization problem could be used \citep{cvxbook}. 
Instead, in Sections \ref{sec:exactH} and \ref{sec:identify}, 
we show that in EiGLasso, because of the special problem structure, 
it is sufficient to solve Eq. \eqref{eq:obj} with Newton's method, ignoring the equality constraint, 
and to adjust the diagonals of $\bm{\Theta}$ and $\bm{\Psi}$ once after convergence to satisfy the constraint.



\subsection{EiGLasso with Exact Hessian} \label{sec:exactH}

We introduce EiGLasso for an efficient estimation of sparse $\bm{\Theta}$ and $\bm{\Psi}$ that form a Kronecker-sum inverse covariance matrix. 
We adopt the framework of QUIC~\citep{quic}, Newton's method for estimating
sparse Gaussian graphical models. In each iteration, QUIC found a Newton direction 
by minimizing the $L_1$-regularized second-order approximation
of the objective, and updated the parameters given this Newton direction and the step size
found by backtracking line search \citep{bertsekas1995,Tseng2009}. 
Using the same strategy, EiGLasso finds Newton directions, $D_{\bm{\Theta}}$ for $\bm{\Theta}$ and $D_{\bm{\Psi}}$ for $\bm{\Psi}$, 
by minimizing the following second-order approximation of the smooth part $g(\bm{\Theta},\bm{\Psi})$ of 
the objective in Eq. \eqref{eq:objrho} with the $L_1$-regularization:
\begin{align}
    (D_{\bm{\Theta}},D_{\bm{\Psi}}) = \argmin_{\Delta_{\bm{\Theta}},\Delta_{\bm{\Psi}}}\; \hat{g}(\Delta_{\bm{\Theta}},\Delta_{\bm{\Psi}}) + h(\bm{\Theta} + \Delta_{\bm{\Theta}}, \bm{\Psi} + \Delta_{\bm{\Psi}}) \label{eq:newton}\\
    \text{subject to } \bm{R}^T
    \begin{bmatrix}
        \tvec(\Delta_{\bm{\Theta}})\\
        \tvec(\Delta_{\bm{\Psi}})
    \end{bmatrix} = 0,
    \nonumber
\end{align}
where
\begin{equation*}\label{eq:g2compact}
    \hat{g}(\Delta_{\bm{\Theta}}, \Delta_{\bm{\Psi}}) = \tvec(\bm{G})^T
    \begin{bmatrix}
        \tvec(\Delta_{\bm{\Theta}})\\
        \tvec(\Delta_{\bm{\Psi}})
    \end{bmatrix} + \frac{1}{2}
    \begin{bmatrix}
        \tvec(\Delta_{\bm{\Theta}})\\
        \tvec(\Delta_{\bm{\Psi}})
    \end{bmatrix}^T
    \bm{H}
     \begin{bmatrix}
        \tvec(\Delta_{\bm{\Theta}})\\
        \tvec(\Delta_{\bm{\Psi}})
    \end{bmatrix}
\end{equation*}
with the gradient $\bm{G}$ and Hessian $\bm{H}$ 
\begin{gather}
    \tvec(\bm{G}) =
        \begin{bmatrix}
            \tvec\left(\bm{G}_{\bm{\Theta}} \right)\\ 
            \tvec\left(\bm{G}_{\bm{\Psi}} \right)
        \end{bmatrix}
        =
        \begin{bmatrix}
            \tvec\left(\nabla_{\bm{\Theta}} g(\bm{\Theta},\bm{\Psi})\right)\\ 
            \tvec\left(\nabla_{\bm{\Psi}} g(\bm{\Theta},\bm{\Psi})\right)
        \end{bmatrix},
   
    \label{eq:compactg}
    \\
    \bm{H} =
        \begin{bmatrix}
            \bm{H}_{\bm{\Theta}}  & \bm{H}_{\bm{\Theta}\bm{\Psi}}\\
            \bm{H}_{\bm{\Theta}\bm{\Psi}}^T  & \bm{H}_{\bm{\Psi}}
        \end{bmatrix}
        =
        \begin{bmatrix}
            \nabla^2_{\bm{\Theta}} g(\bm{\Theta},\bm{\Psi}) & \nabla_{\bm{\Theta}}\nabla_{\bm{\Psi}} g(\bm{\Theta},\bm{\Psi})\\
            \nabla_{\bm{\Theta}}\nabla_{\bm{\Psi}} g(\bm{\Theta},\bm{\Psi})^T & \nabla^2_{\bm{\Psi}} g(\bm{\Theta},\bm{\Psi})
        \end{bmatrix}.
   
    \label{eq:compacth}
\end{gather}
As we will show in Section \ref{sec:identify}, the objective above can be optimized without being
concerned about the equality constraint. Thus, 
the coordinate descent optimization can be used to solve Eq. \eqref{eq:newton} as in QUIC. 
Given the descent directions $D_{\bm{\Theta}}$ and $D_{\bm{\Psi}}$, 
we update the parameters as $\bm{\Theta} \gets \bm{\Theta} + \alpha D_{\bm{\Theta}}$ and $\bm{\Psi} \gets \bm{\Psi} + \alpha D_{\bm{\Psi}}$, 
with the step size $\alpha$ found by the line-search method in Algorithm \ref{alg:linesearch}. 



Two key challenges arise in a direct application of QUIC to Kronecker-sum inverse covariance estimation.
First, as we detail in the next section, the gradient $\bm{G}$ and Hessian $\bm{H}$ 
involves the computation of the $pq\times pq$ matrix $\bm{W}=(\bm{\Theta} \oplus \bm{\Psi})^{-1}$ significantly larger
than a $p\times p$ matrix required for a single graph in QUIC. 
Second, the constraint on the trace ratio in Eq. \eqref{eq:newton} needs to be considered 
to handle the unidentifiability of the diagonals of $\bm{\Theta}$ and $\bm{\Psi}$ 
during the optimization process.
In the rest of Section \ref{sec:eiglasso}, we describe how we address these challenges by leveraging 
the eigen structure of the model, and provide the details of the EiGLasso optimization outlined 
in Algorithm \ref{alg:eiglasso}.


\begin{algorithm}[t!] 
\caption{Line Search}
\label{alg:linesearch}
    \DontPrintSemicolon
    \SetKwInOut{Input}{input}\SetKwInOut{Output}{output}
    \Input{$0<\sigma<0.5$, $0<\beta<1$}
    \Output{step size $\alpha$}
    
    Initialize $\alpha = 1$\;
    
    \For{$i=0,1,\dots$}{
        Check the following conditions:\;
        \begin{enumerate}
            \item $(\bm{\Theta} + \alpha D_{\bm{\Theta}})\oplus(\bm{\Psi} + \alpha D_{\bm{\Psi}})\succ 0$,\;
            \item $f(\bm{\Theta} + \alpha D_{\bm{\Theta}},\bm{\Psi} + \alpha D_{\bm{\Psi}})\le f(\bm{\Theta},\bm{\Psi}) + \alpha\sigma\delta$, where 
        
    $\delta = \tvec(D)^T\tvec(\bm{G}) + h(\bm{\Theta}+D_{\bm{\Theta}},\bm{\Psi}+D_{\bm{\Psi}}) - h(\bm{\Theta},\bm{\Psi})$.
\;
        \end{enumerate}
        If satisfied, break. If not, $\alpha \leftarrow \alpha\beta$.\;
    }
\end{algorithm}

\begin{algorithm}[t] 
\caption{EiGLasso}
\label{alg:eiglasso}
    \DontPrintSemicolon
    \SetKwInOut{Input}{input}\SetKwInOut{Output}{output}
    \Input{{\bf Data} Sample covariances $\bm{S}=\frac{1}{nq}\sum_{i=1}^n\bm{Y}^{iT}\bm{Y}^i$ and $\bm{T}=\frac{1}{np}\sum_{i=1}^n\bm{Y}^{i}\bm{Y}^{iT}$\\
    {\bf Hyperparameters for} regularization $\gamma_{\bm{\Theta}}, \gamma_{\bm{\Psi}}$, Hessian approximation $K$, \\\quad line search $\sigma$, $0<\beta<1$, and trace ratio $\rho=\tr(\bm{\Psi})/\tr(\bm{\Theta})$}
    
    \Output{Parameters $\bm{\Theta}, \bm{\Psi}$}
    
    Initialize $\bm{\Theta} \leftarrow \bm{I}_p, \bm{\Psi} \leftarrow \bm{I}_q$.
    
    \For{$t=0,1,\dots$}{
        
       
        Eigendecompose $\bm{\Theta}$ and $\bm{\Psi}$.\;
        
        Compute $\bm{G}$ and $\bm{H}$.\;
        
        Determine active sets $\mathcal{A}_{\bm{\Theta}}$ and $\mathcal{A}_{\bm{\Psi}}$.\; 
        
        Compute $(D_{\bm{\Theta}},D_{\bm{\Psi}})$ via coordinate descent over the active sets.\;
        
        Use the Armijo rule to compute a step-size $\alpha$ as in Algorithm \ref{alg:linesearch}.
        
      
	Update $\bm{\Theta} \leftarrow \bm{\Theta} + \alpha D_{\bm{\Theta}}$, $\bm{\Psi} \leftarrow \bm{\Psi} + \alpha D_{\bm{\Psi}}$.
        
        Check convergence.\;
   }
   Adjust the diagonal elements of $\bm{\Theta}$ and $\bm{\Psi}$ according to the trace ratio $\rho$.\;
\end{algorithm}



\subsubsection{Efficient Computation of Gradient and Hessian via Eigendecomposition}

In Lemma \ref{lem:gh} below, we provide the form of the gradient $\bm{G}$ in Eq. \eqref{eq:compactg} and 
Hessian $\bm{H}$ in Eq. \eqref{eq:compacth}. 
We show that $\bm{G}$ can be represented in a significantly more compact form 
than what has been previously presented in \citet{biglasso} and \citet{TeraLasso}. 
We exploit this compact representation to further reduce the computation time  
via eigendecomposition of $\bm{\Theta}$ and $\bm{\Psi}$ in Theorem \ref{thm:gheigcomp}.

\begin{lemma}\label{lem:gh}
The gradient $\bm{G}$ in Eq. \eqref{eq:compactg} is given in terms of $\bm{W}=(\bm{\Theta}\oplus\bm{\Psi})^{-1}$
as 
\begin{align}\label{eq:partialmatrix0}
            \bm{G}_{\bm{\Theta}} = q\bm{S} - \bm{W}_{\bm{\Theta}}, 
            \quad
            \bm{G}_{\bm{\Psi}} = p\bm{T} - \bm{W}_{\bm{\Psi}}, 
\end{align}
where
\begin{align}\label{eq:partialmatrixW}
            \bm{W}_{\bm{\Theta}} =  \sum_{i=1}^q(\bm{I}_p\otimes\bm{e}_{q,i})^T\bm{W}(\bm{I}_p\otimes\bm{e}_{q,i}), \quad
            \bm{W}_{\bm{\Psi}} = \sum_{i=1}^p(\bm{e}_{p,i}\otimes\bm{I}_q)^T\bm{W}(\bm{e}_{p,i}\otimes\bm{I}_q).
\end{align}
The Hessian $\bm{H}$ in Eq. \eqref{eq:compacth} is given as
\begin{equation}\label{eq:HcollapseWW}
    \bm{H} = \begin{bmatrix}
            \bm{H}_{\bm{\Theta}}  & \bm{H}_{\bm{\Theta}\bm{\Psi}}\\
            \bm{H}_{\bm{\Theta}\bm{\Psi}}^T  & \bm{H}_{\bm{\Psi}}
        \end{bmatrix}
    = \begin{bmatrix}
    \bm{P}_{\bm{\Theta}}^T(\bm{W}\otimes\bm{W})\bm{P}_{\bm{\Theta}} & \bm{P}_{\bm{\Theta}}^T(\bm{W}\otimes\bm{W})\bm{P}_{\bm{\Psi}}  \\
    \bm{P}_{\bm{\Psi}}^T(\bm{W}\otimes\bm{W})\bm{P}_{\bm{\Theta}} & \bm{P}_{\bm{\Psi}}^T(\bm{W}\otimes\bm{W})\bm{P}_{\bm{\Psi}}
    \end{bmatrix}
    = \bm{P}^T(\bm{W}\otimes\bm{W})\bm{P},
\end{equation}
where $\bm{P}=\left[\bm{P}_{\bm{\Theta}},\bm{P}_{\bm{\Psi}}\right]$ is a $(p^2q^2)\times (p^2+q^2)$ matrix with
$\bm{P}_{\bm{\Theta}}=\sum_{i=1}^q \bm{I}_p\otimes\bm{e}_{q,i}\otimes\bm{I}_p\otimes\bm{e}_{q,i}$ and $ 
    \bm{P}_{\bm{\Psi}} = \sum_{j=1}^p \bm{e}_{p,j}\otimes\bm{I}_q\otimes\bm{e}_{p,j}\otimes\bm{I}_q$.
\end{lemma}
\begin{proof}
    Let $\bm{E}_{p, ij}=\bm{e}_{p,i}\bm{e}_{p,j}^T$ be a $p\times p$ one-hot matrix with $1$ at 
    the $(i,j)$th element and zero elsewhere. Then, for the gradient $\bm{G}$, we have
    \begin{align*}
         [\bm{G}_{\bm{\Theta}}]_{ij}  &= q[\bm{S}]_{ij} - \tr\left(\bm{W}(\bm{E}_{p, ij}\otimes\bm{I}_q)\right) \\
         &= q[\bm{S}]_{ij} - \tr\big(\bm{W}(\bm{E}_{p, ij}\otimes\sum_l\bm{e}_{q,l}\bm{e}_{q,l}^T)\big) \\
         &= q[\bm{S}]_{ij} - \tr\big(\sum_l\bm{W}(\bm{I}_p\otimes\bm{e}_{q,l})(\bm{E}_{p,ij}\otimes1)(\bm{I}_p\otimes\bm{e}_{q,l})^T\big) \\
         &= q[\bm{S}]_{ij} - \tr\big(\sum_l(\bm{I}_p\otimes\bm{e}_{q,l})^T\bm{W}(\bm{I}_p\otimes\bm{e}_{q,l})\bm{E}_{p,ij}\big).
    \end{align*}         
    By using $\tr(\bm{A}\bm{E}_{p,ij}) = \bm{A}_{ij}$, we collect the elements $[\bm{G}_{\bm{\Theta}}]_{ij}$ to obtain Eqs. \eqref{eq:partialmatrix0} and \eqref{eq:partialmatrixW}. The case for $\bm{G}_{\bm{\Psi}}$ can be shown similarly. Now, for the Hessian, we have
    \begin{align*}
        [\bm{H}_{\bm{\Theta}}&]_{jp+i,lp+k}
        = \tr\left(\bm{W}(\bm{E}_{p,ij}\otimes\bm{I}_q)\bm{W}(\bm{E}_{p,kl}\otimes\bm{I}_q)\right)\\
        &=\tr\left(\left[\sum_r  \bm{W}(\bm{I}_p\otimes\bm{e}_{q,r})\bm{E}_{p,ij}(\bm{I}_p\otimes\bm{e}_{q,r})^T \right] \left[\sum_s\bm{W}(\bm{I}_p\otimes\bm{e}_{q,s})\bm{E}_{p,kl}(\bm{I}_p\otimes\bm{e}_{q,s})^T\right]  \right)\\
        &= \tr\left( \sum_{r,s} (\bm{I}_p\otimes\bm{e}_{q,s})^T\bm{W}(\bm{I}_p\otimes\bm{e}_{q,r})\bm{E}_{p,ij}(\bm{I}_p\otimes\bm{e}_{q,r})^T\bm{W}(\bm{I}_p\otimes\bm{e}_{q,s})\bm{E}_{p,kl} \right),\\
        \intertext{by $\tr(\bm{A}\bm{B}\bm{A}\bm{C}) = \tvec(\bm{B})^T(\bm{A}\otimes\bm{A})\tvec(\bm{C})$ for symmetric matrices $\bm{A}$, $\bm{B}$, and $\bm{C}$, and $\tvec(\bm{E}_{p,ij}) = \bm{e}_{p^2,jp+i}$,}
        &= \bm{e}_{p^2,jp+i}^T \left( \sum_{r,s} (\bm{I}_p\otimes\bm{e}_{q,r})^T\bm{W}(\bm{I}_p\otimes\bm{e}_{q,s}) \otimes (\bm{I}_p\otimes\bm{e}_{q,r})^T\bm{W}(\bm{I}_p\otimes\bm{e}_{q,s}) \right) \bm{e}_{p^2,lp+k}\\
        &= \bm{e}_{p^2,jp+i}^T \left( \sum_{r,s} (\bm{I}_p\otimes\bm{e}_{q,r}\otimes\bm{I}_p\otimes\bm{e}_{q,r})^T(\bm{W}\otimes\bm{W})(\bm{I}_p\otimes\bm{e}_{q,s}\otimes\bm{I}_p\otimes\bm{e}_{q,s}) \right) \bm{e}_{p^2,lp+k}.
    \end{align*}
We collect the elements into a matrix 
\begin{align}
 \bm{H}_{\bm{\Theta}} =  \bigg(\sum_{r=1}^q\bm{I}_p\otimes\bm{e}_{q,r}\otimes\bm{I}_p\otimes\bm{e}_{q,r}\bigg)^T(\bm{W}\otimes\bm{W})\bigg(\sum_{s=1}^q \bm{I}_p\otimes\bm{e}_{q,s}\otimes\bm{I}_p\otimes\bm{e}_{q,s}\bigg). \nonumber
\end{align}    
The cases for $\bm{H}_{\bm{\Psi}}$ and $\bm{H}_{\bm{\Theta}\bm{\Psi}}$ can be shown in a similar way.
\begin{comment}
        \begin{equations}
        \begin{alignat}{2}
            &\bm{H}_{\bm{\Psi}} &&= \sum_{i=1}^q\sum_{j=1}^q (\bm{e}_i\otimes\bm{I}_q\otimes\bm{e}_i\otimes\bm{I}_q)^T(\bm{W}\otimes\bm{W})(\bm{e}_j\otimes\bm{I}_q\otimes\bm{e}_j\otimes\bm{I}_q)\label{eq:partialmatrix1d}\\
            &\bm{H}_{\bm{\Theta}\bm{\Psi}} &&= \sum_{i=1}^q\sum_{j=1}^p (\bm{I}_p\otimes\bm{e}_i\otimes\bm{I}_p\otimes\bm{e}_i)^T(\bm{W}\otimes\bm{W})(\bm{e}_j\otimes\bm{I}_q\otimes\bm{e}_j\otimes\bm{I}_q).\label{eq:partialmatrix1e}
        \end{alignat}
    \end{equations}
\end{comment}
\end{proof}




According to Lemma \ref{lem:gh}, the gradient $\bm{G}$ is obtained by collapsing
$\bm{W}$ as in Eq. \eqref{eq:partialmatrixW} and 
the Hessian $\bm{H}$ is obtained by collapsing $\bm{W} \otimes \bm{W}$ as in Eq. \eqref{eq:HcollapseWW}.
This collapse is illustrated in Figure \ref{fig:collapseW}(b) for the gradient. 
While the Kronecker sum inflates $\bm{\Theta}$ and $\bm{\Psi}$ into $\bm{\Omega}$ (Figure \ref{fig:collapseW}(a)), 
this inflated structure gets deflated, 
when computing $\bm{W}_{\bm{\Theta}}$ and $\bm{W}_{\bm{\Psi}}$ in the gradients from $\bm{W}=\bm{\Omega}^{-1}$ 
(Figure \ref{fig:collapseW}(b)). 
This deflation can be viewed as applying the mask $\mathbbm{1}_{p} \oplus \mathbbm{1}_{q}$,
where $\mathbbm{1}_{a}$ is a $a\times a$ matrix of one's, to $\bm{W}$ 
to obtain $\bm{W}^{\text{masked}}$ such that $[\bm{W}^{\text{masked}}]_{ij} =[\bm{W}]_{ij}$  
if $[\mathbbm{1}_{p} \oplus \mathbbm{1}_{q}]_{ij} \neq 0$, otherwise $[\bm{W}^{\text{masked}}]_{ij} = 0 $,
and replacing $\bm{W}$ in Eq. \eqref{eq:partialmatrixW} with $\bm{W}^{\text{masked}}$. 
In other words, only the non-zero elements of $\bm{W}^{\text{masked}}$ 
contribute towards $\bm{H}_{\bm{\Theta}}$ (blue in Figure \ref{fig:collapseW}(b))
and $\bm{H}_{\bm{\Psi}}$ (yellow in Figure \ref{fig:collapseW}(b)). 
The collapse of $\bm{W} \otimes \bm{W}$ in Eq. \eqref{eq:HcollapseWW} to obtain the Hessian
can be viewed as the same type of deflation applied twice to $\bm{W} \otimes \bm{W}$ (details in Appendix A).


Lemma \ref{lem:gh} reveals the challenge in computing $\bm{G}$ and $\bm{H}$ 
in a direct application of QUIC to the Kronecker-sum model: $\bm{G}$
involves computing the large $pq\times pq$ matrix $\bm{W}$ and $\bm{H}$ involves
computing the even larger $p^2q^2\times p^2q^2$ matrix $\bm{W} \otimes \bm{W}$.
In Theorem \ref{thm:gheigcomp} below, we show that $\bm{G}$ and $\bm{H}$ can be obtained
via eigendecomposition of $\bm{\Theta}$ and $\bm{\Psi}$ without explicitly constructing $\bm{W}$ and $\bm{W}\otimes \bm{W}$. 


\begin{theorem}\label{thm:gheigcomp}
Given the eigendecomposition  $\bm{\Theta}=\bm{Q}_{\bm{\Theta}}\bm{\Lambda}_{\bm{\Theta}}\bm{Q}_{\bm{\Theta}}^T$  
and  $\bm{\Psi}=\bm{Q}_{\bm{\Psi}}\bm{\Lambda}_{\bm{\Psi}}\bm{Q}_{\bm{\Psi}}^T$ 
with the $k$th smallest eigenvalue
$\lambda_{\bm{\Theta},k}$ and $\lambda_{\bm{\Psi},k}$ in the $(k,k)$th element of  $\bm{\Lambda}_{\bm{\Theta}}$ and $\bm{\Lambda}_{\bm{\Psi}}$,
the gradient in Eq. \eqref{eq:partialmatrix0} is given as
\begin{equation}\label{eq:partialmatrixeiggrad}
\begin{aligned}
     \bm{G}_{\bm{\Theta}} = q\bm{S} - \bm{Q}_{\bm{\Theta}}\left(\sum_{k=1}^q 
     \bm{\Xi}_{\bm{\Theta},k}
    
     \right)\bm{Q}_{\bm{\Theta}}^T, \quad\quad
     \bm{G}_{\bm{\Psi}} = p\bm{T} - \bm{Q}_{\bm{\Psi}}\left(\sum_{l=1}^p 
     \bm{\Xi}_{\bm{\Psi},l}
    
     \right)\bm{Q}_{\bm{\Psi}}^T,
\end{aligned}
\end{equation}
where
\begin{align}
    \bm{\Xi}_{\bm{\Theta},k} = \left(\bm{\Lambda}_{\bm{\Theta}} + \lambda_{\bm{\Psi},k}\bm{I}_p\right)^{-1}, \quad\quad \bm{\Xi}_{\bm{\Psi},l} = \left(\bm{\Lambda}_{\bm{\Psi}} + \lambda_{\bm{\Theta},l}\bm{I}_q\right)^{-1}. \nonumber
\end{align}
For the Hessian, we have
\begin{equation}
\!\bm{H} = \begin{bmatrix}
            \bm{H}_{\bm{\Theta}} & \!\!\! \bm{H}_{\bm{\Theta}\bm{\Psi}}\\
    \bm{H}_{\bm{\Theta}\bm{\Psi}}^T & \!\!\! \bm{H}_{\bm{\Psi}}
        \end{bmatrix}
    = \begin{bmatrix}
    \bm{Q}_{\bm{\Theta}}\otimes\bm{Q}_{\bm{\Theta}}&\!\!\!\!\bm{0}\\
    \bm{0}& \!\!\!\!\bm{Q}_{\bm{\Psi}}\otimes\bm{Q}_{\bm{\Psi}}
    \end{bmatrix}
    \begin{bmatrix}
    \bm{\Lambda}_{\bm{H}_{\bm{\Theta}}} & \!\!\!\!\! \bm{\Upsilon} \\
    \bm{\Upsilon}^T & \!\!\!\!\! \bm{\Lambda}_{\bm{H}_{\bm{\Psi}}}
    \end{bmatrix}
    \begin{bmatrix}
    \bm{Q}_{\bm{\Theta}}\otimes\bm{Q}_{\bm{\Theta}}&\!\!\!\!\bm{0}\\
    \bm{0}&\!\!\!\!\bm{Q}_{\bm{\Psi}}\otimes\bm{Q}_{\bm{\Psi}}
    \end{bmatrix}^T\!\!\!,
\label{eq:partialmatrixeigH}
\end{equation}
where 
\begin{align}
    \bm{\Lambda}_{\bm{H}_{\bm{\Theta}}} &= \sum_{k=1}^q \bm{\Xi}_{\bm{\Theta},k} \otimes \bm{\Xi}_{\bm{\Theta},k}, \quad
     \bm{\Lambda}_{\bm{H}_{\bm{\Psi}}} = \sum_{l=1}^p \bm{\Xi}_{\bm{\Psi},l} \otimes \bm{\Xi}_{\bm{\Psi},l}, \label{eq:Hcomp} \\
     \bm{\Upsilon} &= \tilde{\bm{I}}_{p^2,p}
        \begin{bmatrix}
        \lambda^2_{\bm{W},11} & \cdots & \lambda^2_{\bm{W},1q} \\
        \vdots & \ddots & \vdots \\
        \lambda^2_{\bm{W},p1} & \cdots & \lambda^2_{\bm{W},pq}
        \end{bmatrix}
        \tilde{\bm{I}}_{q,q^2}. \nonumber
\end{align}
$\tilde{\bm{I}}_{p^2,p}$ and $\tilde{\bm{I}}_{q,q^2}$ above are stretched 
identity matrices, $\tilde{\bm{I}}_{p^2,p} =
        \begin{bmatrix}
        \bm{E}_{p,11}\;
        \bm{E}_{p,22}\;
        \cdots\;
        \bm{E}_{p,pp}
        \end{bmatrix}^T$ and $\tilde{\bm{I}}_{q,q^2} =
        \begin{bmatrix}
        \bm{E}_{q,11}\;
        \bm{E}_{q,22}\;
        \cdots\;
        \bm{E}_{q,qq}
        \end{bmatrix}$, 
where $\bm{E}_{p, ij}$ is a $p\times p$ one-hot matrix with $1$ at the $(i,j)$th element and zero elsewhere, 
and $\lambda_{\bm{W},lk} = (\lambda_{\bm{\Theta},l} + \lambda_{\bm{\Psi},k})^{-1}$.
\end{theorem}
\begin{proof}
We prove for $\bm{G}_{\bm{\Theta}}$, $\bm{H}_{\bm{\Theta}}$, and $\bm{H}_{\bm{\Theta}\bm{\Psi}}$, 
because the case for $\bm{G}_{\bm{\Psi}}$ and $\bm{H}_{\bm{\Psi}}$ can be proved similarly. 
Let $\bm{q}_{\bm{\Theta},i}$ and $\bm{q}_{\bm{\Psi},i}$ be the $i$th eigenvectors of $\bm{\Theta}$ and $\bm{\Psi}$, given as the $i$th columns of $\bm{Q}_{\bm{\Theta}}$ and $\bm{Q}_{\bm{\Psi}}$. Then, we can re-write Eq. \eqref{eq:Weigdecomp} as
\begin{equation}\label{eq:Weigdecompsum}
    \bm{W} = \sum_{l=1}^p \sum_{k=1}^q \lambda_{\bm{W},lk}(\bm{q}_{\bm{\Theta},l}\otimes\bm{q}_{\bm{\Psi},k})(\bm{q}_{\bm{\Theta},l}\otimes\bm{q}_{\bm{\Psi},k})^T.
\end{equation}

For $\bm{G}_{\bm{\Theta}}$, we substitute $\bm{W}$ in $\bm{W}_{\bm{\Theta}}$ in Eq. \eqref{eq:partialmatrixW} with Eq. \eqref{eq:Weigdecompsum}:
\begin{align*}
    \bm{G}_{\bm{\Theta}}
    &= q\bm{S} - \sum_{i=1}^q(\bm{I}_p\otimes\bm{e}_{q,i})^T\left(\sum_{l=1}^p \sum_{k=1}^q \lambda_{\bm{W},lk}(\bm{q}_{\bm{\Theta},l}\otimes\bm{q}_{\bm{\Psi},k})(\bm{q}_{\bm{\Theta},l}\otimes\bm{q}_{\bm{\Psi},k})^T\right)(\bm{I}_p\otimes\bm{e}_{q,i})\\
    &= q\bm{S} - \sum_{i=1}^q\sum_{l=1}^p \sum_{k=1}^q\lambda_{\bm{W},lk}([\bm{q}_{\bm{\Psi},k}]_i\bm{q}_{\bm{\Theta},l})([\bm{q}_{\bm{\Psi},k}]_i\bm{q}_{\bm{\Theta},l})^T \quad\text{since $(\bm{I}\otimes\bm{e}_{q,i})^T(\bm{a}\otimes\bm{b})=\bm{b}_i\bm{a}$}\\
    &= q\bm{S} - \sum_{l=1}^p \sum_{k=1}^q \lambda_{\bm{W},lk}\bm{q}_{\bm{\Theta},l}\bm{q}_{\bm{\Theta},l}^T
    \;\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\text{since $\bm{Q}_{\bm{\Psi}}$ is orthonormal}\\
    &= q\bm{S} - \sum_{k=1}^q \bm{Q}_{\bm{\Theta}}(\bm{\Lambda}_{\bm{\Theta}} + \lambda_{\bm{\Psi},k}\bm{I}_p)^{-1}\bm{Q}_{\bm{\Theta}}^T  .
\end{align*}
This is equivalent to $\bm{G}_{\bm{\Theta}}$ in Eq. \eqref{eq:partialmatrixeiggrad}, 
since the summation can be performed over $\bm{\Xi}_{\bm{\Theta},k}$, 
instead of over $\bm{Q}_{\bm{\Theta}}\bm{\Xi}_{\bm{\Theta},k}\bm{Q}_{\bm{\Theta}}^T$. 

To show $\bm{H}_{\bm{\Theta}}$, we again substitute $\bm{W}$ in $\bm{H}_{\bm{\Theta}}$ of Eq. \eqref{eq:HcollapseWW}
with Eq. \eqref{eq:Weigdecompsum}.
\begin{align*}
    \bm{H}_{\bm{\Theta}} &=  \sum_{i=1}^q \sum_{j=1}^q \bigg(\sum_{l=1}^p \sum_{k=1}^q\lambda_{\bm{W},lk}([\bm{q}_{\bm{\Psi},k}]_i\bm{q}_{\bm{\Theta},l})([\bm{q}_{\bm{\Psi},k}]_j\bm{q}_{\bm{\Theta},l})^T
    \bigg. \\
	& \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \bigg.\otimes\sum_{r=1}^p \sum_{s=1}^q\lambda_{\bm{W},rs}([\bm{q}_{\bm{\Psi},s}]_i\bm{q}_{\bm{\Theta},r})([\bm{q}_{\bm{\Psi},s}]_j\bm{q}_{\bm{\Theta},r})^T \bigg)\\
    &=\sum_{k=1}^q\bigg(\sum_{l=1}^p \lambda_{\bm{W},lk}\bm{q}_{\bm{\Theta},l}\bm{q}_{\bm{\Theta},l}^T\otimes\sum_{r=1}^p\lambda_{\bm{W},rk}\bm{q}_{\bm{\Theta},r}\bm{q}_{\bm{\Theta},r}^T\bigg) \\
    &= \sum_{k=1}^q \bigg(\bm{Q}_{\bm{\Theta}}(\bm{\Lambda}_{\bm{\Theta}} + \lambda_{\bm{\Psi},k}\bm{I}_p)^{-1}\bm{Q}_{\bm{\Theta}}^T\bigg) \otimes \bigg(\bm{Q}_{\bm{\Theta}}(\bm{\Lambda}_{\bm{\Theta}} + \lambda_{\bm{\Psi},k}\bm{I}_p)^{-1}\bm{Q}_{\bm{\Theta}}^T\bigg)\\
    &= (\bm{Q}_{\bm{\Theta}}\otimes\bm{Q}_{\bm{\Theta}}) \bigg(\sum_{k=1}^q (\bm{\Lambda}_{\bm{\Theta}} + \lambda_{\bm{\Psi},k}\bm{I}_p)^{-1}\otimes (\bm{\Lambda}_{\bm{\Theta}} + \lambda_{\bm{\Psi},k}\bm{I}_p)^{-1}\bigg)(\bm{Q}_{\bm{\Theta}}\otimes\bm{Q}_{\bm{\Theta}})^T.
\end{align*}
For $\bm{H}_{\bm{\Theta}\bm{\Psi}}$, we substitute $\bm{W}$ in $\bm{H}_{\bm{\Theta}\bm{\Psi}}$ of Eq. \eqref{eq:HcollapseWW}
with Eq. \eqref{eq:Weigdecompsum}.
\begin{align*}
        \bm{H}_{\bm{\Theta}\bm{\Psi}} &= \sum_{i=1}^q \sum_{j=1}^p \bigg(\sum_{l=1}^p \sum_{k=1}^q\lambda_{\bm{W},lk}([\bm{q}_{\bm{\Psi},k}]_i\bm{q}_{\bm{\Theta},l})([\bm{q}_{\bm{\Theta},l}]_j\bm{q}_{\bm{\Psi},k})^T\bigg. \\
	    &\bigg. \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \otimes\sum_{r=1}^p \sum_{s=1}^q\lambda_{\bm{W},rs}([\bm{q}_{\bm{\Psi},s}]_i\bm{q}_{\bm{\Theta},r})([\bm{q}_{\bm{\Theta},r}]_j\bm{q}_{\bm{\Psi},s})^T \bigg)\\
        &=\sum_{l=1}^p \sum_{k=1}^q \lambda_{\bm{W},lk}^2(\bm{q}_{\bm{\Theta},l}\bm{q}_{\bm{\Psi},k}^T \otimes \bm{q}_{\bm{\Theta},l}\bm{q}_{\bm{\Psi},k}^T) \\
        &= \sum_{l=1}^p (\bm{q}_{\bm{\Theta},l}\otimes\bm{q}_{\bm{\Theta},l})\sum_{k=1}^q\lambda_{\bm{W},lk}^2(\bm{q}_{\bm{\Psi},k}\otimes\bm{q}_{\bm{\Psi},k})^T\\
        &= \begin{bmatrix}
        \bm{q}^T_{\bm{\Theta},1}\otimes\bm{q}^T_{\bm{\Theta},1}\\
        \bm{q}^T_{\bm{\Theta},2}\otimes\bm{q}^T_{\bm{\Theta},2}\\
        \vdots\\
        \bm{q}^T_{\bm{\Theta},p}\otimes\bm{q}^T_{\bm{\Theta},p}
        \end{bmatrix}^T
        \begin{bmatrix}
        \lambda^2_{\bm{W},11} & \cdots & \lambda^2_{\bm{W},1q} \\
        \vdots & \ddots & \vdots \\
        \lambda^2_{\bm{W},p1} & \cdots & \lambda^2_{\bm{W},pq}
        \end{bmatrix}
        \begin{bmatrix}
        \bm{q}^T_{\bm{\Psi},1}\otimes\bm{q}^T_{\bm{\Psi},1}\\
        \bm{q}^T_{\bm{\Psi},2}\otimes\bm{q}^T_{\bm{\Psi},2}\\
        \vdots\\
        \bm{q}^T_{\bm{\Psi},q}\otimes\bm{q}^T_{\bm{\Psi},q}
        \end{bmatrix}\\
        &=(\bm{Q}_{\bm{\Theta}}\otimes\bm{Q}_{\bm{\Theta}})\tilde{\bm{I}}_{p^2,p}
        \begin{bmatrix}
        \lambda^2_{\bm{W},11} & \cdots & \lambda^2_{\bm{W},1q} \\
        \vdots & \ddots & \vdots \\
        \lambda^2_{\bm{W},p1} & \cdots & \lambda^2_{\bm{W},pq}
        \end{bmatrix}
        \tilde{\bm{I}}_{q,q^2}(\bm{Q}_{\bm{\Psi}}\otimes\bm{Q}_{\bm{\Psi}})^T.
    \end{align*}
\end{proof}

Theorem \ref{thm:gheigcomp} provides a significantly more efficient method for computing $\bm{G}$ and $\bm{H}$, 
compared to Lemma \ref{lem:gh}. Since $\bm{G}$ and $\bm{H}$ can be written entirely in terms of 
the eigenvectors and eigenvalues of $\bm{\Theta}$ and $\bm{\Psi}$, the operations that involve
$\bm{W}=(\bm{\Theta} \oplus \bm{\Psi})^{-1}$
is replaced with the cheaper operation of eigendecomposition with time $\mathcal{O}(p^3+q^3)$ and space $\mathcal{O}(p^2 + q^2)$. 
TeraLasso also used the eigendecomposition of $\bm{\Theta}$ and $\bm{\Psi}$ in each iteration of the optimization. However, TeraLasso has done this 
for an efficient projection of the gradient to the re-parameterized space, whereas 
EiGLasso performs the eigendecomposition for an efficient computation of $\bm{G}$ and $\bm{H}$ in the original parameter space. 






It follows from Theorem \ref{thm:gheigcomp} that $(\bm{\Theta},\bm{\Psi})$'s in the same equivalence class 
in Eq. \eqref{eqn:noniddiag} differ only in their eigenvalues, but not in eigenvectors. 
The equivalence class in Eq. \eqref{eqn:noniddiag} can be written equivalently as
 \begin{align}\label{eqn:noniddiagEig}
  S[\bm{\Omega}] = \{ (\bm{\Theta},\bm{\Psi})\in \mathbb{KS}^{p,q}| 
	\left(\bm{\Lambda}_{\bm{\Theta}}-c\bm{I}_p\right)\oplus\left(\bm{\Lambda}_{\bm{\Psi}}+c\bm{I}_q\right)=\bm{\Lambda}_{\bm{\Omega}},
		c\in \mathbb{R}\}, 
 \end{align}
since $\bm{\Theta} + c\bm{I}_p = \bm{Q}_{\bm{\Theta}}\bm{\Lambda}_{\bm{\Theta}}\bm{Q}_{\bm{\Theta}}^T +c\bm{Q}_{\bm{\Theta}}\bm{Q}_{\bm{\Theta}}^T 
= \bm{Q}_{\bm{\Theta}}(\bm{\Lambda}_{\bm{\Theta}}+c\bm{I}_p)\bm{Q}_{\bm{\Theta}}^T$ and similarly for $\bm{\Psi}$.
We exploit this feature 
to handle the unidentifiable parameters during estimation in Section \ref{sec:identify} and 
to analyze the convergence in Section \ref{sec:convg}. 





In the theorem below, we prove several properties of $\bm{H}$ that we will use to prove
the line-search properties, convergence guarantees, and convergence rates 
in Section \ref{sec:convg}.

\begin{theorem}\label{thm:Hpd}
The Hessian $\bm{H}$ in EiGLasso has the following properties:
\begin{itemize}
    \item It is positive semi-definite with the null space 
$\text{null}(\bm{H})=\{[\tvec(\bm{X}_p)^T \tvec(\bm{X}_q)^T]^T| \bm{X}_p\oplus\bm{X}_q = 0,
\bm{X}_p \in \mathbb{R}^{p\times p}, \bm{X}_p \in \mathbb{R}^{q\times q}\}$.
    
    \item On the Kronecker-sum space in Eq. \eqref{def:ksspace},
$\bm{H}$ is positive definite.
The minimum eigenvalue $\lambda_{\bm{H},\min_0}$ outside of the nullspace and the maximum eigenvalue 
$\lambda_{\bm{H},\max}$ are given in terms of the eigenvalues of $\bm{\Theta}$ and $\bm{\Psi}$:
\begin{equation}
\lambda_{\bm{H},\min_0}=\min\{p,q\}(\lambda_{\bm{\Theta},p}+\lambda_{\bm{\Psi},q})^{-2}, \quad
\lambda_{\bm{H},\max}=(p+q)(\lambda_{\bm{\Theta},1}+\lambda_{\bm{\Psi},1})^{-2}. \nonumber
\end{equation}    
\end{itemize}
\end{theorem}
\begin{proof}
To prove the first property, from Lemma \ref{lem:gh}, we notice  that since $\bm{W}$ is positive definite, 
the null space of $\bm{H}$ is given by $\bm{X}_p$ and $\bm{X}_q$ that satisfy
\begin{align*}
    \bm{P}\begin{bmatrix}
        \tvec(\bm{X}_p)\\
        \tvec(\bm{X}_q)
    \end{bmatrix}
    &=0.
\end{align*} 
Since the left-hand side of the above can be written as
\begin{align*}
  \bm{P}\begin{bmatrix}
        \tvec(\bm{X}_p)\\
        \tvec(\bm{X}_q)
    \end{bmatrix}
    &=\sum_{i=1}^q \tvec\left((\bm{I}_p\otimes\bm{e}_{q,i})\bm{X}_p(\bm{I}_p\otimes\bm{e}_{q,i})^T \right) + \sum_{j=1}^p \tvec\left((\bm{e}_{p,j}\otimes\bm{I}_q)\bm{X}_q(\bm{e}_{p,j}\otimes\bm{I}_q)^T \right) \\
    &=\tvec\left( \bm{X}_p\otimes\bm{I}_q\right) +
    \tvec\left( \bm{I}_p\otimes\bm{X}_q\right) \\
    &= \tvec\left(\bm{X}_p\oplus\bm{X}_q\right),
\end{align*}
the null space of $\bm{H}$ is 
given by $\bm{X}_p$ and $\bm{X}_q$ that satisfy $\bm{X}_p\oplus\bm{X}_q = 0$. 


To prove the second property, we notice that the null space of $\bm{H}$ is outside of
the constraint region $\bm{\Theta}\oplus\bm{\Psi}\succ0$;  
thus, $\bm{H}$ is positive definite on the Kronecker-sum space in Eq. \eqref{def:ksspace}.
To find the eigenvalues of $\bm{H}$, we first find the eigenvalues of 
\begin{align*}
    \bm{P}^T\bm{P} =
    \begin{bmatrix}
        \bm{P}_{\bm{\Theta}}^T\bm{P}_{\bm{\Theta}} & \bm{P}_{\bm{\Theta}}^T\bm{P}_{\bm{\Psi}}\\
        \bm{P}_{\bm{\Psi}}^T\bm{P}_{\bm{\Theta}} & \bm{P}_{\bm{\Psi}}^T\bm{P}_{\bm{\Psi}}
    \end{bmatrix}
    =
    \begin{bmatrix}
        q\bm{I}_{p^2} & \tvec(\bm{I}_p)\tvec(\bm{I}_q)^T\\
        \tvec(\bm{I}_q)\tvec(\bm{I}_p)^T & p\bm{I}_{q^2}
    \end{bmatrix},
\end{align*}
by finding the solutions of the characteristic equation $\det(\bm{P}^T\bm{P} - \lambda\bm{I})=0$:
\begin{align*}
    \det(&\bm{P}^T\bm{P} - \lambda\bm{I}_{p^2+q^2})  
    = \det(q\bm{I}_{p^2}-\lambda\bm{I}_{p^2})\det\left(p\bm{I}_{q^2} 
    - \lambda\bm{I}_{q^2} - \frac{p}{q-\lambda}\tvec(\bm{I}_q)\tvec(\bm{I}_q)^T \right)\\
    &=  (q-\lambda)^{p^2} 
    \det\left((p-\lambda)\bm{I}_{q^2} - \frac{p}{q-\lambda}\tvec(\bm{I}_q)\tvec(\bm{I}_q)^T \right)\\
    &=  (q-\lambda)^{p^2-1} (p-\lambda)^{q^2-1}
    \big((p-\lambda)(q-\lambda)-pq\big) \quad\quad \text{from matrix determinant lemma}\\
    &= (p-\lambda)^{q^2-1}(q-\lambda)^{p^2-1}\lambda\big(\lambda-(p+q)\big) = 0.
\end{align*}
Thus, the eigenvalues of $\bm{P}^T\bm{P}$ are $0$, $p$, $q$, and $(p+q)$, each with the algebraic multiplicity
$1$, $q^2-1$, $p^2-1$, and 1. 
From Eq. \eqref{eq:HcollapseWW}, for any unit vector $\bm{u}\in\mathbb{R}^{p^2+q^2}\setminus\text{null}(\bm{H})$,
\begin{align*}
    \lambda_{\bm{W},\min}^2(\bm{u}^T\bm{P}^T\bm{P}\bm{u})
    \le
        \bm{u}^T
    \bm{P}^T(\bm{W}\otimes\bm{W})\bm{P} \bm{u}
    \le \lambda_{\bm{W},\max}^2(\bm{u}^T\bm{P}^T\bm{P}\bm{u}),
\end{align*}
where $\lambda_{\bm{W},\min}= (\lambda_{\bm{\Theta},p}+\lambda_{\bm{\Psi},q})^{-2}$ and 
$\lambda_{\bm{W},\max}= (\lambda_{\bm{\Theta},1}+\lambda_{\bm{\Psi},1})^{-2}$ are the minimum and maximum eigenvalues of $\bm{W}$.
Combining this with the minimum eigenvalue $\min\{p,q\}$ 
and the maximum eigenvalue $(p+q)$ of $\bm{P}^T\bm{P}$ outside of the null space of $\bm{H}$ results in
\begin{equation*}
    \min\{p,q\}\lambda_{\bm{W},\min}^2\le\lambda_{\bm{W},\min}^2(\bm{u}^T\bm{P}^T\bm{P}\bm{u})\quad\text{and}
    \quad \lambda_{\bm{W},\max}^2(\bm{u}^T\bm{P}^T\bm{P}\bm{u})\le (p+q)\lambda_{\bm{W},\max}^2,
\end{equation*}
which completes the proof. 
\end{proof}









\subsection{EiGLasso with Approximate Hessian}\label{sec:approxH}

While EiGLasso pre-computes and stores $\bm{G}$ after eigendecomposing  $\bm{\Theta}$ and $\bm{\Psi}$
at the beginning of the coordinate descent to solve Eq. \eqref{eq:newton}, 
explicitly computing and storing the $(p^2+q^2)\times(p^2+q^2)$ 
matrix $\bm{H}$ is still expensive for large $p$ and $q$ even with the eigendecomposition. 
To further reduce the computation time and memory, we approximate $\bm{H}$ 
based on its form in Theorem \ref{thm:gheigcomp} as follows.
We first drop $\bm{H}_{\bm{\Theta}\bm{\Psi}}$ while keeping only $\bm{H}_{\bm{\Theta}}$ and $\bm{H}_{\bm{\Psi}}$.
Since pre-computing and storing the $p^2\times p^2$ matrix $\bm{H}_{\bm{\Theta}}$ and 
the $q^2\times q^2$ matrix $\bm{H}_{\bm{\Psi}}$ is still expensive,
we approximate $\bm{H}_{\bm{\Theta}}$ and $\bm{H}_{\bm{\Psi}}$ by approximating their eigenvalues $\bm{\Lambda}_{\bm{H}_{\bm{\Theta}}}$ 
and $\bm{\Lambda}_{\bm{H}_{\bm{\Psi}}}$ in Eq. \eqref{eq:partialmatrixeigH}, while keeping the same eigenvectors
$\bm{Q}_{\bm{\Theta}}\otimes\bm{Q}_{\bm{\Theta}}$ and $\bm{Q}_{\bm{\Psi}}\otimes\bm{Q}_{\bm{\Psi}}$. 
The resulting approximate Hessian $\hat{\bm{H}}$ is 
\begin{equation}\label{eq:approxH}
    \hat{\bm{H}} = \begin{bmatrix}
            \hat{\bm{H}}_{\bm{\Theta}} & \!\!\!\bm{0} \\
    \bm{0} & \!\!\! \hat{\bm{H}}_{\bm{\Psi}}
        \end{bmatrix}\\
    = \begin{bmatrix}
    \bm{Q}_{\bm{\Theta}}\otimes\bm{Q}_{\bm{\Theta}}&\!\!\!\bm{0}\\
    \bm{0}&\!\!\!\bm{Q}_{\bm{\Psi}}\otimes\bm{Q}_{\bm{\Psi}}
    \end{bmatrix}
    \begin{bmatrix}
    \hat{\bm{\Lambda}}_{\bm{H}_{\bm{\Theta}}} & \!\!\bm{0} \\
    \bm{0} & \!\!\hat{\bm{\Lambda}}_{\bm{H}_{\bm{\Psi}}}
    \end{bmatrix}
    \begin{bmatrix}
    \bm{Q}_{\bm{\Theta}}\otimes\bm{Q}_{\bm{\Theta}}& \!\!\!\bm{0}\\
    \bm{0}&\!\!\!\bm{Q}_{\bm{\Psi}}\otimes\bm{Q}_{\bm{\Psi}}
    \end{bmatrix}^T\!\!\!,
\end{equation}
where
\begin{equation}\label{eqn:approxHL}
\begin{alignedat}{2}
    \hat{\bm{\Lambda}}_{\bm{H}_{\bm{\Theta}}} &= \sum_{k=1}^{K} \bm{\Xi}_{\bm{\Theta},k}\otimes\bm{\Xi}_{\bm{\Theta},k} 
        + (q-K) \bm{\Xi}_{\bm{\Theta},K} \otimes \bm{\Xi}_{\bm{\Theta},K},\\
    \hat{\bm{\Lambda}}_{\bm{H}_{\bm{\Psi}}} &= \sum_{l=1}^{K} \bm{\Xi}_{\bm{\Psi},l}\otimes\bm{\Xi}_{\bm{\Psi},l} 
        + (p-K) \bm{\Xi}_{\bm{\Psi},K} \otimes \bm{\Xi}_{\bm{\Psi},K}.
\end{alignedat}
\end{equation}
In $\hat{\bm{\Lambda}}_{\bm{H}_{\bm{\Theta}}}$, we keep only the components $\bm{\Xi}_{\bm{\Theta},k}\otimes \bm{\Xi}_{\bm{\Theta},k}$'s 
of $\bm{\Lambda}_{\bm{H}_{\bm{\Theta}}}$ in Eq. \eqref{eq:Hcomp} 
that contain the $K$ smallest eigenvalues $\lambda_{\bm{\Psi},1}, \ldots, \lambda_{\bm{\Psi},K}$. 
These smallest eigenvalues have the largest contribution toward $\hat{\bm{H}}_{\bm{\Theta}}$, 
since the eigenvalues of $\bm{\Psi}$ contribute to the eigenvalues of $\bm{H}_{\bm{\Theta}}$ through their inverses. 
We replace the remaining  $(q-K)$ components with $\bm{\Xi}_{\bm{\Theta},K}\otimes \bm{\Xi}_{\bm{\Theta},K}$, 
because dropping the $(q-K)$ eigenvalue components in ${\bm{\Lambda}}_{\bm{H}_{\bm{\Theta}}}$ would amount to 
assuming that $\bm{\Psi}$ has $(q-K)$ eigenvalues of infinite magnitude. 
This approximation leads to the following form 
\begin{align*}
    	\hat{\bm{H}}_{\bm{\Theta}} &= \sum_{k=1}^K \bm{V}_{\bm{\Theta}, k}\otimes\bm{V}_{\bm{\Theta}, k} 
		+ (q-K) \bm{V}_{\bm{\Theta}, K+1} \otimes \bm{V}_{\bm{\Theta}, K+1}, \\
    	\hat{\bm{H}}_{\bm{\Psi}} &= \sum_{k=1}^K \bm{V}_{\bm{\Psi}, k}\otimes\bm{V}_{\bm{\Psi}, k} 
		+ (q-K) \bm{V}_{\bm{\Psi}, K+1} \otimes \bm{V}_{\bm{\Psi}, K+1},
\end{align*}
where $\bm{V}_{\bm{\Theta},k} = \bm{Q}_{\bm{\Theta}}\bm{\Xi}_{\bm{\Theta},k}\bm{Q}_{\bm{\Theta}}^T$
and $\bm{V}_{\bm{\Psi},k} = \bm{Q}_{\bm{\Psi}}\bm{\Xi}_{\bm{\Psi},k}\bm{Q}_{\bm{\Psi}}^T$.
Then, during the coordinate descent optimization, we pre-compute and store $\bm{V}_{\bm{\Theta},k}$'s 
and $\bm{V}_{\bm{\Psi},k}$'s for $k=1,\ldots,K$, 
and compute the Hessian entries from $\bm{V}_k$'s as needed in each coordinate descent update (details in Appendix B).
We assume the same  $K$ for both $\hat{\bm{H}}_{\bm{\Theta}}$ and $\hat{\bm{H}}_{\bm{\Psi}}$, though this can be relaxed.
In our experiments in Section \ref{sec:exp} we demonstrate that often $K=1$ suffices.



We provide a geometric interpretation of our approximation of $\bm{H}_{\bm{\Theta}}$ with $\hat{\bm{H}}_{\bm{\Theta}}$.
Assuming Eq. \eqref{eq:newton} without the $L_1$-regularization,
we have the  descent direction $-\bm{H}_{\bm{\Theta}}^{-1}\tvec(\bm{G}_{\bm{\Theta}})$ with the exact Hessian
and $-\hat{\bm{H}}^{-1}_{\bm{\Theta}}\tvec(\bm{G}_{\bm{\Theta}})$ with the approximate Hessian. 
Since $\bm{H}_{\bm{\Theta}}$ and $\hat{\bm{H}}_{\bm{\Theta}}$ share the same eigenvectors,
in the coordinate system defined  by these eigenvectors as bases, 
the descent directions become $-\bm{H}_{\bm{\Theta}}^{-1/2} \tvec(\bm{G}_{\bm{\Theta}}) =
-{\bm{\Lambda}_{\bm{H}_{\bm{\Theta}}}^{-\frac{1}{2}}} (\bm{Q}_{\bm{\Theta}}\otimes\bm{Q}_{\bm{\Theta}})^T
    \tvec(\bm{G}_{\bm{\Theta}})$ with the exact Hessian and
$-\hat{\bm{H}}_{\bm{\Theta}}^{-\frac{1}{2}}\tvec(\bm{G}_{\bm{\Theta}}) 
    =  -\hat{\bm{\Lambda}}_{\bm{H}_{\bm{\Theta}}}^{-\frac{1}{2}}
    (\bm{Q}_{\bm{\Theta}}\otimes\bm{Q}_{\bm{\Theta}})^T \tvec(\bm{G}_{\bm{\Theta}})$ with the approximate Hessian
    \citep{cvxbook}. The latter is an element-wise scaled matrix of the former.
Furthermore, since $[\bm{\Lambda}_{\bm{H}_{\bm{\Theta}}}^{-\frac{1}{2}}]_{ii} <= [\hat{\bm{\Lambda}}_{\bm{H}_{\bm{\Theta}}}^{-\frac{1}{2}}]_{ii}$, 
in this new coordinate system,
each element of the descent direction with $\hat{\bm{H}}_{\bm{\Theta}}$ is a convex combination of $0$ and 
the corresponding element of the descent direction with $\bm{H}_{\bm{\Theta}}$.







The following theorem provides the properties of $\hat{\bm{H}}$, which will be used to analyze
the convergence of EiGLasso in Section \ref{sec:convg}.

\begin{theorem}\label{thm:Happroxpd}
The approximate Hessian $\hat{\bm{H}}$ is positive definite. Furthermore, its minimum and maximum eigenvalues are 
\begin{equation} 
\lambda_{\hat{\bm{H}},\min}=\min\big\{\lambda_{\hat{\bm{H}}_{\bm{\Theta}},\min}, \lambda_{\hat{\bm{H}}_{\bm{\Psi}},\min}\big\}, \quad
\lambda_{\hat{\bm{H}},\max}=\max\big\{\lambda_{\hat{\bm{H}}_{\bm{\Theta}},\max},\lambda_{\hat{\bm{H}}_{\bm{\Psi}},\max}\big\},
\nonumber
\end{equation}
where
\begin{align*}
    \lambda_{\hat{\bm{H}}_{\bm{\Theta}},\min} &= \sum_{i=1}^K (\lambda_{\bm{\Theta}, p} + \lambda_{\bm{\Psi},i})^{-2} +  (q-K)(\lambda_{\bm{\Theta}, p} + \lambda_{\bm{\Psi},K})^{-2}\\
    \lambda_{\hat{\bm{H}}_{\bm{\Theta}},\max} &= \sum_{i=1}^K (\lambda_{\bm{\Theta}, 1} + \lambda_{\bm{\Psi},i})^{-2} +  (q-K)(\lambda_{\bm{\Theta}, 1} + \lambda_{\bm{\Psi},K})^{-2}\\
    \lambda_{\hat{\bm{H}}_{\bm{\Psi}},\min} &= \sum_{j=1}^{K} (\lambda_{\bm{\Theta}, j} + \lambda_{\bm{\Psi},q})^{-2} +  (p-K)(\lambda_{\bm{\Theta}, K} + \lambda_{\bm{\Psi},q})^{-2}\\
    \lambda_{\hat{\bm{H}}_{\bm{\Psi}},\max} &= \sum_{j=1}^{K} (\lambda_{\bm{\Theta}, j} + \lambda_{\bm{\Psi},1})^{-2} +  (p-K)(\lambda_{\bm{\Theta}, K} + \lambda_{\bm{\Psi},1})^{-2}.
\end{align*}
\end{theorem}
\begin{proof}
Since a block-diagonal matrix has the same set of eigenvalues as those of its diagonal blocks, 
the eigenvalues of $\hat{\bm{H}}$ are the same as a union of those of $\hat{\bm{H}}_{\bm{\Theta}}$ 
and $\hat{\bm{H}}_{\bm{\Psi}}$ in Eq. \eqref{eqn:approxHL}.
The eigenvalues of $\hat{\bm{H}}_{\bm{\Theta}}$ and $\hat{\bm{H}}_{\bm{\Psi}}$ are positive, because they
are also the eigenvalues of $\bm{W}$ 
in Eq. \eqref{eq:Weigdecomp}, the inverses of which are the eigenvalues of the positive definite 
matrix $\bm{\Omega}$ in Eq. \eqref{eq:omegaeigdecomp}. Thus, $\hat{\bm{H}}$ is positive definite.
\end{proof}




\subsection{Estimation with the Unidentifiable Parameters} \label{sec:identify}



We describe a simple approach to estimating the diagonals of $\bm{\Theta}$ and $\bm{\Psi}$ with EiGLasso.
We show that to solve the optimization problem with the constraint on a fixed trace ratio 
in Eq. \eqref{eq:objrho}, 
it is sufficient to solve the unconstrained problem in Eq. \eqref{eq:obj} and 
to adjust the diagonals once to enforce the given trace ratio after EiGLasso converges. 

In Lemma \ref{lem:invariantGH} below, we show that given the current estimate of $\bm{\Theta}$ and $\bm{\Psi}$, 
the quadratic approximation in Eq. \eqref{eq:newton} for determining descent directions 
is uniquely defined. 



\begin{lemma}\label{lem:invariantGH}
Within the same equivalence class that contains the current estimate of $(\bm{\Theta}, \bm{\Psi})$,
$\bm{G}$, $\bm{H}$, and $\hat{\bm{H}}$ are uniquely defined. Thus,
the second-order approximation in Eq. \eqref{eq:newton} is uniquely defined.
\end{lemma}
\begin{proof}
According to Eq. \eqref{eqn:noniddiagEig}, the unidentifiability in the diagonal elements of $\bm{\Theta}$ and $\bm{\Psi}$ reduces to the unidentifiability in the eigenvalues of $\bm{\Theta}$ and $\bm{\Psi}$. In the gradient in Eq. \eqref{eq:partialmatrixeiggrad} and Hessians in Eqs. \eqref{eq:partialmatrixeigH} and \eqref{eq:approxH},
the eigenvalues of $\bm{\Theta}$ and $\bm{\Psi}$ always appear as a pair, where the shift $c$ cancels out as follows:
\begin{align*}
\big((\lambda_{\bm{\Theta},l}+c) + (\lambda_{\bm{\Psi},k}-c)\big)^{-1} = (\lambda_{\bm{\Theta},l} + \lambda_{\bm{\Psi},k})^{-1},
\forall (l,k)\in\{1,\ldots,p\}\times\{1,\ldots,q\}.
\end{align*}
Thus, $\bm{G}$, $\bm{H}$, and $\hat{\bm{H}}$ are unaffected by the shift $c$ in the diagonals of $\bm{\Theta}$ and $\bm{\Psi}$.
\end{proof}

Lemma \ref{lem:invariantGH} allows us to solve the problem in Eq. \eqref{eq:newton} to obtain
the descent directions $(D_{\bm{\Theta}},D_{\bm{\Psi}})$, ignoring the constraint on the trace ratio. 
However, if $(D_{\bm{\Theta}},D_{\bm{\Psi}})$ is a solution to the problem in Eq. \eqref{eq:newton}, then $(D_{\bm{\Theta}}-c\bm{I},D_{\bm{\Psi}}+c\bm{I})$ 
for $c\in \mathbb{R}$ is also a solution, forming an equivalence class 
\begin{align} \label{eq:equivD}
S[D_{\bm{\Theta}} \oplus D_{\bm{\Psi}}]= 
    \{ (D_{\bm{\Theta}},D_{\bm{\Psi}})| 
	\left(D_{\bm{\Theta}}-c\bm{I}_p\right)\oplus\left(D_{\bm{\Psi}}+c\bm{I}_q\right)= D_{\bm{\Theta}} \oplus D_{\bm{\Psi}}, 
		c\in \mathbb{R}\}.
\end{align}


\begin{lemma}\label{lem:newtonnonid}
Within the equivalence class in Eq. \eqref{eq:equivD}, the line-search method 
in Algorithm \ref{alg:linesearch} identifies a unique step size $\alpha$.
\end{lemma}
\begin{proof}
We only need to show that $\tvec(D)^T\tvec(\bm{G})$, a term that appears in the computation of $\delta$ 
in Algorithm \ref{alg:linesearch}, is invariant within the equivalence class in Eq. \eqref{eq:equivD}, 
since $D_{\bm{\Theta}}$ and $D_{\bm{\Psi}}$ always appear as Kronecker sum of the two in all the other parts 
of the line search such that $c$ and $-c$ in Eq. \eqref{eq:equivD} cancel out.
We have $\tr(\bm{G}_{\bm{\Theta}}) = \tr(\bm{G}_{\bm{\Psi}})$, which
 can be directly verified from Eqs. \eqref{eq:partialmatrix0} and \eqref{eq:partialmatrixW} in Lemma \ref{lem:gh}
as $\tr(\bm{G}_{\bm{\Theta}}) = \sum_{i=1}^n\frac{1}{n}\tr(\bm{Y}^{iT}\bm{Y}^{i}) - \tr(\bm{W}) = \tr(\bm{G}_{\bm{\Psi}})$. 
From this, we have $\tvec(D)^T\tvec(\bm{G})=\tr(D_{\bm{\Theta}}\bm{G}_{\bm{\Theta}}) + \tr(D_{\bm{\Psi}}\bm{G}_{\bm{\Psi}})
= \tr((D_{\bm{\Theta}}-c\bm{I})\bm{G}_{\bm{\Theta}}) + \tr((D_{\bm{\Psi}}+c\bm{I})\bm{G}_{\bm{\Psi}})
= \tr(D_{\bm{\Theta}}) + \tr(D_{\bm{\Psi}})$, which proves $\tvec(D)^T\tvec(\bm{G})$ is unaffected by the
unidentifiable diagonals of $D_{\bm{\Theta}}$ and $D_{\bm{\Psi}}$.
\end{proof}


With Lemmas \ref{lem:invariantGH} and \ref{lem:newtonnonid},
we can minimize the EiGLasso objective in Eq. \eqref{eq:objrho} by updating $\bm{\Theta}$ and $\bm{\Psi}$
with $(D_{\bm{\Theta}}, D_{\bm{\Psi}})$ and step size $\alpha$ 
and adjusting the diagonals of $\bm{\Theta}$ and $\bm{\Psi}$ using Theorem \ref{thm:uniquepair} after each update.
The theorem below shows this procedure can be simplified even further. 



\begin{theorem}\label{thm:adjustonce}
In EiGLasso, given the trace ratio $\rho=\tr(\bm{\Psi})/\tr(\bm{\Theta})$, it is sufficient to identify 
the diagonals of $\bm{\Theta}$ and $\bm{\Psi}$ only once after convergence. This leads to 
the identical estimate obtained by identifying the diagonals of $\bm{\Theta}$ and $\bm{\Psi}$ 
in every iteration to maintain the trace ratio $\rho$. 
At convergence, the diagonal elements of $(\bm{\Theta},\bm{\Psi})$ are adjusted by the scalar factor
\begin{align}\label{eq:c}
\bm{\Theta} &\gets \bm{\Theta} - \frac{\tr({\bm{\Psi}}) - \rho\tr({\bm{\Theta}})}{q + \rho p} \bm{I}_p.
\end{align}
\end{theorem}
\begin{proof}
Regardless of which member of the equivalence class $S[D_{\bm{\Theta}}\oplus D_{\bm{\Psi}}]$ is used to update 
the current estimate of $(\bm{\Theta}, \bm{\Psi})$, we arrive at the same equivalence class $S[\bm{\Theta}\oplus \bm{\Psi}]$ 
for the estimate of $(\bm{\Theta}, \bm{\Psi})$ after the update. From Lemma \ref{lem:invariantGH},
given this equivalence class $S[\bm{\Theta}\oplus \bm{\Psi}]$, the problem in Eq. \eqref{eq:equivD} 
in the next iteration is uniquely defined. Thus, regardless of whether we adjust the diagonals
to meet the constraint on the trace ratio, the sequence of equivalence class
$S[\bm{\Theta}\oplus \bm{\Psi}]$ over iterations is the same, and it is not necessary to identify the diagonals 
in each iteration of Newton's method. The adjustment in Eq. \eqref{eq:c} is found 
by applying Theorem \ref{thm:uniquepair} to ${\bm{\Theta}}\oplus {\bm{\Psi}}$ with the current estimate of $\bm{\Theta}$
and $\bm{\Psi}$.
\end{proof}


When we set $\rho=\frac{q}{p}$, 
the one-time identification of diagonal parameters with Eq. \eqref{eq:c} in EiGLasso becomes 
identical to 
the identification of the parameters that TeraLasso~\citep{TeraLasso} 
performs at the end of every iteration.
At the end of each iteration, TeraLasso evenly distributes $\tau=\frac{\tr(\bm{\Theta}\oplus\bm{\Psi})}{pq}$
between the re-parameterized $\bm{\Theta}$ and $\bm{\Psi}$ in Eq. \eqref{eq:reparam} 
 and performs the update
$\bm{\Theta} \gets \bm{\Theta} - \frac{\tr(\bm{\Theta})}{p}\bm{I}_p + \frac{\tau}{2}\bm{I}_p$
and similarly for $\bm{\Psi}$. It is straightforward to show
that with $\rho=\frac{q}{p}$ and from Eq. \eqref{eq:tromega}, Eq. \eqref{eq:c} reduces to this update in TeraLasso.
Theorem \ref{thm:adjustonce} is analogous to the simple approach to
handling the unidentifiability of the parameters in Kronecker-product inverse covariance estimation~\citep{YIN2012},
where $\bm{\Theta}\otimes\bm{\Psi} = (c \bm{\Theta})\otimes (\frac{1}{c}\bm{\Psi})$ for any positive constant $c$. 
The parameters are identified by rescaling $\bm{\Theta}$ and $\bm{\Psi}$ 
as $\bm{\Theta} \gets c\bm{\Theta}$ and  $\bm{\Psi} \gets \frac{1}{c}\bm{\Psi}$ with some constant $c$ such that 
$\bm{\Theta}_{11}$ is equal to 1 after convergence, similar to the one-time identification in EiGLasso.








\subsection{Active Set and Automatic Detection of Block-diagonal Structure}
\label{sec:block}



In the graphical lasso, a simple strategy for reducing computation time has been introduced that
detects the block-diagonal structure in the inverse covariance parameter from the sample covariance 
matrix prior to estimation~\citep{witten2011,JMLR:v13:mazumder12a}. Then,
only the parameters within the blocks corresponding to the connected components in the graph need to be estimated.
In Theorem \ref{thm:threshold} below, we show that a similar strategy can be applied to 
EiGLasso with both the exact and approximate Hessian,
to detect the block diagonal structures in $\bm{\Theta}$ and $\bm{\Psi}$ from the sufficient statistics $\bm{S}$ and $\bm{T}$. 

\begin{theorem}\label{thm:threshold}
The block-diagonal structure of $\bm{\Theta}$ can be detected by thresholding $\bm{S}$ such that
$[\bm{\Theta}]_{ij}=0$ iff $|[\bm{S}]_{ij}|\le\gamma_{\bm{\Theta}}$.
Similarly, the block-diagonal structure of $\bm{\Psi}$ can be detected by thresholding $\bm{T}$ 
such that $[\bm{\Psi}]_{ij}=0$ iff $|[\bm{T}]_{ij}|\le\gamma_{\bm{\Psi}}$.
\end{theorem}
\begin{proof}
Let $\partial |\cdot|$ denote the subgradient of the $L_1$ norm of a matrix, i.e., $[\partial|\bm{A}|]_{ij}$ is $1$ if $[\bm{A}]_{ij} > 0$, $-1$ if $[\bm{A}]_{ij} < 0$, and
$[\partial|\bm{A}|]_{ij} \in [-1,1]$ if $[\bm{A}]_{ij} = 0$. Then, the Karush-Kuhn-Tucker conditions \citep{cvxbook,witten2011} for $\bm{\Theta}$ in Eq. \eqref{eq:obj} is
\begin{equation}\label{eq:kkttheta}
    \bm{W}_{\bm{\Theta}} - q\bm{S} - q\gamma_{\bm{\Theta}}\partial|\bm{\Theta}| = 0,
\end{equation}
where $\bm{W}_{\bm{\Theta}}$ is given in Eq. \eqref{eq:partialmatrixW}.
If $\bm{\Theta}$ is block-diagonal, $\bm{W}_{\bm{\Theta}}$ is also block-diagonal, since $\bm{\Theta}$ and $\bm{W}_{\bm{\Theta}}$ 
have the same eigenvectors according to Lemma \ref{lem:gh} and Theorem \ref{thm:gheigcomp}.
This implies that if $[\bm{\Theta}]_{ij}=0$
in the off-diagonal blocks,
$|[\bm{S}]_{ij}|\le\gamma_{\bm{\Theta}}$ in Eq. \eqref{eq:kkttheta}.
The case for $\bm{\Psi}$ can be proven similarly.
\end{proof}


QUIC further showed that their active-set strategy amounts to detecting a block-diagonal structure  
in the first iteration, if the parameters are initialized to a diagonal matrix. This strategy
can be extended to EiGLasso, when an approximate Hessian is used.
To reduce the computation time,
as in QUIC, in each Newton iteration, EiGLasso detects the active set of $\bm{\Theta}$ and $\bm{\Psi}$ 
\begin{align*}
\mathcal{A}_{\bm{\Theta}} =\{(i,j)\mid [\bm{\Theta}]_{ij}\neq0 \text{ or } |[\bm{G}_{\bm{\Theta}}]_{ij}| 
    > q\gamma_{\bm{\Theta}} \}, \;\;
\mathcal{A}_{\bm{\Psi}} =\{(i,j)\mid [\bm{\Psi}]_{ij}\neq0 \text{ or } |[\bm{G}_{\bm{\Psi}}]_{ij}| 
    > p\gamma_{\bm{\Psi}}\},
\end{align*}
and update only the parameters in the active sets during the coordinate descent optimization, 
while setting those in the fixed set to zero.
When $\bm{\Theta}$ and $\bm{\Psi}$ are initialized to diagonal matrices, the approximate Hessian $\hat{\bm{H}}^1$ 
in iteration 1 is diagonal, since the eigenvector matrices $\bm{Q}^1_{\bm{\Theta}}$ and $\bm{Q}^1_{\bm{\Psi}}$ 
are diagonal. Then, the optimization problem in Eq. \eqref{eq:newton} 
decouples into 
a set of optimization problems, each of which involves a single element of $(D^1_{\bm{\Theta}},D^1_{\bm{\Psi}})$
and has a closed-form solution for $[D^1_{\bm{\Theta}}]_{ij}$ as
\begin{equation*}
    [D^1_{\bm{\Theta}}]_{ij} = -q\left( \sum_{k=1}^K[\bm{V}^1_{\bm{\Theta}, k}]_{ii}[\bm{V}^1_{\bm{\Theta}, k}]_{jj} 
	+ (q-K)[\bm{V}^1_{\bm{\Theta}, K}]_{ii}[\bm{V}^1_{\bm{\Theta}, K}]_{jj}\right)^{-1}\mathcal{S}\left([\bm{S}]_{ij},
	\gamma_{\bm{\Theta}}\right),\quad\forall i\neq j.
\end{equation*}
The soft-thresholding 
operator $\mathcal{S}(a,b)$ above is defined as $\mathcal{S}(a,b)=\textrm{sign}(a)(|a|-b)_+$,  
where $(c)_+ = 0$ if $c<0$ and $(c)_+ = c$ if $c>0$.
This closed-form solution is $[D^1_{\bm{\Theta}}]_{ij}=0$ if 
$|[\bm{S}]_{ij}| < \gamma_{\bm{\Theta}}$, which is equivalent to the condition for detecting 
the block-diagonal structure in $\bm{\Theta}$ from $\bm{S}$ in Theorem \ref{thm:threshold}. 
The case for $\bm{\Psi}$ can be shown similarly.






\section{Convergence Analysis}\label{sec:convg}

We examine the properties of the line-search method in Algorithm \ref{alg:linesearch}
and analyze the global and local convergence of EiGLasso.

\subsection{Line Search Properties}

EiGLasso inherits some of the line-search properties shown for QUIC \citep{quic}.
This is because our objective in Eq. \eqref{eq:obj} can be written in terms of 
$\bm{\Omega}=\bm{\Theta}\oplus\bm{\Psi}$ in the form that resembles the objective of QUIC
\begin{equation*}\label{eq:objomega}
    f(\bm{\Omega}) =
    \tr(\bm{U} \bm{\Omega})
    - \log|\bm{\Omega}| + \| \bm{\Gamma}\circ\bm{\Omega} \|_{1,\text{off}}, 
\end{equation*}
where $\bm{U}=\frac{1}{n}\sum_{i=1}^n \tvec(\bm{Y}^i)^T\tvec(\bm{Y}^i)$, $\bm{\Gamma}=\gamma_{\bm{\Theta}}\mathbbm{1}_p\oplus\gamma_{\bm{\Psi}}\mathbbm{1}_q$ with an $r\times r$ matrix of one's $\mathbbm{1}_r$, and $\circ$ is an element-wise multiplication operator. 
Then, the EiGLasso's update
$\bm{\Theta} \gets \bm{\Theta} + \alpha D_{\bm{\Theta}}$ and $\bm{\Theta} \gets \bm{\Theta} + \alpha D_{\bm{\Theta}}$ can be written as 
$\bm{\Omega} \gets \bm{\Omega} + \alpha D_{\bm{\Omega}}$, where $D_{\bm{\Omega}} = D_{\bm{\Theta}} \oplus D_{\bm{\Psi}}$, since $\bm{\Omega} + \alpha D_{\bm{\Omega}} 
    =(\bm{\Theta}+\alpha D_{\bm{\Theta}})\oplus(\bm{\Psi}+\alpha D_{\bm{\Psi}})$.
We extend these results for $\bm{\Omega}$, to prove the results for 
individual $\bm{\Theta}$ and $\bm{\Psi}$, where $(\bm{\Theta}, \bm{\Psi}) \in \mathbb{KS}^{p,q}$, 
when  unlike in QUIC the exact Hessian $\bm{H}$ is not positive definite everywhere (Theorem \ref{thm:Hpd}), 
and when the approximate Hessian $\hat{\bm{H}}$ is used.

We begin by showing that both $\bm{H}$ and $\hat{\bm{H}}$ have bounded eigenvalues
in the level set
$\mathcal{U}=\{(\bm{\Theta},\bm{\Psi})\mid \bm{\Theta}\oplus\bm{\Psi}\in\mathbb{S}^{pq}_{++},\;\;f(\bm{\Theta},\bm{\Psi})
\le f(\bm{\Theta}^0,\bm{\Psi}^0)\}$
and are Lipschitz-continuous. This result will be used to prove the line-search properties and global and local convergence.   

\begin{lemma}\label{lemma:eigvalbound}
In the level set $\mathcal{U}$,
$\bm{H}$ and $\hat{\bm{H}}$ are Lipschitz continuous and have bounded eigenvalues:
\begin{align*}
    \min\{p,q\}\bar{\lambda}^{-2}\bm{I}_{p^2+q^2} \preceq
    & \bm{H}\preceq  (p+q)\underline{\lambda}^{-2}\bm{I}_{p^2+q^2},
    \nonumber \\
    \min\{p,q\}\bar{\lambda}^{-2}\bm{I}_{p^2+q^2} \preceq 
    & \hat{\bm{H}} \preceq  \max\{p,q\}\underline{\lambda}^{-2}\bm{I}_{p^2+q^2},
\end{align*}
for some constants $\underline{\lambda},\overline{\lambda}>0$ that depend on $\gamma_{\bm{\Theta}}$, $\gamma_{\bm{\Psi}}$, $f(\bm{\Theta}^0,\bm{\Psi}^0)$, and $\{\bm{Y}^1,\ldots, \bm{Y}^n\}$.
\end{lemma}
\begin{proof}
We bound the eigenvalues of $\bm{\Omega}=\bm{\Theta}\oplus \bm{\Psi}$, $\bm{\Theta}$, and $\bm{\Psi}$, and
use these bounds to bound the eigenvalues of $\bm{H}$ and $\hat{\bm{H}}$.
When the diagonals are not identified,
it directly follows from Lemma 2 in \citet{quic}
that all EiGLasso iterates of $(\bm{\Theta},\bm{\Psi})$ 
are contained in the set with bounded eigenvalues of $\bm{\Theta}\oplus \bm{\Psi}$  
    \begin{equation}\label{eq:compactset_omega}
        \mathcal{C}=\{(\bm{\Theta},\bm{\Psi})\mid \underline{\lambda}\bm{I}\preceq\bm{\Theta}\oplus\bm{\Psi}\preceq\bar{\lambda}\bm{I}\}.
    \end{equation}
Next, given a fixed trace ratio $\rho=\frac{\tr(\bm{\Psi})}{\tr(\bm{\Theta})}$,  
we bound the eigenvalues of $\bm{\Theta}$ and $\bm{\Psi}$.
Since Eq. \eqref{eq:compactset_omega} implies 
$\underline{\lambda}\bm{I}\preceq\bm{\Lambda}_{\bm{\Theta}}\oplus\bm{\Lambda}_{\bm{\Psi}}\preceq\bar{\lambda}\bm{I}$,
we apply Theorem \ref{thm:uniquepair} to $\bm{\Lambda}_{\bm{\Omega}}=\bm{\Lambda}_{\bm{\Theta}}\oplus\bm{\Lambda}_{\bm{\Psi}}$ 
to identify $\bm{\Lambda}_{\bm{\Theta}}$ and
$\bm{\Lambda}_{\bm{\Psi}}$ in the equivalence class in Eq. \eqref{eqn:noniddiagEig} given $\rho$, thus, identifying
the diagonals of $\bm{\Theta}$ and $\bm{\Psi}$. 
We bound each of the two terms in Eq. \eqref{eq:diagmapping} for $\bm{\Theta}$ as
\begin{align*}
    q\underline{\lambda}\bm{I}\preceq \sum_{i=1}^q(\bm{I}_p\otimes\bm{e}_{q,i})^T\bm{\Lambda}_{\bm{\Omega}}(\bm{I}_p\otimes\bm{e}_{q,i})\preceq q\bar{\lambda}\bm{I} \quad\textrm{and}\quad
    \frac{\rho pq}{q+\rho p}\underline{\lambda}\bm{I}\preceq\frac{\rho}{q+\rho p}\tr(\bm{\Lambda}_{\bm{\Omega}})\bm{I}\preceq\frac{\rho pq}{q+\rho p}\bar{\lambda}\bm{I},
\end{align*}
and combine these to obtain the eigenvalue bounds for $\bm{\Theta}$ and similarly for $\bm{\Psi}$
    \begin{align}
        \mathcal{C}_{\rho} = 
        \left\{(\bm{\Theta},\bm{\Psi})\mid
        \Big(\underline{\lambda}-\frac{\rho p}{q + \rho p}\bar{\lambda}\Big)\bm{I}\preceq\bm{\Theta}\preceq\Big(\bar{\lambda}-\frac{\rho p}{q + \rho p}\underline{\lambda}\Big)\bm{I}, 
        \quad\quad\quad\quad 
        \right.\nonumber \\
        \left.\Big(\underline{\lambda}-\frac{q}{q + \rho p}\bar{\lambda}\Big)\bm{I}
        \preceq\bm{\Psi}\preceq\Big(\bar{\lambda}-\frac{q}{q + \rho p}\underline{\lambda}\Big)\bm{I}\right\}. \label{eq:compactset}
    \end{align}
    
From Eq. \eqref{eq:compactset_omega} and Theorem \ref{thm:Hpd}, 
we obtain the bound on the eigenvalues
of $\bm{H}$. From Eq. \eqref{eq:compactset_omega} and Theorem \ref{thm:Happroxpd}, 
we obtain the bound on the eigenvalues
of $\hat{\bm{H}}$, 
since $\lambda_{\hat{\bm{H}}_{\bm{\Theta}},\min} \ge q\bar{\lambda}^{-2}$,
    $\lambda_{\hat{\bm{H}}_{\bm{\Theta}},\max} \le q\underline{\lambda}^{-2}$,
    $\lambda_{\hat{\bm{H}}_{\bm{\Psi}},\min} \ge p\bar{\lambda}^{-2}$,
    $\lambda_{\hat{\bm{H}}_{\bm{\Psi}},\max} \le p\underline{\lambda}^{-2}$.
On the set in Eq. \eqref{eq:compactset},
since the log-determinant is a continuous function of class $C^\infty$, and a continuous function on a compact set is bounded, both $\bm{H}$ and $\hat{\bm{H}}$ are locally Lipschitz-continuous. 
\end{proof}







Now, for EiGLasso with exact and approximate Hessian,
we show that the following three line-search properties hold:
the line search method is guaranteed to terminate  for any symmetric matrices $D_{\bm{\Theta}}$ and $D_{\bm{\Psi}}$,
as the two line search conditions in Algorithm \ref{alg:linesearch} 
are satisfied for some step size $\alpha$ (Lemma \ref{lem:term}); 
the update with the Newton direction is guaranteed to decrease the objective
(Lemma \ref{lem:delta_exact}); and EiGLasso with the exact Hessian
is guaranteed to enter 
pure-Newton phase where the step size $\alpha=1$ is always chosen.
We state and prove the first two properties.
The proof of the last property follows directly from the proof in \citet{Tseng2009} and \citet{quic}. 


\begin{lemma}\label{lem:term}
For any $\bm{\Theta}$ and $\bm{\Psi}$, where $\bm{\Theta}\oplus\bm{\Psi}\succ0$, and symmetric matrices $D_{\bm{\Theta}}$ and $D_{\bm{\Psi}}$ 
for descent directions found with either the exact or inexact Hessian, there exists a step size $\tilde{\alpha}\in(0,1]$ 
such that for all $\alpha<\tilde{\alpha}$ the two conditions in the line search 
in Algorithm \ref{alg:linesearch} are satisfied.
\end{lemma}
\begin{proof}
If $\alpha < \lambda_{\min}(\bm{\Theta}\oplus\bm{\Psi})/\|D_{\bm{\Theta}}\oplus D_{\bm{\Psi}}\|_2$, 
the updated estimates $\bm{\Theta} + \alpha D_{\bm{\Theta}}$ and $\bm{\Psi} + \alpha D_{\bm{\Psi}}$
satisfy $(\bm{\Theta} + \alpha D_{\bm{\Theta}})\oplus(\bm{\Psi} + \alpha D_{\bm{\Psi}}) \succ 0$,
since 
$(\bm{\Theta} + \alpha D_{\bm{\Theta}})\oplus(\bm{\Psi} + \alpha D_{\bm{\Psi}})=\bm{\Theta}\oplus\bm{\Psi} + \alpha(D_{\bm{\Theta}} \oplus D_{\bm{\Psi}})$ 
and we have $\bm{\Theta}\oplus\bm{\Psi}\succ 0$ and $\|\alpha (D_{\bm{\Theta}} \oplus D_{\bm{\Psi}})\|_2  < \lambda_{\min}(\bm{\Theta}\oplus\bm{\Psi})$. 
Thus, the first line-search condition in Algorithm \ref{alg:linesearch} is satisfied.
From Lemma 1 in \citet{Tseng2009} and Proposition 3 in \citet{quic}, it is straightforward
to show the second condition in Algorithm \ref{alg:linesearch} is satisfied.
\end{proof}



\begin{lemma}\label{lem:delta_exact}
Let $\tvec(D)=[\tvec(D_{\bm{\Theta}})^T,\tvec(D_{\bm{\Psi}})^T]^T$ for all symmetric $D_{\bm{\Theta}}$ and $D_{\bm{\Psi}}$.
With the exact Hessian $\bm{H}$, 
$\delta$ in Algorithm \ref{alg:linesearch} is upper bounded if not at the optimum,
\begin{align}\label{eq:deltaexactbound}
\delta \le -\tvec(D)^T\bm{H}\tvec(D) &\le -\lambda_{\bm{H},\min_0}\|\tvec(D)\|_2^2< 0, 
\end{align} 
where $\|\tvec(D)\|_2^2 = \|D_{\bm{\Theta}}\|_F^2 + \|D_{\bm{\Psi}}\|_F^2$, and 
$\lambda_{\bm{H},\min_0}$ is given in Theorem \ref{thm:Hpd}.
With the approximate Hessian $\hat{\bm{H}}$, 
$\delta$ is upper bounded everywhere 
\begin{align} 
\delta \le -\tvec(D)^T\hat{\bm{H}}\tvec(D)   & \le   -\lambda_{\hat{\bm{H}},\min}\|\tvec({D})\|_2^2  <  0, \label{eq:delta_approx}
\end{align} 
where $\lambda_{\hat{\bm{H}},\min}$ is given in Theorem \ref{thm:Happroxpd}.
\end{lemma}
\begin{proof}
    For $\bm{H}$, the first inequality in Eq. \eqref{eq:deltaexactbound} can be shown by a  
    straightforward application of Lemma 1 and Theorem 1 in \citet{Tseng2009} and Proposition 4 
    in \citet{quic}. To prove the second and third inequalities, since our $\bm{H}$ is not 
    positive definite everywhere, we need to show that $\tvec(D)$ is outside of the nullspace of $\bm{H}$ described in Theorem \ref{thm:Hpd}, unless EiGLasso is at the optimum, where $D={0}$.
 The null space of $\bm{H}$, $\{\tvec({D})|\tvec({D}_{\bm{\Theta}}\oplus {D}_{\bm{\Psi}})={0}\}$, 
 is equivalent to $\{\tvec(D)|D_{\bm{\Theta}} = c\bm{I}, D_{\bm{\Psi}} = -c\bm{I} \textrm{ for } c\in \mathbb{R}\}$, 
 which is the equivalence class of the optimality condition $D={0}$.

     This proves the third inequality in Eq. \eqref{eq:deltaexactbound} that holds except when $D=0$.
     For $\hat{\bm{H}}$, the first inequality in Eq. \eqref{eq:delta_approx} can be again shown 
     from Lemma 1 and Theorem 1 in \citet{Tseng2009} and Proposition 4 in \citet{quic}.
    The second and third inequalities hold since $\hat{\bm{H}}$ is positive definite. 
\end{proof}



\subsection{Convergence Analysis}


To show the global convergence of EiGLasso, as in QUIC, we use a more general non-smooth optimization framework, 
the block coordinate descent studied in \citet{Tseng2009}.
EiGLasso satisfies the following two conditions required 
to guarantee that the block coordinate descent algorithm converges to the global optimum.
First, the objective function of EiGLasso has a bounded positive definite 
exact or approximate Hessian in the level set $\mathcal{U}$: 
$a\bm{I} \preceq \bm{H}, \hat{\bm{H}} \preceq b\bm{I}$  for some positive constants $a,b\in\mathbb{R}^+$, 
according to Theorems \ref{thm:Hpd} and \ref{thm:Happroxpd}, and Lemma \ref{lemma:eigvalbound}. 
Second, EiGLasso with the exact or approximate Hessian chooses a subset of variables to be updated in each iteration 
according to the \textit{Gauss-Seidel} rule: with $T=2$, at iteration $t$, 
EiGLasso updates one of the two subsets of variables,
$\mathcal{J}^{2t}=\mathcal{A}^{t}_{\bm{\Theta}}\sqcup\mathcal{A}^{t}_{\bm{\Psi}}$ and
$\mathcal{J}^{2t+1}=\mathcal{J}\setminus\mathcal{A}^{t}_{\bm{\Theta}}\sqcup\mathcal{A}^{t}_{\bm{\Psi}}$,    
where $\mathcal{J}=\mathcal{J}^{2t} \cup \mathcal{J}^{2t+1}$ is the entire set of variables 
and $\mathcal{A}^{t}_{\bm{\Theta}}$ and $\mathcal{A}^{t}_{\bm{\Psi}}$ are active sets in Algorithm \ref{alg:eiglasso}. 
Therefore, EiGLasso is guaranteed to converge to the global optimum according to \citet{Tseng2009}.




Now we analyze the local convergence of EiGLasso.
We adopt a similar strategy used in QUIC \citep{quic}: 
convergence analysis on a smooth function is applied near the global optimum,  
where the $L_1$-regularized non-smooth objective becomes locally smooth. 



\begin{theorem}\label{thm:convgrateexact}
Near the optimum, where step size $\alpha=1$ is chosen,
EiGlasso with the exact Hessian $\bm{H}$
converges to the optimum  quadratically.
EiGLasso with the approximate Hessian $\hat{\bm{H}}$
converges to the optimum linearly.
\end{theorem}
\begin{proof}
Since the exact Hessian $\bm{H}$ is Lipschitz continuous from Lemma \ref{lemma:eigvalbound},
the proof for EiGLasso with the exact Hessian follows from Lemma 2.5 and Theorem 3.1 of
\citet{dunn1980}. 
While the convergence analysis in \citet{dunn1980} assumes that the Hessian is positive definite,
the exact Hessian in EiGLasso is not positive definite everywhere. 
However, according to Theorem \ref{thm:Hpd}, our Hessian 
is positive definite  for all iterates $(\bm{\Theta}^t,\bm{\Psi}^t)\in \mathbb{KS}^{pq}_{++}$,
since  $\mathbb{KS}^{pq}_{++}$ is outside of the null space of $\bm{H}$, so
the analysis in \citet{dunn1980} can be applied to EiGLasso.
The proof for EiGLasso with the approximate Hessian $\hat{\bm{H}}$ follows from the analysis of steepest descent with the quadratic norm~\citep{cvxbook}, 
given the bounded eigenvalues of $\bm{H}$ according to Lemma \ref{lemma:eigvalbound}.
\end{proof}













\begin{comment}
According to \citet{Tseng2009}, EiGLasso with approximate Hessian globally converges at least linearly to the optimum. The local convergence rate of EiGLasso with approximate Hessian can be analyzed in the framework of inexact Newton methods \citep{inexact}. After EiGLasso enters the {\color{red}pure Newton phase}, its convergence rate depends on the residual defined as
\begin{equation}\label{eq:residual}
    \bm{r}^t = (\bm{I} - \bm{H}^t(\hat{\bm{H}}^t)^{-1})\bm{G}^t.
\end{equation}
More specifically, if $\|\bm{r}^t\|_2 \le \zeta_t \|\bm{G}^t\|_2$ with $0\le\zeta_t<1$, then for some $\zeta_t<\bar{\zeta}<1$ for all $t$, the following is satisfied:
\begin{equation*}
    \left\| \bm{H}^*\begin{bmatrix}
\tvec(\bm{\Theta}^{t+1}-\bm{\Theta}^*)\\
\tvec(\bm{\Psi}^{t+1}-\bm{\Psi}^*)
\end{bmatrix}  \right\|_2 \le \bar{\zeta} 
\left\| \bm{H}^*\begin{bmatrix}
\tvec(\bm{\Theta}^{t}-\bm{\Theta}^*)\\
\tvec(\bm{\Psi}^{t}-\bm{\Psi}^*)
\end{bmatrix}  \right\|_2.    
\end{equation*}
Given this and Eq. \eqref{eq:residual}, we can see that how close $\bm{H}^t(\hat{\bm{H}}^t)^{-1}$ is to an identity matrix in Euclidean or spectral norm can be an indirect indicator of the convergence rate since
\begin{align*}
    \|\bm{r}^t\|_2 &= \|(\bm{I} - \bm{H}^t(\hat{\bm{H}}^t)^{-1})\bm{G}^t\|_2 \\
    &\le \|\bm{I} - \bm{H}^t(\hat{\bm{H}}^t)^{-1}\|_2\|\bm{G}^t\|_2\\
    &\le \|\bm{I} - \bm{H}^t(\hat{\bm{H}}^t)^{-1}\|_F\|\bm{G}^t\|_2.
\end{align*}
It implies that
\begin{align*}
    \zeta_t &\le \|\bm{I} - \bm{H}^t(\hat{\bm{H}}^t)^{-1}\|_2\\
    &\le \|\hat{\bm{H}}^t - \bm{H}^t\|_2 \| (\hat{\bm{H}}^t)^{-1} \|_2\\
    &\le \frac{\|\hat{\bm{H}} - \bm{H}\|_2}{\min\{p,q\}\bar{\lambda}^2}.
\end{align*}
\end{comment}




\section{Experiments} \label{sec:exp}

We compare the performance of EiGLasso with that of TeraLasso \citep{TeraLasso} on simulated data 
and on real-world data from genomics and finance. TeraLasso is the state-of-the art method
for Kronecker-sum inverse covariance estimation, and has been shown to be substantially
more efficient than BiGLasso~\citep{biglasso}, so we did not include BiGLasso in our experiments.
We implemented EiGLasso 
in C++ with the sequential version of Intel Math Kernel Library. 
We downloaded the authors' implementation of TeraLasso and modified it to perform more iterations during line search, 
when the safe-step approach suggested by the authors failed 
to find a step-size that satisfies the positive definite condition on $\bm{\Theta}\oplus\bm{\Psi}$.
All experiments were run on a single core of Intel(R) Xeon(R) CPU E5-2630 v3 @ 2.40GHz.
In all of our experiments, we selected the regularization parameters 
$\gamma=\gamma_{\bm{\Theta}}=\gamma_{\bm{\Psi}}$ for EiGLasso 
and used the selected $\gamma$ for TeraLasso,
because EiGLasso is significantly faster than TeraLasso 
and minimizes the same objective as TeraLasso.
To assess convergence, we used the criterion that the decrease in the objective function 
value $f^t$ at iteration $t$
satisfies the condition $\left|\frac{f^t-f^{t-1}}{f^{t}}\right| < \epsilon$ for three consecutive iterations. 


\begin{figure*}[t!]
\centering
\setlength\tabcolsep{0pt}
\begin{tabular}{ccl}
     \begin{subfigure}{0.33\linewidth}
        \includegraphics[width=\linewidth]{plot_simul4/100x100_logiter_marked_legend5.pdf}
        \caption*{(a)}
    \end{subfigure} &
    \begin{subfigure}{0.33\linewidth}
        \includegraphics[width=\linewidth]{plot_simul4/100x100_logiter_marked7.pdf}
        \caption*{(b)}
    \end{subfigure} &
    \multirow{2}{0.3\linewidth}[-1em]{
        \begin{subfigure}{\linewidth}
            \includegraphics[width=\linewidth]{plot_simul4/100x100PR3.pdf}
            \caption*{(e)}
        \end{subfigure}
    }
    \\
    \begin{subfigure}{0.33\linewidth}
        \includegraphics[width=\linewidth]{plot_simul4/100x100_logtime_marked5.pdf}
        \caption*{(c)}
    \end{subfigure} &
    \begin{subfigure}{0.33\linewidth}
        \includegraphics[width=\linewidth]{plot_simul4/100x100_logtime_marked7.pdf}
        \caption*{(d)}
    \end{subfigure} &
\end{tabular}

    \caption{
Comparison of the convergence of EiGLasso and TeraLasso on data simulated from random graphs.
Objective values over iterations are shown for two datasets
in (a) and (b). Objective values over time are shown in (c) and (d) for the same two datasets. 
All methods were run until they reached the objective that EiGLasso with the exact Hessian reached with the convergence criterion $\epsilon=10^{-3}$.
The `$\bullet$', `$\times$', `+', `$\Box$', `$\bigcirc$' and `$\diamondsuit$' mark the points that satisfy 
the convergence criteria $\epsilon=10^{-3},\,10^{-4},\,10^{-5},\,10^{-6},\,10^{-7}$, and $10^{-8}$, respectively. 
Precision-recall curves for EiGLasso at $\epsilon=10^{-3}$ and TeraLasso at different $\epsilon$'s averaged over 10 datasets are shown in (e).
Graph size $p, q=100$ was used.
}
    \label{fig:100r}
\end{figure*}


\subsection{Simulated Data}\label{subsec:simulation}



We compared EiGLasso and TeraLasso on data simulated from the known $\bm{\Theta}$ and $\bm{\Psi}$. 
We used the true $\bm{\Theta}$ and $\bm{\Psi}$ of different sizes ($p,q=$ 100, 200, 500, 1000, 2000, and 5000),
assuming two types of graph structures.
\begin{itemize}
\item
{\bf Random graph}: To set the ground-truth $\bm{\Theta}$, 
we first generated a sparse $p\times p$ matrix $\bm{A}$  by assigning $-1$, 0, or 1 to each element with 
probabilities $\frac{1-\rho}{2}$, $\rho$, and $\frac{1-\rho}{2}$, respectively. We chose $\rho$ such that 
the number of non-zero elements of $\bm{\Theta}$ is $10p$. 
To ensure $\bm{\Theta}$ is positive definite, we set $\bm{\Theta}$ to $\bm{A}\bm{A}^T$ after adding $\sigma + 10^{-4}$ with $\sigma\sim\text{Unif}(0,0.1)$
to each diagonal element of $\bm{A}\bm{A}^T$. 
\item{\bf Graph with clusters}: We set $\bm{\Theta}$ to a block-diagonal matrix such that each block corresponds to a cluster. 
For graphs with $p=100$ and $200$, we assumed five blocks, each with size $\frac{p}{5}\times \frac{p}{5}$. 
For larger graphs with $p=500$, 1000, 2000, and $5000$, we assumed 10 blocks, each with size $\frac{p}{10}\times \frac{p}{10}$. 
Each block was generated as a random graph described above, setting $\rho$ so that we have $p$ nonzero elements in the block. 
\end{itemize}
The ground-truth $\bm{\Psi}$ was set similarly. Given these parameters, we simulated matrix-variate data 
from Gaussian distribution with mean zeros and inverse covariance $\bm{\Theta} \oplus \bm{\Psi}$.

First, we evaluated TeraLasso and EiGLasso on simulated data with graph size $p, q=100$.
The two methods were compared in terms of computation time and the number of iterations required to reach 
the same level of optimality,  which we define as the objective value that
EiGLasso with the exact Hessian converged to with the convergence criterion $\epsilon=10^{-3}$. 
The regularization parameters were selected such that the number of non-zero elements in the estimated 
parameters roughly  matches that of the true parameters.
EiGLasso with approximate Hessian was run with different $K$'s ranging from 
$K=1$ to $K=100$.
Results are shown for four datasets, two from random graphs (Figures \ref{fig:100r}(a)-(d))
and two from graphs with clusters (Figures \ref{fig:100bd}(a)-(d)).
Regardless of the degree of Hessian approximation, EiGLasso required significantly 
fewer iterations than TeraLasso (Figures \ref{fig:100r}(a)-(b) and 
\ref{fig:100bd}(a)-(b)),
as expected since methods that use the second-order information in general converge in fewer iterations 
than the first-order methods. 
With the exact Hessian and approximate Hessian with large $K$,
EiGLasso took longer in each iteration than TeraLasso, but overall required
two to three times less computation time than TeraLasso
(Figures \ref{fig:100r}(c)-(d) and \ref{fig:100bd}(c)-(d)).
As we reduced $K$, the time taken by EiGLasso 
decreased substantially, and EiGLasso with $K=1$ achieved two to three orders-of-magnitude 
speed-up compared to TeraLasso. 



\begin{figure*}[t!]
\centering
\setlength\tabcolsep{0pt}
\begin{tabular}{ccl}
   \begin{subfigure}{0.33\linewidth}
        \includegraphics[width=\linewidth]{plot_simul4/100x100bd_logiter_marked_legend3.pdf}
        \caption*{(a)}
    \end{subfigure}
    &
    \begin{subfigure}{0.33\linewidth}
        \includegraphics[width=\linewidth]{plot_simul4/100x100bd_logiter_marked8.pdf}
        \caption*{(b)}
    \end{subfigure} &
    \multirow{2}{0.3\linewidth}[-1em]{
        \begin{subfigure}{\linewidth}
            \includegraphics[width=\linewidth]{plot_simul4/100x100bdPR3.pdf}
            \caption*{(e)}
        \end{subfigure}
    }
    \\
   \begin{subfigure}{0.33\linewidth}
        \includegraphics[width=\linewidth]{plot_simul4/100x100bd_logtime_marked3.pdf}
        \caption*{(c)}
    \end{subfigure}
    &
    \begin{subfigure}{0.33\linewidth}
        \includegraphics[width=\linewidth]{plot_simul4/100x100bd_logtime_marked8.pdf}
        \caption*{(d)}
    \end{subfigure} &
\end{tabular}

    \caption{
Comparison of the convergence of EiGLasso and TeraLasso on data simulated from graphs with clusters.
Objective values over iterations are shown for two datasets
in (a) and (b). Objective values over time are shown in (c) and (d) for the same two datasets. 
All methods were run until they reached the objective that EiGLasso with the exact Hessian reached with the convergence criterion $\epsilon=10^{-3}$.
The `$\bullet$', `$\times$', `+', `$\Box$', and `$\bigcirc$' mark the points that satisfy 
$\epsilon=10^{-3},\,10^{-4},\,10^{-5},\,10^{-6},$ and $\,10^{-7}$, respectively. 
Precision-recall curves for EiGLasso at $\epsilon=10^{-3}$ and TeraLasso at different $\epsilon$'s averaged over 10 datasets are shown in (e).
Graph size $p, q=100$ was used.
}
    \label{fig:100bd}
\end{figure*}



\begin{table}[h!]
    \centering
    \begin{tabular}{|c|c|rrrrrr|r|}
    \hline
    \multicolumn{1}{|c|}{\multirow{2}{*}{Graph}} & \multicolumn{1}{c}{\multirow{2}{*}{$p,q$}} & \multicolumn{6}{|c|}{EiGLasso} & TeraLasso\\
    && $K=1$ & $K=2$ & $K=5$ & $K=10$ & $K=50$ & $K=100$ & \\
    \hline
    \multirow{6}{*}{\shortstack{Random\\Graphs}} & 100 & 1 & 4 & 6 & 40 & 729 & 1019 & 4194 \\
    & 200 & 7 & 164 & 551 & 3059 & 27049 & 70311 & 16765 \\
    & 500 & 29 & 155 & 2252 & 5994 & 31826 &&\\
    & 1000 & 515 & 11601 & 39568 &&&&\\
    & 2000 & 6361 & $28539$ &&&&&\\
    & 5000 & 10458 & $62685$ &&&&&\\
    \hline
    \multirow{6}{*}{\shortstack{
    Graphs\\with\\Clusters}} & 100 & 2 & 21 & 28 & 629 & 6516 & 12443 & 21115 \\
    & 200 & 11 & 410 & 930 & 1042 & 11281 & 59813 & $67669$ \\
    & 500 & 28 & 5231 & 9252 & 10549 &&&\\
    & 1000 & 2169 & 12614 & 58168 &&&&\\
    & 2000 & 9310 & $41524$ &&&&&\\
    & 5000 & 10461 & $67608$ &&&&&\\
    \hline
    \end{tabular}
    \caption{Comparison of computation time (in seconds) on simulated data. 
    }
    \label{tab:simresult}
\end{table}

We compared TeraLasso and EiGLasso on the accuracy of the recovered non-zero elements in $\bm{\Theta}$ and $\bm{\Psi}$, 
using 10 datasets simulated as above.
The precision-recall curves averaged over the 10 simulated datasets 
(Figure \ref{fig:100r}(e) and \ref{fig:100bd}(e)) show that  
the TeraLasso estimates at the convergence criterion $\epsilon=10^{-3}$ and $10^{-4}$ are significantly inferior to
the EiGLasso estimates at the same $\epsilon$. TeraLasso needed $\epsilon=10^{-6}$ 
to achieve the accuracy similar to EiGLasso at $\epsilon=10^{-3}$,
whereas more stringent criteria such as $\epsilon=10^{-4}$ and $10^{-5}$ 
led to nearly identical results in EiGLasso.
Thus, in all experiments in the rest of Section \ref{sec:exp},
we used $\epsilon=10^{-3}$ as convergence criterion for EiGLasso with all $K$'s, 
but ran TeraLasso until it reached the similar objective value 
that EiGLasso reached with $\epsilon=10^{-3}$, 
which typically required $\epsilon = 10^{-6}$, $10^{-7},$ or $10^{-8}$. 




\begin{figure*}[b!]
\centering
\setlength\tabcolsep{0pt}
\begin{tabular}{c@{}c@{}c}
    \begin{subfigure}{0.33\linewidth}
        \includegraphics[width=\linewidth]{plot_eig4/50x50_eigval_min5.pdf}
        \caption{}
    \end{subfigure} 
    &
   \begin{subfigure}{0.33\linewidth}
        \includegraphics[width=\linewidth]{plot_eig4/50x50_eigval_max5.pdf}
        \caption{}
    \end{subfigure}
    &
    \begin{subfigure}{0.33\linewidth}
        \includegraphics[width=\linewidth]{plot_eig4/50x50_eigval_all_opt5.pdf}
        \caption{}
    \end{subfigure} 
    \\
    \begin{subfigure}{0.33\linewidth}
        \includegraphics[width=\linewidth]{plot_eig4/50x50_logiter_optdiff5.pdf}
        \caption{}
    \end{subfigure} 
    &
   \begin{subfigure}{0.33\linewidth}
        \includegraphics[width=\linewidth]{plot_eig4/50x50_logiter_Hdiff5.pdf}
        \caption{}
    \end{subfigure}
    &
\end{tabular}
\caption{
    Effects of the Hessian approximation on the convergence behavior of EiGLasso 
    for data simulated from a random graph. 
    (a) The summands in $\lambda_{\bm{H}_{\bm{\Theta}},\min}=\sum_{j=1}^q (\lambda_{\bm{\Theta}, p} 
    + \lambda_{\bm{\Psi},j})^{-2}$ and 
    $\lambda_{\hat{\bm{H}}_{\bm{\Theta}},\min}=\sum_{j=1}^K (\lambda_{\bm{\Theta}, p} + \lambda_{\bm{\Psi},j})^{-2} 
    + (q-K)(\lambda_{\bm{\Theta}, p} + \lambda_{\bm{\Psi},K})^{-2}$ for EiGLasso with $K=5$. 
    (b) The summands in $\lambda_{\bm{H}_{\bm{\Theta}},\max}$ and $\lambda_{\hat{\bm{H}}_{\bm{\Theta}},\max}$ 
    as in Panel (a).
    (c) The eigenvalues of $\bm{H}_{\bm{\Theta}}$ and $\hat{\bm{H}}_{\bm{\Theta}}$, sorted 
    with the eigenvalues of $\bm{H}_{\bm{\Theta}}$.
    In Panels (a)-(c), the results are shown 
    at the 5th iteration (top) and at the last iteration (bottom).
    (d) Error measured as 
   $\left\|\protect\begin{bmatrix}
    \tvec(\bm{\Theta}^{t}-\bm{\Theta}^*) \\
    \tvec(\bm{\Psi}^{t}-\bm{\Psi}^*)
    \protect\end{bmatrix}\right\|_2$ 
    over iterations.
    (e) $\|\hat{\bm{H}}-\bm{H}\|_F$ over iterations.
    A dataset simulated from a random graph with $p, q=50$ was used.
}
    \label{fig:50rndemp}
\end{figure*}


\begin{figure*}[b!]
\centering
\setlength\tabcolsep{0pt}
\begin{tabular}{c@{}c@{}c}
    \begin{subfigure}{0.33\linewidth}
        \includegraphics[width=\linewidth]{plot_eig4/50x50bd_eigval_min5.pdf}
        \caption{}
    \end{subfigure} 
    &
   \begin{subfigure}{0.33\linewidth}
        \includegraphics[width=\linewidth]{plot_eig4/50x50bd_eigval_max5.pdf}
        \caption{}
    \end{subfigure}
    &
    \begin{subfigure}{0.33\linewidth}
        \includegraphics[width=\linewidth]{plot_eig4/50x50bd_eigval_all_opt5.pdf}
        \caption{}
    \end{subfigure} 
    \\
    \begin{subfigure}{0.33\linewidth}
        \includegraphics[width=\linewidth]{plot_eig4/50x50bd_logiter_optdiff5.pdf}
        \caption{}
    \end{subfigure} 
    &
   \begin{subfigure}{0.33\linewidth}
        \includegraphics[width=\linewidth]{plot_eig4/50x50bd_logiter_Hdiff5.pdf}
        \caption{}
    \end{subfigure}
    &
\end{tabular}
\caption{
    Effects of Hessian approximation on the convergence behavior of EiGLasso
    for data simulated from a graph with clusters. 
    (a) The summands in $\lambda_{\bm{H}_{\bm{\Theta}},\min}=\sum_{j=1}^q (\lambda_{\bm{\Theta}, p} 
    + \lambda_{\bm{\Psi},j})^{-2}$ and 
    $\lambda_{\hat{\bm{H}}_{\bm{\Theta}},\min}=\sum_{j=1}^K (\lambda_{\bm{\Theta}, p} + \lambda_{\bm{\Psi},j})^{-2} 
    + (q-K)(\lambda_{\bm{\Theta}, p} + \lambda_{\bm{\Psi},K})^{-2}$ for EiGLasso with $K=5$. 
    (b) The summands in $\lambda_{\bm{H}_{\bm{\Theta}},\max}$ and $\lambda_{\hat{\bm{H}}_{\bm{\Theta}},\max}$, 
    as in Panel (a).
    (c) The eigenvalues of $\bm{H}_{\bm{\Theta}}$ and $\hat{\bm{H}}_{\bm{\Theta}}$, sorted 
    with the eigenvalues of $\bm{H}_{\bm{\Theta}}$.
    In Panels (a)-(c), the results are shown 
    at the 5th iteration (top) and at the last iteration (bottom).
    (d) Error measured as $\left\|\protect\begin{bmatrix}
    \tvec(\bm{\Theta}^{t}-\bm{\Theta}^*)\\
    \tvec(\bm{\Psi}^{t}-\bm{\Psi}^*)
    \protect\end{bmatrix}\right\|_2$ over iterations.
    (e) $\|\hat{\bm{H}}-\bm{H}\|_F$ over iterations.
    A dataset simulated from a random graph with $p, q=50$ was used.
}
    \label{fig:50bdemp}
\end{figure*}


Next, we compared the computation time of TeraLasso and EiGLasso on data simulated from larger graphs.
For datasets with size $p,q=100, 200, 500, 1000, 2000$, and $5000$, we ran TeraLasso and EiGLasso with
varying $K$ and recorded the time taken by each method if it converged within 24 hours
(Table \ref{tab:simresult}).
Across the different graph types and sizes, EiGLasso with different $K$'s almost always converged 
faster than TeraLasso. In particular, EiGLasso with $K=1$ and $2$ was 
consistently two to three orders-of-magnitude faster than
TeraLasso. For example, EiGLasso with $K=1$ took 3 hours for graphs with $p,q=5000$, 
while TeraLasso did not converge on the smaller graphs with $p,q=500$ in 24 hours.






We empirically evaluated the effects of Hessian approximation on the convergence behavior of EiGLasso,
using data simulated from a random graph of size $p,q=50$. 
In $\hat{\bm{H}}_{\bm{\Theta}}$ and $\hat{\bm{H}}_{\bm{\Psi}}$ of the approximate Hessian $\hat{\bm{H}}$, 
the components of the eigenvalues 
of $\bm{H}_{\bm{\Theta}}$ and $\bm{H}_{\bm{\Psi}}$ are modified
from the summands in Eq. \eqref{eq:Hcomp} to
the summands in Eq. \eqref{eqn:approxHL}. 
For EiGLasso with $K=5$,
these summands in the eigenvalues before and after the modification
are shown for the minimum and maximum eigenvalues of $\bm{H}_{\bm{\Theta}}$
at the $5$th and last iterations
in Figures \ref{fig:50rndemp}(a) and \ref{fig:50rndemp}(b).
The full spectrums of the same $\bm{H}_{\bm{\Theta}}$ and $\hat{\bm{H}}_{\bm{\Theta}}$
are shown in Figure \ref{fig:50rndemp}(c).
Next, we examined the effects of such Hessian approximation on the convergence 
of EiGLasso (Figures \ref{fig:50rndemp}(d) and \ref{fig:50rndemp}(e)).
While EiGLasso with the exact Hessian converges quadratically,
with the approximate Hessian, 
as the degree of Hessian approximation increases with smaller $K$, 
the convergence of EiGLasso slows down. However,
even with $K=1$, the convergence
is not significantly slower than EiGLasso with the exact Hessian. 
Even if $\|\hat{\bm{H}}-\bm{H}\|_F$ does not approach zero over iterations 
(Figure \ref{fig:50rndemp}(e)), EiGLasso with the approximate Hessian
quickly reaches the same level of optimality as EiGLasso with the exact Hessian. 
We obtained similar results from a graph with clusters 
(Figure \ref{fig:50bdemp}).




\subsection{Mouse Gene Expression Data from RNA-seq}

We compared EiGLasso and TeraLasso using the gene expression levels obtained from RNA sequencing of brain tissues 
for multiple related mice  from the same pedigree \citep{Gonzales2018}. 
The brain tissue samples were collected from generations 50-56 of the LG/J 
$\times$ SM/J advanced intercross line of mice. RNA was extracted 
from three mouse brain tissue types: prefrontal cortex (PFC; 185 mice), striatum (STR; 169 mice), and hippocampus (HIP; 208 mice). 
Several genes from these tissues have been found relevant to psychiatric and  metabolic disorders \citep{Gonzales2018}.  

We analyzed this data with EiGLasso and TeraLasso to estimate $\bm{\Theta}$
for the relationship among the $p$ mice and $\bm{\Psi}$ for the network over $q$ genes.
We selected $10,000$ genes from each tissue type by excluding the genes with low variance in expression 
across the mice. 
As TeraLasso did not converge within 24 hours on this dataset with $q=10,000$, 
to compare the computation time of EiGLasso and TeraLasso,
we formed three smaller datasets for each tissue type 
by performing hierarchical clustering and selecting clusters of genes.
We selected the regularization parameters from the range $[0.1, 1.0]$ 
using Bayesian Information Criterion (BIC). 
EiGLasso was run with $K=1$ with the convergence criterion $\epsilon = 10^{-3}$. 
TeraLasso was run up to 24 hours, until it reached similar objective values to those from EiGLasso, 
which corresponded to $\epsilon = 10^{-6}$ or $10^{-7}$. 

\begin{table}[t!]
\centering
\begin{tabular}{|c|r|r|r|}
    \hline
    Tissue Type & Genes ($q$) & EiGLasso (sec) & TeraLasso (sec) \\
    \hline
    \multicolumn{1}{|c|}{\multirow{4}{*}{\shortstack{PFC\\($p=185$)}}} & 436 & 2 & 887\\
    & 877 & 13 & 18194\\
    & 5108 & 6141 &\\
    & 10000 & 25360 &\\
    \hline
    \multicolumn{1}{|c|}{\multirow{4}{*}{\shortstack{STR\\($p=169$)}}} & 447 & 9 & 6917\\
    & 977 & 77 & 39645\\
    & 5007 & 32977 &\\
    & 10000 & 65797 &\\
    \hline
    \multicolumn{1}{|c|}{\multirow{4}{*}{\shortstack{HIP\\($p=208$)}}} & 508 & 90 & 54454\\
    & 879 & 708 & \\
    & 5529 & 3106 &\\
    & 10000 & 20104 &\\
    \hline
\end{tabular}
\caption{Comparison of the computation time of EiGLasso with $K=1$ and TeraLasso on mouse gene-expression data.}
\label{tab:realruntime}
\end{table}


\begin{figure*}[t!]
\centering
\setlength\tabcolsep{0.1pt}
\begin{tabular}{ccc}
    \begin{subfigure}{0.33\linewidth}
        \includegraphics[width=\linewidth]{plot_real2/pfc1000_logiter.pdf}
    \end{subfigure} &
    \begin{subfigure}{0.33\linewidth}
        \includegraphics[width=\linewidth]{plot_real2/str1000_logiter.pdf}
    \end{subfigure} &
    \begin{subfigure}{0.33\linewidth}
        \includegraphics[width=\linewidth]{plot_real2/hip1000_logiter.pdf}
    \end{subfigure}
    \\
    \begin{subfigure}{0.33\linewidth}
        \includegraphics[width=\linewidth]{plot_real2/pfc1000_logtime.pdf}
        \caption{}
    \end{subfigure} &
    \begin{subfigure}{0.33\linewidth}
        \includegraphics[width=\linewidth]{plot_real2/str1000_logtime.pdf}
        \caption{}
    \end{subfigure} &
    \begin{subfigure}{0.33\linewidth}
        \includegraphics[width=\linewidth]{plot_real2/hip1000_logtime.pdf}
        \caption{}
    \end{subfigure}
\end{tabular}
\caption{Comparison of the convergence of EiGLasso and TeraLasso on mouse gene-expression data. 
(a) PFC with $q=877$ genes. 
(b) STR with $q=977$ genes. 
(c) HIP with $q=879$ genes.
The objective values over iterations (top) and over time (bottom) are shown for 
each of the three brain tissue types. 
The `$\bullet$' and `$\bigcirc$' mark the points that satisfy the convergence criteria $\epsilon=10^{-3}$ and $10^{-7}$, respectively. 
EiGLasso with $K=1$ was used. 
}
    \label{fig:realconvg}
\end{figure*}


\begin{figure*}[t!]
\centering
\setlength\tabcolsep{0.2pt}
\begin{tabular}{ccc}
    \begin{subfigure}{0.33\linewidth}
        \includegraphics[width=\linewidth]{plot_real2/pfc1000_logcv.pdf}
        \caption*{(a)}
    \end{subfigure} &
    \begin{subfigure}{0.33\linewidth}
        \includegraphics[width=\linewidth]{plot_real2/str1000_logcv.pdf}
        \caption*{(b)}
    \end{subfigure} &
    \begin{subfigure}{0.33\linewidth}
        \includegraphics[width=\linewidth]{plot_real2/hip1000_logcv.pdf}
        \caption*{(c)}
    \end{subfigure}
\end{tabular}
\caption{Convergence of EiGLasso on mouse gene-expression data for different regularization parameters. 
For ten different values of the regularization parameter $\gamma=\gamma_{\bm{\Theta}}=\gamma_{\bm{\Psi}}$, 
the objective values over iterations are shown for  (a) PFC, (b) STR, and (c) HIP. 
EiGLasso with $K=1$ was run on the data for the same set of genes used in Figure \ref{fig:realconvg}.
}
    \label{fig:real1000cv}
\end{figure*}





Across all datasets with different tissue types and with different sizes,
EiGLasso was significantly faster than TeraLasso (Table \ref{tab:realruntime}).
On smaller datasets  with less than 1,000 genes,
EiGLasso was two to three orders-of-magnitude faster than TeraLasso. 
On HIP tissue with 879 genes, within 24 hours, TeraLasso was not able to converge to the similar objective value 
that EiGLasso reached, while EiGLasso converged in about 12 minutes. 
On the larger datasets with more than 5,000 genes, TeraLasso was not able to reach
the same objective as EiGLasso with $K=1$ within 24 hours.
On the full dataset with 10,000 genes, EiGLasso took 6 and 7 hours for the HIP and PFC tissue types
and 18 hours for the STR tissue type.


On the second largest dataset from each tissue type in Table \ref{tab:realruntime},
we examined the objective values over iterations and over time for EiGLasso and TeraLasso 
(Figure \ref{fig:realconvg}).
In all three datasets, EiGLasso converged two orders-of-magnitude faster than TeraLasso 
with much fewer iterations. 
To reach the similar objective value as EiGLasso, TeraLasso needed stricter convergence criterion 
$\epsilon=10^{-6}$ or $10^{-7}$. 
The number of iterations required by EiGLasso was not affected significantly
by different regularization  parameters $\gamma$ (Figure \ref{fig:real1000cv}).






\subsection{Financial data}

We applied EiGLasso and TeraLasso to the historical daily stock price data of 
the S\&P 500 constituents
to model the relationship among companies and dependencies across time points.
We obtained the daily stock closing prices for 306 companies 
that remained in the index from 2/16/2010 to 2/13/2020 for 2,523 days.
We excluded the data for six days when the stocks for only a small subset of the constituents were traded. 
We normalized the data by computing the proportion of price change on day $t$
from the previous day as $\frac{\textrm{price}^t - \textrm{price}^{t-1}}{\textrm{price}^{t-1}}$. 
We ran EiGLasso with $K=1$ and $\epsilon=10^{-3}$ and ran TeraLasso until
it reached the similar objective as EiGLasso.
We selected the optimal regularization parameter from
10 different values within the range of $[0.01, 1.0]$ using BIC.


\begin{figure}[t!]
\centering
\begin{tabular}{cc}
\multirow{2}{*}[2.5em]{
        \begin{subfigure}{.45\linewidth}
            \includegraphics[width=.888\linewidth]{plot_real2/Copy_of_untitled.png}
        \end{subfigure}
        } &
        \begin{subfigure}{.33\linewidth}
            \includegraphics[width=\linewidth]{plot_real2/snp500_logiter.pdf}
        \end{subfigure}\\
        &
        \begin{subfigure}{.33\linewidth}
            \includegraphics[width=\linewidth]{plot_real2/snp500_logtime.pdf}
        \end{subfigure}
\end{tabular}
\caption*{\hspace{3em}(a)\hspace{17em}(b)}
\caption{
Results on S\&P 500 data. (a) Graph over S\&P 500 companies estimated by EiGLasso from the stock closing price data.
Data from 306 companies over 2,516 days were used. 
The off-diagonals of the estimated graph are shown. (b) Comparison of the convergence of EiGLasso and TeraLasso on the data from 306 companies over 500 days. The objective values over iterations (top) and time (bottom) are shown. The `$\bullet$' and `$\bigcirc$' mark the points that satisfy 
the convergence criteria $\epsilon=10^{-3}$ and $10^{-7}$, respectively. 
EiGLasso with $K=1$ was used. 
}
    \label{fig:realrecover}
\end{figure}

EiGLasso with $K=1$ took 46 hours, whereas TeraLasso did not 
converge within 48 hours.
The graph over 306 companies estimated by EiGLasso
match the known sectors (Figure \ref{fig:realrecover}(a)).
We compared the computation time of EiGLasso with $K=1$ and TeraLasso
on a small subset over 500 days for the same 306 companies.
On this smaller dataset, EiGLasso took 148 seconds to converge and
TeraLasso took 25,739 seconds, showing that EiGLasso was again two orders-of-magnitude faster than TeraLasso (Figure \ref{fig:realrecover}(b)). 







\section{Conclusion} \label{sec:conclusion}

We introduced EiGLasso, an efficient algorithm for estimating an inverse covariance matrix as
the Kronecker sum of two matrices, each representing a graph over samples and a graph over features. 
Extending the second-order optimization method in QUIC for Gaussian graphical models~\citep{quic}, 
EiGLasso employed the eigendecomposition of the two matrices for individual graphs
during optimization, to avoid an expensive computation of large gradient and Hessian 
matrices with an inflated structure that arises from the Kronecker sum.
EiGLasso approximated the Hessian based on this eigenstructure of the parameters to further 
reduce the computation time. We showed the quadratic convergence of EiGLasso with the exact Hessian
and linear convergence of EiGLasso with the approximate Hessian. In our experiments,
EiGLasso achieved two to three orders-of-magnitude speed-up
over the state-of-the-art method, TeraLasso~\citep{TeraLasso}.

In addition, we introduced a new approach for identifying and estimating
the unidentifiable diagonal elements of the matrices for the feature and sample graphs.
We showed that given the ratio of the traces of the two matrices, the diagonal parameters
can be identified uniquely within the equivalence class of the unidentifiable parameters. 
EiGLasso performed the optimization without identifying 
these parameters, since all quantities involved in the optimization are
invariant within the equivalence class regardless of the trace ratio. 
We further showed that it is sufficient to identify the parameters 
with the fixed trace ratio once,  after the optimization was complete. 
Thus, EiGLasso used a significantly simpler strategy for estimating the
unidentifiable parameters than the previous methods.

There are several possible future directions. One natural extension of EiGLasso
is to adopt Big\&QUIC~\citep{bigquic}, an extension of QUIC to remove the memory 
requirement, within EiGLasso.
Because unlike QUIC, EiGLasso uses the eigendecomposition of the parameter
matrices, which in itself has the space complexity cubic in the graph sizes, 
further research is needed to incorporate the eigendecomposition within
the block-wise update used in Big\&QUIC.  
Another future direction is to develop a Hessian approximation technique with faster
theoretical convergence. 
While we showed that EiGLasso with our approximate Hessian 
converges linearly to the optimum, it may be possible to achieve superlinear convergence
with an approximation Hessian that approaches the exact Hessian over iterations.




\section*{Acknowledgements}
This work was supported by NIH 1R21HG011116 and 1R21HG010948. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by NSF ACI-1548562 and ACI-1445606.

\vskip 0.2in
