# Statistical Analysis of Complex Computer Models in Astronomy

Joshua Lukemire  Department of Biostatistics and Bioinformatics, Emory
University, Atlanta, GA, USA  Qian Xiao  Department of Statistics, University
of Georgia, Athens, GA, USA  Abhyuday Mandal<EMAIL_ADDRESS>Department
of Statistics, University of Georgia, Athens, GA, USA  Weng Kee Wong
Department of Biostatistics, University of California, Los Angeles, CA, USA

###### Abstract

We introduce statistical techniques required to handle complex computer models
with potential applications to astronomy. Computer experiments play a critical
role in almost all fields of scientific research and engineering. These
computer experiments, or simulators, are often computationally expensive,
leading to the use of emulators for rapidly approximating the outcome of the
experiment. Gaussian process models, also known as Kriging, are the most
common choice of emulator. While emulators offer significant improvements in
computation over computer simulators, they require a selection of inputs along
with the corresponding outputs of the computer experiment to function well.
Thus, it is important to select inputs judiciously for the full computer
simulation to construct an accurate emulator. Space-filling designs are
efficient when the general response surface of the outcome is unknown, and
thus they are a popular choice when selecting simulator inputs for building an
emulator. In this tutorial we discuss how to construct these space filling
designs, perform the subsequent fitting of the Gaussian process surrogates,
and briefly indicate their potential applications to astronomy research.

##  1 Introduction

Computer experiments, or simulators, are an increasingly important tool in
many scientific fields. In these experiments, a computer model is defined
relating a set of inputs to an output. Instead of conducting a traditional
experiment, a researcher will provide a set of inputs to the computer model
and obtain the model output. This approach is very appealing in fields such as
physics, where the computer experiment model can be setup using a series of
known relationships/equations and different inputs may consist of unknown
constants in those equations or other properties such as mass or chemical
compositions. These experiments can be effective alternatives to experiments
which may be too expensive or otherwise impossible to perform in a traditional
setting. They differ from standard experiments in several key ways. Most
importantly, computer experiments are generally deterministic; for a set of
input settings the experiment will return the same result every time it is
conducted. Second, the experiments will generally not have an easily described
response surface; for example a standard linear regression model will not
generally describe the outcome accurately.

Figure 1: An example from a simulation examining whether two black holes
merge. Source: https://www.black-holes.org/code/SpEC.html.

Many research areas in astronomy do not easily permit conducting traditional
experiments. For example, researchers may be interested in the formation of
binary black holes. Clearly the researchers will not be able to create
multiple black holes and observe their dynamics over time. Computer
experiments make it possible to study such phenomena by creating computer
models based on the theorized properties of these binary systems and then
comparing the output to what is observed in Nature. For example, Compact
Object Mergers: Population Astrophysics and Statistics (COMPAS) is used to
investigate binary population synthesis. The computer experiment takes input
as initial conditions and simulates the lifespan of stars [50] [56].
Similarly, binary population synthesis code ComBinE has been used to perform
binary population syntheses [32], and the tool UniverseMachine [4] allows
researchers to study galaxy formation.

Computer experiments for many complex systems can be very expensive to perform
(see, for example, [62]). This computational expense can be a significant
problem, especially if a researcher hopes to conduct the experiment for many
sets of inputs. An alternative to directly performing these computer
experiments is to instead create a surrogate or emulator [18]. Surrogates are
popular for computer experiments when it is not realistic to evaluate a fine
grid over the entire input space. Instead, a (relatively) small number of
points are chosen to evaluate under the original computer simulation. Then, a
model is fit to the output from these limited runs. Predictions under this
model for new inputs, as well as uncertainty quantification, can be obtained
from the surrogate without the need to re-run the expensive computer
simulation at the new points. If the model fits well, then the predicted value
will be close to the true value that would have been obtained if the full
computer experiment was used.

The most common tool used to fit the data points and create the surrogate
model is the Gaussian process (GP) [48] [18]. The GP is appealing for creating
surrogates because it interpolates known data to evaluate new data points.
This is especially important when the outcome for a fixed set of inputs is
deterministic, which is frequently the case in computer experiments. This
approach is becoming more popular in the astronomy literature. Some recent
work includes [22], who proposed using Gaussian process emulation to obtain
confidence intervals for the parameter vector of a phase-space distribution
function for dwarf spheroidal galaxies.

Section 2 of this tutorial paper introduces Gaussian process models and
discusses their applications to computer experiments. We provide codes and
examples throughout in the R programming language [45]. Section 3 of this
paper focuses on determining what inputs to use to generate the responses used
to fit the Gaussian process model to obtain an accurate surrogate. We draw
upon the design of experiments statistical literature to discuss design of
computer experiments. In particular we focus on Latin hypercube designs and
discuss several techniques for finding them.

##  2 Surrogates for Computer Models

Simpler surrogates or emulators are often preferred for complex deterministic
computer models. Gaussian Process (GP) models are popular choice for this
purpose [48]. Consider an $n$-run computer experiment with $d$-dimensional
input vectors $\textbf{x}_{i}=(x_{i1},\ldots,x_{id})^{\mbox{\tiny T}}$ and a
deterministic output $y(\textbf{x}_{i})$, for $i=1,2,\ldots,n$. To fix ideas,
assume that we are interested in a 2-dimensional input for a computer
experiment with output given by the Branin function as defined by [6], see
also [14].

$y(x_{1},x_{2})=\left(x_{2}-\frac{5.1}{4\pi^{2}}x_{1}^{2}+\frac{5}{\pi}x_{1}-6\right)^{2}+10\left(1-\frac{1}{8\pi}\right)\cos\left(x_{1}\right)+10,$
(1)

where the design space is given by values of $x_{1}\in[-5,10]$ and
$x_{2}\in[0,15]$. The R code below can be used to evaluate this function.

⬇

1# x is a vector of inputs (length 2)

2branin <- function(x){

3 a <- 1

4 b <- 5.1 / (4 * pi^2)

5 c <- 5 / pi

6 r <- 6

7 s <- 10

8 t <- 1 / (8 * pi)

9 return( a*(x[2] - b*x[1]^2 + c*x[1] -r)^2 + s*(1-t)*cos(x[1]) + s )

10}

The left panel in Figure 2 displays the output for this function over the
entire design space.

Figure 2: Left: the true response under the Branin function. Right: The
estimated response using the surrogate model.

###  2.1 Stationary Gaussian Process - Krigging

The simplest possible GP model, known as ordinary GP or $krigging$, is given
by

$y(\textbf{x}_{i})=\mu+Z(\textbf{x}_{i}),$ (2)

where $\mu$ is the mean and $Z(\textbf{x})$ is a GP, denoted by
$Z(\textbf{x})\sim GP(0,\sigma^{2}R)$. This notation implies that the GP has
zero-mean, and the covariance function
$Cov\left(Z(\textbf{x}_{i}),Z(\textbf{x}_{j})\right)=\sigma^{2}R(\cdot|\mbox{\boldmath${\theta}$})$,
where $\mbox{\boldmath${\theta}$}=(\theta_{1},\ldots,\theta_{d})^{\mbox{\tiny
T}}$ is the vector of unknown correlation parameters with all $\theta_{s}>0$
($s=1,\ldots,d$). The correlation between outputs is determined by a
stationary correlation function $R$ with parameter ${\theta}$. Two of the more
commonly-used correlation functions are the power-exponential and the Gaussian
functions. Under a power-exponential correlation structure the
$(i,j)^{\mbox{th}}$ term is defined as:

$\displaystyle
R(\textbf{x}_{i},\textbf{x}_{j}|\mbox{\boldmath${\theta}$})=\prod_{s=1}^{d}\exp\Bigg{\\{}-\theta_{s}\mid
x_{is}-x_{js}\mid^{p_{s}}\Bigg{\\}}\hskip 14.45377pt\mbox{for all }i,j,$ (3)

where smoothness parameters $p_{1},\ldots,p_{s}$ are all between $0$ and $2$.
Of special importance is $p_{s}=2$, for all $s=1,\ldots,d$, which corresponds
to the popular Gaussian correlation function:

$R(\textbf{x}_{i},\textbf{x}_{j}|\mbox{\boldmath${\theta}$})=\hbox{exp}\left\\{-\sum_{s=1}^{d}\theta_{s}(x_{is}-x_{js})^{2}\right\\}\hskip
14.45377pt\mbox{for all }i,j.$ (4)

The flexibility of the correlation structure is what makes the GP model a
popular surrogate for complex computer models. For any given input
$\textbf{x}^{*}$ in the design space, the fitted GP surrogate gives the
predicted computer model response as,

$\hat{y}(\textbf{x}^{*})=\mu+\mathbf{r}^{\mbox{\tiny
T}}(\textbf{x}^{*})\mathbf{R}^{-1}(\mathbf{y}-\mu\mathbf{1}_{n}),$ (5)

where

$\mathbf{r}(\textbf{x}^{*})=\Bigg{[}\text{corr}\Big{(}Z(\textbf{x}^{*}),Z(\textbf{x}_{1})\Big{)},\text{corr}\Big{(}Z(\textbf{x}^{*}),Z(\textbf{x}_{2})\Big{)},\ldots,\text{corr}\Big{(}Z(\textbf{x}^{*}),Z(\textbf{x}_{n})\Big{)}\Bigg{]}^{\mbox{\tiny
T}},$ (6)

$\mathbf{1}_{n}$ is a vector of ones of length $n$, $\mathbf{R}$ is the
$n\times n$ correlation matrix for
$(Z(\textbf{x}_{1}),...,Z(\textbf{x}_{n}))^{\mbox{\tiny T}}$, $\mathbf{y}$ is
the response vector $(y(\textbf{x}_{1}),\ldots,y(\textbf{x}_{n}))^{\mbox{\tiny
T}}$, and the associated uncertainty estimate is

$s^{2}(\textbf{x}^{*})=\sigma^{2}\Big{(}1-\mathbf{r}(\textbf{x}^{*})^{\mbox{\tiny
T}}\mathbf{R}^{-1}\mathbf{r}(\textbf{x}^{*})\Big{)}.$ (7)

In practice, the parameters $\mu,\sigma^{2}$ and ${\theta}$ in Equations (5)
and (7) are unknown and need to be estimated from the data. The parameters can
be estimated using the mlegp function in R. Assume that we already have a
design with 10 points (details for obtaining this design will be presented in
Section 3). Then we can fit a GP with a Gaussian correlation function as,

⬇

1library(mlegp)

2# Obtaining this design is discussed in Section 3

3design <- matrix(c(-4.25, -1.25, 4.75, 7.75, -2.75, 1.75,

4 6.25, 3.25, 0.25, 9.25, 8.25, 6.75, 11.25, 12.75, 0.75,

5 2.25, 9.75, 5.25, 14.25, 3.75), ncol=2)

6# Obtain the output at the set points in our design

7Yx <- apply(design, 1, branin)

8# Use the observed outputs to construct a surrogate

9branin_surrogate_1 <- mlegp(design, Yx)

Similarly, estimates across the entire design space can be obtained by using
the surrogate model by specifying the inputs on a grid:

⬇

1# Construct a grid of points to obtain predictions at

2x1 <- seq(from = -5, to =10, length.out = 25)

3x2 <- seq(from = 0, to =15, length.out = 25)

4test_points <- expand.grid(x1, x2)

5# Get predictions from the gaussian process

6yhat <- predict(branin_surrogate_1, test_points)

7predictions <- matrix( yhat, nrow = length(x1) )

8# plot the predictions, theta and phi control viewing angle

9persp(x1, x2, predictions, theta = -45, phi=45)

The right panel in Figure 2 displays a plot of the surrogate output. Comparing
this output to the true values in the left panel, it is clear that the
surrogate model is able to obtain a very close approximation to the true
process.

The formulation in equation (2) can be extended to incorporate a global trend
function for the mean $\mu$ [57]. This is known as Universal Kriging:

$y(\textbf{x})=\mu(\textbf{x})+Z(\textbf{x}),$ (8)

with $\mu(\textbf{x})=g(\textbf{x})^{\mbox{\tiny
T}}\beta=\sum_{i=1}^{m}\beta_{i}g_{i}(\textbf{x})$, where $g$ is a
$m$-dimensional known function and
$\beta=(\beta_{1},\ldots,\beta_{m})^{\mbox{\tiny T}}$ is a vector of unknown
parameters. If we assume $g_{1}(\textbf{x})=1$ and let
$\mathbf{G}=(g_{1}(\textbf{x}),\ldots,g_{m}(\textbf{x}))^{\mbox{\tiny T}}$,
then the optimal predictor under model (8) is given by

$\hat{y}(\textbf{x}^{*})=g^{\mbox{\tiny
T}}(\textbf{x}^{*})\hat{\beta}+\mathbf{r}^{\mbox{\tiny
T}}(\textbf{x}^{*})\mathbf{R}^{-1}(\mathbf{y}-\mathbf{G}\hat{\beta}),$ (9)

where $\hat{\beta}=(\mathbf{G}^{\mbox{\tiny
T}}\mathbf{R}^{-1}\mathbf{G})^{-1}(\mathbf{G}^{\mbox{\tiny
T}}\mathbf{R}^{-1}\mathbf{y})$. If the assumed $\mu(\textbf{x})$ is close to
the truth, this formulation will lead to a better prediction than ordinary
krigging. Note that this universal kriging formulation uses $\mu(\textbf{x})$
to capture the known trends, but in most real applications, these trends are
not known, and hence ordinary kriging is commonly used [61].

###  2.2 Non-stationarity

Note that Equation (3) refers to a stationary GP, that is

$Cov\Big{(}Z(\textbf{x}+\textbf{h}),Z(\textbf{x})\Big{)}=\sigma^{2}R(\textbf{h}),$
(10)

where the correlation function $R(\textbf{h})$ is a positive semidefinite
function with $R(\mathbf{0})=1$ and $R(-\textbf{h})=R(\textbf{h})$. These
stationary Gaussian processes are popular surrogates for complex computer
models, since it can be shown that the corresponding predictor of $\mu$ in
equation (2)

$\hat{\mu}=(\mathbf{1}_{n}^{\mbox{\tiny
T}}\mathbf{R}^{-1}\mathbf{1}_{n})^{-1}\mathbf{1}_{n}^{\mbox{\tiny
T}}\mathbf{R}^{-1}\mathbf{y}$ (11)

is the best linear unbiased predictor (BLUP) in the sense that it minimizes
the mean squared prediction error. In reality this assumption of stationarity
may not hold. Under these circumstances, the above predictor is no longer
optimal. Some literature is available to deal with non-stationary Gaussian
processes for emulating computationally expensive functions. For example, [67]
introduced the idea of nonlinear mapping based on a parameterized density
function, and [20] proposed a Bayesian tree structure by dividing the design
space into subregions.

[2] used composite Gaussian process (CGP) models to address the
nonstationarity problem. In their formulation, the model takes the following
form:

$\displaystyle y(\mathbf{x})$ $\displaystyle=$ $\displaystyle Z_{\text{global
}}(\mathbf{x})+Z_{\text{local }}(\mathbf{x}),$ $\displaystyle Z_{\text{global
}}(\mathbf{x})$ $\displaystyle\sim$
$\displaystyle\operatorname{GP}\left(\mu,\tau^{2}R_{1}(\cdot)\right),$ (12)
$\displaystyle Z_{\text{local }}(\mathbf{x})$ $\displaystyle\sim$
$\displaystyle\mathrm{GP}\left(0,\sigma^{2}R_{2}(\cdot)\right).$

Here $Z_{\text{global }}(\mathbf{x})$ and $Z_{\text{local }}(\mathbf{x})$ are
two stationary GPs that are independent of each other. Just as the universal
kriging generalizes the ordinary kriging by adding a trend function
$\mu(\mathbf{x})$, the composite Gaussian process model given in equation
(2.2) is a further extension which adds a more flexible global trend
component. The model was extended to incorporate the non-constant variance
assumption as follows:

$\displaystyle y(\mathbf{x})$ $\displaystyle=$ $\displaystyle Z_{\text{global
}}(\mathbf{x})+\sigma(\textbf{x})Z_{\text{local }}(\mathbf{x}),$
$\displaystyle Z_{\text{global }}(\mathbf{x})$ $\displaystyle\sim$
$\displaystyle\operatorname{GP}\left(\mu,\tau^{2}R_{1}(\cdot)\right),$ (13)
$\displaystyle Z_{\text{local }}(\mathbf{x})$ $\displaystyle\sim$
$\displaystyle\mathrm{GP}\left(0,R_{2}(\cdot)\right).$

The model can be further extended for noisy data by adding a third GP (with
zero correlation) to the model (2.2).

###  2.3 Numeric Considerations - Local GP

Note that the prediction involves the inversion of the $n\times n$ correlation
matrix $\mathbf{R}$, where $n$ is the number of data points (see equation (5)
or (11), for example). This is a big hurdle in implementing GPs. To overcome
this problem, [19] introduced the idea of local Gaussian Process approximation
for large computer models. They provided a family of local sequential design
schemes that dynamically define the support points of a GP predictor based on
a local subset of the data. Their approach is different from that of
$k$-nearest neighbours. The basic idea is simple, under the standard choices
of the covariance structures the correlation between points is dependent on
the distance between those points, with data points far from $\textbf{x}^{*}$
having very little effect on its prediction. Hence it is not a good use of
computational resources to invert the full covariance matrix, as the elements
corresponding to “far away” points will contribute little to the prediction of
$y(\textbf{x}^{*})$. An interested reader should refer to [19] for the
formulas of the GP predictor based on a local subset of data. The end result
is a global predictor that takes advantage of modern multicore parallel
computing tools.

###  2.4 Extension to Qualitative Inputs

The conventional GP models consider quantitative predictor variables only, but
many computer experiments may have both quantitative and qualitative inputs.
In order to construct an emulator with qualitative factors, a naive approach
would be to create distinct GP models for data collected at the different
level combinations of the qualitative factors. Clearly this approach has many
limitations, particularly when there are several qualitative factors. There
are some more advanced techniques to deal with such cases. To fix ideas, for
an $n$-run computer model, denote the $k^{th}$ ($k=1,\ldots,n$) data input as
$\textbf{w}_{k}=(\textbf{x}_{k}^{\mbox{\tiny T}},\textbf{z}_{k}^{\mbox{\tiny
T}})^{\mbox{\tiny T}}$ where
$\textbf{x}_{k}=(x_{k1},\ldots,x_{kp})^{\mbox{\tiny T}}\in\mathbb{R}^{p}$ is
the quantitative part and $\textbf{z}_{k}=(z_{k1},\ldots,z_{kq})^{\mbox{\tiny
T}}\in\mathbb{N}^{q}$ is the qualitative part (coded in levels) of the input.
Note here that previously x denoted the input, which was entirely continuous.
However, now w denotes the entire input, with x referring to the continuous
part. For these kind of problems, a popular GP based model was introduced by
[44], among many others [23], [72], [53], [70] and [71]. Specifically, an
ordinary GP model with a multiplicative covariance function is considered (for
any two inputs $\textbf{w}_{1}$ and $\textbf{w}_{2}$):

$\hbox{Cov}(Z(\textbf{w}_{1}),Z(\textbf{w}_{2}))=\sigma^{2}\prod_{j=1}^{q}\tau^{(j)}_{z_{1j}z_{2j}}R(\textbf{x}_{1},\textbf{x}_{2}|\mbox{\boldmath${\theta}$}),$
(14)

where the parameter $\tau^{(j)}_{z_{1j}z_{2j}}$ represents the correlation
between two levels ($z_{1j}$ and $z_{2j}$) in the $j^{th}$ qualitative factor
$z^{(j)}$, and $R(\textbf{x}_{1},\textbf{x}_{2}|\mbox{\boldmath${\theta}$})$
is given before in equation (4). Different choices of
$\tau^{(j)}_{z_{1j}z_{2j}}$ lead to different types of correlation functions.
For example, an exchangable correlation function is obtained when
$\tau^{(j)}_{z_{1j}z_{2j}}$ is some constant between $0$ and $1$.
Alternatively, an additive GP model was proposed in [11], which adopts the
following covariance function:

$\hbox{Cov}(Z(\textbf{w}_{1}),Z(\textbf{w}_{2}))=\sum_{j=1}^{q}\sigma_{j}^{2}\tau^{(j)}_{z_{1j}z_{2j}}R(\textbf{x}_{1},\textbf{x}_{2}|\mbox{\boldmath${\theta}$}^{(j)}),$
(15)

where $\sigma_{j}^{2}$ and $\mbox{\boldmath${\theta}$}^{(j)}$ ($j=1,\ldots,q$)
are the process variance and correlation parameters, respectively,
corresponding to $z^{(j)}$.

The methods above do not have good physical interpretation of the correlation
structures. Motivated by this, [64] proposed an EzGP method based on ANOVA
(Analysis of Variance) ideas to jointly model the quantitative and qualitative
inputs:

$Y(\textbf{w})=\mu+G_{\textbf{z}}(\textbf{x}),$ (16)

which implies that for any level combination of z, $Y(\textbf{w})$ is a
Gaussian process. In particular, they considered

$G_{\textbf{z}}(\textbf{x})=G_{0}(\textbf{x})+G_{z^{(1)}}(\textbf{x})+\cdots+G_{z^{(q)}}(\textbf{x}),$
(17)

where $G_{0}$ and $G_{z^{(h)}}$ ($h=1,\ldots q$) are independent Gaussian
processes with mean zero and some covariance functions. Here, $G_{0}$ is a
standard GP taking only quantitative inputs x, which can be viewed as the base
GP reflecting the intrinsic relation between $y$ and x. On the other hand,
$G_{z^{(h)}}$’s can be viewed as an adjustment to the base GP by the impact of
the qualitative factor $z^{(h)}$ ($h=1,\ldots q$). This EzGP technique enjoys
some nice theoretical properties and is able to flexibly address heterogeneity
in computer models involving multiple qualitative factors. [64] also developed
two variants of the EzGP model to achieve computational efficiency for data
with high dimensionality and large sizes.

###  2.5 Calibration

The notion of calibration and sensitivity analysis is important in the context
of physical and computer experiments. Instead of observing the real physical
process, $y^{Real}$, we are only able to observe a process $y^{Field}$ as:

$y^{Field}(\textbf{x})=y^{Real}(\textbf{x})+\epsilon,$ (18)

where $\epsilon$ is the usual normal error. This $y^{Real}$ is approximated by
a computer model $y^{Model}$. Note that the computer model $y^{Model}$ not
only has the input variables x, but also some unknown parameters $\theta$,
called calibration parameters which are used to fine tune the model. Note that
these calibration parameters can be, for example, the correlation parameters
discussed above. The field data $y^{Field}$ is used mainly to learn more about
the real phenomenon $y^{Real}$. [30] proposed a Bayesian framework to address
this as follows:

$\displaystyle y^{Real}(\textbf{x})$ $\displaystyle=$ $\displaystyle
y^{Model}(\textbf{x},\theta)+b(\textbf{x})$ (19) $\displaystyle
y^{Field}(\textbf{x})$ $\displaystyle=$ $\displaystyle
y^{Model}(\textbf{x},\theta)+b(\textbf{x})+\epsilon,$

where $b(\textbf{x})$ is a functional discrepancy, called bias. [30] used
Bayesian methods to estimate the bias correction function and unknown
calibration parameter $\theta$ under a GP prior. An alternative to this
Bayesian approach is an iterative history matching algorithm such as the one
proposed by [55] for calibrating a galaxy formation model called GALFORM. This
is actually a hands-on process, which intelligently eliminates the implausible
points from the input (or parameter) space and returns a set of plausible
candidates for the parameters $\theta$. Recently, [5] used this algorithm for
calibrating hydrological time-series models and [46] further extended this
method with a more systematic approach, in which they discretize the target
response series on a few time points, and then iteratively apply the history
matching algorithm with respect to the discretized targets.

##  3 Design of Computer Experiments

The computer experiments under consideration have deterministic outputs, and
thus replicates at a given set of input settings should be avoided, as they do
not provide any further information about the response. Good designs for
computer experiments are then designs that are “space-filling” in some sense,
which make it easier to fit accurate surrogate models. We will next discuss a
few types of space-filling designs and examine techniques which can be used to
construct them.

###  3.1 LHD: Efficient Experimental Designs

Latin hypercube designs (LHDs) are $n\times d$ matrices whose columns are
permutations of numbers 1 to $n$ (or 0 to $n-1$) [39]. They have unique point
projections on every dimension and avoid replications, making them ideal for
determining which inputs to use for computer experiments [13]. For a given
number of runs and input size, an LHD can easily be generated in R:

⬇

1# Load an R library for finding LHDs

2library(LHD)

3# Generate a Latin Hypercube design with 10 runs and 2 factors

4lhd1 <- rLHD(10, 2)

While it is intuitive to favor a design that is space-filling, in practice it
is difficult to identify such designs for experiments with different number of
runs and inputs. One of the more common approaches is to use orthogonal or
nearly-orthogonal LHDs (OLHD). OLHDs minimize the correlations among the input
settings in the design [15] [52]. They can be obtained by minimizing a
correlation-based design criteria. For example, two of the most commonly used
criteria for OLHDs are the average absolute correlation (ave$(|r|)$) and the
maximum absolute correlation (max$|r|$):

$\displaystyle\mathop{\mathrm{ave}}\limits(|r|)$
$\displaystyle=\frac{2\sum_{s=1}^{d-1}\sum_{s^{\prime}=s+1}^{d}|r_{ss^{\prime}}|}{d(d-1)},$
(20) $\displaystyle\mathop{\mathrm{max}}\limits|r|$
$\displaystyle=\mathop{\mathrm{max}}\limits_{s,s^{\prime}}|r_{ss^{\prime}}|,$

where $r_{ss^{\prime}}$ is the correlation between the $s$th and
$s^{\prime}$th columns of the design. If the design is a true orthogonal LHD,
then ave$(|r|)=0$ and max$|r|=0$. For example, to generate an OLHD in R with 8
factors and 32 runs we can write:

⬇

1# Obtain an orthogonal latin hypercube design

2# Need n_factor = r * 2^(C+1)

3OLHD <- OLHD.S2010(C = 3, r = 2, type = ”even”)

The design can easily be verified to be orthogonal by examining:

⬇

1# All off diagonal elements are 0

2t(OLHD) %*% OLHD

However, for many combinations of run size and number of inputs an orthogonal
LHD does not exist, and thus a good design will be one with small ave$(|r|)$
and max$|r|$ values. Many algebraic construction methods have been proposed
for finding OLHDs, and they can also be found via searching algorithms using
ave$(|r|)$ or max$|r|$ as objective functions. Some specific results include
[69], who proposed techniques for constructing orthogonal LHDs with run-size
$n=2^{m}$ or $n=2^{m}+1$ where $m$ is an integer. [3] proposed to rotate the
$2^{d}$ factorial designs for constructing $d$-factor orthogonal LHDs where
$d$ must be some power of 2 and the run-size is $n=2^{d}$. For further
examples, please refer to [7], [49], [9], [35], [51] and [68]; see [59] for a
survey.

Figure 3: Some examples of 9-run 2-factor LHDs

While OLHDs are very commonly used, they are not guaranteed to be space-
filling; see design (a) in Figure 3 for an example [63]. In light of this,
various design optimality criteria have been developed related to measures of
space-filling.

####  3.1.1 Centered $L_{2}$-Discrepancy Criteria

[24] defined several discrepancy based criteria among which the centered
$L_{2}$-discrepancy (CD) is the most popular. The intuition behind the CD
criteria is that a space-filling design should have points spread out
uniformly in the whole design space or any sub-space of the design space. If
this is the case, for any rectangular region of the design space we examine,
the number of design points in that space should be proportional to the volume
of that space. The CD criteria is defined as,

$CD(D_{n})^{2}=\sum_{v\neq\varnothing}\int_{C^{v}}\left|\frac{\\#(D_{n_{v}},J_{x_{v}})}{n}-\text{Volume}(J_{x_{v}})\right|^{2}\mathrm{d}x,$
(21)

where $D_{n}$ is the $n$-run, $d$-factor, $q$-level design, $v$ is some non
empty subset of ${1,2,\ldots,q}$, $C^{v}$ is the subspace defined by the
coordinate indexes selected by $v$, $D_{n_{v}}$ is the projection of $D_{n}$
onto the subspace $C^{u}$, $x_{v}$ is the projection of vector
$x=(x_{1},x_{2},\ldots,x_{q})$ on to the subspace $C^{v}$, $J_{x}$ is the
chosen rectangle space defined by $x$, $J_{x_{u}}$ is the projection of
$J_{x}$ onto the subspace defined by $C^{u}$, $\\#(D_{n_{u}},J_{x_{u}})$ is
the total number of designs points in $D_{n_{u}}$ within the chosen area
defined by $J_{x_{u}}$, and $Volume(J_{x_{u}})$ is the volume of $J_{x_{u}}$.
For more details on the rationale of the CD criteria, see the Chapter 3 in
[13] for a survey.

####  3.1.2 Multi-objective Criteria

Another commonly-used metric for evaluating designs’ space-filling properties
is the maximin distance criterion [27]. This criteria favors designs with
maximum pairwise distances between inputs. Maximin designs are popular due to
their robustness, since the design criteria focuses on optimizing the worst
case scenario $-$ the closest pairwise distance between any two points. [43]
defined a computationally efficient scalar value for evaluating the maximin
distance criterion:

$\phi_{p}=\left(\sum_{i=2}^{n}\sum_{j=1}^{i-1}\frac{1}{u_{i,j}^{p}}\right)^{\frac{1}{p}},$
(22)

where $u_{i,j}$ is the distance between the $i$th and $j$th design points.
Designs with smaller $\phi_{p}$ values are more space-filling. For
sufficiently large $p$ (e.g. $p>15$), the $\phi_{p}$ criterion is
asymptotically identical to the true maximin distance criterion.

Due to the desirability of both the orthogonality and maximin properties, [28]
proposed a multi-objective criterion (denoted OMmcri) to generate orthogonal-
maximin LHDs (OMm LHDs), which act as a compromise between orthogonal and
maximin designs. The OMmcri criteria is given by,

$\text{OMmcri}(x,\omega)=\omega\rho^{2}+(1-\omega)\frac{(\phi_{p}-\phi_{p,lowerbound})}{(\phi_{p,upperbound}-\phi_{p,lowerbound})}.$
(23)

Here, $\phi_{p}$ is the maximin criteria value from Equation (22), $\rho$ is
the $ave(|r|)$ criteria value as defined in Equation (20), $\omega$ is a
weight value reflecting the tradeoff between the orthogonality and maximin
criteria, and $\phi_{p,lowerbound}$ and $\phi_{p,upperbound}$ are given by,

$\displaystyle\phi_{p,lowerbound}$
$\displaystyle=\left\\{{{n}\choose{2}}\left(\frac{\lceil\overline{u}\rceil-\overline{u}}{{\lfloor\overline{u}\rfloor}^{p}}-\frac{\overline{u}-\lfloor\overline{u}\rfloor}{{\lceil\overline{u}\rceil}^{p}}\right)\right\\}^{\frac{1}{p}},\text{
and}$ $\displaystyle\phi_{p,upperbound}$
$\displaystyle=\left(\sum_{i=1}^{n-1}\frac{n-i}{(id)^{p}}\right)^{\frac{1}{p}},$

respectively. Here $\overline{u}$ is the average distance between the design
points and $\lfloor\overline{u}\rfloor$ and $\lceil\overline{u}\rceil$ are the
largest integer smaller than $\overline{u}$ and the smallest integer larger
than $\overline{u}$.

Another popular class of efficient LHDs is the orthogonal array based LHDs
(OALHDs) by [54], where the levels in randomized orthogonal arrays (OAs) are
expanded to form LHDs. The OALHDs have desirable sampling and projection
properties, but they are not necessarily space-filling [63]; see designs (b)
and (c) in Figure 3 for some examples. [33] proposed to use a simulated
annealing algorithm to search for space-filling OALHDs, and [66] further
proposed to consider both level permutation and level expansion for generating
OALHDs. Some algebraic construction methods are also available for
constructing maximin LHDs for certain design sizes [65] [60].

####  3.1.3 Maxpro: Maximum projection designs

Space-filling LHDs, including CD and and maximin distance LHDs, focus on the
design’s properties in the full dimensional spaces. Yet, their space-filling
properties in some sub-spaces (projections) may not be adequate. [29] proposed
the maximum projection LHDs (Maxpro LHDs) that guarantee designs have space-
filling properties in all projections. The maximum projection criterion is
defined as

$\mathop{\mathrm{min}}\limits_{\textbf{X}}\psi(\textbf{X})=\Bigg{\\{}\frac{1}{{n\choose
2}}\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}\frac{1}{\prod_{s=1}^{d}(x_{is}-x_{js})^{2}}\Bigg{\\}}^{1/d}.$
(24)

Here, X is a $n\times d$ matrix where each row is an input to the computer
experiment, and the minimization is over all pairs of rows in X. Clearly a
design minimizing $\psi$ will have every pair of design points apart from each
other in all projections, justifying the name “Maxpro.”

###  3.2 Searching Algorithms for Generating Efficient LHDs with Flexible
Sizes

The metrics discussed above for evaluating designs such as the minimax
criteria provide a way of quantifying how “good” a design is in some sense. It
remains to determine how to actually construct designs that have a good value
of the criterion, which is a challenging problem in many situations. For many
such design problems, it is popular to use metaheuristic optimization
algorithms to find designs. Metaheuristic algorithms are often used to solve
problems in astronomy, see, for example, [8], [16], [41], [40], and [42]. They
can be applied to solve difficult problems such as clustering in complex data
[12], [26].

These algorithms are preferred due to their flexibility - in general they will
work with any objective function. For a more detailed review of metaheuristic
algorithms for finding designs, see [37]. Here we will focus on two of the
more commonly used approaches: Simulated Annealing and the Genetic Algorithm.

####  3.2.1 Simulated Annealing Algorithms

Simulated Annealing (SA) is one of the most widely used general probabilistic
optimization techniques [31]. The algorithm follows the annealing process in
metallurgy, in which materials are heated to a high temperature where their
properties change, and then are allowed to slowly cool. [43] adapted the
classic SA algorithm for finding maximin distance LHDs, and the approach can
easily be modified to search for other types of designs by using the other
optimality criteria defined in Section 3.1.

SA starts with a random LHD and then improves it via an element exchange
method, where two random elements from a random column in the design are
exchanged. If this exchange results in a more efficient design, then the
change is kept. If the exchange does not result in any improvement, the change
is kept with probability controlled by the current temperature (tuning
parameter). Allowing changes that do not improve the design helps the search
algorithm to escape local optima. The SA algorithm will iteratively repeat
this exchange procedure. After a certain number of rounds, the temperature
would be annealed to decrease (cool down) the probability of updating the
current design following the annealing schedule. We summarize a general SA
framework in Algorithm 1, where the target function $\Phi$ to be minimized can
be the optimality criterion defined in (20), (21), (22), (23) and (24) for the
orthogonal, CD, maximin, OMm and Maxpro LHDs, respectively.

In the Algorithm 1, the maximum number of iterations $N$ is recommended to be
around 500 according to the convergence analysis in [59]. The decreasing rate
for the current temperature $T$ is another important tuning parameter. A
larger rate will make $T$ decline faster, and thus lead to a faster stop of
the algorithm. Yet, it may also result in larger probability of missing the
true global optimum. Considering this trade-off, it is recommended to set $T$
between $0.05$ to $0.15$. The tuning parameter $S$ indicates the maximum
consecutive attempts the algorithm will try without improvements before
temperature reduces, and [43] recommends it to be around 5, depending on how
expensive the objective function is to evaluate.

It is straightforward to use Simulated Annealing to find designs in R. For
example:

⬇

1# 10 Runs, 2 inputs, 25 iterations of Simulated Annealing algorithm

2LHD_SA <- SA(n = 10, k = 2, N = 25)

Similarly, designs satisfying the multi-objective approach can be found by:

⬇

1# 10 Runs, 2 inputs, 25 iterations of Simulated Annealing algorithm

2# using multi-objective

3multi_obj_design <- SA2008(n = 10, k = 2, N = 25)

####  3.2.2 Genetic Algorithms

The genetic algorithm (GA) is a metaheuristic algorithm inspired by the
process of natural selection [25] [17]. The GA starts from a population of
randomly generated candidate solutions (designs), called chromosomes. The
population of chromosomes in each iteration is called a generation. For each
generation the objective function will be evaluated for each chromosome, with
the corresponding value being known as the fitness. The more fit chromosomes
will be allowed to survive to the next generation, while the less fit
chromosomes will be replaced by new offspring. These offspring are obtained by
selecting several chromosomes (called parents) and recombining their settings
using crossover and mutation techniques to produce offspring with potentially
better fitness.

[34] adapted the general GA framework for searching for maximin LHDs. Their
approach begins with random LHDs as the initial population. They then perform
a selection step in which the best half of the LHDs are allowed to survive to
the next generation. Then, a crossover step is performed in which random
columns in these survivors are exchanged with other survivors. Additionally,
to encourage diversity in the solutions and prevent the algorithm becoming
stuck in a local optima, a mutation step is performed in which two random
elements in a column are exchanged. Note that the current best chromosome is
excluded from this mutation in order to preserve the best current solution.
Finally, the fitness of the new population of LHDs is calculated, and the
process is repeated until the stopping criteria is satisfied. We include a
detailed description of the GA, along with the tuning parameters, in Algorithm
2.

Algorithm 1 Simulated Annealing for LHD

1: Choose values for the tuning parameters: the starting temperature $T$, the
number of attempts before lowering the temperature $S$, and the maximum number
of iterations $N$.

2: Set the counter index $C=1$.

3: Construct a random starting LHD X.

4: Select a column from X at random.

5: Exchange two randomly selected elements within this chosen column. Denote
the new design by $\textbf{X}_{new}$.

6: If $\Phi(\textbf{X}_{new})<\Phi(\textbf{X})$, then
$\textbf{X}=\textbf{X}_{new}$ (accept the new design). Otherwise, let
$\textbf{X}=\textbf{X}_{new}$ with probability
$\hbox{exp}\left\\{-\frac{\Phi(\textbf{X}_{new})-\Phi(\textbf{X})}{T}\right\\}$.

7: If $S$ attempts have passed since the last improvement, decrease the
temperature $T$ and repeat Steps 4$-$6\.

8: If $C<N$, increment $C$ and repeat Steps 4$-$7; Otherwise, terminate and
return X.

Algorithm 2 Genetic Algorithm for LHD

1: Set the probability of mutation, $p_{mut}$. Suggested setting is $1/(d-1)$
[34]. Set the maximum number of iterations $N$ and the counter index $C=1$.

2: Generate $m$ random $n\times d$ LHDs, denoted by $X_{1},\ldots,X_{m}$,
where $m$ is the population size (number of chromosomes). Here, $m$ must be an
even number.

3: Evaluate the objective function, $\Phi(X_{i})$, for $i=1,\ldots,m$.

4: Select ${\it survivors}$: order the $X_{i}$ by their objective function
values and select the best $\frac{m}{2}$ $X_{i}$ (with the smallest
$\frac{m}{2}$ $\Phi$ values), denoted by $X_{i}^{s}$ for
$i=1,\ldots,\frac{m}{2}$, WLOG.

5: Let $X_{b}^{s}=\mathop{\mathrm{argmin}}\limits_{i}\Phi(X_{i}^{s})$ (i.e.
$X_{b}^{s}$ is the best ${\it survivor}$)

6: for each $X_{i}^{s}$, excluding $X_{b}^{s}$, do

7: Randomly choose a column $j$ from $X_{b}^{s}$, and replace it with the
$j^{th}$ column from $X_{i}^{s}$.

8: end for

9: for each $X_{i}^{s}$, excluding $X_{b}^{s}$, do

10: Randomly choose a column $j$ from $X_{i}^{s}$, and replace it with the
$j^{th}$ column from $X_{b}^{s}$.

11: end for

12: Update $X_{i}$: let $X_{1}=X_{b}^{s}$ and the $X_{2},\ldots,X_{m/2}$ be
the design matrices obtained by steps 6$-$8\. Let $X_{m/2+1}=X_{b}^{s}$ and
$X_{m/2+2},\ldots,X_{m}$ be the design matrices generated by Steps 9$-$11\.

13: for each $X_{i}$ (except $X_{1}$) do

14: for each column $j$ of $X_{i}$ do

15: if $z<p_{mut}$ where $z\sim\text{Uniform}(0,1)$ then

16: Exchange two randomly selected elements in $j$.

17: end if

18: end for

19: end for

20: Calculate $\Phi(X_{i})$ for all $i$.

21: if $C\leq N$, set $C=C+1$ and repeat Steps 4-21; otherwise, stop the
algorithm.

It is also straightforward to use the GA to find space-filling designs in R.
For example:

⬇

1# 10 Runs, 2 inputs, 25 iterations of Genetic algorithm

2# OC is optimality criteria - phi_p is the maximin distance

3LHD_GA <- GA(n = 10, k = 2, N = 25, OC = ”phi_p”)

##  4 Summary and Conclusions

Sophisticated computer simulators allow scientists to test complex systems
which would be too expensive or completely impossible to assess otherwise.
These simulations are usually very time-consuming, and computationally cheap
surrogates are called for to facilitate the analysis and optimization of the
underlying system. Gaussian processes are popular choices for such surrogates
(or emulators). In order to effectively reap the benefits of utilizing the
surrogate, the simulator should be evaluated on a set of points chosen
efficiently. Latin hypercube designs have proven efficient for that purpose.

In this tutorial paper we discussed design criteria and subsequent
metaheuristic optimization strategies for finding designs that allow astronomy
researchers to extract the maximum benefit offered by Gaussian process
surrogate modeling. We provided an overview of model fitting using Gaussian
processes and identification of optimal Latin hypercube designs. Relevant R
codes have been used for illustration. Apart from the libraries discussed in
the paper, there are many other packages in R that can be used. Interested
readers may want to consider the laGP (Local Approximate Gaussian Process
Regression [21]), DiceKriging (Kriging Methods for Computer Experiments [47]),
GPfit (Gaussian Processes Modeling [38]) and SLHD (Maximin-Distance (Sliced)
Latin Hypercube Designs [1]) packages.

One consideration not covered in this tutorial paper is how to best utilize
Gaussian process models when the data sets are astronomically large. Such “big
data” may cause the estimation techniques to become quite slow, requiring
advanced techniques to speed up the estimation. This is a topic of active
research. For further details, see [36].

## References

  * [1] Ba, S. (2015), SLHD: Maximin-Distance (Sliced) Latin Hypercube Designs, URL:https://cran.r-project.org/web/packages/SLHD/index.html, R package version 2.1-1.
  * [2] Ba, S. and Joseph, V. R. (2012), Composite Gaussian process models for emulating expensive functions, _Ann. Appl. Stat._ , 6, 4, 1838–1860.
  * [3] Beattie, S. D. and Lin, D. K. J., (2005), A new class of Latin hypercube for computer experiments, _Contemporary Multivariate Analysis and Designs of Experiments, in Celebration of Prof. Kai-Tai Fang’s 65th Birthday. Singapore: World Scientific_ , 205–226.
  * [4] Behroozi, P., Wechsler, R., Hearin, A., and Conroy, C., (2019), UniverseMachine: The correlation between galaxy growth and dark matter halo assembly from z= 0$-$10, _Monthly Notices of the Royal Astronomical Society_ , 488, 3, 3143–3194, Oxford University Press.
  * [5] Bhattacharjeea, N., Ranjan, P., Mandal, A. and Tollner, E. W. (2019), A history matching approach for calibrating hydrological models, _Environmental and Ecological Statistics_ , 26, 1, 87–105.
  * [6] Bingham, D., Branin Function. _Virtual Library of Simulation Experiments_ , https://www.sfu.ca/~ssurjano/branin.html.
  * [7] Bursztyn, D. and Steinberg, D. M., (2002), Rotation designs: orthogonal first-order designs with higher order projectivity, _Applied Stochastic Models in Business and Industry_ , 18, 3, 197–206, Wiley Online Library.
  * [8] Charbonneau, P., (1995), Genetic algorithms in astronomy and astrophysics, _The Astrophysical Journal Supplement Series_ , 101, 309–334.
  * [9] Cioppa, T. M. and Lucas, T. W., (2007), Efficient nearly orthogonal and space-filling Latin hypercubes, _Technometrics_ , 49, 1, 45–55, Taylor & Francis.
  * [10] Dancik, G. M. (2020), mlegp: Maximum Likelihood Estimates of Gaussian Processes, URL:https://cran.r-project.org/web/packages/mlegp/index.html, R package version 3.1.8.
  * [11] Deng, X., Lin, C. D., Liu, K.-W. and Rowe, R. K., (2017), Additive Gaussian process for computer models with qualitative and quantitative factors, _Technometrics_ , 59, 3, 283–292, Taylor & Francis.
  * [12] Djorgovski, SG., Brunner, R., Mahabal, A., Williams, R., Granat, R., and Stolorz, P., (2003), Challenges for cluster analysis in a virtual observatory, Statistical Challenges in Astronomy, 127–141, Springer.
  * [13] Fang, K. T., Li, R. and Sudjianto, A., (2006), Design and modeling for computer experiments, CRC Press.
  * [14] Forrester, A., Sobester, A., and Keane, A., (2008) Engineering design via surrogate modelling: a practical guide, John Wiley & Sons.
  * [15] Georgiou, S. D., (2009), Orthogonal Latin hypercube designs from generalized orthogonal designs, _Journal of Statistical Planning and Inference_ , 139, 4, 1530–1540, Elsevier.
  * [16] Giuliano, M., and Johnston, M., (2008), _Multi-Objective Evolutionary Algorithms for Scheduling the James Webb Space Telescope._ , ICAPS, 107–115.
  * [17] Goldberg, D. E., (1989), Genetic algorithms in search, _Optimization and MachineLearning_ , Addison Wesley Publishing Co. Inc.
  * [18] Gramacy, R. B. (2020). Surrogates: Gaussian Process Modeling, Design, and Optimization for the Applied Sciences. CRC Press.
  * [19] Gramacy, R. B. and Apley, D. W. (2015), Local Gaussian Process Approximation for Large Computer Experiments, _Journal of Computational and Graphical Statistics_ , 24, 2.
  * [20] Gramacy, R. B. and Lee, H. K. H. (2008), Bayesian treed Gaussian process models with an application to computer modeling. _J. Amer. Statist. Assoc._ , 103, 1119–-1130.
  * [21] Gramacy, R. B., Sun, F. (2019), laGP: Local Approximate Gaussian Process Regression, URL:https://cran.r-project.org/web/packages/laGP/index.html, R package version 1.5-5.
  * [22] Gration, A. and Wilkinson, M., (2019), Dynamical modelling of dwarf spheroidal galaxies using Gaussian-process emulation, _Monthly Notices of the Royal Astronomical Society_ , 485, 4, 4878–4892, Oxford University Press.
  * [23] Han, G., Santner, T. J., Notz, W. I., Bartel, D. L., (2009), Prediction for computer experiments having quantitative and qualitative input variables, _Technometrics_ , 51, 3, 278–288, Taylor & Francis.
  * [24] Hickernell, F., (1998) A generalized discrepancy and quadrature error bound, _Mathematics of Computation of the American Mathematical Society_ , 67, 221, 299–322.
  * [25] Holland, J. H. and others, (1992), Adaptation in natural and artificial systems: an introductory analysis with applications to biology, control, and artificial intelligence, MIT press.
  * [26] Hruschka, E., Campello, R., Freitas, A., and others, (2009) A survey of evolutionary algorithms for clustering, _IEEE Transactions on Systems, Man, and Cybernetics, Part C (Applications and Reviews)_ , 39, 2, 133–155, IEEE
  * [27] Johnson, M. E., Moore, L. M. and Ylvisaker, D., (1990), Minimax and maximin distance designs, _Journal of Statistical Planning and Inference_ , 26, 2, 131–148, Elsevier.
  * [28] Joseph, V. R. and Hung, Y., (2008), Orthogonal-maximin Latin hypercube designs, _Statistica Sinica_ , 171–186, JSTOR.
  * [29] Joseph, V. R., Gul, E. and Ba, S., (2015), Maximum projection designs for computer experiments, _Biometrika_ , 102, 2, 371–380, Oxford University Press.
  * [30] Kennedy, M. and O’Hagan, A. (2002), Bayesian calibration of computer models, _Journal of the Royal Statistical Society Series B (Statistical Methodology)_. 63, 3, 425–464, Wiley Online Library.
  * [31] Kirkpatrick, S., Gelatt, C. D. and Vecchi, M. P., (1983), Optimization by simulated annealing, _Science_ , 220, 4598, 671–680, American Association for the Advancement of Science.
  * [32] Kruckow, M., Tauris, T., Langer, N., Kramer, M., and Izzard, Robert G, (2018), Progenitors of gravitational wave mergers: binary evolution with the stellar grid-based code COMBINE, _Monthly Notices of the Royal Astronomical Society_ , 481,2, 1908–1949, Oxford University Press.
  * [33] Leary, S., Bhaskar, A., and Keane, A., (2003), Optimal orthogonal-array-based latin hypercubes, _Journal of Applied Statistics_ , 30, 5, 585–598, Taylor & Francis.
  * [34] Liefvendahl, M. and Stocki, R., (2006), A study on algorithms for optimization of Latin hypercubes, _Journal of Statistical Planning and Inference_ , 136, 9, 3231–3247, Elsevier.
  * [35] Lin, C. D., Mukerjee, R. and Tang, B., (2009), Construction of orthogonal and nearly orthogonal Latin hypercubes, _Biometrika_ , 243–247, JSTOR.
  * [36] Liu, H., Ong, Y-S., Shen, X. and Cai, J., (2020), When Gaussian process meets big data: A review of scalable GPs, _IEEE transactions on neural networks and learning systems_ , 31, 11, 4405–4423, IEEE.
  * [37] Mandal, A., Wong, W. K., and Yu, Y. (2015), Algorithmic searches for optimal designs, _Handbook of design and analysis of experiments_ ,755–783, CRC Press Boca Raton, FL.
  * [38] MacDoanld, B., Chipman, H., Campbell, C. and Ranjan, P. (2019), GPfit: Gaussian Processes Modeling, URL:https://cran.r-project.org/web/packages/GPfit/index.html, R package version 1.0-8.
  * [39] McKay, M. D. and Beckman, R. J. and Conover, W. J., (1979), Comparison of three methods for selecting values of input variables in the analysis of output from a computer code, _Technometrics_ , 21, 2, 239–245, Taylor & Francis.
  * [40] Misiak, M., and others, (2016) Evolutionary Algorithms in Astrodynamics, _International Journal of Astronomy and Astrophysics_ , 6, 04, 435–439, Scientific Research Publishing.
  * [41] Mohanty, S., (2012) Particle Swarm Optimization and regression analysis–I, _Astronomical Review_ , 7, 2, 29–35, Taylor & Francis
  * [42] Mohanty, S. and Fahnestock, E., (2020), Adaptive spline fitting with particle swarm optimization, _Computational Statistics_ , 36, pages=155–191, Springer.
  * [43] Morris, M. D. and Mitchell, T. J., (1995), Exploratory designs for computational experiments, _Journal of statistical planning and inference_ , 43, 3, 381–402, Elsevier.
  * [44] Qian, P. Z. G. and Wu, H. and Wu, C. F. J., (2008), Gaussian process models for computer experiments with qualitative and quantitative factors, _Technometrics_ , 50, 3, 383–396, Taylor & Francis.
  * [45] R Core Team, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria, 2019, https://www.R-project.org/.
  * [46] Resch, J., Mandal, A. and Ranjan, P. (2021), “Inverse problem for dynamic computer simulators via multiple scalar-valued contour estimation”, https://arxiv.org/abs/2010.08941.
  * [47] Roustant, O., Ginsbourger, D., Deville, Y., Clement, C. and Richet, Y. (2020), DiceKriging: Kriging Methods for Computer Experiments, URL:https://cran.r-project.org/web/packages/DiceKriging/index.html, R package version 1.5.8.
  * [48] Sacks, J., Welch, W. J., Mitchell, T. J., and Wynn, H. P. (1989), Design and analysis of computer experiments. _Statistical science_ , 409–423.
  * [49] Steinberg, D. M. and Lin, D. K. J., (2006), A construction method for orthogonal Latin hypercube designs, _Biometrika_ , 93, 2, 279–288, JSTOR.
  * [50] Stevenson, S., Vigna-Gómez, A., Mandel, I., Barrett, J.W., Neijssel, C.J., Perkins, D. and De Mink, S.E., (2017), Formation of the first three gravitational-wave observations through isolated binary evolution. Nature Communications, 8, 1, 1–7. Vancouver.
  * [51] Sun, F., Liu, M.-Q. and Lin, D. K. J., (2010), Construction of orthogonal Latin hypercube designs with flexible run sizes, _Journal of Statistical Planning and Inference_ , 140, 11, 3236–3242, Elsevier.
  * [52] Sun, F. and Tang, B., (2017), A general rotation method for orthogonal Latin hypercubes, _Biometrika_ , 104, 2, 465–472, Oxford University Press.
  * [53] Swiler, L. P. and Hough, P. D. and Qian, P. Z. G., Xu, X., Storlie, C. and Lee, H., (2014), Surrogate models for mixed discrete-continuous variables, Constraint Programming and Decision Making, 181–202, Springer.
  * [54] Tang, B., (1993), Orthogonal array-based Latin hypercubes, _J. Amer. Statist. Assoc._ , 88, 424, 1392–1397, Taylor & Francis.
  * [55] Vernon, I., Goldstein, M. and Bower, R. G., (2010), Galaxy formation: a Bayesian uncertainty analysis, _Bayesian Analysis_ , 5, 4, 619–669, International Society for Bayesian Analysis.
  * [56] Vigna-Gómez, A., Neijssel, C.J., Stevenson, S., Barrett, J.W., Belczynski, K., Justham, S., de Mink, S.E., Müller, B., Podsiadlowski, P., Renzo, M. and Szécsi, D., 2018. On the formation history of Galactic double neutron stars. _Monthly Notices of the Royal Astronomical Society_ , 481, 3, 4009–4029. Vancouver.
  * [57] Wackernagel H. (2002), Multivariate Geostatistics, Springer.
  * [58] Wang, H., Xiao, Q., and Mandal, A., (2020), LHD: Latin Hypercube Designs (LHDs), URL:https://CRAN.R-project.org/package=LHD, R package version 1.3.1.
  * [59] Wang, H., Xiao, Q. and Mandal, A., (2020) Musings about Constructions of Efficient Latin Hypercube Designs with Flexible Run-sizes, _arXiv preprint arXiv:2010.09154v2_.
  * [60] Wang, L., Xiao, Q. and Xu, H., (2018), Optimal maximin $L_{1}$-distance Latin hypercube designs based on good lattice point designs, _Annals of Statistics_ , 46, 6B, 3741–3766, Institute of Mathematical Statistics.
  * [61] Welch, W. J., Buck, R. J., Sacks, J., Wynn, H. P., Mitchell, T. J. and Morris, M. D. (1992), Screening, predicting, and computer experiments, _Technometrics_ , 34, 15–-25.
  * [62] Williams, D., Heng, I.S., Gair, J., Clark, J.A. and Khamesra, B., (2019), A Precessing Numerical Relativity Waveform Surrogate Model for Binary Black Holes: A Gaussian Process Regression Approach. arXiv preprint arXiv:1903.09204.
  * [63] Xiao, Q., (2017), Constructions and Applications of Space-Filling Designs, Ph.D. Dissertation, University of California Los Angeles.
  * [64] Xiao, Q., Mandal, A., Lin, C. D., and Deng, X. (2021), EzGP: Easy-to-interpret Gaussian Process models for computer experiments with both quantitative and qualitative factors. under revision for _SIAM/ASA Journal on Uncertainty Quantification_.
  * [65] Xiao, Q. and Xu, H., (2017), Construction of maximin distance Latin squares and related Latin hypercube designs, _Biometrika_ , 104, 2, 455–464, Oxford University Press.
  * [66] Xiao, Q. and Xu, H., (2018), Construction of maximin distance designs via level permutation and expansion, _Statistica Sinica_ , 28, 3, 1395–1414, JSTOR.
  * [67] Xiong, Y., Chen, W., Apley, D. W. and Ding, X. (2007), A non-stationary covariance-based kriging method for metamodelling in engineering design. _Internat. J. Numer. Methods Engrg._ , 71, 733-–756.
  * [68] Yang, J. and Liu, M.-Q., (2012), Construction of orthogonal and nearly orthogonal Latin hypercube designs from orthogonal designs, _Statistica Sinica_ , 433–442, JSTOR.
  * [69] Ye, K. Q., (1998), Orthogonal column Latin hypercubes and their application in computer experiments, _Journal of the American Statistical Association_ , 93, 444, 1430–1439, Taylor & Francis.
  * [70] Zhang, Y. and Notz, W. I., (2015), Computer experiments with qualitative and quantitative variables: a review and reexamination, _Quality Engineering_ , 27, 1, 2–13, Taylor & Francis.
  * [71] Zhang, Y., Tao, S., Chen, W. and Apley, D. W., (2019), A latent variable approach to Gaussian process modeling with qualitative and quantitative factors, _Technometrics_ , 1–12, Taylor & Francis.
  * [72] Zhou, Q. and Qian, P. Z. G. and Zhou, S., (2011), A simple approach to emulation for computer models with qualitative and quantitative factors, _Technometrics_ , 53, 3, 266–273, Taylor & Francis.


# Convolutional Learning on Directed Acyclic Graphs

Samuel Rey∗, Hamed Ajorlou† and Gonzalo Mateos†
∗Dept. of Signal Theory and Communications, Rey Juan Carlos University,
Madrid, Spain
Email<EMAIL_ADDRESS>
†Dept. of Electrical and Computer Engineering, University of Rochester, NY,
USA.
Emails: {gmateosb<EMAIL_ADDRESS>This work was supported in part
by the Spanish AEI under Grants PID2019-105032GBI00, TED2021-130347B-I00 and
PID2022-136887NBI00 funded by MCIN/AEI/10.13039/501100011033.

###### Abstract

We develop a novel convolutional architecture tailored for learning from data
defined over directed acyclic graphs (DAGs). DAGs can be used to model causal
relationships among variables, but their nilpotent adjacency matrices pose
unique challenges towards developing DAG signal processing and machine
learning tools. To address this limitation, we harness recent advances
offering alternative definitions of causal shifts and convolutions for signals
on DAGs. We develop a novel convolutional graph neural network that integrates
learnable DAG filters to account for the partial ordering induced by the graph
topology, thus providing valuable inductive bias to learn effective
representations of DAG-supported data. We discuss the salient advantages and
potential limitations of the proposed DAG convolutional network (DCN) and
evaluate its performance on two learning tasks using synthetic data: network
diffusion estimation and source identification. DCN compares favorably
relative to several baselines, showcasing its promising potential.

###### Index Terms:

DAG, graph signal processing, graph neural networks, graph convolution.

##  I Introduction

The vast complexity and heterogeneity of modern networks are propelling the
generation of data represented on non-Euclidean domains. Consequently, fields
such as graph signal processing (GSP) have emerged, which model the underlying
irregular signal structure as a graph and then exploit its topology to process
said relational data [1, 2]. Graph neural networks (GNN) are non-linear,
graph-aware machine learning (ML) models that serve as a prominent example of
this new paradigm [3, 4, 5, 6]. GNNs exploit prior information in the graph
via convolutions [7, 8, 9], attention mechanisms [10], or graph autoencoders
[11, 12], and have attained state-of-the-art performance in a gamut of ML
applications.

Arguably most GNNs and graph-based methods focus on undirected graphs.
However, accounting for directionality can play an important role in
processing information, but this leap comes with well-documented challenges
[13, 14, 15] that are exacerbated when dealing with a directed acyclic graph
(DAG). For instance, DAGs have nilpotent adjacency matrices and hence a
collapsed spectrum that prevents us from directly applying spectral-based
tools [16]. These difficulties notwithstanding, DAGs are prevalent across
domains including causal inference [17], learning Bayesian networks and
structural causal models [18, 19, 20, 21], code parsing [22], or performance
prediction for neural architectures [23], just to name a few.

This work develops a novel _convolutional_ GNN tailored for learning from
signals defined over DAGs. As DAGs impose a partial ordering on the set of
nodes, a key challenge is to incorporate this stronger inductive bias into our
architecture in a principled manner. To that end, we build on a recent
framework that extends GSP tools to partially ordered sets (posets) [24, 16],
to arrive at a DAG convolutional network (DCN). These foundational design
principles can be traced to linear structural equation models (SEMs), making
DCN well-suited to capture causal relations in the data. Moreover, the
resulting architecture admits a spectral representation and is permutation
equivariant. Unlike related alternatives [25, 26], the convolutional layers
primarily involve sparse matrix multiplications, resulting in manageable
computational complexity.

Related work and contributions in context. Developing GNNs to process data
defined over DAGs is an important problem that is starting to draw attention.
Recently, [25] advocated embedding DAGs into a smooth latent space and
introduced a variational autoencoder for DAGs, dubbed D-VAE. The DAG neural
network (DAGNN) from [26] sequentially aggregates the node representations
from predecessors via an attention mechanism, and then uses a gated recurrent
unit (GRU) to encode the causal structure of the DAG during a combination
step. However, the sequential scheme followed by both D-VAE and DAGNN results
in a large computational burden, limiting their applicability to moderately-
sized graphs. Furthermore, DAGNN needs to combine the true and reverse
directions of the DAG in a non-intuitive way to process the data. Later, [27]
proposed encoding the DAG ordering via a reachability-based attention
mechanism, extending graph transformers to DAGs.

All in all, the main contributions of this work are as follows:

  * •

We propose DCN, the first convolutional GNN designed to learn from data
defined over DAGs. We elaborate on the key benefits of the model and discuss
its limitations.

  * •

Different from existing approaches [25, 26], DCN layers rely on formal
convolutions derived in [16] for a signal model over posets, that can be
linked to linear SEMs.

  * •

We evaluate the performance of DCN and show it compares favorably to several
baselines in different graph learning tasks. We share the code to reproduce
our results.

##  II Preliminaries and problem statement

We first define basic concepts and notation about DAGs, as well as the
required background on GSP and GNNs. We then formally state the problem of
learning from DAG signals.

###  II-A Preliminaries: DAGs, graph signals and GNNs

DAGs and graph signals. Let ${\mathcal{D}}=({\mathcal{V}},{\mathcal{E}})$ be a
directed acyclic graph (DAG), where ${\mathcal{V}}$ denotes the set of $N$
nodes and ${\mathcal{E}}\subseteq{\mathcal{V}}\times{\mathcal{V}}$ denotes the
set of edges. An edge $(i,j)\in{\mathcal{E}}$ exists if and only if there is a
link from node $j$ to node $i$. Since ${\mathcal{D}}$ contains no cycles, we
can sort ${\mathcal{V}}$ topologically, ensuring that node $j$ precedes node
$i$ whenever $(i,j)\in{\mathcal{E}}$. Consequently, every DAG induces a unique
partial order on ${\mathcal{V}}$, where $j<i$ if $j$ is a _predecessor_ of
$i$, i.e., if there exists a path from $j$ to $i$ [24]. Note that in posets
not all the elements are comparable. In fact, for every pair
$i,j\in{\mathcal{V}}$ without a path from $i$ to $j$ or vice-versa, we have
that $i\not\leq j$ and $j\not\leq i$.

Let ${\mathbf{A}}\in{\mathbb{R}}^{N\times N}$ denote the (possibly weighted)
adjacency matrix of ${\mathcal{D}}$, where $A_{ij}\neq 0$ if and only if
$(i,j)\in{\mathcal{E}}$. Due to the aforementioned partial ordering of
${\mathcal{V}}$, we can sort the entries of ${\mathbf{A}}$ to obtain a
strictly lower-triangular matrix. The diagonal is null because DAGs cannot
have self-loops.

In addition to the graph ${\mathcal{D}}$, we consider signals defined on the
set of nodes ${\mathcal{V}}$. Formally, graph signals are functions from the
vertex set to the real field $x:{\mathcal{V}}\mapsto{\mathbb{R}}$, which can
be represented as a vector ${\mathbf{x}}\in{\mathbb{R}}^{N}$, where $x_{i}$ is
the signal value at node $i$. A signal processing theory can be constructed
for graph signals supported on DAGs (a.k.a. DAG signals), which was put forth
in [16] and we review in Section III-A.

Convolutional GNNs. A GNN can be represented as a parametric non-linear
function that depends on the underlying graph structure. While several GNN
models have been proposed, most of them build upon an _aggregation function_
driven by the graph topology. In this work we consider a class of GNNs where
the aggregation is accomplished through graph convolutions [7, 9]. More
precisely, we are interested in GNNs employing a bank of learnable
convolutional graph filters, where the output of layer $\ell=1,\ldots,L$ is
given by

${\mathbf{X}}^{(\ell+1)}=\sigma\left(\sum_{r=0}^{R-1}{\mathbf{A}}^{r}{\mathbf{X}}^{(\ell)}\boldsymbol{\Theta}_{r}^{(\ell)}\right).$
(1)

Here, $\boldsymbol{\Theta}_{r}^{(\ell)}\in{\mathbb{R}}^{F_{i}^{(\ell)}\times
F_{o}^{(\ell)}}$ collects the learnable filter coefficients, $F_{i}^{(\ell)}$
and $F_{o}^{(\ell)}$ are the number of input and output features of the
$\ell$-th layer, and the powers ${\mathbf{A}}^{r}$ determine the radius of the
aggregation neighborhood. If necessary, the adjacency matrix could be
substituted with any desired _graph-shift operator_ (GSO) that encodes the
graph structure; see e.g., [1, 9]. The point-wise nonlinear activation
function is often chosen to be a ReLU, i.e., $\sigma(x)=\max(0,x)$.

###  II-B Problem statement

This work addresses the problem of learning from signals defined over a DAG
${\mathcal{D}}$. We are given
${\mathcal{T}}=\\{{\mathbf{X}}_{m},{\mathbf{y}}_{m}\\}_{m=1}^{M}$, the
training set containing $M$ input-output observations from a network process
on ${\mathcal{D}}$. The input graph signals are denoted as
${\mathbf{X}}_{m}\in{\mathbb{R}}^{N\times F}$ and the corresponding outputs as
${\mathbf{y}}_{m}\in{\mathbb{R}}^{N}$. Here outputs are single-feature graph
signals ${\mathbf{y}}_{m}$ for simplicity, and generalizations to multi-
feature signals or nodal/graph labels are straightforward. We use
${\mathcal{T}}$ to learn the map relating ${\mathbf{X}}_{m}$ and
${\mathbf{y}}_{m}$, which is assumed to be accurately represented by a non-
linear parametric function
$f_{\boldsymbol{\Theta}}(\cdot|{\mathcal{D}}):{\mathbb{R}}^{N\times
F}\mapsto{\mathbb{R}}^{N}$. To that end, we estimate the weights
$\boldsymbol{\Theta}$ of $f_{\boldsymbol{\Theta}}(\cdot|{\mathcal{D}})$ by
solving

$\min_{\boldsymbol{\Theta}}\frac{1}{M}\sum_{m=1}^{M}{\mathcal{L}}({\mathbf{y}}_{m},f_{\boldsymbol{\Theta}}({\mathbf{X}}_{m}|{\mathcal{D}})),$
(2)

where ${\mathcal{L}}$ is a loss function, e.g., mean squared error (MSE) for
regression or the cross entropy loss for classification.

##  III Convolutional learning on DAGs

Among the many graph-aware functions available, we confine our selection of
$f_{\boldsymbol{\Theta}}(\cdot|{\mathcal{D}})$ to the class of graph
convolutional neural networks (GCNN); see e.g., [7, 9] for extensive
discussion and justification on their merits. The first step towards designing
a GCNN tailored to DAG signals is to define the convolution operation.
Naively, one could use the graph filters in (1), as it is customary when
graphs have cycles [28, 9]. However, for a DAG all the eigenvalues of
${\mathbf{A}}$ are 0, due to its strict lower triangular structure. This
collapsed spectrum deprives us of a spectral interpretation for such filters.
Moreover, recall that DAGs impose a partial ordering and can potentially
capture causal relations, two properties we wish to encode in our GCNN. To
overcome these challenges and design
$f_{\boldsymbol{\Theta}}(\cdot|{\mathcal{D}})$ in a principled manner (Section
III-B), we bring to bear the novel DAG signal processing framework introduced
in [16], which we briefly outline next.

###  III-A Graph-shift operators and convolution for DAGs

Let ${\mathbf{c}}\in{\mathbb{R}}^{N}$ be a vector collecting the causes or
nodal contributions of a DAG signal ${\mathbf{x}}$. The signal model
introduced in [16] postulates that ${\mathbf{x}}$ can be described by the
causes at predecessor nodes. Specifically, we write
${\mathbf{x}}={\mathbf{W}}{\mathbf{c}}$, where
${\mathbf{W}}\in{\mathbb{R}}^{N\times N}$ is the _weighted transitive closure_
of ${\mathcal{D}}$ with $W_{ij}\neq 0$ if there is a path from $j$ to $i$
(i.e., if $j<i$) or $i=j$. Note that ${\mathbf{W}}$ can be interpreted as the
adjacency matrix of the reachability graph of ${\mathcal{D}}$, with ones on
the diagonal to account for the reflexivity property of posets. Among the
different methods to compute ${\mathbf{W}}$ from ${\mathbf{A}}$ [29], we focus
on ${\mathbf{W}}=({\mathbf{I}}-{\mathbf{A}})^{-1}$. This is equivalent to a
linear SEM for ${\mathbf{x}}$, with exogenous inputs ${\mathbf{c}}$ [16,
Section IV-D].

Given this model, every node $k\in{\mathcal{V}}$ induces a causal GSO on the
signal ${\mathbf{x}}$ given by the matrix
${\mathbf{T}}_{k}\in{\mathbb{R}}^{N\times N}$, such that

$[{\mathbf{T}}_{k}{\mathbf{x}}]_{i}=\sum_{j\leq i\;\mathrm{and}\;j\leq
k}W_{ij}c_{j}.$ (3)

Intuitively, the signal value at node $i$ after shifting ${\mathbf{x}}$ with
respect to node $k$ results in a signal containing only causes from common
predecessors to both $i$ and $k$. This can be compactly written in matrix form
as

${\mathbf{T}}_{k}{\mathbf{x}}={\mathbf{W}}{\mathbf{D}}_{k}{\mathbf{c}}={\mathbf{W}}{\mathbf{D}}_{k}{\mathbf{W}}^{-1}{\mathbf{x}},$
(4)

where ${\mathbf{D}}_{k}$ is a diagonal matrix whose entry
$[{\mathbf{D}}_{k}]_{ii}=1$ if $i\leq k$ (i.e., if node $i$ is a predecessor
of node $k$) and 0 otherwise. The inverse ${\mathbf{W}}^{-1}$ always exists
since ${\mathbf{W}}$ is a full-rank matrix and it can be efficiently computed
using the weighted Moebius inversion [30]. It follows that the columns of
${\mathbf{W}}$ and the diagonal entries of ${\mathbf{D}}_{k}$ are,
respectively, the eigenvectors and eigenvalues of the GSO ${\mathbf{T}}_{k}$.
Matrix ${\mathbf{W}}$ is thus a Fourier basis for DAG signals, while
${\mathbf{W}}^{-1}$ is the corresponding Fourier transform. Because
${\mathbf{x}}={\mathbf{W}}{\mathbf{c}}$, for this DAG signal model the causes
${\mathbf{c}}$ are also spectral coefficients.

In GSP, a graph filter is a polynomial of the (single) GSO. In contrast, now
we have a GSO ${\mathbf{T}}_{k}$ for every node in ${\mathcal{V}}$ and,
moreover, every ${\mathbf{T}}_{k}$ is an idempotent operator. Consequently,
the most general shift-invariant DAG filter ${\mathbf{H}}$ is given by

${\mathbf{H}}=\sum_{k\in{\mathcal{V}}}h_{k}{\mathbf{T}}_{k}={\mathbf{W}}\sum_{k\in{\mathcal{V}}}h_{k}{\mathbf{D}}_{k}{\mathbf{W}}^{-1},$
(5)

resulting in a convolution operation
${\mathbf{h}}*{\mathbf{x}}={\mathbf{H}}{\mathbf{x}}$. Here,
${\mathbf{h}}\in{\mathbb{R}}^{N}$ collects the filter coefficients $h_{k}$,
and the frequency response of ${\mathbf{H}}$ is given by the diagonal of
$\sum_{k\in{\mathcal{V}}}h_{k}{\mathbf{D}}_{k}$. This definition of DAG filter
plays a central role in our convolutional architecture, the subject dealt with
next.

###  III-B DAG Convolutional Network (DCN)

Now, we are ready to introduce our DAG-aware convolutional network, which we
refer to as DCN. The proposed DCN harnesses the definition of causal graph
filters in (5) to implement convolutions, which we compose with a point-wise
nonlinearity to process DAG signals. In the simplest case, the output of layer
$\ell=1,\ldots,L$ can be computed as

${\mathbf{x}}^{(\ell+1)}=\sigma\left(\sum_{k\in{\mathcal{V}}}h_{k}^{(\ell)}{\mathbf{T}}_{k}{\mathbf{x}}^{(\ell)}\right),$
(6)

where $h_{k}^{(\ell)}$ are the learnable filter coefficients.

We refer to (6) as a DAG perceptron, and we can gain intuition from two
different standpoints. In the spectral domain, recall that
${\mathbf{T}}_{k}{\mathbf{x}}^{(\ell)}={\mathbf{W}}{\mathbf{D}}_{k}{\mathbf{c}}^{(\ell)}$.
In words, the $k$-th GSO selects the causes (i.e., the exogenous input in an
SEM) of predecessors of the node $k$, and diffuses them across the
reachability graph via ${\mathbf{W}}$. Clearly, this accounts for the partial
ordering imposed by the DAG. Alternatively, one could look at (6) in the
context of message-passing networks. For every node $i\in{\mathcal{V}}$, each
${\mathbf{T}}_{k}$ forms a message combining the features from predecessors
common to $i$ and $k$, to finally yield $x_{i}^{(\ell+1)}$ as a non-linear
weighted combination of these messages. A complete example and schematic
illustrating these layer operations will be included in the extended journal
version.

While (6) is intuitive, a single filter may not be expressive enough to learn
complex DAG signal representations. Furthermore, input signals typically
comprise several features. We address these limitations following a rationale
similar to that of (1). At each layer, the single filter is replaced with a
bank of learnable DAG filters, resulting in the recursion [cf. (1)]

${\mathbf{X}}^{(\ell+1)}=\sigma\left(\sum_{k\in{\mathcal{V}}}{\mathbf{T}}_{k}{\mathbf{X}}^{(\ell)}\boldsymbol{\Theta}_{k}^{(\ell)}\right).$
(7)

Here, the matrices
$\boldsymbol{\Theta}_{k}^{(\ell)}\in{\mathbb{R}}^{F_{i}\times F_{o}}$ collect
the learnable filter coefficients with $F_{i}$ and $F_{o}$ denoting the number
of input and output features. Then, for a DCN with $L$ layers, the prediction
${\hat{\mathbf{y}}}_{m}={\mathbf{X}}^{(L)}=f_{\boldsymbol{\Theta}}({\mathbf{X}}_{m}|{\mathcal{D}})$
will be the output of the final layer.

###### Remark 1

The DCN whose output is given by the recursion (7) is permutation equivariant.
Due to space limitations, the proof is deferred to the extended journal
version of this work.

Discussion. The DCN convolutional layer (7) is closely related to GCNN
architectures for general (cyclic) graphs. However, relying on causal graph
filters is imperative to implement formal convolutions, and offers some key
advantages. First, the spectrum of ${\mathbf{T}}_{k}$ is well defined,
endowing the DCN with a spectral representation that is fundamental to analyze
properties such as stability, transferability, or denoising capability [9,
31], an interesting future research direction. Furthermore, since the
eigenvalues of ${\mathbf{T}}_{k}$ are the binary matrices ${\mathbf{D}}_{k}$,
no issues stemming from numerical instability are expected when concatenating
several layers.

On the flip side, the number of GSOs potentially involved in the convolution –
hence the number of learnable parameters – grows with the size of the graph.
This may lead to computational and memory limitations. As an engineering
workaround, one can approximate and simplify the convolution in (7) as
$\sum_{k\in{\mathcal{U}}}h_{k}{\mathbf{T}}_{k}$, where
${\mathcal{U}}\subset{\mathcal{V}}$. In principle, this comes at the expense
of reducing the expressiveness of the DCN. Still, this limitation can be
alleviated by concatenating several layers since the cross-product of
different ${\mathbf{T}}_{k}$ and ${\mathbf{T}}_{l}$ can give rise to a new
GSO. As we illustrate in the upcoming numerical evaluation, in practice one
can randomly select a small number of matrices ${\mathbf{T}}_{k}$ to define
the convolution and maintain a competitive performance. Developing a
principled approach to determine the subset of nodes ${\mathcal{U}}$
constitutes an interesting problem that is considered a future research
direction.

| Network Diffusion | Source Identification
---|---|---
| MNSE | Time (s) | Accuracy | Time (s)
DCN | $\mathbf{0.016\pm 0.014}$ | 3.6 | $0.052\pm 0.014$ | 7.5
DCN-30 | $\mathbf{0.029\pm 0.017}$ | 3.5 | $0.052\pm 0.016$ | 7.4
DCN-10 | $0.058\pm 0.021$ | 3.5 | $0.055\pm 0.015$ | 7.2
DCN-T | $0.098\pm 0.024$ | 4.1 | $\mathbf{0.991\pm 0.018}$ | 8.2
DCN-30-T | $0.199\pm 0.030$ | 3.7 | $\mathbf{0.983\pm 0.032}$ | 7.64
DCN-10-T | $0.229\pm 0.030$ | 3.5 | $0.865\pm 0.141$ | 7.38
LS | $0.050\pm 0.022$ | 0.4 | $0.05\pm 0.016$ | 0.36
FB-GCNN | $0.091\pm 0.028$ | 3.4 | $0.739\pm 0.172$ | 7.4
GCN | $0.167\pm 0.037$ | 3.3 | $0.155\pm 0.216$ | 7.1
GAT | $0.649\pm 0.089$ | 13.8 | $0.044\pm 0.081$ | 28.4
GraphSAGE | $0.359\pm 0.039$ | 5.9 | $0.676\pm 0.163$ | 12.5
GIN | $0.402\pm 0.079$ | 6.0 | $0.19\pm 0.163$ | 12.5
MLP | $0.353\pm 0.039$ | 2.2 | $0.050\pm 0.016$ | 4.7
TABLE I: Normalized MSE and accuracy when respectively solving the network
diffusion learning and source identification tasks. We report the mean
performance and standard deviation across 25 realizations, along with the
training time.

##  IV Numerical Evaluation

We assess the performance of DCN in different settings and compare it with
several baselines. Code to reproduce all results is available on
GitHub111https://github.com/reysam93/dag_conv_nn, and the interested reader is
referred there for additional experiments and implementation details.

$0$$0.1$$0.2$$0.3$$0.4$$0.5$$10^{-2}$$10^{-1}$$10^{0}$(a) Normalized noise
power.NMSEDCNLSFBGNN-4GCNNGraphSAGEGIN

0.10.30.50.70.9$0$$0.2$$0.4$$0.6$$0.8$$1$(b) Proportion of unobserved
nodes.AccuracyDCNDCN-30DCN-20DCN-10GraphSAGEGIN

$0.1$$0.3$$0.5$$0.7$$0$$0.2$$0.4$$0.6$$0.8$$1$(c) Edge
ProbabilityAccuracyDCNFBGNN-5GraphSAGEDCN-20FBGNN-2GIN

Figure 1: (a) reports the MNSE in the network diffusion task as the noise in
the observations increases; For the source identification task, (b) and (c)
depict the influence of increasing respectively the proportion of unobserved
nodes and the edge probability. We report the median performance across 25
realizations and values between the first and third quartile in the shaded
area.

Experiment setup. The numerical evaluation is conducted over synthetic data.
Unless specified otherwise, graphs are sampled from an Erdős-Rényi random DAG
model with $N=100$ nodes and an edge probability $p=0.2$. The input-output
signals are related as ${\mathbf{y}}={\mathbf{H}}{\mathbf{x}}$. Here,
${\mathbf{x}}$ is a sparse input signal whose non-zero entries are restricted
to the first nodes, ${\mathbf{H}}$ is a DAG filter as in (5) composed of 25
causal GSOs selected uniformly at random, and ${\mathbf{y}}$ is the output of
the diffusion of ${\mathbf{x}}$ across the DAG. We create $M=2000$ input-
output pairs of DAG signals and split them as 70% for training, 20% for
validation, and 10% for testing.

We consider two learning problems. First, we address the task of learning a
_network diffusion_ process where, given a new sparse input
${\mathbf{x}}_{test}$, we aim at estimating the associated output
${\mathbf{y}}_{test}$. The observed input and output signals are corrupted
with zero-mean additive white Gaussian noise with a normalized power of
$0.05$. This task amounts to solving a regression problem where we consider
the MSE as the loss function ${\mathcal{L}}$, and report the normalized MSE
(NMSE), computed as
$\frac{1}{T}\sum_{t=1}^{T}\|{\mathbf{y}}_{t}-{\hat{\mathbf{y}}}_{t}\|_{2}^{2}/\|{\mathbf{y}}_{t}\|_{2}^{2}$
over $T$ test signals. The second task tackles the _source identification_
problem, where we are given a partially observed output signal
${\mathbf{y}}_{test}$ and the goal is to predict which node was the source
originating the observation. Now, the input ${\mathbf{x}}$ contains a single
non-zero entry constrained to be located within the first 20 nodes. Since the
diffusion ${\mathbf{y}}={\mathbf{H}}{\mathbf{x}}$ does not alter the signal
value at the source nodes, we mask the output ${\mathbf{y}}$ so observations
from possibly seeding nodes are not available. The task is cast as a
classification problem where the graph label is the index of the source node.
We adopt the cross-entropy loss and report the mean accuracy. The reported
results comprise 25 independent realizations.

Baselines. The performance of the DCN is compared with several baselines. The
least-squares (LS) model captures the true linear model, computes an LS
estimate of ${\mathbf{h}}$ from the training data and then predicts test
outputs. As non-linear baselines, we include FB-CGNN, the convolutional GNN
based on graph filters from [9] [cf. (1)], GCN [7], GAT [10], GIN [8],
GraphSAGE [32], and the MLP. We do not include D-VAE [25] and DAGNN [26] here
(although they are GNNs for DAG-based data), because their training time is
prohibitive for moderately large graphs (as here where $N=100$).

Preliminary results. We start by assessing the overall performance of the
models considered in the network diffusion and the source identification
problems. The results collected in Table I showcase the superior performance
of the proposed DCN in both tasks. In the network diffusion setting, we
observe how using a subset of 30 or 10 GSOs randomly selected (DCN-30 and
DCN-10) decouples the DCN from the number of nodes without significantly
hindering the performance; DCN-30 is the second-best alternative. Comparing
DCN and LS, the superior performance of DCN is due to a major resilience to
the noise in the inputs, as will be evident in the upcoming experiments. As
expected, when using the transposed matrices ${\mathbf{T}}_{k}^{\top}$ in
DCN-T (i.e., when following the reverse ordering of the DAG), the performance
drops significantly, especially when using only a subset of the GSOs. In
contrast, moving on to the source identification task we observe that while
the DCN variants are incapable of identifying the source node, the
alternatives DCN-T have an almost perfect accuracy. Intuitively, this reflects
how for identifying the source nodes from the output one needs to navigate the
DAGs in the reverse order as facilitated by the transposed GSOs. In
conclusion, this highlights the importance of properly harnessing the
directionality of the DAG. Note that for the baselines we use the GSO with the
best performance.

Test case 1. The results depicted in Fig. 1a reflect how increasing the power
of the noise in the observed signals impacts the performance of the network
diffusion problem. In the absence of noise, the LS model outperforms the
alternatives, but its performance deteriorates swiftly as the noise power
increases. In contrast, the proposed DCN attains a similar NMSE to that of the
LS in the noiseless setting and its performance is more robust and in the
presence of noise. This experiment showcases that i) the proposed DCN obtains
a performance comparable to that of the optimal solution (LS); and ii)
assuming a non-linear model even though the true generative model is linear
provides a solution more robust to noise.

Test Case 2. Moving on to the source identification problem, Fig. 1b evaluates
the accuracy of the DCN when a different number of GSOs is employed as the
proportion of unobserved nodes increases. Based on the results of Table I, all
the methods use the transposed of their respective GSO. It is worth mentioning
that higher values on the x-axis render a more challenging setting since: i)
fewer nodal observations are available; and ii) every unobserved node is
allowed to be a source node, hence increasing the complexity of the
classification problem. First, we note that the accuracy of DCN remains high
even for large fractions of unobserved nodes, demonstrating its robustness.
From the performance of DCN-30, DCN-20, and DCN-10, the trade-off between
reducing the number of GSOs and maintaining a high accuracy and low
variability is evident. Nonetheless, DCN-30 shows a similar performance to
that of DCN up to observing 50% of the nodes, and DCN-10 outperforms the non-
DAG-based alternatives.

Test Case 3. Finally, we evaluate the impact of the density of the graph on
the source identification problem. To that end, we measure the accuracy as the
edge probability increases (Fig. 1c), where we observe that a denser graph
poses a greater challenge to identify the source. Intuitively, as the graph
gets closer to a fully connected DAG, the signal generated from different
nodes becomes more homogeneous, with the discrepancies being due only to the
distance between seeding nodes. In addition, this experiment shows that both
DCN and DCN-20 consistently outperform FBGNN using filters of order 2 and 5,
emphasizing the benefits stemming from replacing GSP filters with DAG filters.

##  V Conclusions

In this work, we proposed DCN, a GCNN tailored to the particular structure
encoded by DAGs. At its core lies a novel convolutional layer, which is built
upon recent results to process signals defined over DAGs and posets, has ties
to linear SEMs, is permutation equivariant, endows the architecture with a
spectral representation, and exploits the partial order imposed by the DAG.
Moreover, our preliminary numerical evaluation demonstrated a promising
performance and showed that the convolution operation can be detached from the
number of nodes, thus addressing a potential limitation.

## References

  * [1] A. Ortega, P. Frossard, J. Kovačević, J. M. F. Moura, and P. Vandergheynst,  “Graph signal processing: Overview, challenges, and applications,”  Proc. IEEE, vol. 106, no. 5, pp. 808–828, 2018.
  * [2] X. Dong, D. Thanou, L. Toni, M. Bronstein, and P. Frossard,  “Graph signal processing for machine learning: A review and new perspectives,”  IEEE Signal Process. Mag., vol. 37, no. 6, pp. 117–127, 2020\.
  * [3] M. M. Bronstein, J. Bruna, Y. LeCun, A. Szlam, and P. Vandergheynst,  “Geometric deep learning: Going beyond euclidean data,”  IEEE Signal Process. Mag., vol. 34, no. 4, pp. 18–42, July 2017\.
  * [4] S. Rey, A. G. Marques, and S. Segarra,  “An underparametrized decoder architecture for graph signals,”  in IEEE Intl. Wrksp. Computat. Advances Multi-Sensor Adaptive Process. (CAMSAP). IEEE, 2019, pp. 231–235.
  * [5] Z. Wu, S. Pan, F. Chen, G. Long, C. Zhang, and S. Y. Philip,  “A comprehensive survey on graph neural networks,”  IEEE Trans. Neural Netw. Learn. Syst., 2020.
  * [6] V. M. Tenorio, S. Rey, F. Gama, S. Segarra, and A. G. Marques,  “A robust alternative for graph convolutional neural networks via graph neighborhood filters,”  in Conf. Signals, Syst., Computers (Asilomar). IEEE, 2021, pp. 1573–1578.
  * [7] T. N. Kipf and M. Welling,  “Semi-supervised classification with graph convolutional networks,”  arXiv preprint arXiv:1609.02907, 2016.
  * [8] K. Xu, W. Hu, J. Leskovec, and S. Jegelka,  “How powerful are graph neural networks?,”  arXiv preprint arXiv:1810.00826, 2018.
  * [9] L. Ruiz, F. Gama, and A. Ribeiro,  “Graph neural networks: Architectures, stability, and transferability,”  Proc. IEEE, vol. 109, no. 5, pp. 660–682, 2021.
  * [10] P. Veličković, G. Cucurull, A. Casanova, A. Romero, P. Liò, and Y. Bengio,  “Graph attention networks,”  in Intl. Conf. Learn. Repr. (ICLR), 2018.
  * [11] C. Wang, S. Pan, G. Long, X. Zhu, and J. Jiang,  “MGAE: Marginalized graph autoencoder for graph clustering,”  in Assoc. Comput. Mach., 2017, pp. 889–898.
  * [12] S. Rey, V. M. Tenorio, S. Rozada, L. Martino, and A. G. Marques,  “Overparametrized deep encoder-decoder schemes for inputs and outputs defined over graphs,”  in European Signal Process. Conf. (EUSIPCO). IEEE, 2021, pp. 855–859.
  * [13] R. Shafipour, A. Khodabakhsh, G. Mateos, and E. Nikolova,  “Digraph Fourier transform via spectral dispersion minimization,”  in IEEE Int. Conf. Acoustics, Speech Signal Process. (ICASSP), 2018, pp. 6284–6288.
  * [14] R. Shafipour, S. Segarra, A. G. Marques, and G. Mateos,  “Directed network topology inference via graph filter identification,”  in IEEE Data Science Wrksp. (DSW), 2018, pp. 210–214.
  * [15] A. G. Marques, S. Segarra, and G. Mateos,  “Signal processing on directed graphs: The role of edge directionality when processing and learning from network data,”  IEEE Signal Process. Mag., vol. 37, no. 6, pp. 99–116, 2020.
  * [16] B. Seifert, C. Wendler, and M. Püschel,  “Causal Fourier analysis on directed acyclic graphs and posets,”  IEEE Trans. Signal Process., vol. 71, pp. 3805–3820, 2023.
  * [17] J. Peters, D. Janzing, and B. Schölkopf,  Elements of Causal Inference: Foundations and Learning Algorithms,  The MIT Press, 2017.
  * [18] J. Peters and P. Bühlmann,  “Identifiability of Gaussian structural equation models with equal error variances,”  Biometrika, vol. 101, no. 1, pp. 219–228, 2014.
  * [19] X. Zheng, B. Aragam, P. K. Ravikumar, and E. P. Xing,  “DAGs with no tears: Continuous optimization for structure learning,”  Advances Neural Inf. Process. Syst., vol. 31, 2018.
  * [20] D’Acunto G., Di Lorenzo P., and Barbarossa S.,  “Multiscale causal structure learning,”  Trans. Mach. Learning Res., 2023.
  * [21] S. S. Saboksayr, G. Mateos, and M. Tepper,  “CoLiDE: Concomitant linear DAG estimation,”  in Intl. Conf. Learn. Repr. (ICLR), 2024.
  * [22] M. Allamanis, E. T. Barr, P. Devanbu, and C. Sutton,  “A survey of machine learning for big code and naturalness,”  ACM Computing Surveys (CSUR), vol. 51, no. 4, pp. 1–37, 2018.
  * [23] C. Zhang, M. Ren, and R. Urtasun,  “Graph hypernetworks for neural architecture search,”  arXiv preprint arXiv:1810.05749, 2018.
  * [24] M. Püschel, B. Seifert, and C. Wendler,  “Discrete signal processing on meet/join lattices,”  IEEE Trans. Signal Process., vol. 69, pp. 3571–3584, 2021.
  * [25] M. Zhang, S. Jiang, Z. Cui, R. Garnett, and Y. Chen,  “D-VAE: A variational autoencoder for directed acyclic graphs,”  in Advances Neural Inf. Process. Syst., H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, Eds. 2019, vol. 32, Curran Associates, Inc.
  * [26] V. Thost and J. Chen,  “Directed acyclic graph neural networks,”  in Intl. Conf. Learn. Repr. (ICLR), 2021.
  * [27] Y. Luo, V. Thost, and L. Shi,  “Transformers over directed acyclic graphs,”  in Advances Neural Inf. Process. Syst., A. Oh, T. Neumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine, Eds. 2023, vol. 36, pp. 47764–47782, Curran Associates, Inc.
  * [28] S. Segarra, A. G. Marques, and A. Ribeiro,  “Optimal graph-filter design and applications to distributed linear network operators,”  IEEE Trans. Signal Process., vol. 65, no. 15, pp. 4117–4131, 2017\.
  * [29] Daniel J Lehmann,  “Algebraic structures for transitive closure,”  Theoretical Comput. Sci., vol. 4, no. 1, pp. 59–76, 1977.
  * [30] G.-C. Rota,  “On the foundations of combinatorial theory: I. Theory of Möbius functions,”  in Classic Papers Combinatorics, pp. 332–360. Springer, 1964.
  * [31] S. Rey, S. Segarra, R. Heckel, and A. G. Marques,  “Untrained graph neural networks for denoising,”  IEEE Trans. Signal Process., vol. 70, pp. 5708–5723, 2022.
  * [32] W. Hamilton, Z. Ying, and J. Leskovec,  “Inductive representation learning on large graphs,”  Advances Neural Inf. Process. Syst., vol. 30, 2017.