\documentclass[journal,11pt,draftcls,onecolumn]{IEEEtran}
%\documentclass[journal,10pt,twocolumn,twoside]{IEEEtran}

\usepackage{graphicx}
\usepackage{amsmath, amsthm, amssymb, amsfonts, cancel}
\usepackage{algorithm, algorithmic, ifsym, subfigure}


% Alter some LaTeX defaults for better treatment of figures:
    % See p.105 of "TeX Unbound" for suggested values.
    % See pp. 199-200 of Lamport's "LaTeX" book for details.
    %   General parameters, for ALL pages:
    \renewcommand{\topfraction}{0.9}	% max fraction of floats at top
    \renewcommand{\bottomfraction}{0.8}	% max fraction of floats at bottom
    %   Parameters for TEXT pages (not float pages):
    \setcounter{topnumber}{2}
    \setcounter{bottomnumber}{2}
    \setcounter{totalnumber}{4}     % 2 may work better
    \setcounter{dbltopnumber}{2}    % for 2-column pages
    \renewcommand{\dbltopfraction}{0.9}	% fit big float above 2-col. text
    \renewcommand{\textfraction}{0.07}	% allow minimal text w. figs
    %   Parameters for FLOAT pages (not text pages):
    \renewcommand{\floatpagefraction}{0.7}	% require fuller float pages
	% N.B.: floatpagefraction MUST be less than topfraction !!
    \renewcommand{\dblfloatpagefraction}{0.7}	% require fuller float pages

	% remember to use [htp] or [htpb] for placement

%\def\BibTeX{{\rm B\kern-.05em{\sc i\kern-.025em b}\kern-.08em
%    T\kern-.1667em\lower.7ex\hbox{E}\kern-.125emX}}
\def\Var{{\rm Var}}
\def\Cov{{\rm Cov}}
\def\tr{{\rm tr}}
\def\E{{\rm E}}
\def\P{{\rm P}}
\def\bA{ {\mathbf{A}} }
\def\bB{ {\mathbf{B}} }
\def\bC{ {\mathbf{C}} }
\def\bD{ {\mathbf{D}} }
\def\bE{ {\mathbf{E}} }
\def\bF{ {\mathbf{F}} }
\def\bG{ {\mathbf{G}} }
\def\bH{ {\mathbf{H}} }
\def\bI{ {\mathbf{I}} }
\def\bJ{ {\mathbf{J}} }
\def\bK{ {\mathbf{K}} }
\def\bL{ {\mathbf{L}} }
\def\bM{ {\mathbf{M}} }
\def\bN{ {\mathbf{N}} }
\def\bO{ {\mathbf{O}} }
\def\bP{ {\mathbf{P}} }
\def\bQ{ {\mathbf{Q}} }
\def\bR{ {\mathbf{R}} }
\def\bS{ {\mathbf{S}} }
\def\bT{ {\mathbf{T}} }
\def\bU{ {\mathbf{U}} }
\def\bV{ {\mathbf{V}} }
\def\bW{ {\mathbf{W}} }
\def\bX{ {\mathbf{X}} }
\def\bY{ {\mathbf{Y}} }
\def\bZ{ {\mathbf{Z}} }

\def\bTheta{ {\mathbf{\Theta}} }
\def\bSigma{ {\mathbf{\Sigma}} }
\def\bDelta{ {\mathbf{\Delta}} }
\def\bXi{ {\mathbf{\Xi}} }

\def\nn{{ \parallel   }}
\def\RR{{ \mathbb{R}  }}
\def\PP{{ \mathbb{P}  }}
\def\NN{{ \mathbb{N}  }}
\def\vec{{ \text{vec} }}
\def\tr{{ \text{tr}   }}
\def\bdiag{{ \text{bdiag}   }}
\def\diag{{ \text{diag}   }}
\def\corr{{ \text{corr}   }}
\def\bx{{ \mathbf{x}  }}
\def\by{{ \mathbf{y}  }}
\def\bz{{ \mathbf{z}  }}
\def\bv{{ \mathbf{v}  }}
\def\card{{ \text{card} }}
\def\erf{{ \text{erf} }}
\def\sgn{{ \text{sgn} }}

\newtheorem{example}{Example}
\newtheorem{theorem}{Theorem}
\newtheorem{lemma}{Lemma}
\newtheorem{remark}{Remark}
\newtheorem{corollary}{Corollary}
\newtheorem{proposition}{Proposition}
\newtheorem{assumption}{Assumption}
\newtheorem{definition}{Definition}
\newtheorem{conjecture}{Conjecture}
\newtheorem{property}{Property}

\setcounter{page}{1}

\begin{document}
\title{On Convergence of Kronecker Graphical Lasso Algorithms}
\author{Theodoros Tsiligkaridis *, \textit{Student Member, IEEE}, Alfred O. Hero III, \textit{Fellow, IEEE}, Shuheng Zhou, \textit{Member, IEEE}}

\maketitle

\begin{abstract}
This paper presents a thorough convergence analysis of 
Kronecker graphical lasso (KGLasso) algorithms for estimating the covariance of an i.i.d. Gaussian random sample under a sparse Kronecker-product covariance model. The KGlasso model, originally called the transposable regularized covariance model by Allen {\it et al}  \cite{AllenTib10}, implements a pair of $\ell_1$ penalties on each Kronecker factor to enforce sparsity in the covariance estimator.   The KGlasso algorithm generalizes Glasso, introduced by Yuan and Lin \cite{YL07} and Banerjee {\it et al} \cite{ModelSel}, to estimate covariances having Kronecker product form.   It  also generalizes  the unpenalized ML flip-flop (FF) algorithm of Werner {\it et al} \cite{EstCovMatKron} to estimation of sparse Kronecker factors. %We establish that the KGlasso iterates converge pointwise to a local maximum of the penalized likelihood function.
We establish high dimensional rates of convergence to the true covariance as both the number of samples and the number of variables go to infinity. Our results establish that KGlasso has significantly faster asymptotic convergence than Glasso and FF. Simulations are presented that validate the results of our analysis.  For example, for a sparse $10,000 \times 10,000$ covariance matrix equal to the Kronecker product of two $100 \times 100$ matrices, the root mean squared error of the inverse covariance estimate using FF is 3.5 times larger than that obtainable using KGlasso.
\end{abstract} 

\begin{keywords}
	Sparsity, structured covariance estimation, penalized maximum likelihood, graphical lasso, direct product representation.
\end{keywords}


\section{Introduction}
Covariance estimation is a problem of great interest in many different disciplines, including machine learning, signal processing, economics and bioinformatics. In many applications the number of variables is very large, e.g., in the tens or hundreds of thousands, leading to a number of covariance parameters that greatly exceeds the number of observations.  To address this problem constraints are frequently imposed on the covariance to reduce the number of parameters in the model. For example, the Glasso model of  Yuan and Lin \cite{YL07} and Banerjee {\it et al} \cite{ModelSel}  imposes sparsity constraints on the covariance. The Kronecker product  model of Werner {\it et al} \cite{EstCovMatKron} assumes that the covariance can be represented as the Kronecker product  of two lower dimensional covariance matrices. The 
transposable regularized covariance model of Allen {\it et al}  \cite{AllenTib10} imposes a combination of sparsity and Kronecker product form on the covariance. When there is no missing data, an extension of the alternating optimization algorithm of  \cite{EstCovMatKron}, called the flip flop (FF) algorithm, can be applied to estimate the parameters of this combined sparse and Kronecker product model.  In this paper we call this algorithm the Kronecker Glasso (KGlasso) and  we thoroughly analyze convergence of the algorithm in the high dimensional setting.  


As in \cite{EstCovMatKron} we assume that there are $pf$ variables whose covariance $\mathbf{\Sigma}_0$ has the separable positive definite Kronecker product representation:
\begin{equation} \label{factorize}
	\mathbf{\Sigma}_{0} = \mathbf{A}_0 \otimes \mathbf{B}_0
\end{equation}
where $\bA_0$ is a $p\times p$ positive definite matrix and $\bB_0$ is an $f \times f$ positive definite matrix.  This model (\ref{factorize}) is relevant to channel modeling for MIMO wireless communications, where $\bA_0$ is a transmit covariance matrix and $\bB_0$ is a receive covariance matrix \cite{MIMOWerner}. The model is also relevant to other transposable models arising in recommendation systems like NetFlix and in gene expression analysis \cite{AllenTib10}. 
The Kronecker factorization (\ref{factorize}) can easily be generalized to the $k$-fold case, where $\mathbf{\Sigma}_0=\mathbf{A}_1\otimes \mathbf{A}_2 \otimes \dots \otimes \mathbf{A}_k$. 

Under the assumption that the measurements are multivariate Gaussian with covariance having the Kronecker product form (\ref{factorize}), the maximum likelihood (ML) estimator can be formulated \cite{LuZimmerman2}. While the ML estimator has no known closed-form solution, an approximation to the solution can be iteratively computed via an alternating algorithm: the flip-flop (FF) algorithm \cite{LuZimmerman2, EstCovMatKron}.  As compared to the standard saturated (unstructured) covariance model, the number of unknown parameters in (\ref{factorize}) is reduced from order $\Theta(p^2 f^2)$ to order $\Theta(p^2)+\Theta(f^2)$.
This results in a significant reduction in the mean squared error (MSE) and the computational complexity of the maximum likelihood (ML) covariance estimator. This paper establishes that further reductions MSE are achievable when the Kronecker matrix factors are known to have sparse inverses, i.e., the measurements obey a sparse Kronecker structured Gaussian graphical model.

The graphical lasso (Glasso) estimator was originally proposed in \cite{YL07, ModelSel} for estimating a sparse inverse covariance, also called the precision matrix, under an i.i.d. Gaussian observation model. An algorithm for efficiently solving the nonsmooth optimization problem that arises in the Glasso estimator, based on ideas from \cite{ModelSel}, was proposed in \cite{Glasso}. Glasso has been applied to the time-varying coefficients setting in Zhou {\it et al} \cite{TimeVaryingGraphs} using the kernel estimator for covariances at a target time. Rothman {\it et al} \cite{Rothman} derived high dimensional convergence rates for a slight variant of Glasso, i.e., only the off-diagonal entries of the estimated precision matrix were penalized using an $\ell_1$-penalty. The high dimensional convergence rate of Glasso was established by Ravikumar {\it et al}  \cite{RWRY08}. This paper extends their analysis to the case that the covariance has Kronecker structure (\ref{factorize}), showing that significantly higher rates of convergence are achievable.


%The first contribution of this paper is proof that the KGlasso iterations converge to a local maximum of the KGlasso objective function for fixed number of samples $n$ and fixed matrix dimensions $p$ and $f$.
%
The main contribution is the derivation of the high-dimensional MSE convergence rates for KGlasso as $n$, $p$ and $f$ go to infinity.
 When both Kronecker factors are sparse, it is shown that KGlasso \textit{strictly} outperforms FF and Glasso in terms of MSE convergence rate. More specifically, we show KGlasso achieves a convergence rate of $O_P\left(\frac{(p+f)\log \max(p,f,n)}{n}\right)$ and FF achieves a rate of $O_P\left(\frac{(p^2+f^2)\log \max(p,f,n)}{n}\right)$ as $n\to\infty$, while it is known \cite{Rothman, TimeVaryingGraphs} that Glasso achieves a rate of $O_P\left(\frac{(pf+s)\log \max(p,f,n)}{n}\right)$, where $s$ denotes the number of off-diagonal nonzero elements in the true precision matrix $\bTheta_0$. Simulations show that the performance improvements predicted by the high-dimensional analysis continue to hold for small sample size and moderate matrix dimension. For the example studied in Sec. \ref{sec: simulations} the empirical MSE of KGlasso is significantly lower than that of Glasso and FF for $p=f=100$ over the range of $n$  from $10$ to $100$.


The starting point for the MSE convergence analysis is the large-sample analysis of the FF algorithm (Thm. 1 in \cite{EstCovMatKron}). The KGlasso convergence proof uses a large deviation inequality  that shows that the dimension of one estimated Kronecker factor, say $\bA$, acts as a multiplier on the number of independent samples when performing inference on the other factor $\bB$. This result is then used to obtain optimal MSE rates in terms of Frobenius norm error between the KGlasso estimated matrix and the ground truth. The asymptotic MSE convergence analysis is useful since it can be used to guide the selection of sparsity regularization parameters and to determine minimum sample size requirements.

The outline of the paper is as follows. Section \ref{sec: notation} introduces the notation that will be used throughout the paper. In Section \ref{sec: glasso}, the graphical lasso framework is introduced. Section \ref{sec: algorithm} uses this framework to introduce the KGlasso algorithm. Section \ref{sec: convergence_analysis_KGL} shows convergence of KGlasso and characterizes its limit points. The high dimensional MSE convergence rate derivation for the FF algorithm is included in Section \ref{sec: FF_rate}. Section \ref{sec: consistency} presents a high-dimensional MSE rate result that is used to establish the superiority of KGlasso as compared to FF and standard Glasso, under the sparse Kronecker product representation (\ref{factorize}). Section \ref{sec: simulations} presents simulations that empirically validate the theoretical convergence rates obtained in Section \ref{sec: consistency}.


\section{Notation} \label{sec: notation}
For a square matrix $\mathbf{M}$, define $|\mathbf{M}|_1=\nn \vec(\mathbf{M}) \nn_1$ and $|\mathbf{M}|_\infty=\nn \vec(\mathbf{M}) \nn_\infty$, where $\vec(\mathbf{M})$ denotes the vectorized form of $\mathbf{M}$ (concatenation of columns into a vector). $\nn \mathbf{M} \nn_2$ is the spectral norm of $\bM$. $\mathbf{M}_{i,j}$ and $[\mathbf{M}]_{i,j}$ are the $(i,j)$th element of $\mathbf{M}$. Let the inverse transformation (from a vector to a matrix) be defined as: $\vec^{-1}(\mathbf{x})=\mathbf{X}$, where $\mathbf{x}=\vec(\mathbf{X})$. Define the $pf\times pf$ permutation operator $\bK_{p,f}$ such that $\bK_{p,f} \vec(\bN) = \vec(\bN^T)$ for any $p\times f$ matrix $\bN$. For a symmetric matrix $\mathbf{M}$, $\lambda(\mathbf{M})$ will denote the vector of real eigenvalues of $\mathbf{M}$ and define $\lambda_{max}(\mathbf{M})=\nn\bM \nn_2=\max{\lambda_i(\mathbf{M})}$ for p.d. symmetric matrix, and $\lambda_{min}(\mathbf{M}) = \min{\lambda_i(\mathbf{M})}$. Define the sparsity parameter associated with $\bM$ as $s_{M}=\card(\{(i_1,i_2): [\bM]_{i_1,i_2}\neq 0, i_1\neq i_2 \})$. Let $\kappa(\bM):=\frac{\lambda_{max}(\bM)}{\lambda_{min}(\bM)}$ denote the condition number of a symmetric matrix $\bM$.

For a matrix $\bM$ of size $pf\times pf$, let $\{\bM(i,j)\}_{i,j=1}^p$ denote its $f\times f$ block submatrices, where each block submatrix is $\bM(i,j)=[\bM]_{(i-1)f+1:if,(j-1)f+1:jf}$. Also let $\{\overline{\bM}(k,l)\}_{k,l=1}^f$ denote the $p\times p$ block submatrices of the permuted matrix $\overline{\bM}=\bK_{p,f}^T \bM \bK_{p,f}$.

Define the set of symmetric matrices $S^p = \{\mathbf{A}\in \RR^{p\times p}: \mathbf{A}=\mathbf{A}^T\}$, the set of symmetric positive semidefinite (psd) matrices $S_{+}^p$, and the set of symmetric positive definite (pd) matrices $S_{++}^p$. $\bI_d$ is a $d\times d$ identity matrix. It can be shown that $S_{++}^p$ is a convex set, but is not closed \cite{ConvexOpt}. Note that $S_{++}^p$ is simply the interior of the closed convex cone $S_{+}^p$.

Statistical convergence rates will be denoted by the $O_P(\cdot)$ notation, which is defined as follows. Consider a sequence of real random variables $\{X_n\}_{n \in \NN}$ defined on a probability space $(\Omega,\mathcal{F},P)$ and a deterministic (positive) sequence of reals $\{b_n\}_{n\in \NN}$. By $X_n=O_P(1)$ is meant: $\sup_{n\in \NN}{\Pr(|X_n|>K)} \to 0$ as $K\to\infty$, where $\bX_n$ is a sequence indexed by $n$, for fixed $p,f$. %for all $\epsilon>0$, $\exists M=M(\epsilon)<\infty$ such that $\sup_{n\in \NN}{P(|X_n|>M)} < \epsilon$.
The notation $X_n=O_P(b_n)$ is equivalent to $\frac{X_n}{b_n}=O_P(1)$. By $X_n=o_p(1)$ is meant $\Pr(|X_n|>\epsilon) \to 0$ as $n\to\infty$ for any $\epsilon>0$. By $\lambda_n \asymp b_n$ is meant $c_1 \leq \frac{\lambda_n}{b_n} \leq c_2$ for all $n$, where $c_1,c_2>0$ are absolute constants.



\section{Graphical Lasso Framework} \label{sec: glasso}
For simplicity, we assume the number of Kronecker components is $k=2$. Available are $n$ i.i.d. multivariate Gaussian observations $\{\bz_t\}_{t=1}^{n}$, where $\bz_t\in \RR^{pf}$, having zero-mean and covariance equal to $\mathbf{\Sigma} = \mathbf{A}_0 \otimes \mathbf{B}_0$. Then, ignoring irrelevant constants,  the log-likelihood  $l(\mathbf{\Sigma})$ is:
\begin{equation} \label{log_likelihood}
	l(\mathbf{\Sigma}) = \log\det(\mathbf{\Sigma}^{-1}) - \tr(\mathbf{\Sigma}^{-1}\hat{\mathbf{S}}_n),
\end{equation}
where $\mathbf{\Sigma}$ is the positive definite covariance matrix and $\hat{\mathbf{S}}_n=\frac{1}{n}\sum_{t=1}^{n}{\bz_t \bz_t^T}$ is the sample covariance matrix. Recent work \cite{ModelSel, Glasso, RWRY08} has considered $\ell_1$-penalized maximum likelihood estimators for the saturated model where $\bSigma$ belongs to the unrestricted cone of positive definite matrices. These estimators are known as graphical lasso (Glasso) estimators and are obtained as the solution to the $\ell_1$-penalized minimization problem:
\begin{equation} \label{opt_prob}
	\hat{\mathbf{\Sigma}}_n \in \arg \min_{\mathbf{\Sigma} \in S_{++}^p} { \{ -l(\mathbf{\Sigma}) + \lambda |\mathbf{\Sigma}^{-1}|_1 \} },
\end{equation}
where $\lambda \geq 0$ is a regularization parameter. If $\lambda>0$ and $\hat{\mathbf{S}}_n$ is positive definite, then $\hat{\mathbf{\Sigma}}_n$ in (\ref{opt_prob}) is the unique minimizer.

A fast iterative algorithm, based on a block coordinate descent approach, exhibiting a computational complexity $\mathcal{O}((pf)^3)$, was developed in \cite{Glasso} to solve the convex program (\ref{opt_prob}). Under the assumption $\lambda \asymp \sqrt{\frac{\log(pf)}{n}}$ solution of (\ref{opt_prob}) was shown to have high dimensional convergence rate \cite{Rothman, TimeVaryingGraphs}:
\begin{equation} \label{Glasso_rate}
 \nn \bG(\hat{\bS}_n, \lambda)-\bTheta_0 \nn_F = O_P\left( \sqrt{\frac{(pf+s) \log(pf)}{n}} \right)
\end{equation}
where $s$ is an upper bound on the number of non-zero off-diagonal elements of $\bTheta_0$. When $s=O(pf)$, this rate is better than that achieved in the case of the standard sample covariance estimator ($\lambda=0$):
\begin{equation} \label{naive_SCM_rate}
	\nn \hat{\bS}_n - \bSigma_0 \nn_F = O_P\left( \sqrt{\frac{p^2f^2}{n}} \right).
\end{equation}



\section{Kronecker Graphical Lasso} \label{sec: algorithm}
Let $\mathbf{\Sigma}_0=\mathbf{A}_0 \otimes \mathbf{B}_0$ denote the true covariance matrix, where $\mathbf{A}_0=\mathbf{X}_0^{-1}$ and $\mathbf{B}_0=\mathbf{Y}_0^{-1}$ are the true Kronecker factors. Let $\mathbf{A}_{init}$ denote an initial guess of $\mathbf{A}_0=\mathbf{X}_0^{-1}$.

Define $J(\bX,\bY)$ as the negative log-likelihood
\begin{align}
	J(\bX,\bY) &= \tr((\bX \otimes \bY)\hat{\bS}_n) - f \log\det(\bX) \nonumber \\
		& \quad -p \log\det(\bY) \label{J_func_2}
\end{align}
Although the objective function (\ref{J_func_2}) is not jointly convex in $(\bX, \bY)$, it is biconvex. This motivates the flip-flop algorithm \cite{EstCovMatKron}. 
%
Adapting the notation from \cite{EstCovMatKron}, define the mappings $\hat{\mathbf{A}}(\cdot), \hat{\mathbf{B}}(\cdot)$:
\begin{align}
	\underbrace{\hat{\mathbf{A}}(\mathbf{B})}_{p\times p} &= \frac{1}{f} \sum_{k,l=1}^f{[\mathbf{B}^{-1}]_{k,l}\overline{\hat{\mathbf{S}}_n}(l,k)}, \label{A_update} \\
	\underbrace{\hat{\mathbf{B}}(\mathbf{A})}_{f\times f} &= \frac{1}{p} \sum_{i,j=1}^p{[\mathbf{A}^{-1}]_{i,j} \hat{\mathbf{S}}_n(j,i)}, \label{B_update}
\end{align}
where $\overline{\hat{\bS}}_n=\bK_{p,f}^T \hat{\bS}_n \bK_{p,f}$ (see Sec. \ref{sec: notation} for definition of $K_{p,f}$). For fixed $\bB\in S_{++}^f$, $\hat{\bA}(\bB)$ in (\ref{A_update}) is the minimizer of $J(\bA^{-1},\bB^{-1})$ over $\bA\in S_{++}^p$. A similar interpretation holds for (\ref{B_update}). The flip-flop algorithm starts with some arbitrary p.d. matrix $\bA_{init}$ and computes $\bB$ using (\ref{B_update}), then $\bA$ using (\ref{A_update}), and repeats until convergence. This algorithm does not account for sparsity.

If $\bTheta_0=\bX_0\otimes \bY_0$ is a sparse matrix, which implies that at least one of $\bX_0$ or $\bY_0$ is sparse, one can penalize the outputs of the flip-flop algorithm and minimize
\begin{equation} \label{J_lambda_func}
	J_\lambda(\bX,\bY) = J(\bX,\bY) + \bar{\lambda}_X |\bX|_1 + \bar{\lambda}_Y |\bY|_1 .
\end{equation}
where $\bar{\lambda}_X=\lambda_X/f$ and $\bar{\lambda}_Y=\lambda_Y/p$. This leads to an algorithm that we call KGlasso (see Algorithm \ref{alg: algKGL}), which sparsifies the Kronecker factors in proportion to the parameters $\bar{\lambda}_X, \bar{\lambda}_Y >0$. This is the same objective function that was proposed in \cite{AllenTib10} when specialized to the case that there is no missing data.

\begin{algorithm}
\caption{Kronecker Graphical Lasso (KGlasso) }
\label{alg: algKGL}
\begin{algorithmic}[1]
\STATE \textbf{Input:}  {$\hat{\bS}_n$, $p$, $f$, $n$, $\bar{\lambda}_X >0$, $\bar{\lambda}_Y >0$}
\STATE \textbf{Output:} {$\hat{\bTheta}_{KGlasso}$}

	\STATE Initialize $\bA_{init}$ to be positive definite satisfying Assumption \ref{assumption_posdef_unif}.
	\STATE $\check{\bX} \leftarrow \bA_{init}^{-1}$
	
	\REPEAT
	  \STATE	$\hat{\bB} \leftarrow \frac{1}{p} \sum_{i,j=1}^p{[\check{\bX}]_{i,j} \hat{\bS}_n(j,i)}$ (see Eq. (\ref{A_update}))
  	\STATE  $\check{\bY} \leftarrow \bG(\hat{\bB}, \frac{\bar{\lambda}_Y}{p})$, where $\bG(\cdot,\cdot)$ is defined in (\ref{G_operator})
  	\STATE  $\hat{\bA} \leftarrow \frac{1}{f} \sum_{k,l=1}^f{[\check{\bY}]_{k,l} \overline{\hat{\bS}_n}(l,k)}$ (see Eq. (\ref{B_update}))
  	\STATE  $\check{\bX} \leftarrow \bG(\hat{\bA}, \frac{\bar{\lambda}_X}{f})$
	\UNTIL {convergence}
	
	\STATE $\hat{\bTheta}_{KGlasso} \leftarrow \check{\bX} \otimes \check{\bY}$
\end{algorithmic}
\end{algorithm}


The Glasso mapping (\ref{opt_prob}) is written as $\bG(\cdot,\lambda): S^d \to S^d$,
\begin{equation} \label{G_operator}
	\bG(\bT,\lambda) = \arg\min_{\bTheta \in S_{++}^d} \Big\{ \tr(\bTheta \bT)-\log\det(\bTheta)+\lambda|\bTheta|_1 \Big\}.
\end{equation}
As compared to the $\mathcal O(p^3f^3)$ computational complexity of Glasso, KGlasso has a computational complexity of only $\mathcal{O}(p^3+f^3)$.

\section{Convergence of KGlasso Iterations} \label{sec: convergence_analysis_KGL}
In this section, we provide an alternative characterization of the KGlasso algorithm (Algorithm \ref{alg: algKGL}) and show the iterations converge pointwise to a critical point of the objective function. Under a mild condition on the starting point they will converge to a local minimum.


\subsection{Block-Coordinate Reformulation of KGlasso}

The following lemma shows that exploiting the property that the KGlasso algorithm is a block-coordinate optimization of the penalized objective function (\ref{J_lambda_func}), each subproblem takes the form of standard Glasso applied on a compressed version of the SCM that is relevant for inference in each step. %We show this equivalence using duality in the following lemma.

\begin{lemma} \label{dual_lemma}
The  KGlasso objective function (\ref{J_lambda_func}) has the following properties:
\begin{enumerate}
	\item Assume $\bar{\lambda}_X, \bar{\lambda}_Y \geq 0$ and $\bX\in S_{++}^p,\bY\in S_{++}^f$. When one argument of $J_\lambda(\bX,\bY)$ is fixed, the objective function (\ref{J_lambda_func}) is convex in the other argument.
	\item Assume $\hat{\bS}_n$ is positive definite. Consider $J_\lambda(\mathbf{X},\mathbf{Y})$ in (\ref{J_lambda_func}) with matrix $\mathbf{X}\in S_{++}^p$ fixed. Then, the dual subproblem for minimizing $J_\lambda(\mathbf{X},\mathbf{Y})$ over $\mathbf{Y}$ is:
\begin{equation} \label{dualY}
	\max_{ |\mathbf{W}-\frac{1}{p}\sum_{i,j=1}^p{\mathbf{X}_{i,j} \hat{\bS}_n(j,i)}|_{\infty} \leq \lambda_Y }{ \log \det(\mathbf{W})  }
\end{equation}
where $\lambda_Y := \bar{\lambda}_Y/p$.

On the other hand, consider (\ref{J_lambda_func}) with matrix $\mathbf{Y} \in S_{++}^f$ fixed. Then, the dual problem for minimizing $J_\lambda(\mathbf{X},\mathbf{Y})$ over $\bX$ is:
\begin{equation} \label{dualX}
	\max_{ |\mathbf{Z}-\frac{1}{f}\sum_{k,l=1}^f{\mathbf{Y}_{k,l} \overline{\hat{\bS}_n}(l,k)}|_{\infty}\leq \lambda_X  }{ \log \det(\mathbf{Z})  }
\end{equation}
where $\overline{\hat{\bS}_n}:=\bK_{p,f}^T \hat{\bS}_n \bK_{p,f}$ and $\lambda_X := \bar{\lambda}_X/f$.
	
	\item Strong duality holds for (\ref{dualY}) and (\ref{dualX}).
	\item The solutions to (\ref{dualY}) and (\ref{dualX}) are positive definite.
\end{enumerate}
\end{lemma}
\begin{IEEEproof}
	See Appendix.
\end{IEEEproof}


Since the dual subproblems (\ref{dualY}) and (\ref{dualX}) are maximizations of a strictly concave function over a closed convex set they have unique solution attaining the maximum.  Lemma \ref{dual_lemma} is similar to the result obtained in \cite{ModelSel}, but with $(\frac{1}{p}\sum_{i,j=1}^p{\bX_{i,j} \hat{\bS}_n(j,i)}, \lambda_Y)$ playing the role of $(\hat{\bS}_n,\lambda)$, for the fixed $\bX$ subproblem.


% Convergence analysis
\subsection{Limit Point Characterization of KGlasso}
The following theorem establishes that KGlasso converges to a local minimum of the penalized likelihood function (\ref{log_likelihood}).

\begin{theorem} \label{convergence_theorem_critical}
Assume $n>pf$. Then the KGlasso iterations converge to a critical point of the negative penalized likelihood function (\ref{J_lambda_func}).

Assuming $(\mathbf{X}^{(0)},\mathbf{Y}^{(0)})$ is not a local maximum, the KGlasso iterations converge to a local minimum of the negative penalized likelihood function (\ref{J_lambda_func}).
\end{theorem}
\begin{IEEEproof}
	See \cite{TsiligkaridisTSP}. Also see Prop. 2 in \cite{AllenTib10} for an alternative proof of the first part.
\end{IEEEproof}

The proof of Thm. \ref{convergence_theorem_critical} is built on several lemmas (see \cite{TsiligkaridisTSP}). The main line of argument is as follows. For $n>pf$, the SCM is positive definite a.s., which implies that the objective function is bounded below. This can be used to show that the iterates generated by Algorithm 1 converge to a fixed point. Combining this result with the KKT optimality conditions and the strict descent property of the algorithm, we arrive to the result in Thm. \ref{convergence_theorem_critical}.


\section{High Dimensional Consistency of FF} \label{sec: FF_rate}
In this section, we show that the flip-flop (FF) algorithm achieves the optimal (non-sparse) statistical convergence rate of $O_P \left(\sqrt{\frac{(p^2+f^2)\log M}{n}}\right)$. This result (see Thm. \ref{thm: FF_optimal_rate}) will be compared to the statistical convergence rate of KGlasso (see Thm. \ref{thm: KGL_optimal_rate}) to establish that KGlasso has lower asymptotic MSE than FF. We make the following boundedness assumptions on the spectra of the Kronecker factors.

\begin{assumption} \label{assumption_posdef_unif}
Uniformly Bounded Spectra \\
There exist absolute constants $\underline{k}_A, \overline{k}_A, \underline{k}_B, \overline{k}_B, \underline{k}_{A_{init}}, \overline{k}_{A_{init}}$ such that: \\
\indent 1a. $0<\underline{k}_A \leq \lambda_{min}(\bA_0) \leq \lambda_{max}(\bA_0) \leq \overline{k}_A < \infty$ \\
\indent 1b. $0<\underline{k}_B \leq \lambda_{min}(\bB_0) \leq \lambda_{max}(\bB_0) \leq \overline{k}_B < \infty$ \\
\indent 2. $0<\underline{k}_{A_{init}} \leq \lambda_{min}(\bA_{init}) \leq \lambda_{max}(\bA_{init}) \leq \overline{k}_{A_{init}} < \infty$
\end{assumption}

Let $\bSigma_{FF}(3):=\hat{\bA}(\hat{\bB}(\bA_{init})) \otimes \hat{\bB}(\hat{\bA}(\hat{\bB}(\bA_{init})))$ denote the 3-step (noniterative) version of the flip-flop algorithm \cite{EstCovMatKron}. 
%More generally, let $\bSigma_{FF}(k)$ denote the $k$-step version of the flip-flop algorithm, and denote its inverse as $\bTheta_{FF}(k)=(\bSigma_{FF}(k))^{-1}$.

\begin{theorem} \label{thm: FF_optimal_rate}
	Let $\bA_0,\bB_0$, and $\bA_{init}$ satisfy Assumption \ref{assumption_posdef_unif} and define $M=\max(p,f,n)$. Assume $p\geq f\geq 2$ and $p \log M \leq C'' n$ for some finite constant $C''>0$. Finally, assume $n \geq \frac{p}{f} + 1$.
	Then, %for $k\geq 2$ finite,
	\begin{equation} \label{FF_rate_2}
		\nn \bTheta_{FF}(3) - \bTheta_0 \nn_F = O_P\left( \sqrt{\frac{(p^2+f^2) \log M}{n}} \right)
	\end{equation}
	as $n\to\infty$.
\end{theorem}
\begin{IEEEproof}
	See Appendix.
\end{IEEEproof}
\begin{remark}
	The sufficient conditions are symmetric with respect to $p$ and $f$-i.e. for $f\geq p$, the corresponding conditions would become $f\log M\leq C'' n$ for some constant $C''>0$, and $n \geq \frac{f}{p} + 1$.
\end{remark}

For the special case of $p=f$, the sufficient conditions of Thm. \ref{thm: FF_optimal_rate} become $p\log M=O(n)$. The relation (\ref{FF_rate_2}) indicates that the error is asymptotically bounded as long as $n$ is of order $\Omega((p^2+f^2)\log M)$.
The relation (\ref{FF_rate_2}) specifies the rate of reduction of the estimation error for the three step FF algorithm ($k=3$) \cite{EstCovMatKron}. This relation will also hold for the multi-step FF as long as the number of steps are finite. Note that (\ref{FF_rate_2}) specifies a faster rate than that of the ordinary ML sample covariance matrix estimator (\ref{naive_SCM_rate}).
%Furthermore, since the computational complexity for FF is $\mathcal{O}(p^2+f^2)$ which is less than the $\mathcal{O}(p^2f^2)$ complexity of SCM, by exploiting Kronecker structure FF simultaneously achieves improved MSE performance and reduced computational complexity.


\section{High Dimensional Consistency of KGlasso} \label{sec: consistency}
Here a relation like (\ref{FF_rate_2}) is established for KGlasso. Recall that a $p\times p$ matrix is called sparse if its number of nonzero elements is of order $p$. Recall $\lambda_X=\bar{\lambda}_X f$ and $\lambda_Y=\bar{\lambda}_Y p$, as in (\ref{J_lambda_func}).

%Define $\bTheta_{KGlasso}(k)$ as the output of the $k$th compression and sparsification step (two of these steps constitute a full KGlasso iteration).

\begin{theorem} \label{thm: KGL_optimal_rate}
Assume $\bX_0$ and $\bY_0$ are sparse.
Let $\bA_0,\bB_0, \bA_{init}$ satisfy Assumptions \ref{assumption_posdef_unif}.  Let $M=\max(p,f,n)$.	Let $\lambda_Y^{(1)} \asymp \sqrt{\frac{\log M}{np}}$, and $\lambda_X^{(2)},\lambda_Y^{(3)} \asymp \left( \frac{1}{\sqrt{p}} + \frac{1}{\sqrt{f}} \right) \sqrt{\frac{\log M}{n}}$.
%as $p,f,n\to\infty$ for all $k \geq 1$ and $k' \geq 2$.  
Then, if  $\max\left( \frac{p}{f}, \frac{f}{p}\right) \log M=o(n)$, %for $k\geq 2$ finite, we have
	\begin{equation} \label{KGL_perfect_rate}
		\nn \bTheta_{KGlasso}(3) - \bTheta_0 \nn_F = O_P\left( \sqrt{\frac{(p+f) \log M}{n}} \right)
	\end{equation}
	as $n\to\infty$.
\end{theorem}
\begin{IEEEproof}
	See Appendix.
\end{IEEEproof}

Theorem \ref{thm: KGL_optimal_rate} offers a strict improvement over standard Glasso \cite{Rothman, ModelSel} and generalizes Thm. 1 in \cite{Rothman} to the case of sparse Kronecker product structure. Thm. \ref{thm: KGL_optimal_rate} generalizes Thm. \ref{thm: FF_optimal_rate} to the case of sparse Kronecker structure. %The computational complexity of KGlasso is $\mathcal{O}(p^3+f^3)$ which is less than the complexity of standard Glasso $\mathcal{O}(p^3 f^3)$ \cite{Glasso}. 
Comparison between the error expressions (\ref{Glasso_rate}), (\ref{FF_rate_2}) and (\ref{KGL_perfect_rate}) show that, by exploiting both Kronecker structure and sparsity, KGlasso can attain significantly lower estimation error than standard  Glasso \cite{Rothman} and FF \cite{EstCovMatKron}. To achieve accurate covariance estimation for the sparse Kronecker product model, the minimal sample size needed is $n=\Omega((p+f)\log M)$.

The minimal sample size required to achieve accurate covariance estimation is graphically depicted in Fig. \ref{fig: visRate} for the special case $p=f$. The regions below the lines are the MSE convergence regions-i.e., the MSE convergence rate goes to zero as $p,n$ grow together to infinity at a certain growth rate controlled by these regions. It is shown that KGlasso allows the dimension $p$ to grow almost linearly in $n$ and still achieve accurate covariance estimation (see (\ref{KGL_perfect_rate})) and thus, uniformly outperforms FF, Glasso and the naive SCM estimators in the case both Kronecker factors are sparse.
\begin{figure}[htp]
	\centering
		\includegraphics[width=0.50\textwidth]{./Sim/visRates.pdf}
	\caption{Regions of convergence for KGlasso (below upper curve), FF (below second highest curve), Glasso (below third highest curve), and standard sample covariance matrix estimator (SCM) (bottom curve). These regions are obtained from the analytical expressions in equations (\ref{KGL_perfect_rate}),  (\ref{FF_rate_2}), (\ref{Glasso_rate})  and  (\ref{naive_SCM_rate}), respectively. The simulation shown in Fig. \ref{fig: KGlasso_FF_div_conv} establishes that the FF algorithm indeed diverges when the parameters $p$ and $n$ fall inbetween the KGlasso and FF curves in the above figure.}
	\label{fig: visRate}
\end{figure}

Although Thm. \ref{thm: KGL_optimal_rate} shows a rate on the inverse covariance matrix, this asymptotic rate can be shown to hold for the covariance matrix as well (see proof of Thm. \ref{thm: KGL_optimal_rate} in Appendix).


\subsection{Discussion}
Theorem \ref{thm: KGL_optimal_rate} is established using the large deviation bound in Lemma \ref{lemma: large_dev_Ted}. We provide some intuition on this bound below.
Assume that $\bX_{init} = \bX_0$, or $\bA_{init}=\bX_{init}^{-1}=\bA_0$. Define $\bW= \bX_0^{1/2} \otimes \bI_p$ and $\tilde{\bz}_t = \bW \bz_t$, with i.i.d. $\bz_t\sim N(\mathbf{0}, \bA_0\otimes \bB_0)$, $t=1,\dots,n$. Then, $\tilde{\bz}_t$ has block-diagonal covariance
\begin{equation*}
	\Cov(\tilde{\bz}_t) = \bI_p \otimes \bB_0.
\end{equation*}
When $\bW$ is applied to the transformed $pf\times pf$ sample covariance matrix, $\hat{\bS}_n^W := \bW \hat{\bS}_n \bW^T$, the first step of KGlasso produces an iterate $\hat{\bY}_n^{(1)}=\bG(\hat{\bB},\lambda_Y)$ with $\hat{\bB} = \frac{1}{p}\sum_{i=1}^p \hat{\bS}_n^W(i,i)$ (recall (\ref{B_update})). For suitable $\lambda_Y = \lambda_Y^{(1)}$, $\hat{\bY}_n^{(1)}$ converges to $\bY_0$ with respect to maximal elementwise norm at a rate $O_P\left( \sqrt{\frac{\log M}{np}} \right)$. The convergence of $\hat{\bY}_n^{(1)}$ is easily established by applying the Chernoff bound and invoking the jointly Gaussian property of the measurements and the block diagonal structure of $\Cov(\tilde{\bz}_t)$. Lemma \ref{lemma: large_dev_Ted} in the Appendix establishes that this rate holds even if $\bX_{init}\neq \bX_0$ in Assumption \ref{assumption_posdef_unif}. In view of the rate of convergence of $\hat{\bY}^{(1)}$, to achieve a reduction in the MSE of $\bY$, either the sample size $n$ or the dimension $p$ must increase. Lemma \ref{lemma: large_dev_Ted} provides a tight bound that makes the dependence of the convergence rate explicit in $p,f$ and $n$. Theorem \ref{thm: KGL_optimal_rate} uses Lemma \ref{lemma: large_dev_Ted} to show that KGlasso converges to $\bX_0\otimes \bY_0$ with rate $O_P\left( \sqrt{\frac{(p+f)\log M}{n}} \right)$ with respect to Frobenius norm.



\section{Simulation Results} \label{sec: simulations}

In this section, we empirically validate the convergence rates established in previous sections using Monte Carlo simulation. 

Each iteration of the KGlasso involves solving an $\ell_1$ penalized covariance estimation problem of dimension $100 \times 100$ (Step 6 and Step 8 of KGlasso specified by Algorithm \ref{alg: algKGL}). To solve these small sparse covariance estimation problems we used the Glasso algorithm of Hsieh {\it et al} \cite{HsiehNIPS11} where the Glasso stopping criterion was determined by monitoring when the duality gap falls below a threshold of $10^{-3}$.

In each of the simulations the true covariance matrix factors  $\bX_0=\bA_0^{-1}$ and $\bY_0=\bB_0^{-1}$ were unstructured randomly generated positive definite matrices. First, $p$ random nonzero elements were placed on the diagonal of a square $p \times p$ matrix $C$. Then, on average $p$ nonzero elements were placed on the off-diagonal and symmetry was imposed. On average, a total of $3p$ elements were nonzero. The resulting matrix $\tilde{\bC}$ was regularized to produce the sparse positive definite inverse covariance $\bY_0 = \tilde{\bC} + \rho \bI_f$, where $\rho = 0.5-\lambda_{min}(\tilde{\bC})$. A total of $N_{MC}=50$ simulation runs were performed for each sample size $n$, where $n$ ranged from 10 to 100. Performance assessment was based on normalized Frobenius norm error in the covariance and precision matrix estimates. The normalized error was calculated using
\begin{equation*}
	\sqrt{ \frac{1}{N_{MC}} \sum_{i=1}^{N_{MC}}{ \frac{\nn \bSigma_{0}-\hat{\bSigma}(i) \nn_F^2}{\nn \bSigma_0 \nn_F^2} } }
\end{equation*}
where $\hat{\bSigma}(i)$ is the covariance estimate for the $i$-th simulation. The same formula was  used to calculate the normalized error in the precision matrix $\hat{\bTheta}_0$.
In the implementation of KGlasso, the regularization parameters were chosen as follows. The initialization was $\bX_{init}=\bI_p$. The regularization parameters were selected as $\lambda_Y^{(1)}= c_y \sqrt{\frac{\log M}{np}}$, $\lambda_X^{(2)}= c_x \sqrt{\frac{\log M}{nf}} + \lambda_Y^{(1)}$, $\lambda_Y^{(2)} = \lambda_X^{(2)}$, $\lambda_X^{(3)}= \lambda_X^{(2)}$ and so on. We set $c_x=c_y=0.4$.

We considered the setting where $\bX_0$ and $\bY_0$ are large sparse matrices of dimension $p=f=100$ (see Fig. \ref{fig: sim4_matrices}) yeilding a covariance matrix $\bTheta_0$ of dimension $10,000 \times 10,000$. This dimension was too large for implementation of Glasso even when implemented using the state-of-the-art algorithm by Hsieh {\it et al} \cite{HsiehNIPS11}. Figures \ref{fig: sim4_frob_inv} and \ref{fig: sim4_frob_cov} compare the root-mean squared error (RMSE) performance in precision and covariance matrices as a function of $n$. As expected, KGlasso outperforms FF over the range of $n$ for both covariance and inverse covariance estimation problems. KGlasso outperforms FF in the small-sample regime since it exploits sparsity in addition to Kronecker structure.

For $n=10$, there is a $72 \%$ ($\approx 5.53$ dB) RMSE reduction for the precision matrix and $49 \%$ RMSE reduction for the covariance matrix when using KGlasso instead of FF. For $n=100$, there is a $51 \%$ ($\approx 3.10$ dB) RMSE reduction for the precision matrix and $41 \%$ RMSE reduction for the covariance matrix. For the small sample regime, there is approximately a $5.53$ dB reduction for the precision matrix, which is a significant performance gain.
\begin{figure}[htp]
	\centering
%		\includegraphics[width=0.40\textwidth]{./Simulations/KGL_FF_sim_21-Mar-2012-22-15_matrices.png}
		\includegraphics[width=0.40\textwidth]{./Sim/pf100matrices.png}
	\caption{Sparse Kronecker matrix representation. Left panel: left Kronecker factor. Right panel: right Kronecker factor. As the Kronecker-product covariance matrix is of dimension $10,000 \times 10,000$ standard Glasso is not practically implementable for this example.}
	\label{fig: sim4_matrices}
	\vfill
		 \includegraphics[width=0.40\textwidth]{./Sim/KGL_FF_sim_21-Mar-2012-22-15_Frob_inv_std_all.pdf}
	\caption{Normalized RMSE performance for precision matrix as a function of sample size $n$. KGlasso (Kronecker graphical lasso) uniformly outperforms FF (flip-flop) algorithm for all $n$. Here, $p=f=100$ and $N_{MC}=50$. The error bars are centered around the mean with $\pm$ one standard deviation. For $n=10$, there is a $72 \%$ RMSE reduction.  }
	\label{fig: sim4_frob_inv}
	\vfill
		 \includegraphics[width=0.40\textwidth]{./Sim/KGL_FF_sim_21-Mar-2012-22-15_Frob_cov_std_all.pdf}
	\caption{Normalized RMSE performance for covariance matrix as a function of sample size $n$. KGlasso (Kronecker graphical lasso) uniformly outperforms FF (flip-flop) algorithm for all $n$. Here, $p=f=100$ and $N_{MC}=50$. The error bars are centered around the mean with $\pm$ one standard deviation. For $n=10$, there is a $49 \%$ RMSE reduction. }
	\label{fig: sim4_frob_cov}
\end{figure}


\subsection{Empirical Rate Comparison}
Next, we illustrate the rates obtained in for the dimension setting $p(n)=f(n)= \lceil 8 n^\alpha \rceil$, where $\alpha\in\{0.1,0.2,0.3\}$. According to the theory developed, for large $n$, the MSE converges to zero at a certain convergence rate. The predicted rates of FF and KGlasso are fitted on top of the empirical MSE curves by ensuring intersection at $n=1000$. Fig. \ref{fig: KGlasso_FF_curves} shows that the empirical rates match the predicted rates well.
\begin{figure}[htp]
	\centering
		\includegraphics[width=0.50\textwidth]{./Sim/KGlasso_FF_sim_19-Mar-2012-20-19_summary.pdf}
	\caption{ Precision Matrix MSE convergence as a function of sample size $n$ for FF and KGlasso. The dimensions of the Kronecker factor matrices grow as a function of $n$ as: $p(n)=f(n)=\lceil 8 \cdot n^\alpha \rceil$. The true Kronecker factors were set to identity (so their inverses are fully sparse). The predicted MSE curves according to Thm. \ref{thm: FF_optimal_rate} and Thm. \ref{thm: KGL_optimal_rate} are also shown. For both KGlasso and FF, the predicted MSE matches the empirical MSE well, thus verifying the rate expressions (\ref{FF_rate_2}) and (\ref{KGL_perfect_rate}). }
	\label{fig: KGlasso_FF_curves}
\end{figure}

We also show a borderline case $p=f=\lceil n^{0.6} \rceil$. In this case, according to Thm. \ref{thm: FF_optimal_rate} and Thm. \ref{thm: KGL_optimal_rate}, the FF diverges (MSE increases in $n$), while the KGlasso converges (MSE decreases in $n$). This is illustrated in Fig. \ref{fig: KGlasso_FF_div_conv}. Our predicted rates are plotted on top of the empirical curves.
\begin{figure}[htp]
	\centering
		\includegraphics[width=0.50\textwidth]{./Sim/KGlasso_FF_sim_20-Mar-2012-22-54_inv_KGL_FF_pred.pdf}
	\caption{ Precision Matrix MSE as a function of sample size $n$ for FF and KGlasso. The dimensions of the Kronecker factor matrices grow as a function of $n$ as: $p(n)=f(n)=\lceil n^{0.6} \rceil$. The true Kronecker factors were set to identity (so their inverses are fully sparse). The predicted MSE curves according to Thm. \ref{thm: FF_optimal_rate} and Thm. \ref{thm: KGL_optimal_rate} are also shown. As predicted by our theory, and by the predicted convergent regions of $(n,p)$ for FF and KGlasso in Fig. \ref{fig: visRate},  the MSE of the FF diverges while the MSE of the KGlasso converges as $n$ increases. }
	\label{fig: KGlasso_FF_div_conv}
\end{figure}




\section{Conclusion}
We established high dimensional consistency for Kronecker Glasso algorithms that use iterative $\ell_1$-penalized likelihood optimization that exploit both Kronecker structure and sparsity of the covariance. A tight MSE convergence rate was derived for KGlasso, showing significantly better MSE performance than standard Glasso \cite{Rothman, ModelSel} and FF \cite{EstCovMatKron}. Simulations validated our theoretical predictions. 


\section*{Acknowledgement}
The authors thank Prof. Mark Rudelson for very helpful discussions on large deviation theory.


% Appendices with technical proofs
\appendices

\section{Proof of Lemma \ref{dual_lemma}}
\input{appendix_dual_lemma}

%\section{Limit Point Characterization of KGlasso}
%%MOVE THIS ENTIRE SECTION (TO POINT A) TO THE APPENDIX
%%We first show that the KGlasso iterations converge to a fixed point.
%Let $J_\lambda(\bX,\bY)$ be as defined in (\ref{J_lambda_func}) and define $J_\lambda^{(k)}=J_\lambda(\mathbf{X}^{(k)},\mathbf{Y}^{(k)})$ for $k=0,1,2,\dots$.
%
%\begin{theorem} \label{convergence_fixed_point}
%	If $\hat{\bS}_n$ is positive definite, the KGlasso iterations converge to a fixed point. Also, we have $J_\lambda^{(k)} \searrow J_\lambda^{(\infty)}$.
%\end{theorem}
%\begin{IEEEproof}
%	Recall that the basic optimization problem (\ref{opt_prob}) is
%	\begin{equation*}
%		\min_{\bX \in S_{++}^p, \bY \in S_{++}^f} J_\lambda(\bX,\bY)
%	\end{equation*}
%	Let $J^*:=\inf_{\mathbf{X} \in S_{++}^p, \mathbf{Y} \in S_{++}^f}{J_\lambda(\mathbf{X},\mathbf{Y})}$ be the optimal primal value. Note that $J_\lambda^* > -\infty$ when $\hat{\bS}_n \in S_{++}^{pf}$. Consider Algorithm \ref{alg: algKGL}. Assuming $\bX^{(0)}$ and $\hat{\bS}_n$ are p.d., Lemma \ref{dual_lemma} implies that all the iterates are p.d. and
%	\begin{equation} \label{monotonicity}
%		J_\lambda(\bX^{(k)},\bY^{(k)}) \leq	J_\lambda(\bX^{(k-1)},\bY^{(k)}) \leq J_\lambda(\bX^{(k-1)},\bY^{(k-1)})
%	\end{equation}
%	By induction on the number of iterations, we conclude that the iterates yield a nonincreasing sequence of objective functions. Since $\bar{\lambda}_X |\bX|_1, \bar{\lambda}_Y |\bY|_1 \geq 0$, the maximum likelihood objective function provides a lower bound to the optimal primal value
%	\begin{equation} \label{LB_ML}
%		J_\lambda(\bX^{(k)},\bY^{(k)}) \geq J_\lambda^* \geq J_\lambda(\bX_{ML},\bY_{ML}) > -\infty
%	\end{equation}
%	Thus, the sequence $\{J_\lambda^{(k)}:k\geq 0\}$ forms a nonincreasing sequence bounded below \footnote{This follows since for $n>pf$, the log-likelihood function is bounded above by the log-likelihood evaluated at the sample mean and sample covariance matrix.}. The monotone convergence theorem for sequences \cite{IntroRealAnalysis} implies that $\{J_\lambda^{(k)}\}$ converges monotonically to $J_\lambda^{(\infty)}=\inf_k{J_\lambda^{(k)}}$. By the alternating minimization, we conclude that the sequence of iterates $\{(\mathbf{X}^{(k)},\mathbf{Y}^{(k)})\}_k$ converges since the minimizer at each Glasso step is unique.
%\end{IEEEproof}
%
%
%%The following analysis uses Theorem \ref{convergence_fixed_point} to prove convergence of the KGlasso algorithm to a local minimum. To do this, we consider a more general setting. The KGlasso algorithm is a special case of the block-coordinate descent algorithm. Assuming a $k$-fold Kronecker product structure for the covariance matrix, the optimization problem (\ref{J_lambda_func}) can be written as:
%%\begin{equation} \label{obj_func_matrix}
%%	J_\lambda(\bX_1,\dots,\bX_k) = J_0(\bX_1,\dots,\bX_k) + \sum_{i=1}^k{J_i(\bX_i) + \bar{\lambda}_i \eta_1(\bX_i)}
%%\end{equation}
%%where $\bar{\lambda}_i=\lambda_i/d_i$, $\bX_i \in S_{++}^{d_i}$, $\eta_1(\bX_i):=|\bX_i|_1$, $J_0(\bX_1,\dots,\bX_k):=\tr((\bX_1\otimes \bX_2 \otimes \dots \otimes \bX_k)\hat{\bS}_n)$ and $J_i(\bX_i) = - \prod_{i' \neq i}{d_{i'}} \cdot \log\det(\bX_i)$ for $i=1,\dots,k$.
%%
%%Without loss of generality, by reshaping matrices into appropriate vectors, (\ref{obj_func_matrix}) can be rewritten as:
%%\begin{equation} \label{general_f}
%%	J_\lambda(\mathbf{x}_1,\dots,\mathbf{x}_k) = J_0(\mathbf{x}_1,\dots,\mathbf{x}_k) + \sum_{i=1}^k{J_i(\mathbf{x}_i)+\bar{\lambda}_i \eta_i(\mathbf{x}_i)}
%%\end{equation}
%%where the optimization variable is $\mathbf{x}:=[\mathbf{x}_1^T, \mathbf{x}_2^T,\dots,\mathbf{x}_k^T]^T \in \RR^{d'}$, where $\mathbf{x}_i\in \RR^{d_i^2}$ and $d'=\sum_{i=1}^k{d_i^2}$. For example, $\eta_i(\mathbf{X}_i)=|\mathbf{X}_i|_1=\nn \vec(\mathbf{X}_i) \nn_1 = \nn \mathbf{x}_i\nn_1 = \eta_i(\mathbf{x}_i)$. The mapping $\{J_i\}_{i=0}^k$ can be similarly written in terms of the vectors $\mathbf{x}_i$ instead of the matrices $\mathbf{X}_i$.
%%
%%
%%The general optimization problem of interest here is:
%%\begin{equation} \label{gen_func}
%%	\min_{\mathbf{x} \in \RR^{d'}} {J_\lambda(\mathbf{x})} \text{ subject to } \vec^{-1}(\mathbf{x}_i)=\bX_i \in S_{++}^{d_i}, i=1,\dots,k
%%\end{equation}
%%The positive definiteness constraints are automatically taken care of by the construction of the algorithm (see Lemma \ref{dual_lemma}.4). Let the dimension of the covariance matrix be denoted by $d:=\prod_{i=1}^k{d_i}$. We assume $n>d$. To solve (\ref{gen_func}), a block coordinate-descent penalized algorithm can be constructed \cite{TsiligkaridisTSP} that iteratively optimizes (\ref{obj_func_matrix}).
%%
%%%\begin{algorithm}
%%%\caption{Block Coordinate-Descent Penalized Algorithm}
%%%\label{alg: alg2}
%%%\begin{algorithmic}[1]
%%%\STATE \textbf{Input:}  {$\hat{\bS}_n$, $d_i$, $n$, $\epsilon>0$, $\lambda_i >0$}
%%%\STATE \textbf{Output:} {$\hat{\bTheta}$}
%%%
%%%	\STATE Initialize $\mathbf{X}_1^0, \mathbf{X}_2^0, \dots, \mathbf{X}_k^0$ matrices as positive definite matrices, e.g., scaled identity.
%%%	
%%%	\STATE $\hat{\mathbf{\Theta}}_0 \leftarrow \mathbf{X}_1^0\otimes \mathbf{X}_2^0 \otimes \dots \otimes \mathbf{X}_k^0$\;
%%%	
%%%	\STATE $m \leftarrow 0$\;
%%%	\REPEAT {
%%%		\STATE $\hat{\mathbf{\Theta}}_{\text{prev}} \leftarrow \hat{\mathbf{\Theta}}$
%%%		
%%%		\STATE $\mathbf{X}_1^m \leftarrow \arg \min_{\mathbf{A}_1 \succ 0} J_\lambda(\mathbf{A}_1,\mathbf{X}_2^{m-1},\dots,\mathbf{X}_k^{m-1})$
%%%		\STATE $\mathbf{X}_2^m \leftarrow \arg \min_{\mathbf{A}_2 \succ 0} J_\lambda(\mathbf{X}_1^m,\mathbf{A}_2,\dots,\mathbf{X}_k^{m-1})$
%%%		\STATE $\vdots$
%%%		\STATE $\mathbf{X}_k^m \leftarrow \arg \min_{\mathbf{A}_k \succ 0} J_\lambda(\mathbf{X}_1^m,\mathbf{X}_2^m,\dots,\mathbf{A}_k)$
%%%		
%%%		\STATE $\hat{\mathbf{\Theta}} \leftarrow \mathbf{X}_1^m \otimes \mathbf{X}_2^m \otimes \dots \otimes \mathbf{X}_k^m$
%%%		
%%%		\STATE $m \leftarrow m + 1$
%%%	}	\UNTIL {$\frac{\nn \hat{\mathbf{\Theta}}_{\text{prev}}-\hat{\mathbf{\Theta}} \nn}{\nn \hat{\bTheta}_{\text{prev}} \nn} \leq \epsilon$}
%%%	
%%%\end{algorithmic}
%%%\end{algorithm}
%%%
%%%\begin{remark}
%%%	The positive definiteness constraint at each coordinate descent iteration of Algorithms \ref{alg: algKGL} and \ref{alg: alg2} need not be explicit since the objective function $J_\lambda(\cdot)$ acts as a logarithmic barrier function.
%%%\end{remark}
%%
%%Note that Algorithm \ref{alg: algKGL} would be a special case of this algorithm. An extension of Theorem \ref{convergence_fixed_point}, assuming $n>d$ or $J_\lambda^*> -\infty$, based on induction, can be used to show that the limit points of the sequence of iterates $(\mathbf{x}^m)_{m\geq 0}=(\mathbf{x}_1^m,\dots,\mathbf{x}_k^m)_{m \geq 0}$ are fixed points.
%
%\begin{remark}
%	Note that a necessary condition for $\mathbf{x}^*$ to minimize $J_\lambda$ is $0\in \partial J_\lambda(\mathbf{x}^*)$. This is not sufficient however.
%\end{remark}
%
%We next show that the limit point(s) of $(\mathbf{x}^m)_{m\geq 0}$ are nonempty and are local minima.
%
%\begin{theorem} \label{convergence_theorem_critical_manyfactors}
%	Let $(\mathbf{x}^m)=(\mathbf{x}_1^m,\dots,\mathbf{x}_k^m)_{m\geq 0}$ be a sequence generated by the block-coordinate descent algorithm. Assume $n>d$.
%
%\begin{enumerate}
%	%\item The KGlasso iterations converge to a local minimum or a local maximum of the objective function.
%	\item The iterates converge to a local minimum or a local maximum of the objective function.
%	\item If the initial condition $\mathbf{x}^0$ is not a local maximum then the algorithm converges to a local minimum.
%\end{enumerate}
%\end{theorem}
%\begin{IEEEproof}
%	Due to limited space, the proof is included in \cite{TsiligkaridisTSP}.
%\end{IEEEproof}
%%\begin{IEEEproof}
%%	Define the set of critical points $C_J:=\{\mathbf{x}: 0\in \partial J(\mathbf{x})\}=C_{J,min} \cup C_{J,saddle} \cup C_{J,max}$, where $C_{J,min}$ contains all the local minima, $C_{J,saddle}$ contains all the saddle points and $C_{J,max}$ contains all the local maxima.
%%	
%%\begin{enumerate}
%%	\item Let $L(\mathbf{x}^0)=L(\mathbf{x}_1^0,\dots,\mathbf{x}_k^0)$ be the set of all limit points of $(\mathbf{x}^m)_{m\geq 0}$ starting from $\mathbf{x}^0$. By joint continuity of the objective, we have $J_\lambda(\mathbf{x}^{m_j}) \to J_\lambda(\mathbf{x}^*)$ as $j\to\infty$.
%%	%By the continuity of $J_i$, we have $\sum_{i=1}^k{J_i(\mathbf{x}_i^{m_j})} \to \sum_{i=1}^k{J_i(\mathbf{x}_i^*)}$ as $j \to \infty$. Since $J_0(\cdot)$ is jointly continuous, $J_0(\mathbf{x}_1^{m_j},\dots,\mathbf{x}_k^{m_j}) \to J_0(\mathbf{x}_1^*,\dots,\mathbf{x}_k^*)$. By continuity of $\eta_i(\cdot)$, $\sum_{i=1}^k{\bar{\lambda}_i \eta_i(\mathbf{x}_i^{m_j})} \to \sum_{i=1}^k \bar{\lambda}_i \eta_i(\mathbf{x}_i^*)$. Thus, $J_\lambda(\mathbf{x}^{m_j}) \to J_\lambda(\mathbf{x}^*)$ as $j\to\infty$.
%%	Note that $\nabla J_0$ is uniformly continuous on bounded subsets of its domain. Define the sequence $\{(x_i^m)^o\}_m$ as
%%	\begin{align*}
%%		(\mathbf{x}_1^m)^\circ & := \nabla_{\mathbf{x}_1}J_0(\mathbf{x}_1^m,\mathbf{x}_2^m\dots,\mathbf{x}_k^m) \\
%%			&\quad - \nabla_{\mathbf{x}_1}J_0(\mathbf{x}_1^m,\mathbf{x}_2^{m-1}\dots,\mathbf{x}_k^{m-1}) \\
%%		(\mathbf{x}_2^m)^\circ & := \nabla_{\mathbf{x}_2}J_0(\mathbf{x}_1^m,\mathbf{x}_2^m\dots,\mathbf{x}_k^m) \\
%%			&\quad - \nabla_{\mathbf{x}_2}J_0(\mathbf{x}_1^m,\mathbf{x}_2^m,\mathbf{x}_{3}^{m-1}\dots,\mathbf{x}_k^{m-1}) \\
%%		&\vdots \\
%%		(\mathbf{x}_j^m)^\circ & := \nabla_{\mathbf{x}_j}J_0(\mathbf{x}_1^m,\mathbf{x}_2^m\dots,\mathbf{x}_k^m) \\
%%			&\quad - \nabla_{\mathbf{x}_j}J_0(\mathbf{x}_1^m,\dots,\mathbf{x}_j^m,\mathbf{x}_{j+1}^{m-1}\dots,\mathbf{x}_k^{m-1}) \\
%%		&\vdots \\
%%		(\mathbf{x}_k^m)^\circ &:= 0
%%	\end{align*}
%%	Then, it can be shown $((\mathbf{x}_1^m)^\circ,\dots,(\mathbf{x}_k^m)^\circ) \in \partial J_\lambda(\mathbf{x}_1^m,\dots,\mathbf{x}_k^m)$. This implies $((\mathbf{x}^{m_j})^{\circ}) \in \partial J_\lambda(\mathbf{x}^{m^j})$ \cite{TsiligkaridisTSP}. Since the subsequence $(\mathbf{x}^{m_j})_j$ is convergent, we have $(\mathbf{x}^{m_j})^\circ \to 0$ as $j \to \infty$ \cite{TsiligkaridisTSP}. Since $\partial J_\lambda(\mathbf{x}^{m_j})$ is closed (see Theorem 8.6 in \cite{VariationalAnalysis}) for all $j$, we conclude that $\mathbf{x}^* \in C_J$. Thus, $L(\mathbf{x}^0) \subseteq C_J$.
%%	
%%	%Now, it can be shown \cite{TsiligkaridisTSP} that $((\mathbf{x}^{m_j})^{\circ}) \in \partial J_\lambda(\mathbf{x}^{m^j})$. Since the subsequence $(\mathbf{x}^{m_j})_j$ is convergent, we have \cite{TsiligkaridisTSP} $(\mathbf{x}^{m_j})^\circ \to 0$ as $j \to \infty$. As a result, since $\partial J_\lambda(\mathbf{x}^{m_j})$ is closed (see Theorem 8.6 in \cite{VariationalAnalysis}) for all $j$, we conclude that $\mathbf{x}^* \in C_J$. Thus, $L(\mathbf{x}^0) \subseteq C_J$.
%%	
%%	%We have thus proved that limit points are critical points of the objective function.
%%	We can easily rule out convergence to saddle points. Theorem \ref{convergence_fixed_point} implies that $L(\mathbf{x}^0)$ is nonempty and singleton.
%%	
%%	\item Let $\mu(\cdot)$ denote the point-to-point mapping during one iteration step, i.e., $\mathbf{x}^{m+1}=\mu(\mathbf{x}^m)$. We show that if $\mathbf{x}^0 \notin C_J$, then $L(\mathbf{x}^0) \subseteq C_{J,min} \cup C_{J,saddle}$. The result then follows by using the proof of the first part. To this end, let $\mathbf{x}^{'}$ be a fixed point under $\mu$, i.e., $\mu(\mathbf{x}^{'})=\mathbf{x}^{'}$. Then, the subiteration steps of the algorithm yield $0\in \partial_{\mathbf{x}_i} J_\lambda(\mathbf{x}_1^{'},\dots,\mathbf{x}_k^{'})$ for $i=1,\dots,k$, which implies $0\in \partial J_\lambda(\mathbf{x}^{'})$, i.e., $\mathbf{x}^{'}\in C_J$. The contrapositive implies that if $\mathbf{x}\notin C_J$, then $J_\lambda(\mu(\mathbf{x}))<J_\lambda(\mathbf{x})$. A simple induction on the number of iterations then concludes the proof.
%%\end{enumerate}
%%	
%%\end{IEEEproof}%
%
%As a consequence of Theorem \ref{convergence_theorem_critical_manyfactors}, we have the following corollary.
%
%\begin{corollary} \label{cor_localmin_singleton}
%	Assuming $(\mathbf{X}^{(0)},\mathbf{Y}^{(0)})$ is not a local maximum and $n>pf$, the KGlasso algorithm converges to a local minimizer of the objective function (\ref{J_lambda_func}).
%\end{corollary}
%
%%POINT A THIS IS THE END OF WHAT SHOULD GO INTO APPENDIX



%\section{Proof of Theorem \ref{convergence_fixed_point}}
%\input{appendix_convergence_fixed_point}

%\section{Proof of Theorem \ref{convergence_theorem_critical}}
%Define the set of critical points $C_J:=\{\mathbf{x}: 0\in \partial J(\mathbf{x})\}=C_{J,min} \cup C_{J,saddle} \cup C_{J,max}$, where $C_{J,min}$ contains all the local minima, $C_{J,saddle}$ contains all the saddle points and $C_{J,max}$ contains all the local maxima.
%\input{appendix_KGL_convergence_proof}

% statistical analysis

\section{Lemma \ref{lemma: large_dev_Ted}}
The following lemma will be used in the proof of Thm. \ref{thm: FF_optimal_rate} and Thm. \ref{thm: KGL_optimal_rate}.
\input{appendix_large_dev_Ted}

\section{Proposition \ref{prop: Glasso_optimal_rate}}
\input{appendix_Glasso_optimal_rate}

\section{Proof of Theorem \ref{thm: FF_optimal_rate}}
\input{appendix_FF_optimal_rate_achievability_new_detailed_short_arb}

\section{Proof of Theorem \ref{thm: KGL_optimal_rate}}
\input{appendix_KGL_optimal_rate_new_detailed_short_arb}



\bibliographystyle{IEEEtran}
\bibliography{myRefs}


\end{document}