# Directed Acyclic Graph Structure Learning from Dynamic Graphs

Shaohua Fan1, Shuyang Zhang1, Xiao Wang1,2, Chuan Shi1,2 Corresponding author.

###### Abstract

Estimating the structure of directed acyclic graphs (DAGs) of features
(variables) plays a vital role in revealing the latent data generation process
and providing causal insights in various applications. Although there have
been many studies on structure learning with various types of data, the
structure learning on the dynamic graph has not been explored yet, and thus we
study the learning problem of node feature generation mechanism on such
ubiquitous dynamic graph data. In a dynamic graph, we propose to
simultaneously estimate contemporaneous relationships and time-lagged
interaction relationships between the node features. These two kinds of
relationships form a DAG, which could effectively characterize the feature
generation process in a concise way. To learn such a DAG, we cast the learning
problem as a continuous score-based optimization problem, which consists of a
differentiable score function to measure the validity of the learned DAGs and
a smooth acyclicity constraint to ensure the acyclicity of the learned DAGs.
These two components are translated into an unconstraint augmented Lagrangian
objective which could be minimized by mature continuous optimization
techniques. The resulting algorithm, named GraphNOTEARS, outperforms baselines
on simulated data across a wide range of settings that may encounter in real-
world applications. We also apply the proposed approach on two dynamic graphs
constructed from the real-world Yelp dataset, demonstrating our method could
learn the connections between node features, which conforms with the domain
knowledge.

## Introduction

A Bayesian network (BN) is a probabilistic graphical model that represents a
set of variables and their conditional dependencies. It has been widely used
in machine learning applications (Pearl 2011; Ott, Imoto, and Miyano 2003;
Friedman, Geiger, and Goldszmidt 1997). The structure of a BN takes the form
of a directed acyclic graph (DAG) and provides a convenient and interpretable
output which is needed in today’s high-stake applications of artificial
intelligence, such as healthcare, finance and autonomous driving. The edges in
a DAG represent the directed generation relationships between variables (e.g.,
features) in a system. When these edges are not known based on prior
knowledge, one possible solution is to resort to DAG structure learning,
namely, learning the edges in a graphical model from the observed data.

Figure 1: A toy example of feature generation process on dynamic graph.

Existing approaches for DAG learning mostly focus on dealing with tabular
data, i.e., each sample is independently drawn from the same distribution (IID
data) (Spirtes et al. 2000; Neuberg 2003; Spirtes, Meek, and Richardson 2013;
Geiger and Heckerman 1994; Heckerman, Geiger, and Chickering 1995).
Nevertheless, in real-world scenarios, there usually exists associations
between samples, so the generation of features of certain samples may be
influenced by the other samples through the links between samples. Several
pioneer works have proposed constraint-based methods to learn DAGs from static
(i.e., equilibrium) graph data (Maier et al. 2010; Lee and Honavar 2016; Maier
et al. 2013). They usually test the conditional independencies of the
attributes of entities in graph data to form DAGs among those variables.

However, many real-world scenarios exhibit temporal information in graph data.
For example, as shown in Fig. 1, Bob’s friends, Lisa and Mary, went to eat
seafood and posted their recommendation and the meal time on the social media
at timestamp $t-1$. Bob viewed this information. When he has enough available
time at timestamp $t$ (e.g., several days later), he will choose to eat
seafood with high probability. Hence, the generation of Bob’s type value at
timestamp $t$ will be both determined by current available time and
recommendation from social network at previous timestamp. Another example is
that, the risk of ones to be infected by COVID-19 may be both determined by
current protection status (e.g., wearing a mask or keeping social distance)
and the ratio of neighborhood who has been vaccinated or infected at previous
timestamp two weeks ago. Actually, the current methods largely ignore modeling
the temporal interaction, so that the true data generation process could not
be revealed accurately.

When learning DAGs from dynamic graphs, there are two intractable challenges
we need to confront. First, a dynamic graph contains complex temporal
interactions between samples, so what kind of DAGs would reflect the
generation process of features in a dynamic graph? As the samples at each
timestamp will be generated based on the interactions from previous
timestamps, the learned DAG should model the generation process of each new
sample at each timestamp, considering the time-lagged interaction information.
Second, how to efficiently learn a DAG from complex evolutional graph data?
Compared with IID data, a dynamic graph contains both temporal and interaction
information, hence it has a more complicated data generation mechanism. It is
non-trivial to design a DAG learning method for dynamic graphs. Fortunately,
owing to the well-developed optimization techniques, it is possible to develop
a differentiable score function that could measure the validity of candidate
DAGs and resort to blackbox solvers to find the optimal DAG efficiently.

Particularly, to address these two challenges, we propose an effective score-
based approach for learning DAGs that could scale gracefully to dynamic graphs
with high-dimensional node features, called GraphNOTEARS. To solve the first
challenge, we propose to learn an intra-slice matrix to characterize
contemporaneous relationships between variables, and several inter-slice
matrices to characterize multi-step time-lagged graph influence on current
timestamp. Meanwhile, an acyclicity constraint is required to ensure the
acyclicity of the learned whole graph. As for the second challenge, we cast
the problem as a score-based optimization problem, and develop a least-squares
loss based score function. The score function leverages the temporal and
interaction information, as well as two kinds of learnable structural matrices
to reconstruct the data. With the smooth acyclicity constraint, we translate
the original unsolvable constraint problem into an unconstraint augmented
Lagrangian objective. The resulting program could be solved by the standard
second-order optimization schemes efficiently. The main contributions of this
paper are summarized as follows:

  * •

To our best knowledge, we first study the DAG learning problem on dynamic
graphs. Because of the ubiquity of dynamic graph data in real applications,
learning DAGs on such data could reveal the underlying feature generation
process, provide skeletons for possible Bayesian networks, and answer causal
questions, like the effect of various interventions. These applications are
very important for building explainable, robust, and generalized algorithms.

  * •

We develop a score-based learning method for simultaneously estimating the
structure and parameters of a sparse DAG from a dynamic graph. The resulting
method can be used to learn the relationships between variables of arbitrary
time-lagged order in a dynamic graph, without any implicit assumptions on the
underlying DAG topologies.

  * •

We conduct extensive simulation experiments with broad range settings which
may encounter in real world, validating the effectiveness of our approach in
revealing the feature generation mechanism of dynamic graphs. The experiments
on real-world datasets well demonstrate the rationality of the relationships
inferenced by GraphNOTEARS.

## Background and Related Works

A DAG $G$ is faithful with respect to a joint distribution $\mathcal{P}$ of a
set of variables if all and only the conditional independencies of variables
true in $\mathcal{P}$ are entailed by $G$ (Pearl 2014). The faithfulness
assumption enables one to recover $G$ from $\mathcal{P}$. Given samples $D$
from an unknown distribution corresponding to a faithful but unknown DAG,
structure learning refers to recovering the DAG from $D$.

Existing methods for DAG learning can be classified into constraint-based
methods and score-based methods. Most constraint-based DAG learning methods
(Spirtes et al. 2000; Neuberg 2003; Spirtes, Meek, and Richardson 2013) first
use conditional independence tests to find graph skeleton and then determine
the orientations of the edges up to the Markov equivalence class, which
usually contains DAGs that can be structurally diverse and may still have many
unoriented edges. Score-based methods (Geiger and Heckerman 1994; Huang et al.
2018; Hyvärinen and Smith 2013), on the other hand, define a score function to
find the best DAG that fits the given data. Unlike constraint-based methods
that assume faithfulness and identify only the Markov equivalence class, these
methods are able to distinguish between different DAGs in the same equivalence
class, owing to the additional assumptions on data distribution and/or
functional classes. However, due to the acyclicity constraint and
superexponential in the number of nodes of DAGs to search over, score-based
methods are computationally expensive. The recent work, NOTEARS (Zheng et al.
2018), expresses the acyclicity of a DAG by a smooth equality constraint under
the linear assumption, which makes it possible to formulate structure learning
as a smooth minimization problem subject to this equality constraint.
DYNOTEARS (Pamfil et al. 2020) extend NOTEARS to learn the DAG of time-series
data, which incorporates the temporal information into the score function.

The mainstream of existing methods is designed for the tabular samples, which
are independently identically drawn from the same distribution. However, in
many real-world settings, there may exist links between samples and the links
may influence the features generation process of samples. Maier et al. (Maier
et al. 2010) first extend the well-known PC algorithm to the relational
setting for learning causal relationships from relational data, called RPC.
Later, Maier et al. (Maier et al. 2013) demonstrate the lack of completeness
of RPC and introduce a sound and complete algorithm, named RCD. All of these
relational DAG learning algorithms are constraint-based methods and they could
only handle static graphs, ignoring temporal information. In this paper, we
propose an efficient score-based method that could model the temporal
interaction information.

## DAG Structure Learning on Dynamic Graphs

Figure 2: Illustration of intra-slice (solid lines) and inter-slice (dashed
lines) dependencies in a dynamic graph with $n=3$ samples and $d=3$ variables
at each timestamp and time-lagged order $p=2$. For clarity, we ignore the
edges that do not influence the variables $\mathbf{X}_{11}^{(t)}$.

In this section, we introduce the definition of dynamic graphs, the problem of
DAG structure learning on dynamic graphs as well as the proposed model:
GraphNOTEARS.

### Problem Formulation

######  Definition 1 (Dynamic Graph)

A dynamic graph is
$\mathcal{G}=\\{(\mathbf{X}^{(1)},\mathbf{A}^{(1)}),\cdots,(\mathbf{X}^{(T)},\mathbf{A}^{(T)})\\}$,
where $T$ is the total number of timestamps, tuple ($\mathbf{X}^{(T)}$,
$\mathbf{A}^{(T)}$) represents the graph at timestamp $T$,
$\mathbf{X}^{(T)}\in\mathbb{R}^{n\times d}$ is the matrix of node features and
$\mathbf{A}^{(T)}\in\mathbb{R}^{n\times n}$ is the adjacency matrix of nodes,
$n$ is the number of nodes111We assume the number of nodes at each timestamp
remains the same., $d$ is the number of node features (i.e., variables).

At each timestamp $t$, we assume that each feature of nodes is generated based
on the contemporaneous variables and time-lagged neighborhood variables. For
instance, as depicted in Fig. 2, the variable $\mathbf{X}_{11}^{(t)}$ of
sample $\mathbf{X}_{1}^{(t)}$ at timestamp $t$ is determined by the
contemporaneous variables $\mathbf{X}_{12}^{(t)}$ and $\mathbf{X}_{13}^{(t)}$
from the same sample222In this paper, we use node and sample interchangeably.
with coefficients $\mathbf{W}_{21}$ and $\mathbf{W}_{31}$, and the time-lagged
aggregated neighborhood variables $\hat{\mathbf{X}}_{.1}^{(t-1)}$ and
$\hat{\mathbf{X}}_{.1}^{(t-2)}$ from timestamp $t-1$ and $t-2$ with
coefficients $\mathbf{P}^{(t-1)}_{11}$ and $\mathbf{P}^{(t-2)}_{11}$,
respectively. We call these intra-slice and inter-slice dependencies,
respectively. One may argue that the generation of variables of samples at
current timestamp $t$ should also depend on the neighborhood samples at the
same timestamp. In this paper, we assume that neighborhood behaviour needs a
delay to influence ego nodes. The delay influence phenomenon is very common in
real-world applications. One example is that the user’s preference on
restaurant category may be influenced by their friends’ recommendation on Yelp
platform, and the user will go to the restaurant in a few days. Moreover, we
assume the interactions between samples are given at each timestamp, which is
very common in dynamic graph data (Rossi et al. 2020; Sankar et al. 2020). For
example, we could easily get the following relationship between users on Yelp
platform.

Moreover, we assume that in each timestamp, the effects of the time-lagged
influences from neighborhood on current timestamp could be at most $p$
timestamps, where $p\leq T$. We also make the stationary process assumption
(Hamilton 2020) that the generation process is fixed through time and is
identical for generating each timestamp in the time-series, which is very
common assumption in time-series data (Hamilton 2020). Without loss of any
generality, we propose to model the data generation process at timestamp $t$
with the structural vector autoregressive (SVAR) model (Demiralp and Hoover
2003; Kilian 2013):

$\displaystyle\mathbf{X}^{(t)}$
$\displaystyle=\mathbf{X}^{(t)}\mathbf{W}+\mathbf{\hat{A}}^{(t-1)}\mathbf{X}^{(t-1)}\mathbf{P}^{(t-1)}+\cdots$
(1)
$\displaystyle+\mathbf{\hat{A}}^{(t-p)}\mathbf{X}^{(t-p)}\mathbf{P}^{(t-p)}+\mathbf{Z},$

where $\mathbf{W}\in\mathbb{R}^{d\times d}$ and
$\mathbf{P}^{(t-i)}(i\in\\{1,\cdots,p\\})\in\mathbb{R}^{d\times d}$ represents
weighted adjacency matrices with nonzero entries corresponding to the intra-
slice and inter-slice edges, respectively.
$\mathbf{\mathbf{\hat{A}}}^{(t-i)}(i\in\\{1,\cdots,p\\})$ are normalized
adjacency matrices, which is computed by
$\mathbf{D}^{-\frac{1}{2}}(\mathbf{A}+\mathbf{I})\mathbf{D}^{-\frac{1}{2}}$,
$\mathbf{D}_{ii}=\sum_{j}\mathbf{A}_{ij}$, $\mathbf{I}$ is the identity matrix
of $\mathbf{A}$. The reason for $\mathbf{A}+\mathbf{I}$ is that the time-
lagged influence should not only from neighborhood but also from itself at
previous timestamps. And $\mathbf{Z}$ is a centered error matrix, where each
row is independent for each sample. The error variables could be Gaussian,
Exponential or others. The overall intuition of Eq. (1) would be that, the
generation of the value of a new variable $i$ of $\mathbf{X}^{(t)}$ depends on
two parts: parent variables from contemporaneous variables (i.e.,
$\mathbf{X}^{(t)}\mathbf{W}$) and parents variables from time-lagged
neighborhood variables (i.e.,
$\mathbf{\hat{A}}^{(t-1)}\mathbf{X}^{(t-1)}\mathbf{P}^{(t-1)}+\cdots+\mathbf{\hat{A}}^{(t-p)}\mathbf{X}^{(t-p)}\mathbf{P}^{(t-p)}$).
$\mathbf{X}^{(t)}\mathbf{W}$ represents the relationship between
contemporaneous variables. According to SVAR model, contemporaneous variables
usually exhibit a causal order, hence $\mathbf{W}$ represents the causal order
of contemporaneous variables and is acyclic, where $\mathbf{W}_{kj}$
represents the coefficient of the parent variable $k$ at timestamp $t$ on the
variable $j$ at the same timestamp. And $\mathbf{P}^{(t-i)}_{kj}$ represents
the coefficient of $k$-th aggregated node variables at timestamp $t-i$ on
$j$-th variable at timestamp $t$. For the detailed pseudocode of Eq. (1),
please refer to Appendix A.2 333Supplementary
material:https://drive.google.com/file/d/1S1pzEyy
C9kNL6s97yQxRt5lOLdO5L8sz/view?usp=sharing.

Let $\mathbf{M}=[\mathbf{X}^{(t-1)}|\cdots|\mathbf{X}^{(t-p)}]$ be the
$n\times pd$ matrix of time-lagged node features data,
$\mathbf{A}=[\mathbf{\hat{A}}^{(t-1)}|\cdots|\mathbf{\hat{A}}^{(t-p)}]$ be the
$n\times pn$ matrix of time-lagged interaction data, and
$\mathbf{P}=[\mathbf{P}^{(t-1)}|\cdots|\mathbf{P}^{(t-p)}]$ be the $pd\times
d$ matrix of inter-slice weights. We could rewrite Eq. (1) as following
compact form of structural equation model (SEM) (Hox and Bechger 1998):

$\mathbf{X}=\mathbf{X}\mathbf{W}+\mathbf{A}\boxtimes\mathbf{M}\mathbf{P}+\mathbf{Z},$
(2)

where
$\mathbf{A}\boxtimes\mathbf{M}=[\mathbf{\hat{A}}^{(t-1)}\mathbf{X}^{(t-1)}|\cdots|\mathbf{\hat{A}}^{(t-p)}\mathbf{X}^{(t-p)}]\in\mathbb{R}^{n\times
pd}$. This general formulation covers scenarios in which the time-lagged data
matrix $\mathbf{M}$ or $\mathbf{A}$ are not a contiguous sequence of time
slices (i.e., from $t-p$ to $t-1$). For example, someone usually meets their
friends in the weekend, hence one can include the lagged data matrix
$\mathbf{M}$ or $\mathbf{A}$ only those time points that have an impact on the
variables at timestamp $t$. Based on the generation formulation, we formulate
our target problem:

######  Problem 1 (DAG Structure Learning on Dynamic Graph)

Given a dynamic graph $\mathcal{G}$, the goal of DAG structure learning on the
dynamic graph is to estimate the generation process of node features, which
usually forms as a DAG, both considering the contemporaneous effect from each
node itself and neighborhood effects from time-lagged graphs.

Hence, given the data $\mathbf{X}$, $\mathbf{M}$ and $\mathbf{A}$, the goal of
this paper is to estimate weighted adjacency matrices $\mathbf{W}$ and
$\mathbf{P}$, which could characterize the node feature generation process in
a dynamic graph.

### The Proposed Model: GraphNOTEARS

An SEM could be found through minimizing the least-squares (LS) loss (Zheng et
al. 2018). The statistical properties of the LS loss in scoring DAGs have been
extensively studied: The minimizer of the least-squares loss provably recovers
a true DAG with high probability on finite-samples and in high-dimensions
(Aragam, Amini, and Zhou 2015; Loh and Bühlmann 2014). Inspired by this, we
propose to estimate $\mathbf{W}$ and $\mathbf{P}$ by minimizing the following
LS loss:

$\mathcal{L}(\mathbf{W},\mathbf{P})=\frac{1}{2n}||\mathbf{X}-\mathbf{X}\mathbf{W}-\mathbf{A}\boxtimes\mathbf{M}\mathbf{P}||^{2}_{F}.$
(3)

Moreover, the edges in $\mathbf{P}$ go only forward in time and thus they do
not create cycles. To ensure that the whole graph is acyclic, it thus suffices
to require that $\mathbf{W}$ is acyclic. In this work, we utilize acyclic
constraint proposed by NOTEARS (Zheng et al. 2018) to ensure the acyclicity of
learned $\mathbf{W}$, which states that: a directed graph $G$ with binary
adjacency matrix $\mathbf{W}$ is acyclic if and only if:

$h(\mathbf{W}):=\text{trace}(e^{\mathbf{W}\circ\mathbf{W}})-d=0,$ (4)

where $e^{\mathbf{W}\circ\mathbf{W}}$ is the matrix exponential of
${\mathbf{W}\circ\mathbf{W}}$, and $\circ$ denotes the Hadamard product of two
matrices. To enforce the sparsity of $\mathbf{W}$ and $\mathbf{P}$, we also
introduce $\ell_{1}$ penalties in the objective function. The overall
optimization problem is:

$\displaystyle\min_{\mathbf{W},\mathbf{P}}f(\mathbf{W},\mathbf{P})$ (5) s.t.
$\displaystyle
h(\mathbf{W}):=\text{trace}(e^{\mathbf{W}\circ\mathbf{W}})-d=0,$ with
$\displaystyle
f(\mathbf{W},\mathbf{P})=\frac{1}{2n}||\mathbf{X}-\mathbf{X}\mathbf{W}-\mathbf{A}\boxtimes\mathbf{M}\mathbf{P}||^{2}_{F}$
$\displaystyle+\lambda_{\mathbf{W}}||\mathbf{W}||_{1}+\lambda_{\mathbf{A}}||\mathbf{P}||_{1},$

where the $||\cdot||_{1}$ represents the element-wise $\ell_{1}$ norm, and
$\lambda_{\mathbf{W}}$ and $\lambda_{\mathbf{A}}$ represent the coefficients
of $||\mathbf{W}||_{1}$ and $||\mathbf{P}||_{1}$, respectively. As we can see,
given $\mathbf{X}$, $\mathbf{M}$ and $\mathbf{A}$, though minimizing this
objective, our model could simultaneously recover the contemporaneous
dependencies $\mathbf{W}$ and neighborhood dependencies $\mathbf{P}$ from
data, as well as ensure the acyclicity of the learned graph.

### Optimization

The above optimization is a standard equality-constrained program (ECP). We
translate the problem to an unconstrained problem with the following smooth
augmented Lagrangian objective

$\min_{\mathbf{W},\mathbf{P}}f(\mathbf{W},\mathbf{P})+\frac{\rho}{2}h(\mathbf{W})^{2}+\alpha
h(\mathbf{W}).$ (6)

The resulting problem can be solved using efficient solvers such as L-BFGS-B
(Zhu et al. 1997), and the update strategy of $\rho$ and $\alpha$ is the same
as (Zheng et al. 2018). To reduce false discoveries (Zhou 2009), we threshold
the edge weights of $\mathbf{W}$ and $\mathbf{P}$ via hard thresholds
$\tau_{\mathbf{W}}$ and $\tau_{\mathbf{P}}$: After obtaining a stationary
point $\mathbf{W}$ and $\mathbf{P}$, given a fixed threshold
$\mathbf{\tau_{W}}>0$ and $\mathbf{\tau_{P}}>0$, set any weights smaller than
$\mathbf{\tau_{W}}$ and $\mathbf{\tau_{P}}$ in absolute value to zero, termed
as $\widetilde{\mathbf{W}}$ and $\widetilde{\mathbf{P}}$, and their binary
version termed as $\mathbf{\hat{W}}$ and $\mathbf{\hat{P}}$, respectively.

### Discussion

Here we discuss the identifiability, some limitations and possible extensions
of GraphNOTEARS.

Identifiability Identifiability is the key research problem of SVAR models of
the econometrics literature (Kilian 2013). Identifiability of structure
learning on time-series data has been discussed in (Pamfil et al. 2020) by
using the conclusions of SVAR model. As a dynamic graph could be viewed as a
special time-series data, where $\mathbf{A}\boxtimes\mathbf{M}$ in Eq. (2) is
the aggregated time-lagged features, we assume the conditions for the
identifiability in (Pamfil et al. 2020) will be also held in our model, i.e.,
the error $\mathbf{Z}$ could be drawn from non-Gaussian or standard normal
distribution. Thus, $\mathbf{W}$ and $\mathbf{P}$ of our model are
identifiable under reasonable conditions.

Assumptions To be simplified, we have assumed that the structure of the
process of variable generation is fixed across time and is identical for all
timestamps. Based on this assumption, when long time-series data is available,
our model could easily utilize such long time-series in the objective function
Eq. (5) by extending data matrices (i.e., $\mathbf{X}$, $\mathbf{M}$ and
$\mathbf{A}$) to tensors and keeping the parameter matrices $\mathbf{W}$ and
$\mathbf{P}$ unchanged. And all our experiments are based on this extension.
This stationary process assumption is a very common assumption in time-series
data (Hamilton 2020) and could be relaxed in several ways. We could allow the
directed dependency structure of data to vary smoothly over time (Song, Kolar,
and Xing 2009) or have discrete change points (Grzegorczyk and Husmeier 2011).

Nonlinear relationship As a very beginning work of structure learning on
dynamic graphs, we follow previous works on structure learning (Zheng et al.
2018; Pamfil et al. 2020) that first consider the linear scenarios. Note that
linear assumption is made purely for simplicity, so that our paper could focus
on the most salient temporal and network aspects of this problem. Inspired by
the GNNs (Wu et al. 2020; Fan et al. 2019, 2020, 2021, 2022b, 2022a) and the
nonlinear structure learning methods (Zheng et al. 2020; Yu et al. 2019;
Lachapelle et al. 2019), we could model the nonlinear effects of neighbors by
GNNs.

An important feature of GraphNOTEARS is its simplicity, both in terms of
formulating an objective function and optimizing it. Nevertheless, the
proposed model is general enough to be extended to more complicated scenarios.

## Experiments

It is notoriously hard to obtain the ground truth of causal structure because
it is difficult to obtain the underlying data generation process of real-world
problems. To validate the effectiveness of our method, in this section, we
follow the setting in (Zheng et al. 2018; Pamfil et al. 2020; Maier et al.
2010, 2013), which conduct extensive experiments on synthetic data with known
generating mechanisms to simulate real-world scenarios. 444Code and data:
https://github.com/googlebaba/GraphNOTEARS.

Datasets. To validate the effectiveness of GraphNOTEARS against existing
approaches, we simulate data according to the SEM from Eq. (2). We need three
steps to this process: 1) generating the weighted graphs
$\mathcal{G}_{\mathbf{W}}$ and $\mathcal{G}_{\mathbf{P}}$, and adjacency
matrix $\mathbf{A}$; 2) generating data matrices $\mathbf{X}$ and $\mathbf{M}$
based on $\mathcal{G}_{\mathbf{W}}$ and $\mathcal{G}_{\mathbf{P}}$; 3) running
all algorithms on all or partial of $\mathbf{X}$, $\mathbf{M}$ and
$\mathbf{A}$ based on whether the model considers this kind of information and
computing metrics respectively. Particularly, following (Pamfil et al. 2020),
we use either the Erdős-Rényi (ER) model (Newman 2018) or the Barabási-Albert
(BA) model (Barabási and Albert 1999) to generate intra-slice graphs
$\mathcal{G}_{\mathbf{W}}$. And for inter-slice graph
$\mathcal{G}_{\mathbf{P}}$, we use ER model or Stochastic Block Model (SBM)
(Newman 2018). These graph generation models could simulate real-world
variable generation processes, like “rich get richer”, “preferential
attachment” and “cluster effect”. To get weighted intra- and inter-slice
matrices $\mathbf{W}$ and $\mathbf{P}$, we sample weights uniformly at random
from $[-2,-0.5]\cup[0.5,2]$. For the generation of $\mathbf{A}$, we connect
each pair of samples with 0.1 probability. Once given $\mathbf{W}$,
$\mathbf{P}$ and $\mathbf{A}$, we use the SEM from Eq. (2) to generate a data
matrix $\mathbf{X}$ of size $n\times d$. In particular, we first generate the
variables of $\mathbf{X}$ and $\mathbf{M}$ based on the sorted topological
order of $\mathbf{W}$ same as (Zheng et al. 2018). And we generate current
timestamp observation according to:
$\mathbf{X}=\mathbf{X}\mathbf{W}+\mathbf{A}\boxtimes\mathbf{M}\mathbf{P}+\mathbf{Z}$.
For the noise term $\mathbf{Z}$, we utilize Gaussian noise and Exponential
noise. Moreover, to compare GraphNOTEARS against baselines with a wide range
of sample sizes and the number of variables, we vary the sample size
$n\in\\{100,200,300,500\\}$, the number of variables $d\in\\{5,10,20,30\\}$ at
each timestamp, and the length of time-series $T$ is set as 7 for all
experiments. A detailed introduction to the data generation process is in
Appendix A.

Baselines. Because we study a new problem, there is no baseline specially
designed for this problem. We compare the following two alternatives that
could deal with the problem in an indirect way.

  * •

NOTEARS (Zheng et al. 2018)+LASSO: This is a two-step approach. We use static
NOTEARS to estimate $\mathbf{W}$, and use Lasso regression to estimate
$\mathbf{P}$, independently.

NOTEARS:
$\mathcal{L}(\mathbf{W})=\frac{1}{2n}||\mathbf{X}-\mathbf{X}\mathbf{W}||^{2}_{F}+\lambda_{\mathbf{W}}||\mathbf{W}||_{1}.\text{s.t.,}\mathbf{W}\text{
is acyclic.}$

LASSO:
$\mathcal{L}(\mathbf{P})=\frac{1}{2n}||\mathbf{X}-\mathbf{A}\boxtimes\mathbf{M}\mathbf{P}||^{2}_{F}+\lambda_{\mathbf{P}}||\mathbf{P}||_{1}.$

  * •

DYNOTEARS (Pamfil et al. 2020): This is an extension of NOTEARS on time
series. Compared with our method, it ignores the interactions between samples
(i.e., $\mathbf{A}$).

Objective:
$\mathcal{L}(\mathbf{W},\mathbf{P})=\frac{1}{2n}||\mathbf{X}-\mathbf{X}\mathbf{W}-\mathbf{M}\mathbf{P}||^{2}_{F}+\lambda_{\mathbf{W}}||\mathbf{W}||_{1}+\lambda_{\mathbf{P}}||\mathbf{P}||_{1}.\text{s.t.,}\quad\mathbf{W}\text{
is acyclic.}$

For fully utilizing multiple-step graph series, these methods are also
extended to tensor versions.

Figure 3: Example results for data with Gaussian noise, $n=500$ samples, $d=5$
variables at each timestamp, $T=7$ time-series, and $p=2$ time-lagged graph
effect order. Our algorithm recovers weights that are closer to the ground
truth than baselines.

Metrics. We evaluate the learned binary intra-slice $\mathbf{\hat{W}}$ and
inter-slice matrices $\mathbf{\hat{P}}$ separately by two common graph
metrics: F1-score and Structural Hamming Distance (SHD) (Zheng et al. 2018).

Experimental setup. For all methods, we set hyperparameters
$\lambda_{\mathbf{W}}=\lambda_{\mathbf{P}}=0.01$. For the weight thresholds,
following (Zheng et al. 2018), we choose
$\tau_{\mathbf{W}}=\tau_{\mathbf{P}}=0.3$ for all the methods. The relative
ranking of the three methods is not sensitive to the weight thresholds
according to our observation. For utilizing multiple-step graph series, we use
the first $T-1$ timestamps to predict last $T-p$ timestamps, where $p$ time-
lagged influence order is considered at each timestamp. For all experiments,
we utilize 5 different random seeds to generate different datasets and
initialize models, and report the mean value and 95% confidence interval.

Figure 4: F1 scores (higher is better) for different noise models (Gaussian,
Exponential) and different sample sizes ($n\in[100,200,300,500]$). The length
of time-series is 7 and we consider 1-step time-lagged neighbor influence
here. Each panel contains results for two different choices of intra-slice
graphs (columns) and inter-slice graphs (rows). Every marker corresponds to
the mean performance across 5 algorithm runs, where each on a different
simulated dataset, and shade area means the 95% confidence interval.
Continuous and dashed lines represent F1 scores for intra-slice and inter-
slice edges, respectively.

### Performance Evaluation

We start by illustrating the estimated weighted matrices of GraphNOTEARS and
baselines with the ground truth in Fig. 3. The evaluation data is generated
with Gaussian noise, $n=500$ samples, $d=5$ variables, $T=7$ time-series, and
$p=2$ time-lagged graph effect order. And the intra-slice and inter-slice are
both generated with ER graph. As Fig. 3 shows, our estimated weights are much
closer to the true weights for both $\mathbf{W}$, $\mathbf{P}^{(1)}$ and
$\mathbf{P}^{(2)}$ compared with baselines, where the intra-slice $\mathbf{W}$
is recovered perfectly (F1 score =1.0) and the inter-slice $\mathbf{P}^{(1)}$
and $\mathbf{P}^{(2)}$ are also better than baselines with a large margin. The
results validate that our model is a general framework that could estimate the
DAG of a dynamic graph with multiple timestamps and could simultaneously learn
arbitrary time-lagged order influence. However, both NOTEARS+LASSO and
DYNOTEARS could not achieve satisfying results. The reasons could be that
NOTEARS+LASSO is a two-stage method that could not simultaneously consider the
contemporaneous and time-lagged interaction influence to generate variables,
while our model jointly models these two factors into a unified framework.
Furthermore, DYNOTEARS only takes the time-lagged series information into
consideration and ignores the interaction information. As we can see, if the
generation of variables are determined by the neighborhood, it is necessary to
incorporate the graph information into the model.

We present the F1-score and SHD results on the full setting in Fig. 4 and Fig.
5, respectively. Note that, for simplicity, we set time-lagged graph order
$p=1$ here. From the results, we have the following observations: (1)
GraphNOTEARS is the best algorithm. In most cases, for inter-slice graph
$\mathbf{W}$, our model nearly recovers the graph perfectly (F1
score$\approx$1.0 and SHD$\approx$0). For the intra-slice graph $\mathbf{P}$,
our method also achieves satisfying results, outperforming baselines with a
large margin. The phenomenon demonstrates that GraphNOTEARS could learn DAGs
from dynamic graphs, owing to the fact that it could comprehensively consider
complex information. (2) Overall, with the number of variables increasing,
especially in insufficient samples scenario, e.g., $n=100$, all the methods
suffer from the degradation of performance. However, our model performs stable
and still outperforms baselines with a large margin, indicating our model
could handle the challenging high-dimensional scenario well. (3) No matter in
what kind of underlying graphs or noise-type scenarios, our model could all
achieve promising results. It indicates our model has the potential to deal
with various scenarios which may encounter in real-world applications. (4) For
SHD results, GraphNOTEARS requires less modification of edges to reach ground
truth in all settings, further validating the effectiveness of the proposed
method.

Figure 5: SHD scores. Illustrations are the same as Fig. 4.

### Sensitivity of Threshold

We investigate the effect of threshold of $\mathbf{W}$ and $\mathbf{P}$ on the
results in Fig. 6. We present the estimated edge weights of
$\widetilde{\mathbf{W}}$ or $\widetilde{\mathbf{P}}$ in decreasing order in
$n=500$ and $n=100$ scenarios. As the ground-truth graphs are usually very
sparse, we can observe that most of the edge weights are equal or close to
zero as expected. The remaining question is how to choose a threshold that
separates out noise (near zero) from signals (away from zero) so that we can
obtain a DAG with the best performance. With enough samples ($n=500$), one can
often notice a sudden change in the weight distribution. With insufficient
samples ($n=100$), the turning point is less clear, and the optimal choice
that balances between recall and false discovery that depends on the specific
settings. Nonetheless, the predictive performance is less sensitive to the
threshold value as one can see there is a slope of the decrease in the weights
before getting close to zero.

(a) $\widetilde{\mathbf{W}},n=500$

(b) $\widetilde{\mathbf{P}},n=500$

(c) $\widetilde{\mathbf{W}},n=100$

(d) $\widetilde{\mathbf{P}},n=100$

Figure 6: Illustration of the effect of the threshold. In each subfigure, we
plot the sorted weights of $\widetilde{\mathbf{W}}$ and
$\widetilde{\mathbf{P}}$ in decreasing order.

## Application on Real-world Datasets

We consider applying GraphNOTEARS on two real-world dynamic graphs constructed
from Yelp dataset (Luca 2016). (Anderson and Magruder 2012) introduced a toy
SCM that embeds causal knowledge for the Yelp example. That is, there are
three random variables, i.e., the restaurant category $C$, Yelp star rating
$S$, and customer flow $Y$. There exist three directed edges that represent
the three causal relationships between variables: (1) Restaurant category
influences its Yelp rating. For example, the average rating of fast-food
restaurants is lower than that of high-end seafood restaurants. (2) Restaurant
category also influences its customer flow. For example, the average customer
flow of high-end restaurants is lower than fast food. (3) Yelp rating of a
restaurant influences its customer flow. According to this, we construct two
dynamic graphs, i.e., the user graph and business graph, where its node
features are these three variables. Particularly, we construct a user graph
based on whether two users are friends on the Yelp platform. Then we take the
time lag as one month and calculate the average category, the average Yelp
star rating, and the average customer flow of the restaurants they have
visited in this month as the node features. Here we consider 1-step time-
lagged graph information. As the user’s taste may be influenced by their
friends, the generation of users’ features should consider their friends’
influence. For example, if a user’s friend posts a positive review on this
restaurant, the user will have a larger possibility to visit this restaurant.
For the business graph, as the same category of restaurants may have the
“effect of agglomeration” to influence each other, we add edges between the
restaurants which have a similar category and close distance. Then we
calculate the variables of the restaurants same as the user graph.

We apply GraphNOTEARS and DYNOTEARS on the constructed graphs and obtain the
binary DAG via 0.1 threshold, shown in Fig. 7. For both two dynamic graphs,
the estimated relationships of the intra-slice matrix $\hat{\mathbf{W}}$
discovered by GraphNOTEARS coincides with our prior knowledge, i.e., the three
black directed edges among variables in Fig. 7(a)(c). However, DYNOTEARS could
only discover partial edges, e.g., missing $C\rightarrow F$ in Fig. 7(b), and
$C\rightarrow S$ and $C\rightarrow F$ in Fig. 7(d). And we find that there are
strong correlations between the same type of variables as illustrated in
inter-slice $\hat{\mathbf{P}}$ (e.g., the directed edge between the aggregated
neighborhood category (orange S) with self category (blue S)), which could be
explained by the homophily influence of the graph. Overall, our model could
discover an explainable DAG on real-world dynamic graphs.

(a) GraphNOTEARS on user dynamic graph.

(b) DYNOTEARS on user dynamic graph.

(c) GraphNOTEARS on business dynamic graph.

(d) DYNOTEARS on business dynamic graph.

Figure 7: Estimated DAG on real-world Yelp dataset.

## Conclusion

In this paper, we first study a new DAG learning diagram on dynamic graphs,
which plays a vital role in understanding the node features generation
mechanism. To handle such complex data, we propose a score-based DAG method to
learn both intra-slice and inter-slice dependencies between variables
simultaneously, considering both temporal and interaction information. The
resulting method could deal with such a complex problem efficiently and has
the potential for more complicated settings. Extensive experiments on both
simulated and real-world datasets well demonstrate the effectiveness of the
proposed method.

## Acknowledgments

This work is supported in part by the National Natural Science Foundation of
China (No. U20B2045, 62192784, 62172052, 62002029, 62172052, U1936014).

## References

  * Anderson and Magruder (2012) Anderson, M.; and Magruder, J. 2012.  Learning from the crowd: Regression discontinuity estimates of the effects of an online review database.  _The Economic Journal_ , 122(563): 957–989.
  * Aragam, Amini, and Zhou (2015) Aragam, B.; Amini, A. A.; and Zhou, Q. 2015.  Learning directed acyclic graphs with penalized neighbourhood regression.  _arXiv:1511.08963_.
  * Barabási and Albert (1999) Barabási, A.-L.; and Albert, R. 1999.  Emergence of scaling in random networks.  _science_ , 286(5439): 509–512.
  * Demiralp and Hoover (2003) Demiralp, S.; and Hoover, K. D. 2003.  Searching for the causal structure of a vector autoregression.  _Oxford Bulletin of Economics and statistics_ , 65: 745–767.
  * Fan et al. (2022a) Fan, S.; Wang, X.; Mo, Y.; Shi, C.; and Tang, J. 2022a.  Debiasing Graph Neural Networks via Learning Disentangled Causal Substructure.  _NeurIPS_.
  * Fan et al. (2021) Fan, S.; Wang, X.; Shi, C.; Cui, P.; and Wang, B. 2021.  Generalizing Graph Neural Networks on Out-Of-Distribution Graphs.  In _arXiv preprint arXiv:2111.10657_.
  * Fan et al. (2022b) Fan, S.; Wang, X.; Shi, C.; Kuang, K.; Liu, N.; and Wang, B. 2022b.  Debiased Graph Neural Networks with Agnostic Label Selection Bias.  _IEEE Transactions on Neural Networks and Learning Systems_.
  * Fan et al. (2020) Fan, S.; Wang, X.; Shi, C.; Lu, E.; Lin, K.; and Wang, B. 2020.  One2multi graph autoencoder for multi-view graph clustering.  In _WWW_ , 3070–3076.
  * Fan et al. (2019) Fan, S.; Zhu, J.; Han, X.; Shi, C.; Hu, L.; Ma, B.; and Li, Y. 2019.  Metapath-guided heterogeneous graph neural network for intent recommendation.  In _SIGKDD_ , 2478–2486.
  * Friedman, Geiger, and Goldszmidt (1997) Friedman, N.; Geiger, D.; and Goldszmidt, M. 1997.  Bayesian network classifiers.  _Machine learning_ , 29(2): 131–163.
  * Geiger and Heckerman (1994) Geiger, D.; and Heckerman, D. 1994.  Learning gaussian networks.  In _Uncertainty Proceedings 1994_ , 235–243. Elsevier.
  * Grzegorczyk and Husmeier (2011) Grzegorczyk, M.; and Husmeier, D. 2011.  Non-homogeneous dynamic Bayesian networks for continuous data.  _Machine Learning_ , 83(3): 355–419.
  * Hamilton (2020) Hamilton, J. D. 2020.  _Time series analysis_.  Princeton university press.
  * Heckerman, Geiger, and Chickering (1995) Heckerman, D.; Geiger, D.; and Chickering, D. M. 1995.  Learning Bayesian networks: The combination of knowledge and statistical data.  _Machine learning_ , 20(3): 197–243.
  * Hox and Bechger (1998) Hox, J. J.; and Bechger, T. M. 1998.  An introduction to structural equation modeling.
  * Huang et al. (2018) Huang, B.; Zhang, K.; Lin, Y.; Schölkopf, B.; and Glymour, C. 2018.  Generalized score functions for causal discovery.  In _SIGKDD_.
  * Hyvärinen and Smith (2013) Hyvärinen, A.; and Smith, S. M. 2013.  Pairwise likelihood ratios for estimation of non-Gaussian structural equation models.  _Journal of Machine Learning Research_ , 14(Jan): 111–152.
  * Kilian (2013) Kilian, L. 2013.  Structural vector autoregressions.  In _Handbook of research methods and applications in empirical macroeconomics_. Edward Elgar Publishing.
  * Lachapelle et al. (2019) Lachapelle, S.; Brouillard, P.; Deleu, T.; and Lacoste-Julien, S. 2019.  Gradient-based neural dag learning.  _arXiv preprint arXiv:1906.02226_.
  * Lee and Honavar (2016) Lee, S.; and Honavar, V. 2016.  On learning causal models from relational data.  In _AAAI_.
  * Loh and Bühlmann (2014) Loh, P.-L.; and Bühlmann, P. 2014.  High-dimensional learning of linear causal networks via inverse covariance estimation.  _The Journal of Machine Learning Research_ , 15(1): 3065–3105.
  * Luca (2016) Luca, M. 2016.  Reviews, reputation, and revenue: The case of Yelp. com.  _Harvard Business School NOM Unit Working Paper_ , (12-016).
  * Maier et al. (2013) Maier, M.; Marazopoulou, K.; Arbour, D.; and Jensen, D. 2013.  A sound and complete algorithm for learning causal models from relational data.  _arXiv preprint arXiv:1309.6843_.
  * Maier et al. (2010) Maier, M.; Taylor, B.; Oktay, H.; and Jensen, D. 2010.  Learning causal models of relational domains.  In _AAAI_ , volume 24.
  * Neuberg (2003) Neuberg, L. G. 2003.  Causality: models, reasoning, and inference, by judea pearl, cambridge university press, 2000.  _Econometric Theory_ , 19(4): 675–685.
  * Newman (2018) Newman, M. 2018.  _Networks_.  Oxford university press.
  * Ott, Imoto, and Miyano (2003) Ott, S.; Imoto, S.; and Miyano, S. 2003.  Finding optimal models for small gene networks.  In _Biocomputing 2004_ , 557–567. World Scientific.
  * Pamfil et al. (2020) Pamfil, R.; Sriwattanaworachai, N.; Desai, S.; Pilgerstorfer, P.; Georgatzis, K.; Beaumont, P.; and Aragam, B. 2020.  Dynotears: Structure learning from time-series data.  In _AISTATS_.
  * Pearl (2011) Pearl, J. 2011.  Bayesian networks.
  * Pearl (2014) Pearl, J. 2014.  _Probabilistic reasoning in intelligent systems: networks of plausible inference_.  Elsevier.
  * Rossi et al. (2020) Rossi, E.; Chamberlain, B.; Frasca, F.; Eynard, D.; Monti, F.; and Bronstein, M. 2020.  Temporal graph networks for deep learning on dynamic graphs.  _arXiv preprint arXiv:2006.10637_.
  * Sankar et al. (2020) Sankar, A.; Wu, Y.; Gou, L.; Zhang, W.; and Yang, H. 2020.  Dysat: Deep neural representation learning on dynamic graphs via self-attention networks.  In _Proceedings of the 13th international conference on web search and data mining_ , 519–527.
  * Song, Kolar, and Xing (2009) Song, L.; Kolar, M.; and Xing, E. 2009.  Time-varying dynamic bayesian networks.  In _NeurIPS_.
  * Spirtes et al. (2000) Spirtes, P.; Glymour, C. N.; Scheines, R.; and Heckerman, D. 2000.  _Causation, prediction, and search_.  MIT press.
  * Spirtes, Meek, and Richardson (2013) Spirtes, P. L.; Meek, C.; and Richardson, T. S. 2013.  Causal inference in the presence of latent variables and selection bias.
  * Wu et al. (2020) Wu, Z.; Pan, S.; Chen, F.; Long, G.; Zhang, C.; and Philip, S. Y. 2020.  A comprehensive survey on graph neural networks.  _IEEE transactions on neural networks and learning systems_ , 32(1): 4–24.
  * Yu et al. (2019) Yu, Y.; Chen, J.; Gao, T.; and Yu, M. 2019.  DAG-GNN: DAG structure learning with graph neural networks.  In _ICML_.
  * Zheng et al. (2018) Zheng, X.; Aragam, B.; Ravikumar, P.; and Xing, E. P. 2018.  Dags with no tears: Continuous optimization for structure learning.  In _NeurIPS_.
  * Zheng et al. (2020) Zheng, X.; Dan, C.; Aragam, B.; Ravikumar, P.; and Xing, E. 2020.  Learning sparse nonparametric dags.  In _International Conference on Artificial Intelligence and Statistics_ , 3414–3425. PMLR.
  * Zhou (2009) Zhou, S. 2009.  Thresholding procedures for high dimensional variable selection and statistical estimation.  _NeurIPS_.
  * Zhu et al. (1997) Zhu, C.; Byrd, R. H.; Lu, P.; and Nocedal, J. 1997.  Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound-constrained optimization.  _ACM Transactions on mathematical software (TOMS)_ , 23(4): 550–560.


i.i.d. noise draws $\epsilon(\xi)$ are additive, zero mean, with bounded
variance $\sigma^{2}$ for all $\xi$. We assume we have fixed estimates
$\hat{L_{1}}$ and $\hat{\sigma}$ with $K_{1}>0$ and $K_{2}>0$ as in (22). By
fixing the step size as $\hat{h}_{\tilde{\mathcal{A}}}$ in (97), the
approximately active smoothing parameter $\hat{\mu}_{\tilde{\mathcal{A}}}^{*}$
in (97) minimizes the error between the gradient oracle in Algorithm 5 and the
true directional derivative of $f$ in the direction
$u_{\tilde{\mathcal{A}}}^{(k)}$ in the $\tilde{j}$-dimensional AS
$\tilde{\mathcal{A}}$. That is, $\mathcal{E}(\mu)$ in (21) (with
$u=u_{\tilde{\mathcal{A}}}^{(k)}$) is minimized by the choice
$\mu=\hat{\mu}_{\tilde{\mathcal{A}}}^{*}$. In particular, we have the bound

(98)
$\mathbb{E}_{u_{\tilde{\mathcal{A}}}^{(k)},\xi_{1},\xi_{2}}\left(\mathcal{E}^{\tilde{\mathcal{A}}}\left(\mu^{*}_{\tilde{\mathcal{A}}}\right)\right)\leq\frac{K_{1}+K_{2}}{\sqrt{2K_{1}K_{2}}}\sigma
L_{1}\sqrt{\tilde{j}(\tilde{j}+6)^{3}}.$

Proof: The proof is identical to the proof of STARS Theorem 4.3 (Modified) and
ASTARS Corollary 2 (making the usual substitutions), until we formulate (as in
(32))

(99)
$\mathbb{E}_{u,\xi_{1},\xi_{2}}(\mathcal{E}^{\tilde{\mathcal{A}}}(\mu))\leq\frac{1}{\mu^{2}}\mathbb{E}_{u}\left(2\sigma^{2}||u||^{2}+\frac{L_{1}^{2}}{4}\mu^{4}||u||^{6}\right)$

and proceed to bound the right hand side. Applying ASTARS Proposition 1, we
have $(u_{\tilde{\mathcal{A}}}^{(k)})_{p}\sim N(0,1)$ for $p=1,\ldots,P$.
Taking $u=u_{\tilde{\mathcal{A}}}^{(k)}$ and noting
$u_{\tilde{\mathcal{A}}}^{(k)}\in\tilde{\mathcal{A}}$ (and
$\dim\tilde{\mathcal{A}}=\tilde{j}$), we apply the bounds on the moments
$M_{p}$ of $||u||^{p}$ given in (11). Using $p=2$ and $p=6$ – but replacing
the AS dimension $j$ with $\tilde{j}$ – we have

(100)
$\mathbb{E}_{u_{\tilde{\mathcal{A}}}^{(k)},\xi_{1},\xi_{2}}(\mathcal{E}(\mu))\leq(2\sigma^{2}\tilde{j})\frac{1}{\mu^{2}}+\left(\frac{L_{1}^{2}(\tilde{j}+6)^{3}}{4}\right)\mu^{2}.$

We again observe that the right-hand side of the above inequality is uniformly
convex for $\mu>0$ with minimizer
$\mu^{*}_{\tilde{\mathcal{A}}}:=\left(\frac{8\sigma^{2}\tilde{j}}{L_{1}^{2}(\tilde{j}+6)^{3}}\right)^{1/4}$.
Analogously to the proof of STARS Theorem 4.3 (Modified), our optimal choice
of smoothing, given the information we have available, will require us to swap
out $L_{1}$ and $\sigma^{2}$ in $\mu^{*}_{\mathcal{A}}$ with their estimates,
$\hat{L_{1}}$ and $\hat{\sigma}^{2}$, and to swap out $j$ for $\tilde{j}$,
recovering $\hat{\mu}^{*}_{\tilde{\mathcal{A}}}$ (97). This particular choice
$\mu=\hat{\mu}^{*}_{\tilde{\mathcal{A}}}$ can be plugged into (100), which
gives us the bound in (98), our main result. $\blacksquare$

Next, we redefine $\mathcal{Q}_{k}:=\\{\xi_{k}\\}_{k=1}^{M_{\mathcal{A}}}$ and
$\mathcal{U}_{k}:=\\{u^{(k)}\\}_{k=1}^{M_{\mathcal{A}}}$, which are two sets
containing all random variables that appear in FAASTARS’ regular STARS burn-in
phase, iterations $k=1,\ldots,M_{\mathcal{A}}$. Likewise, we extend the
definitions of each set so that
$\mathcal{Q}_{k}=\\{\xi_{k}\\}_{k=M_{\mathcal{A}}+1}^{M}$ and
$\mathcal{U}_{k}^{\tilde{\mathcal{A}}}=\\{u^{(k)}_{\tilde{\mathcal{A}}}\\}_{k=M_{\mathcal{A}}+1}^{M}$
also contain all random variables that appear in FAASTARS’ approximate ASTARS
phase, iterations $k=M_{\mathcal{A}}+1,\ldots,M$.

Let $\phi_{0}:=f(x^{(0)})$,
$\phi_{k}:=\mathbb{E}_{\mathcal{Q}_{k-1},\mathcal{U}_{k-1}}(f(x^{(k)}))$,
$1\leq k\leq M_{\mathcal{A}},$ and
$\phi_{k}:=\mathbb{E}_{\mathcal{Q}_{k-1},\mathcal{U}_{k-1}^{\tilde{\mathcal{A}}}}(f(x^{(k)}))$,
$M_{\mathcal{A}}+1\leq k\leq M$ where the $x^{(k)}$’s are now STARS iterates
for $1\leq k\leq M_{\mathcal{A}},$ and FAASTARS iterates for
$M_{\mathcal{A}}+1\leq k\leq M$.

FAASTARS Corollary 2 (FAASTARS, approximate ASTARS phase result): For all
FAASTARS iterates in phase 3, $k=M_{\mathcal{A}}+1,\ldots,M$, let the vectors
$u_{\tilde{\mathcal{A}}}^{(k)}$ denote those drawn using Algorithm 5. Also,
let $f\in\mathcal{C}^{1,1}(\Lambda)$ with $f$ convex, and let i.i.d. noise
draws $\epsilon(\xi)$ be additive, zero mean, with bounded variance
$\sigma^{2}$ for all appearing $\xi$. We assume we have fixed estimates
$\hat{L_{1}}$ and $\hat{\sigma}$ with $K_{1}>0$ and $K_{2}>0$ as in (22). Let
$\tilde{\mathcal{A}}$ denote the approximated $\tilde{j}$-dimensional AS of
$f$. Let $x^{(M_{\mathcal{A}})}$ be fixed and given and let
$\\{x^{(k)}\\}_{k=M_{\mathcal{A}}+1}^{M}$ denote a sequence of FAASTARS
iterates formed using the approximate active hyperparameters in (97)).
Finally, we require $0<K_{1}<4$ and $K_{2}>0$, the values defined in (22).
Then for any $M-M_{\mathcal{A}}$ total number of approximate ASTARS iterations
within FAASTARS, $k=M_{\mathcal{A}}+1,\ldots,M$,

(101)
$\sum_{k=M_{\mathcal{A}}}^{M}\frac{\phi_{k}-f^{*}_{\tilde{\mathcal{A}}}}{M-M_{\mathcal{A}}}\leq\frac{4L_{1}(\tilde{j}+4)||P_{\tilde{\mathcal{A}}}(x^{(M_{\mathcal{A}})}-x^{*})||^{2}}{\sqrt{K_{1}}(2-\sqrt{K_{1}})(M-M_{\mathcal{A}}+1)}+\frac{4\sigma(\tilde{j}+4)}{\sqrt{2K_{2}}(2-\sqrt{K_{1}})}C_{5},$

where $C_{5}:=\sqrt{K_{1}}\cdot 0.036+\frac{3K_{1}+K_{2}}{16}\cdot 1.034$.

Proof: The proof is almost identical to the proof of ASTARS Corollary 3, but
with $\tilde{j}$’s replacing the roles of $j$’s (since we take steps with
$\tilde{j}$-dimensional $u^{(k)}_{\tilde{\mathcal{A}}}$ vectors and not
$u^{(k)}_{\mathcal{A}}\in\mathcal{A}$), $\hat{\mu}^{*}_{\tilde{\mathcal{A}}}$
replacing $\mu^{*}_{\mathcal{A}}$, and $x^{*}_{\tilde{\mathcal{A}}}$ replacing
$x^{*}_{\mathcal{A}}$ (since we are converging to the minimum of
$f_{\tilde{\mathcal{A}}}$). Here, $K_{1}$ and $K_{2}$ are not necessarily
equal to $1$ as before, since we are assuming inexact active hyperparameters,
formed with the estimates $\hat{L_{1}}$ and $\hat{\sigma}^{2}$. We outline the
required changes one must make to the proofs of STARS Lemma 4.4 (Modified) and
STARS Theorem 4.5 (Modified) to obtain our desired result. These changes
essentially amount to keeping the logic from STARS Lemma 4.4 (Modified) and
STARS Theorem 4.5 (Modified) to account for $K_{1}$ and $K_{2}$ but to
replacing $j$ with $\tilde{j}$.

We begin by obtaining a bound analogous to that of STARS Lemma 4.4 (Modified),
but note that the approximate ASTARS steps are taken with random vectors
$u^{(k)}_{\tilde{\mathcal{A}}}\in\tilde{\mathcal{A}}$, and
$\dim\tilde{\mathcal{A}}=\tilde{j}$. First, we replace $u^{(k)}$ in (37) with
$u^{(k)}_{\tilde{\mathcal{A}}}$, so $g_{0}(x^{(k)}):=\langle\nabla
f(x^{(k)}),u^{(k)}_{\tilde{\mathcal{A}}}\rangle
u^{(k)}_{\tilde{\mathcal{A}}}$. Similarly, note that we set
$u^{(k)}=u^{(k)}_{\tilde{\mathcal{A}}}$ in (20). Also, let
$s_{\hat{\mu}^{*}_{\tilde{\mathcal{A}}}}=s_{\hat{\mu}^{*}_{\tilde{\mathcal{A}}}}^{\tilde{\mathcal{A}}}$
for cleaner notation. Then, using FAASTARS Corollary 1 – instead of ASTARS
Corollary 2 – (38) becomes

(102)
$\mathbb{E}\left(||s_{\hat{\mu}^{*}_{\tilde{\mathcal{A}}}}||^{2}-2\langle
s_{\hat{\mu}^{*}_{\tilde{\mathcal{A}}}},g_{0}(x^{(k)})\rangle+||g_{0}(x^{(k)})||^{2}\right)\leq\frac{K_{1}+K_{2}}{\sqrt{2K_{1}K_{2}}}\sigma
L_{1}\sqrt{\tilde{j}(\tilde{j}+6)^{3}}.$

Continuing with $u^{(k)}_{\tilde{\mathcal{A}}}$ replacing $u^{(k)}$, we
proceed identically, noting
$C_{1}=\frac{K_{1}+K_{2}}{\sqrt{2K_{1}K_{2}}}\sigma
L_{1}\sqrt{\tilde{j}(\tilde{j}+6)^{3}}$. Modifying (19) for
$x^{(k)}\in\tilde{\mathcal{A}}$ with $\dim\tilde{\mathcal{A}}=\tilde{j}$, we
have

(103)
$\mathbb{E}\left(||s_{\hat{\mu}^{*}_{\tilde{\mathcal{A}}}}||^{2}\right)\leq
2(\tilde{j}+4)||\nabla
f(x^{(k)})||^{2}+\frac{(\hat{\mu}^{*}_{\tilde{\mathcal{A}}})^{2}L_{1}^{2}}{2}(\tilde{j}+6)^{3}+C_{1},$

Plugging in the value of $\hat{\mu}^{*}_{\tilde{\mathcal{A}}}$, we obtain the
bound

(104)
$\mathbb{E}\left(||s_{\hat{\mu}^{*}_{\tilde{\mathcal{A}}}}||^{2}\right)\leq
2(\tilde{j}+4)||\nabla f(x^{(k)})||^{2}+C_{2},$

where
$C_{2}:=\frac{3K_{1}+K_{2}}{\sqrt{2K_{1}K_{2}}}L_{1}\sigma\sqrt{\tilde{j}(\tilde{j}+6)^{3}}=\frac{3K_{1}+K_{2}}{\sqrt{2}}\hat{L_{1}}\hat{\sigma}\sqrt{\tilde{j}(\tilde{j}+6)^{3}}$.
The bound in (104) is the analogous result to STARS Lemma 4.4 (Modified) in
the case of ASTARS performed in the estimated $\tilde{\mathcal{A}}$ and with
inexact hyperparameters.

We now proceed to proving the analogous result to STARS Theorem 4.5 (Modified)
in our case. We redefine $r_{k}:=||x^{(k)}-x^{*}_{\tilde{\mathcal{A}}}||$ for
approximate ASTARS iterates $x^{(k)}$. The first three equations appearing in
the proof of STARS Theorem 4.5 (Modified), (48) through (50), are nearly
identical for this proof – one must replace $x^{*}$ with
$x^{*}_{\tilde{\mathcal{A}}}$, replace $\hat{h}$ with
$\hat{h}_{\tilde{\mathcal{A}}}$, and note that here we have
$s_{\mu^{(k)}}=s^{\tilde{\mathcal{A}}}_{\hat{\mu}^{*}_{\tilde{\mathcal{A}}}}$.
Then, invoking (104), we rewrite (51) in our case as

(105) $\mathbb{E}(r_{k+1}^{2})\leq
r_{k}^{2}-2\hat{h}_{\tilde{\mathcal{A}}}\langle\nabla
f_{\mu}(x^{(k)}),x^{(k)}-x^{*}_{\tilde{\mathcal{A}}}\rangle+\hat{h}_{\tilde{\mathcal{A}}}^{2}\left(2(\tilde{j}+4)||\nabla
f(x^{(k)})||^{2}+C_{2}\right),$

where we recall
$C_{2}=\frac{3K_{1}+K_{2}}{\sqrt{2}}\hat{L_{1}}\hat{\sigma}\sqrt{\tilde{j}(\tilde{j}+6)^{3}}$
here. Now, again, equations (52) through (56) are identical for this proof as
long as $x^{*}_{\tilde{\mathcal{A}}}$ replaces $x^{*}$ and
$\hat{h}_{\tilde{\mathcal{A}}}$ replaces $\hat{h}$ throughout, and similarly
for equations (57) through (59) with the additional needed replacement of
$\tilde{j}$ for $P$ and using $C_{2}$ as we have defined it here. Taking
$C_{3}=\hat{h}_{\tilde{\mathcal{A}}}\hat{\mu}^{*}_{\tilde{\mathcal{A}}}L_{1}\tilde{j}+\hat{h}_{\tilde{\mathcal{A}}}^{2}C_{2}$,
(60) holds (with the usual replacments) and plugging
$\hat{h}_{\tilde{\mathcal{A}}}$ in (105):

(106) $\mathbb{E}(r_{k+1}^{2})\leq
r_{k}^{2}-\frac{1}{4L_{1}(\tilde{j}+4)}(f(x^{(k)})-f(x^{*}_{\tilde{\mathcal{A}}}))+C_{3}.$

Now plugging $\hat{h}_{\tilde{\mathcal{A}}}$ into $C_{3}$, we obtain
$C_{3}\leq\frac{\sqrt{K_{1}}}{\sqrt{K_{2}}}\cdot\frac{\sigma}{\sqrt{2}L_{1}}\left(\sqrt{K_{1}}\cdot
0.036+\frac{3K_{1}+K_{2}}{16}\cdot 1.034\right)=:C_{4}$, so

(107) $\mathbb{E}(r_{k+1}^{2})\leq
r_{k}^{2}-\frac{1}{4L_{1}(\tilde{j}+4)}(f(x^{(k)})-f(x^{*}_{\tilde{\mathcal{A}}}))+C_{4}.$

We may now apply the expectation over $\mathcal{P}_{k}$ and
$\mathcal{U}_{k}^{\tilde{\mathcal{A}}}$, $k=M_{\mathcal{A}},\ldots,M$,
rearrange, and sum over $k=M_{\mathcal{A}},\ldots,M$, similarly to before. We
obtain a modification to (68) with

(108)
$\sum_{k=M_{\mathcal{A}}}^{M}\frac{\phi_{k}-f^{*}_{\tilde{\mathcal{A}}}}{M-M_{\mathcal{A}}+1}\leq\frac{4L_{1}(\tilde{j}+4)}{(M-M_{\mathcal{A}}+1)}\left(r_{M_{\mathcal{A}}}-\mathbb{E}_{\mathcal{U}_{k},\mathcal{P}_{k}}(r_{k+1}^{2})\right)+4L_{1}(\tilde{j}+4)C_{4}.$

We drop the strictly negative term (again involving
$\mathbb{E}_{\mathcal{P}_{k},\mathcal{U}_{k}^{\tilde{\mathcal{A}}}}$) and plug
in the definition of
$r_{M_{\mathcal{A}}}=||x^{(M_{\mathcal{A}})}-x^{*}_{\tilde{\mathcal{A}}}||^{2}$.
Then, recalling
$x^{*}_{\tilde{\mathcal{A}}}=P_{\tilde{\mathcal{A}}}(x^{*})+P_{\tilde{\mathcal{I}}}(x^{(M_{\mathcal{A}})})$,
writing
$x^{(M_{\mathcal{A}})}=P_{\tilde{\mathcal{A}}}(x^{(M_{\mathcal{A}})})+P_{\tilde{\mathcal{I}}}(x^{(M_{\mathcal{A}})})$,
and noting that $P_{\tilde{\mathcal{A}}}$ is linear, we have
$r_{M_{\mathcal{A}}}=||P_{\tilde{\mathcal{A}}}(x^{(M_{\mathcal{A}})}-x^{*})||^{2}$.
Plugging in the value of $C_{4}$, we obtain (101). $\blacksquare$

FAASTARS Corollary 2 shows that its iterates during its third phase
(approximate ASTARS) converge to $f^{*}_{\tilde{\mathcal{A}}}$; however, we
need a result for the convergence to $f^{*}$ like we did in the last section.
In particular, we need a result analogous to ASTARS Corollary 4, generally
regarding the distance between evaluations of $f$ and
$f_{\tilde{\mathcal{A}}}$ (instead of $f_{\mathcal{A}}$) at their respective
minimizers, $x^{*}$ and $x^{*}_{\tilde{\mathcal{A}}}$. Recall, by
$f_{\tilde{\mathcal{A}}}$, we mean
$f_{\mathcal{A}}(\lambda):=f(V_{\tilde{\mathcal{A}}}V_{\tilde{\mathcal{A}}}^{\top}\lambda)$,
where the application of
$V_{\tilde{\mathcal{A}}}V_{\tilde{\mathcal{A}}}^{\top}$ is a linear
transformation of $\lambda$ into $\tilde{\mathcal{A}}$, as in [ConstantineK],
which we again rely on heavily for the following result.

FAASTARS Corollary 3: Let $x^{(M_{\mathcal{A}})}$ denote any initial iterate
for phase 3 of FAASTARS. Let $\mathcal{A}$ continue to denote the true
$j$-dimensional AS of $f$ and let $\tilde{\mathcal{A}}$ denote the approximate
$\tilde{j}$-dimensional AS of $f$ and the true minimizer of
$f_{\tilde{\mathcal{A}}}$ with $x^{*}_{\tilde{\mathcal{A}}}$ and corresponding
noiseless function evaluation $f^{*}_{\tilde{\mathcal{A}}}$. Assume
$||x^{*}-x^{*}_{\tilde{\mathcal{A}}}||_{\Lambda}<\infty$, where
$||\cdot||_{\Lambda}$ denotes a norm. Assume $\tilde{j}=j$ and for $\delta>0$,
$||V-\tilde{V}||_{2}<\delta$, where $||\cdot||_{2}$ here denotes the matrix
2-norm induced by the 2-norm in $\Lambda$. Also, assume the sign of
$(\tilde{V}_{\tilde{\mathcal{I}}})_{i}$, the $i$-th column of
$\tilde{V}_{\tilde{\mathcal{I}}}$ to be chosen so that
$||(\tilde{V}_{\tilde{\mathcal{I}}})_{i}-(V_{\mathcal{I}})_{i}||_{2}$ is
minimized for $i=j+1,\ldots,P$, where $||\cdot||_{2}$ denotes the vector
2-norm. Then the difference between $|f^{*}-f^{*}_{\tilde{\mathcal{A}}}|$ is
bounded by

(109)
$|f^{*}-f^{*}_{\tilde{\mathcal{A}}}|\leq\sqrt{\tilde{a}_{1}}\left(\delta\sqrt{q_{1}+\cdots+q_{j}}+\sqrt{q_{j+1}+\cdots+q_{P}}\right),$

where $0\leq\tilde{a}_{1}<\infty$ is an eigenvalue of the positive-semi
definite matrix
$(x^{*}-x^{*}_{\tilde{\mathcal{A}}})(x^{*}-x^{*}_{\tilde{\mathcal{A}}})^{\top}$
and the $q$’s are our notations for the exact eigenvalues associated with the
eigendecomposition of the sensitivity matrix $W$. Also,
$\tilde{a}_{1}=||P_{\tilde{\mathcal{I}}}(x^{(M_{\mathcal{A}})}-x^{*})||_{\Lambda}^{2}.$

Proof: We bound the quantity $(f^{*}-f^{*}_{\tilde{\mathcal{A}}})^{2}$. First,
we expand $f*=f(x^{*})$ around $x^{*}_{\tilde{\mathcal{A}}}$ by applying
Extended Mean Value Theorem analogously to ASTARS Corollary 4. For $c\in[0,1]$
and $\tilde{z}:=cx_{\tilde{\mathcal{A}}}^{*}+(1-c)x^{*}$ we have
$f^{*}=f(x^{*})=f(x_{\tilde{\mathcal{A}}}^{*})+\nabla
f(\tilde{z})^{\top}(x^{*}-x^{*}_{\tilde{\mathcal{A}}})$. Note that since the
components of $x_{\tilde{\mathcal{A}}}^{*}$ and rotated
$V_{\tilde{\mathcal{A}}}V_{\tilde{\mathcal{A}}}^{\top}x^{*}$ must match for
indices $i=1,\ldots,j$, the point $\tilde{z}$ varies along
$\tilde{\mathcal{I}}$ only; that is, $\tilde{z}$ is fixed in
$\tilde{\mathcal{A}}$. Thus, $\nabla
f(\tilde{z})=\nabla_{\tilde{\mathcal{I}}}f(\tilde{z})$, the gradient taken in
the approximate inactive subspace only.

The expansion of $f^{*}$ around $x^{*}_{\tilde{\mathcal{A}}}$ writes
$(f^{*}-f^{*}_{\tilde{\mathcal{A}}})^{2}=\left(\nabla_{\tilde{\mathcal{I}}}f(\tilde{z})^{\top}(x^{*}-x^{*}_{\tilde{\mathcal{A}}})\right)^{2}=\nabla_{\tilde{\mathcal{I}}}f(\tilde{z})^{\top}\tilde{A}\nabla_{\tilde{\mathcal{I}}}f(\tilde{z})$,
where
$\tilde{A}:=(x^{*}-x^{*}_{\tilde{\mathcal{A}}})(x^{*}-x^{*}_{\tilde{\mathcal{A}}})^{\top}$
is a $P\times P$ matrix. Note that $\tilde{A}$ is a square, positive semi-
definite, rank 1 matrix. Hence, it has 1 eigenvalue that is positive or zero,
which we denote with $\tilde{a}_{1}\geq 0$; all other $P-1$ eigenvalues are
$0$. We find
$\tilde{a}_{1}=||x^{*}-x^{*}_{\tilde{\mathcal{A}}}||^{2}_{\Lambda}$ so by
definition,
$\tilde{a}_{1}=||P_{\tilde{\mathcal{I}}}(x^{(M_{\mathcal{A}})}-x^{*})||_{\Lambda}^{2}.$
Observe $\tilde{a}_{1}<\infty$ because
$||x^{*}-x^{*}_{\tilde{\mathcal{A}}}||_{\Lambda}<\infty$ and also
$\tilde{a}_{1}=0\iff x^{*}=x^{*}_{\tilde{\mathcal{A}}}$, which case
$f^{*}=f^{*}_{\tilde{\mathcal{A}}}$ so that
$|f^{*}-f^{*}_{\tilde{\mathcal{A}}}|=0$. Then
$(f^{*}-f^{*}_{\tilde{\mathcal{A}}})^{2}\leq\tilde{a}_{1}\nabla_{\tilde{\mathcal{I}}}f(\tilde{z})^{\top}\nabla_{\tilde{\mathcal{I}}}f(\tilde{z})$.

Recall $\tilde{j}=j$ and for $\delta>0$, $||V-\tilde{V}||_{2}<\delta$. Then we
have comparable partitions, in the sense that the submatrices
$V_{\mathcal{A}}$ and $\tilde{V}_{\tilde{\mathcal{A}}}$ (of $V$ and
$\tilde{V}$ respectively) are both $j$-dimensional, and likewise the
submatrices $V_{\mathcal{I}}$ and $V_{\tilde{\mathcal{I}}}$ are both
$P-j$-dimensional. We shall need the results from Lemma 3.4 in [ConstantineK]
which state $||V_{\mathcal{I}}^{\top}\tilde{V}_{\tilde{\mathcal{I}}}||_{2}\leq
1$ and
$||V_{\mathcal{A}}^{\top}\tilde{V}_{\tilde{\mathcal{I}}}||_{2}\leq\delta$,
where $||\cdot||_{2}$ denotes the matrix 2-norm. Note that Lemma 3.4 requires
that the sign of $(\tilde{V}_{\tilde{\mathcal{I}}})_{i}$, the $i$-th column of
$\tilde{V}_{\tilde{\mathcal{I}}}$ to be chosen so that
$||(\tilde{V}_{\tilde{\mathcal{I}}})_{i}-(V_{\mathcal{I}})_{i}||_{2}$ is
minimized for $i=j+1,\ldots,P$, where $||\cdot||_{2}$ denotes the vector
2-norm. Now the chain rule provides
$\nabla_{\tilde{\mathcal{I}}}f=V_{\mathcal{I}}^{\top}\tilde{V}_{\tilde{\mathcal{I}}}\nabla_{\mathcal{I}}f+V_{\mathcal{A}}^{\top}\tilde{V}_{\tilde{\mathcal{I}}}\nabla_{\mathcal{A}}f$
([ConstantineK], pp. A1510).

Applying the expectation over both $\mathcal{I}$ and $\mathcal{A}$ to both
sides (where the left-hand side is constant with respect to this expectation)
we have
$(f^{*}-f^{*}_{\tilde{\mathcal{A}}})\leq\tilde{a}_{1}\mathbb{E}_{\mathcal{A},\mathcal{I}}\left(\nabla_{\tilde{\mathcal{I}}}f(\tilde{z})^{\top}\nabla_{\tilde{\mathcal{I}}}f(\tilde{z})\right)$.
Using the chain rule and
$||V_{\mathcal{I}}^{\top}\tilde{V}_{\tilde{\mathcal{I}}}||_{2}\leq 1$ and
$||V_{\mathcal{A}}^{\top}\tilde{V}_{\tilde{\mathcal{I}}}||_{2}\leq\delta$
gives

(110)
$(f^{*}-f^{*}_{\tilde{\mathcal{A}}})^{2}\leq\tilde{a}_{1}\left(\mathbb{E}_{\mathcal{I}}\left(\nabla_{\mathcal{I}}f(\tilde{z})^{\top}\nabla_{\mathcal{I}}f(\tilde{z})\right)+2\delta\,\mathbb{E}_{\mathcal{A},\mathcal{I}}\left(\nabla_{\mathcal{I}}f(\tilde{z})^{\top}\nabla_{\mathcal{A}}f(\tilde{z})\right)+\delta^{2}\mathbb{E}_{\mathcal{A}}\left(\nabla_{\mathcal{A}}f(\tilde{z})^{\top}\nabla_{\mathcal{A}}f(\tilde{z})\right)\right),$

where linearity of $\mathbb{E}$ allows for breaking the expectation of the
appearing sum into a sum of expectations taken over only $\mathcal{A}$, only
$\mathcal{I}$, or both $\mathcal{A}$ and $\mathcal{I}$, depending on which
types of gradients appear. The Cauchy-Schwarz inequality may be applied to the
second term on the right-hand side of (110) to write

(111)
$2\delta\mathbb{E}_{\mathcal{A},\mathcal{I}}\left(\nabla_{\mathcal{I}}f(\tilde{z})^{\top}\nabla_{\mathcal{A}}f(\tilde{z})\right)\leq
2\delta\,\mathbb{E}_{\mathcal{I}}\left(\nabla_{\mathcal{I}}f(\tilde{z})^{\top}\nabla_{\mathcal{I}}f(\tilde{z})\right)\mathbb{E}_{\mathcal{A}}\left(\nabla_{\mathcal{A}}f(\tilde{z})^{\top}\nabla_{\mathcal{I}}f(\tilde{z}))\right),$

where we can split up the two terms inside
$\mathbb{E}_{\mathcal{A},\mathcal{I}}(\cdot)$ into two separate expectations
taken over $\mathcal{I}$ and $\mathcal{A}$ individually, due to their
independence. Substituting this bound into (86) and then factoring the
resulting terms will yield

(112)
$(f^{*}-f^{*}_{\tilde{\mathcal{A}}})^{2}\leq\tilde{a}_{1}\left(\mathbb{E}_{\mathcal{I}}\left(\nabla_{\mathcal{I}}f(\tilde{z})^{\top}\nabla_{\mathcal{I}}f(\tilde{z})\right)^{1/2}+\delta\,\mathbb{E}_{\mathcal{A}}\left(\nabla_{\mathcal{A}}f(\tilde{z})^{\top}\nabla_{\mathcal{A}}f(\tilde{z})\right)^{1/2}\right)^{2},$

Citing [ConstantineK] Lemma 2.2, we have
$\mathbb{E}_{\mathcal{A}}\left(\nabla_{\mathcal{A}}f(\tilde{z})^{\top}\nabla_{\mathcal{A}}f(\tilde{z})\right)<q_{1}+\cdots+q_{j}$,
where $q_{i}$, $i=1,\ldots,j$ are the first $j$ eigenvalues of the sensitivity
matrix $W$ and
$\mathbb{E}_{\mathcal{I}}\left(\nabla_{\mathcal{I}}f(\tilde{z})^{\top}\nabla_{\mathcal{I}}f(\tilde{z})\right)<q_{j+1}+\cdots+q_{P}$,
where $q_{i}$, $i=j+1,\ldots,P$ are the last $P-j$ eigenvalues of the
sensitivity matrix $W$. Hence,
$(f^{*}-f^{*}_{\mathcal{A}})^{2}\leq\tilde{a}_{1}\left(\delta\,\sqrt{q_{1}+\cdots+q_{j}}+\sqrt{q_{j+1}+\cdots+q_{P}}\right)^{2}$.
Applying a square root to both sides, we obtain (109). $\blacksquare$

We now present a statement about the convergence of FAASTARS to $f^{*}$ by
combining the previous two results, along with STARS Theorem 4.5 (Modified).

Recall for $k=1,\ldots,M_{\mathcal{A}}$ iterates correspond to STARS with
estimated hyperparameters and for $k=M_{\mathcal{A}}+1,\ldots,M$ the iterates
correspond to ASTARS with estimated hyperparameters and estimated
$\mathcal{A}$.

FAASTARS Theorem 4: For $k=M_{\mathcal{A}},\ldots,M$, let the assumptions from
FAASTARS Corollaries 3 and 4 hold. For any total number of FAASTARS iterations
$M\geq 1$,

(113)
$\begin{split}\sum_{k=M_{\mathcal{A}}}^{M}\frac{\phi_{k}-f^{*}}{M-M_{\mathcal{A}}+1}\leq&\frac{4L_{1}(j+4)||P_{\tilde{\mathcal{A}}}(x^{(M_{\mathcal{A}})}-x^{*})||^{2}}{\sqrt{K_{1}}(2-\sqrt{K_{1}})(M-M_{\mathcal{A}}+1)}+\frac{4\sigma(j+4)}{\sqrt{2K_{2}}(2-\sqrt{K_{1}})}C_{5}\\\
&+\sqrt{\tilde{a}_{1}}\left(\delta\sqrt{q_{1}+\cdots+q_{j}}+\sqrt{q_{j+1}+\cdots+q_{P}}\right),\end{split}$

where $C_{5}:=\sqrt{K_{1}}\cdot 0.036+\frac{3K_{1}+K_{2}}{16}\cdot 1.034$.

Proof: Using the triangle inequality, we have
$|\phi_{k}-f^{*}|\leq|\phi_{k}-f^{*}_{\tilde{\mathcal{A}}}|+|f^{*}-f^{*}_{\tilde{\mathcal{A}}}|.$
Thus, noting that the quantity $\phi_{k}-f^{*}_{\tilde{\mathcal{A}}}$ is
always nonnegative, we use the triangle inequality to rewrite the summation on
the left-hand side of (113) as

(114)
$\sum_{k=M_{\mathcal{A}}}^{M}\frac{\phi_{k}-f^{*}}{M-M_{\mathcal{A}}+1}\leq\sum_{k=M_{\mathcal{A}}}^{M}\frac{\phi_{k}-f^{*}_{\tilde{\mathcal{A}}}}{M-M_{\mathcal{A}}+1}+|f^{*}-f^{*}_{\tilde{\mathcal{A}}}|.$

Now the second sum on the right-hand side of (115) is in a more useful form
for us, since the final phase of iterates, $k=M_{\mathcal{A}},\ldots,M$ are
(approximate) ASTARS iterates converging to $f^{*}_{\tilde{\mathcal{A}}}$ (not
$f^{*}$); thus, we recall that we already proved a bound for this sum in
FAASTARS Corollary 2. Note that since we assume in FAASTARS Corollary 3 that
$\tilde{j}=j$, we replace $\tilde{j}$ with $j$ in the result from FAASTARS
Corollary 2. Otherwise, we just invoke FAASTARS Corollary 3, bounding the last
term on the right-hand side of (114) to obtain (113). $\blacksquare$

We use the results above to analyze the complexity of FAASTARS in the
following remark.

Remark: We again mimic the complexity analyses we performed in the preceding
sections for FAASTARS. Let $||\cdot||$ denote a norm in $\Lambda$ throughout.
Define $R_{\tilde{\mathcal{A}}}^{2}$ as a bound
$||P_{\tilde{\mathcal{A}}}(x^{(M_{\mathcal{A}})}-x^{*})||^{2}\leq
R_{\tilde{\mathcal{A}}}^{2}$. We recall
$\tilde{a}_{1}=||P_{\tilde{\mathcal{I}}}(x^{(M_{\mathcal{A}})}-x^{*})||_{\Lambda}^{2}$
and define a bound $\tilde{a}_{1}\leq R_{\tilde{\mathcal{I}}}^{2}$.

Define
$x^{\dagger}:=\text{argmin}_{x}\\{f(x):x\in\\{x^{(0)},\ldots,x^{M}\\}\\}$ and
$\phi_{\dagger}:=\mathbb{E}_{\mathcal{U}_{k-1},\mathcal{P}_{k-1}}(f(x^{\dagger}))$
as before. Then the value $\phi_{\dagger}-f^{*}$ must be less than or equal to
the average improvement for any given run of STARS; that is, (113), along with
our new definitions, implies that

(115)
$\begin{split}\phi_{\dagger}-f^{*}\leq&\frac{4L_{1}(j+4)}{\sqrt{K_{1}}(2-\sqrt{K_{1}})(M-M_{\mathcal{A}}+1)}R_{\tilde{\mathcal{A}}}^{2}+\frac{4\sigma(j+4)}{\sqrt{2K_{2}}(2-\sqrt{K_{1}})}C_{5}\\\
&+R_{\tilde{\mathcal{I}}}\left(\delta\sqrt{q_{1}+\ldots+q_{j}}+\sqrt{q_{j+1}+\cdots+q_{P}}\right),\end{split}$

Now in this analysis, we are only taking enough approximate STARS steps to
learn a surrogate for $f$ to obtain $\tilde{\mathcal{A}}$ Thus, we have
$M_{\mathcal{A}}\sim\mathcal{O}(L(P))$, where $L$ is a function of $P$ that
depends on the surrogate method used to learn $\mathcal{A}$. For instance, if
we use a linear surrogate or RBF’s – which begin by fitting a linear
surrogate, and then later a quadratic surrogate once enough ASTARS steps are
taken – then $L(P)=P$. If we use quadratic surrogates, then $L(P)=P^{2}$.
(Typically we use quadratic surrogates for higher-quality active subspaces.)

The terms involving Phase 3 of FAASTARS (approximate ASTARS phase) are the
three terms on the right-hand side of (115). We assume that we wish to achieve
a final accuracy of $\epsilon_{\text{tol}}>0$. Then we will need
$\phi_{\dagger}-f^{*}\leq\epsilon_{\text{tol}}$. If we take

(116)
$\frac{4\sigma(j+4)}{\sqrt{2K_{2}}(2-\sqrt{K_{1}})}C_{5}\leq\frac{\epsilon_{\text{tol}}}{3},$

then we must require that the noise not exceed the following threshold:

(117)
$\sigma\leq\frac{\sqrt{2K_{2}}(2-\sqrt{K_{1}})\epsilon_{\text{tol}}}{12\sigma(j+4)}.$

If we satisfy (117) and also have

(118)
$R_{\tilde{\mathcal{I}}}\left(\delta\sqrt{q_{1}+\ldots+q_{j}}+\sqrt{q_{j+1}+\cdots+q_{P}}\right)\leq\frac{\epsilon_{\text{tol}}}{3},$

then we can achieve $\epsilon_{\text{tol}}$ accuracy as long as

(119)
$\frac{4L_{1}(j+4)}{\sqrt{K_{1}}(2-\sqrt{K_{1}})(M-M_{\mathcal{A}}+1)}R_{\tilde{\mathcal{A}}}^{2}\leq\frac{\epsilon_{\text{tol}}}{3}\iff
M\geq
M_{\mathcal{A}}+\frac{12L_{1}(j+4)R_{\tilde{\mathcal{A}}}^{2}}{\epsilon_{\text{tol}}\sqrt{K_{1}}(2-\sqrt{K_{1}})}-1.$

Hence, we achieve $\epsilon_{\text{tol}}$ accuracy as long as: the noise is
small enough; the eigenvalues of the inactive subspace and distance from
$x^{*}$ to $x^{*}_{\tilde{\mathcal{A}}}$ is small enough; and $M$ is large
enough. Details of those required bounds given by (117), (118), and (119),
respectively. With these assumptions, we achieve $\epsilon_{\text{tol}}$
accuracy in

(120)
$M\sim\mathcal{O}\left(\max\left\\{L(P),\,\frac{L_{1}jR_{\tilde{\mathcal{A}}}^{2}}{\epsilon_{\text{tol}}\sqrt{K_{1}}(2-\sqrt{K_{1}})}\right\\}\right).$

This analysis also shows that given a particular variance in the noise,
$\sigma$, as well as the term involving the eigenvalues of the inactive
subspace and distance from $x^{*}$ to $x^{*}_{\tilde{\mathcal{A}}}$, the
achievable accuracy can be no better (i.e., less) than the value given below:

(121)
$\epsilon_{\text{tol}}\geq\max\left\\{\frac{12(j+4)C_{5}}{\sqrt{K_{1}}(2-\sqrt{K_{1}})}\sigma,\,3R_{\tilde{\mathcal{I}}}\left(\delta\sqrt{q_{1}+\cdots+q_{j}}+\sqrt{q_{j+1}+\cdots+q_{P}}\,\right)\right\\},$

showing that the achievable accuracy will either be limited (mainly) by
hyperparameter approximations or (mainly) by the error in
$\tilde{\mathcal{A}}$.

We find similar complexity results to that of ASTARS, but also pay a price for
approximations to $\tilde{\mathcal{A}}$, especially when the approximation is
poor, usually do to both an insufficient number of samples and samples that do
not explore $\Lambda$ sufficiently. (Recall that since $\tilde{\mathcal{A}}$
is formed using information from a surrogate trained on STARS samples in the
FAASTARS routine, it is really poor surrogates that hurt convergence.)