text:
\section{Appendix}

\iffalse
\subsection{Coordinate Descent Updates for Newton Coordinate Descent Method}
In this section, we provide details of the coordinate-descent update equations
for the Newton coordinate descent algorithm.
In the Newton coordinate descent method, the coordinate-descent method for
computing a generalized Newton direction iteratively updates
each element of $\Delta_\mathbf{\Lambda}$ as follows:
\begin{align*}
(\Delta_\mathbf{\Lambda})_{ij} \leftarrow& (\Delta_\mathbf{\Lambda})_{ij} - c_\mathbf{\Lambda} + S_{\lambda_\mathbf{\Lambda}/a_\mathbf{\Lambda}}(c_\mathbf{\Lambda}-\frac{b_\mathbf{\Lambda}}{a_\mathbf{\Lambda}}),
\end{align*}
where
$S_{r}(w)=\textrm{sign}(w)\max(|w|-r,0)$ is the soft-thresholding operator and
\begin{align*}
a_\mathbf{\Lambda} =& \mathbf{\Sigma}_{ij}^2+\mathbf{\Sigma}_{ii}\mathbf{\Sigma}_{jj}+\mathbf{\Sigma}_{ii}\mathbf{\Psi}_{jj}+2\mathbf{\Sigma}_{ij}\mathbf{\Psi}_{ii} \\
b_\mathbf{\Lambda} =& (\mathbf{S_{yy}})_{ij}-\mathbf{\Sigma}_{ij}-\Psi_{ij} + (\mathbf{\Sigma}\Delta_\mathbf{\Lambda}\mathbf{\Sigma})_{ij}
-\bPhi_{ij} - \bPhi_{ji} + (\mathbf{\Psi}\Delta_\mathbf{\Lambda} \mathbf{\Sigma})_{ij} +
(\mathbf{\Psi}\Delta_\mathbf{\Lambda} \mathbf{\Sigma})_{ji} \\
c_\mathbf{\Lambda} =& \mathbf{\Lambda}_{ij} + (\Delta_\mathbf{\Lambda})_{ij}
\end{align*}
given $\mathbf{\Sigma} = \mathbf{\Lambda}^{-1}$, $\mathbf{\Psi} := \mathbf{\Sigma}\mathbf{\Theta}^T \mathbf{S_{xx}}\mathbf{\Theta}\mathbf{\Sigma}$,
and $\bPhi := \mathbf{\Sigma}\mathbf{\Theta}^T \mathbf{S_{xx}} \Delta_\mathbf{\Theta} \mathbf{\Sigma}$.

Similarly, in each coordinate-descent iteration,
each element of $\Delta_\mathbf{\Theta}$ is updated as follows:
\begin{align*}
(\Delta_\mathbf{\Theta})_{ij} \leftarrow& (\Delta_\mathbf{\Theta})_{ij} - c_\mathbf{\Theta} + S_{\lambda_\mathbf{\Theta}/a_\mathbf{\Theta}}(c_\mathbf{\Theta}-\frac{b_\mathbf{\Theta}}{a_\mathbf{\Theta}}),
\end{align*}
where
\begin{align*}
a_\mathbf{\Theta} =& \mathbf{\Sigma}_{jj}(\mathbf{S_{xx}})_{ii} \\
b_\mathbf{\Theta} =& 2(\mathbf{S_{xy}})_{ij} + 2(\mathbf{S_{xx}}\mathbf{\Theta}\mathbf{\Sigma})_{ij}
+ 2(\mathbf{S_{xx}}\Delta_\mathbf{\Theta} \mathbf{\Sigma})_{ij} - 2(\mathbf{S_{xx}}\mathbf{\Theta}\mathbf{\Sigma}\Delta_\mathbf{\Lambda} \mathbf{\Sigma})_{ij} \\
c_\mathbf{\Theta} =& \mathbf{\Theta}_{ij} + (\Delta_\mathbf{\Theta})_{ij}.
\end{align*}
\fi

\subsection{Coordinate Descent Updates for Alternating Newton Coordinate Descent Method}

In our alternating Newton coordinate descent algorithm,
each element of $\Delta_\mathbf{\Lambda}$ is updated as follows:
\begin{align*}
(\Delta_\mathbf{\Lambda})_{ij} \leftarrow& (\Delta_\mathbf{\Lambda})_{ij} - c_\mathbf{\Lambda} + S_{\lambda_\mathbf{\Lambda}/a_\mathbf{\Lambda}}(c_\mathbf{\Lambda}-\frac{b_\mathbf{\Lambda}}{a_\mathbf{\Lambda}}),
\end{align*}
where $S_{r}(w)=\textrm{sign}(w)\max(|w|-r,0)$ is the soft-thresholding operator and
\begin{align*}
a_\mathbf{\Lambda} =& \mathbf{\Sigma}_{ij}^2+\mathbf{\Sigma}_{ii}\mathbf{\Sigma}_{jj}+\mathbf{\Sigma}_{ii}\mathbf{\Psi}_{jj}+2\mathbf{\Sigma}_{ij}\mathbf{\Psi}_{ii} \\
b_\mathbf{\Lambda} =& (\mathbf{S_{yy}})_{ij}-\mathbf{\Sigma}_{ij}-\mathbf{\Psi}_{ij} + (\mathbf{\Sigma}\Delta_\mathbf{\Lambda}\mathbf{\Sigma})_{ij}
 + (\mathbf{\Psi}\Delta_\mathbf{\Lambda} \mathbf{\Sigma})_{ij} +(\mathbf{\Psi}\Delta_\mathbf{\Lambda} \mathbf{\Sigma})_{ji} \\
c_\mathbf{\Lambda} =& \mathbf{\Lambda}_{ij} + (\Delta_\mathbf{\Lambda})_{ij}.
\end{align*}
For $\mathbf{\Theta}$, we perform coordinate-descent updates directly on $\mathbf{\Theta}$ without
forming a second-order approximation of the log-likelihood to find a Newton direction,
as follows:
\begin{align*}
\mathbf{\Theta}_{ij} \leftarrow& \mathbf{\Theta}_{ij} - c_\mathbf{\Theta} + S_{\lambda_\mathbf{\Theta}/a_\mathbf{\Theta}}(c_\mathbf{\Theta}-\frac{b_\mathbf{\Theta}}{a_\mathbf{\Theta}}),
\end{align*}
where
\begin{align*}
a_\mathbf{\Theta} =& \mathbf{\Sigma}_{jj}(\mathbf{S_{xx}})_{ii} \\
b_\mathbf{\Theta} =& 2(\mathbf{S_{xy}})_{ij} + 2(\mathbf{S_{xx}}\mathbf{\Theta}\mathbf{\Sigma})_{ij} \\
c_\mathbf{\Theta} =& \mathbf{\Theta}_{ij}.
\end{align*}

\iffalse
\subsection{Time Complexity Analysis of Alternating Newton Block Coordinate Descent}

In this section we describe the time complexity
of the alternating Newton block coordinate descent method.
The active set sizes for $\mathbf{\Lambda}$ and $\mathbf{\Theta}$
are $m_\mathbf{\Lambda}$ and $m_\mathbf{\Theta}$, respectively.
Also, $K \sim 10$ is the average number of conjugate gradient iterations.

\subsubsection{Cost of updating $\mathbf{\Lambda}$}

The cost of computing $\mathbf{\Sigma}_{C_z}$ and $\mathbf{\Psi}_{C_z}$ for all $z \in \{1,\dots,k_\mathbf{\Lambda}\}$
is $O(qm_\mathbf{\Lambda} K)$ and $O(nq^2)$, respectively, regardless of sparsity in $\mathbf{\Lambda}$.
However, the cost of computing $\mathbf{\Sigma}_{C_r}$ and $\mathbf{\Psi}_{C_r}$
for all $q \ne z$ for all $z \in \{1,\dots,k_\mathbf{\Lambda}\}$
depends on the block-diagonal structure in $\mathbf{\Lambda}$.
Let $B_{zr} = |\{j | j \in C_r \text{ and } \exists i \in C_z \text{ s.t. } (i,j) \in \mathcal{S}_\mathbf{\Lambda}\}|$
be the number of columns in $C_r$ with at least one active variable in $(C_z, C_r)$.
Then $B = \sum_{z\ne r}{|B_{zr}|}$ is the number of times
we need to recompute a column in $\mathbf{\Sigma}$ and $\mathbf{\Psi}$ while traversing off-diagonal blocks.
This implies that coordinate descent for $\mathbf{\Lambda}$ costs $O((q+B) m_\mathbf{\Lambda} K)$ for computing $\mathbf{\Sigma}$,
and $O((q+B)nq)$ for computing $\mathbf{\Psi}$.
As for alternating Newton coordinate descent,
the cost of each individual coordinate update is $O(q)$.
The overall cost is then $O((q+B)(m_\mathbf{\Lambda} K + nq) + m_\mathbf{\Lambda} q)$.

The dependence on $B$ highlights how recognizing diagonal structure improves performance.
In the best case where $\mathbf{\Lambda}$ becomes perfectly block-diagonal,
$B = 0$ and the overall cost of updating $\mathbf{\Lambda}$ reduces to
$O(q m_\mathbf{\Lambda} K + n q^2 + m_\Lambda q)$.
In the worst case where we need to recompute all $q/k_\mathbf{\Lambda}$ columns of
$\mathbf{\Sigma}_{C_r}$ and $\mathbf{\Psi}_{C_r}$ in each of the off-diagonal blocks,
$B \sim q k_\mathbf{\Lambda}$ and the overall cost of updating all blocks is
$O(k_\mathbf{\Lambda} q m_\mathbf{\Lambda} K + k_\mathbf{\Lambda} nq^2 + m_\mathbf{\Lambda} q)$.

\subsubsection{Time Cost of Updating $\mathbf{\Theta}$}

To sweep through all the blocks in $\mathbf{\Theta}$, we compute $\mathbf{\Sigma}$ only once
while potentially computing columns of $\mathbf{S_{xx}}$ multiple times.
Thus, the cost of computing $\mathbf{\Sigma}$ is $O(qm_\mathbf{\Lambda} K)$
and does not depend on sparsity in $\mathbf{\Theta}$.
The cost of computing $\mathbf{V} = \mathbf{\Theta} \mathbf{\Sigma}$ is $O(q m_\mathbf{\Theta})$
and also does not depend on sparsity or the number of blocks.
For the cost of computing $\mathbf{S_{xx}}$, we assume that
$\mathbf{\Theta}$ has $\tilde{p} \le p$ non-empty rows, so that
the number of blocks with at least one active variable is at most $\tilde{p} k_\mathbf{\Theta}$.
Let $B = I_\phi(\mathcal{S}_{(i,C_r)})$ be the number of blocks
$(i,C_r)$ where we need to compute $\mathbf{S_{xx}}_i$.
The cost of computing $\mathbf{S_{xx}}_i$ for each cache miss is $O(n\tilde{p})$.
Thus, computing $\mathbf{S_{xx}}$ for updating all blocks in $\mathbf{\Theta}$ is
$O(B n \tilde{p})$ which is upper-bounded by $O(k_\mathbf{\Theta} n \tilde{p}^2)$.
The cost of each individual coordinate update is $O(\tilde{p}+ q/k_\mathbf{\Theta})$.
The overall cost is then $O(B n\tilde{p} + qm_\mathbf{\Lambda} K + q m_\mathbf{\Theta}
+ m_\mathbf{\Theta}(\tilde{p} + q/k_\mathbf{\Theta}))$.

For large problems the time needed to compute $\mathbf{S_{xx}}$ dominates.
In the worst case where $\tilde{p} = p$ and
graph clustering is ineffectual so that $B = k_\mathbf{\Theta} p$,
updating all blocks requires computing $\mathbf{S_{xx}}$ a total of
$k_\mathbf{\Theta}$ times with cost $O(k_\mathbf{\Theta} n p^2)$.
However, because we store only one row of $\mathbf{S_{xx}}$ at a time,
$k_\mathbf{\Theta}$ will be small, even for problems with large $p$.
In fact, for large sparse problems,
the effect of sparsity can outweigh the additional $k_\mathbf{\Theta}$ factor
so that $Bn\tilde{p}$ can be smaller than $np^2$.
\fi

\begin{figure*}[h!]
\centering
\begin{tabular}{@{}c@{\hspace{10pt}}c}
\includegraphics[scale=0.45]{figures/scaling-N.pdf}
& \includegraphics[scale=0.45]{figures/accuracy-N.pdf} \\
\vspace{-3pt}
(a) & (b)
\end{tabular}
\vspace{-0.05in}
\caption{Results from varying sample size $n$ on chain graph with $p=q=10{,}000$.
(a) Comparison of computation time of different methods.
(b) Comparison of edge recovery accuracy as measured by $F_1$-score.
}
\label{fig:scaling-n}
\end{figure*}

\subsection{Additional Results from Synthetic Data Experiments}

We compare the performance of the different algorithms on synthetic datasets
with different sample sizes $n$, using a chain graph structure
with $p=q=10{,}000$.
Figure \ref{fig:scaling-n}(a) shows that our methods run significantly faster
than the previous method across all sample sizes.
In Figure \ref{fig:scaling-n}(b) we measure the accuracy in recovering
the true chain graph structure in terms of $F_1$-score for different sample sizes $n$.
At convergence, $F_1$-score was the same for all methods to three significant digits.
As expected, the accuracy improves as the sample size increases.

\subsubsection*{References}%
}

\DeclarePairedDelimiter\ceil{\lceil}{\rceil}
\DeclarePairedDelimiter\floor{\lfloor}{\rfloor}
\DeclareMathOperator{\tr}{tr}
\DeclareMathOperator{\vect}{vec}
\DeclareMathOperator*{\argmin}{\arg\!\min}

\newcommand{\genfrac{}{}{}1}{\genfrac{}{}{}1}
\newcommand{\mathbf{y}}{\mathbf{y}}
\newcommand{\mathbf{x}}{\mathbf{x}}
\newcommand{\mathbf{Y}}{\mathbf{Y}}
\newcommand{\mathbf{X}}{\mathbf{X}}
\newcommand{\mathbf{\Lambda}}{\mathbf{\Lambda}}
\newcommand{\mathbf{\Theta}}{\mathbf{\Theta}}
\newcommand{\mathbf{\Sigma}}{\mathbf{\Sigma}}
\newcommand{\mathbf{\Psi}}{\mathbf{\Psi}}
\newcommand{\mathbf{\Gamma}}{\mathbf{\Gamma}}

\newcommand{\mathbf{S_{xx}}}{\mathbf{S_{xx}}}
\newcommand{\mathbf{S_{xy}}}{\mathbf{S_{xy}}}
\newcommand{\mathbf{S_{yy}}}{\mathbf{S_{yy}}}

\newcommand{\mathbf{D}_\bLambda}{\mathbf{D}_\mathbf{\Lambda}}
\newcommand{\mathbf{D}_\bTheta}{\mathbf{D}_\mathbf{\Theta}}

\newcommand{\mathbf{U}}{\mathbf{U}}
\newcommand{\mathbf{V}}{\mathbf{V}}
\newcommand{\mathbf{R}}{\mathbf{R}}

\begin{document}
\runningtitle{Large-Scale Optimization Algorithms for Sparse CGGMs}
\twocolumn[
\aistatstitle{Large-Scale Optimization Algorithms for \\ Sparse Conditional Gaussian Graphical Models}

\aistatsauthor{ Calvin McCarter \And Seyoung Kim }

\aistatsaddress{ Machine Learning Department \\ Carnegie Mellon University \And Computational Biology Department \\ Carnegie Mellon University } ]

\begin{abstract}
This paper addresses the problem of scalable optimization for $l_1$-regularized
conditional Gaussian graphical models.  Conditional Gaussian graphical models
generalize the well-known Gaussian graphical models to conditional distributions
to model the output network influenced by conditioning input variables.
While highly scalable optimization methods exist for sparse
Gaussian graphical model estimation, state-of-the-art methods for
conditional Gaussian graphical models are not efficient enough
and more importantly, fail due to memory constraints for very large problems.
In this paper, we propose a new optimization procedure based on a Newton method
that efficiently iterates over two sub-problems, leading to drastic improvement
in computation time compared to the previous methods.
We then extend our method to scale to large problems under memory constraints,
using block coordinate descent to limit memory usage while achieving fast convergence.
Using synthetic and genomic data, we show that our methods can solve problems with millions of variables and tens of billions of parameters to high accuracy on a single machine.
\end{abstract}

\section{\uppercase{Introduction}}
\label{sec:introduction}

Sparse Gaussian graphical models (GGMs)~\citep{glasso} have been extremely popular
as a tool for learning a network structure over a large
number of continuous variables in many different application
domains including neuroscience~\citep{ng2011} and biology~\citep{glasso}.
A sparse GGM can be estimated as a
sparse inverse covariance matrix by
minimizing the convex function of $l_1$-regularized negative log-likelihood.
Highly scalable learning algorithms such as graphical lasso~\citep{glasso},
QUIC~\citep{Hsieh:2011}, and BigQUIC~\citep{Hsieh:2013}
have been proposed to learn the model.

In this paper, we address the problem of scaling up the optimization of sparse conditional GGM (CGGM),
a model closely related to sparse GGM, to very large problem sizes without requiring excessive time or memory.
Sparse CGGMs have been introduced as a discriminative extension of sparse GGMs
to model a sparse network over outputs conditional on input
variables~\cite{Sohn:2012,Wytock:2013}.
CGGMs can be viewed as a Gaussian analogue of conditional random field~\citep{crf},
while GGMs are a Gaussian analogue of Markov random field.
Sparse CGGMs have recently been applied to various settings
including energy forecasting~\citep{Wytock:2013:energy}, finance~\citep{yuan2014}, and biology~\citep{Zhang:2014},
where the goal is to predict structured outputs influenced by inputs.
A sparse CGGM can be estimated by minimizing a convex
function of $l_1$-regularized negative log-likelihood.
This optimization problem is closely related to
that for sparse GGMs because CGGMs also model the network over outputs.
However, the presence of the additional parameters in CGGMs for the functional mapping
from inputs to outputs makes the optimization significantly more complex than in sparse GGMs.

Several different approaches have been previously proposed to estimate
sparse CGGMs, including OWL-QN~\citep{Sohn:2012},
accelerated proximal gradient method~\citep{yuan2014},
and Newton coordinate descent algorithm~\citep{Wytock:2013}.
In particular, the Newton coordinate descent algorithm
extends the QUIC algorithm~\cite{Hsieh:2011} for sparse GGM estimation to the case of CGGMs,
and has been shown to have superior computational speed and convergence.
This approach finds in each iteration a descent direction
by minimizing a quadratic approximation of the original negative
log-likelihood function along with $l_1$ regularization.
Then, the parameter estimate is updated
with this descent direction and a step size found by line search.

Although the Newton coordinate descent method~\cite{Wytock:2013} is state-of-the-art
for its scalability and fast convergence, it is still not efficient enough
to be applied to many real-world problems even with tens
of thousands of variables. More importantly,
it suffers from a large space requirement, because
for very high-dimensional problems,
several large dense matrices need to be precomputed and stored during optimization.
For a CGGM with $p$ inputs and $q$ outputs,
the algorithm requires storing several $p \times p$ and $q \times q$ dense matrices,
which cannot fit in memory for large $p$ and $q$.

We propose new algorithms for learning $l_1$-regularized CGGMs
that significantly improve the computation time
of the previous Newton coordinate descent algorithm and also remove
the large memory requirement.
We first propose an optimization method,
called an alternating Newton coordinate descent algorithm,
for improving computation time.
Our algorithm is based on the key observation that
the computation simplifies drastically, if
we alternately optimize the two sets of
parameters for output network and for mapping inputs to outputs,
instead of updating all parameters at once as in the previous approach.
The previous approach updated all parameters simultaneously
by forming a second-order approximation of the objective on all parameters, which
requires an expensive computation of the large Hessian
matrix of size $(p+q)\times (p+q)$ in each iteration.
Our approach of alternate optimization forms a second-order approximation
only on the network parameters, which requires the Hessian of size $q \times q$,
as the other set of parameters can be updated easily using a simple coordinate descent.

In order to overcome the constraint on the space requirement,
we then extend our algorithm to an alternating Newton block coordinate descent method
that can be applied to problems of unbounded size on a machine with limited memory.
Instead of recomputing each element of the large matrices on demand,
we divide the parameters into blocks for block-wise updates such that the results of computation
can be reused within each block.
Block-wise parameter updates were previously used in BigQUIC~\citep{Hsieh:2013}
for learning a sparse GGM, where the block sparsity pattern of
the network parameters was leveraged to overcome the space limitations.
We propose an approach for block-wise update of the output network parameters in CGGMs that
extends their idea. We then propose a new block-wise update strategy
for the parameters for mapping inputs to outputs.
In our experiments, we show that we can solve problems with a million inputs and hundreds of
thousands of outputs on a single machine.

\iffalse
We first propose an optimization method for $l_1$-regularized CGGMs,
called an alternating Newton coordinate descent algorithm,
that significantly speeds up the previous Newton coordinate descent algorithm.
Our approach is based on the observation that
it is not necessary to form a second-order approximation
on all parameters as in the Newton coordinate descent algorithm,
because given the parameters for output network structure,
we can easily optimize over the other set of the parameters
for mapping inputs to outputs. Thus, instead of performing a descent
on all parameters simultaneously, we alternately optimize over
each of the two sets of parameters.
We show that our approach simplifies many aspects of the computation,
significantly reducing computation time.

We further extend our approach to handle the space requirement,
introducing an alternating Newton block coordinate descent method
that can be applied to problems of unbounded size
on a machine with limited memory.
Block-wise parameter updates were previously used in BigQUIC~\citep{Hsieh:2013}
for learning a sparse GGM, where the block sparsity pattern of
the network parameters
was leveraged to overcome the space limitations.
We adopt this key idea and propose a strategy for dividing the
parameters into blocks for block-wise updates such that the computation
can be grouped and reused within each block.
\fi

The rest of the paper is organized as follows.
In Section \ref{sec:background}, we provide a brief review
of sparse CGGMs and  the current state-of-the-art Newton coordinate descent
algorithm~\cite{Wytock:2013} for learning the models.
In Section \ref{sec:alternating-cd}, we propose
an alternating Newton coordinate descent algorithm that
significantly reduces computation time compared to the previous method.
In Section \ref{sec:alternating-bcd}, we further extend our algorithm
to perform block-wise updates in order to scale up to very large problems
on a machine with bounded memory.
In Section \ref{sec:experiments}, we demonstrate our proposed algorithms
on synthetic and real-world genomic data.

\iffalse
\subsection{Notation}
The number of non-zeros in $\Theta$ is $m_\Theta$ and for $\Lambda$ is $m_\Lambda$.
We assume that Conjugate Gradient (CG) method runs at most $C$ iterations.
We center and scale the data matrices so that their inner products equal their covariances.
Formally $\tilde{X} = \frac{1}{\sqrt{n}}(I_n - \frac{1}{n}\mathbf{1}\mathbf{1}^T)X$
and $\tilde{Y} = \frac{1}{\sqrt{n}}(I_n - \frac{1}{n}\mathbf{1}\mathbf{1}^T)Y$.
We define $\tilde{p}$ as the number of non-zero rows in $\Theta$.
\fi

\section{\uppercase{Background}}
\label{sec:background}

\subsection{The Conditional Gaussian Graphical Model}
A CGGM~\citep{Sohn:2012,Wytock:2013} models the conditional probability density of $\mathbf{x} \in \mathbb{R}^{p}$
given $\mathbf{y} \in \mathbb{R}^{q}$ as follows:
\begin{align*}
p(\mathbf{y}|\mathbf{x}; \mathbf{\Lambda},\mathbf{\Theta}) = \exp\{-\mathbf{y}^T\mathbf{\Lambda} \mathbf{y} - 2\mathbf{x}^T\mathbf{\Theta} \mathbf{y}\}/Z(\mathbf{x}),
\end{align*}
where $\mathbf{\Lambda}$ is a $q \times q$ matrix for modeling the network
over $\mathbf{y}$ and $\mathbf{\Theta}$ is a $p \times q$ matrix for modeling the mapping
between the input variables $\mathbf{x}$ and output variables $\mathbf{y}$.
\iffalse
The normalization constant is given as
$
        Z(\mathbf{x})
        = \int \exp
        \big(
                - \frac{1}{2} \, {\mathbf{y}}^T \mathbf{\Lambda} \mathbf{y}
                - {\mathbf{x}}^T \mathbf{\Theta} \mathbf{y}
        \big)
        d\mathbf{y}
        =\sqrt{\frac{(2\pi)^q}{|\mathbf{\Lambda}|}}
        \exp\big(\frac{1}{2} {\mathbf{x}}^T\mathbf{\Theta} \mathbf{\Lambda}^{-1} \mathbf{\Theta}^T\mathbf{x}\big)
$.
\fi
The normalization constant is given as
$Z(\mathbf{x}) = c |\mathbf{\Lambda}|^{-1} \exp(\mathbf{x}^T\mathbf{\Theta}\mathbf{\Lambda}^{-1}\mathbf{\Theta}^T\mathbf{x})$.
Inference in a CGGM gives
$
p(\mathbf{y}|\mathbf{x}) = \mathcal{N}(\mathbf{B} \mathbf{x},\mathbf{\Lambda}^{-1})
$, where $\mathbf{B}\!=\!-\mathbf{\Theta}\mathbf{\Lambda}^{-1}$,
showing the connection between a CGGM and multivariate linear regression.

Given a dataset of $\mathbf{X} \in \mathbb{R}^{n \times p}$
and $\mathbf{Y} \in \mathbb{R}^{n \times q}$ for $n$ samples,
and their covariance matrices
$\mathbf{S_{xx}} = \frac{1}{n}\mathbf{X}^T \mathbf{X}, \mathbf{S_{xy}} = \frac{1}{n}\mathbf{X}^T \mathbf{Y}, \mathbf{S_{yy}} = \frac{1}{n}\mathbf{Y}^T \mathbf{Y}$,
a sparse estimate of CGGM parameters can be obtained
by minimizing $l_1$-regularized negative log-likelihood:
\begin{align}
\min_{\mathbf{\Lambda} \succ 0,\mathbf{\Theta}} & f(\mathbf{\Lambda},\mathbf{\Theta})
= g(\mathbf{\Lambda},\mathbf{\Theta}) + h(\mathbf{\Lambda},\mathbf{\Theta}), \label{eq:opt}
\end{align}
where
$g(\mathbf{\Lambda}, \mathbf{\Theta}) \!\!=\!\! -\!\log|\mathbf{\Lambda}| \!+\!
\tr(\mathbf{S_{yy}}\mathbf{\Lambda} \!+\! 2\mathbf{S_{xy}}^T\mathbf{\Theta} \!+\!
\mathbf{\Lambda}^{-1}\mathbf{\Theta}^T \mathbf{S_{xx}} \mathbf{\Theta})$ for the smooth negative log-likelihood
and
$h(\mathbf{\Lambda},\mathbf{\Theta}) \!\!=\!\! \lambda_{\mathbf{\Lambda}} \|\mathbf{\Lambda}\|_1
+ \lambda_{\mathbf{\Theta}} \|\mathbf{\Theta}\|_1$
for the non-smooth elementwise $l_1$ penalty.
$\lambda_\mathbf{\Lambda}, \lambda_\mathbf{\Theta} > 0$ are regularization parameters.
\iffalse
We do not penalize the diagonal entries of $\mathbf{\Lambda}$,
following the common practice for sparse inverse covariance estimation.
\fi
As observed in previous work~\citep{Sohn:2012,yuan2014,Wytock:2013}, this objective is convex.

\subsection{Optimization}
The current state-of-the-art method for solving Eq. (\ref{eq:opt})
for $l_1$-regularized CGGM
is the Newton coordinate descent algorithm~\citep{Wytock:2013}
that extends QUIC~\citep{Hsieh:2011} for $l_1$-regularized GGM estimation.
In each iteration, this algorithm found a generalized Newton descent direction
by forming a second-order approximation of the smooth part of the objective
and minimizing this along with the $l_1$ penalty.
Given this Newton direction, the parameter estimates were updated
with a step size found by line search using Armijo's rule~\citep{Armijo:1966}.

In each iteration, the Newton coordinate descent algorithm
found the Newton direction as follows:
\begin{align}
\mathbf{D}_\bLambda, \mathbf{D}_\bTheta = \argmin_{\Delta_\mathbf{\Lambda},\Delta_\mathbf{\Theta}} & \
\bar{g}_{\mathbf{\Lambda},\mathbf{\Theta}}(\Delta_\mathbf{\Lambda},\Delta_\mathbf{\Theta}) \nonumber \\
& + h(\mathbf{\Lambda}\!+\!\Delta_\mathbf{\Lambda},\mathbf{\Theta}\!+\!\Delta_\mathbf{\Theta}),
\label{eq:nd}
\end{align}
where $\bar{g}_{\mathbf{\Lambda},\mathbf{\Theta}}$ is the second-order approximation of
$g$ given by Taylor expansion:
\begin{eqnarray*}
\bar{g}_{\mathbf{\Lambda},\mathbf{\Theta}}(\Delta_\mathbf{\Lambda},\Delta_\mathbf{\Theta})
= \vect(\nabla g(\mathbf{\Lambda},\mathbf{\Theta}))^T \vect([\Delta_\mathbf{\Lambda} \ \Delta_\mathbf{\Theta}])
\nonumber \\
\quad\quad +\frac{1}{2}\vect([\Delta_\mathbf{\Lambda} \ \Delta_\mathbf{\Theta}])^T
\nabla^2 g(\mathbf{\Lambda},\mathbf{\Theta})
\vect([\Delta_\mathbf{\Lambda} \ \Delta_\mathbf{\Theta}]).
\end{eqnarray*}
The gradient and Hessian matrices above are given as:
\begin{eqnarray}
\nabla g(\mathbf{\Lambda},\mathbf{\Theta}) =&
[\nabla_\mathbf{\Lambda} g(\mathbf{\Lambda},\mathbf{\Theta})\quad\quad  \nabla_\mathbf{\Theta} g(\mathbf{\Lambda},\mathbf{\Theta}) ]
\nonumber \\
=& [\mathbf{S_{yy}} \!-\! \mathbf{\Sigma} \!-\! \mathbf{\Psi} \quad 2\mathbf{S_{xy}} + 2\mathbf{\Gamma}]
\label{eq:grad} \\
\nabla^2 g(\mathbf{\Lambda},\mathbf{\Theta}) =&
\begin{bmatrix}
\nabla^2_\mathbf{\Lambda} g(\mathbf{\Lambda},\mathbf{\Theta}) & \nabla_\mathbf{\Lambda} \nabla_\mathbf{\Theta} g(\mathbf{\Lambda},\mathbf{\Theta}) \\
\nabla_\mathbf{\Lambda} \nabla_\mathbf{\Theta} g(\mathbf{\Lambda},\mathbf{\Theta})^T & \nabla^2_\mathbf{\Theta} g(\mathbf{\Lambda},\mathbf{\Theta})
\end{bmatrix}
\nonumber \\
=& \begin{bmatrix}
\mathbf{\Sigma} \otimes (\mathbf{\Sigma} + 2 \mathbf{\Psi})
& -2\mathbf{\Sigma} \otimes \mathbf{\Gamma}^T \\
-2\mathbf{\Sigma} \otimes \mathbf{\Gamma}
& 2\mathbf{\Sigma} \otimes \mathbf{S_{xx}}
\end{bmatrix}, \label{eq:hess}
\end{eqnarray}
where $\mathbf{\Sigma} = \mathbf{\Lambda}^{-1}$, $\mathbf{\Psi} = \mathbf{\Sigma}\mathbf{\Theta}^T \mathbf{S_{xx}} \mathbf{\Theta} \mathbf{\Sigma}$,
 and $\mathbf{\Gamma} = \mathbf{S_{xx}}\mathbf{\Theta}\mathbf{\Sigma}$.
Given the Newton direction in Eq. (\ref{eq:nd}), the parameters can be updated as
$\mathbf{\Lambda} \gets \mathbf{\Lambda} + \alpha \mathbf{D}_\bLambda$ and $\mathbf{\Theta} \gets \mathbf{\Theta} + \alpha \mathbf{D}_\bTheta$,
where step size $0 < \alpha \le 1$ ensures
sufficient decrease in Eq. (\ref{eq:opt}) and positive definiteness of $\mathbf{\Lambda}$.

The \textit{Lasso} problem~\citep{lasso} in Eq. (\ref{eq:nd})
was solved using coordinate descent.
Despite the efficiency of coordinate descent for \textit{Lasso},
applying coordinate updates repeatedly to all $q^2+pq$ variables in
$\mathbf{\Lambda}$ and $\mathbf{\Theta}$ is costly.
So, the updates were restricted to an active set of variables given as:
\begin{align}
\mathcal{S}_\mathbf{\Lambda} \ &= \quad \{(\Delta_\mathbf{\Lambda})_{ij} :
|(\nabla_\mathbf{\Lambda} g(\mathbf{\Lambda},\mathbf{\Theta}))_{ij}| > \lambda_\mathbf{\Lambda} \vee \mathbf{\Lambda}_{ij} \ne 0\}
\nonumber \\
\mathcal{S}_\mathbf{\Theta} \ &= \quad \{(\Delta_\mathbf{\Theta})_{ij} :
|(\nabla_\mathbf{\Theta} g(\mathbf{\Lambda},\mathbf{\Theta}))_{ij}| > \lambda_\mathbf{\Theta} \vee \mathbf{\Theta}_{ij} \ne 0\}.
\nonumber
\end{align}
Because the active set sizes
$m_\mathbf{\Lambda} = |\mathcal{S}_\mathbf{\Lambda}|, m_\mathbf{\Theta} = |\mathcal{S}_\mathbf{\Theta}|$
approach the number of non-zero entries in the sparse solutions for
$\mathbf{\Lambda}^*$ and $\mathbf{\Theta}^*$ over iterations,
this strategy yields a substantial speedup.

To further improve the efficiency of coordinate descent,
intermediate results were stored for
the large matrix products that need to be computed repeatedly.
At the beginning of the optimization for Eq. (\ref{eq:nd}),
$\mathbf{U} := \Delta_\mathbf{\Lambda} \Sigma$ and $\mathbf{V} := \Delta_\mathbf{\Theta} \Sigma$
were computed and stored.
Then, after a coordinate descent update to $(\Delta_\mathbf{\Lambda})_{ij}$,
row $i$ and $j$ of $\mathbf{U}$ were updated. Similarly,
after an update to $(\Delta_\mathbf{\Theta})_{ij}$, row $i$ of $\mathbf{V}$ was updated.

\iffalse
\begin{algorithm}[tb]
\caption{Newton Coordinate Descent}
\label{alg:Matt-alg}
   \DontPrintSemicolon
   \SetKwInOut{Input}{input}\SetKwInOut{Output}{output}
   \Input{Inputs $X \in \mathbb{R}^{n \times p}$
and outputs $Y \in \mathbb{R}^{n \times q}$;
regularization parameters $\lambda_{\mathbf{\Lambda}}, \lambda_{\mathbf{\Theta}}$}
   \Output{Parameters $\mathbf{\Lambda}, \mathbf{\Theta}$}
   Initialize $\mathbf{\Theta} \leftarrow 0$, $\mathbf{\Lambda} \leftarrow I_q$\;
   \For{$t=0,1,\dots$}{
       Compute active set $\mathcal{S} = (\mathcal{S}_\mathbf{\Lambda}, \mathcal{S}_\mathbf{\Theta})$\;
       Solve via coordinate descent:
       $D_\mathbf{\Lambda}, D_\mathbf{\Theta} =
       \underset{\Delta_\mathbf{\Lambda},\Delta_\mathbf{\Theta},\Delta_{\overline{\mathcal{S}}=0}}{\argmin}
       \bar{g}_{\mathbf{\Lambda},\mathbf{\Theta}}(\Delta_\mathbf{\Lambda},\Delta_\mathbf{\Theta})
       + h(\mathbf{\Lambda}\!+\!\Delta_\mathbf{\Lambda},\mathbf{\Theta}\!+\!\Delta_\mathbf{\Theta})$\;
   Update $\mathbf{\Lambda}^+ = \mathbf{\Lambda} + \alpha D_\mathbf{\Lambda}$ and $\mathbf{\Theta}^+ = \mathbf{\Theta}+\alpha D_\mathbf{\Lambda}$, where step size $\alpha$ is found by line search.\;
   }
\end{algorithm}
\fi

\subsection{Computational Complexity and Scalability}

Although the Newton coordinate descent method is computationally more efficient than
other previous approaches, it does not scale to problems even with tens of thousands of
variables. The main computational cost of the algorithm comes from computing
the large $(p+q) \times (p+q)$ Hessian matrix in Eq. (\ref{eq:hess})
in each application of Eq. (\ref{eq:nd}) to find the Newton direction.
At the beginning of the optimization in Eq. (\ref{eq:nd}),
large dense matrices $\mathbf{\Sigma}$, $\mathbf{\Psi}$, and $\mathbf{\Gamma}$, for computing
the gradient and Hessian in Eqs. (\ref{eq:grad}) and (\ref{eq:hess}), are precomputed
and reused throughout the coordinate descent iterations.
Initializing $\mathbf{\Sigma}=\mathbf{\Lambda}^{-1}$ via Cholesky decomposition costs up to $O(q^3)$ time,
although in practice, sparse Cholesky decomposition
exploits sparsity to invert $\mathbf{\Lambda}$ in much less than $O(q^3)$.
Computing $\mathbf{\Psi}=\mathbf{R}^T \mathbf{R}$, where $\mathbf{R}=\mathbf{X}\mathbf{\Theta}\mathbf{\Sigma}$,
requires $O(nm_\mathbf{\Theta} + nq^2)$ time,
and computing $\mathbf{\Gamma}$ costs $O(npq + nq^2)$.
After the initializations, the cost of coordinate descent update
per each active variable $(\Delta_\mathbf{\Lambda})_{ij}$ and $(\Delta_\mathbf{\Theta})_{ij}$ is $O(p+q)$.
During the coordinate descent for solving Eq. (\ref{eq:nd}),
the entire $(p+q) \times (p+q)$ Hessian matrix in Eq. (\ref{eq:hess})
needs to be evaluated, whereas for the gradient in Eq. (\ref{eq:grad})
only those entries corresponding to the parameters in active sets
are evaluated.

A more serious obstacle to scaling up to problems with large $p$ and $q$
is the space required to store dense matrices $\mathbf{\Sigma}$ (size $q\!\times\!q$),
$\mathbf{\Psi}$ (size $q\!\times\!q$), and $\mathbf{\Gamma}$ (size $p\!\times\!q$).
In our experiments on a machine with 104 Gb RAM,
the Newton coordinate descent method exhausted memory when $p+q$ exceeded $80{,}000$.

In the next section, we propose a modification
of the Newton coordinate descent algorithm
that significantly improves the computation time.
Then, we introduce block-wise update strategies to our algorithm
to remove the memory constraint and to scale to arbitrarily large problem sizes.

\section{\uppercase{Alternating Newton Coordinate Descent}}
\label{sec:alternating-cd}

In this section, we introduce our alternating Newton coordinate descent algorithm
for learning an $l_1$-regularized CGGM
that significantly reduces computation time compared to the previous method.
Instead of performing Newton descent
for all parameters $\mathbf{\Lambda}$ and $\mathbf{\Theta}$ simultaneously,
our approach alternately updates $\mathbf{\Lambda}$ and $\mathbf{\Theta}$,
optimizing Eq. (\ref{eq:opt}) over $\mathbf{\Lambda}$ given $\mathbf{\Theta}$
and vice versa until convergence.

Our approach is based on the key observation that with $\mathbf{\Lambda}$ fixed,
the problem of solving Eq. (\ref{eq:opt}) over $\mathbf{\Theta}$ becomes
simply minimizing a quadratic function with $l_1$ regularization.
Thus, it can be solved efficiently using a coordinate descent method,
without the need to form a second-order approximation or to perform line search.
On the other hand, optimizing Eq. (\ref{eq:opt}) for $\mathbf{\Lambda}$ given $\mathbf{\Theta}$ still requires forming a quadratic
approximation to find a generalized Newton direction and performing line search
to find the step size. However, this computation involves only $q \times q$ Hessian matrix and
is significantly simpler than
performing the same type of computation on both $\mathbf{\Lambda}$ and $\mathbf{\Theta}$ jointly
as in the previous approach.

\begin{algorithm}[tb]
\caption{Alternating Newton Coordinate Descent}
\label{alg:Fast-alg}
   \DontPrintSemicolon
   \SetKwInOut{Input}{input}\SetKwInOut{Output}{output}
   \Input{Inputs $X \in \mathbb{R}^{n \times p}$ and $Y \in \mathbb{R}^{n \times q}$;
regularization parameters $\lambda_{\mathbf{\Lambda}}, \lambda_{\mathbf{\Theta}}$}
   \Output{Parameters $\mathbf{\Lambda}, \mathbf{\Theta}$}
   Initialize $\mathbf{\Theta} \leftarrow 0, \mathbf{\Lambda} \leftarrow I_q$\;
   \For{$t=0,1,\dots$}{
     Determine active sets $\mathcal{S}_\mathbf{\Lambda}, \mathcal{S}_\mathbf{\Theta}$\;
     Solve via coordinate descent:\;
     \quad $D_\mathbf{\Lambda} \!=\!

     \underset{\Delta_\mathbf{\Lambda},\Delta_{\overline{\mathcal{S}_\mathbf{\Lambda}}=0}}{\argmin}
     \bar{g}_{\mathbf{\Lambda},\mathbf{\Theta}}(\mathbf{\Lambda}+\Delta_\mathbf{\Lambda},\!\mathbf{\Theta})
     + h(\mathbf{\Lambda}\!+\!\Delta_\mathbf{\Lambda},\!\mathbf{\Theta})$\;
     Update $\mathbf{\Lambda}^+ = \mathbf{\Lambda} + \alpha D_\mathbf{\Lambda}$,
     where step size $\alpha$ is found\;
     \quad with line search\;
     Solve via coordinate descent:\;
     \quad $\mathbf{\Theta}^+ = \arg\!\min_{\mathbf{\Theta}_{\mathcal{S}_\mathbf{\Theta}}} \ g_\mathbf{\Lambda}(\mathbf{\Theta}) + \lambda_{\mathbf{\Theta}} \|\mathbf{\Theta}\|_1$\;

   }
\end{algorithm}

\subsection{Coordinate Descent Optimization for $\mathbf{\Lambda}$}
Given fixed $\mathbf{\Theta}$,
the problem of minimizing the objective in Eq. (\ref{eq:opt})
with respect to $\mathbf{\Lambda}$ becomes
\begin{eqnarray}
\argmin_{\mathbf{\Lambda} \succ 0} g_\mathbf{\Theta}(\mathbf{\Lambda})
+ \lambda_\mathbf{\Lambda} \|\mathbf{\Lambda}\|_{1},
\label{eq:opt_lambda}
\end{eqnarray}
where $g_\mathbf{\Theta}(\mathbf{\Lambda})= -\log|\mathbf{\Lambda}| + \tr(\mathbf{S_{yy}}\mathbf{\Lambda}
+ \mathbf{\Lambda}^{-1} \mathbf{\Theta}^T \mathbf{S_{xx}} \mathbf{\Theta})$.
In order to solve this, we first find a generalized Newton direction
that minimizes the $l_1$-regularized quadratic
approximation of $g_\mathbf{\Theta}(\mathbf{\Lambda})$:
\begin{eqnarray}
\mathbf{D}_\bLambda = \argmin_{\Delta_\mathbf{\Lambda}}
\bar{g}_{\mathbf{\Lambda},\mathbf{\Theta}}(\Delta_\mathbf{\Lambda})
+ \lambda_\mathbf{\Lambda} \|\mathbf{\Lambda}+\Delta_\mathbf{\Lambda}\|_{1},
\label{eq:ancd_lambda}
\end{eqnarray}
where $\bar{g}_{\mathbf{\Lambda},\mathbf{\Theta}}(\Delta_\mathbf{\Lambda})$ is obtained from
a second-order Taylor expansion and is given as
\begin{eqnarray*}
\bar{g}_{\mathbf{\Lambda},\mathbf{\Theta}}(\Delta_\mathbf{\Lambda})
= \vect(\nabla_\mathbf{\Lambda} g(\mathbf{\Lambda}, \mathbf{\Theta}))^T \vect(\Delta_\mathbf{\Lambda})
\nonumber \\
+ \frac{1}{2}\vect(\Delta_\mathbf{\Lambda})^T
\nabla^2_\mathbf{\Lambda} g(\mathbf{\Lambda},\mathbf{\Theta}) \vect(\Delta_\mathbf{\Lambda}).
\end{eqnarray*}
The $\nabla_\mathbf{\Lambda} g(\mathbf{\Lambda}, \mathbf{\Theta})$ and
$\nabla^2_\mathbf{\Lambda} g(\mathbf{\Lambda},\mathbf{\Theta})$ above are components of the gradient
and Hessian matrices corresponding to $\mathbf{\Lambda}$ in Eqs. (\ref{eq:grad}) and (\ref{eq:hess}).
We solve the \textit{Lasso} problem in Eq. (\ref{eq:ancd_lambda}) via coordinate descent.
Similar to the Newton coordinate descent method,
we maintain $\mathbf{U} := \Delta_\mathbf{\Lambda} \mathbf{\Sigma}$ to reuse intermediate results
of the large matrix-matrix product.
Given the Newton direction for $\mathbf{\Lambda}$,
we update $\mathbf{\Lambda} \leftarrow \mathbf{\Lambda} + \alpha \Delta_\mathbf{\Lambda}$,
where $\alpha$ is obtained by line search.

Restricting the generalized Newton descent to $\mathbf{\Lambda}$
simplifies the computation significantly for coordinate descent updates,
compared to the previous approach \citep{Wytock:2013}
that applies it to both $\mathbf{\Lambda}$ and $\mathbf{\Theta}$ jointly.
Our updates only involve $\nabla_\mathbf{\Lambda} g(\mathbf{\Lambda},\mathbf{\Theta})$ and
$\nabla_\mathbf{\Lambda}^2 g(\mathbf{\Lambda},\mathbf{\Theta})$, and
no longer involve $\nabla_\mathbf{\Theta} g(\mathbf{\Lambda},\mathbf{\Theta})$ and
$\nabla_\mathbf{\Lambda}\nabla_\mathbf{\Theta} g(\mathbf{\Lambda},\mathbf{\Theta})$, eliminating
the need to compute the large $p \times q$ dense matrix $\mathbf{\Gamma}$ in $O(npq + nq^2)$ time.
Our approach also reduces the computational cost for the coordinate descent update
of each element of $\Delta_\mathbf{\Lambda}$ from $O(p+q)$ to $O(q)$.

\subsection{Coordinate Descent Optimization for $\mathbf{\Theta}$}
With $\mathbf{\Lambda}$ fixed, the optimization problem in Eq. (\ref{eq:opt}) with
respect to $\mathbf{\Theta}$ becomes
\begin{eqnarray}
\argmin_\mathbf{\Theta} \  g_\mathbf{\Lambda}(\mathbf{\Theta}) + \lambda_{\mathbf{\Theta}} \|\mathbf{\Theta}\|_1,
\label{eq:ancd_theta}
\end{eqnarray}
where $g_\mathbf{\Lambda}(\mathbf{\Theta}) = \tr(2\mathbf{S_{xy}}^T\mathbf{\Theta} + \mathbf{\Lambda}^{-1}\mathbf{\Theta}^T \mathbf{S_{xx}} \mathbf{\Theta})$.
Since $g_\mathbf{\Lambda}(\mathbf{\Theta})$ is a quadratic function itself,
there is no need to form its second-order Taylor expansion
or to determine a step size via line search.
Instead, we solve
Eq. (\ref{eq:ancd_theta}) directly with coordinate descent method,
storing and maintaining $\mathbf{V} := \mathbf{\Theta} \mathbf{\Sigma}$.
Our approach reduces the computation time for updating $\mathbf{\Theta}$
compared to the corresponding computation in the previous algorithm.
We avoid computing the large $p\!\times \! q$ matrix $\mathbf{\Gamma}$,
which had dominated overall computation time with $O(npq)$.
Our approach also eliminates the need for line search for updating $\mathbf{\Theta}$.
Finally, it reduces the cost
for each coordinate descent update in $\mathbf{\Theta}$ to
$O(p)$, compared to $O(p + q)$ for the corresponding computation
for $\Delta_{\mathbf{\Theta}}$ in the previous method.

Our approach is summarized in Algorithm \ref{alg:Fast-alg}.
We provide the details of the coordinate descent update equations in Appendix.
For approximately solving Eqs. (\ref{eq:ancd_lambda}) and (\ref{eq:ancd_theta}),
we warm-start $\mathbf{\Lambda}$ and $\mathbf{\Theta}$ from the results of the previous iteration
and make a single pass over the active set,
which ensures decrease in the objective in Eq. (\ref{eq:opt})
and reduces the overall computation time in practice.

\section{\uppercase{Alternating Newton Block Coordinate Descent}}
\label{sec:alternating-bcd}

The alternating Newton coordinate descent algorithm in the previous section
improves the computation time of the previous state-of-the-art method, but
is still limited by the space required to store large matrices during coordinate
descent computation. Solving Eq. (\ref{eq:ancd_lambda}) for updating $\mathbf{\Lambda}$
requires precomputing and storing $q \times q$ matrices, $\mathbf{\Sigma}$ and $\mathbf{\Psi}$,
whereas solving Eq. (\ref{eq:ancd_theta}) for updating $\mathbf{\Theta}$ requires
$\mathbf{\Sigma}$ and a $p \times p$ matrix $\mathbf{S_{xx}}$.
A naive approach to reduce the memory footprint would be to recompute
portions of these matrices on demand
for each coordinate update, which would be very expensive.

In this section we describe how our algorithm in the previous section
can be combined with block coordinate descent to scale up the optimization
to very large problems on a machine with limited memory.
During coordinate descent optimization, we update blocks of
$\mathbf{\Lambda}$ and $\mathbf{\Theta}$
so that within each block, the computation of the large matrices
can be cached and re-used, where
these blocks are determined automatically by exploiting the sparse stucture.
For $\mathbf{\Lambda}$,
we extend the block coordinate descent approach in BigQUIC \citep{Hsieh:2013}
developed for GGMs to take into account the conditioning variables
in CGGMs.
For $\mathbf{\Theta}$, we describe a new approach for block coordinate descent
update.
Our algorithm can, in principle, be applied to problems of any size on a
machine with limited memory.

\begin{algorithm}[tb]
\caption{Alternating Newton Block Coordinate Descent}
\label{alg:Huge-alg}
   \DontPrintSemicolon
   \SetKwInOut{Input}{input}\SetKwInOut{Output}{output}
   \SetKwComment{Comment}{$\triangleright$\ }{}
   \Input{Inputs $X \in \mathbb{R}^{n \times p}$
and outputs $Y \in \mathbb{R}^{n \times q}$;
regularization parameters $\lambda_{\mathbf{\Lambda}}, \lambda_{\mathbf{\Theta}}$}
   \Output{Parameters $\mathbf{\Lambda}, \mathbf{\Theta}$}
   Initialize $\mathbf{\Theta} \leftarrow 0, \mathbf{\Lambda} \leftarrow I_q$\;
   \For{$t=0,1,\dots$}{
     Determine active sets $\mathcal{S}_\mathbf{\Lambda}, \mathcal{S}_\mathbf{\Theta}$\;
     Partition columns of $\mathbf{\Lambda}$ into $k_\mathbf{\Lambda}$ blocks\;
     \Comment*[c]{Minimize over $\mathbf{\Lambda}$}
     Initialize $\Delta_\mathbf{\Lambda} \leftarrow 0$\;
     \For{$z=1$ \KwTo $k_\mathbf{\Lambda}$}{
       Compute $\mathbf{\Sigma}_{C_z}, \mathbf{U}_{C_z}$, and $\mathbf{\Psi}_{C_z}$\;
       \For{$r=1$ \KwTo $k_\mathbf{\Lambda}$}{
         \If{$z \ne r$}{
           Identify columns $B_{zr} \subset C_r$ with active\;
           \quad elements in $\mathbf{\Lambda}$\;
           Compute $\mathbf{\Sigma}_{B_{zr}}, \mathbf{U}_{B_{zr}}$, and $\mathbf{\Psi}_{B_{zr}}$\;
         }
         Update all active $(\Delta_\mathbf{\Lambda})_{ij}$ in $(C_z,C_r)$\;
       }
     }
     Update $\mathbf{\Lambda}^+ \leftarrow \mathbf{\Lambda} + \alpha\Delta_\mathbf{\Lambda}$,
     	where step size $\alpha$ is found\;
        \quad with line search.\;
     Partition columns of $\mathbf{\Theta}$ into $k_\mathbf{\Theta}$ blocks\;
     \Comment*[c]{Minimize over $\mathbf{\Theta}$}
     \For{$r=1$ \KwTo $k_\mathbf{\Theta}$}{
       Compute $\mathbf{\Sigma}_{C_r}$, and initialize $\mathbf{V} \leftarrow \mathbf{\Theta}\mathbf{\Sigma}_{C_r}$\;
       \For{{\upshape row} $i \in \{1,\dots,p\}$ {\upshape if} $I_\phi(\mathcal{S}_{(i,C_r)})$}{
         Compute $(\mathbf{S_{xx}})_{ij}$ for non-empty rows $j$ in $V_{C_r}$\;
         Update all active $\mathbf{\Theta}_{ij}$ in $(i,C_r)$\;
       }
     }
   }
\end{algorithm}

\begin{figure}[t]
\centering
\raisebox{-0.5\height}{\includegraphics[scale=1.2]{lambda_sweep.pdf}}
\ \ \
\raisebox{-0.5\height}{\includegraphics[scale=0.7]{lambda_figure_grid.pdf}}
\vspace{-10pt}
\caption{
Schematic of block coordinate descent for $\mathbf{\Lambda}$.
The $\mathbf{\Lambda}$ of size $q=9$ is updated for each of the $k_\mathbf{\Lambda}^2$ blocks in turn with $k_\mathbf{\Lambda}=3$.
Filled elements denote the parameters in the active set.
The green arrows denote rows of $\mathbf{\Sigma}$ and $\mathbf{\Psi}$ that are computed once and reused
while sweeping through a row of blocks. The red arrows denote cache misses and the corresponding columns of
$\mathbf{\Sigma}$ and $\mathbf{\Psi}$ need to be recomputed.
}
\vspace{-5pt}
\label{fig:lambda-schema-grid}
\end{figure}

\begin{figure}[t]
\centering
\begin{tabular}{c@{\hspace{25pt}}c@{\hspace{10pt}}c}
\raisebox{-0.5\height}{\includegraphics[scale=1.5]{theta_sweep.pdf}}
&
\raisebox{-0.5\height}{\includegraphics[scale=0.65]{theta_figure_grid1.pdf}}
&
\raisebox{-0.5\height}{\includegraphics[scale=0.65]{theta_figure_grid2.pdf}}
\end{tabular}
\vspace{-5pt}
\caption{
Schematic of block coordinate descent for $\mathbf{\Theta}$.  The $\mathbf{\Theta}$ of size $p\!=\!8,q\!=\!6$
is updated for each of the $p \times k_\mathbf{\Theta}$ blocks with $k_\mathbf{\Theta}\!=\!2$.
Filled elements denote the parameters in the active set.
Green arrows denote columns of $\mathbf{\Sigma}$ that are computed once and reused while sweeping
through the column of $p$ blocks. The red arrows denote cache misses for $(\mathbf{S_{xx}})_i$.
}
\vspace{-5pt}
\label{fig:theta-schema-grid}
\end{figure}

\iffalse
\begin{figure*}[t!]
\centering
\raisebox{-0.5\height}{\includegraphics[scale=1.0]{lambda_sweep.pdf}}
\qquad\quad
\raisebox{-0.5\height}{\includegraphics[scale=0.8]{lambda_figure.pdf}}
\vspace{-4pt}
\caption{Schematic of block coordinate descent for $\mathbf{\Lambda}$ where $q=9$ and $k_\mathbf{\Lambda}=3$.
Filled elements denote coordinates in the active set.
The green arrows denote columns of $\mathbf{\Sigma}$ and $\mathbf{\Psi}$ that are computed
before sweeping through a row of blocks, while the red arrows denote cache misses.
With 2 cache misses, columns of both $\mathbf{\Sigma}$ and $\mathbf{\Psi}$ are computed 11 times.
This compares to computing 18 columns from both $\mathbf{\Sigma}$ and $\mathbf{\Psi}$
without exploiting block-diagonal sparsity.
}
\label{fig:lambda-schema}
\end{figure*}
\begin{figure*}
\centering
\raisebox{-0.5\height}{\ \ \includegraphics[scale=1.0]{theta_sweep.pdf}}
\qquad\quad
\raisebox{-0.5\height}{\includegraphics[scale=0.75]{theta_figure.pdf}}
\vspace{-4pt}
\caption{
Schematic of block coordinate descent for $\mathbf{\Theta}$ where $p=8,q=6$, and $k_\mathbf{\Theta}=2$.
The green arrows denote columns of $\mathbf{\Sigma}$ that are computed before sweeping
through blocks of $\mathbf{\Theta}$, while red arrows denote cache misses for $(\mathbf{S_{xx}})_i$.
Columns of $\mathbf{S_{xx}}$ are computed only 7 times, and for each column only 5 out of 8 covariances are computed.
This compares to computing 128 covariances without exploiting row-wise sparsity.
Each column of $\mathbf{\Sigma}$ is computed once.
}
\label{fig:theta-schema}
\end{figure*}
\fi

\subsection{Blockwise optimization for $\mathbf{\Lambda}$}

\subsubsection{Block coordinate descent method}

A coordinate-descent update of $(\Delta_\mathbf{\Lambda})_{ij}$
requires the $i$th and $j$th columns of $\mathbf{\Sigma}$ and $\mathbf{\Psi}$.
If these columns are in memory, they can be
reused. Otherwise, it is a cache miss and we should compute them on demand.
$\mathbf{\Sigma}_i$ for the $i$th column of $\mathbf{\Sigma}$ can be obtained
by solving linear system $\mathbf{\Lambda} \mathbf{\Sigma}_i = \mathbf{e}_i$ with
conjugate gradient method in $O(m_\mathbf{\Lambda} K)$ time,
where $K$ is the number of conjugate gradient iterations.
Then, $\mathbf{\Psi}_i$ can be obtained from
$\mathbf{R}^T \mathbf{R}_i$, where $\mathbf{R} = \mathbf{X} \mathbf{\Theta} \mathbf{\Sigma}$
in $O(nq)$ time.

In order to reduce cache misses, we perform block coordinate descent,
where within each block, the columns of $\mathbf{\Sigma}$ are cached and re-used.
Suppose we partition $\mathcal{N} = \{1,\dots,q\}$ into $k_\mathbf{\Lambda}$ blocks,
$C_1,\dots,C_{k_\mathbf{\Lambda}}$. We apply this partitioning to the rows and columns of
$\Delta_\mathbf{\Lambda}$ to obtain $k_\mathbf{\Lambda} \times k_\mathbf{\Lambda}$ blocks.
We perform coordinate-descent updates in each block,
updating all elements in the active set within that block.
Let $\mathbf{A}_{C_r}$ denote a $q$ by $|C_r|$ matrix
containing columns of $\mathbf{A}$ that corresponds to the subset $C_r$.
In order to perform coordinate-descent updates on $(C_r, C_z)$ block of $\Delta_\mathbf{\Lambda}$,
we need $\mathbf{\Sigma}_{C_r}$, $\mathbf{\Sigma}_{C_z}$, $\mathbf{\Psi}_{C_r}$, and $\mathbf{\Psi}_{C_r}$.
Thus, we pick the smallest possible $k_\mathbf{\Lambda}$ such that we can store
$2q/k_\mathbf{\Lambda}$ columns of $\mathbf{\Sigma}$
and $2q/k_\mathbf{\Lambda}$ columns of $\mathbf{\Psi}$
in memory.
When updating the variables within block $(C_z, C_r)$ of $\Delta_\mathbf{\Lambda}$,
there are no cache misses
once $\mathbf{\Sigma}_{C_z}$, $\mathbf{\Sigma}_{C_z}$, $\mathbf{\Psi}_{C_z}$, and $\mathbf{\Psi}_{C_r}$
are computed and stored.
After updating each $(\Delta_\mathbf{\Lambda})_{ij}$ to $(\Delta_\mathbf{\Lambda})_{ij} + \mu$,
we maintain ${\mathbf{U}}_{C_z}$ and $\mathbf{U}_{C_r}$ by
$\mathbf{U}_{it} \leftarrow \mathbf{U}_{it} + \mu \mathbf{\Sigma}_{jt}, \mathbf{U}_{jt} \leftarrow
\mathbf{U}_{jt} + \mu\mathbf{\Sigma}_{it}, \forall t \in \{C_z \cup C_r\}$.

To go through all blocks, we update blocks $(C_z, C_1), \dots, (C_z, C_k)$
for each $z \in \{1,\dots,k_\mathbf{\Lambda}\}$.
Since all of these blocks share $\mathbf{\Sigma}_{C_z}$ and $\mathbf{\Psi}_{C_z}$,
we precompute and store them in memory.
When updating an off-diagonal block $(C_z, C_r), z \neq r$,
we compute $\mathbf{\Sigma}_{C_r}$ and $\mathbf{\Psi}_{C_r}$.
In the worst case, overall
$\mathbf{\Sigma}$ and $\mathbf{\Psi}$ will be computed $k_\mathbf{\Lambda}$ times.

\subsubsection{Reducing computational cost using graph clustering}
In typical real-world problems, the graph structure of $\mathbf{\Lambda}$
will exhibit clustering, with an approximately block diagonal structure.
We exploit this structure by choosing a partition
$\{C_1,\dots, C_{k_\mathbf{\Lambda}}\}$ that reduces cache misses.
Within diagonal blocks $(C_z, C_z)$'s, once $\mathbf{\Sigma}_{C_z}$ and $\mathbf{\Psi}_{C_z}$
are computed, there are no cache misses.
For off-diagonal blocks $(C_z, C_r)$'s, $r \neq z$,
we have a cache miss only if some variable in
$\{\Delta_{ij} | i \in C_z, j \in C_r \}$ lies in the active set.
We thus minimize the active set in off-diagonal blocks
via clustering, following the strategy for sparse GGM estimation
in \citep{Hsieh:2013}.
In the best case, if all parameters in the active set appear in the diagonal blocks,
$\mathbf{\Sigma}$ and $\mathbf{\Psi}$ are computed only once with no cache misses.
We use the METIS \citep{Karypis:1995} graph clustering library.
Our method for updating $\mathbf{\Lambda}$ is illustrated in Figure \ref{fig:lambda-schema-grid}.

\iffalse
The time complexity of each $\mathbf{\Lambda}$ update is dominated by
the cost of precomputing columns of $\mathbf{\Sigma}$ and $\mathbf{\Psi}$.
The cost of these precomputations is $O\Big(\big[1+\frac{B_\mathbf{\Lambda}}{q}\big][m_\mathbf{\Lambda} Kq + nq^2]\Big)$
where $B_\mathbf{\Lambda}=\sum_{z \ne r}{|\{j|i \in C_z, j \in C_r, (i,j) \in \mathcal{S}_\mathbf{\Lambda}\}|}$
is the number of cache misses.
Although the worst-case of $B_\mathbf{\Lambda} = k_\mathbf{\Lambda} q$
requires computing $\mathbf{\Sigma}$ and $\mathbf{\Psi}$ a total of $k_\mathbf{\Lambda}$ times,
in practice, graph clustering dramatically
reduces this additional cost of block-wise optimization.
In the best case,
when graph clustering identifies perfect block-diagonal structure
in the active set,
the number of cache misses $B_\mathbf{\Lambda}=0$ and we incur no runtime penalty
from limited memory.
\fi

\subsection{Blockwise Optimization for $\mathbf{\Theta}$}
\label{subsec:bcd-bTheta}

\subsubsection{Block coordinate descent method}

The coordinate descent update of $\mathbf{\Theta}_{ij}$
requires $(\mathbf{S_{xx}})_i$ and $\mathbf{\Sigma}_{j}$ to compute $(\mathbf{S_{xx}})_i^T \mathbf{V}_j$,
where $\mathbf{V}_{j} = \mathbf{\Theta} \mathbf{\Sigma}_{j}$.
If $(\mathbf{S_{xx}})_i$ and $\mathbf{\Sigma}_{j}$ are not already in the memory, it is a cache miss.
Computing $(\mathbf{S_{xx}})_i$ takes $O(np)$, which is expensive
if we have many cache misses.

We propose a block coordinate descent approach for solving Eq. (\ref{eq:ancd_theta})
that groups these computations to reduce cache misses.
Given a partition $\{1,\dots,q\}$ into
$k_\mathbf{\Theta}$ subsets, $C_1,\dots,C_{k_\mathbf{\Theta}}$,
we divide $\mathbf{\Theta}$ into $p \times k_\mathbf{\Theta}$ blocks,
where each block comprises a portion of a row of $\mathbf{\Theta}$.
We denote each block $(i, C_r)$, where $i \in \{1,\dots,p\}$.
Since updating block $(i,C_r)$ requires $(\mathbf{S_{xx}})_i$ and $\mathbf{\Sigma}_{C_r}$,
we pick smallest possible $k_\mathbf{\Theta}$ such that we can store
$q/k_\mathbf{\Theta}$ columns of $\mathbf{\Sigma}$ in memory.
While performing coordinate descent updates within block $(i,C_r)$ of $\mathbf{\Theta}$,
there are no cache misses, once $(\mathbf{S_{xx}})_i$ and $\mathbf{\Sigma}_{C_r}$ are in memory.
After updating each $\mathbf{\Theta}_{ij}$ to $\mathbf{\Theta}_{ij}+\mu$,
we update $\mathbf{V}_{C_r}$ by
$\mathbf{V}_{it} \leftarrow \mathbf{V}_{it} + \mu \mathbf{\Sigma}_{jt}, \forall t \in C_r$.

In order to sweep through all blocks,
each time we select a $q \in \{1,\dots,k_\mathbf{\Theta}\}$ and update blocks $(1,C_r), \dots, (p, C_r)$.
Since all of these $p$ blocks with the same $C_r$ share the computation of $\mathbf{\Sigma}_{C_r}$,
we compute and store $\mathbf{\Sigma}_{C_r}$ in memory.
Within each block, the computation of $(\mathbf{S_{xx}})_i$ is shared, so we pre-compute and store
it in memory, before updating this block.
The full matrix of $\mathbf{\Sigma}$ will be computed once while sweeping through the full $\mathbf{\Theta}$,
whereas in the worst case $\mathbf{S_{xx}}$ will be computed $k_\mathbf{\Theta}$ times.

\subsubsection{Reducing computational cost using row-wise sparsity}

We further reduce cache misses for $(\mathbf{S_{xx}})_i$
by strategically selecting partition $C_1,\dots,C_{k_\mathbf{\Theta}}$,
based on the observation that
if the active set is empty in block $(i,C_r)$,
we can skip this block and forgo computing $(\mathbf{S_{xx}})_i$.
We therefore choose a partition
where the active set variables are clustered into as few blocks as possible.
Formally, we want to minimize $\sum_{i,q}{|I_\phi(\mathcal{S}_{(i,C_r)})|}$,
where $I_\phi(\mathcal{S}_{(i,C_r)})$ is an indicator function that outputs
1 if the active set $\mathcal{S}_{(i,C_r)}$ within block $(i,C_r)$
is not empty
We therefore perform graph clustering over
the graph $G=(V,E)$ defined from the active set in $\mathbf{\Theta}$,
where $V = \{1, \ldots, q\}$ with one node for each column of $\mathbf{\Theta}$,
and $E = \{(j,k)|
\mathbf{\Theta}_{ij} \in \mathcal{S}_{\mathbf{\Theta}},
\mathbf{\Theta}_{ik} \in \mathcal{S}_\mathbf{\Theta}
\textrm{ for } i = 1, \ldots, p\}$,
connecting two nodes $j$ and $k$ with an edge
if both $\mathbf{\Theta}_{ij}$ and $\mathbf{\Theta}_{ik}$ are in the active set.
This edge set corresponds to the non-zero elements of $\mathbf{\Theta}^T \mathbf{\Theta}$,
so the graph can be computed quickly in $O(m_\mathbf{\Theta} q)$.

We also exploit row-wise sparsity in $\mathbf{\Theta}$ to
reduce the cost of each cache miss.
Every empty row in $\mathbf{\Theta}$ corresponds to an empty row in $\mathbf{V} = \mathbf{\Theta} \mathbf{\Sigma}$.
Because we only need elements in $(\mathbf{S_{xx}})_i$ for the dot product $(\mathbf{S_{xx}})^T_i \mathbf{V}_j$,
we skip computing the $k$th element of $(\mathbf{S_{xx}})_{i}$
if the $k$th row of $\mathbf{\Theta}$ is all zeros.
Our strategy for updating $\mathbf{\Theta}$ is illustrated in Figure \ref{fig:theta-schema-grid}.

Our method is summarized in Algorithm \ref{alg:Huge-alg}.
See Appendix for analysis of the computational cost.

\iffalse
The overall runtime is dominated by the cost of
precomputing columns of $\mathbf{S_{xx}}$ and $\mathbf{\Sigma}$.
The complexity of these operations is
$O(m_\mathbf{\Lambda} Kq\!+\!m_\mathbf{\Theta} q\!+\!n\tilde{p}B_\mathbf{\Theta})$,
where $\tilde{p}$ is the number of non-empty rows in $\mathbf{\Theta}$ and
$B_\mathbf{\Theta}=\sum_{i,r}{|\{i|j \in C_r, (i,j) \in \mathcal{S}_\Theta\}|}$
is the number of cache misses.
Without any row-wise sparsity we have $\tilde{p}\!=\!p$ and $B_\mathbf{\Theta}\!=\!k_\mathbf{\Theta} p$,
so the worst-case is that $\mathbf{S_{xx}}$ is computed a total of $k_\mathbf{\Theta}$ times.
This additional cost is substantially reduced in real datasets,
where most input variables influence few or none of the outputs.
\fi

\subsection{Parallelization}

The most expensive computations in our algorithm are embarrassingly parallelizable,
allowing for further speedups on machines with multiple cores.
Throughout the algorithm, we parallelize matrix and vector multiplications.
In addition, for block-wise $\mathbf{\Lambda}$ updates, we
compute multiple columns of $\mathbf{\Sigma}_{C_z}$ and $\mathbf{\Psi}_{C_z}$
as well as multiple columns of $\mathbf{\Sigma}_{C_r}$ and $\mathbf{\Psi}_{C_r}$ for multiple cache misses in parallel,
running multiple conjugate gradient methods in parallel.
For block-wise $\mathbf{\Theta}$ updates, we compute
multiple columns of $\mathbf{\Sigma}$ in parallel before sweeping through blocks
and perform a parallel compuation within each cache miss, computing elements within each $\mathbf{S_{xx}}_i$ in parallel.

\section{\uppercase{Experiments}}
\label{sec:experiments}
We compare the performance of our proposed methods with the existing state-of-the-art
Newton coordinate descent algorithm, using synthetic and real-world genomic
datasets. All methods were implemented in C++ with parameters represented in sparse matrix format.
All experiments were run on 2.6GHz Intel Xeon E5 machines with
8 cores and 104 Gb RAM, running Linux OS.
We run the Newton coordinate descent and alternating Newton coordinate descent
algorithms as a single thread job on a single core.
For our alternating Newton block coordinate descent method,
we run it on a single core and with parallelization on 8 cores.

\subsection{Synthetic Experiments}

\begin{figure*}
\centering
\begin{tabular}{@{}c@{}c@{}c}
\includegraphics[scale=0.47]{figures/scaling-PQ.pdf}
& \includegraphics[scale=0.47]{figures/scaling-iPQ.pdf}
& \includegraphics[scale=0.47]{figures/chain-convergence-40k-20k.pdf}
\vspace{-3pt}
\\
\vspace{-5pt}
(a) & (b) & (c)
\end{tabular}
\caption{Comparison of scalability on chain graphs.
We vary $p$ and $q$, where (a) $p=q$ and (b) $p=2q$.
The Newton coordinate descent and alternating Newton coordinate
descent methods could not be run beyond the problem sizes
shown due to memory constraint.
(c) Convergence when $q=20{,}000$ and $p=40{,}000$.
}
\label{fig:chain-results}
\vspace{-4pt}
\end{figure*}
\begin{figure*}
\centering
\begin{tabular}{@{}c@{}c@{}c}
\includegraphics[scale=0.47]{figures/scaling-P.pdf}
& \includegraphics[scale=0.47]{figures/scaling-Q.pdf}
& \includegraphics[scale=0.47]{figures/cluster-active-20k-10k.pdf}
\vspace{-4pt}
\\
\vspace{-4pt}
(a) & (b) & (c)
\end{tabular}
\vspace{-4pt}
\caption{Comparison of scalability on random graphs with clustering.
(a) Varying $p$ with $q$ fixed at $10{,}000$.
(b) Varying $q$ with $p$ fixed at $40{,}000$.
(c) Active set size versus time with $p=20{,}000$ and $q=10{,}000$.
}
\label{fig:cluster-results}
\vspace{-4pt}
\end{figure*}

We compare the different methods on two sets
of synthetic datasets, one for chain graphs
and another for random graphs with clustering for $\mathbf{\Lambda}$, generated as follows.
For chain graphs, the true sparse parameters $\mathbf{\Lambda}$ is set
with $\mathbf{\Lambda}_{i,i-1}=1$ and $\mathbf{\Lambda}_{i,i}=2.25$ and
the ground truth $\mathbf{\Theta}$ is set with $\mathbf{\Theta}_{i,i}=1$.
We perform one set of chain graph experiments where the number of inputs $p$
equals the number of outputs $q$,
and another set of experiments with an additional $q$ irrelevant features
unconnected to any outputs, so that $p$ equals $2q$.
For random graphs with clustering,
following the procedure in \citep{Hsieh:2013} for generating a GGM,
we set the true $\mathbf{\Lambda}$ to a graph with
clusters of nodes of size $250$ and with $90\%$ of edges connecting
randomly-selected nodes within clusters.
We set the number of edges so that
the average degree of each node is 10, with edge weights set to 1.
We then set the diagonal values
so that $\mathbf{\Lambda}$ is positive definite.
To set the sparse patterns for $\mathbf{\Theta}$, we randomly select
$100\sqrt{p}$ input variables as having edges to at least one output
and distribute total $10q$ edges among those selected inputs to influence
randomly selected outputs.
We set the edge weights of $\mathbf{\Theta}$ to 1.

Then, we draw samples from the CGGM defined by these true $\mathbf{\Lambda}$ and $\mathbf{\Theta}$.
We generate datasets with $n=100$ samples for the chain graphs and
$n=200$ samples for random graphs with clustering.
We choose $\lambda_\mathbf{\Lambda}$ and $\lambda_\mathbf{\Theta}$ so that the number of
estimated edges in $\mathbf{\Lambda}$ and $\mathbf{\Theta}$ is close to ground truth.
Following the strategy used in GGM estimation \citep{Hsieh:2013},
we use the minimum-norm subgradient of the objective as our stopping criterion:
$\|\text{grad}^S(\mathbf{\Lambda}_t,\mathbf{\Theta}_t)\|_1 < 0.01(\|\mathbf{\Lambda}\|_1 + \|\mathbf{\Theta}\|_1)$.

We compare the scalability of the different methods for chain graph experiments,
as we vary the problem sizes.
We show the computation times for datasets with $p=q$
in Figure \ref{fig:chain-results}(a) and for datasets with $p=2q$
with $q$ additional irrelevant features in Figure \ref{fig:chain-results}(b).
For large problems, computation times are not shown for
Newton coordinate descent and alternating Newton coordinate descent methods
because they could not complete the optimization with limited memory.
Also, for large problems, alternating block coordinate descent was terminated
after 60 hours of computation.
We provide more results on varying the sample size $n$ in Appendix.
In Figure \ref{fig:chain-results}(c), using the dataset with $p=40{,}000$ and $q=20{,}000$,
we plot the suboptimality in the objective $f - f^*$,
where $f^*$ is obtained by running alternating Newton
coordinate descent algorithm to numerical precision.
Our new methods converge substantially faster than the previous approach,
regardless of desired accuracy level.
We notice that as expected from the convexity of the optimization problem,
all algorithms converge to the global optimum and find nearly identical parameter estimates.

\begin{table*}[t!]
\caption{
Computation time in hours on genomic dataset.
`*' indicates running out of memory.
}
\centering
\label{table:real}
\begin{tabular}{@{}c@{\hskip 0.15in}c@{\hskip 0.15in}c@{\hskip 0.15in}c@{\hskip 0.15in}c@{\hskip 0.15in}c@{\hskip 0.15in}c}
\hline
$p$ & $q$ & $\|\mathbf{\Lambda}^*\|_0$ & $\|\mathbf{\Theta}^*\|_0$ & Newton CD & Alternating Newton CD & Alternating Newton BCD \\
\hline
34{,}249 & 3{,}268 & 34{,}914 & 28{,}848 & 22.0 & 0.51 & \textbf{0.24} \\
34{,}249 & 10{,}256 & 86{,}090 & 103{,}767 & $>50$ & 2.4 & \textbf{2.3} \\
442{,}440 & 3{,}268 & 26,232 & 30,482 & * & * & \textbf{11} \\
\hline
\vspace{-15pt}
\end{tabular}
\end{table*}

In Figure \ref{fig:cluster-results}, we compare scalability
of different methods for random graphs with clustering.
In Figure \ref{fig:cluster-results}(a), we vary $p$, while setting $q$ to $10{,}000$.
In Figure \ref{fig:cluster-results}(b),
we vary $q$, fixing $p$ to $40{,}000$.
Similar to the results from chain graph,
for larger problems,
Newton coordinate descent and alternating Newton coordinate descent methods
ran out of memory and alternating block coordinate descent was terminated
after 60 hours.
For all problem sizes, our alternating Newton coordinate descent algorithm
significantly reduces the computation time of the previous method,
the Newton coordinate descent algorithm. This gap in the computation time
increases for larger problems.
In Figure \ref{fig:cluster-results}(c),
we compare the convergence in sparsity pattern for the different methods
as measured by the active set size,
for problem size $p=20{,}000$ and $q=10{,}000$.
All our methods recover the optimal sparsity pattern
much more rapidly than the previous approach.

Figures \ref{fig:chain-results} and \ref{fig:cluster-results}
show that our alternating Newton block coordinate descent
can run on much larger problems than any other methods, while those methods without block
coordinate descent run out of memory.
For example, in Figure \ref{fig:cluster-results}(a) alternating Newton block coordinate
descent could handle problems with one million inputs,
while without block-wise optimization it ran out of memory when $p>100{,}000$.
We also notice that on a single core,
the alternating Newton block coordinate
\begin{figure}
\centering
\includegraphics[scale=0.45]{figures/compara-200-40000-20000.pdf}
\vspace{-9pt}
\caption{Speedup with parallelization for alternating Newton block coordinate descent.}
\vspace{-5pt}
\label{fig:parallel}
\end{figure}
descent is slighly slower
than the same method without block-wise
optimization because of
the need to recompute $\mathbf{\Sigma}$ and $\mathbf{S_{xx}}$. However, it is still
substantially faster than the previous method.

Finally, we evaluate the parallelization scheme for
our alternating Newton block coordinate descent method
on multi-core machines.
Given a dataset generated from chain graph with $p=40{,}000$ and $q=20{,}000$,
in Figure \ref{fig:parallel},
we show the folds of speedup for different numbers of cores with respect to a single core.
We obtained about 7 times speed up with 8 cores on a machine with 104Gb RAM,
and about 12 times speedup with 16 cores on a machine with 28Gb RAM.
In general, we observe greater speedup on larger problem sizes and
also for random graphs,
because such problems tend to have more cache misses that can be computed in parallel.

\subsection{Genomic Data Analysis}

We compare the different methods on a genomic dataset.
The dataset consists of genotypes for $442{,}440$ single nucleotide polymorphisms (SNPs)
and $10{,}256$ gene expression levels for 171 individuals with asthma,
after removing genes with variance under $0.01$.
We fit a sparse CGGM using SNPs as inputs
and expressions as outputs to model
a gene network influenced by SNPs.
We also compared the methods on a smaller dataset of
$34{,}249$ SNPs from chromosome 1 and
$3{,}268$ genes with variance $>0.1$.
As typically sparse model structures are of interests
in this type of analysis,
we chose regularization parameters
so that the number of non-zero entries in
each of $\mathbf{\Theta}$ and $\mathbf{\Lambda}$ at convergence
was approximately $10$ times the number of genes.

\iffalse
The computation time of different methods are provided in Table \ref{table:real}.
On the smallest problem with $p=34{,}249$ and $q=3{,}268$, both of our methods
are more than 50 times faster than the previous approach.
For the problem with $p=34{,}249$ and $q=10{,}256$,
both of our methods finished in under 3 hours,
whereas the previous approach did not converge in 50 hours.
On the large problem with $p=442{,}440$ and $q=3{,}268$,
our block coordinate descent approach converged in around 11 hours.
\fi

The computation time of different methods are provided in Table \ref{table:real}.
On the largest problem, the previous approach
could not run due to memory constraint, whereas our block coordinate
descent
\begin{figure}
\centering
\begin{tabular}{@{}c@{}}
\includegraphics[scale=0.47]{figures/real-convergence.pdf} \vspace{-5pt} \\
(a) \\
\vspace{-1pt} \\
\includegraphics[scale=0.47]{figures/real-active-time.pdf} \vspace{-5pt} \\
(b) \vspace{-6pt}
\end{tabular}
\caption{Convergence results on genomic dataset.
(a) Suboptimality and (b) active set size over time.
}
\label{fig:real-convergence}
\vspace{-4pt}
\end{figure}
converged in around 11 hours.
We also compare the convergence of
the different methods on the dataset
with $34{,}249$ SNPs and $3{,}268$ gene expressions in Figure \ref{fig:real-convergence},
and find that our methods provide vastly superior convergence
than the previous method.

\section{\uppercase{Conclusion}}
\label{sec:conclusion}
In this paper, we addressed the problem of large-scale optimization
for sparse CGGMs.
We proposed a new optimization procedure, called alternating Newton
coordinate descent, that reduces computation time
by alternately optimizing for the two sets of parameters $\mathbf{\Lambda}$ and $\mathbf{\Theta}$.
Further, we extended this with block-wise optimization
so that it can run on any machine with limited memory.

\iffalse
In this paper, we addressed the problem of large-scale optimization
for sparse CGGMs on a single machine.
The previously proposed Newton coordinate descent method is state-of-the-art
but scales poorly in both runtime and memory usage.
We proposed a new optimization procedure, called alternating Newton
coordinate descent, that reduces computation time
by alternately optimizing for the two sets of parameters $\mathbf{\Lambda}$ and $\mathbf{\Theta}$.
Further, we proposed the alternating Newton block coordinate descent method
that can run on any machine with limited memory.
Using both synthetic and genomic data, we demonstrated
that our methods achieve orders of magnitude speedup on large problems
and can handle ultra-large-scale problems
which were previously unsolvable on a single machine.
\fi

\newpage
\clearpage

redpajama_set_name: RedPajamaArXiv


text:
Brendan Joseph Byrne
priest theology educator
Lecturer, Corpus Christi College, Werribee, Australia, 1971; professor new testament united faculty theology, Jesuit Theological College, Parkville, Australia, since 1977; principal, Jesuit Theological College, Parkville, 1992-1997. Rector Campion College, Kew, Australia, 1979-1984. Vice president Melbourne College Division, since 1998.
src="/web/img/loading.gif" data-src="/web/show-photo.jpg?id=1100668&cache=false" alt="Other photo of Brendan Joseph Byrne" class="gallery__img" height="167"
Other photo of Brendan Joseph Byrne
data-src="/web/show-photo.jpg?id=1100668&cache=false" src="" alt="" class="gallery-slider__content__img" height=""
The Hospitality of God: A Reading of Luke's Gospel
( Luke portrays the life and ministry of Jesus as a divin...)
Luke portrays the life and ministry of Jesus as a divine visitation" to the world, seeking hospitality. The One who comes as visitor and guest becomes host and offers a hospitality in which the entire world can become truly human, be at home, and know salvation in the depths of their hearts. In The Hospitality of God Brendan Byrne, SJ, provides to-the-point commentary on those parts of Luke's Gospel that bring home to people a sense of the extravagance of God's love for them. The Hospitality of God approaches Luke's Gospel through the interpretive key of "hospitality." It looks at the Gospel as a whole, yet lingers upon scenes where the theme of "hospitality" is particularly prominent, such as the infancy stories, Jesus at Nazareth, Jesus in the house of Simon, the Good Samaritan, Martha and Mary, the banquet in 14:1-35, the Prodigal Son, Jesus' visit to the tax collector Zacchaeus, the institution of the Eucharist, and the Emmaus event. Byrne stresses that those in Luke's Gospel who readily offer hospitality - chiefly the marginalized and the poor - find themselves drawn into a much deeper sphere of hospitality, the hospitality of God. Those who have difficulty sharing are challenged by Jesus to conversion so that they, too, may not be left out of the banquet of life to which God calls all human beings. Luke's Gospel, ever interested in the process of human transformation, explores this resistance to God's gift and the ways in which individuals need to be converted if they are to come to "knowledge of salvation." Luke's Gospel is written out of faith in the risen Lord. Its primary intent is to bring people of al subsequent Christian generations into saving encounter with the Lord who was raised from the dead and lives among us in the power of the Spirit. Reading and hearing the Gospel in faith allows individuals and communities to access the hospitality of God brought by Jesus. In this work, Byrne offers an invitation for all to join in the life-giving "hospitality of God." Chapters are "The Prologue: 1:1-4 and the Lukan Time Frame of Salvation," "The Infancy Stories I: Before the Birth of Jesus: 1:5-80," "The Infancy Stories II: The Birth and Childhood of Jesus: 2:1-52," "Prelude to the Ministry of Jesus: 3:1-4:13," "Hospitality and Inhospitality at Nazareth: Jesus Inaugurates His Mission: 4:16-30," "The Early Galilean Ministry: 4:31-6:11," "The Community of the Kingdom: 6:12-49," "Response to the Ministry of Jesus: 7:1-8:3," "Later Galilean Ministry: 8:4-9:17," "Climax of the Galilean Ministry: 9:18-50," "The Journey to Jerusalem Begins: 9:51-10:24," "The Way to Eternal Life: 10:25-11:13," "The Prophet Continues His Way I: 11:14-12:53," "The Prophet Continues His Way II: 12:54-14:35," "Celebrating God's Acceptance: 15:1-32," "The Hospitality of the Poor: Two Responses to Wealth: 16:1-31," "Towards the End of the Journey: 17:1-18:30," "Hospitality in Jericho: 18:31-19:27," "Jesus in Jerusalem: 19:28-21:4," "Hope and Endurance: Discourse on the Future: 21:5-38," "Jesus' Journey to God: The Passion Narrative I: 22:1-53," "Jesus' Journey to God: The Passion Narrative II: 22:54-23:56,""Hospitality Continues: The Community of the Risen Lord: 24:1-53." Brendan Byrne, S.J., is professor of New Testament at Jesuit Theological College, within the United Faculty of Theology, Parkville, Melbourne, Australia. He is the author of the commentary Romans from the Sacra Pagina series, and Pal and the Christian Woman and is a contributor to the New Jerome Bible Handbook, published by The Liturgical Press. He was Catholic Biblical Association Visiting Professor at the Pontifical Biblical Institute, Rome, and has also taught in Africa and East Asia. Pope John Paul II named him a member of the Pontifical Biblical Commission in 1990. "
Luke portrays the life and ministry of Jesus as a divine "visitation" to the world, seeking hospitality. The One who comes as visitor and guest becomes host and offers a hospitality in which the entire world can become truly human, be at home, and know salvation in the depths of their hearts. In this new edition of The Hospitality of God, Brendan Byrne, SJ, provides to-the-point commentary on those parts of Luke's gospel that bring home to people a sense of the extravagance of God's love for them. The Hospitality of God approaches Luke's gospel through the interpretive key of "hospitality." It looks at the gospel as a whole, yet lingers upon scenes where the theme of hospitality is particularly prominent, such as the infancy stories, Jesus at Nazareth, Jesus in the house of Simon, the Good Samaritan, Martha and Mary, the banquet in 14:1-35, the Prodigal Son, Jesus' visit to the tax collector Zacchaeus, the institution of the Eucharist, and the Emmaus event. This revised edition includes a new introductory chapter titled "Features of Luke," which includes distinctive themes and motifs of the gospel that invite all to join in the life-giving hospitality of God. This revised edition is an excellent resource for students, teachers, and pastors interested in Luke's gospel or in search of a helpful perspective for pastoral and homiletic reflection and preparation.
Sons of God-Seed of Abraham: A Study of the Idea of the Sonship of God of All Christians in Paul Against the Jewish Background
Paul and the Christian Woman by Brendan S. J. Byrne (1989-08-03)
Reckoning with Romans: Contemporary Reading of Paul's Gospel
Oliver Rackham
Kenneth John Falconer
Meredith Brook Handspicker
John Bolt
Emanuel Tov
Brendan Joseph Byrne	Edit Profile
Brendan Joseph Byrne, Australian theology educator, priest. Joined Society of Jesus (Jesuit), 1957. Member Society for New Testament Studies, Catholic Biblical Association American, Australian Catholic Biblical Association (president 1997-1998).
Byrne, Brendan Joseph was born on October 12, 1939 in Melbourne, Victoria, Australia. Son of Francis Patrick and Mary Alice (Flanagan) Byrne.
Bachelor with honors, University Melbourne, Australia, 1966. Master of Arts, University Melbourne, 1972. Bachelor's Degree with honors, Melbourne College Divinity, 1972.
Doctor of Philosophy, Oxford (England) University, 1977.
Brendan Joseph Byrne has been listed as a noteworthy Theology educator, priest by Marquis Who's Who.
Member Society for New Testament Studies, Catholic Biblical Association American, Australian Catholic Biblical Association (president 1997-1998).
Avocations: reading, cycling.
Francis Patrick Byrne
Mary Alice (Flanagan) Byrne
University Melbourne , Bachelor of Arts with honors
University Melbourne , Master of Arts
Melbourne College Divinity , Bachelor's Degree with honors
Oxford (England) University , Doctor of Philosophy
Lecturer , Corpus Christi College
Werribee, Victoria, Australia
professor new testament united faculty theology , Jesuit Theological College
Parkville, Maryland, United States

redpajama_set_name: RedPajamaCommonCrawl


text:
February 8, 2019 10:34AM ET
Watch Ariana Grande Party in 'Break Up With Your Girlfriend I'm Bored' Video
Pop star dropped her 'Sweetener' follow-up 'Thank U, Next' early Friday morning
Brittany Spanos
@ohheybrittany	Follow
Brittany Spanos's Most Recent Stories
King Princess Announces Debut Album 'Cheap Queen,' New Song 'Prophet'
Watch Korn's Desert-Set, Sci-Fi 'You'll Never Find Me' Video
Ariana Grande makes a simple request in the "Break Up With Your Girlfriend I'm Bored" video. The singer dropped her Sweetener follow-up Thank U, Next last night.
Grande begins the video by entering a club and spotting a boy from across the dancefloor. Soon she sees his girlfriend, and ends up joining them on the floor. She begins being friendly with the girlfriend, playing with her hair and eventually following the couple to a house party. From there, she is a third wheel following them around, rolling her eyes and eventually joining them in a hot tub. At the end, she ends up going in for a kiss with the girlfriend instead of the beau you presume Grande is aiming to seduce.
"Break Up With Your Girlfriend I'm Bored" is the latest single from Thank U, Next and features a sample of 'N Sync's "It Makes Me Ill." It was preceded by the album's title track and the trap-pop "7 Rings," both of which were Number One singles.
Grande won't be attending or performing at the Grammys on Sunday. In an interview published on Thursday, the award show's producer, Ken Ehrlich, claimed Grande would be absent because she "felt it was too late for her to pull something together."
The singer took to Twitter on Thursday to correct the record. "I've kept my mouth shut but now you're lying about me," Grande wrote. "I can pull together a performance over night and you know that, Ken. It was when my creativity & self expression was stifled by you that I decided not to attend." She described the Grammys as "just a game" before adding, "I'm sorry, but that's not what music is to me."
In This Article: 'N Sync,	Ariana Grande


text:
Home Golf Tiger Woods Withdraws from the Ryder Cup
Tiger Woods Withdraws from the Ryder Cup
The Ryder Cup will not be the same this time around as Tiger Woods withdraws from the Ryder Cup as per news just coming in. The next edition of the Ryder Cup is going to be played in Scotland and with Tiger Woods withdrawal, the happiest person on the earth will be Tom Watson as now he has nothing to worry about his 'captain's pick' for the USA Ryder Cup team.
As per the statement issued and published by Tiger Woods himself on his official website, the doctors have advised Tiger Woods to take a rest because his back muscles need to be rehabilitated and healed. Woods has already talked over the phone with Watson on this issue as per this statement. He has also assured the USA Golf fans that he will cheer the US team during the tournament.
It is going to be a long wait for the Tiger Woods fans as he has declared that he will return to the professional golf only on the 1st of December this year when he will play the World Challenge Tournament at Isleworth in Orlando, Florida.
There was a rumor prior to this announcement making rounds for the last couple of weeks that Tom Watson was considering other players then Woods for the Ryder Cup team and a few of them are even unknown. It would have been very embarrassing for Woods had that happened. So does this mean that to avoid the embarrassment Tiger Woods has withdrawn his name for the consideration in the USA Ryder Cup team?
Well, it doesn't seem to be the case as he has declared his return in December and not immediately after the Ryder Cup gets over. Questions are also arising that if Tiger Woods is undoubtedly not fit, then why he was allowed to play the Bridgestone Invitational Tournament and the PGA Championship which were played side by side in last fortnight! Were his doctors sleeping?
Something, somewhere is definitely missing here, but then this decision by Woods will only help the USA Ryder Cup team as it has avoided a controversy which was looming large on them.
The 2014 Ryder Cup will be the 40th Edition and it will be hosted by PGA Centenary Course at the Gleneagles Hotel near Auchterarder in the Perthshire of Scotland. The tournament be played between 26th September and 28th September.
Previous articleIs Dhoni becoming the worst captain ever overseas?
Next articleTiger and Mcllory to appear together on The Tonight Show
The Rod Laver Arena, Where Australian Open is played
Five cricketers who can hit the six at will
Player Profile – Lasith Malinga
Best Surfing destinations: Check Top 5 Surfing destinations of the world
Your Quick Guide To Group F Teams of the FIFA World Cup 2014
Tiger Woods to Design New Golf Course in Dubai
RR Stuns RCB In Their Own Backyard
Player Profile – Dale Steyn
Important Team Tournaments in Golf
All about Carrom Board Game


These Forgefix Modesty Blocks are versatile joining blocks. Manufactured from durable polyamide.
Suitable for 6 and 8 gauge (3.5 - 4.5 mm) screws.
Used when joining timber and carcasses assembly in many joinery applications. Widely used in flat-pack furniture assembly.


text:
The Foundation supports healthy families through traditional grantmaking and New Haven Healthy Start, a program of The Foundation and supported by a federal grant from Health Resources and Services Administration (HRSA) Maternal and Child Health.
Giving Babies a Healthy Start
New Haven Healthy Start is a program of The Community Foundation for Greater New Haven and an outgrowth of the work of the Commission on Child and Infant Health, convened by The Foundation in the mid 1980s. It is one of two Healthy Start programs housed at a community foundation in the country.
New Haven Healthy Start Programs
Back to School Strategies and Techniques ›
Tools that help parents and children center, ground and get in touch with feelings
New Haven Healthy Start RFP 2020 ›
This request for proposal process has ended and a vendor has been selected
Family Support Magistrates Offer Guidance for Court-involved Fathers ›
New Haven Healthy Start hosts a judicial panel to help court-involved fathers and families better understand how the court system works and advocate for themselves.
Healthy Families Need Healthy Fathers
While health care spending is at an all-time high, the health system is failing young African American men. A March, 2015 commentary published in The Journal of American Medical Association co-authored by former New Haven Healthy Start Director Kenn Harris and Boston Medical Center research physicians, offers analysis and recommendations for how the medical field can better meet the needs of African American men.
During Black History Month, we celebrate the many people who are committed to New Haven Healthy Start and to securing a strong Maternal and Child Health system in our community and beyond. Learn more about these local champions and their outstanding work below #blackchampions4health:
Black Champions For Health: Dr. Marcella Nunez-Smith ›
Adult Primary Care Physician and Researcher at Yale Schools of Medicine and Public Health.
Black Champions For Health: Whitney Woodhouse ›
Care Coordinator for New Haven Healthy Start at Cornell Scott Hill Health Center.
Black Champions For Health: Dr. Byron Kennedy ›
Director of Health at the New Haven Health Department.
Natasha Ray, MS
Program Director New Haven Healthy Start
Celebrating Latino Leaders: Yolanda Caldera-Durant
"An important part of my continued success is giving back to young people, helping them understand that they deserve success as much as anyone else."
Celebrating Latino Leaders: Jordy Padilla
"It's a duty of mine, as an educated individual, to fight for the people who didn't have the same opportunity as me."
Celebrating Latino Leaders: Ricardo Henriquez
"You do not always have to understand what someone is asking you to do if there is trust in the middle."


In this project we integrated payments of Discount offers website with Stripe Payment Gateway. We built a module that facilitates payment through Stripe Payment Gateway and parses the response received backThis is an eCommerce portal built over popular Learning Management Systems like Moodle, Totara, Bridge and Blackboard in Heroku infrastructure in CakePHP. We have integrated it with Paypal payment gateway to facilitate purchases.


After a couple months of topics for our South Asian women in sports community I decided to change it up and bring a different format to the chat. An individual who is in the industry and well versed for a number of years was brought in to bring a different dynamic. In addition, it gave folks to connect with yet another South Asian female in the industry. She's been working in the sports industry for almost 10 years!
April's chat was with Pooja Kobawala Van Dyke who is a Senior Manager at ESPN CreativeWorks. The conversation started off with learning more about who and what ESPN CreativeWorks is and does.
The path Pooja took to get to where she is today wouldn't seem so "traditional" in a South Asian household. She has been in the sports industry since college, started off by working for the school paper.
And when asked what sent her down this path, she highlighted something I thought was very important, having family support. It's very tough for South Asians to explain why they want to pursue a career in this field. Though, the more successful "Poojas" we have in the world the more success we'll have as a community.
Not only did her success come from having support at home but also having a platform to prepare her for the real world and learn about the barriers minorities deal with. That platform was the South Asian Journalists Association.
There is also the very powerful espnW Summit hosted every year for the last four years where females in the sports industry come together to share stories, ideas, and empower one another to become successful. This chat is doing the same on a smaller scale, it's about looking after one another, sharing experiences and ideas on a platform everyone can go to.
And finally, Pooja talked about how she attended SXSW this year and how creative the conference has become.
As these chats continue the idea is to get many folks engaged, sharing stories and ideas with one another, and most importantly for everyone to get to know what each other is doing in the industry.
It was delightful chatting with Pooja about her thoughts and experiences in the sports industry! The next chat is scheduled for Wednesday, May 14 from 9 to 10 PM EST. Join in using #SAWSports!
Next > Chat Recap – Can women have it all? A sports perspective.


VGP – A total number of 8,465 firms re-operated in January, 2019, up 84.5% compared to the same period last year, according to the General Statistics Office.
Retails, sales take the lead in the number of new firms in January, 2019.
As many as 10,079 enterprises were established in the first month of 2019 with the total registered capital of VND151,100 billion, down 7% in the number and up 53.8% in capital against the same period last year. The average capital for each enterprise was estimated at VND15 billion, up 64.8%.
As much as VND484,000 billion added to the economy in January, raising the total registered and additional capital to VND635,100 billion.
The number of newly-established enterprises in all regions nationwide decreases compared to the same period last year, including Southeast (down 11.8%) and Red River Delta (down 9.8%).
The total laborers working for new firms in the month hit 107,900, up 26.5%.
In January, the number of dissolved enterprises reached 1,802, up 16%.


text:
Casey by one by Arnold in the hunt
Cronulla's Scott Arnold is two shots back from Englishman Paul Casey after a morning of hot scoring on the second day of the Emirates Australian Open.
Arnold carded a tidy four under 67 to go with his opening 68 in the milder morning conditions.
Casey, the highest-ranked player in the field, sits at nine under par for the tournament after taking advantage of vastly improved conditions for the second round.
Scott Arnold stares down a shot during round two.
The world No.14 bagged seven birdies – including four straight to close his round – in a six-under-par 65 to vault to the top of a jam-packed leaderboard on Friday.
"I've not been able to get it close to any flag for about a day and a half and then suddenly the back nine for me, the last half a dozen holes, I had some really good birdie opportunities," Casey said.
"I don't know if it was fully in the zone but it was really nice stuff and it felt really good.
"So suddenly from a case of sort of bumbling around, just a couple under par, challenging for lead, probably in the lead right now; so good stuff."
Cameron Tringale also collected seven birdies in a 65 that could have been even better had he not made a mess of the par-5 fifth hole, where the American took bogey.
Former British Open Champion Louis Oosthuizen was another to make a major move, with the South African climbing from three to eight to be one stroke adrift of Casey.
Daniel Ieremia plays his approach to the last.
Kiwi Daniel Ieremia was the leading Australasian, also at eight under, after posting an equal-morning-best 65.
Overnight joint leader Chun-An Yu was just two shots off the pace after carding a second-round 69.
Former Champion Matt Jones has made a bright start to his second round, adding two birdies to his overnight total to be at six-under. Dimi Papadatos has also grabbed two early birdies to go to six under alongside Jones.
Adam Scott clawed his way back with a four under 67 in round two.
After being in danger of missing the cut following his opening 75, pre-tournament favourite Adam Scott clawed his way back to even par with a round of four-under 67 to all but ensure himself a weekend start.
"It was much, much better than yesterday," Scott said.
"Hopefully I didn't shoot myself in the foot with a bad round yesterday and I'd like to play two more rounds and keep working on this.
"You never know, I can have a decent result here and then be in good form going into (the Presidents Cup) next week and the (Australian) PGA.
Friday's low-scoring came after the bushfire-driven haze that blanketed The Australian Golf Club on Thursday afternoon cleared significantly for the second-round early starters.
Australian Open Golf, Paul Casey, Scott Arnold
© Golf NSW | ALL RIGHTS RESERVED


Luciano Frigerio was born in Desio in 1928.
Her father, John begins, since 1889, a workshop crafts woodworking.
Over the years, he collaborated with architects and famous Italian designers such as Franco Albini, Sergio Asti and later Paolo Portoghesi.
Permanent exhibitions in Milan Gallery in the cemetery 1959-1963, Galleria Vittorio Emanuele permanent in 1964 * 1967 Permanent exhibition "Art Gallery" Via Montenapoleone 1968-1989.
He died in Sanremo in 1999.
Luciano Frigerio, pair of table lamps, bronze, Frigerio di Desio production, signed under the base, circa 1970, Italy. Measures: Height 46 cm, diameter 21 cm.
Frigerio, pair of table lamps, gilded bronze, circa 1970, Italy. Measure: Height 60 cm, width 12 cm, depth 12 cm.
Piedestal by Frigerio Di Desio (1928 - 1999) Piedestal.
Frigerio Di Desio ( 1928 - 1999 ) Piedestal,lacquered wood and glass top,Italy, circa 1970. Height: 75 cm, diameter: 130 cm.


Satan is alive and active on planet earth – and George Soros is doing his will, destroying one patriotic nation after another, according to a Hungarian politician.
The mega-billionaire's war on countries that refuse to accept millions of Muslim immigrants from the Middle East is really "Satan's Soros Plan," he said.
Andras Aradszki tore the lid off the plot on the floor of Hungary's national parliament in Budapest.
"Yes, I speak of an attack by Satan," said Aradszki, the nation's energy undersecretary.
Aradszki, who is a member of the Christian Democratic People's Party, said George Soros wants the European Union to siphon away the national sovereignty of its member nations, bribe them – or force them – to accept waves of Islamist "migrants," and destroy the traditional Christian civilization of each and every country in the world.
Soros admitted parts of the "Soros Plan" in an op-ed in 2015, calling for Europe to do nothing about the uncontrolled flood of Islamist immigrants. EU nations would just receive these people, and their languages and customs, whether they like it or not.
"The EU has to accept at least a million (Muslim) asylum-seekers annually for the foreseeable future," Soros wrote.
He said the EU needed to develop a "comprehensive plan" that will "reassert effective governance" over its 28 member nations.
The most important thing, he wrote, is that the EU "put an end to the panic," because too many nations were opposed to being flooded with Muslim fundamentalists coming from terrorist-sponsors.
Why not? The jet-setting Soros can always flee to another one of his mansions, behind massive gates secured by teams of armed guards.
This this is all part of his plan to destroy traditional religious culture and create a new world government, Aradszki stated.
"Soros and his comrades want to destroy the independence and values of nation states for the purpose of watering down the Christian spirit of Europe with the forced settlement of tens of millions of migrants," Aradszki said.
Aradszki said, though he was in Hungary at the time, he prayed in his heart with millions of Polish people during the "Rosary on the Borders." Millions circled Poland, surrounded it in prayer, and made a statement that their nation was closed to Islamic terrorism.
Until then, he praised the government for mailing response forms to every home in Hungary, asking their opinion on "the Soros Plan" earlier this month.
Hungary has come out swinging against Soros' financial – and political – war against the West.
Prime Minister Viktor Orban is Soros' Public Enemy Number One. He has exposed Soros' plans and massive influence on the governments of national, and international, bodies. He has opposed Soros opening his own university in Budapest.
Hungary's deputy prime minister said that Soros' war is part of a long, proto-Communist hatred of God, freedom and patriotism.
"Fundamentally, the root of this hundreds of years ago started from Freemasonary's inspiration, which later had a Jacobin version and a Bolshevik version, and one of its many branches is the Soros-type extreme-liberal thing," said Zsolt Semjén on national radio the day after Aradszki's speech.
Media outlets have dismissed these warnings about Satan's plan to import terrorists as anti-Semitic hatred.
Soros, who was born Jewish, allegedly collaborated with the Nazis before fleeing his home country.
A Soros spokesperson told the Financial Times that normal European citizens have nothing to worry about. "There is no such thing as a global conspiracy against Hungary orchestrated by George Soros," he said.
Soros has made a fortune betting on the stock market – but the Hungarian people aren't buying his bull on this.
"We want a Hungarian Hungary and a European Europe," said Orban.
Hungarians overwhelmingly want to preserve their own culture, language, and borders – and they're reward politicians who are willing to do something about it. More than 80 percent of the Hungarian people support the two toughest parties on border security: Viktor Orban's Fidesz Party and the even more nationalist Jobbik Party. Orban seems to be a shoo-in to win a third term as prime minister in April.
Meanwhile, Soros is ginning up his operations in the U.S. He refilled his Open Society Foundation with $18 billion of his personal wealth to wage political warfare against America.
Maybe it isn't just political – maybe it truly is demonic.
Be sure to read The Horn News as we expose Soros' anti-American political maneuvers around the world.
And say a prayer for America this Sunday – and everyday – that Soros' money won't convince us to destroy our country's foundations.


Designed in a micro floral print, this top showcases a ruffled button-front and Victorian-inspired puffed sleeves.
Designed in a micro floral print, this top showcases a ruffled button-front and Victorian-inspired puffed sleeves. The soft cotton silhouette pairs especially well with denim.


Lead Grey Roofers' Seal is a single component, high performance elastomeric sealant for most roofing applications and repairs.
This lead grey coloured sealant is suitable for most roofing applications and repairs; it adheres to various roofing substrates in both wet and dry conditions. It is ideal for sealing around gutters, downpipes, lead flashing and corrugated roofing.