# “Generalized” generalized linear models:
Convex estimation and online bounds

Anatoli Juditsky LJK, Université Grenoble Alpes, B.P. 53, 38041 Grenoble Cedex
9, France<EMAIL_ADDRESS>Arkadi Nemirovski Georgia Institute of
Technology, Atlanta, Georgia 30332, USA<EMAIL_ADDRESS>Yao Xie
Georgia Institute of Technology, Atlanta, Georgia 30332, USA,
<EMAIL_ADDRESS>Chen Xu Georgia Institute of Technology, Atlanta,
Georgia 30332, USA<EMAIL_ADDRESS>
The work of Xu and Xie is partially supported by an NSF CAREER CCF-1650913,
and NSF DMS-2134037, CMMI-2015787, CMMI-2112533, DMS-1938106, and DMS-1830210.

###### Abstract

We introduce a new computational framework for estimating parameters in
generalized generalized linear models (GGLM), a class of models that extends
the popular generalized linear models (GLM) to account for dependencies among
observations in spatio-temporal data. The proposed approach uses a monotone
operator-based variational inequality method to overcome non-convexity in
parameter estimation and provide guarantees for parameter recovery. The
results can be applied to GLM and GGLM, focusing on spatio-temporal models. We
also present online instance-based bounds using martingale concentrations
inequalities. Finally, we demonstrate the performance of the algorithm using
numerical simulations and a real data example for wildfire incidents.

##  1 Introduction

We aim to perform parameter estimation for generalized generalized linear
models (GGLM) and provide methods for quantifying the uncertainty in the
estimated parameters. GGLM can be thought of as an extension of the widely
used generalized linear models (GLM) [5, 4, 1], which offer a versatile set of
models for statistical analysis of data when we have a set of measurements or
counts accompanied by relevant contextual information. GLM encompasses a broad
range of models that map input/response variables to predictors, such as logit
and probit models for quantal response data and log-linear models for count
response data. However, standard GLM models assume that the observations are
independent, which precludes modeling of ”data exhibiting the autocorrelations
of time series and spatial processes,” as noted in [4]. GGLM extends GLM to
account for dependencies among observations in spatio-temporal data.

This paper presents a computational framework for the “generalized”
generalized linear model (GGLM) using a new convex formulation by solving the
monotone operator-based variational inequalities (VI) [3]. In particular, we
emphasize the use of GGLM for spatio-temporal models, although the results can
still be applied to GLM if we remove the spatial and temporal dependence. We
aim to build the bridge between the generalized general linear model and the
variational inequality in optimization. This leads to efficient computational
procedures that overcome the non-convexity in directly fitting the problem via
least-square formulations. In some cases, the new VI approach will also enable
performance guarantees for maximum likelihood estimates (MLE) for model
parameters.

The conventional approach for estimating parameters in GLM models involves
using weighted least-squares (WLS) [5, 4, 1]. Although WLS is commonly used in
practice and offers a convenient estimation strategy, it does not provide
precise guarantees for parameter recovery. In contrast, the current study
proposes a VI monotone operator recovery method that offers guaranteed
parameter recovery and online concentration bounds for GGLM. Notably, these
results are applicable to the classic GLM and its spatio-temporal extension.

Moreover, our approach avoids the non-convexity caused by non-linear link
functions. In order to construct a convex VI problem, the structural
assumptions to be imposed on the model are essentially weaker than those
resulting in convex ML-based or LS-based problems. The VI problem we solve may
or may not be equivalent to a convex minimization problem. For instance, we
will show that the VI problem we construct will reduce to the LS problem under
linear link functions. Moreover, the proposed GGLM emphasizes spatio-temporal
dependencies. Thus, the observations are no longer i.i.d. and can have complex
dependencies. Therefore, martingale concentrations are developed to quantify
the recovery error.

The rest of the paper is organized as follows. Section 2 introduces the GGLM
model, recovery algorithm, and performance metrics. Section 4 presents
performance guarantees for recoveries based on solving monotone operator-based
variational inequalities. Section 5 discusses how to use the framework in a
special case: spatio-temporal processes, followed by two representative
examples for Bernoulli processes (Section 5.1) and Poisson processes (Section
5.2). Finally, we conclude the paper with discussions.

Notation: The notations of the paper is standard. All linear spaces we work in
are $\ bR^{n}$, with different dimensions $n$, equipped with the standard
Euclidean structure. Sometimes the vectors from ${\mathbf{R}}^{n}$ in question
have additional structure – they are represented as one-, two-, or even three-
dimensional arrays of blocks of common size. In all cases, we assume that the
entries in a vector from ${\mathbf{R}}^{n}$ under consideration independently
of the additional block structure, if any, are assigned serial numbers. This
allows us to write inner products as $x^{T}y$, to identify a linear mapping
from ${\mathbf{R}}^{n}$ into ${\mathbf{R}}^{m}$ with $m\times n$ matrix, its
conjugate – with the transpose of this matrix, etc. Let $[x]_{i}$ denote the
$i$th entry of a vector $x$, and $[X]_{ij}$ denote the $(i,j)$th entry of a
matrix $X$; $I_{m}$ denotes an identity matrix of size $m\times m$.

##  2 “Generalized” Generalized Linear Models

We start by presenting the main modeling framework, which we will demonstrate
using examples of its wide applicability. In statistics, the generalized
linear model (GLM) [4] is a flexible generalization of the Gaussian model:
where the response variable is related to the predictor by a linear function
contaminated by additive Gaussian noise. The GLM generalizes linear regression
by allowing the response variable to be related to the predictors via a non-
linear link function, which determines the mean and variance (dispersion) of
the response.

###  2.1 Situation and goal

Consider the “Generalized” Generalized Linear Models (GGLM) as follows

> At time $t=1,2,...$ we observe random pair $\upsilon_{t}=(\zeta_{t},H_{t})$,
> where random vector $\zeta_{t}$ is the response that takes values in a known
> subset ${{\cal Z}}_{t}$ of some ${\mathfrak{Z}}_{t}={\mathbf{R}}^{n_{t}}$,
> $H_{t}$ is random $\kappa\times n_{t}$ matrix that contains the predictor
> variables taking values in a known subset ${{\cal H}}_{t}$ of the space
> ${\mathbf{R}}^{\kappa\times n_{t}}$ of $\kappa\times n_{t}$ matrices.
> Similar to the GLM model, we specify the (conditional) mean of the response
> variable. For every $t=1,2,...$ the conditional,
> $\upsilon^{t-1}=(\upsilon_{1},...,\upsilon_{t-1})$ (assuming
> $\upsilon^{0}=0$) and $H_{t}$ given, distribution of $\zeta_{t}$ has
> expectation
>
>
> ${\mathbf{E}}_{|\upsilon^{t-1},H_{t}}\\{\zeta_{t}\\}=\Phi_{t}(H_{t}^{T}\beta),$
> (1)
>
> where
>
>   * •
>
> $\Phi_{t}(\cdot):{\mathrm{Dom}\,}\Phi_{t}\to{\mathfrak{Z}}_{t}$ — the $t$-th
> known link function — is a continuous monotone mapping defined on a closed
> convex domain ${\mathrm{Dom}\,}\Phi_{t}\subset{\mathfrak{Z}}_{t}$,
> monotonicity meaning that
>
> $[\Phi_{t}(x)-\Phi_{t}(y)]^{T}[x-y]\geq 0\,\,\forall
> x,y\in{\mathrm{Dom}\,}\Phi_{t};$
>
>   * •
>
> $\beta$ is the unknown vector of parameters taking values in a convex
> compact subset ${\cal B}$ of ${\mathfrak{B}}={\mathbf{R}}^{\kappa}$;
> $\kappa$ is the number of parameters.
>
>   * •
>
> ${{\cal H}}_{t}$ and ${\mathfrak{B}}$ are such that
> $H_{t}^{T}x\in{\mathrm{Dom}\,}\Phi_{t}$ for all $H_{t}\in{{\cal H}}_{t}$,
> $x\in{\mathfrak{B}}$, and $t$.
>
>

>
> We call the situation stationary, if $n_{t}$, ${\mathfrak{Z}}_{t}$, ${{\cal
> Z}}_{t}$, ${{\cal H}}_{t}$ and $\Phi_{t}$ are independent of $t$. Our goal
> is to recover $\beta$ from observations
> $\upsilon^{N}=(\upsilon_{1},...,\upsilon_{N})$, where $N$ is a given time
> horizon. A special case is when the mapping (or link function) $\Phi(\cdot)$
> is the same for each entry of $\zeta_{t}$, for instance, the sigmoid
> function $\phi(z)=1/(1+e^{-z})$, such that
> $\Phi_{t}(\zeta)=[\phi([\zeta]_{1}),\ldots,\phi([\zeta]_{n})]$ for
> $\zeta\in{\mathbf{R}}^{n}$ and $[\zeta]_{i}$ is the $i$th entry (see [4] for
> examples of link functions). Later on, we will also further assume that the
> entries of $\zeta_{t}$ are conditionally independent, given $\upsilon_{t}$
> and $H_{t}$.

Figure 1: Illustration of the Generalized Generalized Linear Model (GGLM),
where $\zeta_{t}$ and $H_{t}$ are a sequence of predictors and response
variables, respectively. The link function $\Phi_{t}$ is known and the
parameters for predictors $\beta$ are unknown.

###  2.2 Connection with GLM

The stationary version of this model resembles the generalized linear model
(GLM) [4], when the observations $\upsilon_{t}=(H_{t},\zeta_{t})$,
$t=1,2,\ldots$ are i.i.d., $H_{t}\in{\mathbf{R}}^{\kappa}$ corresponds to the
predictors (regressors), and $\zeta_{t}\in{\mathbf{R}}$ corresponds to the
response variable (in GLM it is typically a scalar.) GLM specifies the model
using a systematic component that ${\mathbf{E}}[\zeta]=\phi(H_{t}^{T}\beta)$
for a monotonic differentiable link function $\phi$ (induced by the
distribution of $\zeta$), and a random component assuming the response follows
a distribution that usually comes from an exponential family. Thus,
observations are realizations of a random variable $\zeta$, whose distribution
$\psi$ takes the form of

$\zeta\sim\psi(\zeta;\theta,\varrho)=\exp\\{(\zeta\theta-B(\theta))/A(\varrho)+D(\zeta,\varrho)\\},$
(2)

for some specific function $A(\cdot)$, $B(\cdot)$, and $D(\cdot)$; the natural
parameter is $\theta$ and the nuisance parameter is $\varrho$. If $\varrho$ is
known, this is a single-parameter exponential-family model.

It can be shown that the mean ${\mathbf{E}}[\zeta]=B^{\prime}(\theta)$ and
variance $\mbox{Var}(\zeta)=B^{\prime\prime}(\theta).$ Except for the Gaussian
distribution, the variance is related to the mean for all other exponential
family distributions. Thus, one can define GLM induced by an exponential
family distribution by specifying the conditional mean in the form of (1). GLM
typically assumes that $\theta$ linearly depends on the predictors
(regressors), and the link function can be chosen according to the assumed
distribution. Some examples are as follows: (i) Gaussian distribution
$\zeta\sim\mathcal{N}(\mu,\sigma^{2})$, $\theta=\mu$, $\varrho=\sigma^{2}$,
$A(\varrho)=\varrho$, $B(\theta)=\theta^{2}/2$,
$D(\zeta,\varrho)=-\frac{1}{2}\\{\zeta^{2}/\sigma^{2}+\log(2\pi\sigma^{2})\\}$.
The mean ${\mathbf{E}}[\zeta]=B^{\prime}(\theta)=\theta=\mu$ and the link
function is an identity (GLM reduces to the linear model in this case). (ii)
Bernoulli distribution $\zeta\in\\{0,1\\}$ with parameter $p\in[0,1]$. The
density function is given by

$\psi(\zeta;p)=\zeta^{p}(1-\zeta)^{1-p}=\exp\\{\zeta\theta-\log(1+e^{\theta})\\},$

where the natural parameter $\theta=\log\frac{p}{1-p}$, $\varrho=1$,
$A(\varrho)=1$, $B(\theta)=\log(1+e^{\theta})$;
${\mathbf{E}}[\zeta]=B^{\prime}(\theta)=1/(1+e^{-\theta})=p$ and
$\mbox{Var}[\zeta]=B^{\prime\prime}(\theta)=e^{\theta}/(1+e^{\theta})^{2}=p(1-p)$;
the link function $\phi(z)$ is a sigmoid function. (iii) Poisson distribution
$\zeta=0,1,2,\ldots$ with parameter $\lambda$. The density function is given
by

$\psi(\zeta;\lambda)=\lambda^{\zeta}e^{-\lambda}/\zeta!=\exp\\{\zeta\log\lambda-\lambda-\log(\zeta!)\\}$

with natural parameter $\theta=\log\lambda$, $B(\theta)=e^{\theta}$. Thus,
${\mathbf{E}}[\zeta]=B^{\prime}(\theta)=e^{\theta}=\lambda$ and
$\mbox{Var}[\zeta]=B^{\prime\prime}(\theta)=e^{\theta}=\lambda$, and the link
function is the log function.

##  3 Convex parameter recovery

The problem we address here is recovering the parameters $\beta$ of the GGLM
model from observations. To this end, we intend to use a Variational
Inequality (VI) oriented procedure; when the link functions are monotone (as
it was assumed when presenting a GGLM model), the resulting procedure is
computationally efficient. In contrast, the computation is not always
efficient for more traditional recovery procedures using classical loss
functions or based on Maximum Likelihood, which requires the concavity of the
log-likelihood function.

Monotone vector fields and variational inequalities. We start by introducing
the necessary background information.

######  Definition 3.1 (Monotone vector field.).

A vector field $F:{\cal X}\rightarrow{\mathbf{R}}^{N}$ defined on a nonempty
convex subset ${\cal X}$ of ${\mathbf{R}}^{N}$ is called monotone, if $\langle
F(x)-F(y),x-y\rangle\geq 0$ whenever $x,y\in{\cal X}$.

The set ${\cal X}$ and a monotone on ${\cal X}$ vector field $F$ specify
Variational Inequality

$\hbox{Find $z\in{\cal X}:$}\ \langle F(w),w-z\rangle\geq 0,\quad\forall
w\in{\cal X}.$ $None$

Vectors $z$ satisfying the requirement just defined are called weak solutions
to ${\mathrm{VI}}[F,{\cal X}]$. A vector $z\in{\cal X}$ is called a strong
solution to ${\mathrm{VI}}[F,{\cal X}]$ when

$\langle F(z),w-z\rangle\geq 0,\quad\forall w\in{\cal X}.$

From the monotonicity of $F$ it follows that every strong solution is a weak
one. When the monotone vector field $F$ is continuous, the weak solutions to
${\mathrm{VI}}[F,{\cal X}]$ are the same as strong solutions. Weak solutions
definitely exist whenever ${\cal X}$ is a nonempty convex compact set.

Some basic examples of monotone vector fields are
— the gradient field $\nabla f(x)$ of continuously differentiable convex
function $f:{\cal X}\to{\mathbf{R}}$;
— the vector field $F(u,v)=[\nabla_{u}f(u,v);-\nabla_{v}f(u,v)]$ of
continuously differentiable convex in $u$ and concave in $v$ function
$f(u,v):{\cal X}:=U\times V\to{\mathbf{R}}$.

In these cases, the weak$\equiv$strong solutions to ${\mathrm{VI}}[F,{\cal
X}]$ are, respectively, the minimizers of $f$ on ${\cal X}$ and the saddle
points ($\min$ in $u$, $\max$ in $v$) of $f$ on $U\times V$.

Application to GGLM case. Given a GGLM from Section 2, we associate with an
observation

$\upsilon^{N}=\\{\upsilon_{\tau}=(\zeta_{\tau},H_{\tau}),1\leq\tau\leq
N\\}\in\overbrace{[{{\cal Z}}_{1}\times{{\cal
H}}_{1}]}^{\Upsilon_{1}}\times\ldots\times\overbrace{[{{\cal
Z}}_{N}\times{{\cal H}}_{N}]}^{\Upsilon_{N}},$

a collection of scalar weights $\lambda_{t}>0$, $1\leq t\leq N,$ which should
be deterministic functions of $(\upsilon^{t-1},H_{t})$:
$\lambda_{t}=\Lambda_{t}(\upsilon^{t-1},H_{t})$ and two vector fields. First,
the observable field

$F_{\upsilon^{N}}(z)={1\over
N}\sum_{t=1}^{N}\Lambda_{t}^{-1}(\upsilon^{t-1},H_{t})H_{t}\left[\Phi_{t}\left(H_{t}^{T}z\right)-\zeta_{t}\right]:{\cal
B}\to{\mathfrak{B}},$ (3)

and second, the unobservable field which is a “population” version of
$F_{\upsilon^{N}}(\cdot)$,

$\overline{F}_{\upsilon^{N}}(z)={1\over
N}\sum_{t=1}^{N}\Lambda_{t}^{-1}(\upsilon^{t-1},H_{t})H_{t}\left[\Phi_{t}\left(H_{t}^{T}z\right)-\Phi_{t}\left(H_{t}^{T}\beta\right)\right]:{\cal
B}\to{\mathfrak{B}}.$ (4)

Note that in our situation, both of these vector fields are well-defined and
continuous on ${\cal B}$; they clearly are monotone along with
$\Phi_{t}$.111We have used the elementary observation that if $x\mapsto Ax+b$
is affine mapping which maps a convex domain $X\subset{\mathbf{R}}^{N}$ into
the domain ${\mathrm{Dom}\,}G\subset{\mathbf{R}}^{M}$ of a monotone vector
field $G$, then the vector field $A^{T}G(Ax+b)$ is monotone on $X$. In
addition, ${\cal B}$ is compact.

Consequently, the VI is defined using the observed field (3):

$\hbox{Find\ }z\in{\cal B}:\langle F_{\upsilon^{N}}(x),x-z\rangle\geq
0,\,\,\forall x\in{\cal B}$ $None$

has solutions, and these weak solutions are the same as the strong solutions
to the VI. Our estimate of $\beta$, denoted by $\widehat{\beta}$, will be the
weak$\equiv$strong solution to the variational inequality
${\mathrm{VI}}(F_{\upsilon^{N}},{\cal B})$.

Special case. When the monotone vector mappings
$\Phi_{t}:{\mathrm{Dom}\,}\Phi_{t}\to{\mathfrak{Z}}_{t}$ are the gradient
fields of convex functions
$\mathfrak{F}_{t}(\cdot):{\mathrm{Dom}\,}\Phi_{t}\to{\mathbf{R}}$, the vector
field $F_{\upsilon^{N}}$ is the gradient field of the convex function

$\mathfrak{F}_{\upsilon^{N}}(x)={1\over
N}\sum_{t=1}^{N}\Lambda^{-1}_{t}(\upsilon^{t-1},H_{t})\left[\mathfrak{F}_{t}(H_{t}^{T}x)-x^{T}H_{t}\zeta_{t}\right]$
(5)

so that the solutions to ${\mathrm{VI}}(F_{\upsilon^{N}},{\cal B})$ are
exactly the optimal solutions to the convex minimization problem

$\min_{x\in{\cal B}}\mathfrak{F}_{\upsilon^{N}}(x).$ (6)

In particular, in the case of the identity link
$\Phi_{t}(\eta_{t})\equiv\eta_{t}$, $\eta_{t}\in{\mathfrak{Z}}_{t}$, we have
$\Phi_{t}(\eta_{t})=\nabla\left[{1\over 2}\eta_{t}^{T}\eta_{t}\right]$, thus
$\mathfrak{F}_{t}(x)\propto\frac{1}{2}x^{T}x$,
$\Lambda^{-1}_{t}(\upsilon^{t-1},H_{t})=1$, and it can be verified that the
proposed recovery becomes just the Least Squares:

$\min_{x\in{\cal
B}}{1\over{2N}}\sum_{t=1}^{N}\left\|H_{t}^{T}x-\zeta_{t}\right\|^{2}.$ (7)

Maximum Likelihood estimation. Assuming that the density, taken w.r.t. some
reference measure, of the conditional, $\upsilon^{t-1}$ and $H_{t}$ given,
distribution of $\zeta_{t}$ is a known function of the form
$\psi_{t}(\zeta_{t};H_{t}^{T}\beta)$ as specified in Section 2.2, we could try
to recover $\beta$ by maximizing the likelihood of what we have observed, that
is, by solving the optimization problem

$\max_{x\in{\cal B}}{1\over
N}\sum_{t=1}^{N}\ln(\psi_{t}(\zeta_{t};H_{t}^{T}x)).$

In order for this approach to be computation-friendly, the problem should be
convex. When the link function is derived from the exponential family, the
maximum likelihood problems are usually convex. In such cases, we can use the
gradient of the log-likelihood function as the monotone operator in VI, and it
can be shown that the weights $\Lambda_{t}^{-1}(\upsilon^{t-1},H_{t})$ in (3)
can be related to terms in the gradient of the log-likelihood function.

Remark. For GLM, the maximum likelihood parameter estimation is typically done
by Iterative Weighted Least Squares (IWLS) [4]. One may wonder whether it
relates to the proposed VI-based approach using (3). However, they are not the
same. IWLS is a heuristic iterative procedure motivated by Fisher’s scoring
method (a variant of Newton’s step using the expected value of the Hessian
matrix). IRLS iteratively fits the adjusted response variable with respect to
predictors using the weights that depend on the fitted parameter values from
the previous iterations. The parameter recovery guarantee for IWLS is unclear.
$\Box$

##  4 Performance Guarantees

We now establish results underlying theoretical guarantees for parameter
recovery. We focus on the results for parameters recovered using VI.

###  4.1 Concentration of empirical vector field

An important observation is as follows:

$\beta$ is the root of $\overline{F}_{\upsilon^{N}}(\cdot)$ in ${\cal B}$:
$\overline{F}_{\upsilon^{N}}(\beta)=0$,

and

$F_{\upsilon^{N}}(\beta)=F_{\upsilon^{N}}(\beta)-\overline{F}_{\upsilon^{N}}(\beta)={1\over
N}\sum_{t=1}^{N}\Lambda_{t}^{-1}(\upsilon^{t-1},H_{t})H_{t}\left[\Phi_{t}\left(H_{t}^{T}\beta\right)-\zeta_{t}\right].$

What is under the summation is martingale-difference, since

${\mathbf{E}}_{|\upsilon^{t-1},H_{t}}[F_{\upsilon^{N}}(\beta)-\overline{F}_{\upsilon^{N}}(\beta)]={1\over
N}\Lambda_{t}^{-1}(\upsilon^{t-1},H_{t})H_{t}\cdot{\mathbf{E}}_{|\upsilon^{t-1},H_{t}}\left[\Phi_{t}\left(H_{t}^{T}\beta\right)-\zeta_{t}\right]=0;$

we have used (1). We first derive the following concentration inequality,
which is useful in establishing our results for parameter recovery and online
bounds for GGLM. We will illustrate their usage for spatial-temporal processes
in Section 5.

######  Proposition 4.1 (Concentration of empirical vector field).

Assume that for every $t$, every
$\upsilon^{t-1}\in\Upsilon_{1}\times\ldots\times\Upsilon_{t-1}$ and every
$H_{t}\in{{\cal H}}_{t}$ the conditional, $\upsilon^{t-1},H_{t}$ given,
distribution $P_{|\upsilon^{t-1},H_{t}}$ of $\zeta_{t}$ satisfies

$\ln\left({\mathbf{E}}_{\zeta_{t}\sim
P_{|\upsilon^{t-1},H_{t}}}\left\\{\exp\\{h^{T}\zeta_{t}\\}\right\\}\right)\leq
h^{T}\Phi_{t}\left(H_{t}^{T}\beta\right)+\Psi_{t,\upsilon^{t-1},H_{t}}(\|P_{t}h\|_{1}),\,\,\forall
h\in{\mathfrak{Z}}_{t}$ (8)

where the rate functions
$\Psi_{t,\upsilon^{t-1},H_{t}}(r):{\mathbf{R}}_{+}\to{\mathbf{R}}$ are
continuous nondecreasing functions, and $P_{t}$ are known deterministic
matrices.

Let $E_{t}$, $t=1,2,...,$ be $m\times n_{t}$ matrices which are deterministic
functions of $(\upsilon^{t-1},H_{t})$: $E_{t}=E_{t}(\upsilon^{t-1},H_{t})$,
and let

$G_{\upsilon^{N}}={1\over
N}\sum_{t=1}^{N}E_{t}[\Phi_{t}(H_{t}^{T}\beta)-\zeta_{t}].$

Then

(i) For $\upsilon_{\tau}\in\Upsilon_{\tau}$, $1\leq\tau\leq N$, $1\leq j\leq
m$, let $\Theta_{sj}(\upsilon^{s-1},H_{s})$ be the $\|\cdot\|_{1}$-norm of
$j$-th column in $P_{s}E^{T}_{s}$, $s=1,2,\ldots,N$. Given a finite set
$\Gamma=\\{\alpha_{i}>0,1\leq i\leq K\\}$ on the positive ray, and a tolerance
$\epsilon\in(0,1)$, for every $j\leq m$ the probability of the event

$\bigg{\\{}\upsilon^{N}:|[G_{\upsilon^{N}}]_{j}|>\underbrace{\min\limits_{1\leq
i\leq
K}\bigg{[}\alpha_{i}\ln(2K/\epsilon)+\alpha_{i}\sum_{t=1}^{N}\Psi_{t,\upsilon^{t-1},H_{t}}(\alpha_{i}^{-1}N^{-1}\Theta_{tj}(\upsilon^{t-1},H_{t}))}_{\delta_{j}(\upsilon^{N})}\bigg{]}\bigg{\\}}$
(9)

is at most $\epsilon$.

(ii) Assume that $\Psi_{t,\upsilon^{t-1},H_{t}}(\cdot)\leq\Psi(\cdot)$ for all
$t$, where $\Psi(\cdot):{\mathbf{R}}_{+}\to{\mathbf{R}}_{+}$ is a continuous
nondecreasing “worst case” rate function, and that for some $\Theta<\infty$,

$\max_{i}\|\mathrm{Col}_{i}[P_{t}E_{t}^{T}]\|_{1}\leq\Theta,\mbox{ whenever
}\upsilon_{s}\in\Upsilon_{s},1\leq s\leq t,\forall t.$

(from now on, $\mathrm{Col}_{j}[A]$ stand for $j$-to column, and
$\mathrm{Row}_{i}[A]$ – for the transpose of $i$-th row in matrix $A$). Then
for every $\delta\geq 0$ one has

${\mathrm{Prob}}\\{\|G_{\upsilon^{N}}\|_{\infty}>\delta\\}\leq
2m\exp\\{-N\Psi_{*}(\delta/\Theta)\\},$ (10)

where for $s\geq 0$

$\Psi_{*}(s)=\sup_{\gamma\geq 0}\left[\gamma s-\Psi(\gamma)\right]\geq 0,$
(11)

is the Fenchel-Legendre transformation of $\Psi(\cdot)$.

Remark. Assume that we are in the situation of the item (ii) and are given
$\varepsilon\in(0,1)$, and let us look at the upper bound of the
$(1-\varepsilon)$-quantile of $\|G_{\upsilon^{N}}\|_{\infty}$ we could get
from (10). In order for (10) to imply that
${\mathrm{Prob}}\\{\|G_{\upsilon^{N}}\|_{\infty}>\delta\\}\leq\varepsilon$ we
need to ensure that

$\sup_{\gamma>0}\left[\gamma[\delta/\Theta]-\Psi(\gamma)\right]\geq
N^{-1}\ln(2m/\varepsilon),$

or, which is the same (set $\gamma=\Theta/(N\alpha)$), that

$\sup_{\alpha>0}\left[\delta-\alpha
N\Psi(\Theta/(N\alpha))-\alpha\ln(2m/\varepsilon)\right]\geq 0.$

The smallest $\delta=\widehat{\delta}$ satisfying this requirement is

$\widehat{\delta}=\inf_{\alpha>0}\left[\alpha\ln\left(2m/\varepsilon\right)+\alpha
N\Psi(N^{-1}\Theta/\alpha)\right].$

On the other hand, we are in the situation when
$\Theta_{tj}(\upsilon^{t-1},H_{t})\leq\Theta$ for all $t$, $j$ and all
$\upsilon^{t-1}\in\Upsilon^{t-1}$, $H_{t}\in{{\cal H}}_{t}$. It follows that
given $\delta>\widehat{\delta}$, selecting $\bar{\alpha}>0$ in such a way that

$\left[\bar{\alpha}\ln\left(2m/\varepsilon\right)+\bar{\alpha}N\Psi(N^{-1}\Theta/\bar{\alpha})\right]\leq\delta$

and specifying $\Gamma=\\{\bar{\alpha}\\}$, $\epsilon=\varepsilon/m$, the
right-hand side in the inequality in (9) would, for every $j\leq m$ be
$\leq\delta$, implying, via the union bound, that the probability of the event

$\|G_{\upsilon^{N}}\|_{\infty}>\delta$

is $\leq\varepsilon$. Thus, in the case of (ii), the result of (i) with
properly selected singleton $\Gamma$ would result in the upper bound on the
$(1-\varepsilon)$-quantile of $\|G_{\upsilon^{N}}\|_{\infty}$ which can be
made arbitrarily close to the bound yielded by (ii). ∎

The rationale beyond (i) is that it results in online bounds
$\delta_{j}(\upsilon^{N})$ on the $(1-\epsilon)$-quantile of
$\|G_{\upsilon^{N}}\|_{\infty}$, meaning that when
$\Psi_{t,\upsilon^{t-1},H_{t}}$ is observable, the bound is an observable
function of $\upsilon^{N}$ such that the probability for
$|[G_{\upsilon^{N}}]_{j}|$ to exceed this bound is at most $\epsilon$. This
online bound adjusts itself to the observed rate functions
$\Psi_{t,\upsilon^{t-1},H_{t}}(\cdot)$ and observed values of
$\Theta_{tj}(\upsilon^{t-1},H_{t})$, in contrast to the deterministic bound
$\widehat{\delta}$ oriented on the worst case rate function and a priori upper
bound $\Theta$ on $\Theta_{tj}(\upsilon^{t-1},H_{t})$’s. Note that the price
for cardinality $K$ of $\Gamma$ (the larger $K$, the stronger “adaptive
abilities” of the online bound $\delta_{j}(\cdot)$) is low – $K$ is under the
log.

Proof of Proposition 4.1. Let us set

$\begin{array}[]{rcl}S_{t}(\upsilon^{t})&=&\sum_{s=1}^{t}{1\over
N}E_{s}\left[\Phi_{s}\left(H_{s}^{T}\beta\right)-\zeta_{s}\right],\\\
R_{t}(\upsilon^{t};g)&=&\sum_{s=1}^{t}\Psi_{s,\upsilon^{s-1},H_{s}}(N^{-1}\|P_{s}E_{s}^{T}g\|_{1})\quad[g\in{\mathbf{R}}^{m}],\\\
\end{array}$

and let ${\mathbf{E}}_{t}$ stand for expectation over $\upsilon^{t}$, and
${\mathbf{E}}^{+}_{t}$ stand for expectation over $(\upsilon^{t},H_{t+1})$.

(i): Let us fix $j$, and let $g=\chi e_{j}$, where $e_{j}$ is $j$-th standard
basic orth in ${\mathbf{R}}^{m}$ and $\chi=\pm 1$, and $\gamma\geq 0$, and let

$\begin{array}[]{rcl}\Sigma_{t}(\upsilon^{t})&=&\sum_{s=1}^{t}\bigg{[}{1\over
N}\gamma
g^{T}E_{s}\left[\Phi_{s}\left(H_{s}^{T}\beta\right)-\zeta_{s}\right]-\Psi_{s,\upsilon^{s-1},H_{s}}(\|\gamma
N^{-1}P_{s}E_{s}^{T}g\|_{1})\bigg{]}\\\ &=&\gamma
g^{T}S_{t}(\upsilon^{t})-R_{t}(\upsilon^{t};\gamma g).\end{array}$ (12)

We have

$\begin{array}[]{l}{\mathbf{E}}_{t+1}\left\\{\exp\\{\Sigma_{t+1}(\upsilon^{t+1})\\}\right\\}\\\
={\mathbf{E}}_{t}^{+}\bigg{\\{}{\mathbf{E}}_{\zeta_{t+1}\sim
P_{|\upsilon^{t},H_{t+1}}}\bigg{\\{}\exp\bigg{\\{}\Sigma_{t}(\upsilon^{t})-\Psi_{t+1,\upsilon^{t},H_{t+1}}(\|\gamma
N^{-1}P_{t+1}E_{t+1}^{T}g\|_{1})\\\ \lx@intercol\hfil+\gamma
N^{-1}[E_{t+1}^{T}g]^{T}[{\Phi}_{t+1}\left(H_{t+1}^{T}\beta\right)-\zeta_{t+1}]\bigg{\\}}\bigg{\\}}\bigg{\\}}\lx@intercol\\\
={\mathbf{E}}_{t}^{+}\bigg{\\{}\exp\bigg{\\{}\Sigma_{t}(\upsilon^{t})-\Psi_{t+1,\upsilon^{t},H_{t+1}}(\|\gamma
N^{-1}P_{t+1}E_{t+1}^{T}g\|_{1})\bigg{\\}}\\\
\lx@intercol\hfil\times\underbrace{{\mathbf{E}}_{\zeta_{t+1}\sim
P_{|\upsilon^{t},H_{t+1}}}\bigg{\\{}\exp\big{\\{}\gamma
N^{-1}[E_{t+1}^{T}g]^{T}[{\Phi}_{t+1}\left(H_{t+1}^{T}\beta\right)-\zeta_{t+1}]\big{\\}}\bigg{\\}}}_{\leq\exp\\{\Psi_{t+1,\upsilon^{t},H_{t+1}}(\|\gamma
N^{-1}P_{t+1}E_{t+1}^{T}g\|_{1})\\}\hbox{\tiny\ by
(\ref{eq2})}}\bigg{\\}}\lx@intercol\\\
\leq{\mathbf{E}}_{t}^{+}\left\\{\exp\\{\Sigma_{t}(\upsilon^{t})\\}\right\\}={\mathbf{E}}_{t}\left\\{\exp\\{\Sigma_{t}(\upsilon^{t})\\}\right\\}.\\\
\end{array}$

This reasoning works for all $t\geq 0$, with the empty sum
$\Sigma_{0}(\upsilon^{0})$, as usual, identically equal to 0. Consequently,
${\mathbf{E}}_{\upsilon^{N}}\left\\{\exp\\{\Sigma_{N}(\upsilon^{N})\\}\right\\}\leq
1,$ or, which is the same (note that $G_{\upsilon^{N}}=S_{N}(\upsilon^{N})$):

${\mathbf{E}}_{\upsilon^{N}}\left\\{\exp\\{\gamma
g^{T}G_{\upsilon^{N}}-R_{N}(\upsilon^{N};\gamma g)\\}\right\\}\leq 1,$

whence for every $\delta\geq 0$ one has by Markov inequality that

${\mathrm{Prob}}\left\\{\gamma g^{T}G_{\upsilon^{N}}-R_{N}(\upsilon^{N};\gamma
g)>\gamma\delta\right\\}\leq\exp\\{-\gamma\delta\\}.$

Setting $\gamma=1/\alpha$ and recalling what $R_{N}$ is, we get

${\mathrm{Prob}}\left\\{g^{T}G_{\upsilon^{N}}>\alpha\sum_{t=1}^{N}\Psi_{t,\upsilon^{t-1},H_{t}}(\|\alpha^{-1}N^{-1}P_{t}E_{t-1}^{T}g\|_{1})+\delta\right\\}\leq\exp\\{-\delta/\alpha\\},$

whence, due to monotonicity of $\Psi_{t,\upsilon^{t-1},H_{t}}(\cdot)$ and the
fact that $\|P_{t}E_{t}^{T}g\|_{1}\leq\Theta_{tj}(\upsilon^{t-1},H_{t})$, we
have

${\mathrm{Prob}}\left\\{g^{T}G_{\upsilon^{N}}>\alpha\sum_{t=1}^{N}\Psi_{t,\upsilon^{t-1},H_{t}}(\alpha^{-1}N^{-1}\Theta_{tj}(\upsilon^{t-1},H_{t}))+\delta\right\\}\leq\exp\\{-\delta/\alpha\\}.$

Using the union bound, we get

${\mathrm{Prob}}\left\\{|[G_{\upsilon^{N}}]_{j}|>\delta+\alpha\sum_{t=1}^{N}\Psi_{t,\upsilon^{t-1},H_{t}}(\alpha^{-1}N^{-1}\Theta_{tj}(\upsilon^{t-1},H_{t}))\right\\}\leq
2\exp\\{-\delta/\alpha\\}.$ (13)

For $1\leq i\leq K$, let us set $\delta_{i}=\alpha_{i}\ln(2K/\epsilon)$, so
that by (13) the probability of the event

${\cal
E}_{i}=\\{\upsilon^{N}:|[G_{\upsilon^{N}}]_{j}|>\alpha_{i}\ln(2K/\epsilon)+\alpha_{i}\sum_{t=1}^{N}\Psi_{t,\upsilon^{t-1},H_{t}}(\alpha_{i}^{-1}N^{-1}\Theta_{tj}(\upsilon^{t-1},H_{t}))\\}$

is at most $\epsilon/K$. Consequently, the probability of the event
$\cup_{i\leq K}{\cal E}_{i}$, which is exactly the event (9), is at most
$\epsilon$, as claimed in (i).

(ii): From the definition of $\Theta$ it follows that
$\|P_{t}E_{t}^{T}g\|_{1}\leq\Theta\|g\|_{1}$ for every $g\in{\mathbf{R}}^{m}$
and every $t$ and $\upsilon^{t}\in\Upsilon^{t}$. Now let
$g\in{\mathbf{R}}^{m}$ with $\|g\|_{1}=1$, and let $\gamma\geq 0$. We have

$\begin{array}[]{l}{\mathbf{E}}_{t+1}\left\\{\exp\\{\gamma
g^{T}S_{t+1}(\upsilon^{t+1})\\}\right\\}\\\
={\mathbf{E}}_{t}^{+}\left\\{\exp\\{\gamma
g^{T}S_{t}(\upsilon^{t})\\}{\mathbf{E}}_{\zeta_{t+1}\sim
P_{|\upsilon^{t},H_{t+1}}}\left\\{\exp\\{[\gamma
N^{-1}[E_{t+1}^{T}g]^{T}[{\Phi}_{t+1}\left(H_{t+1}^{T}\beta\right)-\zeta_{t+1}]\\}\right\\}\right\\}\\\
\leq{\mathbf{E}}_{t}^{+}\left\\{\exp\\{\gamma
g^{T}S_{t}(\upsilon^{t})+\Psi(\|\gamma
N^{-1}P_{t+1}E_{t+1}^{T}g\|_{1})\\}\right\\}\hbox{\ [by (\ref{eq2}) and due to
$\Psi_{t,\upsilon^{t-1},H_{t}}(\cdot)\leq\Psi(\cdot)$]}\\\
\leq{\mathbf{E}}_{t}^{+}\left\\{\exp\\{\gamma
g^{T}S_{t}(\upsilon^{t})\\}\right\\}\exp\\{\Psi(\gamma\Theta/N)\\},\\\
\lx@intercol\hfil\hbox{\ [since $\Psi$ is nondecreasing and
$\|P_{t+1}E_{t+1}^{T}g\|_{1}\leq\Theta\|g\|_{1}=\Theta$]}\lx@intercol\\\
={\mathbf{E}}_{t}\left\\{\exp\\{\gamma
g^{T}S_{t}(\upsilon^{t})\\}\right\\}\exp\\{\Psi(\gamma\Theta/N)\\},\end{array}$

whence

$\ln\left({\mathbf{E}}_{\upsilon^{N}}\left\\{\gamma
g^{T}S_{N}(\upsilon^{N})\right\\}\right)\leq N\Psi(\gamma\Theta/N).$

Therefore, for every $\delta>0$ one has, using Markov inequality

${\mathrm{Prob}}\\{g^{T}S_{N}(\upsilon^{N})>\delta\\}\leq\exp\\{N\Psi(\gamma\Theta/N)-\gamma\delta\\},$

whence finally

$\displaystyle{\mathrm{Prob}}\\{g^{T}S_{N}(\upsilon^{N})>\delta\\}$
$\displaystyle\leq$ $\displaystyle\exp\\{-N\Psi_{*}(\delta/\Theta)\\},$ (14)
$\displaystyle\Psi_{*}(s)$ $\displaystyle=$ $\displaystyle\sup_{\alpha\geq
0}\left[\alpha s-\Psi(\alpha)\right].$ (15)

Taking into account that $G_{\upsilon^{N}}=S_{N}(\upsilon^{N})$, selecting $g$
as $\pm$ basic orths in ${\mathbf{R}}^{m}$ and using the union bound, we
arrive at (10). (ii) is proved. ∎

###  4.2 Quality of parameter recovery

Now we present a model parameter recovery performance bound, as a consequence
of the results in Section 4.1. Assume that we are under the premise of
Proposition 4.1.ii and that the conditional by the past distributions of
$\zeta_{t}$ are “light tail” ones, so that the worst case rate function
$\Psi(\cdot)$ is such that $\Psi(s)\leq Cs$ for all $s\in[0,c]$, with properly
selected $c>0$ and $C<\infty$. Specifying

$m=\kappa,\,\,E_{t}=\Lambda_{t}^{-1}(\upsilon^{t-1},H_{t})H_{t},$ (16)

we get

$G_{\upsilon^{N}}=F_{\upsilon^{N}}(\beta).$

In this case, under the premise of Proposition 4.1.ii, the typical
$\|\cdot\|_{\infty}$ norm of $F_{\upsilon^{N}}(\beta)$ for large $N$ is of
order of $\Theta\sqrt{\ln(\kappa)N}$. Thus, when there are reasons to believe
that $F_{\upsilon^{N}}$ is, typically, strongly monotone on ${\cal B}$ with
parameter of strong monotonicity not deteriorating as $N$ grows, we may hope
that the solution to ${\mathrm{VI}}(F_{\upsilon^{N}},{\cal B})$ will be at the
distance $O(N^{-1/2})$ from $\beta$.

To formulate the precise statement, let us define the $p$-modulus of strong
monotonicity $\theta_{p}(\upsilon^{N})$ of $F_{\upsilon^{N}}(\cdot)$ on ${\cal
B}$, $p\in[1,\infty]$, as

$\theta_{p}(\upsilon^{N})=\max\left\\{\theta:\langle
F_{\upsilon^{N}}(x)-F_{\upsilon^{N}}(y),x-y\rangle\geq\theta\|x-y\|_{p}^{2},\,\,\forall
x,y\in{\cal B}\right\\}.$

Note that the strict positivity of $\theta_{p}(\upsilon^{N})$ is independent
of the value of $p$.

######  Proposition 4.2 (Parameter recovery guarantee).

Let $\upsilon^{N}$ be such that $\theta_{p}(\upsilon^{N})>0$. Then the
weak$\equiv$strong solution $\widehat{\beta}(\upsilon^{N})$ to the variational
inequality ${\mathrm{VI}}(F_{\upsilon^{N}},{\cal B})$ is unique, and

$\|\widehat{\beta}(\upsilon^{N})-\beta\|_{p}\leq\|F_{\upsilon^{N}}(\beta)\|_{\infty}/\sqrt{\theta_{p}(\upsilon^{N})\theta_{1}(\upsilon^{N})}.$
(17)

Proof. We know that in the situation of Section 2.1 weak solution
$\widehat{\beta}=\widehat{\beta}(\upsilon^{N})$ to
${\mathrm{VI}}(F_{\upsilon^{N}},{\cal B})$ exists and is a strong solution;
besides this, it is well known that under the premise of Proposition, this
solution is unique. Since $\widehat{\beta}$ is a strong solution, we have
$\langle F_{\upsilon^{N}}(\widehat{\beta}),x-\widehat{\beta}\rangle\geq 0$ for
all $x\in{\cal B}$, and in particular $\langle
F_{\upsilon^{N}}(\widehat{\beta}),\beta-\widehat{\beta}\rangle\geq 0$. As a
result,

$\|F_{\upsilon^{N}}(\beta)\|_{\infty}\|\beta-\widehat{\beta}\|_{1}\geq\langle
F_{\upsilon^{N}}(\beta),\beta-\widehat{\beta}\rangle=\underbrace{\langle
F_{\upsilon^{N}}(\beta)-F_{\upsilon^{N}}(\widehat{\beta}),\beta-\widehat{\beta}\rangle}_{\geq\sqrt{\theta_{p}(\upsilon^{N})\theta_{1}(\upsilon^{N})}\|\beta-\widehat{\beta}\|_{p}\|\beta-\widehat{\beta}\|_{1}}+\underbrace{\langle
F_{\upsilon^{N}}(\widehat{\beta}),\beta-\widehat{\beta}\rangle}_{\geq 0},$

and (17) follows. ∎

To extract from Propositions 4.1 information on performance guarantees of our
estimate, we need to understand how to lower-bound the modulus of strong
monotonicity. The related theoretical analysis would require heavy assumptions
on inter-dependence of the subsequent “regressors” $H_{t}$, $t=1,2,\ldots$.
This is an issue we do not want to touch here. Instead, let us focus on how to
lower-bound this modulus “online.” An immediate observation is that if
$\Phi_{t}$ is strongly monotone, with positive moduli $\gamma_{t}$, on their
domains:

$\langle\Phi_{t}(u)-\Phi_{t}(v),u-v\rangle\geq\gamma_{t}\langle
u-v,u-v\rangle,\,\forall u,v\in{\mathrm{Dom}\,}\Phi_{t},$

then

$\begin{array}[]{ll}&\langle
F_{\upsilon^{N}}(x)-F_{\upsilon^{N}}(y),x-y\rangle\\\ =&{1\over
N}\sum_{t=1}^{N}\Lambda_{t}^{-1}(\upsilon^{t-1},H_{t})\langle
H_{t}\Phi_{t}(H_{t}^{T}x)-H_{t}\Phi_{t}(H_{t}^{T}y),x-y\rangle\\\ \geq&{1\over
N}\sum_{t=1}^{N}\gamma_{t}\Lambda_{t}^{-1}(\upsilon^{t-1},H_{t})\langle
H_{t}^{T}(x-y),H_{t}^{T}(x-y)\rangle\\\ =&(x-y)^{T}\underbrace{\left[{1\over
N}{\sum}_{t=1}^{N}\gamma_{t}\Lambda_{t}^{-1}(\upsilon^{t-1},H_{t})H_{t}H_{t}^{T}\right]}_{\Gamma(\upsilon^{N})}(x-y),\\\
\end{array}$

whence

$\theta_{2}(\upsilon^{N})\geq\lambda_{\min}(\Gamma(\upsilon^{N})),$ (18)

where $\lambda_{\min}(Q)$ is the minimal eigenvalue of a symmetric matrix $Q$.
We see that $\theta_{2}(\upsilon^{N})$ admits an “online observable” lower
bound. We have also

$\theta_{p}(\upsilon^{N})\geq\min_{h:\|h\|_{p}=1}h^{T}\Gamma(\upsilon^{N})h=\left[\max_{g:\|g\|_{p_{*}}\leq
1}g^{T}\Gamma^{-1}(\upsilon^{N})g\right]^{-1}$ $None$

When $1\leq p\leq 2$, the right-hand side in this relation can be efficiently
lower-bounded, the lower bound being tight within absolute constant [6] (equal
to $\pi/2$ when $p=1$).

####  4.2.1 The case of affine $\Phi_{t}$’s

Assume that the links $\Phi_{t}(\cdot)$ are affine; in this situation, the
vector field $F_{\upsilon^{N}}(\cdot)$ is affine as well. Assume also that the
functions $\Psi_{t,\upsilon^{t-1},H_{t}}$ are observable. In this case
observations $\upsilon^{N}$ imply $(1-\kappa\epsilon)$-reliable upper bound
$\delta(\upsilon^{N})$ on $\|F_{\upsilon^{N}}(\beta)\|_{\infty}$,
specifically, the bound $\delta(\upsilon^{N}):=\max\limits_{1\leq
j\leq\kappa}\delta_{j}(\upsilon^{N})$ with $\delta_{j}(\upsilon^{N})$ given by
(9) as applied with

$\widetilde{E}_{t}=\Lambda^{-1}_{t}(\upsilon^{t-1},H_{t})H_{t}.$

(Recall that this choice of $\widetilde{E}_{t}$ results in
$G_{\upsilon^{N}}=F_{\upsilon^{N}}(\beta)$). As a result, the observable under
our assumptions convex set

$\Delta(\upsilon^{N})=\\{x\in{\cal
B}:\|F_{\upsilon^{N}}(x)\|_{\infty}\leq\delta(\upsilon^{N})\\}$ (19)

is “$(1-\kappa\epsilon)$-confidence interval” for $\beta$ – this set contains
$\beta$ with probability $\geq 1-\kappa\epsilon$. In particular, given a norm
$\|\cdot\|$ on ${\mathbf{R}}^{\kappa}$, the quantity

$d(\upsilon^{N})=\max\limits_{x\in\Delta(\upsilon^{N})}\|x-\widehat{\beta}\|$

is an $(1-\kappa\epsilon)$-reliable upper bound on the $\|\cdot\|$-error of
our estimate $\widehat{\beta}$. Whether this online error bound can or cannot
be efficiently computed, depends on the norm $\|\cdot\|$; e.g., for the
uniform norm $\|x\|_{\infty}=\max_{i}|x_{i}|$, computing the bound is easy -
it reduces to $2\kappa$ maximizations of linear forms over
$\Delta(\upsilon^{N})$.

Note that under our present assumptions every policy ${\cal R}$ for generating
row vectors $\widetilde{E}_{t}\in{\mathbf{R}}^{1\times n_{t}}$ as
deterministic functions of $\upsilon^{t-1},H_{t}$, $t=1,...,N$, induces an
observable, given $\upsilon^{N}$, affine real-valued function

$G^{{\cal R}}_{\upsilon^{N}}(x)={1\over
N}\sum\limits_{t=1}^{N}\widetilde{E}_{t}\left[\Phi_{t}(H_{t}^{T}x)-\zeta_{t}\right]:{\cal
B}\to{\mathbf{R}}$

along with observable quantity $\delta^{{\cal R}}_{\epsilon}(\upsilon^{N})$,
given by (9), such that

${\mathrm{Prob}}\\{|G^{{\cal R}}_{\upsilon^{N}}(\beta)|\geq\delta^{{\cal
R}}(\upsilon^{N})\\}\leq\epsilon.$

Thus, the pair of observable linear inequalities $|G^{{\cal
R}}_{\upsilon^{N}}(\beta)|\leq\delta^{{\cal R}}(\upsilon^{N})$ on $\beta$ do
hold true with probability $\geq 1-\epsilon$, and we can add these
inequalities to the description of $\Delta(\upsilon^{N}))$ as given by (19),
thus refining the confidence interval, and we can implement several
refinements of this type stemming from several policies ${\cal R}$. For an
instructive application example, see Section 5.3.

##  5 Illustrations: Spatio-Temporal Processes

Assume that there are $L$ spatial locations; at discrete time $t$, the state
of location $k\in\\{1,\ldots,L\\}$ is described by:

realizations $\omega_{t,k}$ of random variable taking values in a set
$Z\subset{\mathbf{R}}^{\mu}$,

and $\mu$ is the number of possible states when the observation is non-zero.
Our observation at time $t$ is the block vector

$\zeta_{t}=[\omega_{t,1};\ldots;\omega_{t,L}]\in{{\cal Z}}=Z\times...\times
Z\subset{\mathfrak{Z}}={\mathbf{R}}^{\mu L}.$

We assume that the conditional (by what happened prior to time $t$)
expectation of $\zeta_{t}$ depends solely on the collection

$(\zeta_{t-d},\zeta_{t-d+1},\ldots,\zeta_{t-1})=\\{\omega_{s,\ell}:t-d\leq
s\leq t-1,1\leq\ell\leq L\\},$

where $d\geq 1$ is the “memory depth” of our process. Specifically, setting

$\zeta_{\tau}^{t}=(\zeta_{\tau},\zeta_{\tau+1},\ldots,\zeta_{t})=\\{\omega_{s,\ell}\in{\mathbf{R}}^{\mu}:\tau\leq
s\leq t,1\leq\ell\leq L\\},\zeta^{t}=\zeta_{-d+1}^{t}.$

Assume that our observations start at time $-d+1$ and that for every $t\geq 1$
the conditional, given $\zeta^{t-1}$, expectation
${\mathbf{E}}_{|\zeta^{t-1}}$ of $\zeta_{t}$ is

${\mathbf{E}}_{|\zeta^{t-1}}\\{\zeta_{t}\\}=\Phi(\eta^{T}(\zeta_{t-d}^{t-1})\beta),$
(20)

where

  * •

$\beta$ is the vector of parameters of the process. This vector is obtained by
writing down in once for every prescribed order the entries of the elements of
the collection

$\\{\beta^{0}\in{\mathbf{R}}^{\mu L},\beta^{s}\in{\mathbf{R}}^{\mu L\times\mu
L},1\leq s\leq d\\}$

of “actual parameters,” and

$\eta^{T}(\zeta_{t-d}^{t-1})\beta=\beta^{0}+\sum_{s=1}^{b}\beta^{s}\zeta_{t-1}\in{\mathbf{R}}^{\mu
L}.$

In other words, denoting by

$\kappa=\mu L+d\mu^{2}L^{2}$

the dimension of $\beta$, $\eta(\gamma)$ is matrix-valued function, taking
values in ${\mathbf{R}}^{\kappa\times\mu L}$, of $L\times d$ array
$\xi=\\{\xi_{s\ell},1\leq s\leq d,1\leq\ell\leq L\\}$ of vectors
$\xi_{s\ell}\in{\mathbf{R}}^{\mu}$, and this function is uniquely defined by
the requirement that

$\forall
x\in{\mathfrak{B}}:\eta^{T}(\xi)x=x^{0}+\sum_{s=1}^{d}x^{s}[\xi_{t-s,1};\ldots;\xi_{t-s,L}],$

where ${\mathfrak{B}}={\mathbf{R}}^{\kappa}$ is the space of column vectors
$x$ obtained by arranging in a prescribed order the entries in a collection
$\\{x^{0}\in{\mathbf{R}}^{\mu L},x^{s}\in{\mathbf{R}}^{\mu L\times\mu L},1\leq
s\leq d\\}$.

We always assume that $\beta$ resides in a known in advance convex compact
subset ${\cal B}$ of ${\mathfrak{B}}$;

  * •

$\Phi(\cdot):{\mathrm{Dom}\,}\Phi\to{\mathfrak{Z}}$ is the link function – a
continuous and monotone vector field on a closed convex domain
${\mathrm{Dom}\,}\Phi\subset{\mathfrak{Z}}$. We always assume that whenever
$\xi\in{{\cal Z}}\times\ldots\times{{\cal Z}}$ and $x\in{\cal B}$, we have
$\eta^{T}(\xi)x\in{\mathrm{Dom}\,}\Phi.$

A spatio-temporal process we just have described can be modeled by stationary
GGLM, where matrix component $H_{t}$ of observations at time $t$ is known
deterministic function of $\zeta_{t-d},\zeta_{t-d+1},\ldots,\zeta_{t-1}$:

$H_{t}=\eta(\zeta_{t-d}^{t-1})\in{\mathbf{R}}^{\kappa\times\mu L},$

and all other components of the model are readily given by the above
description. Consequently, we can apply the machinery from Section 3 to
recover the parameters $\beta$ from observations $\zeta_{t}$ on time horizon
$-d+1\leq t\leq N$.

We are especially interested in the case of categorial and Poisson spatio-
temporal processes.

###  5.1 Categorial spatio-temporal process

The spatio-temporal categorial process with $\mu+1$ states models the
situation when at time $t$ the $k$-th location can be in one of $\mu+1$ states
— ground state 0 and nontrivial states $1,\ldots,\mu$. We specify the model as
follows:

  * •

encode these states by vectors from ${\mathbf{R}}^{\mu}$, with the ground
state encoded by the origin and nontrivial state $i$ encoded by $i$-th basic
orth in ${\mathbf{R}}^{\mu}$. As a result, $Z$ becomes the finite set – the
vertices of the simplex $\\{z\in{\mathbf{R}}^{\mu}_{+}:\sum_{i}z_{i}\leq
1\\}$, ${{\cal Z}}$ becomes the set of Boolean block-vectors with $L$ blocks
of dimension $\mu$ each, and at most one nonzero entry in every block. A
special case is when $\mu=1$, then the model becomes the Bernoulli process
where the observations are binary vectors.

  * •

assume that the states $\omega_{tk}$ are conditionally on $\zeta^{t-1}$
independent across $k$.

  * •

the values of $\Phi$ are block vectors with $L$ blocks of dimension $\mu$
each; denoting by $\Phi_{ik}(z)$ the $i$-th entry of $k$-th block in
$\Phi(z)$.

In order to meet (20), the quantity
$\Phi_{ik}(\eta^{T}(\zeta_{t-d}^{t-1})\beta)$ should be the conditional on
$\zeta^{t-1}$ probability for $k$-th location at time $t$ to be at active
state $i$. Thus, it should be nonnegative, and the sum of these quantities
should be $\leq 1$. Consequently, when speaking about categorial spatio-
temporal processes, we always assume that we are given a closed convex set
$\overline{Z}\subset{\mathrm{Dom}\,}\Phi$ such that

$\forall z\in\overline{Z}:\Phi(z)\geq 0\ \&\sum_{i=1}^{\mu}\Phi_{ik}(z)\leq
1,\,1\leq k\leq L,$

and that

$\eta^{T}(\zeta_{1}^{d})x\in\overline{Z}\,\,\forall x\in{\cal
B},\zeta_{1}^{d}\in{{\cal Z}}\times\cdots\times{{\cal Z}}.$

Note that the conditional on $\zeta^{t-1}$ probability for $k$-th location at
time $t$ to be in ground state is

$1-\sum_{i=1}^{\mu}\Phi_{ik}(\eta^{T}(\zeta_{t-d}^{t-1})\beta).$

Categorial spatio-temporal process of the outlined type and VI-based
techniques for parameter recovery in these processes were considered in [2].

Observe that in the resulting stationary GGLM for every $t\geq 1$ and
$h\in{\mathfrak{Z}}$ setting

$\nu=\nu_{t}(\upsilon^{t-1},H_{t})=\Phi(H_{t}^{T}\beta)$, we get $\nu\geq
0,\|\nu\|_{\infty}\leq 1,$

and

$\begin{array}[]{rcl}{\mathbf{E}}_{|\upsilon^{t-1},H_{t}}\left\\{\exp\\{h^{T}\zeta_{t}\\}\right\\}&=&\prod_{k=1}^{L}\left(1+\sum_{i=1}^{\mu}\nu_{ik}[{\rm
e}^{h_{ik}}-1]\right),\\\ \lx@intercol\Rightarrow\hfil\lx@intercol&&\\\
\ln\left({\mathbf{E}}_{|\upsilon^{t-1},H_{t}}\left\\{\exp\\{h^{T}\zeta_{t}\\}\right\\}\right)&=&h^{T}\nu+\left[\sum_{k=1}^{L}\ln\left(1+\sum_{i=1}^{\mu}\nu_{ik}[{\rm
e}^{h_{ik}}-1]\right)-h^{T}\nu\right].\end{array}$ (21)

The function

$f_{\nu}(g)=\sum_{k=1}^{L}\ln\left(1+\sum_{i=1}^{\mu}\nu_{ik}[{\rm
e}^{g_{ik}}-1]\right)-g^{T}\nu$

of $g\in{\mathfrak{Z}}$ is convex, and therefore its maximum over the set
$\\{g:\|g\|_{1}\leq r:=\|h\|_{1}\\}$ is achieved at a vertex, where all but
one entry of $g$ are zero, and the remaining entry is $\pm r$.222The unit
$\ell_{1}$ ball is also known as the cross-polytope; the vertices of a cross-
polytope can be the canonical basis, i.e., all the vectors as permutations of
$(\pm 1,0,0,\ldots,0)$. As a result, setting

$\begin{array}[]{l}\alpha_{t}(\upsilon^{t-1},H_{t})=\max\limits_{x\in{\cal
B}}\|\Phi(H_{t}^{T}x)\|_{\infty}\in[\|\nu_{t}\|_{\infty},1],\\\
\Psi^{\mathrm{\tiny
cat}}_{t,\upsilon^{t-1},H_{t}}(s):=\max\limits_{\alpha}\left\\{\max\left[\ln\left(1+\alpha[{\rm
e}^{s}-1]\right)-\alpha s,\ln\left(1+\alpha[{\rm e}^{-s}-1]\right)+\alpha
s\right]:0\leq\alpha\leq\alpha_{t}(\upsilon^{t-1},H_{t})\right\\}\\\
\leq\Psi^{\mathrm{\tiny
cat}}(s):=\max\limits_{\alpha}\left\\{\max\left[\ln\left(1+\alpha[{\rm
e}^{s}-1]\right)-\alpha s,\ln\left(1+\alpha[{\rm e}^{-s}-1]\right)+\alpha
s\right]:0\leq\alpha\leq 1\right\\}\\\ \end{array}$

we get continuous convex and even function of $s\in{\mathbf{R}}$, with
$\Psi_{t,\upsilon^{t-1},H_{t}}(\cdot)$ observable at time $t$, such that

$\forall(h\in{\mathfrak{Z}},t\geq 1,\upsilon_{s}\in\Upsilon_{s},s\leq
t):{\mathbf{E}}_{|\upsilon^{t-1},H_{t}}\left\\{\exp\\{h^{T}\zeta_{t}\\}\right\\}\leq
h^{T}\Phi(H_{t}^{T}\beta)+\Psi{}{}^{\mathrm{\tiny cat}}(\|h\|_{1}),$

that is, relation (8) holds true with
$\Psi_{t,\upsilon^{t-1},H_{t}}=\Psi^{\mathrm{\tiny
cat}}_{t,\upsilon^{t-1},H_{t}}$ and the identity matrix in the role of
$P_{t}$.

Observe also that, as it is immediately seen, whenever $\zeta_{s}$, $1\leq
s\leq d$, are Boolean vectors (which definitely is the case when
$\zeta_{s}\in{{\cal Z}}$), $s\leq d$, the matrix
$\eta^{T}({}{\color[rgb]{1,0,0}\zeta_{1}^{d}})$ is Boolean, and every column
of it has at most one nonzero entry, that is, the columns in all realizations
of $H_{t}$ are of $\ell_{1}$-norm not exceeding 1. As a result, when applying
Proposition 4.1 to a spatio-temporal categorial process, we can use the $\mu
L\times\mu L$ identity matrix in the role of $P_{t}$’s, and set
$\Psi_{t,\upsilon^{t-1},H_{t}}\equiv\Psi^{\mathrm{\tiny
cat}}_{t,\upsilon^{t-1},H_{t}}$. Further, when specifying $E_{t}$ according to
(16), resulting in $G_{\upsilon^{N}}=F_{\upsilon^{N}}$, we can set
$\Theta_{tj}(\upsilon^{t-1},H_{t})=\Lambda_{t}^{-1}(\upsilon^{t-1},H_{t})$ in
(i); this choice works also whenever the $\|\cdot\|_{1}$-norms of all columns
in $E_{t}$ are bounded by 1. When, in addition, $\Lambda_{t}(\cdot)\equiv 1$,
the premise of item (ii) is satisfied with $\Theta=1$ and
$\Psi=\Psi^{\mathrm{\tiny cat}}$.

###  5.2 Poisson spatio-temporal process

Spatio-temporal Poisson process models the situation where the state
$\omega_{tk}$ of the $k$-th location is a non-negative integer (that is,
$\mu=1$ and $Z=\\{0,1,2,\ldots\\}$), and the conditional (by the past)
distribution of $\zeta_{t}=[\omega_{t1};\omega_{t2};\ldots;\omega_{tL}]$ is
$k$ independent Poisson components $\omega_{tk}$ with parameters specified by
$\Phi_{k}(\eta^{T}(\zeta_{t-d}^{t-1})\beta)$. In order to meet (20) such that
the expectation of the Poisson random variable is nonnegative, we should
assume now that we are given a closed convex set
$\overline{Z}\subset{\mathrm{Dom}\,}\Phi$ such that

$\Phi(z)\geq 0\,\forall z\in\overline{Z}$

and that for all nonnegative integral vectors $z_{s}\in{\mathbf{R}}^{L}$,
$1\leq s\leq d$, it holds

$\eta^{T}(z_{1},\ldots,z_{d})x\in\overline{Z}\,\,\forall x\in{\cal B}.$

Observe that in the resulting stationary GGLM for all $t\geq 1$ and all
$h=[h_{1};\ldots;h_{L}]\in{\mathfrak{Z}}$, setting
$\nu_{t}=\Phi(H_{t}^{T}\beta)$, based on the conditional independence across
locations, and applying the Poisson moment generating function for each
location

${\mathbf{E}}\\{\exp\\{t\xi\\}\\}=\exp\\{\lambda
e^{t}-1\\},\quad\xi\sim\mbox{Poisson}(\lambda),$

we have

$\begin{array}[]{rcl}{\mathbf{E}}_{|\upsilon^{t-1},H_{t}}\left\\{\exp\\{h^{T}\zeta_{t}\\}\right\\}&=&\exp\left\\{\sum_{k=1}^{L}[\nu_{t}]_{k}(\exp\\{h_{k}\\}-1)\right\\},\\\
\lx@intercol\Rightarrow\hfil\lx@intercol&&\\\
\ln\left({\mathbf{E}}_{|\upsilon^{t-1},H_{t}}\left\\{\exp\\{h^{T}\zeta_{t}\\}\right\\}\right)&=&h^{T}\nu_{t}+\sum_{k=1}^{L}[\nu_{t}]_{k}\left[\exp\\{h_{k}\\}-h_{k}-1\right]\\\
&=&h^{T}\Phi(H_{t}^{T}\beta)+\sum_{k=1}^{L}[\nu_{t}]_{k}\left[\exp\\{h_{k}\\}-h_{k}-1\right].\\\
\end{array}$ (22)

Now let us set

$\chi_{t}(\upsilon^{t-1},H_{t})=\max_{u\in{\cal
B}}\|\Phi(H^{T}_{t}u)\|_{\infty},\,\,f(s)=\exp\\{s\\}-s-1.$

Then $\|\nu_{t}\|_{\infty}\leq\chi_{t}(\upsilon^{t-1},H_{t})$ and

$\begin{array}[]{l}\sum_{k=1}^{L}[\nu_{t}]_{k}\left[\exp\\{h_{k}\\}-h_{k}-1\right]\\\
\leq\|\nu_{t}\|_{\infty}\|[f(h_{1});f(h_{2});\ldots;f(h_{L})]\|_{1}\leq\chi_{t}(\upsilon^{t-1},H_{t})\max_{g}\left\\{E(g):\|g\|_{1}\leq\|h\|_{1}\right\\},\\\
\end{array}$

where $E(g)=\|[f(g_{1});...;f(g_{L})]\|_{1}$. Since $f$ is nonnegative and
convex on ${\mathbf{R}}$, the function $E(g)$ is convex, and therefore its
maximum over the $\ell_{1}$-ball $\\{g:\|g\|_{1}\leq\|h\|_{1}\\}$ is achieved
at a vertex, implying that

$\max_{g}\\{E(g):\|g\|_{1}\leq\|h\|_{1}\\}=f(\|h\|_{1}).$

Thus, setting

$\Psi_{t,\upsilon^{t-1},H_{t}}(r)=\chi_{t}(\upsilon^{t-1},H_{t})\Psi(r),\,\Psi(r):=\exp\\{r\\}-r-1,\,\,P_{t}=I_{L},$

and invoking (22), we conclude that

${\mathbf{E}}_{|\upsilon^{t-1},H_{t}}\left\\{\exp\\{h^{T}\zeta_{t}\\}\right\\}\leq
h^{T}\Phi_{t}(H_{t}^{T}\beta)+\Psi_{t,\upsilon^{t-1},H_{t}}(\|P_{t}h\|_{1}),$

as required by (8). Note that $\Psi_{t,\upsilon^{t-1},H_{t}}(\cdot)$, is,
modulo computational aspects, observable at time $t$.

###  5.3 Numerical illustration for Poisson case

Now we present some numerical examples of the Poisson spatio-temporal process
case to demonstrate the computational procedure for estimation and online
bounds for the recovery error. In the experiments to be reported, we dealt
with the identity link function (thus, the non-negativity of the Poisson
intensity is ensured by constraints on the feasible parameters.) The process
was

$[\zeta_{t}]_{k}\sim\mathrm{Poisson}\left(\beta^{0}_{k}+\sum_{s=1}^{d}\sum_{\ell=1}^{L}\beta^{s}_{k\ell}[\zeta_{t-s}]_{\ell}\right),\,t\geq
1,1\leq k\leq L.$

Recall that $d$ is the memory depth, $L$ is the number of locations, and
$[\zeta_{t}]_{k}$ is the $k$-th coordinate location of response at time $t$.

##### Generating parameter $\beta$.

Since $\Phi$ is the identity, selection of $\beta$ should ensure nonnegativity
of the quantities

$\beta^{0}_{k}+\sum_{s=1}^{d}\sum_{\ell=1}^{L}\beta^{s}_{k\ell}[\zeta_{t-s}]_{\ell},\,k\leq
L,$

whatever be nonnegative integral $\zeta_{s}\in{\mathbf{R}}^{L}$. To this end,
we restrict $\beta$ to be nonnegative. Besides this, with positive
$\beta^{0}_{i}$, in order for the process not to explode, we should have

$\sum_{s=1}^{d}\sum_{\ell=1}^{L}\beta^{s}_{k\ell}<1,\,k\leq L.$

Our generation was governed by two parameters $a>0$ and $b\leq 1$,
specifically, we ensured that

$\beta\geq 0\ \&\ \beta^{0}_{k}=a\ \&\
\sum_{s=1}^{d}\sum_{\ell=1}^{L}\beta^{s}_{k\ell}=b,\,\,k\leq L.$ (23)

Generation was organized as follows: given $d$, $L$, we selected $\beta$
satisfying (23) at random

The set ${\cal B}$ – a priori localized for $\beta$ – also was given by
parameters $a,b$ according to

${\cal B}=\\{x=\\{x^{0}_{k},x^{s}_{k\ell},1\leq s\leq d,1\leq k,\ell\leq
L\\}:x\geq 0,x^{0}_{k}\leq
1.1a,\sum_{s=1}^{d}\sum_{\ell=1}^{L}x^{s}_{k\ell}\leq 1.1b\,\forall k.\\}$

##### Building estimate $\widehat{\beta}$.

Our estimate was exactly as explained in the special case of Section 3 with
$\Lambda_{t}\equiv 1$, that is, $\widehat{\beta}(\upsilon^{N})$ was an optimal
solution to the Least Squares problem

$\min\limits_{x\in{\cal B}}\left[{1\over
2N}\sum_{t=1}^{N}\|H_{t}^{T}x-\zeta_{t}\|_{2}^{2}:={1\over
2N}\sum_{t=1}^{N}\sum_{k=1}^{L}\left([\zeta_{t}]_{k}-x^{0}_{k}-\sum_{s=1}^{d}\sum_{\ell=1}^{L}x^{s}_{k\ell}[\zeta_{t-s}]_{\ell}\right)^{2}\right].$

##### On-line bounds.

The major goal of our experimentation was to get online upper bounds on the
recovery errors. As it turned out, bounds stemming from Proposition 4.2 under
the circumstances were quite poor – the values of the resulting error bounds
typically were larger than the corresponding sizes of our a priori parameter’s
localizer ${\cal B}$. We used the approach described in Section 4.2.1 to get
meaningful error bounds. The choice of $M$ and $\alpha_{i}$’s will be
explained later. Specifically, we

  1. 1.

selected $M$ of policies ${\cal R}^{\iota}$, $1\leq\iota\leq M$, of generating
$L$-dimensional row vectors $\widetilde{E}_{t}$ as deterministic functions of
$\upsilon^{t-1}$ and $H_{t}$, resulting in $M$ affine functions

$G^{\iota}_{\upsilon^{N}}(x)={1\over
N}\sum_{t=1}^{N}\widetilde{E}_{t}^{\iota}[H_{t}x-\zeta_{t}];$

($\widetilde{E}_{t}^{\iota}$ are yielded by ${\cal R}^{\iota}$.)

  2. 2.

used Proposition 4.1.i and the preceding results of this section to build
online $(1-\epsilon/M)$-reliable upper bounds on the quantities
$|G^{\iota}_{\upsilon^{N}}(\beta)|$, specifically, the bounds

$\begin{array}[]{c}\delta^{\iota}(\upsilon^{N})=\min_{i}\left[\alpha_{i}\ln(2KM/\epsilon)+\alpha_{i}\sum_{t=1}^{N}\Psi_{t}(\alpha_{i}^{-1}N^{-1}\|\widetilde{E}^{\iota}_{t}\|_{1})\right]\\\
\Psi_{t}(r)=\chi_{t}[\exp\\{r\\}-r-1],\\\ \chi_{t}=\max\limits_{x\in{\cal
B}}\|H_{t}^{T}x\|_{\infty},\end{array}$ (24)

on $|G^{\iota}_{\upsilon^{N}}(\beta)|$.

  3. 3.

specified $(1-\epsilon)$-confidence interval

$\Delta(\upsilon^{N})=\\{x\in{\cal
B}:|G^{\iota}_{\upsilon^{N}}(x)|\leq\delta^{\iota}(\upsilon^{N}),1\leq\iota\leq
M\\}$

– under the circumstances, this is a convex compact set that contains $\beta$
with probability at least $1-\epsilon$.

  4. 4.

built induced confidence intervals $\Delta^{0}_{k}$, $\Delta^{s}_{k\ell}$ on
the entries in $\beta$, the endpoints of these intervals being the minima and
the maxima of the corresponding entries in $x$ taken over
$x\in\Delta(\upsilon^{N})$.
By construction, the probability for all entries in $\beta$ to be covered by
the corresponding intervals is at least $1-\epsilon$, so that the maximal,
over points from the intervals, deviations of the points from the
corresponding entries in $\widehat{\beta}$ form a vector which is an
$(1-\epsilon)$-reliable upper bound on the vector of magnitudes of the
recovery errors.

The essence of the matter is how to specify the policies ${\cal R}^{\iota}$.
The first $\kappa$ of our policies resulting in basic online bounds were to
take, as $\widetilde{E}^{\iota}_{t}$, the $\iota$’th row in $H_{t}$, resulting
in

$[G^{1}_{\upsilon^{N}}(x);\ldots;G^{\kappa}_{\upsilon^{N}}(x)]=F_{\upsilon^{N}}(x).$

The remaining policies were more sophisticated and inspired by the following
consideration.

> Given index $k\leq\kappa$ of an entry in $\beta$ and a “scale parameter”
> $\theta>0$, let us impose on $\widetilde{E}_{t}^{\iota}$’s the restriction
> $\|\widetilde{E}_{t}^{\iota}\|_{1}\leq\theta$, thus imposing upper bounds
> $\Psi_{t}(\alpha_{i}^{-1}N^{-1}\theta)$ on the terms
> $\Psi_{t}(\alpha_{i}^{-1}N^{-1}\|\widetilde{E}^{\iota}_{t}\|_{1})$ in (24),
> and under this assumption let us try to select $\widetilde{E}_{t}^{\iota}$’s
> in order to make the slope of the associated affine function
> $G_{\upsilon^{N}}^{\iota}(\cdot)$ as close as possible to $k$-th basic orth
> in ${\mathbf{R}}^{\kappa}$. Assume for a moment that we fully succeeded and
> the slope is exactly this basic orth. Then we get at our disposal inequality
> an inequality for $s$:
>
> $|s-c|\leq\delta=\delta^{\iota}(\upsilon^{N}),\,c:={1\over
> N}\sum_{t=1}^{N}\widetilde{E}_{t}^{\iota}H_{t}\zeta_{t},$
>
> in variable $s$ with hopefully small $\delta$ with is with overwhelming
> probability satisfied by $k$-th entry of $\beta$.
>
> The outlined methodology reduces to finding $L$-dimensional row vectors
> $\widetilde{E}_{t}$ of $\|\cdot\|_{1}$-norm not exdceeding a given $\theta$
> such that the vector
>
> ${1\over N}\sum_{t=1}^{T}H_{t}\widetilde{E}_{t}^{T}$
>
> is as close as possible to $k$-th basic orth $e_{k}$ in
> ${\mathbf{R}}^{\kappa}$. The difficulty here is that $\widetilde{E}_{t}$
> should be “non-anticipating” – they should be specified in terms of
> $\nu^{t-1}$ and $H_{t}$. The simplest policy of this type is the advanced
> policy
>
> $\begin{array}[]{rcl}f_{1}&=&e_{k},\\\
> \widetilde{E}_{t}&\in&\mathrm{Argmin}_{g}\left\\{\|f_{t}-N^{-1}H_{t}g\|:\|g\|_{1}\leq\theta\right\\},\\\
> f_{t+1}&=&f_{t}-N^{-1}\widetilde{E}_{t}.\\\ \end{array}$
>
> In our experiments, the only one of the standard norms which worked in this
> advanced policy was $\|\cdot\|_{2}$, and these were the policies we used. To
> implement such a policy, we need to specify $\theta$. This parameter is
> responsible for the tradeoff between the closeness of
> $e_{k}(\upsilon^{N}):={1\over N}\sum_{t=1}^{T}H_{t}\widetilde{E}_{t}^{T}$ to
> our target $e_{k}$ and online upper bound
>
>
> $\delta=\min_{i}\left[\alpha_{i}\ln(2KM/\epsilon)+\sum_{t=1}^{N}\alpha_{i}\Psi_{t}(\alpha_{i}^{-1}N^{-1}\theta)\right]$
>
> on $|e_{k}^{T}(\upsilon^{n})\beta-{1\over
> N}\sum_{t=1}^{N}\widetilde{E}_{t}H_{t}\zeta_{t}|$ we end up with. Resolving
> this tradeoff analytically seems to be an impossible task, but we can use
> the same computation-oriented approach as with optimizing in the scale
> parameter $\alpha$, specifically, select in advance a finite grid
> $\Theta=\\{\theta_{k}:k\leq K\\}$ of values of $\theta$ and run $K$ advanced
> policies parameterized by $\theta\in\Theta$.

##### “Advanced bounds.”

In our implementation, we selected a $K$-element set $\Theta$ on the positive
ray and, on top of already explained $\kappa$ policies ${\cal R}^{\iota}$,
$\iota\leq\kappa$, given by $\widetilde{E}_{t}=\mathrm{Row}_{\iota}(H_{t})$
used $\kappa K$ policies more, associating with every $k\leq\kappa$, the $K$
aforementioned advanced policies with common target $e_{k}$ and different
values of $\theta$ running through $\Theta$.

##### Implementation and numerical results.

In the experiments to be reported, we used $d=L=5$, resulting in $\kappa=130$,
$N=100,000$, $\epsilon=0.01$, and $M=780=5\times 130+130$ \- according to the
above explanation, $M$ is $\kappa$ plus $\kappa$ times the cardinality of
$\Theta$.

$\Gamma=\\{\alpha_{i}:=10^{0.25i-4},0\leq i\leq
36\\},\,\,\Theta=\\{0.5,0.75,1.00,1.25,2\\}.$

Numerical results. We present numerical results for two (in fact, quite
representative) experiments. In the first (Experiment A) we used $a=1$, $b=1$,
resulting in self-exciting process with growing with time $t$ total, over
locations, number $\omega^{\hbox{\tiny total}}_{t}:=\sum_{k=1}^{L}\omega_{tk}$
of events at time $t$. In the second (Experiment B), we used $a=1$, $b=0.5$,
resulting in a stable process, as shown in Fig. 2.

(a) Experiment A, ${1\over LN}\sum_{t}\omega^{\hbox{\tiny
total}}_{t}\approx$1.7e4

(b) Experiment B, ${1\over LN}\sum_{t}\omega^{\hbox{\tiny total}}_{t}=2.00$

Figure 2: Sample trajectories of $\omega^{\hbox{\tiny total}}_{t}$ versus $t$.

##### Recovery errors.

Now, we report the graphical and numerical data on actual recovery errors and
online bounds on these errors. On the plots and in the tables, “basic” online
bounds stem from the $\kappa$ initial policies ${\cal R}^{\iota}$,
$\iota\leq\kappa$, while “advanced” bounds stem from all $\kappa(K+1)$
policies described above.

(a) Experiment A

(b) Experiment B

Figure 3: Coordinate-wise recovery errors and online upper bounds on these
errors. Red dots: actual coordinate-wise recovery errors; blue dots: basic
online bounds; black dots: advanced online bounds. Ordering of coordinates
makes the errors non-decreasing. Table 1: Recovery error and bounds for
simulated examples.

| $\min$ | median | mean | $\max$
---|---|---|---|---
actual errors | 0.0000/0.0000 | 0.0005/0.0004 | 0.0402/0.0018 | 1.0000/0.0312
bounds, basic | 0.0962/0.0962 | 1.0775/1.0725 | 1.0045/1.0007 | 1.100/1.100
bounds, advanced | 0.0041/0.0041 | 0.0983/0.0932 | 0.1215/0.0824 | 1.100/0.2053

(a) Experiment A, coordinate-wise errors and their online bounds

| $\min$ | median | mean | $\max$
---|---|---|---|---
actual errors | 0.0000/0.0000 | 0.0007/0.0006 | 0.0024/0.0013 | 0.0413/0.0064
bounds, basic | 0.1216/0.1991 | 0.5485/0.5487 | 0.5140/0.5293 | 0.5500/0.5500
bounds, advanced | 0.0018/0.0018 | 0.0161/0.0158 | 0.0183/0.0152 | 0.1155/0.0249

(b) Experiment B, coordinate-wise errors and their online bounds

In Table 1(b), the first number in a cell represents the “column quantity,”
say, the median, taken over all $\kappa=130$ coordinate-wise errors and their
upper bounds, while the second number corresponds to recovery errors/bounds on
the errors for 125 “location to location influence” parameters
$\beta^{s}_{k\ell}$. This distinction is because in Experiment A, the “birth
rates” $\beta^{0}_{k}$ do not admit nontrivial estimates. Indeed, with our
setup, these birth rates are equal to 1, while the typical observation in the
unstable case is of the order of 104; of course, there is absolutely no way to
recover birthrates of the order of 1 on the “background” of magnitude 104. In
Experiment A, the estimates of birth rates were nearly zero; this quite
natural fact is responsible for the “large” – equal to 1 – uniform recovery
error in observed in this experiment. As the corresponding table shows, the
uniform norm of recovering all but 5 “birth rate” coefficients in $\beta$ is
0.03 - 30 times smaller than the uniform norm of recovering the entire
$\beta$. A similar albeit less profound phenomenon takes place in Experiment
B.

From the above data, it is clear that utilizing advanced policies improved
significantly the quality of online error bounds. While even after this
improvement, the advanced bounds are much larger than the actual errors, these
bounds seem to yield some information.

##### Upper-bounding $\|F_{\upsilon^{N}}(\beta)\|_{\infty}$.

Our experiments show that the online upper bound on
$\|F_{\upsilon^{N}}(\beta)\|_{\infty}$ is within factor 5-10 of the actual
value of the quantity.

##### Predictive power.

One way to utilize an estimate $\widehat{\beta}$ of $\beta$ is to use the
observations obtained so far to predict the states in the locations in the
future. Specifically, given observations $\omega^{t}:=\\{\omega_{\tau
k}:\tau\leq t,k\leq L\\}$, and prediction step $p\geq 1$, we can write a
simple recurrence specifying the conditional, $\omega^{t}$ given, expectations
$\bar{\omega}_{t+p,k}$ of the states $\omega_{t+p,k}$ and use these
expectations to build confidence intervals for $\omega_{t+p,k}$ as if the
conditional, $\omega^{t}$ given, distributions of these random variables were
$\mathrm{Poisson}(\bar{\omega}_{t+p,k})$ (“as if” reflects the fact that the
actual conditional distributions in question probably are
$\mathrm{Poisson}(\bar{\omega}_{t+p,k})$ only when $p=1$). To write down the
aforementioned recurrence, we need to know the process’s true parameters
$\beta$; by replacing these parameters with their estimates $\widehat{\beta}$,
we get “empirical” confidence intervals for future states. Table 2(b) shows
how the prediction worked in the testing sample.

Table 2: Average, over 100,000-element testing sample, frequency for $p$-step-
ahead empirical 0.95-confidence intervals not to cover the actual states of
locations $k$, $1\leq k\leq 5$. Last row: trivial prediction with
$\bar{\omega}_{t+p,k}$ set to $\omega_{tk}$ (“tomorrow will be same as
today”).

$p$ | $k=1$ | $k=2$ | $k=3$ | $k=4$ | $k=5$
---|---|---|---|---|---
1 | 0.049 | 0.050 | 0.050 | 0.048 | 0.050
2 | 0.049 | 0.051 | 0.050 | 0.050 | 0.050
3 | 0.088 | 0.054 | 0.054 | 0.054 | 0.068
4 | 0.088 | 0.056 | 0.057 | 0.058 | 0.113
5 | 0.135 | 0.063 | 0.063 | 0.061 | 0.116
6 | 0.137 | 0.069 | 0.068 | 0.080 | 0.142
7 | 0.166 | 0.075 | 0.070 | 0.084 | 0.162
8 | 0.170 | 0.080 | 0.077 | 0.093 | 0.170
9 | 0.194 | 0.087 | 0.081 | 0.106 | 0.191
10 | 0.198 | 0.093 | 0.088 | 0.112 | 0.205
| 0.322 | 0.156 | 0.182 | 0.184 | 0.223

(a) Experiment A

$p$ | $k=1$ | $k=2$ | $k=3$ | $k=4$ | $k=5$
---|---|---|---|---|---
1 | 0.015 | 0.015 | 0.014 | 0.014 | 0.014
2 | 0.016 | 0.015 | 0.014 | 0.015 | 0.014
3 | 0.016 | 0.015 | 0.015 | 0.015 | 0.020
4 | 0.017 | 0.015 | 0.015 | 0.016 | 0.021
5 | 0.021 | 0.016 | 0.015 | 0.017 | 0.021
6 | 0.021 | 0.017 | 0.018 | 0.018 | 0.023
7 | 0.021 | 0.018 | 0.018 | 0.019 | 0.024
8 | 0.023 | 0.018 | 0.018 | 0.019 | 0.025
9 | 0.025 | 0.019 | 0.019 | 0.019 | 0.026
10 | 0.025 | 0.019 | 0.019 | 0.019 | 0.026
| 0.202 | 0.183 | 0.180 | 0.181 | 0.203

(b) Experiment B

---
Experiment A
Experiment B
Figure 4: Prediction on time interval $50,000\leq t\leq 50,050$. Top to
bottom: locations. Blue crosses: actual states; 10 orange vertical segments:
in-between the crosses: 0.95-confidence intervals with prediction steps 1, 2,
…, 10; black vertical segments: confidence intervals yielded by trivial
prediction.

##  6 Real Data illustrations

###  6.1 Dataset

We apply the Poisson point process model to the wildfire dataset. The dataset
contains wildfire occurrences in California from 06/2014 to 12/2019 at various
locations. The dataset is retrieved from the California Public Utilities
Commission333Website: https://www.cpuc.ca.gov/wildfires/. We discretize the
space containing the fire incidents into 26 disjoint rectangle cells (i.e.,
$L=26$); see Figure 5 for the raw data (grey circles) and discretized cells
(blue rectangles). Fire incidents are grouped by months, so each
$[\zeta_{t}]_{k}$ denotes the number of wildfire incidents at location $k$ in
month $t=1,\ldots,67$, after space discretization. Data from all locations and
all except the last year (i.e., 55 months) are used for training, and data in
the last year (i.e., 12 months) are for testing. We choose the memory depth
$d=12$, assuming observations up to 12 months in the past may influence the
present intensities.

(a) Raw data and discretized space

(b) Total number of fires per month

Figure 5: Figure (a) visualizes the wildfire data for training and test set,
shown as scatters on the terrain map and spatial discretization for the
Poisson process model. Figure (b) visualizes the total number of fire
incidents over all locations at each unit time (month) $t$ (including both
training and test time).

Figure 6: Visualization of estimated interaction parameters
$\hat{\beta}^{s}_{kl}$ over neighbors $l=1,\ldots,L$ and lags $s=1,\ldots,d$.
We select two locations $k$, where the left (resp. right) figure shows the
heatmap of these interaction parameters at location 13 (resp. 7), with the
highest (resp.) lowest fire incidents in training data. The estimated baseline
intensity $\hat{\beta}^{0}_{k}$ at $k=13$ and $k=7$ are 0.0088 and 6.2$\times
10^{-5}$, respectively.

(a) Lag $s=1$

(b) Lag $s=6$

(c) Lag $s=12$

Figure 7: Overlaying base rates $\hat{\beta}^{0}_{k}$ and interaction
parameters $\hat{\beta}^{s}_{kl}$ on a terrain map over different lag $s$. In
each subfigure, the base rates are plotted as red solid circles, with color
and size proportional to the value of $\hat{\beta}^{0}_{k}$. The interactions
$\hat{\beta}^{s}_{kl}$ are plotted as blue-directed edges from location $l$ to
location $k$, with width proportional to the magnitude.

###  6.2 Model formulation and estimation

Following the simulation examples, we assume the underlying point process is

$[\zeta_{t}]_{k}\sim\mathrm{Poisson}\left(\beta^{0}_{k}+\sum_{s=1}^{d}\sum_{\ell=1}^{L}\beta^{s}_{k\ell}[\zeta_{t-s}]_{\ell}\right),\,t\geq
1,1\leq k\leq L.$ (25)

The parameters
$\beta=\\{\beta^{0}_{k},\beta^{s}_{kl}\\},\beta\in\mathbb{R}^{\kappa}$ are
estimated as an optimal solution to the unconstrained least-square objective

$\min\limits_{x\in\mathbb{R}^{\kappa}}\left[{1\over
2N}\sum_{t=1}^{N}\sum_{k=1}^{L}\left([\zeta_{t}]_{k}-x^{0}_{k}-\sum_{s=1}^{d}\sum_{\ell=1}^{L}x^{s}_{k\ell}[\zeta_{t-s}]_{\ell}\right)^{2}\right].$
(26)

During prediction, in order to ensure non-negativity of estimated rates in
(25), we use
$\max\\{0,\beta^{0}_{k}+\sum_{s=1}^{d}\sum_{\ell=1}^{L}\beta^{s}_{k\ell}[\zeta_{t-s}]_{\ell}\\}$
as the estimated rates. We note that unlike a constrained problem, which may
require non-negativity in the entries of $\beta$, the unconstrained problem in
(26) allows for greater flexibility. Specifically, it can capture possibly
negative spatial-temporal influences through the parameters $\beta^{s}_{kl}$.

###  6.3 Evaluation metrics

We evaluate the model’s performance on average and in specific instances.

We consider several metrics for average performance. Below, $r_{k}$ denotes
the observed average count (frequency) for a location $k$, and $\hat{r}_{k}$
denotes the conditional expectation for location $k$ based on the parameters
estimated using training data from time 1 to $N$. We also consider a
simulation-based test, $\hat{r}_{k}^{\text{simul}}$, which has been used in
[2, Table V] for Bernoulli processes; the rationale is we model the
probability distribution of counts, and thus, can check whether or not our
model can generate random instances with an empirical frequency that matches
the observed frequency of the data. Specifically, we sample from a Poisson
distribution with predicted intensity to generate replicas
$[\hat{\zeta}_{t}]_{k}^{(i)},i=1,\ldots,I$ for $I=100$. Finally, we also
compare with a baseline using the seasonal prediction: in (30), the estimates
$[\tilde{\zeta}_{t}]_{k}$ are seasonal averages (with a period of 12 months)
of past samples $\\{[\zeta_{t^{\prime}}]_{k}:(t-t^{\prime})\equiv 0\
(\text{mod}~{}12),t^{\prime}\leq N\\}$. In our case, such predictions are
meaningful due to seasonality in raw data, as shown in Figure 5(b). We thus
compute the following quantities:

$\displaystyle r_{k}$
$\displaystyle:=\frac{1}{T-N}\sum_{t=N+1}^{T}[\zeta_{t}]_{k}$ (27)
$\displaystyle\hat{r}_{k}$
$\displaystyle:=\frac{1}{T-N}\sum_{t=N+1}^{T}\left(\hat{\beta}^{0}_{k}+\sum_{s=1}^{d}\sum_{\ell=1}^{L}\hat{\beta}^{s}_{k\ell}[\zeta_{t-s}]_{\ell}\right)$
(28) $\displaystyle\hat{r}_{k}^{\text{simul}}$
$\displaystyle:=\frac{1}{(T-N)I}\sum_{t=N+1}^{T}\sum_{i=1}^{I}[\hat{\zeta}_{t}]_{k}^{(i)}$
(29) $\displaystyle\hat{r}_{k}^{\text{seasonal}}$
$\displaystyle:=\frac{1}{T-N}\sum_{t=N+1}^{T}[\tilde{\zeta}_{t}]_{k}.$ (30)

Table 3: For wildfire data: the comparison of true average counts as in (27)
versus estimated average counts as in (28)–(30).

Location $r_{k}$ $\hat{r}_{k}$ $\hat{r}_{k}^{\text{simul}}$
$\hat{r}_{k}^{\text{seasonal}}$ 0 0.333 0.706 0.738 0.729 1 1.833 1.003 0.994
0.988 2 1.917 1.441 1.528 1.371 3 0.083 0.220 0.201 0.204 4 1.167 1.076 1.092
1.029 5 2.083 1.084 1.097 1.021 6 1.917 1.258 1.302 1.254 7 0.000 0.047 0.052
0.017 8 4.417 3.308 3.212 3.354 9 1.167 1.217 1.200 1.221 10 4.917 3.060 3.062
2.971 11 0.500 0.484 0.508 0.442 12 1.250 1.438 1.444 1.400 (a) First 13
locations Location $r_{k}$ $\hat{r}_{k}$ $\hat{r}_{k}^{\text{simul}}$
$\hat{r}_{k}^{\text{seasonal}}$ 13 5.333 4.000 4.123 3.996 14 2.167 2.334
2.357 2.392 15 0.917 0.763 0.757 0.746 16 0.083 0.063 0.065 0.038 17 1.167
1.567 1.658 1.529 18 2.583 2.527 2.514 2.388 19 2.667 3.832 3.899 3.596 20
0.000 0.055 0.057 0.050 21 0.167 0.274 0.278 0.208 22 2.250 2.691 2.734 2.683
23 0.333 0.196 0.177 0.204 24 0.083 0.049 0.054 0.021 25 0.333 0.367 0.386
0.350 (b) Rest 13 locations $\hat{r}_{k}$ $\hat{r}_{k}^{\text{simul}}$
$\hat{r}_{k}^{\text{seasonal}}$ MAE 0.421 0.424 0.425 (c) Mean Absolute Error
(MAE) between $r_{k}$ and the estimates

###  6.4 Results

We first visualize estimated parameters, both as heatmaps over lag
$s=1,\ldots,d,$ in Figure 6 and edges overlaid on a terrain map in Figure 7.
Figure 6 shows that at location $k=13$ with the most fire incidents, the
estimated $\beta^{s}_{kl}$ are much higher in magnitude than those at location
$k=7$ (with the least fire incidents). The large magnitude parameters indicate
significant spatial-temporal influence received by location $k$ from its
neighboring and past fire incidents. Such strong influences are also shown in
Figure 7, visualized as directed edges. Meanwhile, from both figures, we see
that values of $\hat{\beta}^{s}_{kl}$ tend to be smaller as the lag $s$
increases (especially at $s=1$ versus $s=6$), indicating a temporally decaying
influence from historical fire incidents.

Figure 8: Absolute error of metrics shown in Table 3. The locations on the
x-axis are sorted by the observed number of fire incidents in the training
data. Table 4: For wildfire data, empirical 0.95-confidence interval for
$p$-step-ahead prediction. We report top-2, middle-2, and bottom-2 locations
with the fire incidents level that are high, medium, and low in training data.
The second-to-bottom row labeled with persistent prediction refers to the
prediction using the previous observation$[\zeta_{t-1}]_{k}$, assuming
tomorrow will be the same as today. The bottom row predicts using seasonal
average $[\tilde{\zeta}_{t}]_{k}$ defined in (30).

$p$ | $k=13$ | $k=19$ | $k=4$ | $k=1$ | $k=7$ | $k=24$
---|---|---|---|---|---|---
1 | 0.167 | 0.000 | 0.000 | 0.167 | 0.0 | 0.000
2 | 0.240 | 0.000 | 0.000 | 0.182 | 0.0 | 0.000
3 | 0.246 | 0.002 | 0.001 | 0.200 | 0.0 | 0.000
4 | 0.242 | 0.002 | 0.002 | 0.222 | 0.0 | 0.000
5 | 0.267 | 0.000 | 0.000 | 0.250 | 0.0 | 0.000
6 | 0.301 | 0.000 | 0.007 | 0.285 | 0.0 | 0.000
7 | 0.333 | 0.000 | 0.004 | 0.333 | 0.0 | 0.000
8 | 0.400 | 0.000 | 0.012 | 0.400 | 0.0 | 0.000
9 | 0.500 | 0.000 | 0.014 | 0.250 | 0.0 | 0.000
10 | 0.333 | 0.000 | 0.018 | 0.333 | 0.0 | 0.000
11 | 0.000 | 0.000 | 0.029 | 0.000 | 0.0 | 0.000
12 | 0.000 | 0.000 | 0.103 | 0.000 | 0.0 | 0.000
Persistent pred, $p=1$ | 0.417 | 0.167 | 0.167 | 0.250 | 0.0 | 0.083
Seasonal average pred, $p=1$ | 0.083 | 0.083 | 0.000 | 0.167 | 0.0 | 0.083

Figure 9: Instance-based “online” prediction intervals on test data at
locations selected in Table 4: observations are blue dots, and intervals are
in magenta; we note that the prediction interval achieves good coverage
probability with a reasonable width.

We then compare the performance of our model against baselines defined in
(30). Table 3(c) shows that, on average, the Mean Absolute Error (MAE) of
estimated rates by our model is lower than that of the baselines. Upon sorting
all locations by descending number of fire incidents, we see in Figure 8 that
our method yields very similar performance in most locations.

We lastly compare on an instance level the non-coverage of prediction
intervals as considered in simulation (i.e., Table 2(b)). For instance-based
performance, we report the $p$-step ahead ($p=1,\ldots,12$) interval non-
coverage (similar to Table 2(b) and Figure 4 for simulation examples). Here,
Table 4 shows $p$-step ahead coverage at different prediction horizons for the
wildfire data, where the non-coverage tends to increase as the model predicts
further into the future. We observe that our model incurs a lower non-coverage
than the persistent prediction and the seasonal average prediction at $p=1$
(i.e., one-step ahead prediction). Specifically, our model incurs a non-
coverage exceeding the target 0.05 only at 2 out of 6 locations, whereas the
other baselines incur non-coverage at 5 or 4 out of 6 locations. Figure 9
visualizes prediction intervals when $p=1$ at corresponding locations, where
intervals are also adaptive in widths at different locations and test times
with a reasonable width.

##  7 Conclusion

We have introduced a computational method to estimate and quantify the
uncertainty in the Generalized Generalized Linear Model (GGLM) for modeling
spatio-temporal data. Our method provides a parameter recovery guarantee as
well as uncertainty quantification through a set of concentration
inequalities. To assess the performance of our method, we conducted numerical
simulations and applied the method to a real-world problem of wildfire
prediction.

## References

  * [1] Alan Agresti.  Categorical data analysis, volume 482.  John Wiley & Sons, 2003.
  * [2] Anatoli Juditsky, Arkadi Nemirovski, Liyan Xie, and Yao Xie.  Convex parameter recovery for interacting marked processes.  IEEE Journal on Selected Areas in Information Theory, 2020.
  * [3] Anatoli B Juditsky and AS Nemirovski.  Signal recovery by stochastic optimization.  Automation and Remote Control, 80(10):1878–1893, 2019.
  * [4] Peter McCullagh.  Generalized linear models.  Routledge, 2018.
  * [5] John Ashworth Nelder and Robert WM Wedderburn.  Generalized linear models.  Journal of the Royal Statistical Society: Series A (General), 135(3):370–384, 1972.
  * [6] Yu Nesterov.  Semidefinite relaxation and nonconvex quadratic optimization.  Optimization Methods and Software, 9(1-3):141–160, 1998.