11affiliationtext:  Department of Statistics and Applied Probability,
University of California, Santa Barbara, Santa Barbara, CA22affiliationtext:
Department of Statistics, Oregon State University, Corvallis, OR

# Sparse estimation of parameter support sets for generalized vector
autoregressions by resampling and model aggregation

Trevor D. Ruiz Corresponding Author<EMAIL_ADDRESS>Sharmodeep Bhattacharyya
and Sarah C. Emerson

###### Abstract

The central problem we address in this work is estimation of the parameter
support set $S$, the set of indices corresponding to nonzero parameters, in
the context of a sparse parametric likelihood model for count-valued
multivariate time series. We develop a computationally-intensive algorithm
that performs the estimation by aggregating support sets obtained by applying
the LASSO to data subsamples. Our approach is to identify several well-fitting
candidate models and estimate $S$ by the most frequently-used parameters, thus
aggregating candidate models rather than selecting a single candidate deemed
optimal in some sense. While our method is broadly applicable to any selection
problem, we focus on the generalized vector autoregressive model class, and in
particular the Poisson case, due to (i) the difficulty of the support
estimation problem due to complex dependence in the data, (ii) recent work
applying the LASSO in this context, and (iii) interesting applications in
network recovery from discrete multivariate time series. We establish
benchmark methods based on the LASSO and present empirical results
demonstrating the superior performance of our method. Additionally, we present
an application estimating ecological interaction networks from
paleoclimatology data.

Keywords: model selection; computationally-intensive methods; LASSO; vector
time series

##  1 Introduction

Vector autoregressive (VAR) models describe the serial dependence in
multivariate time series, including both inter- and intra-series correlations,
but are densely parametrized: the number of parameters scales quadratically
with the process dimension. Bach and Jordan (2004); Arnold et al. (2007) and
Hsu et al. (2008) used the LASSO (Tibshirani, 1996) to apply VAR models to
high-dimensional data, and more recently Basu et al. (2015) established
theoretical support for the method. Variants on the LASSO specifically adapted
to the VAR setting were introduced by Shojaie and Michailidis (2010); Song and
Bickel (2011), and Fan et al. (2011), and alternative sparse estimation
methods that sought to reduce overfitting were proposed by Han et al. (2015)
and Davis et al. (2016). Among applications, a central theme in estimating VAR
models with sparsity is learning causal relationships between the
observational units corresponding to the components of the multivariate
process (Valdés-Sosa et al., 2005; Arnold et al., 2007; Lozano et al., 2009).
Causality is narrowly defined in this context based on forecasting accuracy in
the sense proposed by Granger (1969): one process is causal for another if the
former improves forecasts of the latter. Granger causality can be identified
directly from the support set of VAR model parameters, a fact which underpins
the explicitly causal interpretation of a VAR process as a graphical model
(Dahlhaus and Eichler, 2003). Thus, in practice, sparse estimation is utilized
specifically for identifying a causal network.

Network estimation is also a central application of generalized vector
autoregressive (GVAR) models for discrete data (Mark et al., 2017; Hall et
al., 2018). These models — discrete-time stochastic processes
$\\{X_{t}\\}_{t\in\mathbb{Z}}$ specified through a parametric model for
$X_{t}|\mathcal{H}_{t}$ in the exponential family in which the parameter is a
function of the process history $\mathcal{H}_{t}$ — provide a VAR-like
framework for non-Gaussian and discrete data (Cox et al., 1981; Fahrmeir and
Tutz, 1994; Karlis, 2016; Fokianos, 2021). For a subclass of GVAR models,
LASSO estimates exhibit consistency and support recovery properties under
certain conditions (Pandit et al., 2020; Han et al., 2023), and can be
computed efficiently using methods for generalized linear models (Friedman et
al., 2010). Although LASSO estimation of GVAR models has only recently
attracted attention in the statistical community, it has been in use in
applications for some time in the sciences (Brown et al., 2004; Pillow et al.,
2008; Paul et al., 2008).

The parameter support set itself, apart from the values of model parameters,
is thus an object of particular interest in the application of GVAR models,
and consequently, the selection behavior of sparse estimators in practice with
limited data is of special importance. However, the conditions under which
theoretical guarantees hold are difficult to verify in practice, and moreover,
the use of ad hoc approaches for fixing regularization hyperparameters can
severely compromise performance. There is little guidance available on how to
directly resolve these challenges.

In other contexts, data resampling has been used to formulate selection and
sparse estimation methods that achieve strong empirical performance.
Resampling allows for the possibility of obtaining collections of models with
comparable sparsity levels that can be considered as candidates and leveraged
in various ways to formulate more robust strategies for searching large model
spaces. For example, if many resamples are considered, candidate models can be
used to estimate the frequency with which parameters are selected across data
subsets at different sparsity levels. Bach (2008) used these frequencies to
reduce overfitting across the regularization path by returning only common
parameter supports across all resamples. (Meinshausen and Bühlmann, 2010)
performed selection without fixing regularization hyperparameters by choosing
parameters having a maximum selection frequency above a certain threshold.
Ruiz et al. (2020) considered estimating VAR models by pooling supports
meeting optimality criteria across resamples.

In this work we propose a sparse estimation method based on simple aggregation
operations applied to multiple estimates obtained from data resampling, and
demonstrate the superiority of the approach compared with simpler
alternatives. The method is developed in the context of first-order Poisson
GVAR processes, but generalizes easily to any support set estimation or model
selection problem. A high-level presentation of the methodology is given in
Section 2. We provide empirical support for the method in Section 3 along with
a discussion of special challenges in selecting parameter values for
simulating process realizations. We give a data application in 4 illustrating
the use of the methodology and its advantage over alternatives in practice.
Section 5 closes the paper with a summary of findings.

##  2 Methods

$\\{X_{t}\\}_{t\in\mathbb{Z}}$ is a generalized vector autoregressive (GVAR)
process if for some exponential family probability density $p$ and for every
$t\in\mathbb{Z}$:

$(X_{t,m}\;|\;X_{s}=x_{s},s<t)\stackrel{{\scriptstyle
indep.}}{{\sim}}p\left(\cdot\;;\theta_{t,m}=\theta_{m}(x_{t-1},x_{t-2},\dots)\right)\text{,
}m\in\\{1,\dots,M\\}$

Given an initial distribution $\nu_{0}$, the joint density of any finite-time
process realization $\\{x_{t}\\}_{t=0}^{T}$ is:

$\text{Pr}(X_{0}=x_{0},\dots,X_{T}=x_{T})=\nu_{0}(x_{0})\prod_{t=1}^{T}\prod_{m=1}^{M}p(x_{t,m};\theta_{t,m})$

The Poisson GVAR(D) process has the Poisson density as $p$ with a canonical
parameter $\theta_{t,m}$ linear in $D$ time lags:

$p(x_{t,m};\theta_{t,m})=\frac{\exp\\{x_{t,m}\theta_{t,m}-\exp(\theta_{t,m})\\}}{x_{t,m}!}\quad\text{and}\quad\theta_{t,m}=\nu_{m}+\sum_{d=1}^{D}a_{dm}^{\prime}x_{t-d}$

This process is identifiable for any values of the parameters
$\nu,A_{1},\dots,A_{D}$ provided $D<T$ but moments do not exist without
suitable constraints (Ruiz et al., 2022). The conditional means are log-linear
in past values:

$\ln\left\\{\mathbb{E}(X_{t}\;|\;X_{s},s<t)\right\\}=\nu+\sum_{d=1}^{D}A_{d}X_{t-d}$

The log-likelihood
$\ell(\nu,A_{1},\dots,A_{d};x)=\ln\\{\text{Pr}(X_{0}=x_{0},\dots,X_{T}=x_{T})\\}$
can be written directly from the process definition above. Maximum likelihood
estimates can be computed using standard techniques for generalized linear
models. The basic approach is to arrange the multivariate time series into a
single long pseudo-response vector, form a block-diagonal pseudo-covariate-
matrix containing the lagged values, and estimate the regression pseudo-
parameter:

$\beta=\text{vec}(\nu\;A_{1}\;\cdots\;A_{D})^{\prime}$

We express the estimand $S$ as the set of indices designating the nonzero
entries in $\beta$:

$S=\left\\{j:\beta_{j}\neq 0\right\\}$

The LASSO estimator of $S$ is:

$\hat{S}^{\lambda}=\\{j:\hat{\beta}_{j}^{\lambda}\neq
0\\}\quad\text{where}\quad\hat{\beta}_{j}^{\lambda}=\operatorname{argmin}_{\beta}\left\\{-\ell(\beta;x)+\lambda\|\text{vec}(A_{1}\;\cdots\;A_{D})\|_{1}\right\\}$

The details of the regression formulation and computations are given as
supplementary material for brevity. Here we simply note: (i) some care must be
taken when applying regularization constraints so as not to penalize the
intercept terms, which are distributed throughout the parameter vector unlike
in a standard regression problem; and (ii) although the parameters $A_{d}$
must be estimated simultaneously, the sparsity of the pseudo-covariate-matrix
should be leveraged for computational efficiency.

Our methods repeatedly calculate parameter estimates on data subsamples with
and without regularization. To simplify the presentation, we introduce the
shorthand notation in Table 1 for key operations: (i) resampling/subsampling;
(ii) LASSO estimation of support sets; (iii) maximum likelihood estimation of
parameters constrained to a support set; and (iv) calculation of an arbitrary
error metric or optimality criterion of the analyst’s choosing.

Table 1: Generic notations used throughout the paper. Notation |  Meaning | Example
---|---|---
$\text{resample}(\;\cdot\;)$ |  produce test and training subsamples/resamples | $\text{resample}(\text{data})=(\text{train},\text{test})$
$\hat{S}(\;\cdot\;;\lambda)$ |  compute LASSO supports for tuning hyperparameter $\lambda$ | $\hat{S}(\text{train};0.012)=\\{1,3\\}$
$\hat{\beta}(\;\cdot\;;S)$ |  compute MLE of parameter vector $\beta$ constrained to a specific support set $S$ | $\hat{\beta}(\text{train};\\{1,3\\})=(6.1\;00.76\;0\cdots\;0)^{T}$
$e(\;\cdot\;;\beta)$ |  compute a model-based error metric on some data for parameter value $\beta$ | $e(\text{test};(6.1\;0\;0.76\;0\;\cdots\;0)^{T})=153.8$

The standard approach to sparse estimation using the LASSO is to select the
regularization hyperparameter $\lambda$ by optimizing a criterion such as
goodness of fit or forecast accuracy. In other words, an analyst in practice
would optimize $e\\{\text{test}\;;\hat{\beta}^{\lambda}(\text{train})\\}$ by a
grid search over values of $\lambda$, where the data partitioning is employed
to reduce the bias of the error estimate. If sample size is limited, the data
may be repeatedly partitioned and average error across partitions chosen as
the quantity to optimize. However, we consider this to be a naive approach —
optimizing fit and/or forecasts for the LASSO estimator does not, in general,
produce a good estimator for the support set $S$, especially for vector
autoregressive models (Arnold et al., 2007; Davis et al., 2016; Shojaie and
Michailidis, 2010). Due to the dense parametrization — the number of
parameters $DM(M+1)$ scales quadratically in the process dimension $M$ and
linearly in the model order $D$ — often the model best approximating the true
data-generating process is extremely sparse. If so, strong regularization is
needed to accurately recover $S$, which induces heavy bias in the LASSO
estimator. In this circumstance, the hyperparameter $\lambda$ corresponding to
the best estimate of $S$ produces a poorly-fitting model that yields strongly
biased predictions. We therefore consider, to start, a modification of the
standard approach wherein regularization is used to determine candidate
supports but the selection of an optimal model is based on unpenalized
estimates. In short: optimize
$e[\text{test};\hat{\beta}\\{\text{train};\hat{S}(\text{train};\lambda)\\}]$
by a grid search instead of
$e\\{\text{test}\;;\hat{\beta}^{\lambda}(\text{train})\\}$. We consider this
as a benchmark method.

We propose two innovations. First we implement the ‘bootstrapped LASSO’
(BoLASSO) of Bach (2008), replacing the LASSO supports with the intersection
of LASSO support estimates calculated from a specified number of subsamples of
input data. In other words, we replace $\hat{S}(\cdot;\lambda)$ by:

$\hat{S}_{\text{agg}}(\;\cdot\;;\lambda)=\hat{S}(\text{resample}(\cdot);\lambda)\cap\cdots\cap\hat{S}(\text{resample}(\cdot);\lambda)$

However, we relax the strict intersection and instead introduce a threshold on
the proportion of sets in the collection that contain a given index. Next, we
propose selecting an optimal support set $\hat{S}_{j}^{*}$ on each data
partition $j$ and aggregating the collection by returning those indices that
occur in at least $100\times\gamma$% of optimal supports. That is, if
$\\{\hat{S}_{j}^{*}\\}$ are the optimal supports identified on each of the
data partitions, then we take as a final estimate of $S$:

$\hat{S}^{*}=\left\\{i:J^{-1}\sum_{j}\mathbf{1}\\{i\in\hat{S}^{*}_{j}\\}>\gamma\right\\}$

The model aggregation method is depicted in detail by the schematic in Fig. 1.

data $x$ $x^{*}_{11}$ train $x^{*}_{12}$ test
$\hat{S}\left(\lambda_{1};{\color[rgb]{1,0.3,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.3,0.3}x^{*}_{11}}\right)$$\hat{\beta}\left(\cdot\;;{\color[rgb]{1,0.3,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.3,0.3}x^{*}_{11}}\right)$$e\left(\cdot\;;{\color[rgb]{0.79,0.3,0.51}\definecolor[named]{pgfstrokecolor}{rgb}{0.79,0.3,0.51}x^{*}_{12}}\right)$$\hat{S}\left(\lambda_{2};{\color[rgb]{1,0.3,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.3,0.3}x^{*}_{11}}\right)$$\hat{\beta}\left(\cdot\;;{\color[rgb]{1,0.3,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.3,0.3}x^{*}_{11}}\right)$$e\left(\cdot\;;{\color[rgb]{0.79,0.3,0.51}\definecolor[named]{pgfstrokecolor}{rgb}{0.79,0.3,0.51}x^{*}_{12}}\right)$$\vdots$$\vdots$$\vdots$$\hat{S}\left(\lambda_{K};{\color[rgb]{1,0.3,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.3,0.3}x^{*}_{11}}\right)$$\hat{\beta}\left(\cdot\;;{\color[rgb]{1,0.3,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.3,0.3}x^{*}_{11}}\right)$$e\left(\cdot\;;{\color[rgb]{0.79,0.3,0.51}\definecolor[named]{pgfstrokecolor}{rgb}{0.79,0.3,0.51}x^{*}_{12}}\right)$$\hat{S}^{*}_{1}$min
$x^{*}_{21}$ train $x^{*}_{22}$ test
$\hat{S}\left(\lambda_{1};{\color[rgb]{1,0.65,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.65,0.3}x^{*}_{21}}\right)$$\hat{\beta}\left(\cdot\;;{\color[rgb]{1,0.65,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.65,0.3}x^{*}_{21}}\right)$$e\left(\cdot\;;{\color[rgb]{0.79,0.545,0.51}\definecolor[named]{pgfstrokecolor}{rgb}{0.79,0.545,0.51}x^{*}_{22}}\right)$$\hat{S}\left(\lambda_{2};{\color[rgb]{1,0.65,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.65,0.3}x^{*}_{21}}\right)$$\hat{\beta}\left(\cdot\;;{\color[rgb]{1,0.65,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.65,0.3}x^{*}_{21}}\right)$$e\left(\cdot\;;{\color[rgb]{0.79,0.545,0.51}\definecolor[named]{pgfstrokecolor}{rgb}{0.79,0.545,0.51}x^{*}_{22}}\right)$$\vdots$$\vdots$$\vdots$$\hat{S}\left(\lambda_{K};{\color[rgb]{1,0.65,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.65,0.3}x^{*}_{21}}\right)$$\hat{\beta}\left(\cdot\;;{\color[rgb]{1,0.65,0.3}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.65,0.3}x^{*}_{21}}\right)$$e\left(\cdot\;;{\color[rgb]{0.79,0.545,0.51}\definecolor[named]{pgfstrokecolor}{rgb}{0.79,0.545,0.51}x^{*}_{22}}\right)$$\hat{S}^{*}_{2}$min$\cdots$
$x^{*}_{J1}$ train $x^{*}_{J2}$ test
$\hat{S}\left(\lambda_{1};{\color[rgb]{0.82,0.82,0.4}\definecolor[named]{pgfstrokecolor}{rgb}{0.82,0.82,0.4}\pgfsys@color@cmyk@stroke{0}{0}{0.42}{0.18}\pgfsys@color@cmyk@fill{0}{0}{0.42}{0.18}x^{*}_{J1}}\right)$$\hat{\beta}\left(\cdot\;;{\color[rgb]{0.82,0.82,0.4}\definecolor[named]{pgfstrokecolor}{rgb}{0.82,0.82,0.4}\pgfsys@color@cmyk@stroke{0}{0}{0.42}{0.18}\pgfsys@color@cmyk@fill{0}{0}{0.42}{0.18}x^{*}_{J1}}\right)$$e\left(\cdot\;;{\color[rgb]{0.79,0.79,0.51}\definecolor[named]{pgfstrokecolor}{rgb}{0.79,0.79,0.51}\pgfsys@color@cmyk@stroke{0.21}{0.21}{0.49}{0}\pgfsys@color@cmyk@fill{0.21}{0.21}{0.49}{0}x^{*}_{J2}}\right)$$\hat{S}\left(\lambda_{2};{\color[rgb]{0.82,0.82,0.4}\definecolor[named]{pgfstrokecolor}{rgb}{0.82,0.82,0.4}\pgfsys@color@cmyk@stroke{0}{0}{0.42}{0.18}\pgfsys@color@cmyk@fill{0}{0}{0.42}{0.18}x^{*}_{J1}}\right)$$\hat{\beta}\left(\cdot\;;{\color[rgb]{0.82,0.82,0.4}\definecolor[named]{pgfstrokecolor}{rgb}{0.82,0.82,0.4}\pgfsys@color@cmyk@stroke{0}{0}{0.42}{0.18}\pgfsys@color@cmyk@fill{0}{0}{0.42}{0.18}x^{*}_{J1}}\right)$$e\left(\cdot\;;{\color[rgb]{0.79,0.79,0.51}\definecolor[named]{pgfstrokecolor}{rgb}{0.79,0.79,0.51}\pgfsys@color@cmyk@stroke{0.21}{0.21}{0.49}{0}\pgfsys@color@cmyk@fill{0.21}{0.21}{0.49}{0}x^{*}_{J2}}\right)$$\vdots$$\vdots$$\vdots$$\hat{S}\left(\lambda_{K};{\color[rgb]{0.82,0.82,0.4}\definecolor[named]{pgfstrokecolor}{rgb}{0.82,0.82,0.4}\pgfsys@color@cmyk@stroke{0}{0}{0.42}{0.18}\pgfsys@color@cmyk@fill{0}{0}{0.42}{0.18}x^{*}_{J1}}\right)$$\hat{\beta}\left(\cdot\;;{\color[rgb]{0.82,0.82,0.4}\definecolor[named]{pgfstrokecolor}{rgb}{0.82,0.82,0.4}\pgfsys@color@cmyk@stroke{0}{0}{0.42}{0.18}\pgfsys@color@cmyk@fill{0}{0}{0.42}{0.18}x^{*}_{J1}}\right)$$e\left(\cdot\;;{\color[rgb]{0.79,0.79,0.51}\definecolor[named]{pgfstrokecolor}{rgb}{0.79,0.79,0.51}\pgfsys@color@cmyk@stroke{0.21}{0.21}{0.49}{0}\pgfsys@color@cmyk@fill{0.21}{0.21}{0.49}{0}x^{*}_{J2}}\right)$$\hat{S}^{*}_{J}$min$\hat{S}^{*}=\left\\{i:\frac{1}{J}\sum_{j=1}^{J}\mathbf{1}\\{i\in\hat{S}^{*}_{j}\\}\geq\gamma\right\\}$
Figure 1: Schematic diagram of the model aggregation method. The functions
$\hat{S}(\;\cdot\;;\;\cdot\;)$, $\hat{\beta}(\;\cdot\;;\;\cdot\;)$, and
$e(\;\cdot\;;\;\cdot\;)$ are written with the order of arguments reversed to
better indicate that the arrows mean passing the output of one function to the
first argument of the next, e.g.,
$\hat{S}(\lambda_{1};x^{*}_{11})\rightarrow\hat{\beta}(\cdot\;;x^{*}_{11})$
means $\hat{\beta}(\hat{S}(\lambda_{1};x^{*}_{11});x^{*}_{11})$.

We refer to the first modification as ‘support aggregation’ and the second as
‘model aggregation’. Both are, essentially, thresholded intersection
operations – applied once to generate candidate models at the initial support
estimation step and again to combine the best candidates in place of selecting
a single model.

##  3 Simulations

We present the results of a simulation study designed to compare method
performance in the context of estimating the $M$-dimensional Poisson GVAR(1)
model:

$\displaystyle X_{t,m}|X_{t-1}\stackrel{{\scriptstyle
indep.}}{{\sim}}\text{Poisson}(\lambda_{t,m})\;,\quad m\in\\{1,\dots,M\\}$
$\displaystyle\log\lambda_{t}=\nu+AX_{t-1}$

The parameter of interest is $A\in\mathbb{R}^{M\times M}$ and it is assumed to
be $s$-sparse ($\|\text{vec}(A)\|_{0}=sM^{2}$). For the study, synthetic data
of length $T$ were generated and each method was used to estimate the
parameter $A$. The process dimension $M$, parametric sparsity $s$, and
realization length $T$ were varied systematically to generate the distinct
simulation settings. At each setting, several values of the parameter $A$ were
generated, and for each distinct value, several data realizations were
simulated.

Although the main objective is to compare our proposed method, the combination
of support aggregation and model aggregation, with the benchmark of model
selection via grid search, we also employed algorithms with each modification
separately in order to assess the marginal effects of each algorithm feature.
As such, four distinct estimation methods were applied to each dataset in the
simulation study.

### Study design

The study is a factorial design generated by three levels of process
dimensions, $M=10,15,20$, and three levels of sparsity, $s=0.01,0.02,0.05$.
For each factorial combination of sparsity and dimension, $10$ different $A$
matrices were generated. Further, for each distinct parameter matrix $A$ three
data realization lengths $T=500,1000,2000$ were considered; $5$ process
realizations of each length were simulated. In total, the number of datasets
simulated is $5\times 3\times 10\times 3\times 3=1350$ (dataset replicates
$\times$ realization lengths $\times$ parameter replicates $\times$ sparsity
levels $\times$ dimensions).

### Parameter generation

Poisson GVAR(1) process parameters cannot be drawn at random due to the
existence of regimes in the parameter space that are either (i) unstable
(unbounded means and/or variances) or (ii) unidentifiable in practice (produce
data realizations containing no information about the parameters with high
probability). Therefore some care is required in generating parameter values
for simulating data.

In our study, parameters are generated by first simulating a stable proposal
and then computing a ‘recoverability’ statistic to determine whether to accept
proposals. Proposals are constructed by the following steps:

  1. 1.

draw nonzero parameter magnitudes from a uniform distribution on $(0,1)$

  2. 2.

allocate a sign to each parameter

  3. 3.

fix the positions of nonzero parameters to satisfy stability constraints

To ensure that proposals are neither too difficult nor too easy to estimate,
we compute an expected ‘recoverability’ score based on the expression for the
signal-to-noise ratio in linear regression and similar to that used by Czanner
et al. (2015) to quantify signal-to-noise for GLMs. We characterize parameter
recoverability by the statistic:

$\text{recoverability}(A,\nu;x)=\frac{\text{dev}(0,\nu;x)-\text{dev}(A,\nu;x)}{\text{dev}(A,\nu;x)}$

This statistic is the relative change in deviance attributable to the
transition matrix $A$. (It can also be viewed as a generalization of the
$F$-test for overall significance in regression.) We approximate its expected
value by simulation: repeatedly simulate realizations using the parameters $A$
and $\nu$ and compute the average relative change in deviance across
realizations. Empirical experiments indicated a modest number ($n=10$
replicates) of long ($T=100,000$) simulated realizations produces a stable
estimate with relatively low variance. For our simulation study, intercepts
were fixed at a constant value for each setting and proposals for $A$ were
accepted if the estimated expected recoverability was at least 0.5 but did not
exceed 1.5 – that is, if the null model has an expected increase in deviance
of 50%-150% relative to the true data-generating model.

### Results

Figure 2: Average selection errors observed in simulation, reported as the
average proportion of matrix entries that are incorrectly classified as zero
or nonzero, for each combination of support estimation and support selection
method. Average selection errors are plotted against simulated realization
length. Separate panels are shown for each combination of sparsity $s$ (row)
and process dimension $M$ (column).

Figure 2 shows the average selection error for each method — the percentage of
erroneously-identified parameters in the estimated support set. Average error
rates are plotted against realization length $T$ for each combination of
parameter sparsity $s$ and process dimension $M$. The lowest-dimensional and
sparsest case ($M=10$ and $s=0.01$) represents selection and estimation of a
single nonzero autoregressive parameter out of a possible $100$; and the
highest-dimensional and densest case ($M=20$ and $s=0.05$) represents
selection and estimation of $20$ nonzero parameters out of a possible $400$.
For all methods, the average error rates range from $0$ to $10\%$ across
combinations of $M,s,T$, and selection errors decrease with $T$ and increase
with $M$ and $s$. Only in the lowest-dimensional and sparsest setting do all
methods perform equally well on average. Otherwise, support aggregation
outperforms the LASSO and model aggregation outperforms model selection; the
performance gaps widen as $s$ and $M$ increase. The combination performs
favorably, maintaining a $0-2\%$ error rate across cases. These improvements
are driven largely by better control of false positive rates, as shown in
Figure 3.

Figure 3: Average false positives (left) and false negatives (right) for each
method as a function of simulated series length. Separate panels are shown for
each combination of sparsity $s$ (row) and process dimension $M$ (column).

Figure 4 shows the marginal effect of each methodological innovation on a
dataset-by-dataset basis. At left, the error rate from the support aggregation
method is plotted against the error rate from the LASSO support estimation
method keeping other algorithm features fixed; these comparisons are shown for
every dataset generated in the simulation study and arranged in panels
according to the combination of sparsity $s$ and dimension $M$ in the same
layout as Figure 2. At right, the error rates for model aggregation are
compared with those for cross validation in the same manner. Simple regression
lines are superimposed to help visualize trends, along with a $y=x$ reference
that indicates equal performance. The trends, and in fact the majority of the
points themselves, are below the reference line, indicating that each
innovation outperforms its benchmark independently of other algorithm
features.

Figure 4: Marginal effect of algorithm features: support estimation method
(left) and selection vs. aggregation (right). Separate panels are shown for
each combination of sparsity $s$ (row) and process dimension $M$ (column).

##  4 Application

We demonstrate application of our methodology using paleoclimatology data from
Barron et al. (2005b); data are publicly available online through the World
Data Service for Paleoclimatology managed by the National Centers for
Environmental Information at NOAA Barron et al. (2005a). The data are
multivariate counts of individuals in eight diatom taxa over time, shown in
Fig. 5. Data were collected from samples taken at evenly-spaced depths in
marine sediment cores from the Guayamas basin in Baja California; the depths
sampled yield approximately homogeneous time intervals of 40-60 years between
observations, and in all, the data span 16,000 years. The observation window
covers the transition from the Pleistocene epoch to the Holocene epoch around
11,700 before present. Sea surface temperature reconstructions from the
sampling location reflect a shift toward warmer average temperatures occurring
at approximately that time.

Figure 5: Counts of the numbers of individuals over time identified at the
taxon level for eight diatom taxa present in marine sediment core samples from
the Guayamas basin in Baja California; sea surface temperature reconstructions
show a dramatic change following the end of the Pleistocene epoch, indicated
by the dashed red line.

We used our method to estimate a interaction networks describing serial
interdependence between diatom taxa. By partitioning the data according to
geological epoch and fitting separate models to each partition, we identified
a clear and apparent change in the connectivity structure before and after the
transition between geological epochs. This change was not identifiable using
alternative estimation methods. We used log-linear Poisson GVAR(1) models to
represent the dependence structure for each epoch and fit the two models using
three methods each: a naive direct application of the LASSO; the benchmark
method from our simulation study; and the aggregation method. Results are
shown in Fig. 6. Due to its sparsity, only the aggregation method clearly
supports the conclusion that there is a change in connectivity structure.
Although there are some differences in connectivity as estimated by the
alternative methods, there also remain inter-taxon links that persist across
epochs that introduce ambiguity in interpreting the results; there is no such
ambiguity using our method.

Figure 6: Interaction graphs representing inter-taxon dependencies before and
after the climate change event as estimated by three methods: a ‘naive’ method
that consists in directly applying the LASSO at left; the benchmark considered
in our simulation study in center; and our proposed method at right. Only the
proposed method identifies a clear and apparent change in connectivity.

##  5 Discussion

This paper presents novel methodology for estimation of parameter supports
using resampling and aggregation operations. The methodology is fully general
but developed in the context of GVAR models for multivariate count-valued time
series. We provide empirical validation demonstrating a consistent improvement
over benchmark methods across process and parameter dimension, sparsity level,
and series length or sample size; not only is there a consistent average
effect, but the method outperforms alternatives for the vast majority of
simulated datasets. The simulation study reveals that most of the advantage of
our method lies in better control of false positive rates in estimating
parameter supports. We illustrate the practical advantages of this behavior
with a data analysis involving estimating ecological interaction networks and
comparing the connectivity structures for differences. While there is no
ground truth in this context, the greater sparsity of estimates resulting from
our method makes apparent differences in network structure that are otherwise
obscured by additional parameters that reflect either errors or weak
correlations.

It is worth highlighting that all of the benchmark methods we consider in our
simulation study involve using maximum likelihood refits of LASSO estimates
rather than the LASSO estimates themselves. Although we exclude ‘raw’ LASSO
estimates from consideration in the simulations, we illustrate the effect of
refitting in our data analysis by including the ‘naive’ method of direct LASSO
estimation; the naive method produces a nearly-fully-dense model. Anecdotally,
similar behavior was observed with simulated data in the course of developing
the simulation study — the naive LASSO estimate that optimized deviance
produced an unacceptably high selection error rate. Thus, the benchmarks in
this paper would appear to represent a considerable improvement over the
methods of Mark et al. (2017); Hall et al. (2018); Pandit et al. (2020) in the
small-sample and low-dimensional setting.

The resampling and aggregation approach overcomes the problem of choosing
exactly one regularization hyperparameter. However, the method still requires
a good choice of regularization path at an appropriate resolution, and it
introduces other hyperparameters that require some amount of tuning: the
numbers of resamples and the threshold parameter for the aggregation
operations. In practice, we found that these are not too challenging to
specify, but do require some trial and error. We found that relatively low
numbers of resamples — around 10 — were sufficient for stable performance
across settings, and we cross-validated the entire estimation procedure(s) to
select an optimal threshold parameter for model aggregation.

Although it is under-emphasized throughout the paper, the resampling of
multiple time series can be accomplished in a number of ways. Ideally, if
enough data are available, non-overlapping subsamples could be used to
minimize the correlation between estimates on each subsample due to long-range
dependence in the data. In our implementations, we subdivided data evenly into
blocks and omitted one block at a time, concatenating the remaining data and
introducing one discontinuity point across which the subsampled process is no
longer stationary; the concatenated data was used for training and the held-
out block was used for testing. Thus, for our method, increasing the number of
resamples taken diminishes the length of series used for validation. One might
instead use a block bootstrap or a residual bootstrap method such as described
by Kreiss and Lahiri (2012).

It is an important feature of our method that support sets are re-estimated
across the regularization paths for each resample, rather than fixed based on
a single estimate from the full dataset or some other heuristic. This is a key
difference between our approach and the method by which Ruiz et al. (2020)
generate and select model candidates for pooling. We found that conditioning
our procedure on a fixed collection of candidate supports, although
dramatically less computationally burdensome, produces a noteworthy
deterioration of selection performance.

Finally, we note that the simulation of stable process realizations without
artificial thresholds proved to be challenging, as there is no obvious
necessary and sufficient condition on the parameters that ensures moments
exist. Furthermore, some process parametrizations are extremely challenging to
recover from data because the likelihood varies little across large
neighborhoods with high probability. In the course of developing empirical
validation of our methodology, we therefore also developed heuristics for
selecting process parameters that are both stable and recoverable.

## References

  * Arnold et al. [2007] Andrew Arnold, Yan Liu, and Naoki Abe.  Temporal causal modeling with graphical granger methods.  In _Proceedings of the 13th ACM SIGKDD international conference on Knowledge discovery and data mining_ , pages 66–75, 2007.
  * Bach [2008] Francis R Bach.  Bolasso: model consistent lasso estimation through the bootstrap.  In _Proceedings of the 25th international conference on Machine learning_ , pages 33–40, 2008.
  * Bach and Jordan [2004] Francis R Bach and Michael I Jordan.  Learning graphical models for stationary time series.  _IEEE transactions on signal processing_ , 52(8):2189–2199, 2004.
  * Barron et al. [2005a] J.A. Barron, D. Bukry, and J.L. Bischoff.  High resolution guaymas basin geochemical, diatom, and silicoflagellate data.  IGBP PAGES/World Data Center for Paleoclimatology, Data Contribution Series # 2005-022, 2005a.  NOAA/NGDC Paleoclimatology Program, Boulder CO, USA.
  * Barron et al. [2005b] John A Barron, David Bukry, and Walter E Dean.  Paleoceanographic history of the guaymas basin, gulf of california, during the past 15,000 years based on diatoms, silicoflagellates, and biogenic sediments.  _Marine Micropaleontology_ , 56(3-4):81–102, 2005b.
  * Basu et al. [2015] Sumanta Basu, George Michailidis, et al.  Regularized estimation in sparse high-dimensional time series models.  _The Annals of Statistics_ , 43(4):1535–1567, 2015.
  * Brown et al. [2004] Emery N Brown, Robert E Kass, and Partha P Mitra.  Multiple neural spike train data analysis: state-of-the-art and future challenges.  _Nature neuroscience_ , 7(5):456–461, 2004.
  * Cox et al. [1981] David R Cox, Gudmundur Gudmundsson, Georg Lindgren, Lennart Bondesson, Erik Harsaae, Petter Laake, Katarina Juselius, and Steffen L Lauritzen.  Statistical analysis of time series: Some recent developments.  _Scandinavian Journal of Statistics_ , pages 93–115, 1981.
  * Czanner et al. [2015] Gabriela Czanner, Sridevi V Sarma, Demba Ba, Uri T Eden, Wei Wu, Emad Eskandar, Hubert H Lim, Simona Temereanca, Wendy A Suzuki, and Emery N Brown.  Measuring the signal-to-noise ratio of a neuron.  _Proceedings of the National Academy of Sciences_ , 112(23):7141–7146, 2015.
  * Dahlhaus and Eichler [2003] Rainer Dahlhaus and Michael Eichler.  Causality and graphical models in time series analysis.  _Oxford Statistical Science Series_ , pages 115–137, 2003.
  * Davis et al. [2016] Richard A Davis, Pengfei Zang, and Tian Zheng.  Sparse vector autoregressive modeling.  _Journal of Computational and Graphical Statistics_ , 25(4):1077–1096, 2016.
  * Fahrmeir and Tutz [1994] Ludwig Fahrmeir and Gerhard Tutz.  _Multivariate Statistical Modeling Based on Generalized Linear Models_.  Springer-Verlag: New York, 1994.
  * Fan et al. [2011] Jianqing Fan, Jinchi Lv, and Lei Qi.  Sparse high-dimensional models in economics.  _Annu. Rev. Econ_ , 3:291–317, 2011.
  * Fokianos [2021] Konstantinos Fokianos.  Multivariate count time series modelling.  _Econometrics and Statistics_ , 2021.
  * Friedman et al. [2007] Jerome Friedman, Trevor Hastie, Holger Höfling, and Robert Tibshirani.  Pathwise coordinate optimization.  _The Annals of Applied Statistics_ , 1(2):302–332, 2007.
  * Friedman et al. [2010] Jerome Friedman, Trevor Hastie, and Rob Tibshirani.  Regularization paths for generalized linear models via coordinate descent.  _Journal of statistical software_ , 33(1):1, 2010\.
  * Granger [1969] Clive WJ Granger.  Investigating causal relations by econometric models and cross-spectral methods.  _Econometrica: journal of the Econometric Society_ , pages 424–438, 1969.
  * Hall et al. [2018] Eric C Hall, Garvesh Raskutti, and Rebecca M Willett.  Learning high-dimensional generalized linear autoregressive models.  _IEEE Transactions on Information Theory_ , 65(4):2401–2422, 2018.
  * Han et al. [2015] Fang Han, Huanran Lu, and Han Liu.  A direct estimation of high dimensional stationary vector autoregressions.  _Journal of Machine Learning Research_ , 2015.
  * Han et al. [2023] Yuefeng Han, Ruey S Tsay, and Wei Biao Wu.  High dimensional generalized linear models for temporal dependent data.  _Bernoulli_ , 29(1):105–131, 2023.
  * Hsu et al. [2008] Nan-Jung Hsu, Hung-Lin Hung, and Ya-Mei Chang.  Subset selection for vector autoregressive processes using lasso.  _Computational Statistics & Data Analysis_, 52(7):3645–3657, 2008.
  * Karlis [2016] Dimitris Karlis.  Models for multivariate count time series.  _Handbook of discrete-valued time series_ , pages 407–424, 2016.
  * Kreiss and Lahiri [2012] Jens-Peter Kreiss and Soumendra Nath Lahiri.  Bootstrap methods for time series.  In _Handbook of statistics_ , volume 30, pages 3–26. Elsevier, 2012\.
  * Lozano et al. [2009] Aurélie C Lozano, Naoki Abe, Yan Liu, and Saharon Rosset.  Grouped graphical granger modeling for gene expression regulatory networks discovery.  _Bioinformatics_ , 25(12):i110–i118, 2009.
  * Mark et al. [2017] Ben Mark, Garvesh Raskutti, and Rebecca Willett.  Network estimation via poisson autoregressive models.  In _2017 IEEE 7th International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP)_ , pages 1–5, 2017.
  * Meinshausen and Bühlmann [2010] Nicolai Meinshausen and Peter Bühlmann.  Stability selection.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 72(4):417–473, 2010.
  * Pandit et al. [2020] Parthe Pandit, Mojtaba Sahraee-Ardakan, Arash A Amini, Sundeep Rangan, and Alyson K Fletcher.  Generalized autoregressive linear models for discrete high-dimensional data.  _IEEE Journal on Selected Areas in Information Theory_ , 1(3):884–896, 2020.
  * Paul et al. [2008] Michaela Paul, Leonhard Held, and André M Toschke.  Multivariate modelling of infectious disease surveillance data.  _Statistics in medicine_ , 27(29):6250–6267, 2008\.
  * Pillow et al. [2008] Jonathan W Pillow, Jonathon Shlens, Liam Paninski, Alexander Sher, Alan M Litke, EJ Chichilnisky, and Eero P Simoncelli.  Spatio-temporal correlations and visual signalling in a complete neuronal population.  _Nature_ , 454(7207):995–999, 2008.
  * Ruiz et al. [2020] Trevor Ruiz, Sharmodeep Bhattacharyya, Mahesh Balasubramanian, and Kristofer Bouchard.  Sparse and low-bias estimation of high dimensional vector autoregressive models.  In _Learning for Dynamics and Control_ , pages 55–64. PMLR, 2020\.
  * Ruiz et al. [2022] Trevor Ruiz, Sharmodeep Bhattacharyya, and Sarah C Emerson.  A graphical sufficient condition for the stability of first-order log-linear poisson generalized vector autoregressive processes.  Preprint available at SSRN: https://ssrn.com/abstract=4283477 or http://dx.doi.org/10.2139/ssrn.4283477, 2022.
  * Shojaie and Michailidis [2010] Ali Shojaie and George Michailidis.  Discovering graphical granger causality using the truncating lasso penalty.  _Bioinformatics_ , 26(18):i517–i523, 2010.
  * Song and Bickel [2011] Song Song and Peter J Bickel.  Large vector auto regressions.  _arXiv preprint arXiv:1106.3915_ , 2011.
  * Tibshirani [1996] Robert Tibshirani.  Regression shrinkage and selection via the lasso.  _Journal of the Royal Statistical Society. Series B (Methodological)_ , pages 267–288, 1996.
  * Valdés-Sosa et al. [2005] Pedro A Valdés-Sosa, Jose M Sánchez-Bornot, Agustín Lage-Castellanos, Mayrim Vega-Hernández, Jorge Bosch-Bayard, Lester Melie-García, and Erick Canales-Rodríguez.  Estimating brain functional connectivity with sparse multivariate autoregression.  _Philosophical Transactions of the Royal Society of London B: Biological Sciences_ , 360(1457):969–981, 2005.
  * van der Kooij [2007] Anita J van der Kooij.  _Prediction accuracy and stability of regression with optimal scaling transformations_.  PhD thesis, Lieden University, 2007.
  * Zou and Hastie [2005] Hui Zou and Trevor Hastie.  Regularization and variable selection via the elastic net.  _Journal of the royal statistical society: series B (statistical methodology)_ , 67(2):301–320, 2005.

## Appendix

##  Appendix A Algorithms

The algorithms used to produce estimates for numerical experiments and data
analysis in the main body of the paper are presented in detail in this
appendix. Essentially, we apply the coordinate descent method of Friedman et
al. [2010] after manipulating the likelihood optimization problem into a form
similar to univariate regression and partitioning coordinates into blocks for
computational efficiency. The core technique is based on coordinate descent
updates for the weighted least squares regression estimate with an elastic net
penalty. We first present this method, with full derivations, then the
approach of Friedman et al. [2010], and lastly the details of our
implementation for GVAR(D) estimation.

Weighted least squares with elastic net. Consider estimation of the regression
model

$z=X\beta+\epsilon$

where $z$ is the response, $X$ is the matrix of fixed covariate information,
$\beta$ is a vector of fixed but unknown parameters, and $\epsilon$ is a
vector of uncorrelated random errors with mean zero.

The weighted least squares estimate of $\beta$ is given by

$\hat{\beta}=\operatorname{argmin}_{\beta}\left\\{\frac{1}{2}(z-X\beta)^{\prime}W(z-X\beta)\right\\}$

where $W$ is a diagonal matrix of known weights. Denote the objective function
by

$f(\beta)=\frac{1}{2}(z-X\beta)^{\prime}W(z-X\beta)$ (1)

The weighted least squares estimate with an elastic net penalty is obtained by
adding a penalty term $P(\beta)$ to the objective function

$\hat{\beta}=\operatorname{argmin}_{\beta}\left\\{\frac{1}{2}(z-X\beta)^{\prime}W(z-X\beta)+\lambda
P(\beta)\right\\}$ (2)

where $P(\beta)=\frac{1-\alpha}{2}\|\beta\|_{2}^{2}+\alpha\|\beta\|_{1}$ and
$\lambda$ is a hyperparameter that controls the strength of the penalty. The
hyperparameter $\alpha$ controls the relative balance of the $L_{2}$ (ridge)
and $L_{1}$ (LASSO) terms in the penalty.

The coordinate descent approach to computing $\hat{\beta}$ in Equation (2) is
to iteratively minimize the objective function in each $\beta_{j}$ conditional
on the current values of the other parameters $\beta_{k}$ for $k\neq j$ until
convergence. Each iteration can be accomplished by a simple update rule that
can be equivalently derived by a number of approaches [Zou and Hastie, 2005,
van der Kooij, 2007, Friedman et al., 2007, 2010]. Here the rule is derived
from the update given by Newton’s method.

For the purpose of taking derivatives in $\beta_{j}$, the first portion of the
objective function $f(\beta)$ can be rewritten as

$f(\beta)=\frac{1}{2}\left(z-\sum_{k\neq
j}x_{k}\beta_{k}-x_{j}\beta_{j}\right)^{\prime}W\left(z-\sum_{k\neq
j}x_{k}\beta_{k}-x_{j}\beta_{j}\right)$

where $x_{j}$ denotes the $j$th column of $X$. The first and second partial
derivatives of this portion with respect to $\beta_{j}$ are

$\displaystyle\frac{\partial}{\partial\beta_{j}}f(\beta)$
$\displaystyle=-\left(z-\sum_{k\neq
j}x_{k}\beta_{k}-x_{j}\beta_{j}\right)^{\prime}Wx_{j}$
$\displaystyle=\beta_{j}x_{j}^{\prime}Wx_{j}-\left(z-\sum_{k\neq
j}x_{k}\beta_{k}\right)^{\prime}Wx_{j}$
$\displaystyle=\beta_{j}x_{j}^{\prime}Wx_{j}-\left(z-z^{(j)}\right)^{\prime}Wx_{j}$
$\displaystyle\frac{\partial^{2}}{\partial^{2}\beta_{j}}f(\beta)$
$\displaystyle=x_{j}^{\prime}Wx_{j}$

where $z^{(j)}$ denotes $\sum_{k\neq j}x_{k}\beta_{k}$. The quantity
$\left(z-z^{(j)}\right)$ is known as the $j$th partial residual.

Now, the derivatives of the penalty portion of the objective function are

$\displaystyle\frac{\partial}{\partial\beta_{j}}\left(\lambda\frac{1-\alpha}{2}\|\beta\|_{2}^{2}+\lambda\alpha\|\beta\|_{1}\right)$
$\displaystyle=\lambda(1-\alpha)\beta_{j}+\lambda\alpha\frac{\partial}{\partial\beta_{j}}|\beta_{j}|$
$\displaystyle=\begin{cases}\lambda(1-\alpha)\beta_{j}+\lambda\alpha&\text{ ,
if }\beta_{j}>0\\\ \lambda(1-\alpha)\beta_{j}-\lambda\alpha&\text{ , if
}\beta_{j}<0\\\
\lambda(1-\alpha)\beta_{j}+\lambda\alpha\partial|\beta_{j}|&\text{ , if
}\beta_{j}=0\end{cases}$
$\displaystyle\frac{\partial^{2}}{\partial^{2}\beta_{j}}\left(\lambda\frac{1-\alpha}{2}\|\beta\|_{2}^{2}+\lambda\alpha\|\beta\|_{1}\right)$
$\displaystyle=\lambda(1-\alpha)$

where $\partial|\beta_{j}|$ denotes a subgradient. Supposing the current
estimate is $\hat{\beta}^{(k)}$, the Newton update for minimization in the
$j$th coordinate is

$\hat{\beta}^{(k+1)}_{j}\longleftarrow\hat{\beta}^{(k)}_{j}-\frac{\frac{\partial}{\partial\beta_{j}}\left(f(\beta)+\lambda
P(\beta)\right)\rvert_{\beta=\hat{\beta}^{(k)}}}{\frac{\partial^{2}}{\partial^{2}\beta_{j}}\left(f(\beta)+\lambda
P(\beta)\right)\rvert_{\beta=\hat{\beta}^{(k)}}}$ (3)

When articulated casewise according to whether $\hat{\beta}^{(k)}_{j}>0$ or
$\hat{\beta}^{(k)}_{j}<0$ or $\hat{\beta}^{(k)}_{j}=0$, Equation (3) gives the
update rule (after minor algebraic simplification)

$\hat{\beta}_{j}^{(k+1)}\longleftarrow\begin{cases}\frac{\left(z-z^{(j)}\right)^{\prime}Wx_{j}-\lambda\alpha}{x_{j}^{\prime}Wx_{j}+\lambda(1-\alpha)}&\text{
, if }\hat{\beta}_{j}^{(k)}>0\\\
\frac{\left(z-z^{(j)}\right)^{\prime}Wx_{j}+\lambda\alpha}{x_{j}^{\prime}Wx_{j}+\lambda(1-\alpha)}&\text{
, if }\hat{\beta}_{j}^{(k)}<0\\\
\frac{\left(z-z^{(j)}\right)^{\prime}Wx_{j}-\lambda\alpha
c}{x_{j}^{\prime}Wx_{j}+\lambda(1-\alpha)}&\text{ , if
}\hat{\beta}_{j}^{(k)}=0\\\ \end{cases}$

where $c\in[-1,1]$. Now if
$|\left(z-z^{(j)}\right)^{\prime}Wx_{j}|\leq\lambda\alpha$, the (sub)gradient
in Equation (3) is not a descent direction, and in this case, the optimal
value of the objective function is achieved by setting
$\hat{\beta}_{j}^{(k+1)}$ to zero. Rewriting the conditions in terms of
$\left(z-z^{(j)}\right)^{\prime}Wx_{j}$ and dropping the iteration count $k$
yields the update rule

$\hat{\beta}_{j}\longleftarrow\begin{cases}\frac{\left(z-z^{(j)}\right)^{\prime}Wx_{j}-\lambda\alpha}{x_{j}^{\prime}Wx_{j}+\lambda(1-\alpha)}&\text{
, if }\left(z-z^{(j)}\right)^{\prime}Wx_{j}>0\text{ and
}|\left(z-z^{(j)}\right)^{\prime}Wx_{j}|>\lambda\alpha\\\
\frac{\left(z-z^{(j)}\right)^{\prime}Wx_{j}+\lambda\alpha}{x_{j}^{\prime}Wx_{j}+\lambda(1-\alpha)}&\text{
, if }\left(z-z^{(j)}\right)^{\prime}Wx_{j}<0\text{ and
}|\left(z-z^{(j)}\right)^{\prime}Wx_{j}|>\lambda\alpha\\\ 0&\text{ , if
}\text{ and
}|\left(z-z^{(j)}\right)^{\prime}Wx_{j}|\leq\lambda\alpha\end{cases}$

Finally, the update can be re-expressed in terms of the soft-threshold
function

$S(x,y)=\text{sign}(x)(|x|-y)_{+}=\begin{cases}x-y&\text{ , if }x>0\text{ and
}|x|>y\\\ x+y&\text{ , if }x<0\text{ and }|x|>y\\\ 0&\text{ , if }|x|\leq
y\end{cases}$

compactly as the coordinate update

$\hat{\beta}_{j}\longleftarrow\frac{S\left(\left(z-z^{(j)}\right)^{\prime}Wx_{j},\lambda\alpha\right)}{x_{j}^{\prime}Wx_{j}+\lambda(1-\alpha)}$

The coordinate descent algorithm using this update rule and a gradient norm
stopping criterion is given as Algorithm A.1.

Algorithm A.1 Coordinate descent algorithm for computing the weighted least
squares estimate with an elastic net penalty.

0:

Data $z,X$; hyperparameters $\lambda,\alpha$; and convergence threshold
$\delta$

Initialize $\hat{\beta}$

while $\|\nabla f(\hat{\beta})\|_{2}>\delta$ do

for $j=1,\dots,J$ do

$z^{(j)}\longleftarrow\sum_{k\neq j}x_{k}\hat{\beta}_{k}$

$\hat{\beta}_{j}\longleftarrow\frac{S\left(\left(z-z^{(j)}\right)^{\prime}Wx_{j},\lambda\alpha\right)}{x_{j}^{\prime}Wx_{j}+\lambda(1-\alpha)}$

end for

end while

$\hat{\beta}$

Penalized IRLS for univariate generalized linear models. Consider the
univariate generalized linear model: a response vector $y$ and covariate
information $X$ where

$y_{i}\stackrel{{\scriptstyle
indep.}}{{\sim}}p(\cdot\;;\theta_{i})\quad\text{and}\quad
g(\mathbb{E}y)=g(\mu)=\eta=X\beta$

where $p$ is an exponential family density of the form

$p(y_{i};\theta_{i})\propto\exp\left\\{\frac{y_{i}\theta_{i}-b(\theta_{i})}{a(\phi_{i})}\right\\}$

Assume that $g$ is the canonical link function that results from equating
$\eta=\theta$. The log-likelihood as a function of the linear predictor $\eta$
is, up to a constant,

$\ell(\eta;x,y)=\sum_{i=1}^{n}\left(\frac{y_{i}\eta_{i}-b(\eta_{i})}{a(\phi_{i})}\right)$

The classical technique for unpenalized maximum likelihood estimation
maximizes $\ell(\eta;x,y)$ iteratively by:

  1. 1.

computing an approximation $\tilde{\ell}$ to the log-likelihood from current
estimates;

  2. 2.

maximizing the approximation $\tilde{\ell}$ to get new estimates;

  3. 3.

repeating 1 - 2 until the estimates converge.

Conveniently, the approximation $-\tilde{\ell}$ takes the form of the loss
function for a weighted least squares problem for which a closed form solution
is available, which eases calculations in step 2.

Consider now computing the penalized MLE

$\hat{\beta}=\operatorname{argmin}_{\beta}\left\\{-\ell(\eta;x,y)+\lambda
P(\beta)\right\\}$ (4)

The same strategy for unpenalized estimation is applicable with the penalty,
though generally iterative methods will replace closed form solutions for the
optimization step (step 2). The objective function $-\ell(\eta;x,y)+\lambda
P(\beta)$ can be approximated by the loss function for a penalized weighted
least squares problem. For the elastic net penalty, the approximation can be
maximized with respect to $\beta$ via Algorithm A.1. The penalized MLE in
Equation (4) with the elastic net penalty can therefore be found iteratively
by:

  1. 1.

computing the approximation $-\tilde{\ell}(\eta;x,y)+\lambda P(\beta)$ to the
(negative) penalized log-likelihood from current estimates;

  2. 2.

maximizing the approximation $-\tilde{\ell}(\eta;x,y)+\lambda P(\beta)$ with
respect to one coordinate via the update rule in Algorithm A.1 to get new
estimates;

  3. 3.

repeating 1 - 2 for each coordinate and cycling through the coordinates until
the estimates converge.

This procedure is given as Algorithm A.2.

Algorithm A.2 Coordinatewise IRLS algorithm for estimation of generalized
linear models with an elastic net penalty.

0: data $x,y$; hyperparameters $\lambda,\alpha$; convergence threshhold
$\delta$

Initialize $\hat{\beta}$

while $\|\nabla\ell(\hat{\eta};x,y)\|_{2}>\delta$ do

for $j=1,\dots,J$ do

$W\longleftarrow\left(\text{diag}(g^{\prime}(\hat{\mu})\right)^{-1}$

$z\longleftarrow x\hat{\beta}-(y-\hat{\mu})W^{-1}$

$z^{(j)}\longleftarrow\sum_{k\neq j}x_{k}\hat{\beta}_{k}$

$\hat{\beta}_{j}\longleftarrow\frac{S\left(\left(z-z^{(j)}\right)^{\prime}Wx_{j},\lambda\alpha\right)}{x_{j}^{\prime}Wx_{j}+\lambda(1-\alpha)}$

end for

end while

$\hat{\beta}$

The update rules for $z$ and $W$ in Algorithm A.2 are derived from a Taylor
expansion around current estimates. Let $\hat{\eta}$ denote the linear
predictor from the current estimates $\hat{\beta}$, $\ell(\eta)$ denote the
log-likelihood (suppressing the data arguments), $v$ and $W$ denote the
gradient and Hessian of $\ell(\eta)$. A second-order Taylor expansion of the
log-likelihood about $\hat{\eta}$ gives

$\displaystyle\ell(\eta)$
$\displaystyle\approx\ell(\hat{\eta})+(\eta-\hat{\eta})^{\prime}\ell^{\prime}(\hat{\eta})-\frac{1}{2}(\eta-\hat{\eta})^{\prime}\ell^{\prime\prime}(\hat{\eta})(\eta-\hat{\eta})$
$\displaystyle=C+(\eta-\hat{\eta})^{\prime}\ell^{\prime}(\hat{\eta})-\frac{1}{2}(\eta-\hat{\eta})^{\prime}\ell^{\prime\prime}(\hat{\eta})(\eta-\hat{\eta})$
$\displaystyle=C+(\eta-\hat{\eta})^{\prime}v-\frac{1}{2}(\eta-\hat{\eta})^{\prime}W(\eta-\hat{\eta})$
$\displaystyle=C^{*}-\frac{1}{2}(\eta-\hat{\eta}-vW^{-1})^{\prime}W(\eta-\hat{\eta}-vW^{-1})$
$\displaystyle=C^{*}-\frac{1}{2}(\eta-z)^{\prime}W(\eta-z)$

with $z=\hat{\eta}-vW^{-1}$. Now $W$ and $z$ depend on the current estimates
and can be derived in general for any GLM. Since by hypothesis
$\mu_{i}=b^{\prime}(\theta_{i})=g^{-1}(\eta_{i})$,

$\displaystyle\frac{\partial\ell}{\partial\eta_{i}}$
$\displaystyle=\left(\frac{y_{i}-b^{\prime}(\theta_{i})}{a_{i}(\phi)}\right)\frac{\partial\theta_{i}}{\partial\eta_{i}}$
$\displaystyle=\left(\frac{y_{i}-b^{\prime}(\theta_{i})}{a_{i}(\phi)}\right)\frac{\partial\theta_{i}}{\partial\mu_{i}}\frac{\partial\mu_{i}}{\partial\eta_{i}}$
$\displaystyle=\left(\frac{y_{i}-b^{\prime}(\theta_{i})}{a_{i}(\phi)}\right)\frac{1}{g^{\prime}(\mu_{i})}\frac{1}{b^{\prime\prime}(\theta_{i})}$
$\displaystyle=(y_{i}-\mu_{i})\left(a_{i}(\phi)g^{\prime}(\mu_{i})b^{\prime\prime}(\theta_{i})\right)^{-1}$

When $a_{i}(\phi)=1$ and $g$ is the canonical link,
$g^{\prime}(\mu_{i})=\frac{\partial\eta_{i}}{\partial\mu_{i}}$ and
$b^{\prime\prime}(\theta_{i})=\frac{\partial\mu_{i}}{\partial\theta_{i}}$, so

$\frac{\partial\ell}{\partial\eta_{i}}=y_{i}-\mu_{i}$

and

$\frac{\partial^{2}\ell}{\partial\eta_{i}\partial\eta_{j}}=\begin{cases}-\frac{1}{g^{\prime}(\mu_{i})}&\qquad
i=j\\\ 0&\qquad i\neq j\end{cases}$

So $W=\text{diag}(1/g^{\prime}(\mu))$ and $v=(y_{1}-\mu_{1}\cdots
y_{N}-\mu_{N})$. Therefore, given current estimates, updating the
approximation amounts to computing

$\displaystyle z$ $\displaystyle=X^{\prime}\hat{\beta}-(y-\hat{\mu})W^{-1}$
(5) $\displaystyle W$ $\displaystyle=\left(g^{\prime}(\hat{\mu})\right)^{-1}$
(6)

Equations (5) and (6) give the updates that appear in the outer while loop of
Algorithm A.2.

GVAR(D) estimation. Algorithm A.2 can be modified for efficient estimation of
GVAR(D) models with an elastic net penalty using a blockwise update strategy
for the inner while loop.

The GVAR(D) model for data $\\{x_{t}\in\mathbb{R}^{M}\\}_{t=0}^{T}$ is the
random process $\\{X_{t}\in\mathbb{R}^{M}\\}_{t=0}^{T}$ characterized by the
conditional mean structure

$g\left(\mathbbm{E}\left(X_{t}\;|\;X_{t-1}=x_{t-1},\dots,X_{t-D}=x_{t-D}\right)\right)=\nu+\sum_{d=1}^{D}A_{d}x_{t-d}$

under the assumption that for each $m=1,\dots,M$ and each $t$

$X_{t,m}\;|\;X_{t-1}=x_{t-1},\dots,X_{t-D}=x_{t-D}\stackrel{{\scriptstyle
indep.}}{{\sim}}p\left(\;\cdot\;\Big{|}\;\theta_{t,m}=\nu_{m}+\sum_{d=1}^{D}a_{d,m}^{\prime}x_{t-d}\right)$

where $p(\cdot|\theta)$ is an exponential family density and $g$ is the
canonical link so that
$\eta_{t,m}=\theta_{t,m}=g(\mu_{t,m})=\mathbbm{E}(X_{t,m}|X_{t-1},\dots,X_{t-D})$).
Denoting the data by $X=(X_{0}\;\cdots\;X_{T})$ log-likelihood by $\ell(B;X)$,
the penalized maximum likelihood estimate of $B$ is

$\hat{B}=\arg\min_{B}\left\\{-\ell(B;X)+\lambda P(B)\right\\}$

where the explicit form of the likelihood depends only on the underlying
density $p$. Here let $P$ be the elastic net penalty applied only to the
elements of $A_{1},\dots,A_{D}$

$P(B)=\frac{1-\alpha}{2}\sum_{d,i,j}a_{dij}^{2}+\alpha\sum_{d,i,j}|a_{dij}|$

An iterative estimation algorithm based on Algorithm A.2 is derived below by
vectorizing the problem and developing a blockwise update strategy.

The GVAR(D) model can be expressed as a multivariate generalized regression
model

$g\left(\mathbb{E}(Y|U)\right)=g\left(\mu\right)=\eta=UB$

expressed in terms of $Y,U,B$ where

$Y=\left[\begin{array}[]{c}X_{T}^{\prime}\\\ X_{T-1}^{\prime}\\\ \vdots\\\
X_{D}^{\prime}\end{array}\right]\quad
U=\left[\begin{array}[]{cccc}1&X_{T-1}^{\prime}&\cdots&X_{T-D}^{\prime}\\\
1&X_{T-2}^{\prime}&\cdots&X_{T-D-1}^{\prime}\\\ \vdots&\vdots&\ddots&\vdots\\\
1&X_{D-1}^{\prime}&\cdots&X_{0}^{\prime}\\\ \end{array}\right]\quad
B=\left[\begin{array}[]{c}\nu^{\prime}\\\ A_{1}^{\prime}\\\ \vdots\\\
A_{D}^{\prime}\end{array}\right]$

The core model relationship can be expressed in univariate terms by applying
the vectorization transformations

$\mathcal{Y}\stackrel{{\scriptstyle\text{def}}}{{=}}\text{vec}Y\;,\quad\mathcal{X}\stackrel{{\scriptstyle\text{def}}}{{=}}U\otimes
I_{M}\;,\quad\text{and
}\beta\stackrel{{\scriptstyle\text{def}}}{{=}}\text{vec}B$

with the result

$\mathcal{Y}=\left[\begin{array}[]{c}x_{1}\\\ \hdashline x_{2}\\\
\hdashline\vdots\\\ \hdashline
x_{M}\end{array}\right]\qquad\mathcal{X}=\left[\begin{array}[]{c:c:c:c}U&0&\cdots&0\\\
\hdashline 0&U&\cdots&0\\\ \hdashline\vdots&\vdots&\ddots&\vdots\\\ \hdashline
0&0&\cdots&U\end{array}\right]\qquad\beta=\text{vec}(B)$

where $x_{m}=(x_{0,m}\;\cdots\;x_{T,m})^{\prime}$ denote the univariate series
in each component so that

$g\left(\mathbb{E}(\mathcal{Y}|\mathcal{X})\right)=g\left(\text{vec}(\mu)\right)=\text{vec}(\eta)=\mathcal{X}\beta$

It is straightforward to verify by inspection that if the the elements of
$\mathcal{Y}$ are assumed to be Poisson conditional on $\mathcal{X}$ with the
corresponding mean, the likelihoods of $\beta$ and $B$ are identical, i.e.,

$\ell(\beta;\mathcal{Y},\mathcal{X})=\ell(B;Y,U)$

since the likelihood contributions from each $x_{t,m}$ are exactly the same.
Therefore, the estimate $\hat{B}$ can be recovered from the estimate

$\hat{\beta}=\operatorname{argmin}_{\beta}\left\\{-\ell(\beta;\mathcal{Y},\mathcal{X})+\lambda
P^{*}(\beta)\right\\}$

where $P^{*}(\beta)=P(B)$.

Due to the sparsity pattern in $\mathcal{X}$, the computations involved in the
update rule for the $l$th coordinate in the vectorized problem can be
articulated as operations involving only $U$ and columns of $Y$. Specifically,
denoting the row and column dimensions of $U$ as $N=T-D$ and $P=DM+1$ and the
row and column dimensions of $\mathcal{X}$ as $K=NM$ and $L=M(DM+1)$, direct
application of the calculations for the coordinate update for the $l$th
position in Algorithm A.2 to the vectorized problem in terms of
$\mathcal{Y},\mathcal{X}$ gives

$\displaystyle\left(\mathcal{Z}-\mathcal{Z}^{(l)}\right)^{\prime}\mathcal{W}\mathcal{X}_{l}$
$\displaystyle=\left(x_{m}-x_{m}^{(p)}\right)^{\prime}W_{m}u_{p}$
$\displaystyle\mathcal{X}_{l}^{\prime}\mathcal{W}\mathcal{X}_{l}$
$\displaystyle=u_{p}^{\prime}W_{m}u_{p}$

where $\mathcal{Z}$ is the working response vector, $\mathcal{W}$ is the
diagonal matrix of update weights, and the index transformation is given by

$\displaystyle(n,m)$
$\displaystyle=\left(k\text{mod}M+M\mathbbm{1}\\{k\text{mod}M=0\\},\frac{k-k\text{mod}M}{M}+\mathbbm{1}\\{k=k\text{mod}M\\}\right)$
$\displaystyle(p,m)$
$\displaystyle=\left(l\text{mod}P+P\mathbbm{1}\\{l\text{mod}P=0\\},\frac{l-l\text{mod}P}{P}+\mathbbm{1}\\{l=l\text{mod}P\\}\right)$

Note that the computation of $x_{mn}^{(p)}=\sum_{j\neq p}u_{nj}b_{jm}$ depends
only on the $m$th column of $B$.

Consequently, a more efficient strategy is to divide the coordinate iterations
into blocks corresponding to the blocks of $B$ and compute the update rules
using data stored in the formats specified by $U$ and $Y$. A slight ridge
penalty is included for numerical stability. This is given as Algorithm A.3
and that which is used in the simulation experiments and data analyses
presented in the main paper.

Algorithm A.3 Coordinatewise IRLS algorithm for computing LASSO GVAR(D)
estimates with blockwise updates.

0:

Data $Y,U$, hyperparameters $\lambda,\alpha$, and convergence threshold
$\delta$

Initialize $\hat{B}$

while $\|\nabla\ell(\hat{B};X)\|_{2}>\delta$ do

for $m=1,\dots,M$ do

for $p=1,\dots,DM+1$ do

$W_{m}\longleftarrow\left(\text{diag}(g^{\prime}(\hat{\mu}_{m})\right)^{-1}$

$z_{m}\longleftarrow U\hat{b}_{m}-(x_{m}-\hat{\mu}_{m})W_{m}^{-1}$

$z_{mn}^{(p)}\longleftarrow\sum_{j\neq p}u_{nj}\hat{b}_{jm}$ for $n=1,\dots,N$

$\hat{b}_{pm}\longleftarrow\begin{cases}\frac{\left(z_{m}-z_{m}^{(p)}\right)^{\prime}W_{m}u_{p}}{u_{p}^{\prime}W_{m}u_{p}}&\text{,
if }p=1\\\
\frac{S\left(\left(z_{m}-z_{m}^{(p)}\right)^{\prime}W_{m}u_{p},\lambda\alpha\right)}{u_{p}^{\prime}W_{m}u_{p}+\lambda(1-\alpha)}&\text{,
if }p>1\end{cases}$

end for

end for

end while

$\hat{B}$


[cite=DBLP:conf/wdag/LeinweberH23]Publisher
Editionhttps://doi.org/10.4230/LIPIcs.DISC.2023.43 Rotem Oshman 1 37th
International Symposium on Distributed Computing (DISC 2023) DISC 2023 DISC
2023 October 10-12, 2023 L’Aquila, Italy 281 40 Institute of Information
Security and Dependability (KASTEL), Karlsruhe Institute of Technology (KIT),
Germanymarc.leinweber@kit.eduhttps://orcid.org/0000-0002-9638-8526 Institute
of Information Security and Dependability (KASTEL), Karlsruhe Institute of
Technology (KIT),
Germanyhannes.hartenstein@kit.eduhttps://orcid.org/0000-0003-3441-3180 Marc
Leinweber and Hannes Hartenstein Security and privacy Distributed systems
security This work was supported by funding from the topic Engineering Secure
Systems of the Helmholtz Association (HGF).

# Let It TEE: Asynchronous Byzantine Atomic Broadcast with $n\geq 2f+1$

Marc Leinweber  Hannes Hartenstein

###### Abstract

Asynchronous Byzantine Atomic Broadcast (ABAB) promises simplicity in
implementation as well as increased performance and robustness in comparison
to partially synchronous approaches. We adapt the recently proposed DAG-Rider
approach to achieve ABAB with $n\geq 2f+1$ processes, of which $f$ are faulty,
with only a constant increase in message size. We leverage a small Trusted
Execution Environment (TEE) that provides a unique sequential identifier
generator (USIG) to implement Reliable Broadcast with $n>f$ processes and show
that the quorum-critical proofs still hold when adapting the quorum size to
$\lfloor\frac{n}{2}\rfloor+1$. This first USIG-based ABAB preserves the
simplicity of DAG-Rider and serves as starting point for further research on
TEE-based ABAB.

###### category:

Brief Announcement

###### keywords:

Byzantine Fault Tolerance, Trusted Execution Environments, Asynchrony

##  1 Introduction

Atomic Broadcast primitives play a crucial role for Byzantine-fault tolerant
(BFT) State Machine Replication (SMR). A prominent example for BFT SMR in the
partially synchronous model is PBFT [6]. By use of small Trusted Execution
Environments (TEE) that generate and sign unique sequential identifiers on
each process, called USIGs, Veronese et al. [13] showed that PBFT’s
communication complexity can be reduced and the fault tolerance can be
increased to $n\geq 2f+1$ while still tolerating Byzantine faults. The
authenticity/integrity of TEEs can be verified remotely and, thus, TEEs are
assumed to only fail by crashing. However, as shown by Miller et al. [12],
Asynchronous Byzantine Atomic Broadcast (ABAB) outperforms approaches based on
the partially synchronous model particularly under faults and tends to show a
simpler design. While it is known that any asynchronous crash fault-tolerant
algorithm can be compiled to withstand Byzantine faults using TEEs [3, 7], the
proposed compilers show either a polynomial or an exponential overhead in
runtime. We are interested in a simple and straightforward design of a USIG-
enhanced ABAB that does not add further message rounds and only adds a
constant number of bits to each message (essentially a counter value and a
signature). To this end, we adapt DAG-Rider [10] to provide ABAB with $n\geq
2f+1$ processes. We give a quick recap on DAG-Rider and explain the adaption
_TEE-Rider_. Besides using TEE-based Reliable Broadcast and changing the
required quorums from $2f+1$ to $\lfloor\frac{n}{2}\rfloor+1$, we leave DAG-
Rider unchanged. We show that the quorum-based arguments of DAG-Rider still
hold for TEE-Rider.

##  2 TEE-Rider: Transforming DAG-Rider to $n\geq 2f+1$

We make use of the following definition of Atomic Broadcast for a set of
processes $P,n:=|P|$. The processes communicate over authenticated point-to-
point links with eventual delivery.

######  Definition 2.1 (Atomic Broadcast).

Each process $p_{i}\in P$ receives client transactions $t$ via events
$\mathrm{clientRequest}(t)$. Correct processes deliver tuples $(t,r,p_{i})$,
where $t$ is a client transaction, $r\in\mathbb{N}_{0}$ a round number, and
$p_{i}\in P$ the process that initially received $t$, satisfying the following
properties:
Agreement: If a correct process $p_{i}\in P$ delivers $(t,r,p_{j})$, then
every other correct process $p_{k}\in P,k\not=i$ eventually delivers
$(t,r,p_{j})$ with probability 1.
Integrity: For each round $r\in\mathbb{N}_{0}$ and process $p_{j}\in P$, a
correct process $p_{i}\in P$ delivers $(t,r,p_{j})$ at most once.
Validity: If a correct process $p_{i}\in P$ receives an event
$\mathrm{clientRequest}(t)$, then every correct process $p_{k}\in P$
eventually delivers $(t,r,p_{i})$ with probability 1.
Total Order: Let $m_{1}$ and $m_{2}$ be any two valid tuples that are
delivered by any two correct processes $p_{i},p_{j}\in P$. If $p_{i}$ delivers
$m_{1}$ before $m_{2}$, then $p_{j}$ delivers $m_{1}$ before $m_{2}$.

###  2.1 Changes in Assumptions, Building Blocks, and Setup

In addition to the assumptions of DAG-Rider, we assume that each process is
equipped with a USIG [13] that may only fail by crashing. It implements a
signature service that binds a unique counter value to each signature it
produces. The USIG is used for Reliable Broadcast with a fault tolerance of
$n>f$ as, e.g., implemented in [8, Algorithm 1]: it is a ‘single echo’
algorithm with USIG-signed messages and attached counter value. Correct
processes relay a message once and reject messages with invalid USIG
signatures or with counter values already received which prevents equivocating
messages for the same counter. An instance of the Reliable Broadcast
abstraction has two functions: $\mathrm{broadcast}(r,m)$ to reliably broadcast
exactly one arbitrary message $m$ for round $r$ to all processes in $P$, and
$\mathrm{delivered}()$ which returns all messages that were received by the
instance since the last call to $\mathrm{delivered}()$. We expect the Reliable
Broadcast abstraction to fullfil the following properties:

######  Definition 2.2 (Reliable Broadcast).

A sender $p_{s}\in P,n:=|P|$ can call $\mathrm{broadcast}(m)$. Correct
processes deliver $(c,m)$ where $c\in\mathbb{N}_{0}$ and $m$ an arbitrary
message satisfying the following properties:
RB-Agreement. If a correct process $p_{i}\in P$ delivers $(c,m)$, then every
other correct process $p_{k}\in P,k\not=i$ eventually delivers the same
$(c,m)$.
RB-Integrity. For each $c\in\mathbb{N}_{0}$, a correct process $p_{k}\in P$
delivers $(c,m)$ at most once.
RB-Validity. If a correct sender calls $\mathrm{broadcast}(m)$, then every
correct process $p_{i}\in P$ eventually delivers $(c,m)$.

Additionally, we assume an asynchronous common coin, e.g. as defined by Cachin
et al. [5], that produces a uniformly distributed common random number $p$ out
of $\\{p\mid p\in\mathbb{N}_{0}\colon p<n\\}$ for all correct processes and a
name $i\in\mathbb{N}_{0}$ as soon as $f+1$ processes invoked
$\mathrm{toss}(i)$; repetitive calls with same the $i$ yield the same $p$. We
further assume a trusted setup of the common coin and the USIGs using a public
key infrastructure (to set up the common coin’s threshold signature scheme,
dealerless variants [4] and those with an asynchronous setup [1] exist).

###  2.2 The Algorithm

Figure 1: Example DAG for $n=5$ processes of which at maximum $2$ may be
faulty. Shown is the ‘global’ state of the graph, i.e., after every process
eventually received every vertex. For simplicity, all vertices are valid and
weak edges are left out. The direct commit rule is not fulfilled for any
process for wave $w$; it is fulfilled for processes $p_{3}$, $p_{4}$, and
$p_{5}$ for wave $w+1$. The green coloring highlights the effect of the direct
commit rule as proven in Lemma 3.1. The direct commit rule ensures that a
correct process can commit a wave retrospectively if it was not able to commit
when it finished the wave. Since the leader vertex of wave $w+1$ has a strong
path to the leader vertex of wave $w$, wave $w$ will be committed
retrospectively.

DAG-Rider uses $n$ Reliable Broadcast instances to disseminate process
messages and to construct a directed acyclic graph (DAG) that captures the
communication history of all processes. In a second step, each process derives
consensus on the order of transactions using a Common Coin based on the graph
structure. The adapted DAG-Rider algorithm executed by a correct process
$p_{i}\in P$ is shown in Algorithm 1. Utility functions are listed in
Algorithm 2. The core of the approach is the construction and interpretation
of a (local) DAG that captures received transactions and the observed
communication sequence between processes. The DAG is structured in rounds and
a round contains at maximum one vertex per process, i.e., $n$ vertices. Rounds
are addressed in an array style and the local view of a process $p_{i}$ on the
DAG is indicated by an index $i$. The very first round $\mathit{DAG}_{i}[0]$
is initialized with $n$ hard-coded “genesis” vertices. A vertex in round $r$
has two types of edges: strong edges point to vertices of round $r-1$ and weak
edges point to vertices of any round $r^{\prime}\leq r-2$. As soon as $p_{i}$
received $\lfloor\frac{n}{2}\rfloor+1$ _valid_ vertices for a round $r$, i.e.,
$\lfloor\frac{n}{2}\rfloor+1$ vertices referencing
$\lfloor\frac{n}{2}\rfloor+1$ vertices of round $r-1$ as strong edges
($v.\mathit{strong}$, l. 11), for which it also knows its predecessors (l.
15), $p_{i}$ will _complete_ round $r$ and transition to round $r+1$. Now, as
soon as $p_{i}$ receives a client transaction, it will become the payload of a
vertex $v$ which is created and broadcast by $p_{i}$ for round $r+1$ (ll. 44
and 21-27). The vertex $v$ connects to all vertices $p_{i}$ received for round
$r$ (l. 23). If $p_{i}$ received vertices $u$ for older rounds that are not
reachable from the newly created vertex using the transitive closure of strong
and weak edges (a ‘path’), $u$ will become a weak edge of $v$
($v.\mathit{weak}$, l. 26). The new vertex is broadcast using Reliable
Broadcast instance $i$ to all processes (l. 27). Every fourth round a so-
called wave, consisting of four rounds, is completed (l. 19) and the DAG
structure is used to derive a total order on the transactions (ll. 28-42).
Each wave $w$ has exactly one wave leader $v$ which is chosen calling
$\mathit{coin}.\mathrm{toss}(w)$ from the vertices of $w$’s first
$\mathrm{round}(w,1)$. The random number is used to select the process whose
vertex is to be used as wave leader. If $v$ was not (yet) received or there
are no $\lfloor\frac{n}{2}\rfloor+1$ vertices in the $w$’s fourth
$\mathrm{round}(w,4)$ that have $v$ in their transitive closure of strong
edges (a ‘strong path’), i.e. the _direct commit rule_ is not fulfilled, the
wave cannot be committed (l. 30). If wave $w$ can be committed, process
$p_{i}$ checks first if there are wave leaders of waves $w^{\prime}$ between
the last wave that was committed (variable $\mathit{decidedWave}$) and the
current wave $w$ that were received in the meantime and are connected to the
leader of the wave $w^{\prime}+1$ (retrospective commit, ll. 33-36). The wave
leaders are used as the root for a deterministic graph traversal to determine
the total order of transactions (ll. 38-42). An example for a resulting graph
with $n=5$ processes, i.e. $f\leq 2$, is shown in Figure 1.

Algorithm 1 TEE-Rider pseudocode for process $p_{i}\in
P,n:=|P|,\boldsymbol{n\geq 2f+1}$

1:state $\mathit{DAG}\colon$ array of sets of vertices, $\mathit{DAG}[0]$
initialized with “genesis” vertices

2:state $r\colon\mathbb{N}$, initialized with 1

3:state $\mathit{decidedWave}\colon\mathbb{N}_{0}$, initialized with 0

4:state $\mathit{transactionsToPropose}\colon$ queue of client transactions
$t$, initialized empty

5:state $\mathit{buffer}\colon$ set of vertices, initialized empty

6:state $\mathit{rb}\colon$ array of $n$ Reliable Broadcast instances with
delivered() and broadcast($r,m$)

7:state $\mathit{coin}\colon$ common coin instance with toss($w$)

8:while True do

9: for $k\leftarrow 0$ up to $n-1$ do

10: for $m=(r^{\prime},v)\in\mathit{rb}[k].\mathrm{delivered}()$ do

11: if $|v.strong|<\lfloor\frac{n}{2}\rfloor+1$ then continue

12: $v.\mathit{source}\leftarrow p_{k};v.\mathit{round}\leftarrow
r^{\prime};v.\mathit{delivered}\leftarrow\mathrm{False}$

13: $\mathit{buffer}.\mathrm{add}(v)$

14: for $v\in\mathit{buffer}$ do

15: if $v.\mathit{round}>r\lor\exists u\in v.\mathit{strong}\cup
v.\mathit{weak}\colon u\not\in\cup_{r^{\prime}\geq 0}\mathit{DAG}[r^{\prime}]$
then continue

16: $\mathit{DAG}[v.\mathit{round}].\mathrm{add}(v)$

17: $\mathit{buffer}.\mathrm{remove}(v)$

18: if $|\mathit{DAG}[r]|<\lfloor\frac{n}{2}\rfloor+1$ then continue

19: if $r\mod 4=0$ then $\mathrm{waveReady}(\frac{r}{4})$

20: $r\leftarrow r+1$

21: wait until $\neg\mathit{transactionsToPropose}.\mathrm{isEmpty}($)

22: $v\leftarrow$ new vertex

23:
$v.\mathit{block}\leftarrow\mathit{transactionsToPropose}.\mathrm{dequeue}();v.\mathit{strong}\leftarrow\mathit{DAG}[r-1]$

24: for $r^{\prime}\leftarrow r-2$ down to $1$ do

25: for $u\in\mathit{DAG}[r^{\prime}]$ do

26: if $\neg\mathrm{path}(v,u)$ then $v.\mathit{weak}.\mathrm{add}(u)$

27: $\mathit{rb}[i].\mathrm{broadcast}(r,v)$

28:function _waveReady_($w$)

29: $v\leftarrow\mathrm{getWaveLeader}(w)$ $\triangleright$ $\bot$ if
$\mathrm{round}(w,1)$ vertex of chosen process is not in $\mathit{DAG}$

30: if $v=\bot\lor|\\{u\mid
u\in\mathit{DAG}[\mathrm{round}(w,4)]\colon\mathrm{strongPath}(u,v)\\}|<\lfloor\frac{n}{2}\rfloor+1$
then return

31: $\mathit{leadersStack}\leftarrow$ new stack;
$\mathit{leadersStack}.\mathrm{push}(v)$

32: for $w^{\prime}\leftarrow w-1$ down to $\mathit{decidedWave}+1$ do

33: $u\leftarrow\mathrm{getWaveLeader}(w^{\prime})$

34: if $u\not=\bot\land\mathrm{strongPath}(v,u)$ then

35: $\mathit{leadersStack}.\mathrm{push}(u);v\leftarrow u$

36: $\mathit{decidedWave}\leftarrow w$

37: while $\neg\mathit{leadersStack}.\mathrm{isEmpty}()$ do

38: $v\leftarrow\mathit{leadersStack}.\mathrm{pop}()$

39: $\mathit{verticesToDeliver}\leftarrow\\{u\mid
u\in\cup_{r^{\prime}>0}\mathit{DAG}[r^{\prime}]\colon\mathrm{path}(v,u)\land\neg
u.\mathit{delivered}\\}$

40: for $u\in\mathit{verticesToDeliver}$ in deterministic order do

41: $u.\mathit{delivered}\leftarrow\mathrm{True}$; deliver
$(u.\mathit{block},u.\mathit{round},u.\mathit{source})$

42:upon clientRequest ($t$)

43: $\mathit{transactionsToPropose}.\mathrm{enqueue}(t)$

Algorithm 2 Utility functions pseudocode

1:function _path_($v,u$) : boolean

2: return exists a sequence of vertices
$(v_{1},v_{2},...,v_{k})\in\cup_{r^{\prime}\geq 0}DAG[r^{\prime}]$ such that

3: $v_{1}=v\land v_{k}=u\land\forall i\in[2,k]\colon v_{i}\in
v_{i-1}.strong\cup v_{i-1}.weak$

4:function _strongPath_($v,u$) : boolean

5: return exists a sequence of vertices
$(v_{1},v_{2},...,v_{k})\in\cup_{r^{\prime}\geq 0}DAG[r^{\prime}]$ such that

6: $v_{1}=v\land v_{k}=u\land\forall i\in[2,k]\colon v_{i}\in v_{i-1}.strong$

7:function _getWaveLeader_($w$) : vertex or $\bot$

8: $p_{j}\leftarrow coin.\mathrm{toss}(w)$

9: if $\exists v\in DAG[\mathrm{round}(w,1)]\colon v.source=p_{j}$ then return
$v$

10: return $\bot$

11:function _round_($w,i$) : $\mathbb{N}$

12: return $4(w-1)+i$

##  3 Correctness Argument

Lemmas 1 and 2 of the original DAG-Rider publication [10] are crucial for
Total Order and Agreement and rely on quorum intersection arguments. The
following Lemmas 3.1 and 3.3 show the corresponding results for a quorum size
of $\lfloor\frac{n}{2}\rfloor+1$. Results for Integrity and Validity simply
follow from the original paper.

######  Lemma 3.1.

If a correct process $p_{i}\in P$ commits the wave leader $v$ of a wave $w$
when it completes wave $w$ in $\mathrm{round}(w,4)$, then any valid vertex
$v^{\prime}$ of any process $p_{j}\in P$ broadcast for a round
$r\geq\mathrm{round}(w+1,1)$ will have a strong path to $v$.

######  Proof 3.2.

Since $p_{i}$ commits $v$ in $\mathrm{round}(w,4)$, the direct commit rule is
fulfilled (l. 30): $\exists
U\subseteq\mathit{DAG}_{i}[\mathrm{round}(w,4)]\colon|U|\geq\lfloor\frac{n}{2}\rfloor+1\land\forall
u\in U\colon\mathrm{strongPath}(u,v)$. A valid vertex must reference at least
$\lfloor\frac{n}{2}\rfloor+1$ distinct vertices of the previous round with a
strong edge (l. 11). Thus, a process $p_{j}\in P$ broadcasting a valid vertex
$v_{j}$ for $\mathrm{round}(w+1,1)$ selected at least
$\lfloor\frac{n}{2}\rfloor+1$ vertices of $\mathrm{round}(w,4)$ as strong
edges for $v_{j}$. Any two subsets of size $\lfloor\frac{n}{2}\rfloor+1$ of a
superset of size $n$ intersect at least in one element. Thus, every valid
vertex of a process broadcast for $\mathrm{round}(w+1,1)$ must have at least
one edge to a vertex of $U$, and, via $U$ to $v$. As every valid vertex of
$\mathrm{round}(w+1,1)$ has a strong path to $v$, and every valid vertex of
$\mathrm{round}(w+1,2)$ connects to at least $\lfloor\frac{n}{2}\rfloor+1$
vertices of $\mathrm{round}(w+1,1)$, by induction, any valid vertex
$v^{\prime}$ of any process $p_{j}\in P$ broadcast for a round
$r\geq\mathrm{round}(w+1,1)$ has a strong path to $v$.

######  Lemma 3.3.

When a correct process $p_{i}\in P$ completes $\mathrm{round}(w,4)$ of wave
$w$, then $\exists
V_{1}\subseteq\mathit{DAG}_{i}[\mathrm{round}(w,1)],V_{4}\subseteq\mathit{DAG}_{i}[\mathrm{round}(w,4)]\colon|V_{1}|\geq\lfloor\frac{n}{2}\rfloor+1\land|V_{4}|\geq\lfloor\frac{n}{2}\rfloor+1\land(\forall
v_{1}\in V_{1},\forall v_{4}\in V_{4}\colon\mathrm{strongPath}(v_{4},v_{1}))$.

######  Proof 3.4.

By use of Reliable Broadcast and validity checks in ll. 11 and 15, faulty
processes are limited to omission faults. Thus, the _get-core_ argument of
Attiya and Welch [2, Sec. 14.3.1] still holds [2, Sec. 14.3.3]: Let
$A\in\\{0,1\\}^{n\times n}$ be a matrix that contains a row for each possible
vertex of $\mathrm{round}(w,3)$ and a column for each possible vertex of
$\mathrm{round}(w,2)$. Let $A[j,k]=1$ if the vertex of process $p_{j}$ of
$\mathrm{round}(w,3)$ has a strong edge to the vertex of process $p_{k}$ of
$\mathrm{round}(w,2)$ or $p_{j}$ sends no vertex (or an invalid one) but
$p_{k}$ sends a valid vertex for $\mathrm{round}(w,2)$. As there are at least
$\lfloor\frac{n}{2}\rfloor+1\leq n-f$ correct processes, each row of $A$
contains at least $\lfloor\frac{n}{2}\rfloor+1$ ones and $A$ contains at least
$n(\lfloor\frac{n}{2}\rfloor+1)$ ones. Since there are $n$ columns, there must
be a column $l$ with at least $\lfloor\frac{n}{2}\rfloor+1$ ones. This implies
there is a vertex $v_{l}$ by process $p_{l}$ in $\mathrm{round}(w,2)$ s.t.
$\exists
V_{3}\subseteq\mathit{DAG}_{i}[\mathrm{round}(w,3)]\colon|V_{3}|\geq\lfloor\frac{n}{2}\rfloor+1\land\forall
v_{3}\in V_{3}\colon\mathrm{strongPath}(v_{3},v_{l})$. As at most $f$ vertices
in $V_{3}$ belong to faulty processes that may commit send omission faults for
$\mathrm{round}(w,3)$ and $\lfloor\frac{n}{2}\rfloor+1\geq f+1$, by quorum
section at least one vertex of $V_{3}$ is received by any correct process
$p_{j}\in P$ before it sends its vertex for $\mathrm{round}(w,4)$. Thus, every
valid vertex in $\mathit{DAG}_{i}[\mathrm{round}(w,4)]$ has at least one
strong edge to a vertex of $V_{3}$. Since $v_{l}$ must be valid and thus has a
strong edge to each vertex of a set
$V_{1}\subseteq\mathit{DAG}_{i}[\mathrm{round}(w,1)],|V_{1}|\geq\lfloor\frac{n}{2}\rfloor+1$,
any valid vertex of rounds $r\geq\mathrm{round}(w,4)$ has a strong path to
every vertex, including $V_{1}$, reached by $v_{l}$ via strong paths. Please
note that the construction of the set $V_{1}$ is valid for all correct
processes that complete the wave and, thus, represents the ‘common core’.

##  4 Discussion and Conclusion

DAG-Rider shows the power of causal order broadcast to implement consensus.
The adaption for TEEs preserves the simplicity of DAG-Rider while increasing
its fault tolerance and reducing the communication effort (i.e., from ‘double
echo‘ to ‘single echo‘ Reliable Broadcast). The ease of adaption of DAG-Rider
for TEEs make it a perfect textbook example for TEE-based ABAB. Follow-up work
to DAG-Rider, Tusk [9], addresses a major deployability issue, namely garbage
collection, shortens the wave length, and replaces the underlying Reliable
Broadcast with a communication scheme that leverages the graph structure to
achieve linear communication complexity in the happy case. Additionally, to
the best of our knowledge, there exists no TEE-based, dealerless, and
asynchronous common coin primitive. In summary, investigating a TEE-based
dealerless setup as well as transforming the follow-ups of DAG-Rider for
empirical studies to investigate the assumed superiority of asynchronous TEE-
based approaches, e.g., in comparison to MinBFT [13], is a promising line of
research.

## References

  * [1] Ittai Abraham, Marcos Kawazoe Aguilera, and Dahlia Malkhi.  Fast asynchronous consensus with optimal resilience.  In Nancy A. Lynch and Alexander A. Shvartsman, editors, Distributed Computing, volume 6343 of Lecture Notes in Computer Science, pages 4–19. Springer, 2010.  doi:10.1007/978-3-642-15763-9\\_3.
  * [2] Hagit Attiya and Jennifer L. Welch.  Distributed computing - fundamentals, simulations, and advanced topics (2. ed.).  Wiley, 2004.  doi:10.1002/0471478210.
  * [3] Naama Ben-David, Benjamin Y. Chan, and Elaine Shi.  Revisiting the power of non-equivocation in distributed protocols.  In Alessia Milani and Philipp Woelfel, editors, PODC ’22: ACM Symposium on Principles of Distributed Computing, Salerno, Italy, July 25 - 29, 2022, pages 450–459. ACM, 2022.  doi:10.1145/3519270.3538427.
  * [4] Alexandra Boldyreva.  Threshold signatures, multisignatures and blind signatures based on the gap-diffie-hellman-group signature scheme.  In Yvo Desmedt, editor, Public Key Cryptography — PKC 2003, volume 2567 of Lecture Notes in Computer Science, pages 31–46. Springer, 2003.  doi:10.1007/3-540-36288-6\\_3.
  * [5] Christian Cachin, Klaus Kursawe, and Victor Shoup.  Random oracles in constantinople: Practical asynchronous byzantine agreement using cryptography.  J. Cryptol., 18(3):219–246, 2005.  doi:10.1007/s00145-005-0318-0.
  * [6] Miguel Castro and Barbara Liskov.  Practical byzantine fault tolerance and proactive recovery.  ACM Trans. Comput. Syst., 20(4):398–461, 2002.  doi:10.1145/571637.571640.
  * [7] Allen Clement, Flavio Junqueira, Aniket Kate, and Rodrigo Rodrigues.  On the (limited) power of non-equivocation.  In Darek Kowalski and Alessandro Panconesi, editors, ACM Symposium on Principles of Distributed Computing, PODC ’12, Funchal, Madeira, Portugal, July 16-18, 2012, pages 301–308. ACM, 2012.  doi:10.1145/2332432.2332490.
  * [8] Miguel Correia, Giuliana Santos Veronese, and Lau Cheuk Lung.  Asynchronous byzantine consensus with 2f+1 processes.  In Sung Y. Shin, Sascha Ossowski, Michael Schumacher, Mathew J. Palakal, and Chih-Cheng Hung, editors, Proceedings of the 2010 ACM Symposium on Applied Computing (SAC), Sierre, Switzerland, March 22-26, 2010, pages 475–480. ACM, 2010.  doi:10.1145/1774088.1774187.
  * [9] George Danezis, Lefteris Kokoris-Kogias, Alberto Sonnino, and Alexander Spiegelman.  Narwhal and Tusk: a DAG-based mempool and efficient BFT consensus.  In Yérom-David Bromberg, Anne-Marie Kermarrec, and Christos Kozyrakis, editors, EuroSys ’22: Seventeenth European Conference on Computer Systems, Rennes, France, April 5 - 8, 2022, pages 34–50. ACM, 2022\.  doi:10.1145/3492321.3519594.
  * [10] Idit Keidar, Eleftherios Kokoris-Kogias, Oded Naor, and Alexander Spiegelman.  All you need is DAG.  In Avery Miller, Keren Censor-Hillel, and Janne H. Korhonen, editors, PODC ’21: ACM Symposium on Principles of Distributed Computing, Virtual Event, Italy, July 26-30, 2021, pages 165–175. ACM, 2021\.  doi:10.1145/3465084.3467905.
  * [11] Marc Leinweber and Hannes Hartenstein.  Brief Announcement: Let It TEE: Asynchronous Byzantine Atomic Broadcast with n $\geq$ 2f+1.  In Rotem Oshman, editor, 37th International Symposium on Distributed Computing, DISC 2023, October 10-12, 2023, L’Aquila, Italy, volume 281 of LIPIcs, pages 43:1–43:7. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2023.  doi:10.4230/LIPICS.DISC.2023.43.
  * [12] Andrew Miller, Yu Xia, Kyle Croman, Elaine Shi, and Dawn Song.  The honey badger of BFT protocols.  In Edgar R. Weippl, Stefan Katzenbeisser, Christopher Kruegel, Andrew C. Myers, and Shai Halevi, editors, Proc. 2016 ACM SIGSAC Conf. on Computer and Communications Security, Vienna, Austria, 2016, pages 31–42. ACM, 2016.  doi:10.1145/2976749.2978399.
  * [13] Giuliana Santos Veronese, Miguel Correia, Alysson Neves Bessani, Lau Cheuk Lung, and Paulo Veríssimo.  Efficient byzantine fault-tolerance.  IEEE Trans. Computers, 62(1):16–30, 2013.  doi:10.1109/TC.2011.221.