# Solution Path of Time-varying Markov Random Fields with Discrete
Regularization††thanks: Salar Fattahi is supported, in part, by grant 2152776
from the National Science Foundation, grant N00014-22-1-2127 from the Office
of Naval Research, and a MICDE Catalyst Grant. Andrés Gómez is supported, in
part, by grant 2152777 from the National Science Foundation.

Salar Fattahi
Industrial and Operations Engineering
University of Michigan
<EMAIL_ADDRESS>Andrés Gómez
Industrial and Systems Engineering
University of Southern California
<EMAIL_ADDRESS>
###### Abstract

We study the problem of inferring sparse time-varying Markov random fields
(MRFs) with different discrete and temporal regularizations on the parameters.
Due to the intractability of discrete regularization, most approaches for
solving this problem rely on the so-called maximum-likelihood estimation (MLE)
with relaxed regularization, which neither results in ideal statistical
properties nor scale to the dimensions encountered in realistic settings. In
this paper, we address these challenges by departing from the MLE paradigm and
resorting to a new class of constrained optimization problems with exact,
discrete regularization to promote sparsity in the estimated parameters.
Despite the nonconvex and discrete nature of our formulation, we show that it
can be solved efficiently and parametrically for all sparsity levels. More
specifically, we show that the entire solution path of the time-varying MRF
for all sparsity levels can be obtained in $\mathcal{O}(pT^{3})$, where $T$ is
the number of time steps and $p$ is the number of unknown parameters at any
given time. The efficient and parametric characterization of the solution path
renders our approach highly suitable for cross-validation, where parameter
estimation is required for varying regularization values. Despite its
simplicity and efficiency, we show that our proposed approach achieves
provably small estimation error for different classes of time-varying MRFs,
namely Gaussian and discrete MRFs, with as few as one sample per time.
Utilizing our algorithm, we can recover the complete solution path for
instances of time-varying MRFs featuring over 30 million variables in less
than 12 minutes on a standard laptop computer. Our code is available at
https://sites.google.com/usc.edu/gomez/data. 111An earlier version of this
paper [Fattahi and Gomez, 2021], which concentrated solely on sparsely-
changing Gaussian MRFs with a fixed parameter, was presented at the Neural
Information Processing Systems (NeurIPS) conference in 2021. The current
submission is a substantial expansion of our work on multiple fronts. Firstly,
we have shown that the complete solution path of time-varying MRFs can be
recovered for all regularization parameters. Secondly, our new approach
enables us to tackle a broader class of temporal changes beyond sparsity.
Lastly, our statistical guarantees apply to a more extensive class of MRFs,
specifically discrete MRFs.

##  1 Introduction

Most modern systems are massive-scale with a hierarchy of unknown and ever-
changing topologies that must be estimated in real-time. For instance, the
real-time inference of time-varying stock correlation networks is crucial for
the identification of sharp changes in market conditions, and a need to
rebalance the portfolio [Talih and Hengartner, 2005, Polson and Tew, 2000,
Hallac et al., 2017]. Another example is the inference of temporal gene
regulatory networks based on gene expression data, which has immediate
applications in early diagnosis of cancer and other disease processes
[Hartemink et al., 2000, Karlebach and Shamir, 2008, Ravikumar et al., 2023].

The behavior of time-varying networks, such as those mentioned above, can be
captured by time-varying Markov Random Fields (MRF). Time-varying MRFs are
associated with a temporal sequence of undirected Markov graphs
$\mathcal{G}_{t}(V,E_{t})$, where $V$ is the set of nodes and $E_{t}$ is the
set of edges in the graph at time $t$. The node set $V$ represents the random
variables in the model, while the edge set $E_{t}$ captures the conditional
dependency between these variables that may change over time. In other words,
if two nodes $i$ and $j$ are not connected by an edge, then the random
variables $i$ and $j$ are independent at time $t$, conditioned on the rest of
the variables; a feature known as Markov property. The identification of such
conditional independence properties would lead to simpler and more
interpretable models, which are invaluable for improving human understanding
of different phenomena.

###  1.1 Maximum Likelihood Estimation

A popular approach for the inference of time-varying MRFs is based on the so-
called maximum likelihood estimation (MLE): to find a probabilistic graphical
model, based on which the observed data is most probable to occur [Wainwright
and Jordan, 2008]. Despite being known as theoretically powerful tools—a fact
noted as early as 1960s [Kalman, 1960, Shellenbarger, 1966]—MLE-based methods
face several fundamental challenges.

First, MLE-based approaches have overwhelmingly high computational complexity,
which limits their application to small- and medium-scale problems. This
challenge is further exacerbated in the time-varying setting where a graphical
model needs to be inferred for each time, thereby leading to a dramatic
increase in the number of parameters to be estimated. For instance, in order
to obtain an $\epsilon$-accurate solution, typical MLE-based methods have
complexity ranging from $\mathcal{O}(Tp^{6}\log(1/\epsilon))$ (via general
interior-point methods) [Mohan et al., 2014, Potra and Wright, 2000] to
$\mathcal{O}(Tp^{3}/\epsilon)$ (via tailored first-order methods, such as
ADMM) [Hallac et al., 2017, Ravikumar et al., 2010, Banerjee et al., 2008].
Here, $p$ is the number of unknown parameters, and $T$ is the number of time
steps. Solvers with such computational complexity fall short of practical use
in settings where the dimension $p$ is large, and the desired $\epsilon$ is
small. This prohibitive complexity of MLE-based methods is also exemplified in
their practical performance [Fattahi and Sojoudi, 2019, Zhang et al., 2018,
Fattahi et al., 2019].

Second, MLE-based methods fail to efficiently incorporate prior structural and
temporal information into their estimation procedure. For instance, it is
well-known that different time-varying networks exhibit sparse topologies,
which can be captured via a combinatorial or discrete regularizer like the
“$\ell_{0}$-norm”. However, due to the combinatorial nature of the
$\ell_{0}$-norm, MLE-based methods inevitably resort to relaxed or weaker
variants of such regularization, such as $\ell_{1}$-norm, thereby suffering
from inferior statistical properties [Bertsimas et al., 2019]. In particular,
a regularizer based on $\ell_{1}$-norm results in shrinkage of the entries of
the precision matrix, and the bias due to the penalty can dominate the
variance from the likelihood, ultimately leading to poor estimates. In §2.2,
we provide an example to illustrate the performance of such relaxation.

Finally, in most applications, one needs to obtain a statistically meaningful
regularization parameter via cross-validation [Cox and Hinkley, 1979], which
amounts to solving the MLE repeatedly or parametrically for different values
of the regularization parameter. However, obtaining parametric solutions to
optimization or inference problems is a highly non-trivial task even if the
fixed-parameter problem can be easily tackled. For example, while the well-
known lasso regression can be solved efficiently, parametric algorithms such
as LARS [Efron et al., 2004] have in general exponential complexity [Mairal
and Yu, 2012].

###  1.2 Our contributions

To address the aforementioned challenges, our approach departs from the
conventional wisdom in statistics and machine learning that the inference with
discrete regularizers is intractable, and convex surrogates should be used
instead. In fact, we prove the contrary: time-varying MRFs with discrete
regularizers can be inferred efficiently and parametrically; that is, for all
sparsity levels. Our contributions are as follows:

  * -

More tractable formulation: We introduce a much simpler formulation for the
inference of time-varying MRFs that can leverage and promote different
sparsity and smoothness structures. Unlike most MLE-based methods which
typically rely on convexity of the regularizer to ensure tractability, our
proposed formulation can directly incorporate the nonconvex and discrete
$\ell_{0}$ regularizer to promote sparsity in the estimated parameters and/or
their differences.

  * -

Recovering the solution path via dynamic programming: We propose a dynamic
programming approach that can recover the entire solution path of time-varying
MRFs for all sparsity levels, i.e., for all values of the regularization
coefficients controlling the sparsity, with complexity $\mathcal{O}(pT^{3})$
(Theorem 1), which scales linearly with the dimension $p$. This makes our
approach particularly suitable for cross-validation, where the goal is to find
the regularization parameter that leads to the best generalization performance
of the model.

  * -

Statistical guarantees: We show that our proposed estimation method enjoys a
strong statistical guarantee, provided that an approximate backward mapping of
the underlying distribution is available (Theorem 2). For two classes of
Gaussian and Discrete time-varying MRFs, we derive sharp non-asymptotic
guarantees for the estimated parameters, showing that they can be accurately
inferred with as few as one sample per time (Propositions 5 and 6). In the
case where the underlying Markov graphs and their temporal differences are
sparse, we establish high probability guarantees on the sparsistency of the
obtained solutions and their differences.

  * -

Implementation: Using our algorithm, we can recover the entire solution path
of an instance of time-varying MRFs with more than 30M variables in less than
12 minutes on a normal laptop computer. Our code, as well as our test cases,
are available at https://sites.google.com/usc.edu/gomez/data.

Notations. The $i^{th}$ element of a time-series vector $v_{t}$ is denoted as
$v_{t;i}$; the $(i,j)^{th}$ element of a time-indexed matrix $V_{t}$ is
denoted as $V_{t;ij}$. For a vector $v$, the notation $v_{i:j}$ is used to
denote the subvector of $v$ from index $i$ to $j$. The $\ell_{q}$-ball with
radius $\rho$ centered at $\bar{\mu}$ is defined as
$\mathcal{B}_{q}(\bar{\mu},\rho):=\\{\mu:\|\mu-\bar{\mu}\|_{q}\leq\rho\\}$.
For a vector $v$, the notations $\|v\|_{\infty}$, $\|v\|_{2}$, $\|v\|_{0}$
denote the $\ell_{\infty}$ norm, $\ell_{2}$ norm, and “$\ell_{0}$-norm”, i.e.,
the number of nonzero elements, respectively. Moreover, for a matrix $M$, the
notations $\|M\|_{2}$, $\|M\|_{\infty}$, $\|M\|_{1/1}$,
$\|M\|_{\infty/\infty}$ refer to the induced 2-norm, induced $\infty$-norm,
$\ell_{1}/\ell_{1}$ norm, and $\ell_{\infty}/\ell_{\infty}$ norm,
respectively. For a $d\times d$ matrix $M$, We define
$\|M\|_{\mathrm{off}}=\|M\|_{1/1}-\sum_{i=1}^{d}|M_{ii}|$. For a vector $v$
and matrix $M$, the notations $\mathrm{supp}(v)$ and $\mathrm{supp}(M)$ are
defined as the location of their nonzero elements. Given two sequences $f(n)$
and $g(n)$, the notation $f(n)\lesssim g(n)$ implies that there exists a
constant $C<\infty$ that satisfies $f(n)\leq Cg(n)$. Moreover, $f(n)\asymp
g(n)$ implies that $f(n)\lesssim g(n)$ and $g(n)\lesssim f(n)$. Given two
scalars $a$ and $b$, the symbols $a\wedge b$ and $a\vee b$ are used to denote
their minimum and maximum, respectively.

##  2 Problem Formulation

In this work, we consider a class of time-varying MRFs that can be expressed
as families of exponential distributions, defined as:

$\displaystyle\mathbb{P}(x_{t};\theta_{t}^{\star})=\exp\left\\{\langle\theta^{\star}_{t},\phi(x_{t})\rangle-A(\theta^{\star}_{t})\right\\}\qquad\text{for}\quad
t=0,\dots,T.$ (1)

Here, $x_{t}\in\mathbf{R}^{n}$ is the random variable at time $t$,
$\theta^{\star}_{t}$ is the (unknown) true canonical parameter from the domain
$\mathcal{D}\subseteq\mathbb{R}^{p}$, the function
$\phi:\mathbb{R}^{n}\to\mathbb{R}^{p}$ is the sufficient statistics, and
$A:\mathbb{R}^{p}\to\mathbb{R}$ is the log-partition function, which
normalizes the distribution. Special classes of time-varying MRFs include
time-varying Gaussian MRFs (GMRFs) and time-varying Discrete MRFs,
respectively corresponding to multivariate Gaussian and discrete random
processes. As will be explained later, due to the equivalence between MRFs and
exponential families, Markov graphs can be systematically obtained from their
canonical parameters $\\{\theta^{\star}_{t}\\}_{t=0}^{T}$ [Wainwright and
Jordan, 2008].

An alternative parameterization of exponential families is via mean or moment
parameters, i.e., the moments of the sufficient statistics
$\mu_{t}(\theta^{\star}_{t})=\mathbb{E}_{\theta^{\star}_{t}}[\phi(x_{t})]\in\mathcal{R}$,
where $\mathcal{R}$ is the so-called moment polytope. Given the true canonical
parameters $\\{\theta^{\star}_{t}\\}_{t=0}^{T}$, the mean parameters
$\\{\mu_{t}(\theta^{\star}_{t})\\}_{t=0}^{T}$ can be obtained via the forward
mapping $F:\mathcal{D}\to\mathcal{R}$, where
$\mu_{t}(\theta_{t})=F(\theta_{t})=\nabla A(\theta_{t})$. The conjugate (or
Fenchel) duality can be used to define the _backward mapping_
$F^{*}:\mathcal{R}\to\mathcal{D}$ with
$\theta_{t}(\mu_{t})=F^{*}(\mu_{t})=\nabla A^{*}(\mu_{t})$, where $A^{*}$ is
the conjugate dual of the log-partition function [Wainwright and Jordan, 2008,
Chapter 3]. For well-known classes of MRFs, the sufficient statistics,
canonical parameters, and backward and forward mappings can be characterized
based on the underlying probability distribution.

Figure 1:  (left) The Markov graph of a time-varying GMRF, (right) the
sparsity pattern of the precision matrix.

##### Gaussian setting.

The canonical parameters in a Gaussian time-varying MRF (GMRF) correspond to
the tuple of inverse covariance matrix $\Theta_{t}\in\mathbb{R}^{n\times n}$
(also known as the precision matrix) and mean vector
$\eta_{t}\in\mathbb{R}^{n}$. In this setting, the exponential distribution (1)
reduces to

$\displaystyle\mathbb{P}(x_{t})\\!=\\!\exp\left\\{-\frac{1}{2}\langle\Theta^{\star}_{t},x_{t}x_{t}^{\top}\rangle\\!-\\!A(\Theta^{\star}_{t})\right\\},$
(2)

with sufficient statistics $x_{t}x_{t}^{\top}$, where without loss of
generality we assumed that $\eta_{t}=\mathbb{E}[x_{t}]=0$. For GMRFs, the edge
set of the Markov graph $\mathcal{G}_{t}$ coincides with the off-diagonal
nonzero elements of $\Theta^{\star}_{t}$ [Weiss and Freeman, 2000]. Figure 1
shows the sparsity pattern of a precision matrix for a GMRF and its
corresponding Markov graph. The moment parameter $\mu_{t}$ in this setting
corresponds to the covariance matrix
$\Sigma_{t}^{\star}=\mathbb{E}[x_{t}x_{t}^{\top}]$. Accordingly, the forward
and backward mappings $F:\mathcal{D}\to\mathcal{R}$ and
$F^{*}:\mathcal{R}\to\mathcal{D}$ can be defined as
$F(\Theta_{t})=\Sigma_{t}=\Theta_{t}^{-1}$ and $F(\Sigma_{t})=\Sigma_{t}^{-1}$
with $\mathcal{D}=\\{\Theta_{t}:\Theta_{t}=\Theta_{t}^{\top},\Theta_{t}\succ
0\\}$ and
$\mathcal{R}=\\{\Sigma_{t}:\Sigma_{t}=\Sigma_{t}^{\top},\Sigma_{t}\succ 0\\}$.

##### Discrete setting.

In discrete time-varying MRFs (DMRF), each element of the random variable
$x_{t}$ takes a value from the discrete set $\mathcal{K}$ with cardinality
$|\mathcal{K}|=K$. For simplicity, we assume that $\mathcal{K}\subseteq[0,1]$.
In this case, (1) reduces to

$\displaystyle\mathbb{P}(x_{t})\\!=\\!\exp\left\\{\sum_{i\in[p],k\in\mathcal{K}}\theta^{\star}_{t;ik}\mathbb{I}[x_{t;i}=k]+\sum_{i,j\in[p],k,l\in\mathcal{K}}\theta^{\star}_{t;ijkl}\mathbb{I}[x_{t;i}=k,x_{t;j}=l]-A(\theta^{\star}_{t})\right\\},$
(3)

where $\theta^{\star}_{t;i\bm{\cdot}}\in\mathbb{R}^{K}$ is a vector parameter
for the node-wise indicator function $\mathbb{I}[x_{t;i}=k]$, and
$\theta^{\star}_{t;ij\bm{\cdot}\bm{\cdot}}\in\mathbb{R}^{K\times K}$ is a
matrix parameter for the edge-wise indicator function
$\mathbb{I}[x_{t;i}=k,x_{t;j}=l]$. The node- and edge-wise indicator functions
form an (over-complete) sufficient statistics. In this case, the edge set of
the Markov graph $\mathcal{G}_{t}$ corresponds to the nonzero matrix
parameters $\theta^{\star}_{t;ij\bm{\cdot}\bm{\cdot}}$; in other words,
$(i,j)\in E_{t}$ if and only if
$\theta^{\star}_{t;ij\bm{\cdot}\bm{\cdot}}\not=0_{K\times K}$. The moment
parameter, in this case, coincides with the node- and edge-wise marginal
probabilities, defined as
$\mu^{\star}_{t;ik}=\mathbb{E}[\mathbb{I}[x_{t;i}=k]]=\mathbb{P}(x_{t;i}=k)$
and
$\mu^{\star}_{t;ijkl}=\mathbb{E}[\mathbb{I}[x_{t;i}=k,x_{t;j}=l]]=\mathbb{P}(x_{t;i}=k,x_{t;j}=l)$.
Accordingly, the forward mapping $F:\mathcal{D}\to\mathcal{R}$ is a function
that computes the node- and pair-wise marginal probabilities given the
canonical parameters of the joint distribution. Moreover, the backward mapping
$F^{*}:\mathcal{R}\to\mathcal{D}$ recovers the set of canonical parameters
under which the joint probability distribution has node- and pair-wise
marginal distributions that agree with the mean parameters. In this setting,
$\mathcal{D}=\\{\\{(\theta_{i\bm{\cdot}},\theta_{jv\bm{\cdot}\bm{\cdot}})\\}_{i,j,v}|\theta_{i\bm{\cdot}}\in\mathbb{R}^{K},\theta_{jv\bm{\cdot}\bm{\cdot}}\in\mathbb{R}^{K\times
K},\theta_{jv\bm{\cdot}\bm{\cdot}}=\theta_{jv\bm{\cdot}\bm{\cdot}}^{\top},\theta_{jv\bm{\cdot}\bm{\cdot}}=\theta_{vj\bm{\cdot}\bm{\cdot}},\forall
i,j,v\in[p]\\}$, and $\mathcal{R}$ corresponds to the convex hull of all valid
marginal distributions (also known as marginal polytope). Unfortunately,
obtaining the exact and compact representation of the marginal polytope is NP-
hard [Roughgarden and Kearns, 2013], which leads to the intractability of
computing the backward mapping for time-varying DMRFs [Bresler et al., 2014,
Montanari, 2015]. To circumvent this difficulty, a common approach is to
resort to variational approximations of the backward mapping via tree-
reweighted entropy [Wainwright et al., 2003, Wainwright and Jordan, 2008],
which will be discussed in detail in §6.2.

###  2.1 Estimating the Canonical Parameters

In practice, the true mean parameters are rarely available and should be
replaced by their empirical analogs $\\{\widehat{\mu}_{t}\\}_{t=0}^{T}$, where
$\widehat{\mu}_{t}=\frac{1}{N_{t}}\sum_{i=1}^{N_{t}}\phi\left(x^{(i)}_{t}\right)$
and $\left\\{x^{(i)}_{t}\right\\}_{i=1}^{N_{t}}$ are the samples at time $t$
[Jordan et al., 1999]. Given the empirical mean parameters and the backward
mapping, one can estimate the canonical parameters from which the Markov
graphs can be extracted. However, the backward mapping may be ill-defined or
hard to derive explicitly, especially in the time-varying regime where the
number of available samples at each time may be extremely small [Wainwright
and Jordan, 2008, Zhou et al., 2010]. For example, for the time-varying GMRF,
the empirical moment parameter coincides with the sample covariance matrix
$\widehat{\Sigma_{t}}=\frac{1}{N_{t}}\sum_{i=1}^{N_{t}}x^{(i)}_{t}{x^{(i)}_{t}}^{\top}$.
Recalling the explicit form of the backward mapping, the precision matrix can
be estimated as $\widehat{\Theta_{t}}=\widehat{\Sigma_{t}}^{-1}$. In time-
varying settings, $N_{t}$ is often significantly smaller than $p$, which in
turn implies that the sample covariance matrix is rank deficient and,
therefore, not invertible. Under such circumstances, a consistent estimation
of the canonical parameters is only possible by leveraging additional prior
information on the structural or temporal properties of the MRF. For instance,
many real-world problems have sparsely- [Barabasi and Oltvai, 2004, Montanari
et al., 2012, Ravikumar et al., 2023] and/or smoothly-changing [Sporns et al.,
2004, Yu et al., 2008] structures. These structures can be promoted via
different regularizers.

For example, consider a sparsely-changing GMRF, where the time-varying
precision matrices are sparse and change sparsely over time. The corresponding
$\ell_{1}$-regularized MLE, also known as time-varying Graphical Lasso [Hallac
et al., 2017], is defined as:

$\displaystyle\\{\widehat{\Theta}_{t}\\}_{t=0}^{T}=\arg\\!\sup_{\\{{\Theta}_{t}\\}_{t=0}^{T}}$
$\displaystyle\
\sum_{t=0}^{T}\left(\langle\Theta_{t},\widehat{\Sigma}_{t}\rangle-\log\det(\Theta_{t})\right)+\gamma_{1}\sum_{t=0}^{T}\|\Theta_{t}\|_{\mathrm{off}}+\gamma_{2}\sum_{t=0}^{T}\|\Theta_{t}-\Theta_{t-1}\|_{1/1}$
(4a) s.t. $\displaystyle\ \ \Theta_{t}\succ 0\qquad\quad t=0,1,\dots,T$ (4b)

where $\left\lVert\Theta_{t}\right\rVert_{1}$ and
$\left\lVert\Theta_{t}-\Theta_{t-1}\right\rVert_{1/1}$ are convex relaxations
of the “$\ell_{0}$-norm” that promote sparsity in the precision matrices and
their differences. However, $\ell_{1}$-regularized MLE suffers from inferior
statistical properties, which is shown in our next example.

(a)

(b)

Figure 2: (a) The heatmap of the mismatch error. (b) The true and estimated
nonzero elements of the precision matrix.

##### Example (regularized MLE for time-varying GMRF).

Consider an instance of time-varying GMRF with $T=5$ and $n=25$, where
$\\{\Theta^{\star}_{t}\\}_{t=0}^{4}\in\mathbb{R}^{25\times 25}$ are randomly
generated symmetric and sparse matrices. At each time $t=0,\dots,4$, the true
precision matrix $\Theta^{\star}_{t}$ has exactly 30 off-diagonal elements
with value one (while maintaining its symmetry), and the remaining off-
diagonal entries are set to zero. Moreover, the diagonal entries
$\Theta^{\star}_{t;ii}$ are set to $1+\sum_{j\not=i}\Theta^{\star}_{t;ij}$ to
ensure $\Theta^{\star}_{t}\succ 0$ for every $t=0,\dots,4$. At every time, 5
nonzero off-diagonal elements are changed to zero, and 5 zero elements are set
to one. The sample covariance $\widehat{\Sigma}_{t}$ for every $t=0,\dots,4$
is obtained by collecting 500 samples from the Gaussian distribution with the
constructed precision matrices. Figure 2a illustrates a heatmap of the
mismatch error, i.e., the total number of false positives and false negatives
in the sparsity patterns of the true and estimated precision matrices and
their differences, for different values of the regularization coefficients. It
can be seen that after an exhaustive search over the regularization
coefficient space, the best achievable mismatch error is in the order of 50.
Note that simply predicting a fully sparse precision matrix has a mismatch
error of 190. In light of this, the estimated parameters reveal little
information about the true structure of the time-varying GMRF. Moreover,
Figure 2b shows the concatenation of the nonzero elements in the true
precision matrices (dashed red line), and their corresponding estimated values
via the regularized MLE (blue) and our proposed method (green). It can be seen
that, even when the sparsity pattern of the elements is correctly recovered,
the estimated nonzero entries via regularized MLE are “shrunk” toward zero due
to the shrinking effect of $\ell_{1}$ regularizer, thereby incurring a
substantial bias. In contrast, our proposed method circumvents this
undesirable bias by directly employing the $\ell_{0}$ regularizer instead.

##  3 Related Work

##### Time-varying MRF.

The problem of inferring time-varying MRFs can be traced back to Kalman
filters, where the goal is to predict a random signal evolving over time by
filtering out the observational noise [Kalman, 1960, Brown and Hwang, 1992].
More recent results have studied the non-asymptotic inference of time-varying
MRF with side information, such as sparsity and smoothness. In addition to the
time-varying Grahical Lasso introduced in (4), a recent line of works has
studied the inference of smoothly-changing MRFs [Kolar and Xing, 2011,
Greenewald et al., 2017, Zhou et al., 2010]. These methods rely on kernel
methods, where the empirical mean parameter at any given time is estimated as
a weighted average of the samples over time, and the weights are collected
from a predefined kernel. However, these methods do not leverage the prior
information about the sparsity of the parameter differences. With the goal of
addressing this deficiency, several works have studied the inference of
sparsely-changing MRF (also known as sparse differential networks) [Wang et
al., 2018, Zhao et al., 2014, Liu et al., 2017]. However, the main drawback of
these methods is that they only estimate the parameter differences, and their
theoretical guarantees are restricted to problems with two time steps ($T=1$).
Similarly, regression-based approaches have been proposed for change point
detection problems [Kolar and Xing, 2012, Roy et al., 2017] with MRFs and two
time periods, assuming the sparsity pattern of all entries of the precision
matrices change at the same time. In contrast, Wang and Kolar [2014] have
studied the inference of sparse MRFs given an index variable under the
assumption that the sparsity pattern is invariant, whereas Geng et al. [2019]
have assumed that the precision matrix is a linear function of the index
variable.

##### Sparsity-promoting optimization.

Unfortunately, solving statistical inference problems with sparsity is often
an NP-hard discrete optimization problem [Chickering, 1996, Natarajan, 1995],
and exact methods requiring exhaustive search [Lauritzen, 1996] are
impractical. Perhaps the most popular approaches are Fused Graphical Lasso and
Group Graphical Lasso [Danaher et al., 2014], which propose to relax
$\|\theta_{t}\|_{0}$ and $\|\theta_{t}-\theta_{t-1}\|_{0}$ with their
$\ell_{1}$-approximations. The $\ell_{1}$-norm proxy for sparsity is by now a
standard approach in inference with graphical models [Banerjee et al., 2008,
Friedman et al., 2008, Meinshausen and Bühlmann, 2006]. However, such
relaxation results in inferior statistical performance than its exact, non-
convex counterpart [Bertsimas et al., 2019]. Moreover, MLE-based estimators
with $\ell_{1}$ regularization are statistically consistent only under
stringent conditions, while in some cases, their nonconvex counterparts do not
require such strong restrictions [Lam and Fan, 2009].

Also closely related to our setting is the class of optimization problems of
the form

$\displaystyle\min_{\theta\in[\ell,u]^{p}}\;$
$\displaystyle\|\theta\|_{0}+\sum_{i=1}^{p}\sum_{j=i+1}^{p}g_{ij}(\theta_{i}-\theta_{j}),$
(5)

for given one-dimensional functions $g_{ij}:\mathbf{R}\to\mathbf{R}$,. If
functions $g_{ij}$ are convex, then problem (5) admits pseudo-polynomial time
algorithm [Ahuja et al., 2004, Bach, 2019]. Moreover, convex relaxations that
deliver near-optimal solutions for (5) were proposed for the special case of
convex quadratic $g$ functions [Atamtürk et al., 2018]. If, additionally, we
have $\ell=0$ and $u=\infty$, then problem (5) is, in fact, solvable in
strongly polynomial time [Atamtürk and Gómez, 2018]. As another special case,
if $\ell=-\infty$, $u=\infty$, and $g_{ij}(x)=0$ whenever $j>i+1$, then the
problem (5) with convex quadratic $g$ functions can be solved in quadratic
time [Liu et al., 2022]. On the other hand, problem (5) is much more
challenging for non-convex $g$: if $g(x)=\mathbbm{1}\left\\{x\neq 0\right\\}$,
as is the case for sparsely-changing MRFs, then problem (5) is NP-hard even if
the term $\|\theta\|_{0}$ is dropped from the objective [Hochbaum, 2001].
Nonetheless, as we show in this paper, problem (5) can be solved efficiently
in the context of time-varying MRFs, where $g_{ij}(x)=0$ whenever $j>i+1$.

##  4 Proposed Method

We propose the following class of constrained optimization problems:

$\displaystyle\\{\widehat{\theta}_{t}\\}_{t=0}^{T}=\arg\min_{\\{{\theta}_{t}\\}_{t=0}^{T}}$
$\displaystyle\gamma\\!\\!\\!\\!\underbrace{\sum_{t=0}^{T}\left\lVert\theta_{t}\right\rVert_{0}}_{\textbf{absolute
regularizer}}+(1-\gamma)\underbrace{\sum_{t=1}^{T}\left\lVert\theta_{t}-\theta_{t-1}\right\rVert^{q}_{q}}_{\textbf{temporal
regularizer}},$ (ProxGL) $\displaystyle\mathrm{s.t.}$
$\displaystyle\underbrace{\left\lVert\theta_{t}-\widetilde{F}^{*}(\widehat{\mu}_{t})\right\rVert_{\infty}}_{\textbf{backward
mapping deviation}}\leq\lambda_{t}\qquad\forall t=0,\dots,T,$

where $q\geq 0$ in the temporal $\ell_{q}$-regularizer is chosen to match the
temporal changes in the canonical parameters. For example, $q=0$ is a natural
choice for sparsely-changing MRF where only a few elements of the canonical
parameters change in consecutive times, while $q=2$ is suitable for smoothly-
changing MRFs where canonical parameters change smoothly over time. The
parameter $0\leq\gamma\leq 1$ is a regularization parameter that is used to
balance the sparsity level and the temporal changes. Indeed, a larger value
for $\gamma$ would lead to sparser parameters, whereas a smaller $\gamma$
would result in smaller temporal changes in the parameters. In the above
optimization, the backward mapping deviation controls the deviation of the
canonical parameters from an approximate backward mapping
$\widetilde{F}^{*}(\widehat{\mu}_{t})$, which we assume can be obtained
directly from the data. In §6, we provide a detailed discussion on how to
obtain such approximate backward mapping for different classes of time-varying
MRFs. Despite its nonconvexity, our next theorem shows that ProxGL can be
solved efficiently. Recall that $T$ is the number of time steps, and $p$ is
the dimension of the canonical parameter at each time.

######  Theorem 1 (Efficient algorithm).

For all values of the regularization parameter $\gamma\in[0,1]$ and any
$q\in\\{0\\}\cup[1,\infty)$, the entire solution path of ProxGL can be
obtained in at most $\mathcal{O}(pT^{3})$ time on a single thread.

The proof of Theorem 1 is provided in §5 together with the detailed
implementation of our proposed algorithm. To the best of our knowledge,
Theorem 1 is the first result showing that the entire solution path of time-
varying MRFs can be estimated in strongly polynomial time. The availability of
the solution path facilitates the automated fine-tuning of the sparsity level
and/or smoothness of the obtained canonical parameters via cross-validation or
by solving a bilevel optimization, where the estimated canonical parameters on
the training set are parameterized as functions of $\gamma$ and optimized on
the validation set [Sinha et al., 2014, MacKay et al., 2019]. Our algorithm
for solving ProxGL relies on two key properties of our proposed formulation:
First, ProxGL is decomposable into $p$ independent subproblems over different
coordinates of $\\{\theta_{t}\\}_{t=0}^{T}$. This guarantees that the
complexity of our algorithm depends linearly on the dimension $p$. Moreover,
we propose a dynamic programming approach that can solve each subproblem in
$O(T^{3})$ for all values of $\gamma$.

Given the efficient solvability of ProxGL, our next goal is to show that the
obtained solution indeed enjoys a small estimation error. To this goal, we
formally define two classes of temporal structures for the canonical
parameters, namely sparsely-changing MRFs and smoothly-changing MRFs.

######  Definition 1 (Sparsely-changing MRF).

A time-varying MRF is sparsely-changing with parameter $D_{0}\geq 0$ if
$\left\lVert\theta^{\star}_{t}-\theta^{\star}_{t-1}\right\rVert_{0}\leq D_{0}$
for every $t=1,\dots,T$.

Intuitively, we say that “a time-varying MRF is sparsely-changing” if it is
sparsely-changing with a parameter that satisfies $D_{0}\ll p$. We note that,
if $\left\lVert\theta^{\star}_{t}\right\rVert_{0}\leq S_{0}$ for every $t$,
then $D_{0}\leq 2S_{0}$. In other words, time-varying MRFs with sparse
parameters are also sparsely-changing with parameter $2S_{0}$.

######  Definition 2 (Smoothly-changing MRF).

A time-varying MRF is smoothly-changing with parameters $(q,D)$ if
$\underset{1\leq i\leq
n}{\max}\left\\{\sum_{t=1}^{T}|\theta^{\star}_{t;i}-\theta^{\star}_{t-1;i}|^{q}\right\\}\leq
D$.

Similarly, We say “a time-varying MRF is smoothly-changing” if it is smoothly-
changing with $q\geq 1$ and $D\ll T$. Indeed, a time-varying MRF can be
simultaneously sparsely-changing and smoothly-changing. In such settings, only
a few elements of the canonical parameters change over time and the changes
are smooth. However, to streamline the presentation, we study these two
scenarios separately.

Let $\mathcal{S}_{t}$ denote the set of indices corresponding to the nonzero
elements of $\theta_{t}^{\star}$ for every $t=0,\dots,T$. Similarly, let
$\mathcal{D}_{t}$ denote the set of indices corresponding to the nonzero
elements of $\theta_{t}^{\star}-\theta_{t-1}^{\star}$ for every $t=1,\dots,T$.
Our next theorem provides conditions under which an optimal solution of ProxGL
enjoys a small estimation error and correct sparsity pattern.

######  Theorem 2 (Estimation error and sparsistency).

Suppose that the parameter $\lambda_{t}$ in ProxGL satisfies

$\displaystyle\left\lVert\theta_{t}^{\star}-\widetilde{F}^{*}(\widehat{\mu}_{t})\right\rVert_{\infty}\leq\lambda_{t},\quad\forall
t=0,\dots,T.$ (A1)

Then the following statements hold:

  * •

Estimation error. We have

$\displaystyle\left\lVert\widehat{\theta}_{t}-\theta^{\star}_{t}\right\rVert_{\infty}\leq
2\lambda_{t},\quad\forall t=0,\dots,T.$ (6)

  * •

Sparsistency for smoothly-changing MRF. Suppose that the time-varying MRF is
smoothly-changing with parameters $(q,D)$ for some $q\geq 1,D\geq 0$.
Additionally, suppose that $0<\gamma<1/(1+D)$ and
$2\lambda_{t}\leq\min_{i\in\mathcal{S}_{t}}|\theta_{t;i}^{\star}|$. Then, the
optimal solution of ProxGL with temporal $\ell_{q}$-regularizer satisfies:

$\displaystyle\mathrm{supp}(\widehat{\theta}_{t})=\mathrm{supp}({\theta}^{\star}_{t}),\qquad\forall
t=0,\dots,T.$ (7)

  * •

Sparsistency for sparsely-changing MRF. Suppose that the time-varying MRF is
sparsely-changing with parameter $D_{0}\geq 0$. Moreover, suppose that
$2\lambda_{t}\leq\min_{i\in\mathcal{S}_{t}}|\theta_{t;i}^{\star}|$ and
$2\lambda_{t}+2\lambda_{t-1}\leq\min_{i\in\mathcal{D}_{t}}|\theta_{t;i}^{\star}-\theta_{t-1;i}^{\star}|$.
Then, with any choice of $0<\gamma<1$, the optimal solution of ProxGL with
temporal $\ell_{0}$-regularizer satisfies:

$\displaystyle\mathrm{supp}(\widehat{\theta}_{t})=\mathrm{supp}({\theta}^{\star}_{t}),$
$\displaystyle\qquad\forall t=0,\dots,T,$ (8)
$\displaystyle\mathrm{supp}(\widehat{\theta}_{t}-\widehat{\theta}_{t-1})=\mathrm{supp}({\theta}^{\star}_{t}-{\theta}^{\star}_{t-1}),$
$\displaystyle\qquad\forall t=1,\dots,T.$

The above theorem shows that if the true canonical parameters remain feasible
for ProxGL with sufficiently small $\lambda_{t}$, then the estimated
parameters have small errors and correct sparsity patterns. Moreover, if the
time-varying MRF is sparsely-changing, then the solution of ProxGL with a
temporal $\ell_{0}$-regularizer recovers the correct sparsity pattern of the
individual parameters, as well as their temporal changes. In §6, we use the
above theorem to characterize the sample complexity of our proposed framework
for two classes of time-varying MRFs, namely GMRFs and DMRFs. Next, we present
the proof of Theorem 2.

Proof of Theorem 2. We first present the proof of estimation error bound (6).

##### Estimation error bound.

The solution of ProxGL satisfies

$\displaystyle\left\lVert\widehat{\theta}_{t}-\theta^{\star}_{t}\right\rVert_{\infty}\leq\left\lVert\widehat{\theta}_{t}-\widetilde{F}^{*}(\widehat{\mu}_{t})\right\rVert_{\infty}+\left\lVert\theta^{\star}_{t}-\widetilde{F}^{*}(\widehat{\mu}_{t})\right\rVert_{\infty}\leq
2\lambda_{t},$

where in the first inequality, we used the triangle inequality, and in the
second inequality, we used the assumption that the true canonical parameter
$\\{\theta^{\star}_{t}\\}_{t=0}^{T}$ is feasible for ProxGL.

##### Sparsistency.

Next, we show the sparsistency of the estimated parameters for any
$q\in\\{0\\}\cup[1,\infty)$. First, we prove
$\mathrm{supp}(\widehat{\theta}_{t})=\mathrm{supp}(\theta^{\star}_{t})$ for
any $q\geq 1$. To this goal, we show
$\mathrm{supp}(\widehat{\theta}_{t})\subseteq\mathrm{supp}(\theta^{\star}_{t})$
and
$\mathrm{supp}(\theta^{\star}_{t})\subseteq\mathrm{supp}(\widehat{\theta}_{t})$.
For any $i\in\mathcal{S}_{t}$, we have

$\displaystyle\left|\widehat{\theta}_{t;i}\right|$
$\displaystyle=\left|\widehat{\theta}_{t;i}-\theta_{t;i}^{\star}+\theta_{t;i}^{\star}\right|$
$\displaystyle\geq\left|\theta_{t;i}^{\star}\right|-\left|\widehat{\theta}_{t;i}-\theta_{t;i}^{\star}\right|$
$\displaystyle\geq\left|\theta_{t;i}^{\star}\right|-2\lambda_{t}$
$\displaystyle>0,$

where in the second inequality we used the upper bound on the estimation
error, and in the last inequality we used the assumption
$2\lambda_{t}\leq\min_{i\in\mathcal{S}_{t}}|\theta^{\star}_{t;i}|$. This
implies that
$\mathrm{supp}(\theta^{\star}_{t})\subseteq\mathrm{supp}(\widehat{\theta}_{t})$.
To prove
$\mathrm{supp}(\widehat{\theta}_{t})\subseteq\mathrm{supp}(\theta^{\star}_{t})$,
we first rewrite ProxGL as follows:

$\displaystyle\\{\widehat{\theta}_{t}\\}_{t=0}^{T}=\arg\min_{\\{{\theta}_{t}\\}_{t=0}^{T}}$
$\displaystyle\gamma\sum_{t=0}^{T}\sum_{i=1}^{n}\mathbb{I}[\theta_{t;i}\neq
0]+(1-\gamma)\sum_{t=1}^{T}\sum_{i=1}^{n}\left|\theta_{t;i}-\theta_{t-1;i}\right|^{q}$
(9) $\displaystyle\mathrm{s.t.}$
$\displaystyle\left|\theta_{t;i}-\left[\widetilde{F}^{*}(\widehat{\mu}_{t})\right]_{i}\right|\leq\lambda_{t}\qquad\forall
t=0,\dots,T,\;i=1,\dots,p.$

We make the following observation. First, due to the use of the infinity norm
for the backward mapping deviation, problem (9) decomposes into $p$
independent subproblems, one for each coordinate of $\theta_{t}$. On the other
hand, due to the optimality of $\\{\widehat{\theta}_{t}\\}_{t=0}^{T}$ and the
feasibility of $\\{\theta^{\star}_{t}\\}_{t=0}^{T}$, we have for every $1\leq
i\leq p$

$\displaystyle(1\\!-\\!\gamma)\sum_{t=0}^{T}\mathbb{I}\\{\widehat{\theta}_{t;i}\not=0\\}\\!+\\!\gamma\sum_{t=1}^{T}{\left|\widehat{\theta}_{t;i}\\!-\\!\widehat{\theta}_{t-1;i}\right|^{q}}\\!\leq\\!(1-\gamma)\sum_{t=0}^{T}\mathbb{I}\\{\theta^{\star}_{t;i}\not=0\\}\\!+\\!\gamma\sum_{t=1}^{T}{\left|{\theta}^{\star}_{t;i}-{\theta}^{\star}_{t-1;i}\right|^{q}}$
$\displaystyle\implies$
$\displaystyle(1-\gamma)\sum_{t=0}^{T}\mathbb{I}\\{\widehat{\theta}_{t;i}\not=0\\}\leq(1-\gamma)\sum_{t=0}^{T}\mathbb{I}\\{\theta^{\star}_{t;i}\not=0\\}+\gamma
D$ $\displaystyle\implies$
$\displaystyle(1\\!-\\!\gamma)\\!\sum_{t=0}^{T}\mathbb{I}\\{\widehat{\theta}_{t;i}\\!\not=\\!0\\}\mathbb{I}\\{i\\!\in\\!\mathcal{S}_{t}\\}\\!+\\!(1\\!-\\!\gamma)\\!\sum_{t=0}^{T}\mathbb{I}\\{\widehat{\theta}_{t;i}\\!\not=\\!0\\}\mathbb{I}\\{i\\!\in\\![n]\backslash\mathcal{S}_{t}\\}\\!\leq\\!(1\\!-\\!\gamma)\\!\sum_{t=0}^{T}\\!\left|{\theta}^{\star}_{t;i}\right|_{0}\\!+\\!\gamma
D_{q}$ $\displaystyle\implies$
$\displaystyle(1-\gamma)\sum_{t=0}^{T}\mathbb{I}\\{\widehat{\theta}_{t;i}\not=0\\}\mathbb{I}\\{i\in[n]\backslash\mathcal{S}_{t}\\}\leq\gamma
D$ $\displaystyle\implies$
$\displaystyle\sum_{t=0}^{T}\mathbb{I}\\{\widehat{\theta}_{t;i}\not=0\\}\mathbb{I}\\{i\in[n]\backslash\mathcal{S}_{t}\\}<1$

where the second inequality follows from the assumption
$\sum_{t=1}^{T}\sum_{i\in[p]}|\theta_{t;i}^{*}-\theta_{t-1;i}^{*}|^{q}\leq D$.
Moreover, the fourth inequality follows from
$\mathrm{supp}(\theta^{\star}_{t})\subseteq\mathrm{supp}(\widehat{\theta}_{t})$.
The last inequality follows from the assumption $0<\gamma<1/(1+D)$. This
implies that
$\mathbb{I}\\{\widehat{\theta}_{t;i}\not=0\\}\mathbb{I}\\{i\in[n]\backslash\mathcal{S}_{t}\\}=0$
for every $t$ and $i$, and hence, we have
$\mathrm{supp}(\widehat{\theta}_{t})\subseteq\mathrm{supp}(\theta^{\star}_{t})$.

Finally, we show that for sparsely-changing MRFs, ProxGL with temporal
$\ell_{0}$-regularizer satisfies
$\mathrm{supp}(\widehat{\theta}_{t})=\mathrm{supp}({\theta}^{\star}_{t})$ and
$\mathrm{supp}(\widehat{\theta}_{t}-\widehat{\theta}_{t-1})=\mathrm{supp}({\theta}^{\star}_{t}-{\theta}^{\star}_{t-1})$.
An argument identical to $q\geq 1$ can be invoked to show that
$\mathrm{supp}(\theta^{\star}_{t})\subseteq\mathrm{supp}(\widehat{\theta}_{t})$.
Similarly, given any $i\in\mathcal{D}_{t}$, one can write

$\displaystyle\left|\widehat{\theta}_{t;i}-\widehat{\theta}_{t-1;i}\right|$
$\displaystyle=\left|\widehat{\theta}_{t;i}-\theta_{t;i}^{*}+\theta_{t;i}^{*}-\theta_{t-1;i}^{*}+\theta_{t-1;i}^{*}-\widehat{\theta}_{t-1;i}\right|$
$\displaystyle\geq\left|\theta_{t;i}^{*}-\theta_{t-1;i}^{*}\right|-\left|\widehat{\theta}_{t;i}-\theta_{t;i}^{*}\right|-\left|\widehat{\theta}_{t-1;i}-\theta_{t-1;i}^{*}\right|$
$\displaystyle\geq\left|\theta_{t;i}^{*}-\theta_{t-1;i}^{*}\right|-2\lambda_{t}-2\lambda_{t-1}$
$\displaystyle>0,$

where the last inequality is due to our assumption
$2\lambda_{t}+2\lambda_{t-1}\leq\min_{i\in\mathcal{D}_{t}}|\theta_{t;i}^{\star}-\theta_{t-1;i}^{\star}|$.
This implies that
$\mathrm{supp}({\theta}^{\star}_{t}-{\theta}^{\star}_{t-1})\subseteq\mathrm{supp}(\widehat{\theta}_{t}-\widehat{\theta}_{t-1})$.
On the other hand, due to the optimality of
$\\{\widehat{\theta}_{t}\\}_{t=0}^{T}$ and the feasibility of
$\\{\theta^{\star}_{t}\\}_{t=0}^{T}$, we have

$\displaystyle(1-\gamma)\sum_{t=0}^{T}\|\widehat{\theta}_{t}\|_{0}+\gamma\sum_{t=1}^{T}\|\widehat{\theta}_{t}-\widehat{\theta}_{t-1}\|_{0}\leq(1-\gamma)\sum_{t=0}^{T}\|\theta_{t}^{*}\|_{0}+\gamma\sum_{t=1}^{T}\|\theta_{t}^{*}-\theta_{t-1}^{*}\|_{0}$
$\displaystyle\implies$
$\displaystyle(1\\!-\\!\gamma)\sum_{t=0}^{T}\\!\left(\sum_{i\in[n]\backslash\mathcal{S}_{t}}\mathbb{I}\\{\widehat{\theta}_{t;i}\not=0\\}\\!+\\!\sum_{i\in\mathcal{S}_{t}}\mathbb{I}\\{\widehat{\theta}_{t;i}\not=0\\}\right)$
$\displaystyle\\!+\\!\gamma\sum_{t=1}^{T}\left(\sum_{i\in[n]\backslash\mathcal{D}_{t}}\\!\mathbb{I}\\{\widehat{\theta}_{t;i}-\widehat{\theta}_{t-1;i}\not=0\\}\\!+\\!\sum_{i\in\mathcal{D}_{t}}\\!\mathbb{I}\\{\widehat{\theta}_{t;i}-\widehat{\theta}_{t-1;i}\not=0\\}\right)$
$\displaystyle\leq(1-\gamma)\sum_{t=0}^{T}\sum_{i\in\mathcal{S}_{t}}\mathbb{I}\\{\theta^{\star}_{t;i}\not=0\\}+\gamma\sum_{t=1}^{T}\sum_{i\in\mathcal{D}_{t}}\mathbb{I}\\{\theta^{\star}_{t;i}-\theta^{\star}_{t-1;i}\not=0\\}$
$\displaystyle\implies$
$\displaystyle(1-\gamma)\sum_{t=0}^{T}\sum_{i\in[n]\backslash\mathcal{S}_{t}}\mathbb{I}\\{\widehat{\theta}_{t;i}\not=0\\}+\gamma\sum_{t=1}^{T}\sum_{i\in[n]\backslash\mathcal{D}_{t}}\mathbb{I}\\{\widehat{\theta}_{t;i}-\widehat{\theta}_{t-1;i}\not=0\\}\leq
0$

where the last inequality follows from
$\mathrm{supp}(\theta^{*}_{t})\subseteq\mathrm{supp}(\widehat{\theta}_{t})$
and
$\mathrm{supp}(\theta_{t}^{*}-\theta_{t-1}^{*})\subseteq\mathrm{supp}(\widehat{\theta}_{t}-\widehat{\theta}_{t-1})$.
Due to $0<\gamma<1$, the above inequality implies that
$\widehat{\theta}_{t;i}=0$ for every $t=0,\dots,T$ and
$i\in[n]\backslash\mathcal{S}_{t}$, and
$\widehat{\theta}_{t;i}-\widehat{\theta}_{t-1;i}=0$ for every $t=1,\dots,T$
and $i\in[n]\backslash\mathcal{D}_{t}$. Therefore, we have
$\mathrm{supp}(\widehat{\theta}_{t})\subseteq\mathrm{supp}(\theta^{*}_{t})$
and
$\mathrm{supp}(\widehat{\theta}_{t}-\widehat{\theta}_{t-1})\subseteq\mathrm{supp}(\theta_{t}^{*}-\theta_{t-1}^{*})$.
This completes the proof.$\hfill\square$

##  5 Efficient Algorithm

In this section, we describe the proposed algorithm for solving ProxGL,
assuming that the approximate backward mapping
$\widetilde{F}^{*}(\widehat{\mu}_{t})$ is known. In the next section, we study
efficient ways for obtaining $\widetilde{F}^{*}(\widehat{\mu}_{t})$. Recall
that, due to the use of the infinity norm for the backward mapping deviation,
problem (9) decomposes into $p$ independent subproblems, one for each
coordinate of $\theta_{t}$. Moreover, if $\gamma=1$, then problem (9) can be
solved trivially, by setting $\theta_{t;i}=0$ whenever that solution is
feasible and setting $\theta_{t;i}=\left[F^{*}(\widehat{\mu}_{t})\right]_{i}$
otherwise. Thus, we assume without loss of generality that $p=1$ (so we omit
subscripts $i$) and $\gamma<1$, and focus on solving the problem

$\displaystyle\min_{\\{{\theta}_{t}\\}_{t=0}^{T}}$
$\displaystyle\bar{\gamma}\sum_{t=0}^{T}\mathbb{I}[\theta_{t}\neq
0]+\sum_{t=1}^{T}\left|\theta_{t}-\theta_{t-1}\right|^{q}$ (10)
$\displaystyle\mathrm{s.t.}$ $\displaystyle\ell_{t}\leq\theta_{t}\leq
u_{t}\qquad\forall t=0,\dots,T,$

where $\bar{\gamma}=\gamma/(1-\gamma)$,
$\ell_{t}=-\lambda_{t}+\widetilde{F}^{*}(\widehat{\mu}_{t})$ and
$u_{t}=\lambda_{t}+\widetilde{F}^{*}(\widehat{\mu}_{t})$. First, in §5.1, we
discuss how to solve (10) for a given value of the parameter $\bar{\gamma}$.
We note that Fattahi and Gomez [2021] proposed a method for solving (10) with
$q=0$ by casting it as a shortest path problem, which is a special case of our
proposed dynamic programming approach. Then, in §5.2, we give an algorithm to
solve problem (10) for all values of $\bar{\gamma}$. Finally, in §5.3 we
discuss further refinements of the algorithm. Overall, the main result of this
section is an algorithm that solves (10) in strongly polynomial time,
summarized in Proposition 1 below.

######  Proposition 1.

The entire solution path of (10) for all values of
$0\leq\bar{\gamma}\leq\infty$ and any fixed $q\in\\{0\\}\cup[1,\infty)$ can be
computed in $\mathcal{O}(T^{3})$.

###  5.1 Fixed parameter solution

A key insight to solve (10) is that if $\theta_{t}=0$ for any index $t$, then
(10) decomposes into two independent problems. This key observation gives rise
to a dynamic programming approach for solving (10). Given $0\leq a<b\leq T+1$,
define function $v(a,b)$ to be the best objective value of (10) from time
period $a$ to $b-1$ assuming that $\theta_{a-1}=\theta_{b}=0$ (if such
variables exist), that is,

$\displaystyle v(a,b)=\min_{\\{{\theta}_{t}\\}_{t=a}^{b-1}}$
$\displaystyle\bar{\gamma}\sum_{t=a}^{b-1}\mathbb{I}[\theta_{t}\neq
0]\\!+\\!\\!\sum_{t=a+1}^{b-1}\left|\theta_{t}\\!-\\!\theta_{t-1}\right|^{q}\\!+\\!\mathbb{I}[a>0](\theta_{a}-0)^{q}\\!+\\!\mathbb{I}[b<T+1]|\theta_{b-1}-0|^{q}$
(11) $\displaystyle\mathrm{s.t.}$ $\displaystyle\ell_{t}\leq\theta_{t}\leq
u_{t}\qquad\forall t=a,\dots,b-1.$

Clearly, the optimal objective value of problem (10) is $v(0,T+1)$. Moreover,
let $f(a,b)$ denote the best objective value of (10) from indexes $a$ to $b-1$
while ignoring the absolute regularizer terms, that is,

$\displaystyle f(a,b)=\min_{\\{{\theta}_{t}\\}_{t=a}^{b-1}}$
$\displaystyle\sum_{t=a+1}^{b-1}\left|\theta_{t}-\theta_{t-1}\right|^{q}+\mathbb{I}[a>0]\left|\theta_{a}-0\right|^{q}+\mathbb{I}[b<T+1]\left|\theta_{b-1}-0\right|^{q}$
(12) $\displaystyle\mathrm{s.t.}$ $\displaystyle\ell_{t}\leq\theta_{t}\leq
u_{t}\qquad\forall t=a,\dots,b-1.$

By convention, we also let $f(a,a)=v(a,a)=0$ for all $a=0,\dots,T+1$. In §5.3,
we discuss the complexity of solving (12).

Functions $v$ and $f$ are related through the following recursion:

######  Lemma 1.

For any $0\leq a<b\leq T+1$, we have

$v(a,b)=\min\Big{\\{}\underbrace{\min_{t\in\\{a,\dots,b-1\\}:\;0\in[\ell_{t},u_{t}]}v(a,t)+v(t+1,b)}_{\text{Cost
of setting }\theta_{t}=0},\;\underbrace{f(a,b)+\bar{\gamma}(b-a)}_{\text{Cost
of setting }\theta_{t}\neq 0,\;a\leq t\leq b-1}\Big{\\}}.$ (13)

###### Proof.

The optimal solution corresponding to $v(a,b)$ either satisfies
$\theta_{t}\not=0$ for all $a\leq t\leq b-1$, or $\theta_{t}=0$ for some
$a\leq t\leq b-1$. If $\theta_{t}\not=0$ for all $a\leq t\leq b-1$, then we
have $v(a,b)=f(a,b)+\bar{\gamma}\sum_{t=a}^{b-1}\mathbb{I}[\theta_{t}\neq
0]=f(a,b)+\bar{\gamma}(b-a)$. Otherwise, the problem (11) decomposes at the
break-point $t$ with $\theta_{t}=0$ which results in $v(a,b)=v(a,t)+v(t+1,b)$.
Therefore, $v(a,b)$ can be obtained by taking the minimum of
$f(a,b)+\bar{\gamma}(b-a)$ and $v(a,t)+v(t+1,b)$ for the best choice of the
break-point $t$. ∎

In our subsequent arguments, $EO$ stands for “evaluation oracle” and denotes
the complexity of solving (12). Moreover, $EO_{2}$ stands for the complexity
of solving (12) for all combinations of parameters $0\leq a<b\leq T+1$.
Clearly, $EO_{2}=\mathcal{O}(T^{2}EO)$, by simply calling the evaluation
oracle $T^{2}$ times, although more efficient methods may exist.

######  Lemma 2.

For a given value of $\bar{\gamma}$, (10) can be solved with complexity
$\mathcal{O}\big{(}EO_{2}+T^{3}\big{)}$.

###### Proof.

From the recursion (13), it follows that if all values of $v(a,t)$ and
$v(t+1,b)$ for all $t\in\\{a,\dots,b-1\\}$ are available, and $f(a,b)$ is
available as well, then computing $v(a,b)$ reduces to performing
$\mathcal{O}(T)$ comparisons. Since $v$ needs to be computed
$\mathcal{O}(T^{2})$ for all $a<b$, to ensure $v(a,t)$ and $v(t+1,b)$ are
available, we obtain the complexity of $\mathcal{O}(T^{3})$. The term
$\mathcal{O}(EO_{2})$ corresponds to computing function $f$ as well. ∎

###  5.2 Parametric scheme for solution path

We now show that (10) can be solved _for all values of $\bar{\gamma}$_ with
the same theoretical complexity as the dynamic programming algorithm discussed
in §5.1. Observe that the optimal solution $\\{\hat{\theta}_{t}\\}_{t=0}^{T}$
of (10) is also optimal for the optimization problem

$\displaystyle\min_{\\{{\theta}_{t}\\}_{t=0}^{T}}$
$\displaystyle\sum_{t=1}^{T}\left|\theta_{t}-\theta_{t-1}\right|^{q}$ (14)
$\displaystyle\mathrm{s.t.}$
$\displaystyle\sum_{t=0}^{T}\mathbb{I}[\theta_{t}\neq 0]\leq k$
$\displaystyle\ell_{t}\leq\theta_{t}\leq u_{t}\qquad\forall t=0,\dots,T$

with $k=\sum_{t=0}^{T}\mathbb{I}[\hat{\theta}_{t}\neq 0]$. Thus, to solve (10)
for all values of $\bar{\gamma}$, one can instead solve (14) for all values of
$k\in\\{0,\dots,T\\}$ and select the solution with best objective value.
Observe that this approach produces optimal solutions for all values of the
regularization parameter $\bar{\gamma}$. Also note that since problems (10)
and (14) are non-convex, the regularized and constrained problems are not
equivalent: some optimal solutions of (14) for a given value of $k$ may not be
optimal for (10) for any value of $\bar{\gamma}$.

Given $0\leq b\leq T+1$ and $k\in\mathbb{Z}$, define function $\nu(b,k)$ to be
the best objective value of (14) from indexes $0$ to $b-1$ assuming
$\theta_{b}=0$ (if $b\leq T$), that is,

$\displaystyle\nu(b,k)=\min_{\\{{\theta}_{t}\\}_{t=0}^{b-1}}$
$\displaystyle\sum_{t=1}^{b-1}\left|\theta_{t}-\theta_{t-1}\right|^{q}+\mathbb{I}[b<T+1]\left|\theta_{b-1}-0\right|^{q}$
$\displaystyle\mathrm{s.t.}\;$
$\displaystyle\sum_{t=0}^{b-1}\mathbb{I}[\theta_{t}\neq 0]\leq k$
$\displaystyle\ell_{t}\leq\theta_{t}\leq u_{t}\qquad\forall t=a,\dots,b-1,$

where we adopt the convention that $\nu(b,k)=\infty$ if the optimization is
not feasible (e.g., if $k<0$), and $\nu(0,k)=0$ for all $k\in\mathbb{Z}_{+}$.
Our goal is to compute $\nu(T+1,k)$ for all values of $k\in\\{0,\dots,T+1\\}$.
Note that for all $k\geq b$, we have that $\nu(b,k)=f(0,b)$, where $f$ is
given by (12). For $k<b$, functions $\nu$ and $f$ are related through the
following recursion:

######  Lemma 3.

For any $0\leq k<b\leq T+1$, we have

$\nu(b,k)=\min_{t\in\\{0,\dots,b-1\\}:\;0\in[\ell_{t},u_{t}]}\nu(t,k-(b-t-1))+f(t+1,b).$
(15)

###### Proof.

The optimal solution corresponding to $\nu(b,k)$ satisfies $\theta_{t}=0$ for
some $0\leq t\leq b-1$ (due to the assumption that $k<b$). Letting $\bar{t}$
be the largest such index, and noting that the problem decomposes at the
breakpoint $\bar{t}$ with $\theta_{\bar{t}}=0$, it follows that

$\nu(b,k)=\underbrace{\nu(\bar{t},k-(b-\bar{t}-1))}_{\text{Cost from 0 to
$\bar{t}-1$ with remaining budget}}+\underbrace{f(\bar{t}+1,b)}_{\text{Cost
from $\bar{t}+1$ to $b-1$ with $\theta_{t}\neq 0$.}}.$

Therefore, $\nu(b,k)$ can be obtained by selecting the break-point $\bar{t}$
with least cost. ∎

######  Lemma 4.

The problem (14) can be solved for all values of $\bar{\gamma}$ with
complexity $\mathcal{O}(T^{3}+EO_{2})$, where $EO_{2}$ is the complexity of
solving (12) for all combinations of parameters $a<b$.

###### Proof.

Computing all values of $\nu(b,k)$ amounts to computing $f(a,b)$ for all
$\mathcal{O}(T^{2})$ combinations of parameters $a<b$ and doing
$\mathcal{O}(T^{3})$ comparisons, corresponding to $\mathcal{O}(T^{2})$
choices of $b$ and $k$, with each computation of $v(b,k)$ requiring
$\mathcal{O}(T)$ comparisons. ∎

Observe that if solving (12) for all combinations of $a<b$ can be done in
$\mathcal{O}(T^{3})$, then the overall complexity of the method is cubic in
$T$. We now show that this is indeed the case.

###  5.3 Solving the subproblems

Solving (10), either for a fixed value of $\bar{\gamma}$ or parametrically for
all values of the regularization parameter, requires solving problem (12) for
all $0\leq a<b\leq T+1$. We now discuss how to accomplish this task
efficiently, depending on the value of $q$. In particular, we show that (12)
for all $0\leq a<b\leq T+1$ can be solved in $O(T^{2})$ for $q=0$, and in
$O(T^{3})$ for any $q\geq 1$.

####  5.3.1 Case $q=0$

We now provide an efficient algorithm (Algorithm 1) for solving (12) for the
case $q=0$, which we restate for convenience:

$\displaystyle f(a,b)=\min_{\\{{\theta}_{t}\\}_{t=a}^{b-1}}$
$\displaystyle\sum_{t=a+1}^{b-1}\mathbb{I}[\theta_{t}\neq\theta_{t-1}]+\mathbb{I}[a>0]\mathbb{I}[\theta_{a}\neq
0]+\mathbb{I}[b<T+1]\mathbb{I}[\theta_{b-1}\neq 0]$ (16)
$\displaystyle\mathrm{s.t.}$ $\displaystyle\ell_{t}\leq\theta_{t}\leq
u_{t}\qquad\forall t=a,\dots,b-1.$

For simplicity, we can introduce a dummy variable $\theta_{a-1}$ with bounds
$\ell_{a-1}=u_{a-1}=0$ if $a>0$, and bounds $-\ell_{a-1}=u_{a-1}=\infty$ if
$a=0$. Similarly, we can introduce a variable $\theta_{b}$ with bounds
$\ell_{b}=u_{b}=0$ if $b<T+1$, and bounds $-\ell_{a-1}=u_{a-1}=\infty$ if
$b=T+1$. With these additional variables, we can rewrite (16) as

$\displaystyle f(a,b)=\min_{\\{{\theta}_{t}\\}_{t=a-1}^{b}}$
$\displaystyle\sum_{t=a}^{b}\mathbb{I}[\theta_{t}\neq\theta_{t-1}]$ (17)
$\displaystyle\mathrm{s.t.}$ $\displaystyle\ell_{t}\leq\theta_{t}\leq
u_{t}\qquad\forall t=a-1,\dots,b.$

The next lemma provides a method to compute lower bounds on the objective
value of (17).

######  Lemma 5.

Given any indexes $a\leq\tau_{1}<\tau_{2}\leq b$, if $\max_{\tau_{1}\leq
t\leq\tau_{2}}\ell_{t}>\min_{\tau_{1}\leq t\leq\tau_{2}}u_{t}$, then

$\sum_{t=\tau_{1}+1}^{\tau_{2}}\mathbb{I}[\theta_{t}\neq\theta_{t-1}]\geq 1,$
(18)

for any feasible solution $\theta_{t}$ of (17).

###### Proof.

The left hand side of (18) is zero if and only if
$\theta_{\tau_{1}}=\theta_{\tau_{1}+1}=\dots=\theta_{\tau_{2}}$. However, the
condition $\max_{\tau_{1}\leq t\leq\tau_{2}}\ell_{t}>\min_{\tau_{1}\leq
t\leq\tau_{2}}u_{t}$ ensures that such solutions are not feasible. ∎

Algorithm 1 below solves (17) to optimality. At a high level, starting from
time period $\tau_{1}=a$, the algorithm greedily finds the largest time period
$\tau_{2}$ such that setting
$\theta_{\tau_{1}}=\theta_{\tau_{1}+1}=\dots=\theta_{\tau_{2}}$ is feasible.
In the algorithm, we maintain a set $\Gamma$ corresponding to the first and
last time periods, as well as breakpoints $\tau$ such that
$\theta_{\tau_{1}}\neq\theta_{\tau_{1}-1}$. We use $\texttt{last}(\Gamma)$ to
denote the last element added to set $\Gamma$.

Algorithm 1 $\ell_{0}$ greedy

1:Input: $\hat{\theta}\in\mathbf{R}^{b-a+2}$: array to store the optimal
solution

2:Output: The optimal objective value $f(a,b)$

3:$\bar{f}\leftarrow 0$

4:$\bar{\ell}\leftarrow\ell_{a-1}$, $\bar{u}\leftarrow u_{a-1}$
$\triangleright$ Current bounds on constant interval

5:$\Gamma\leftarrow\\{a-1\\}$ $\triangleright$ Set of breakpoints

6:for $t=a,\dots,b$ do

7: if $\max\\{\ell_{t},\bar{\ell}\\}<\min\\{u_{t},\bar{u}\\}$ then
$\triangleright$ Possible to extend constant interval

8: $\bar{\ell}\leftarrow\max\\{\ell_{t},\bar{\ell}\\}$,
$\bar{u}\leftarrow\min\\{u_{t},\bar{u}\\}$

9: else$\triangleright$ Breakpoint found

10: $\bar{f}\leftarrow\bar{f}+1$ $\triangleright$ Updates objective

11:
$\hat{\theta}_{\tau}\leftarrow(\bar{\ell}+\bar{u})/2\qquad\forall\tau=\texttt{last}(\Gamma),\dots,t-1$
$\triangleright$ Updates solution

12: $\Gamma\leftarrow\Gamma\cup\\{t\\}$

13: $\bar{\ell}\leftarrow\ell_{t}$, $\bar{u}\leftarrow u_{t}$ $\triangleright$
Reset bounds on interval

14: end if

15:end for

16:if $\texttt{last}(\Gamma)\neq b$ then $\triangleright$ Either $b=T+1$ or
$0\in[\bar{\ell},\bar{u}]$

17:
$\hat{\theta}_{\tau}\leftarrow(\bar{\ell}+\bar{u})/2\qquad\forall\tau=\texttt{last}(\Gamma),\dots,b$

18: $\Gamma\leftarrow\Gamma\cup\\{b\\}$

19:end if

20:Return $\bar{f}$;

######  Proposition 2.

Algorithm 1 runs in $\mathcal{O}(b-a)$ time and memory, and returns a solution
$\\{\hat{\theta}_{t}\\}_{a-1}^{b}$ that is optimal for (17). Moreover, given
any $b_{0}<b$, the truncated solution $\\{\hat{\theta}_{t}\\}_{a-1}^{b_{0}}$
is optimal for the problem associated with function $f(a,b_{0})$.

###### Proof.

The runtime readily follows by counting the number of iterations in Algorithm
1. We first prove the first sentence of the proposition. Let
$a-1=\tau_{0}<\tau_{1}<\tau_{2}<\ldots<\tau_{m}=b$ be the elements of the set
$\Gamma$ from Algorithm 1. By construction, $\max_{\tau_{i}\leq
t\leq\tau_{i+1}}\ell_{t}>\min_{\tau_{i}\leq t\leq\tau_{i+1}}u_{t}$ for all
$i=0,\dots,m-2$ (since elements are added to $\Gamma$ in line 12 only when the
condition in line 7 fails). It follows from Lemma 5 that

$m-1\leq\sum_{i=0}^{m-2}\sum_{t=\tau_{i}+1}^{\tau_{i+1}}\mathbb{I}[\theta_{t}\neq\theta_{t-1}]=\sum_{t=a-1}^{\tau_{m-1}}\mathbb{I}[\theta_{t}\neq\theta_{t-1}].$

Moreover, if $\ell_{b}=u_{b}=0$ and $0\not\in[\max_{\tau_{m-1}\leq t\leq
b-1}\ell_{t},\min_{\tau_{m-1}\leq t\leq b-1}u_{t}]$, then

$1\leq\sum_{t=\tau_{m-1}+1}^{b}\mathbb{I}[\theta_{t}\neq\theta_{t-1}]$

since $\theta_{t}$ needs to be nonzero for some $\tau_{m-1}\leq
t\leq\tau_{b-1}$ and $\theta_{b}=0$. Therefore, we find the lower bound on the
optimal objective value of (17) given by

$m-1+\mathbb{I}\left[\text{$\ell_{b}=u_{b}=0$ and
$0\not\in[\max_{\tau_{m-1}\leq t\leq b-1}\ell_{t},\min_{\tau_{m-1}\leq t\leq
b-1}u_{t}]$}\right]\leq\sum_{t=a-1}^{b}\mathbb{I}[\theta_{t}\neq\theta_{t-1}].$

We now check the objective value of Algorithm 1, given by $\bar{f}$. Note that
set $\bar{\Gamma}$ contains $m+1$ elements. The first element added (line 5)
does not increase the objective value. The next $m-1$ elements added (line 12)
increase the objective function by one. Finally, the last element added
increases the objective value by one unless $b=T+1$ (in which case
$-\ell_{b}=u_{b}=\infty$) or $0\in[\bar{\ell},\bar{u}]$ (line 18). Thus, the
value $\bar{f}$ returned by Algorithm 1 matches the lower, and is thus
optimal.

The second sentence of the proposition follows since Algorithm 1 is greedy:
the first $b_{0}$ steps used when computed $f(a,b)$ are exactly the same as
those used to compute $f(a,b_{0})$. ∎

From Proposition 2 we see that solving (16) for a fixed $a$ and all $b>a$ can
be done in linear time. It follows that solving (16) for all $a<b$ is possible
in $\mathcal{O}(T^{2})$ time.

####  5.3.2 Case $q\geq 1$

We now discuss how to solve (12) for $q\geq 1$. We assume for simplicity that
$0<a<b<T+1$, although the arguments can easily be extended to the cases with
$a=0$ or $b=T+1$.

######  Lemma 6.

If $q\geq 1$, there exists an optimal solution $\theta^{*}$ of (12) where:
$\bullet$ for all $a<t<b$, either
$\theta_{t}^{*}=(\theta_{t-1}^{*}+\theta_{t+1}^{*})/2$, or
$\theta_{t}^{*}\in\\{\ell_{t},u_{t}\\}$; $\bullet$ either
$\theta_{a}^{*}=\theta_{a+1}^{*}$, or $\theta_{a}^{*}\in\\{\ell_{a},u_{a}\\}$;
$\bullet$ either $\theta_{b-1}^{*}=\theta_{b-2}^{*}$, or
$\theta_{b-1}^{*}\in\\{\ell_{b-1},u_{b-1}\\}.$

###### Proof.

We only prove the first statement in the lemma, since the other two follow
from identical arguments. Given any index $a<\tau<b-1$, consider optimization
(12) where all values except for $\theta_{\tau}$ are fixed to any feasible
value (we omit terms that do not depend on $\theta_{\tau}$)

$\displaystyle\min_{\theta_{\tau}}$
$\displaystyle\left|\theta_{\tau}-\theta_{\tau-1}\right|^{q}+\left|\theta_{\tau+1}-\theta_{\tau}\right|^{q}$
(19) $\displaystyle\mathrm{s.t.}$
$\displaystyle\ell_{\tau}\leq\theta_{\tau}\leq u_{\tau}.$

Observe that if an optimal solution $\theta_{\tau}^{*}$ of (19) satisfies
$\theta_{\tau}^{*}>\max\\{\theta_{\tau-1},\theta_{\tau+1}\\}$, then
necessarily $\theta_{\tau}^{*}=\ell_{\tau}$, since otherwise it would be
possible to decrease $\theta_{\tau}^{*}$, improving the objective value.
Similarly the case
$\theta_{\tau}^{*}<\min\\{\theta_{\tau-1},\theta_{\tau+1}\\}$ implies that
$\theta_{\tau}^{*}=u_{\tau}$.

To prove the first statement, assume that
$\theta_{\tau}^{*}\not\in\\{\ell_{t},u_{t}\\}$, which can only happen if
$\min\\{\theta_{\tau-1},\theta_{\tau+1}\\}\leq\theta_{\tau}^{*}\leq\max\\{\theta_{\tau-1},\theta_{\tau+1}\\}$.
If $\theta_{\tau-1}\leq\theta_{\tau+1}$, then it follows that $\theta_{\tau}$
is optimal for

$\displaystyle\min_{\theta_{\tau}\in\mathbf{R}}\left(\theta_{\tau}-\theta_{\tau-1}\right)^{q}+\left(\theta_{\tau+1}-\theta_{\tau}\right)^{q}.$
(20)

First, suppose that $q>1$. Then, by convexity and differentiability, the
derivative of (20) at the optimal solution should be zero. This implies

$q(\theta_{t}-\theta_{t-1})^{q-1}-q(\theta_{t+1}-\theta_{t})^{q-1}=0\Leftrightarrow\theta_{t}-\theta_{t-1}=\theta_{t+1}-\theta_{t}\Leftrightarrow\theta_{\tau}=(\theta_{t+1}+\theta_{t-1})/2.$

Moreover, if $q=1$, it can be easily seen that any feasible value
$\theta_{\tau}\in[\theta_{\tau-1},\theta_{\tau+1}]$ is optimal. Thus, at least
one of the following five cases holds: (i)
$\theta_{\tau}^{*}=(\theta_{t+1}+\theta_{t-1})/2\in[\theta_{\tau-1},\theta_{\tau+1}]$
is optimal; (ii)
$\theta_{\tau}^{*}=\ell_{t}\in[\theta_{\tau-1},\theta_{\tau+1}]$ is optimal;
(iii) $\theta_{\tau}^{*}=u_{t}\in[\theta_{\tau-1},\theta_{\tau+1}]$; (iv)
$\theta_{\tau}^{*}=\ell_{t}\geq\theta_{\tau+1}$ is optimal; or (v)
$\theta_{\tau}^{*}=u_{t}\leq\theta_{\tau-1}$ is optimal. Since all fives cases
satisfy
$\theta_{\tau}^{*}\in\\{\ell_{\tau},(\theta_{t+1}+\theta_{t-1})/2,u_{\tau}\\}$,
the statement of the lemma holds.

The case $\theta_{\tau-1}\geq\theta_{\tau+1}$ is handled identically,
concluding the proof. ∎

From Lemma 6, we find that $\theta_{t}$ (when not fixed to a bound) is
obtained by interpolation of the adjacent values. In other words, in an
optimal solution to (12), there is a set of indexes
$a\leq\tau_{1}<\tau_{2}<\dots<\tau_{m}\leq b-1$ corresponding to variables
fixed to a bound, and for a variable $t$ who is not at a bound with
$\tau_{j}<t<\tau_{j+1}$, we find its optimal value given by

$\theta_{t}^{*}=\theta_{\tau_{j}}^{*}+\frac{t-\tau_{j}}{\tau_{j+1}-\tau_{j}}(\theta_{\tau_{j+1}}^{*}-\theta_{\tau_{j}}^{*}).$
(21)

Moreover, $\theta_{1}>a$ (thus $\theta_{a}^{*}\not\in\\{\ell_{a},u_{a}\\}$),
then it also follows from Lemma 6 that
$\theta_{a}^{*}=\theta_{a_{+}1}^{*}=\theta_{\tau_{1}}$; similarly, if
$\tau_{m}<b-1$, then
$\theta_{\tau_{m}}^{*}=\theta_{\tau_{m}+1}^{*}=\theta_{b-1}$.

We now propose a dynamic programming approach for solving (12) with $q\geq 1$.
Let set $S=\left\\{\texttt{lb},\texttt{mid},\texttt{ub}\right\\}$, where lb
stands for “lower bound”, ub stands for “upper bound”, and mid stands for
“neither lower or upper bound”. Given $s_{a},s_{b}\in S$, define
$f^{0}(a,s_{a},b,s_{b})$ as the objective value of (12) when $\theta_{a}$ and
$\theta_{b-1}$ are set to their bounds indicated by $s_{a}$ and $s_{b}$,
respectively, and the remaining variables are set according to (21). Note that
$f^{0}(a,s_{a},b,s_{b})$ can be computed in $\mathcal{O}(T)$ time, and
$f^{0}(a,s_{a},b,s_{b})=\infty$ if the sequence indicated by (21) is
infeasible.

Given $s_{a},s_{b}\in S$, define $f^{1}(a,s_{a},b,s_{b})$ as the optimal
solution of (12), and observe that $f^{1}$ can be computed via the recursion
(where $a\leq t<b-1$ is the largest index of a variable set to a bound)

$f^{1}(a,s_{a},b,s_{b})=\min_{t\in\\{a,\dots,b-1\\},s\in\\{\texttt{lb},\texttt{ub}\\}}f^{1}(a,s_{1},t+1,s)+f^{0}(t,s,b,s_{b}).$
(22)

Thus, assuming that values $f^{1}(a,s_{a},t,s)$ with $s\in S$ and $t<b$ have
been computed, (22) can be computed in $\mathcal{O}(T)$ time, and computing
$f^{1}$ for all $\mathcal{O}(T^{2})$ combinations of parameters requires
$\mathcal{O}(T^{3})$ time. Figure 3 illustrates an example of the recursion
(22).

##  6 Statistical Guarantees

(a) $f^{0}(0,\texttt{ub},5,\texttt{lb})=\infty$ (b)
${f^{1}(0,\texttt{mid},5,\texttt{lb})=f^{1}(0,\texttt{mid},3,\texttt{ub})+f^{0}(2,\texttt{ub},5,\texttt{lb})}$

Figure 3: Example of recursion (22). (Left) Evaluation of the cost setting
$\theta_{0}=u_{0}$, $\theta_{4}=\ell_{4}$ and remaining variables by
interpolation (21). Since this solution is infeasible due to
$\theta_{2}\not\in[\ell_{2},u_{2}]$, the cost is infinity. (Right) Optimal
solution, which sets $\theta_{2}=u_{2}$.

In the previous section, we showed that ProxGL can be solved efficiently for
all values of the regularization parameter $0\leq\gamma\leq 1$. However, the
efficient solvability of ProxGL is contingent upon the availability of the
approximate backward mapping $\widetilde{F}^{*}$. On the other hand, according
to Theorem 2, the accuracy of this backward mapping and the choice of
$\lambda_{t}$ directly controls the estimation error of the obtained
parameters. In this section, our goal is to provide statistically and
computationally efficient ways to obtain approximate backward mappings for two
classes of time-varying MRFs, namely time-varying GMRFs and time-varying
DMRFs. We then use our results to provide end-to-end sample complexity
guarantees for the solution of ProxGL.

To achieve this goal, we first present a general statistical bound on the
backward mapping deviation. Recall that $\mathcal{R}$ is the domain of the
backward mapping and $\mathcal{B}_{\infty}(\mu^{\star},\alpha)$ is the
$\ell_{\infty}$-ball with radius $\rho$ centered at $\mu^{\star}$.

######  Definition 3 (Local lipschitzness of
$\boldsymbol{\widetilde{F}^{*}}$).

$\widetilde{F}^{*}$ is $(L,\alpha)$-locally Lipschitz for some $L,\alpha>0$ if
$\|\widetilde{F}^{*}(\mu)-\widetilde{F}^{*}(\mu^{\prime})\|_{\infty}\leq
L\|\mu-\mu^{\prime}\|_{\infty}$ for every
$\mu,\mu^{\prime}\in\mathcal{R}\cap\mathcal{B}_{\infty}(\mu^{\star},\alpha)$.

The local Lipschitzness of the approximate backward mapping will play a
crucial role in our subsequent statistical analysis. Roughly speaking, it
entails that if the mean parameters
$\mu,\mu^{\prime}\in\mathcal{R}\cap\mathcal{B}^{\infty}_{\mu^{\star}}(\alpha)$
are close, so are their images under $\widetilde{F}^{*}$.

######  Assumption 1 (Exponential tail bound).

Let
$\widehat{\mu}_{t}=\frac{1}{N_{t}}\sum_{i=1}^{N_{t}}\phi\left(x^{(i)}_{t}\right)$
be the empirical mean parameter and
$\mu^{\star}_{t}=\mathbb{E}_{\theta_{t}^{\star}}[\phi(x_{t})]$ its true
counterpart. There exist $c,c^{\prime},v>0$ such that

$\displaystyle\mathbb{P}(|\widehat{\mu}_{t;i}-{\mu}^{\star}_{t;i}|\geq\delta)\leq
c^{\prime}\exp(-cN_{t}\delta^{2})\quad\text{for every}\quad\delta\in(0,v],\ \
i=1,\dots,p,\ \ t=0,\dots,T$

The above assumption implies that the empirical mean parameters are
concentrated around their expectation with a tail bound that decreases
exponentially fast. As will be shown later, this assumption is satisfied for a
wide range of MRFs, including GMRFs and DMRFs.

Recall that
$\left\lVert\theta_{t}^{\star}-\widetilde{F}^{*}(\widehat{\mu}_{t})\right\rVert_{\infty}$
plays an important role in Theorem 2. Under the local Lipschitzness of
$\widetilde{F}^{*}$ and Assumption 1, we will next show that the backward
mapping deviation
$\left\lVert\theta_{t}^{\star}-\widetilde{F}^{*}(\widehat{\mu}_{t})\right\rVert_{\infty}$
can be controlled with high probability.

######  Proposition 3 (Backward mapping deviation).

Suppose that Assumption 1 is satisfied and $\widetilde{F}^{*}$ is
$(L,\alpha)$-locally Lipschitz with some $\alpha,L>0$. Then, for any constant
$\tau\geq 0$ and sample size
$N_{t}\geq\frac{1+\tau}{c\min\\{v^{2},\alpha^{2}\\}}\cdot\log p$, we have

$\displaystyle\left\lVert\theta_{t}^{\star}-\widetilde{F}^{*}(\widehat{\mu}_{t})\right\rVert_{\infty}\leq\underbrace{\left\lVert\theta_{t}^{\star}-\widetilde{F}^{*}({\mu}^{\star}_{t})\right\rVert_{\infty}}_{:=\Delta_{t}}+\sqrt{\frac{(1+\tau)L^{2}}{c}\cdot\frac{\log
p}{N_{t}}}\\!\\!\\!\\!,\quad\text{with probability of }1-p^{-\tau}.$

###### Proof.

Applying Assumption 1 with the choice of
$\delta=\sqrt{\frac{1+\tau}{c}\cdot\frac{\log p}{N_{t}}}$ implies that

$\displaystyle\mathbb{P}\left(|\widehat{\mu}_{t;i}-{\mu}^{\star}_{t;i}|\leq\sqrt{\frac{1+\tau}{c}\cdot\frac{\log
p}{N_{t}}}\right)\geq 1-c^{\prime}p^{-(1+\tau)},$

where $\delta\leq v$ is automatically satisfied due to the lower bound on
$N_{t}$. Then, a simple union bound over index $i$ leads to

$\left\lVert\widehat{\mu}_{t}-\mu^{\star}_{t}\right\rVert_{\infty}\leq\sqrt{\frac{1+\tau}{c}\cdot\frac{\log
p}{N_{t}}}\leq\alpha,\quad\text{with probability of $1-c^{\prime}p^{-\tau}$.}$

Under this event, one can write

$\displaystyle\left\lVert\theta_{t}^{\star}-\widetilde{F}^{*}(\widehat{\mu}_{t})\right\rVert_{\infty}$
$\displaystyle\leq\left\lVert\theta_{t}^{\star}-\widetilde{F}^{*}({\mu}^{\star}_{t})\right\rVert_{\infty}+\left\lVert\widetilde{F}^{*}(\widehat{\mu}_{t})-\widetilde{F}^{*}({\mu}^{\star}_{t})\right\rVert_{\infty}$
$\displaystyle\leq\left\lVert\theta_{t}^{\star}-\widetilde{F}^{*}({\mu}^{\star}_{t})\right\rVert_{\infty}+L\left\lVert\widehat{\mu}_{t}-{\mu}^{\star}_{t}\right\rVert_{\infty}$
$\displaystyle\leq\left\lVert\theta_{t}^{\star}-\widetilde{F}^{*}({\mu}^{\star}_{t})\right\rVert_{\infty}+\sqrt{\frac{(1+\tau)L^{2}}{c}\cdot\frac{\log
p}{N_{t}}}$

where in the second inequality we used the $(L,\alpha)$-local Lipschitzness of
$\widehat{F}^{*}$. ∎

Proposition 3 shows that in order to guarantee a small backward mapping
deviation, $\widetilde{F}^{*}$ should satisfy two properties: first, it must
incur a small error at the true mean parameter (captured by $\Delta_{t}$), and
second, it must be $(\alpha,L)$-Lipschitz for some $\alpha$ and $L$. However,
there is often a trade-off between these two properties. To see this trade-
off, consider $\widetilde{F}^{*}(\Sigma)=\Sigma^{-1}$ as a natural choice of
the approximate backward mapping for GMRFs. Evidently, we have $\Delta_{t}=0$
for $\widetilde{F}^{*}(\Sigma)=\Sigma^{-1}$. However, this choice of
approximate backward mapping loses its local Lipschitzness as $\Sigma$
approaches singularity (i.e., when the sample size falls below the dimension).
To alleviate this issue, we propose a proxy backward mapping that strikes a
balance between $\Delta_{t}$ and the local Lipschitzness.

###  6.1 Approximate Backward Mapping for GMRFs

Suppose that at any given time $t$, a sequence of data samples
$\left\\{x_{t}^{(i)}\right\\}_{i=1}^{N_{t}}\subset\mathbb{R}^{n}$ are
collected from a time-varying GMRF with distribution (2). Our next lemma is
borrowed from [Ravikumar et al., 2011] and shows that the empirical mean
parameters (i.e., sample covariance matrices) satisfy Assumption 1.

######  Lemma 7 (Lemma 1, Ravikumar et al. [2011]).

Suppose the samples $\left\\{x_{t}^{(i)}\right\\}_{i=1}^{N_{t}}$ are collected
from a time-varying GMRF. Then, the empirical mean parameter defined as
$\widehat{\Sigma}_{t}=({1}/{N_{t}})\sum_{i=1}^{N_{t}}x_{t}^{(i)}{x_{t}^{(i)}}^{\top}$
satisfies Assumption 1 with $c^{\prime}=4$,
$c=3200\max_{i}\\{\Sigma^{\star}_{t;ii}\\}^{2}$, and
$v=40\max_{i}\\{\Sigma^{\star}_{t;ii}\\}$.

Recall that the canonical parameter corresponds to the inverse covariance
matrix for GMRFs. Hence, we have $p=n^{2}$. To address the singularity of the
backward mapping when $N_{t}<n$, we use the approximate backward mapping
introduced by Yang et al. [2014]: consider the soft-thresholding operator
$\mathrm{ST}_{\nu}:\mathbb{R}^{n\times n}\to\mathbb{R}^{n\times n}$ defined as
$\texttt{ST}_{\nu}(M)_{ij}=M_{ij}-\mathrm{sign}(M_{ij})\min\\{|M_{ij}|,\nu\\}$
if $i\not=j$, and $\texttt{ST}_{\nu}(M)_{ij}=M_{ij}$ if $i=j$. The proposed
approximate backward mapping is then given by
$\widetilde{F}^{*}(\widehat{\Sigma}_{t})=[\texttt{ST}_{\nu}(\widehat{\Sigma}_{t})]^{-1}$.
Clearly, we have $\Delta_{t}=0$ for the choice of $\nu=0$. However, we will
later show that $\nu>0$ is essential for guaranteeing the $(L,\alpha)$-local
Lipschitzness of $\widetilde{F}^{*}$ when $N_{t}<n$.

######  Assumption 2 (Bounded norm).

There exist constant numbers $\kappa_{1}<\infty$, $\kappa_{2}>0$, and
$\kappa_{3}<\infty$ such that

$\displaystyle\|\Theta^{\star}_{t}\|_{\infty}\leq\kappa_{1},\;\inf_{w:\|w\|_{\infty}=1}\|\Sigma^{\star}_{t}w\|_{\infty}\geq\kappa_{2},\;\|\Sigma^{\star}_{t}\|_{\infty/\infty}\leq\kappa_{3},\quad\text{for
every $t=0,\dots,T$.}$

The above assumption is fairly mild and implies that the true covariance
matrices and their inverses have bounded norms. Without loss of generality and
to streamline the presentation, we assume that $\kappa_{1},\kappa_{3}\geq 1$
and $\kappa_{2}\leq 1$.

######  Definition 4 (Weak sparsity).

We say $\Sigma\in\mathbb{R}^{n\times n}$ is $(s,r)$-weakly sparse for some
$s\geq 0$ and $0\leq r<1$ if $\max_{i}\sum_{j=1}^{n}|\Sigma_{ij}|^{r}\leq
s_{r}$.

Intuitively, we say that “the true covariance matrices are weakly sparse” if
$\\{\Sigma^{\star}_{t}\\}_{t=0}^{T}$ are $(s,r)$-weakly sparse with $s_{r}\ll
n$ for some $0\leq r<1$. The notion of weak sparsity extends the classical
notion of sparsity for matrices. Indeed, except for a few special cases, a
sparse matrix does not have a sparse inverse. Consequently, a sparse precision
matrix may not lead to a sparse covariance matrix. However, a large class of
sparse precision matrices have weakly sparse inverses:

  * •

If $\Theta^{\star}_{t}$ has a banded structure with small bandwidth, then it
is known that the elements of $\Sigma^{\star}_{t}=\Theta_{t}^{-1}$ enjoy
exponential decay away from the main diagonal [Demko et al., 1984, Kershaw,
1970]. Under such condition, one can verify that $s\leq\frac{C}{1-\rho^{r}}$
for some constant $C>0$ and $\rho<1$.

  * •

A similar statement holds for a class of inverse covariance matrices whose
support graphs have large average path length [Benzi and Razouk, 2007, Benzi
and Simoncini, 2015]; a large class of inverse covariance matrices with row-
and column-sparse structures satisfy this condition.

  * •

Suppose that $\Theta^{\star}_{t}$ is diagonally dominant with
$\Theta^{\star}_{t,ii}-\sum_{j\not=i}|\Theta^{\star}_{t,ij}|>0$ for
$i=1,2,\dots,d$. Then, a simple derivation reveals that $\Sigma^{\star}_{t}$
is $(s,r)$-weakly sparse with $s\leq
d^{1-r}/\left(\min_{i}\left\\{\Theta^{\star}_{t,ii}-\sum_{j\not=i}|\Theta^{\star}_{t,ij}|\right\\}\right)^{r}$
for any $0\leq r<1$. The diagonally dominant structures naturally arise in the
context of graphical model inference; see [Egilmez et al., 2016] for various
applications, including those with graph Laplacian structures.

Under Assumption 2 and $(s,r)$-weak sparsity of the true covariance matrices,
we show that the proposed approximate backward mapping
$\widetilde{F}^{*}(\widehat{\Sigma}_{t})=[\texttt{ST}_{\nu}(\widehat{\Sigma}_{t})]^{-1}$
is $(L,\alpha)$-locally Lipschitz and has small $\Delta_{t}$.

######  Lemma 8.

Suppose that Assumption 2 is satisfied and $\Sigma^{\star}_{t}$ is
$(s,r)$-weakly sparse for some $s\geq 0$ and $0\leq r<1$. Then, the following
statements hold for the approximate backward mapping
$\widetilde{F}^{*}(\widehat{\Sigma}_{t})=[\texttt{ST}_{\nu}(\widehat{\Sigma}_{t})]^{-1}$:

  * •

We have $\Delta_{t}\leq({2\kappa_{1}}/{\kappa_{2}})\nu$;

  * •

Assume that $\nu\leq(\kappa_{2}/(20s))^{1/(1-r)}$. Then,
$\widetilde{F}^{*}(\Sigma)$ is $(L,\alpha)$-locally Lipschitz with
$L=4/\kappa_{2}^{2}$ and $\alpha=\min\\{\kappa_{2}\nu^{r}/(24s),\nu/2\\}$.

To provide the proof of the above lemma, we will rely on the following
intermediate lemma adapted from [Yang et al., 2014, Lemma 1].

######  Lemma 9 (Yang et al. [2014]).

We have

$\displaystyle\left\lVert\texttt{ST}_{\nu}(\Sigma)-\Sigma^{\star}_{t}\right\rVert_{\infty}\leq
5s_{r}\nu^{1-r}+3s_{r}\nu^{-r}\alpha,$

for every $\Sigma\in\mathcal{B}_{\Sigma_{t}^{\star}}(\alpha)$ and
$\alpha\leq\nu/2$.

##### Proof of Lemma 8

We first start with the proof of the first statement. One can write

$\displaystyle\Delta_{t}=\left\lVert\Theta_{t}^{\star}-[\texttt{ST}_{\nu}(\Sigma_{t}^{\star})]^{-1}\right\rVert_{\infty/\infty}$
$\displaystyle\leq\left\lVert\Theta^{\star}_{t}\right\rVert_{\infty}\left\lVert[\texttt{ST}_{\nu}(\Sigma_{t}^{\star})]^{-1}\right\rVert_{\infty}\left\lVert\texttt{ST}_{\nu}(\Sigma_{t}^{\star})-\Sigma_{t}^{\star}\right\rVert_{\infty/\infty}$
$\displaystyle\leq\kappa_{1}\nu\left\lVert[\texttt{ST}_{\nu}(\Sigma_{t}^{\star})]^{-1}\right\rVert_{\infty}$
(23)

where the last inequality follows from
$\left\lVert\texttt{ST}_{\nu}(\Sigma_{t}^{\star})-\Sigma_{t}^{\star}\right\rVert_{\infty/\infty}\leq\nu$
and our assumption
$\left\lVert\Theta^{\star}_{t}\right\rVert_{\infty}\leq\kappa_{1}$. Next, we
provide an upper bound for
$\left\lVert[\texttt{ST}_{\nu}(\Sigma_{t}^{\star})]^{-1}\right\rVert_{\infty}$.
One can write

$\displaystyle\left\lVert\texttt{ST}_{\nu}(\Sigma_{t}^{\star})w\right\rVert_{\infty}$
$\displaystyle\geq\left\lVert\Sigma_{t}^{\star}w\right\rVert_{\infty}-\left\lVert\left(\texttt{ST}_{\nu}(\Sigma_{t}^{\star})-\Sigma_{t}^{\star}\right)w\right\rVert_{\infty}$
$\displaystyle\geq\left(\kappa_{2}-\left\lVert\left(\texttt{ST}_{\nu}(\Sigma_{t}^{\star})-\Sigma_{t}^{\star}\right)\right\rVert_{\infty}\right)\left\lVert
w\right\rVert_{\infty}$ $\displaystyle\geq(\kappa_{2}/2)\left\lVert
w\right\rVert_{\infty}$

where the last inequality is due to our choice of $\nu$. This implies that
$\left\lVert[\texttt{ST}_{\nu}(\Sigma_{t}^{\star})]^{-1}\right\rVert_{\infty}\leq
2/\kappa_{2}$. Combining this inequality with (23) leads to
$\Delta_{t}\leq(2\kappa_{1}/\kappa_{2})\nu$.

Moving on to the proof of the second statement, we invoke Lemma 9 to write

$\displaystyle\left\lVert\texttt{ST}_{\nu}(\Sigma)-\Sigma^{\star}_{t}\right\rVert_{\infty}\leq
5s_{r}\nu^{1-r}+3s_{r}\nu^{-r}\alpha\qquad\text{for every
$\Sigma\in\mathcal{B}_{\Sigma_{t}^{\star}}(\alpha)$ and $\alpha\leq\nu/2$.}$

Based on our assumption on $\nu$, we have $5s_{r}\nu^{1-r}\leq\kappa_{2}/4$.
Similarly, based on our assumed upper bound on $\alpha$, we have
$3s_{r}\nu^{-r}\alpha\leq\kappa_{2}/4$. This implies that
$\left\lVert\texttt{ST}_{\nu}(\Sigma)-\Sigma^{\star}_{t}\right\rVert_{\infty}\leq\kappa_{2}/2$
for every $\Sigma\in\mathcal{B}_{\Sigma_{t}^{\star}}(\alpha)$. On the other
hand, for every
$\Sigma,\Sigma^{\prime}\in\mathcal{B}_{\Sigma_{t}^{\star}}(\alpha)$, we have

$\displaystyle\left\lVert[\texttt{ST}_{\nu}(\Sigma)]^{-1}-[\texttt{ST}_{\nu}(\Sigma^{\prime})]^{-1}\right\rVert_{\infty/\infty}\leq\left\lVert[\texttt{ST}_{\nu}(\Sigma)]^{-1}\right\rVert_{\infty}\left\lVert[\texttt{ST}_{\nu}(\Sigma^{\prime})]^{-1}\right\rVert_{\infty}\left\lVert\Sigma-\Sigma^{\prime}\right\rVert_{\infty/\infty}$
(24)

Similarly, one can write

$\displaystyle\|\texttt{ST}_{\nu}({\Sigma})w\|_{\infty}$
$\displaystyle\geq\|\Sigma^{\star}_{t}w\|_{\infty}-\left\|(\texttt{ST}_{\nu}({\Sigma})-\Sigma^{\star}_{t})w\right\|_{\infty}$
$\displaystyle\geq\left(\kappa_{2}-\left\|\texttt{ST}_{\nu}({\Sigma})-\Sigma^{\star}_{t}\right\|_{\infty}\right)\|w\|_{\infty}$
$\displaystyle\geq(\kappa_{2}/2)\|w\|_{\infty}$

where in the last inequality, we used
$\left\lVert\texttt{ST}_{\nu}(\Sigma)-\Sigma^{\star}_{t}\right\rVert_{\infty}\leq\kappa_{2}/2$.
This in turn implies
$\left\lVert[\texttt{ST}_{\nu}(\Sigma)]^{-1}\right\rVert_{\infty}\leq
2/\kappa_{2}$. Similarly, we have
$\left\lVert[\texttt{ST}_{\nu}(\Sigma^{\prime})]^{-1}\right\rVert_{\infty}\leq
2/\kappa_{2}$. Combining these inequalities with (24) leads to

$\displaystyle\left\lVert[\texttt{ST}_{\nu}(\Sigma)]^{-1}-[\texttt{ST}_{\nu}(\Sigma^{\prime})]^{-1}\right\rVert_{\infty/\infty}\leq(4/\kappa_{2}^{2})\left\lVert\Sigma-\Sigma^{\prime}\right\rVert_{\infty/\infty},\qquad\text{for
every $\Sigma,\Sigma^{\prime}\in\mathcal{B}_{\Sigma_{t}^{\star}}(\alpha)$.}$

This completes the proof. $\hfill\square$

Lemma 8 highlights the trade-off in choosing the threshold $\nu$: to keep the
error $\Delta_{t}$ small, it is desirable to pick a small value for $\nu$.
However, $\nu$ cannot be too small, as it would shrink the radius of the
neighborhood within which $\widetilde{F}^{*}$ is locally Lipschitz. Our next
lemma shows that it suffices to pick $\nu$ large enough so that
$\left\lVert\widehat{\Sigma}_{t}-\Sigma_{t}^{\star}\right\rVert_{\infty/\infty}\leq\alpha$
with high probability, where $\alpha$ is a function of $\nu$ defined in Lemma
8. To streamline the presentation, we first define the following quantities
which will be used extensively in our subsequent arguments:

$\displaystyle\mathcal{N}_{\mathrm{g}}:=\left(\frac{4s_{r}}{\kappa_{3}\kappa_{2}}\right)^{\frac{2}{1-r}},\qquad\mathcal{E}_{\mathrm{g}}:=\frac{10\kappa_{1}+2\kappa_{3}^{-1}}{5\kappa_{2}^{2}},$
$\displaystyle\bar{\Theta}_{t}:=\\!\\!\min_{(i,j)\in\mathcal{S}_{t}}|\Theta^{\star}_{t;ij}|,\
$
$\displaystyle\Delta\bar{\Theta}_{t}:=\\!\\!\min_{(i,j)\in\mathcal{D}_{t}}|\Theta^{\star}_{t;ij}-\Theta^{\star}_{t-1;ij}|,\
\Theta_{t}^{\min}:=\begin{cases}\min\left\\{\bar{\Theta}_{t},\frac{\Delta\bar{\Theta}_{t}}{2},\frac{\Delta\bar{\Theta}_{t-1}}{2}\right\\}&t\geq
1\\\
\min\left\\{\bar{\Theta}_{t},\frac{\Delta\bar{\Theta}_{t}}{2}\right\\}&t=0\end{cases}$

The quantities $\mathcal{N}_{\mathrm{g}}$ and $\mathcal{E}_{\mathrm{g}}$ will
be used in our subsequent error bounds, while $\bar{\Theta}_{t}$,
$\Delta\bar{\Theta}_{t}$, and $\Theta_{t}^{\min}$ are defined to capture the
minimum nonzero elements of the true parameter and their differences.

######  Lemma 10.

Suppose that Assumption 2 is satisfied and $\Sigma^{\star}_{t}$ is
$(s,r)$-weakly sparse for some $s\geq 0$ and $0\leq r<1$. Suppose that
$N_{t}\geq\mathcal{N}_{\mathrm{g}}\log n$. Then, with the choice of
$\nu_{t}=\sqrt{\log n/N_{t}}$, we have

$\left\lVert\Theta_{t}^{\star}-[\texttt{ST}_{\nu_{t}}(\widehat{\Sigma}_{t})]^{-1}\right\rVert_{\infty/\infty}\leq\mathcal{E}_{\mathrm{g}}\sqrt{\frac{\log
n}{N_{t}}},\qquad\text{with probability of $1-4n^{-15}$.}$

###### Proof.

In order to apply Proposition 3, we need to verify that
$N_{t}\geq(1+\tau)/(c\min\\{v^{2},\alpha^{2}\\})\log n$. This is readily
implied by setting $\tau=15$, and invoking the derived values for $c$, $v$,
and $\alpha$ in Lemmas 7 and 8 together with our choice of $\nu_{t}$ and the
assumed lower bound on $N_{t}$. The result then follows from Proposition 3 and
the provided upper bounds on $\Delta_{t}$ and $L$ in Lemma 8. ∎

Given Lemma 10, we are ready to provide the sample complexity of ProxGL for
GMRFs.

######  Theorem 3 (Sample Complexity of GMRFs).

Suppose that a sequence of data samples
$\left\\{x_{t}^{(i)}\right\\}_{i=1}^{N_{t}}$ are collected from the
distribution (2) at every $t=0,\dots,T$. Suppose that Assumption 2 is
satisfied and $\Sigma^{\star}_{t}$ is $(s,r)$-weakly sparse for some $s\geq 0$
and $0\leq r<1$. Consider ProxGL with
$\widetilde{F}^{*}(\widehat{\Sigma}_{t})=[\texttt{ST}_{\nu_{t}}(\widehat{\Sigma}_{t})]^{-1}$
and parameters $\nu_{t}\asymp\sqrt{\log n/N_{t}}$ and
$\lambda_{t}\asymp\mathcal{E}_{\mathrm{g}}\sqrt{\log n/N_{t}}$. The following
statements hold with probability of $1-4Tn^{-15}$:

  * -

Estimation error: Suppose that $N_{t}\gtrsim\mathcal{N}_{\mathrm{g}}\log p$
for every $t=0,1,\dots,T$. For any $q\in\\{0\\}\cup[1,\infty)$, the solution
of ProxGL with a temporal $\ell_{q}$-regularizer satisfies

$\displaystyle\left\lVert\widehat{\Theta}_{t}-\Theta^{\star}_{t}\right\rVert_{\infty/\infty}\lesssim\mathcal{E}_{\mathrm{g}}\sqrt{\frac{\log
n}{N_{t}}},\qquad\text{ for every }t=0,\dots,T.$ (25)

  * -

Sparsistency for smoothly-changing GMRFs. Suppose that the GMRF is smoothly
changing with parameters $(q,D),q\geq 1,D\geq 0$. Suppose that
$N_{t}\gtrsim\left(\mathcal{N}_{\mathrm{g}}\vee(\mathcal{E}_{\mathrm{g}}/\bar{\Theta}_{t})^{2}\right)\log
n$ for every $t=0,\dots,T$. Then, the solution of ProxGL with
$0<\gamma<1/(1+D)$ and a temporal $\ell_{q}$-regularizer satisfies

$\displaystyle\mathrm{supp}(\widehat{\Theta}_{t})=\mathrm{supp}({\Theta}^{\star}_{t}),\qquad\forall
t=0,\dots,T.$ (26)

  * -

Sparsistency for sparsely-changing GMRFs. Suppose that the GMRF is sparsely
changing with parameter $D_{0}\geq 0$. Suppose that
$N_{t}\gtrsim\left(\mathcal{N}_{\mathrm{g}}\vee(\mathcal{E}_{\texttt{g}}/{\Theta}^{\min}_{t})^{2}\right)\log
n$ for every $t=0,\dots,T$. Then, with any choice of $0<\gamma<1$, the optimal
solution of ProxGL with temporal $\ell_{0}$-regularizer satisfies:

$\displaystyle\mathrm{supp}(\widehat{\Theta}_{t})=\mathrm{supp}({\Theta}^{\star}_{t}),$
$\displaystyle\qquad\forall t=0,\dots,T,$ (27)
$\displaystyle\mathrm{supp}(\widehat{\Theta}_{t}-\widehat{\Theta}_{t-1})=\mathrm{supp}({\Theta}^{\star}_{t}-{\Theta}^{\star}_{t-1}),$
$\displaystyle\qquad\forall t=1,\dots,T.$

###### Proof.

The proof is a direct consequence of Lemma 10 and Theorem 2. The details are
omitted for brevity. ∎

######  Remark 1.

We note that Theorem 3 does not impose any condition on the weak sparsity
parameter $s$ of the true covariance matrices; instead, it uses this parameter
to control the sample complexity of our method. In other words, our
theoretical result still holds for large values of $s$, provided that the
number of samples scales accordingly.

###  6.2 Approximate Backward Mapping for DMRFs

Suppose that at any given time $t$, a sequence of discrete data samples
$\left\\{x_{t}^{(i)}\right\\}_{i=1}^{N_{t}}$ are collected from a DMRF with
distribution (3). Our next lemma shows that the empirical mean
parameters—corresponding to the empirical node- and edge-wise marginal
probabilities—satisfy Assumption 1.

######  Lemma 11 (Hoeffding’s inequality).

For every $1\leq i,j\leq p$ and $k,l\in\mathcal{K}$, define
$\widehat{\mu}_{t;ik}=(1/N_{t})\sum_{s=1}^{N_{t}}\mathbb{I}[x^{(s)}_{t;i}=k]$
and
$\widehat{\mu}_{t;ijkl}=(1/N_{t})\sum_{s=1}^{N_{t}}\mathbb{I}[x^{(s)}_{t;i}=k,x^{(s)}_{t;j}=l]$
as the empirical node- and edge-wise marginal probabilities. Then, for all
$\widehat{\mu}_{t;ik}$ and $\widehat{\mu}_{t;ijkl}$, Assumption 1 holds with
$c^{\prime}=2$, $c=2$, and $v=\infty$.

###### Proof.

Note that both $\widehat{\mu}_{t;ik}$ and $\widehat{\mu}_{t;ijkl}$ are
empirical averages of binary random variables. Therefore, the proof follows
from a direct application of Hoeffding’s inequality; see e.g. [Wainwright,
2019, Proposition 2.5]. ∎

Bresler et al. [2014] showed that, unlike GMRFs, obtaining the backward
mapping of DMRFs is NP-hard, even if random variables are restricted to binary
values. Therefore, one can only hope to provide an approximation of this
backward mapping. A promising candidate for such approximation is the so-
called tree-reweighted entropy mapping $\widetilde{F}^{*}$ introduced by
Wainwright et al. [2003]:

$\displaystyle\left[\widetilde{F}^{*}_{\mathrm{trw}}(\widehat{\mu}_{t})\right]_{\phi}=\begin{cases}\log\widehat{\mu}_{t;ik},&\text{if}\quad\phi=ik\\\
\rho_{t;ij}\log\left(\frac{\widehat{\mu}_{t;ijkl}}{\widehat{\mu}_{t;ik}\widehat{\mu}_{t;jl}}\right),&\text{if}\quad\phi=ijkl\end{cases}$
(28)

The weights $0\leq\rho_{t;ij}\leq 1$ are selected from the so-called spanning
tree polytope. For any spanning tree $\mathcal{T}_{t}$ of the Markov graph
$\mathcal{G}_{t}$, one can set $\rho_{t;ij}=1$ if $(i,j)\in\mathcal{T}_{t}$
and $\rho_{t;ij}=0$ otherwise. However, due to the unknown nature of
$\mathcal{G}_{t}$, its spanning trees are also unknown. To circumvent this
issue, it is common to select $\rho_{t;lk}=1$ for all tuples $(t,l,k)$
[Wainwright, 2002].

Due to the intractability of the backward mapping, in general one cannot
expect the defined approximate backward mapping to be tight even if
$\widehat{\mu}_{t}=\mu^{\star}_{t}$. This implies that, unlike GMRFs, it may
not be possible to push $\Delta_{t}$ towards zero. Nonetheless, we will
establish the local Lipschitzness of the above approximate backward mapping
under the following mild condition.

######  Assumption 3.

There exists $\mu_{\min}>0$ such that:

  * •

We have $\mu^{\star}_{t;ik}=P(X_{i}=k)\geq\mu_{\min}$ for every $i,k$.

  * •

We have $\mu^{\star}_{t;ijkl}=P(X_{i}=k,X_{j}=l)\geq\mu_{\min}$ for every
$i,j,k,l$.

The above assumption requires the density of the joint probability
distribution to be strictly positive on every atom of the distribution. Given
the above assumption, we next prove the local Lipschitzness of our defined
approximate backward mapping.

######  Lemma 12.

Suppose that Assumption 3 is satisfied. Then, the approximate backward mapping
(28) is $(L,\alpha)$-locally Lipschitz with $L=6/\mu_{\min}$ and
$\alpha=\mu_{\min}/2$.

###### Proof.

Suppose that $\phi=ijkl$ is an arbitrary index. For any
$\mu,\mu^{\prime}\in\mathcal{B}_{\mu^{\star}_{t}}(\alpha)$, we have

$\displaystyle\left|\left[\widetilde{F}^{*}_{\mathrm{trw}}({\mu})\right]_{\phi}-\left[\widetilde{F}^{*}_{\mathrm{trw}}({\mu^{\prime}})\right]_{\phi}\right|$
$\displaystyle=\left|\rho_{t;ij}\log\left(\frac{{\mu}_{ijkl}}{{\mu}_{ik}{\mu}_{jl}}\right)-\rho_{t;ij}\log\left(\frac{{\mu}^{\prime}_{ijkl}}{{\mu}^{\prime}_{ik}{\mu}^{\prime}_{jl}}\right)\right|$
$\displaystyle\leq\left|\log\left(\frac{\mu^{\prime}_{ijkl}}{\mu_{ijkl}}\right)\right|+\left|\log\left(\frac{\mu^{\prime}_{ik}}{\mu_{ik}}\right)\right|+\left|\log\left(\frac{\mu^{\prime}_{jl}}{\mu_{jl}}\right)\right|$

On the other hand, one can write

$\displaystyle\left|\log\left(\frac{\mu^{\prime}_{ijkl}}{\mu_{ijkl}}\right)\right|$
$\displaystyle\leq\left|\log\left(\frac{\min\\{\mu^{\prime}_{ijkl},\mu_{ijkl}\\}+|\mu^{\prime}_{ijkl}-\mu_{ijkl}|}{\min\\{\mu^{\prime}_{ijkl},\mu_{ijkl}\\}}\right)\right|$
$\displaystyle\leq\left|\log\left(1+\frac{|\mu^{\prime}_{ijkl}-\mu_{ijkl}|}{\min\\{\mu^{\prime}_{ijkl},\mu_{ijkl}\\}}\right)\right|$
$\displaystyle\leq\frac{|\mu^{\prime}_{ijkl}-\mu_{ijkl}|}{\min\\{\mu^{\prime}_{ijkl},\mu_{ijkl}\\}}$
$\displaystyle\leq\left({2}/{\mu_{\min}}\right)|\mu^{\prime}_{ijkl}-\mu_{ijkl}|$

where the third inequality follows from $\log(1+x)\leq x$ for any $x>-1$ and
the last inequality is due to
$\mu,\mu^{\prime}\in\mathcal{B}_{\mu^{\star}_{t}}(\alpha)$ and our choice of
$\alpha$. Similarly, one can show that
$\left|\log\left({\mu^{\prime}_{ik}}/{\mu_{ik}}\right)\right|\leq\left({2}/{\mu_{\min}}\right)|\mu^{\prime}_{ik}-\mu_{ik}|$
and
$\left|\log\left({\mu^{\prime}_{jl}}/{\mu_{jl}}\right)\right|\leq\left({2}/{\mu_{\min}}\right)|\mu^{\prime}_{jl}-\mu_{jl}|$.
This implies that

$\left|\left[\widetilde{F}^{*}_{\mathrm{trw}}({\mu})\right]_{ijkl}-\left[\widetilde{F}^{*}_{\mathrm{trw}}({\mu^{\prime}})\right]_{ijkl}\right|\leq(6/\mu_{\min})\cdot\max\\{|\mu^{\prime}_{ijkl}-\mu_{ijkl}|,|\mu^{\prime}_{ik}-\mu_{ik}|,|\mu^{\prime}_{jl}-\mu_{jl}|\\}$

Similarly, for an arbitrary index $\phi=ik$, we have
$\left|\left[\widetilde{F}^{*}_{\mathrm{trw}}({\mu})\right]_{\phi}-\left[\widetilde{F}^{*}_{\mathrm{trw}}({\mu^{\prime}})\right]_{\phi}\right|\leq(2/\mu_{\min})|\mu^{\prime}_{ik}-\mu_{ik}|$.
This implies that

$\left\lVert\widetilde{F}^{*}_{\mathrm{trw}}({\mu})-\widetilde{F}^{*}_{\mathrm{trw}}({\mu^{\prime}})\right\rVert_{\infty}\leq(6/\mu_{\min})\left\lVert\mu-\mu^{\prime}\right\rVert_{\infty},$

thereby completing the proof. ∎

Lemma 12 combined with Proposition 3 leads to the following concentration
bound on the backward mapping deviation.

######  Lemma 13.

Suppose that Assumption 3 is satisfied. Moreover, suppose that
$N_{t}\geq(16/\mu_{\min}^{2})\log p$. Then, we have

$\displaystyle\left\lVert\theta_{t}^{\star}-\widetilde{F}^{*}_{\mathrm{trw}}(\widehat{\mu}_{t})\right\rVert_{\infty}\leq\Delta_{t}+\frac{12}{\mu_{\min}}\sqrt{\frac{\log
p}{N_{t}}},\qquad\text{ with probability of $1-2p^{-8}$.}$ (29)

###### Proof.

The proof immediately follows from Lemma 12 and Proposition 3. ∎

Lemma 13 implies that the backward mapping deviation can be upper bounded with
$2\Delta_{t}$, so long as
$N_{t}\geq\left(12/(\mu_{\min}\Delta_{t})\right)^{2}\log p$. The above result
implies that, with logarithmic sample size, the backward mapping deviation
stays in the order of the approximation error $\Delta_{t}$. Equipped with the
above lemma, we are ready to present the sample complexity of ProxGL for
DMRFs. Similar to GMRFs, we define the following quantities:

$\displaystyle\bar{\theta}_{t}:=\min_{(i,j)\in\mathcal{S}_{t}}|\theta^{\star}_{t;ij}|,\quad$
$\displaystyle\Delta\bar{\theta}_{t}:=\min_{(i,j)\in\mathcal{D}_{t}}|\theta^{\star}_{t;ij}-\theta^{\star}_{t-1;ij}|,\quad\theta_{t}^{\min}:=\begin{cases}\min\left\\{\bar{\theta}_{t},\frac{\Delta\bar{\theta}_{t}}{2},\frac{\Delta\bar{\theta}_{t-1}}{2}\right\\}&t\geq
1\\\
\min\left\\{\bar{\theta}_{t},\frac{\Delta\bar{\theta}_{t}}{2}\right\\}&t=0\end{cases}$

######  Theorem 4 (Sample Complexity of DMRFs).

Suppose that a sequence of data samples
$\left\\{x_{t}^{(i)}\right\\}_{i=1}^{N_{t}}$ are collected from the
distribution (3) at every $t=0,\dots,T$. Suppose that Assumption 3 is
satisfied. Consider ProxGL with the backward mapping
$\widetilde{F}_{\mathrm{trw}}^{*}$ defined as (28) and parameter
$\lambda_{t}\asymp\Delta_{t}+\frac{1}{\mu_{\min}}\sqrt{\frac{\log p}{N_{t}}}$.
The following statements hold with probability of $1-2Tp^{-8}$:

  * -

Estimation error: Suppose that $N_{t}\gtrsim\mu_{\min}^{-2}\log p$ for every
$t=0,\dots,T$. For any $q\in\\{0\\}\cup[1,\infty)$, the solution of ProxGL
with a temporal $\ell_{q}$-regularizer satisfies

$\displaystyle\left\lVert\widehat{\theta}_{t}-\theta^{\star}_{t}\right\rVert_{\infty}\lesssim\Delta_{t}+\frac{1}{\mu_{\min}}\sqrt{\frac{\log
p}{N_{t}}},\qquad\text{ for every }t=0,\dots,T.$

  * -

Sparsistency for smoothly-changing DMRFs. Suppose that the time-varying DMRF
is smoothly changing with parameters $(q,D),q\geq 1,D\geq 0$. Suppose that
$\Delta_{t}\leq\bar{\theta}_{t}/2$ and
$N_{t}\gtrsim(\mu_{\min}\cdot\bar{\theta}_{t})^{-2}\log p$ for every
$t=0,\dots,T$. Then, the solution of ProxGL with $0<\gamma<1/(1+D)$ and a
temporal $\ell_{q}$-regularizer satisfies

$\displaystyle\mathrm{supp}\left(\widehat{\theta}_{t}\right)=\mathrm{supp}\left({\theta}^{\star}_{t}\right),\qquad\forall
t=0,\dots,T.$ (30)

  * -

Sparsistency for sparsely-changing DMRFs. Suppose that the time-varying DMRF
is sparsely changing with parameter $D_{0}\geq 0$. Suppose that
$\Delta_{t}\leq\theta^{\min}_{t}/2$ and
$N_{t}\gtrsim(\mu_{\min}\cdot\theta^{\min}_{t})^{-2}\log n$ for every
$t=0,1,\dots,T$. Then, with any choice of $0<\gamma<1$, the optimal solution
of ProxGL with temporal $\ell_{0}$-regularizer satisfies:

$\displaystyle\mathrm{supp}\left(\widehat{\theta}_{t}\right)=\mathrm{supp}({\theta}^{\star}_{t}),$
$\displaystyle\qquad\forall t=0,\dots,T,$ (31)
$\displaystyle\mathrm{supp}\left(\widehat{\theta}_{t}-\widehat{\theta}_{t-1}\right)=\mathrm{supp}({\theta}^{\star}_{t}-{\theta}^{\star}_{t-1}),$
$\displaystyle\qquad\forall t=1,\dots,T.$

###### Proof.

The proof is a direct consequence of Lemma 13 and Theorem 2. The details are
omitted for brevity. ∎

According to Theorem 4, the exact sparsity recovery is only possible if the
approximation error of the backward mapping $\Delta_{t}$ is small.

###  6.3 Kernel Averaging

In most realistic time-varying MRFs, the underlying graphical model changes
continuously over time as the samples continue to arrive. For instance, in
financial markets, the underlying stock correlation network may change quickly
in response to global events. Therefore, a stock holder needs to identify the
sharp changes in the market and rebalance their portfolio “on the go” [Talih
and Hengartner, 2005, Hallac et al., 2017]. Evidently, in these applications,
the sample size $N_{t}$ may be significantly smaller than what is required for
ProxGL to provide a reliable estimation of the canonical parameters.

To address the scarcity of data in this setting, we leverage and combine the
information provided by the samples over time. In particular, we propose to
replace the empirical mean parameters with their weighted averages over time,
where the weights are obtained from a nonparametric kernel
$K:\mathbb{R}\to\mathbb{R}_{+}$. Without loss of generality, suppose that
$N_{t}=1$222If $N_{t}>1$, the sufficient statistics $\phi(x_{s})$ in (32) can
be replaced by $\frac{1}{N_{t}}\sum_{i=1}^{N_{t}}\phi(x^{(i)}_{s})$. for every
$t=0,\dots,T$ and consider the weighted mean parameter

$\displaystyle\widehat{\mu}^{\mathrm{ker}}_{t}=\sum_{s=0}^{T}w(s,t)\phi(x_{s}),\quad\text{where
}\ w(s,t)=\frac{1}{Th}K\left(\frac{s-t}{Th}\right).$ (32)

Here, $K(\cdot)$ is a symmetric nonnegative kernel that satisfies a set of
mild conditions which will be delineated later. Moreover, the parameter $h$ is
the bandwidth of the kernel, controlling the decay rate of the weights. The
key insight behind kernel averaging is simple: at any given time $t$, we
estimate the mean parameters by taking the weighted average of the samples
over time, where the weights are obtained from a kernel that assigns smaller
weights to samples that are temporally farther away from $t$. The idea of
kernel averaging (also known as kernel smoothing) dates back to late 60s
[Epanechnikov, 1969, Altman, 1992, Wand and Jones, 1994], and has been
recently used in graphical model inference [Greenewald et al., 2017, Zhou et
al., 2010].

To ensure that such weighted averaging is statistically consistent, we must
assume that its true counterpart changes slowly over time to ensure that the
samples collected over time can be used to reveal partial information about
the true mean parameter at time $t$. More precisely, our hope is to guarantee
that
$\sum_{s=0}^{T}w(s,t)\phi(x_{s})\approx\mathbb{E}[\phi(x_{t})]=\mu^{\star}_{t}$.
To formalize this intuition, we assume that $\mu^{\star}_{t}=\mu(t/T)$, for
some twice continuously differentiable $\mu(x):[0,1]\to\mathbb{R}^{n}$ that
satisfies the following assumption:333A function can be differentiable only on
open sets. Without loss of generality and with a slight abuse of notation, we
replace the derivatives of $\mu(x)$ at $x=0$ and $x=1$ with their right and
left derivatives, respectively.

######  Assumption 4 (Bounded derivatives).

There exist constants $0\leq\Gamma_{0},\Gamma_{1},\Gamma_{2}<\infty$ such that
$\left\lVert\mu(x)\right\rVert_{\infty}\leq\Gamma_{0}$,
$\left\lVert\mu^{\prime}(x)\right\rVert_{\infty}\leq\Gamma_{1}$, and
$\left\lVert\mu^{\prime\prime}(x)\right\rVert_{\infty}\leq\Gamma_{2}$, for
every $0\leq x\leq 1$.

The above assumption implies that the function $\mu$ and its first derivative
change smoothly.

######  Assumption 5 (Kernel properties).

The kernel $K(x)$ satisfies the following conditions:

  * -

$\int_{-1}^{1}K(x)dx=1$;

  * -

There exist constants $K_{l},K_{u},K_{1},K_{2}<\infty$ such that:

$\displaystyle 0<K_{l}\leq K(x/h)\leq K_{u},\quad|K^{\prime}(x/h)|\leq
K_{1}\cdot h^{-1},\quad|K^{\prime\prime}(x/h)|\leq K_{2}\cdot h^{-2},$

for every $-h\leq x\leq h$ and $0<h\leq 1$.

The above assumption holds for most standard kernels. For instance, it is easy
to verify that the uniform kernel $K(x)=1/2$ with domain $[-1,1]$ satisfies
Assumption 5 with parameters $K_{l}=K_{u}=1/2$ and $K_{1}=K_{2}=0$. As another
example, consider the truncated Gaussian kernel $K(x)=Ce^{-x^{2}/2}$, where
$C=(\Phi(1)-\Phi(-1))^{-1}\approx 1.465$ is the normalizer of the kernel and
$\Phi$ is the CDF of the normal distribution. Simple calculation reveals that
this kernel satisfies Assumption 5 with parameters $K_{l}=Ce^{-1/2}$,
$K_{u}=K_{1}=C$ and $K_{2}=2C$.

Given the above two assumptions, our next goal is to derive a generalization
of Proposition 3 to the kernel averaged empirical mean parameters. Recall that
Proposition 3 relies on an upper bound on
$\left\lVert\widehat{\mu}_{t}-{\mu}^{\star}_{t}\right\rVert_{\infty}$. With
kernel averaging, the empirical mean parameter $\widehat{\mu}_{t}$ is replaced
with $\widehat{\mu}_{t}^{\mathrm{ker}}$. As a result, we need to obtain a
concentration bound on the deviation of $\widehat{\mu}_{t}^{\mathrm{ker}}$
from $\mu^{\star}_{t}$. To this goal, we write

$\displaystyle\left\lVert\widehat{\mu}_{t}^{\mathrm{ker}}-\mu^{\star}_{t}\right\rVert_{\infty}\leq\underbrace{\left\|\mu^{\star}_{t}-\mathbb{E}[\widehat{\mu}_{t}^{\mathrm{ker}}]\right\|_{\infty}}_{\text{kernel
error}}+\underbrace{\left\|\widehat{\mu}_{t}^{\mathrm{ker}}-\mathbb{E}[\widehat{\mu}_{t}^{\mathrm{ker}}]\right\|_{\infty}}_{\text{deviation
from kernel mean}}$

Our next two lemmas provide upper bounds on the kernel error and the deviation
of kernel mean.

######  Lemma 14 (Kernel error).

Under Assumptions 4 and 5, we have

$\displaystyle\left\|\mu^{\star}_{t}-\mathbb{E}[\widehat{\mu}_{t}^{\mathrm{ker}}]\right\|_{\infty}\leq\frac{\Gamma_{2}}{4}\cdot
h^{2}+{\Gamma_{1}}\cdot h+\frac{\Gamma_{2}K_{u}}{24T^{2}}\cdot
h^{-1}+\frac{\Gamma_{1}K_{1}}{12T^{2}}\cdot
h^{-3}+\frac{\Gamma_{0}K_{2}}{24T^{2}}\cdot h^{-5}$ (33)

###### Proof.

One can write

$\displaystyle\left\|\mu^{\star}_{t}-\mathbb{E}[\widehat{\mu}_{t}^{\mathrm{ker}}]\right\|_{\infty}\leq$
$\displaystyle\underbrace{\left\|\mu^{\star}_{t}-\int_{0}^{1}\frac{1}{h}K\left(\frac{x-t/T}{h}\right)\mu(x)dx\right\|_{\infty}}_{\mathcal{A}}$
$\displaystyle+\underbrace{\left\|\int_{0}^{1}\frac{1}{h}K\left(\frac{x-t/T}{h}\right)\mu(x)dx-\mathbb{E}[\widehat{\mu}_{t}^{\mathrm{ker}}]\right\|_{\infty}}_{\mathcal{B}}$

According to Zhou et al. [2010, Lemma 14] we have
$\mathcal{A}\leq\frac{\Gamma_{2}}{4}\cdot h^{2}+{\Gamma_{1}}\cdot h$.
Therefore, it remains to control $\mathcal{B}$. To this goal, define
$f(x):=\frac{1}{h}K\left(\frac{x-t/T}{h}\right)\mu(x)$. Note that

$\displaystyle\mathbb{E}[\widehat{\mu}_{t}^{\mathrm{ker}}]=\frac{1}{T}\sum_{s=0}^{T}\frac{1}{h}K\left(\frac{s-t}{Th}\right)\mu\left(\frac{s}{T}\right)=\frac{1}{T}\sum_{s=0}^{T}f\left(\frac{s}{T}\right).$

This implies that $\mathbb{E}[\widehat{\mu}_{t}^{\mathrm{ker}}]$ is the
Riemann approximation of $\int_{0}^{1}f(x)dx$. Therefore, we have

$\displaystyle\mathcal{B}=\left\|\int_{0}^{1}f(x)dx-\frac{1}{T}\sum_{s=0}^{T}f\left(\frac{s}{T}\right)\right\|_{\infty}\leq\frac{\max_{0\leq
x\leq 1}\|f^{\prime\prime}(x)\|_{\infty}}{24T^{2}}$ (34)

where the second inequality follows from a standard bound on the Riemann sum
[Hughes-Hallett et al., 2020]. Next, we provide an upper bound on $\max_{0\leq
x\leq 1}\|f^{\prime\prime}(x)\|_{\infty}$. It is easy to verify that

$\displaystyle f^{\prime\prime}(x)=h^{-3}\cdot
K^{\prime\prime}\left(\frac{x-t/T}{h}\right)\mu(x)+2h^{-2}\cdot
K^{\prime}\left(\frac{x-t/T}{h}\right)\mu^{\prime}(x)+h^{-1}K\left(\frac{x-t/T}{h}\right)\mu^{\prime\prime}(x),$
(35)

Therefore, with Assumption 5, we have

$\displaystyle\|f^{\prime\prime}(x)\|_{\infty}\leq h^{-5}\cdot
K_{2}\Gamma_{0}+2h^{-3}\cdot K_{1}\Gamma_{1}+h^{-1}\cdot K_{u}\Gamma_{2}.$

The above inequality combined (34) completes the proof. ∎

######  Lemma 15 (Deviation from kernel mean).

Suppose that Assumptions 4 and 5 are satisfied. Moreover, suppose that
$T\geq(\tau+1)(K_{u}/K_{l})^{4}\log p$ for an arbitrary $\tau>0$. Then, we
have

$\displaystyle\left\lVert\widehat{\mu}_{t}^{\mathrm{ker}}-\mathbb{E}[\widehat{\mu}_{t}^{\mathrm{ker}}]\right\rVert_{\infty}\leq\sqrt{\frac{(\tau+1)K_{u}^{2}\log
p}{cTh^{2}}},\quad\text{with probability}\ 1-c^{\prime}p^{-\tau}$

###### Proof.

Let us define $g_{s}=(1/h)K((s-t)/(Th))$ and $Y_{s}=g_{s}\phi(x_{s})$. We have
$\widehat{\mu}_{t}^{\mathrm{ker}}=(1/T)\sum_{s=0}^{T}Y_{s}$. According to
Assumption 1, $Y_{s}$ satisfies

$\displaystyle\mathbb{P}(|Y_{s;i}-\mathbb{E}[Y_{s;i}]|\geq\delta)\leq
c^{\prime}\exp(-(c/g_{s}^{2})\delta^{2})\quad\text{for
every}\quad\delta\in(0,g_{s}v],\ \ i=1,\dots,n.$

Therefore, $Y_{s;i}$ is a sub-exponential random variable with parameters
$(g_{s}^{2}/(2c),1/(g_{s}v))$; see [Wainwright, 2019, Chapter 2] for the
definition of sub-exponential random variables). On the other hand, it is easy
to see that
$\mathbb{E}[\widehat{\mu}^{\mathrm{ker}}_{t;i}]=\frac{1}{T}\sum_{s=0}^{T}g_{s}\mu_{s}=\frac{1}{T}\sum_{s=0}^{T}\mathbb{E}[Y_{s;i}]$.
Therefore, one can invoke the concentration bound on the sum of sub-
exponential random variables [Wainwright, 2019, Proposition 2.9] to obtain

$\displaystyle\mathbb{P}\left(\left|\frac{1}{T}\sum_{s=0}^{T}Y_{s;i}-\mathbb{E}[\widehat{\mu}^{\mathrm{ker}}_{t;i}]\right|\geq\delta\right)\leq
c^{\prime}\exp\left(-\frac{cT^{2}}{\sum_{s=0}^{T}g_{s}^{2}}\cdot\delta^{2}\right)\quad\text{for
every}\quad\delta\in\left(0,\frac{\sum_{s=0}^{T}g_{s}^{2}}{T\max_{s}g_{s}}\right].$

A simple union bound leads to

$\mathbb{P}\left(\left\|\frac{1}{T}\sum_{s=0}^{T}Y_{s}-\mathbb{E}[\widehat{\mu}^{\mathrm{ker}}_{t}]\right\|_{\infty}\geq\delta\right)\leq
c^{\prime}\exp\left(\log
p-\frac{cT^{2}}{\sum_{s=0}^{T}g_{s}^{2}}\cdot\delta^{2}\right)\\\ \text{for
every}\quad\delta\in\left(0,\frac{\sum_{s=0}^{T}g_{s}^{2}}{T\max_{s}g_{s}}\right].$

Upon choosing
$\delta=\sqrt{\frac{(\tau+1)\left(\sum_{s=0}^{T}g_{s}^{2}\right)\log
p}{T^{2}c}}$, we have

$\displaystyle\left\|\frac{1}{T}\sum_{s=0}^{T}Y_{s}-\mathbb{E}[\widehat{\mu}^{\mathrm{ker}}_{t}]\right\|_{\infty}\leq\sqrt{\frac{(\tau+1)\left(\sum_{s=0}^{T}g_{s}^{2}\right)\log
p}{T^{2}c}},$

with probability $1-c^{\prime}p^{-\tau}$, provided that

$\displaystyle\sqrt{\frac{(\tau+1)\left(\sum_{s=0}^{T}g_{s}^{2}\right)\log
p}{T^{2}c}}\leq\frac{\sum_{s=0}^{T}g_{s}^{2}}{T\max_{s}g_{s}}.$ (36)

To verify the last inequality, recall that according to Assumption 5, we have
$K_{l}\leq K(\cdot)\leq K_{u}$, which in turn implies that $h^{-1}K_{l}\leq
g_{s}\leq h^{-1}K_{u}$ for every $s$. This results in

$\displaystyle\sqrt{\frac{(\tau+1)\left(\sum_{s=0}^{T}g_{s}^{2}\right)\log
p}{T^{2}c}}\leq\sqrt{\frac{(\tau+1)K_{u}^{2}\log
p}{Th^{2}c}},\qquad\frac{K_{l}^{2}}{hK_{u}}\leq\frac{\sum_{s=0}^{T}g_{s}^{2}}{T\max_{s}g_{s}}.$

Therefore, to satisfy (36), it suffices to have

$\displaystyle\sqrt{\frac{(\tau+1)K_{u}^{2}\log
p}{Th^{2}c}}\leq\frac{K_{l}^{2}}{hK_{u}}$

which is satisfied with $T\geq(\tau+1)(K_{u}/K_{l})^{4}\log p$. This completes
the proof. ∎

Combining Lemmas 14 and 15 leads to an upper bound on
$\left\lVert\widehat{\mu}_{t}^{\mathrm{ker}}-\mu_{t}^{\star}\right\rVert_{\infty}$.
To simplify our subsequent arguments, we define
$\mathcal{C}^{\mathrm{ker}}=\Gamma_{2}/4+\Gamma_{1}+\Gamma_{2}K_{u}/24+\Gamma_{1}K_{1}/12+\Gamma_{0}K_{2}/24$.

######  Lemma 16.

Suppose that Assumptions 4 and 5 are satisfied, and
$T\geq(\tau+1)(K_{u}/K_{l})^{4}\log p$ for any arbitrary $\tau\geq 0$. Then,
upon choosing $h=T^{-1/4}$, we have

$\displaystyle\left\lVert\widehat{\mu}_{t}^{\mathrm{ker}}-\mu_{t}^{\star}\right\rVert_{\infty}\leq\frac{\mathcal{C}^{\mathrm{ker}}}{T^{1/4}}+\frac{K_{u}\sqrt{(\tau+1)\log
p}}{\sqrt{c}T^{1/4}},\quad\text{with probability of $1-c^{\prime}p^{-\tau}$.}$
(37)

###### Proof.

The proof follows from Lemmas 14 and 15 after replacing $h=T^{-1/4}$. ∎

Given Lemma 16, we are now ready to provide an upper bound for the backward
mapping deviation of the kernel averaged empirical mean parameter.

######  Proposition 4 (Backward mapping deviation with kernel averaging).

Suppose that Assumptions 4 and 5 are satisfied and $\widetilde{F}^{*}$ is
$(L,\alpha)$-locally Lipschitz for some $\alpha,L>0$. Then, for any $\tau\geq
0$ and
$T\geq\max\left\\{16(\mathcal{C}^{\mathrm{ker}}/\alpha)^{4},\left(\frac{16K_{u}^{4}(\tau+1)^{2}}{c^{2}K_{l}^{4}\alpha^{4}}\cdot\log^{2}p\right)\right\\}$,
we have

$\displaystyle\left\lVert\theta_{t}^{\star}-\widetilde{F}^{*}(\widehat{\mu}^{\mathrm{ker}}_{t})\right\rVert_{\infty}\leq\underbrace{\left\lVert\theta_{t}^{\star}-\widetilde{F}^{*}({\mu}^{\star}_{t})\right\rVert_{\infty}}_{:=\Delta_{t}}+\frac{\mathcal{C}L}{T^{1/4}}+\frac{K_{u}L\sqrt{(\tau+1)\log
p}}{\sqrt{c}T^{1/4}},$

with probability of $1-c^{\prime}p^{-\tau}$.

###### Proof.

Due to Lemma 16 and
$T\geq\max\left\\{16(\mathcal{C}^{\mathrm{ker}}/\alpha)^{4},\left(\frac{16K_{u}^{4}(\tau+1)^{2}}{c^{2}K_{l}^{4}\alpha^{4}}\cdot\log^{2}p\right)\right\\}$,
we have

$\displaystyle\left\lVert\widehat{\mu}_{t}^{\mathrm{ker}}-\mu_{t}^{\star}\right\rVert_{\infty}\leq\frac{\mathcal{C}}{T^{1/4}}+\frac{K_{u}\sqrt{(\tau+1)\log
p}}{\sqrt{c}T^{1/4}}\leq\alpha$

with probability of $1-c^{\prime}p^{-\tau}$. Conditioned on this event, we
have

$\displaystyle\left\lVert\theta_{t}^{\star}-\widetilde{F}^{*}(\widehat{\mu}^{\mathrm{ker}}_{t})\right\rVert_{\infty}$
$\displaystyle\leq\left\lVert\theta_{t}^{\star}-\widetilde{F}^{*}({\mu}^{\star}_{t})\right\rVert_{\infty}+\left\lVert\widetilde{F}^{*}(\widehat{\mu}^{\mathrm{ker}}_{t})-\widetilde{F}^{*}({\mu}^{\star}_{t})\right\rVert_{\infty}$
$\displaystyle\leq\left\lVert\theta_{t}^{\star}-\widetilde{F}^{*}({\mu}^{\star}_{t})\right\rVert_{\infty}+L\left\lVert\widehat{\mu}^{\mathrm{ker}}_{t}-{\mu}^{\star}_{t}\right\rVert_{\infty}$
$\displaystyle\leq\left\lVert\theta_{t}^{\star}-\widetilde{F}^{*}({\mu}^{\star}_{t})\right\rVert_{\infty}+\frac{\mathcal{C}L}{T^{1/4}}+\frac{K_{u}L\sqrt{(\tau+1)\log
p}}{\sqrt{c}T^{1/4}}$

where in the second inequality we used the $(L,\alpha)$-local Lipschitzness of
$\widetilde{F}^{*}$. ∎

Relying on Proposition 4, we next present the sample complexity of ProxGL for
GMRF and DMRF with kernel averaging. To this goal, we first provide the
explicit form of the kernel empirical mean parameter (32) in both settings:

  * •

For GMRF, the kernel mean parameters are the weighted sample covariance
matrices defined as
$\widehat{\Sigma}^{\mathrm{ker}}_{t}=\sum_{s=0}^{T}w(s,t)x_{s}x_{s}^{\top}$,
for every $t=0,\dots,T$.

  * •

For DMRF, the kernel mean parameters are the weighted empirical node- and
edge-wise marginal probabilities respectively defined as
$\widehat{\mu}^{\mathrm{ker}}_{t;ik}=\sum_{s=0}^{T}w(s,t)\mathbb{I}[x_{s;i}=k]$
and
$\widehat{\mu}^{\mathrm{ker}}_{t;ijkl}=\sum_{s=0}^{T}w(s,t)\mathbb{I}[x_{s;i}=k,x_{s;j}=l]$,
for every $i,j,k,l$.

To streamline the presentation, we define the following quantities:

$\displaystyle\mathcal{E}^{\mathrm{ker}}_{\mathrm{g}}:=\frac{\mathcal{C}^{\mathrm{ker}}+\kappa_{1}+K_{u}\kappa_{3}^{-1}}{\kappa_{2}^{2}},\quad\mathcal{N}^{\mathrm{ker}}_{\mathrm{g}}:=\max\left\\{\left(\frac{\mathcal{C}^{\mathrm{ker}}s}{\kappa_{2}}\right)^{\frac{4}{1-r}},\left(\frac{K_{u}s}{K_{l}\sqrt{\kappa_{3}}\kappa_{2}}\right)^{\frac{4}{1-r}}\right\\}$

where $\kappa_{1},\kappa_{2},\kappa_{3},r,s$ are defined in Assumptions 2 and
Definition 4.

######  Proposition 5 (Sample complexity of GMRF with kernel averaging).

Suppose that a sequence of data samples $\left\\{x_{t}\right\\}_{t=0}^{T}$ are
collected from the distribution (2). Moreover, suppose that Assumptions 2, 4,
and 5 are satisfied and $\Sigma^{\star}_{t}$ is $(s,r)$-weakly sparse for some
$s\geq 0$ and $0\leq r<1$. Consider ProxGL with
$\widetilde{F}^{*}(\widehat{\Sigma}^{\mathrm{ker}}_{t})=[\texttt{ST}_{\nu_{t}}(\widehat{\Sigma}^{\mathrm{ker}}_{t})]^{-1}$
and parameters $\nu_{t}\asymp{\sqrt{\log n}}/{T^{1/4}}$,
$\lambda_{t}\asymp{\mathcal{E}_{g}\log n}/{T^{1/4}}$, and $h\asymp T^{-1/4}$.
Then, with probability of $1-4Tn^{-15}$, the following statements hold:

  * -

Estimation error: Suppose that
$T\gtrsim\mathcal{N}^{\mathrm{ker}}_{\mathrm{g}}\log^{2}n$. For any
$q\in\\{0\\}\cup[1,\infty)$, the solution of ProxGL with a temporal
$\ell_{q}$-regularizer satisfies

$\displaystyle\left\lVert\widehat{\Theta}_{t}-\Theta^{\star}_{t}\right\rVert_{\infty/\infty}\lesssim\mathcal{E}^{\mathrm{ker}}_{\mathrm{g}}{\frac{\sqrt{\log
n}}{T^{1/4}}},\qquad\text{ for every }\ t=0,1,\dots,T.$

  * -

Sparsistency for smoothly-changing GMRF. Suppose that the GMRF is smoothly
changing with parameters $(q,D),q\geq 1,D\geq 0$. Suppose that
$T\gtrsim\left(\mathcal{N}^{\mathrm{ker}}_{\mathrm{g}}\vee(\mathcal{E}^{\mathrm{ker}}_{\mathrm{g}}/\bar{\Theta}_{t})^{4}\right)\log^{2}n$.
Then, the solution of ProxGL with $0<\gamma<1/(1+D)$ and a temporal
$\ell_{q}$-regularizer satisfies

$\displaystyle\mathrm{supp}(\widehat{\Theta}_{t})=\mathrm{supp}({\Theta}^{\star}_{t}),\qquad\forall
t=0,\dots,T.$ (38)

  * -

Sparsistency for sparsely-changing GMRF. Suppose that the GMRF is sparsely
changing with with sparsity parameter $D_{0}\geq 0$. Suppose that
$T\gtrsim\left(\mathcal{N}^{\mathrm{ker}}_{\mathrm{g}}\vee(\mathcal{E}^{\mathrm{ker}}_{\mathrm{g}}/{\Theta}^{\min}_{t})^{4}\right)\log^{2}n$.
Then, with any choice of $0<\gamma<1$, the optimal solution of ProxGL with
temporal $\ell_{0}$-regularizer satisfies

$\displaystyle\mathrm{supp}(\widehat{\Theta}_{t})=\mathrm{supp}({\Theta}^{\star}_{t}),$
$\displaystyle\qquad\forall t=0,\dots,T,$ (39)
$\displaystyle\mathrm{supp}(\widehat{\Theta}_{t}-\widehat{\Theta}_{t-1})=\mathrm{supp}({\Theta}^{\star}_{t}-{\Theta}^{\star}_{t-1}),$
$\displaystyle\qquad\forall t=1,\dots,T.$

###### Proof.

The proof follows that of Theorem 3. First, due to Lemma 8, we have
$\Delta_{t}\leq(2\kappa_{1}/\kappa_{2})\nu_{t}$ and
$\widetilde{F}^{*}(\Sigma)$ is $(L,\alpha)$-locally Lipschitz with
$L=4/\kappa_{2}^{2}$ and
$\alpha=\min\\{\kappa_{2}\nu_{t}^{r}/(24s),\nu_{t}/2\\}$, provided that
$\nu_{t}\leq(\kappa_{2}/(40s))^{1/(1-r)}$. This combined with Proposition 4
and the selected values for $\nu_{t}$ and $T$ leads to

$\left\lVert\Theta_{t}^{\star}-[\texttt{ST}_{\nu_{t}}(\widehat{\Sigma}^{\mathrm{ker}}_{t})]^{-1}\right\rVert_{\infty/\infty}\leq\mathcal{E}^{\mathrm{ker}}_{\mathrm{g}}\frac{\sqrt{\log
n}}{T^{1/4}},\quad\text{for every $t$ with probability of $1-4Tn^{-15}$.}$

Therefore, upon choosing
$\lambda_{t}\asymp{\mathcal{E}^{\mathrm{ker}}_{\mathrm{g}}\log n}/{T^{1/4}}$,
Theorem 2 can be invoked to write
$\left\lVert\widehat{\Theta}_{t}-\Theta^{\star}_{t}\right\rVert_{\infty/\infty}\leq\lambda_{t}$,
which completes the proof of the first statement. To prove the correct
sparsity recovery, according to Theorem 2 we additionally need to satisfy
$\lambda_{t}\leq\bar{\Theta}_{t}/2$ for smoothly-changing GMRFs, and
$\lambda_{t}\leq{\Theta}^{\min}_{t}/2$ for sparsely-changing GMRFs. These are
guaranteed to hold with the above error bound and the assumed lower bounds on
$T$. ∎

Finally, we provide the sample complexity if inferring DMRFs with kernel
averaging. Similar to GMRF, we define the following quantities:

$\displaystyle\mathcal{E}^{\mathrm{ker}}_{\mathrm{d}}:=\frac{\mathcal{C}^{\mathrm{ker}}+K_{u}}{\mu_{\min}},\quad\mathcal{N}^{\mathrm{ker}}_{\mathrm{d}}:=\max\left\\{\left(\frac{\mathcal{C}^{\mathrm{ker}}}{\mu_{\min}}\right)^{4},\left(\frac{K_{u}}{K_{l}\mu_{\min}}\right)^{4}\right\\}$

######  Proposition 6 (Sample complexity of DMRF with kernel averaging).

Suppose that a sequence of data samples $\left\\{x_{t}\right\\}_{t=0}^{T}$ are
collected from the distribution (3). Moreover, suppose that Assumptions 3, 4,
and 5 are satisfied. Consider ProxGL with parameter
$\lambda_{t}\asymp\Delta_{t}+\mathcal{E}^{\mathrm{ker}}_{\mathrm{d}}\sqrt{\log
p}/T^{1/4}$ and the backward mapping $\widetilde{F}_{\mathrm{trw}}^{*}$
defined as (28). Then, with probability of $1-2Tp^{-8}$, the following
statements hold:

  * -

Estimation error: Suppose that
$T\gtrsim\mathcal{N}^{\mathrm{ker}}_{\mathrm{d}}\log^{2}p$. For any
$q\in\\{0\\}\cup[1,\infty)$, the solution of ProxGL with a temporal
$\ell_{q}$-regularizer satisfies

$\displaystyle\left\lVert\widehat{\theta}_{t}-\theta^{\star}_{t}\right\rVert_{\infty}\lesssim\mathcal{E}^{\mathrm{ker}}_{\mathrm{d}}{\frac{\sqrt{\log
p}}{T^{1/4}}},\qquad\text{ for every }\ t=0,\dots,T.$

  * -

Sparsistency for smoothly-changing DMRF. Suppose that the DMRF is smoothly
changing with parameters $(q,D),q\geq 1,D\geq 0$. Suppose that
$\Delta_{t}\leq\bar{\theta}_{t}/2$ and
$T\gtrsim\left(\mathcal{N}^{\mathrm{ker}}_{\mathrm{d}}\vee(\mathcal{E}^{\mathrm{ker}}_{\mathrm{d}}/\bar{\theta}_{t})^{4}\right)\log^{2}p$.
Then, the solution of ProxGL with $0<\gamma<1/(1+D)$ and a temporal
$\ell_{q}$-regularizer satisfies

$\displaystyle\mathrm{supp}\left(\widehat{\theta}_{t}\right)=\mathrm{supp}\left({\theta}^{\star}_{t}\right),\qquad\forall
t=0,\dots,T.$ (40)

  * -

Sparsistency for sparsely-changing DMRF. Suppose that the time-varying DMRF is
sparsely changing with sparsity parameter $D_{0}\geq 0$. Suppose that
$\Delta_{t}\leq\theta_{t}^{\min}/2$ and
$T\gtrsim\left(\mathcal{N}^{\mathrm{ker}}_{\mathrm{d}}\vee(\mathcal{E}^{\mathrm{ker}}_{\mathrm{d}}/{\theta}^{\min}_{t})^{4}\right)\log^{2}p$.
Then, with any choice of $0<\gamma<1$, the optimal solution of ProxGL with
temporal $\ell_{0}$-regularizer satisfies

$\displaystyle\mathrm{supp}\left(\widehat{\theta}_{t}\right)=\mathrm{supp}({\theta}^{\star}_{t}),$
$\displaystyle\qquad\forall t=0,\dots,T,$ (41)
$\displaystyle\mathrm{supp}\left(\widehat{\theta}_{t}-\widehat{\theta}_{t-1}\right)=\mathrm{supp}({\theta}^{\star}_{t}-{\theta}^{\star}_{t-1}),$
$\displaystyle\qquad\forall t=1,\dots,T.$

###### Proof.

The proof is similar to those of Theorem 3 and Proposition 5. The details are
omitted for brevity. ∎

The above propositions show that the previously imposed lower bounds on the
number of samples per time can be relaxed via kernel averaging. However, this
comes at the expense of increasing the error bound from
$\mathcal{O}(N_{t}^{-1/2})$ to $\mathcal{O}(T^{-1/4})$ (in this case, $T$
plays the role of sample size $N_{t}$).

##  7 Computational experiments

In this section, we test the proposed method in both synthetic and real
datasets. We point out that the single-parameter version of the proposed
algorithm was showcased in Fattahi and Gomez [2021], thus in this section we
seek to demonstrate the parametric version of the algorithm.

In §7.1 we consider synthetically generated instances of time-varying GMRFs
and compare the performance of our proposed method with two other well-known
techniques. In §7.2, we test our approach for time-varying DMRFs, using stock
market data for our case study.

###  7.1 Experiments with synthetic data

We consider randomly generated instances of sparsely-changing GMRFs. Given the
dimension $n$, the number of time periods $T$ and the number of observations
per time period $N_{t}$, the instances are constructed as follows. Initially,
for time period $t=0$, we construct a “true” precision matrix
$\bar{\Theta}_{t}\in\mathbb{R}^{n\times n}$ with exactly $3n$ off-diagonal
non-zero elements. Each non-zero off-diagonal element satisfies
$(\bar{\Theta}_{t})_{ij}=-0.4$, while

$(\bar{\Theta}_{t})_{ii}=1+\sum_{j\neq i}|(\bar{\Theta}_{t})_{ij}|$

to guarantee positive definiteness. Then, for all subsequent time periods, 4%
of the non-zero off-diagonal elements are set to $0$, and the same number of
previously zero elements are set to the value $-0.4$. Once the precision
matrices have been constructed, we generate for each period $N_{t}$ iid
samples $x~{}\sim\mathcal{N}(0,\bar{\Theta}_{t}^{-1})$ that are used to train
the models, and another $N_{t}$ samples to use for validation. Thus, at the
end of this process, we obtain two datasets, each containing $T\times N_{t}$
samples of dimension $n$. Finally, for each combination of parameters, we
repeat this process five times and report the averages over the five instances
generated with identical parameters.

####  7.1.1 Methods

We test three approaches:
$\bullet$ ProxGL with cross-validation We use the parametric algorithm
described in §5, which solves ProxGL for all values of the parameter $\gamma$.
We then use cross-validation to determine the best value of $\gamma$ as
follows. Given any fixed value of $\gamma$, the method produces a sequence of
estimated precision matrices $\\{\hat{\Theta}_{t}\\}_{t=0}^{T}$. Then, given
the sequence of samples in the validation set for time period $t$,
$\\{x_{t}^{(i)}\\}_{i=1}^{N_{t}}$, the negative log-likelihood that the
samples where generated from the estimated precision matrices is proportional
to

$-\frac{N_{t}}{2}\log(\text{det}(\hat{\Theta}_{t}))+\frac{1}{2}\sum_{i=1}^{N_{t}}\left(x_{t}^{(i)}\right)^{\top}\hat{\Theta}_{t}x_{t}^{(i)},$

and the overall negative log-likelihood can be obtained by summing across all
time periods. We select the value of $\gamma$ that minimizes the negative log-
likelihood. In our computations, we set the $\ell_{\infty}$ parameter
$\lambda=0.2$ ( this parameter naturally matches the data-generation process,
where parameters change by multiples of $0.4$), and the shrinkage parameter
$\nu=\nu_{0}\sqrt{(\log n)/(T\cdot N_{t})}$ for
$\nu_{0}\in\\{0.0,0.2,0.5,0.8,2.0\\}$. Observe that we could select $\nu_{0}$
using cross-validation as well. Instead, we show that while the two extreme
choices $\nu_{0}\in\\{0,2\\}$ (corresponding to no shrinkage, or excessive
shrinkage) are unsatisfactory, all other choices of $\nu_{0}$ outperform the
alternatives.
$\bullet$ TVGL with cross-validation The time varying Graphical Lasso (TVGL),
computed as (4), is a well-known regularized MLE approach for estimating the
sparsely-changing GMRFs [Hallac et al., 2017, Cai et al., 2018]. Similar to
ProxGL, we pick $\gamma_{1}$ and $\gamma_{2}$ via cross-validation, i.e., by
selecting the parameters that minimize the negative log-likelihood.
$\bullet$ L1E with cross-validation Consider an $\ell_{1}$ relaxation of our
proposed estimator ProxGL, where the $\ell_{0}$ regularizer in the objective
function is replaced with its $\ell_{1}$ relaxation. The resulted estimator
reduces to that of Yang et al. [2014] for $T=0$, and to that of Wang et al.
[2018] for $T=1$ and $\gamma=1$. Similar to ProxGL and TVGL, we fine-tune the
parameter $\gamma$ via cross-validation.

####  7.1.2 Metrics

We test the quality of the approaches in terms of the estimation error with
respect to the true parameters, as well as the ability to correctly recover
the sparsity pattern of the sequence of precision matrices and their changes.
The estimation error is computed as

$\texttt{error}=\sqrt{\frac{\sum_{t=0}^{T}\sum_{i=1}^{n}\sum_{j=i}^{n}\left(\Theta^{\star}_{t;ij}-\hat{\Theta}_{t;ij}\right)^{2}}{\sum_{t=0}^{T}\sum_{i=1}^{n}\sum_{j=i}^{n}{\Theta^{\star}_{t;ij}}^{2}}},$

where $\\{\Theta^{\star}_{t}\\}_{t=0}^{T}$ are the “true” precision matrices
of the process that generated the data, and $\\{\hat{\Theta}_{t}\\}_{t=0}^{T}$
are the estimates delivered by a given method. The ability to recover the true
sparsity pattern is computed using the

$\texttt{F1-score}=2\times\frac{\texttt{Recall}\times\texttt{Precision}}{\texttt{Recall}+\texttt{Precision}},$

where $\texttt{Recall}=\texttt{TP}/(\texttt{TP}+\texttt{FP})$,
$\texttt{Precision}=\texttt{TP}/(\texttt{TP}+\texttt{FN})$, and TP, FP and FN
denote respectively the number of true positives, false positives and false
negatives in the sequence of estimated precision matrices. We also use the
F1-score to evaluate the detection of changes in the sparsity patterns, where
the formulas are identical but TP, FP and FN refer to changes in the sparsity
pattern of the matrices from a time period to the next.

####  7.1.3 Performance

We test the performance of all methods by varying the number of samples per
time period. Specifically, we fix $n=50$, $T=10$ and let $N_{t}=n\kappa$ for
integer $1\leq\kappa\leq 20$. Table 3 (at the end of the paper) presents the
results for $\kappa=\\{1,5,10,15,20\\}$ and varying levels of threshold in the
approximate backward mapping. Note that $\nu_{0}=0$ corresponds to no
thresholding, whereas $\nu_{0}=2$ sets every off-diagonal entry of the sample
covariance matrix to zero. As can be seen in the table, ProxGL with
$\nu_{0}\in\\{0,2\\}$ is clearly inferior to other choices of the parameters,
so these values are excluded from our subsequent experiments. Figure 4 depicts
the performance of ProxGL, TVGL, and L1E for more values of $N_{t}/n$. We note
that both ProxGL and TVGL tend to perform better as the number of samples
$N_{t}$ increases (as expected). While TVGL results in a better estimation
error with few samples ($N_{t}=2n$), ProxGL performs better for all values of
$N_{t}\geq 8n$ and all parameters $\nu_{0}\in\\{0.2,0.5,0.8\\}$. In
particular, for the largest value $N_{t}=40n$, ProxGL decreases the estimation
error from 7.4% of TVGL to up to 5.3% (if $\nu_{0}=0.2$). We also observe that
both TVGL and ProxGL have good performance in terms of detecting the true
sparsity pattern of $\hat{\Theta}_{t}$. However, TVGL fails to correctly
identify changes in the underlying graphical model (with a F1-score of $0.2$
for most values of $N_{t}$), whereas ProxGL is able to do so accurately (with
F1-score above to $0.7$ for the largest values of $N_{t}$ considered). Method
L1E performs poorly according to the three metrics considered, and is inferior
to the other two approaches. Finally, we point out that no value of $\nu_{0}$
for ProxGL clearly dominates the others, with larger values resulting in
better performance when $N_{t}$ is small and smaller values performing best
with a large number of samples per time period.

We now briefly discuss computational times. TVGL is solved via the solver
MOSEK [ApS, 2019], which is one of the fastest available solvers for general-
purpose convex optimization problems. We highlight the fact that while faster
algorithms for static Graphical Lasso (such as QUIC [Hsieh et al., 2014] and
BIG-QUIC [Hsieh et al., 2013]) exist, these algorithms do not readily extend
to the time-varying settings. Moreover, L1E reduces to linear programming
(LP), for which fast algorithms (such as GUROBI [Gurobi Optimization, LLC,
2023], CPLEX [Cplex, 2009], and MOSEK) exist. However, these solvers are
incapable of solving LPs parametrically, thereby leading to larger runtimes
with cross-validation. In our experiments, we use MOSEK to solve L1E.

Our implementation of ProxGL is very fast: for the large-sample instance with
$N_{t}=40n$, it requires, on average, 11 seconds per instance on a single
thread. We point out that this time includes computing the backward mappings,
solving ProxGL for all values of $\gamma$, and printing to disc (a CSV file)
the estimated precision matrices for all possible values of $\gamma$. We also
highlight that the most expensive component is, in fact, printing to file,
which accounts for 9 out of the 11 seconds on average, while computing the
backward mappings and solving ProxGL is done in less than two seconds. On the
other hand, for the same large-sample instance with $N_{t}=40n$, TVGL with
cross-validation over 16 choices of the parameters $(\gamma_{1},\gamma_{2})$
requires on average 1600 seconds per instance on a single thread, which is at
least 145 times slower than ProxGL. Finally, L1E with cross-validation over 10
choices of the parameter $\gamma$ requires on average 61 seconds per instance
on a single thread, which is at least 5 times slower than ProxGL.

(a) F1-score: parameters (b) F1-score: differences between parameters (c)
Estimation error

Figure 4: Performance metrics between ProxGL, TVGL and L1E. The top row
depicts the ability of each method to recover the sparsity pattern and
identify changes in the underlying graph (top left: F1-score of the
parameters, top right: F1-score of changes in the parameters), and bottom
graphs depict the actual estimation error.

In summary, we see that while TVGL has a reasonable performance and can
identify the sparsity pattern of the precision matrices, ProxGL is able to
decrease further the estimation error and is substantially better at
identifying changes in the underlying graphs, and does so with only a fraction
of the computational cost.

####  7.1.4 Further experiments on computational times

We further discuss the computational times of solving ProxGL for larger
instances of GMRFs. To do so, we fix $T=10$ and $N_{t}=n/2$, and vary the
dimension $n$ (and thus the number of parameters to be estimated). The
computational times for DMRF and with varying $T$ are discussed in the next
section. The results are summarized in Table 1. It shows the dimension $n$ and
the total number of parameters to be estimated—computed as $n(n+1)T/2$.
Moreover, it shows the time (in seconds) required to compute the backward