# Bayesian Geographically Weighted Regression
using Fused Lasso Prior

Toshiki Sakai
Graduate School of Culture and Information Science, Doshisha University
<EMAIL_ADDRESS>
&Jun Tsuchida
Faculty of Data Science, Kyoto Women’s University
Hiroshi Yadohisa
Faculty of Culture and Information Science, Doshisha University

###### Abstract

A main purpose of spatial data analysis is to predict the objective variable
for the unobserved locations. Although Geographically Weighted Regression
(GWR) is often used for this purpose, estimation instability proves to be an
issue. To address this issue, Bayesian Geographically Weighted Regression
(BGWR) has been proposed. In BGWR, by setting the same prior distribution for
all locations, the coefficients’ estimation stability is improved. However,
when observation locations’ density is spatially different, these methods do
not sufficiently consider the similarity of coefficients among locations.
Moreover, the prediction accuracy of these methods becomes worse. To solve
these issues, we propose Bayesian Geographically Weighted Sparse Regression
(BGWSR) that uses Bayesian Fused Lasso for the prior distribution of the BGWR
coefficients. Constraining the parameters to have the same values at adjacent
locations is expected to improve the prediction accuracy at locations with a
low number of adjacent locations. Furthermore, from the predictive
distribution, it is also possible to evaluate the uncertainty of the predicted
value of the objective variable. By examining numerical studies, we confirmed
that BGWSR has better prediction performance than the existing methods (GWR
and BGWR) when the density of observation locations is spatial difference.
Finally, the BGWSR is applied to land price data in Tokyo. Thus, the results
suggest that BGWSR has better prediction performance and smaller uncertainty
than existing methods.

_K_ eywords Hierarchical Bayesian Model $\cdot$ Laplace Prior $\cdot$ Spatial
Heterogeneity

##  1 Introduction

Spatial data denote data that comprise individual values and the spatial
location information where they were observed. The application of spatial
data, such as mapping to understand visual features and creating hazard maps
for disasters (Fischer and Getis, 2010), is used in various research fields.
When we obtain spatial data, we sometimes set objective variables and
covariates for the variables obtained. The purposes of analyzing such spatial
data include visualization of spatial data, classification and division of
space, and prediction of the value of the objective variable (e.g., Jena et
al., 2023; Wang et al., 2023). The prediction of the objective variable’s
value is important in analyzing spatial data because there are many unobserved
locations for the objective variables in many cases.

There are two main ways to predict the objective variable. The first applies
to the spatial data containing only the objective variable. The methods used
in this case include Kriging (Krige, 1951), inverse distance weighting (Bailey
et al., 1995), and spline (Cline, 1973). For example, Kriging predicts the
value of the objective variable at the prediction location using a weighted
average of the available objective variable data. The second applies to the
spatial data containing covariates in addition to the objective variable. The
methods used in this case include Kriging with an external drift (Ahmed and De
Marsily, 1987). One of the most frequently used methods is Geographically
Weighted Regression (GWR) (Brunsdon et al., 1996; Fotheringham et al., 2003).
GWR is based on linear regression and assumes a multiple regression model with
different coefficients for each location. Further, it assumes that the values
of the coefficients are similar for locations close to each other—the method
assumes spatial autocorrelation in the coefficients. In Yu et al. (2007), GWR
is used to investigate the relationship among housing prices, floor space, and
age of the house in Milwaukee, Wisconsin, in the United States. In Koh et al.
(2020), GWR is used to investigate the relationship among Nitrate
concentrations, land use, and precipitation in South Korea.

Notably, when applying GWR to real data, the estimation of coefficients
becomes unstable when the number of observation locations is small (LeSage,
2004; Subedi et al., 2018) because of the small number of observations used to
estimate the coefficients. In extreme cases, for estimation coefficients, we
could use only one or two locations, and the solution of the weighted least
squared method may not be uniquely determined. To solve this issue, Bayesian
Geographically Weighted Regression (BGWR) has been proposed (LeSage, 2004).
BGWR is a method for estimating parameters of GWR within the framework of
Bayesian statistics. Specifically, parameters such as coefficients are assumed
to follow a certain probability distribution, and the estimator is the
expected value of the posterior distribution calculated by using the data.
Setting a common prior distribution for all locations allows for a stable
estimation of coefficients when the number of locations is small, as
information from the prior distribution becomes more weighted in the
estimation. Additionally, BGWR allows us to evaluate the instability of the
objective variable’s predicted value because the predicted value’s credibility
is calculated by the predictive distribution.

However, these methods sometimes estimate considerably different coefficients
for adjacent locations. For example, two close yet adjacent locations do not
share the same location for estimation. In such a situation, the similarity of
the coefficients at adjacent locations, as assumed by spatial autocorrelation,
cannot be fully captured because the coefficients are estimated using
completely different data (Fotheringham et al., 2017).

To resolve this issue, we propose a novel method by using the idea of Fused
Lasso (Tibshirani et al., 2005). Fused Lasso penalizes the absolute value of
the difference between the coefficients at adjacent locations, making it
easier to obtain the same estimated value of the coefficient at adjacent
locations. For example, Li and Sang (2019) and Zhong et al. (2023) use Fused
Lasso penalties to consider the similarity of the coefficient of regression
and show improved prediction performance compared to the GWR. Hence, the
proposed method—Bayesian Geographically Weighted Sparse Regression
(BGWSR)—combines BGWR and Bayesian Fused Lasso Kyung et al. (2010) and is
expected to obtain similar estimates of coefficients between adjacent
locations. This is because the Laplace distribution is set for the difference
in coefficients between adjacent locations, and the difference between
adjacent locations is easily estimated as $0$. Moreover, as BGWSR can consider
the similarity between adjacent locations, it is expected to improve the
prediction accuracy of coefficients and objective variables compared to
existing methods.

In BGWR and GWR, when observation locations are obtained from a spatially
heterogeneous environment (i.e., there are areas whose density is sparse and
dense), the estimation of coefficients and predictions for the observed
locations may be less accurate. For example, if adjacent locations are defined
as extremely close in distance, the number of observation locations used to
estimate the coefficients will not be sufficient in areas where the density of
observation locations is sparse, and the resulting estimation will be
unstable. Therefore, for stable estimation, it is necessary to have excessive
information in the prior distribution. Nevertheless, when adjacent locations
are defined as distant locations, the estimated values of the coefficients may
be the same for most locations because the same data are shared by many
locations and the penalty of Fused Lasso. Furthermore, in areas with dense
observation locations, it is necessary to define adjacent locations as those
close to each other to consider rapid changes adequately. By introducing a
prior distribution for setting the width of each location to be an adjacent
location, it is expected that the appropriate adjacent location will be
determined for all locations. This allows for a trade-off between the breadth
of the definition of the adjacency points and the strength of the to-be-
determined Fused Lasso penalty so that the estimation performance is better
for each location. Therefore, introducing this prior is expected to improve
the prediction accuracy of coefficients and objective variables even when
locations are spatially heterogeneous.

This paper is organized as follows. In Section 2, we introduce the proposed
method and other predictive algorithms. In Section 3, we evaluate the
numerical performance of the proposed methods together with some existing
methods through numerical studies. In Section 4, we demonstrate the proposed
method through spatial regression modeling of the land prices in the Tokyo
metropolitan area. Finally, Section 5 summarizes our main conclusions and
identifies areas where further research is necessary.

##  2 Bayesian Geographically Weighted Sparse Regression (BGWSR)

###  2.1 Model

This section describes the model for Bayesian estimation of the BGWSR. Let
$D\subset\mathbb{R}^{2}$ be the region to be analyzed, and let
$\bm{s}_{i}=(s_{1i},\,s_{2i})^{\prime}\in D$ be the $i$-th location. Suppose a
set of spatial data ${(x(\bm{s}_{i}),\,y(\bm{s}_{i}))}$ is observed locations
$\bm{s}_{1},\,\ldots,\,\bm{s}_{n}\in D$, where the objective variable
$y(\bm{s}_{i})$ is assumed to be spatially correlated, and
$\bm{x}(\bm{s}_{i})=(x_{1}(\bm{s}_{i}),\,\ldots,\,x_{p}(\bm{s}_{i}))^{\prime}$
is the $p$-dimensional vector of covariates at the location $\bm{s}_{i}$.

$\displaystyle
y(\bm{s}_{i})=\bm{x}(\bm{s}_{i})^{\prime}\bm{\beta}(\bm{s}_{i})+\varepsilon(\bm{s}_{i}).$

We assume that $\bm{\beta}(\bm{s}_{i})$ is the coefficient and that the error
terms $\varepsilon(\bm{s}_{i})$ follow a Normal distribution
${\rm{N}}(0,\,\sigma^{2})(\sigma^{2}>0)$ independent of each other.

To estimate the coefficients, we consider that the coefficients are described
by $n$-dimensional vectors. Let $\bm{\beta}_{k}$ be a vector of coefficients
whose $i$-th elements are the coefficient
$\beta_{k}(\bm{s}_{i})(k=1,\,2,\,\ldots,\,p)$, and be defined as
$\bm{y}_{k}=\bm{X}_{k}\bm{\beta}_{k}+\bm{\varepsilon}_{k}$, where
$\bm{\varepsilon}_{k}\sim{\rm{MN}}(\bm{0}_{n},\,\sigma^{2}\bm{I}_{n}/p)$, and
$\bm{I}_{n}$is $n$-dimensional identity matrix.
${\rm{MN}}(\bm{\mu},\,\bm{\Sigma})$ denotes the Multivariate Normal
distribution of the mean vector $\bm{\mu}$ and covariance matrix
$\bm{\Sigma}$. The $\bm{X}_{k}$ is an $n$-dimensional diagonal matrix and
$\bm{X}_{k}={\rm{diag}}\\{x_{k}(\bm{s}_{1}),\,x_{k}(\bm{s}_{2}),\,\ldots,\,x_{k}(\bm{s}_{n}))\\}\in\mathbb{R}^{n\times
n}$. Using these terms, the $\bm{y}=(y(\bm{s}_{i}))$ is represented as
$\bm{y}=\sum_{k=1}^{p}\bm{y}_{k}=\sum_{k=1}^{p}\bm{X}_{k}\bm{\beta}_{k}+\bm{\varepsilon}_{k}$.
The model for estimating the parameters at location $\bm{s}_{i}$ is given as
follows based on BGWR (LeSage, 2004).

$\displaystyle\bm{W}(\bm{s}_{i})\bm{y}_{k}=\bm{W}(\bm{s}_{i})\bm{X}_{k}\bm{\beta}_{k}+\bm{\varepsilon}_{k}.$
(1)

$\bm{W}(\bm{s}_{i})$ is an $n$-dimensional diagonal matrix, and
$\bm{W}(\bm{s}_{i})={\rm{diag}}\\{w_{1}(\bm{s}_{i}),\,w_{2}(\bm{s}_{i}),\,\ldots,\,w_{n}(\bm{s}_{i})\\}$
and $w_{j}(\bm{s}_{i})\,(\geq 0;j=1,\,2,\,\ldots,\,n)$ are weights determined
by the distance between two locations $\bm{s}_{i}$ and $\bm{s}_{j}$.

For handling the spatial autocorrelation of coefficients, the Fused Lasso
Prior is selected as the prior distribution of the coefficients in BGWSR. The
following distribution is set as the prior distribution of the coefficient
$\bm{\beta}_{k}\,(k=1,\,2,\,\ldots,\,p)$ for the $k$-th covariates.

$\displaystyle\begin{split}\bm{\beta}_{k}|\sigma^{2},\,S,\,N,\,C,\,\lambda_{k,\,1},\,\lambda_{k,\,2}\sim(\sigma^{2})^{-\frac{2n+\sum_{i=1}^{n}n_{i}}{4}}&\prod_{i=1}^{n}{\rm{Laplace}}\left(\frac{\beta_{k}(\bm{s}_{i})}{\sqrt{\sigma^{2}}}\Big{|}0,\lambda_{k,\,1}\right)\\\
&\quad\times\prod_{(i,j)\in\bm{C}}{\rm{Laplace}}\left(\frac{\beta_{k}(\bm{s}_{i})-\beta_{k}(\bm{s}_{j})}{\sqrt{\sigma^{2}}}\Big{|}0,\,\lambda_{k,\,2}\right),\end{split}$
(2)

where ${\rm{Laplace}}(x|\mu,\,b)$ denotes the Laplace distribution of mean
$\mu$ and scale parameter $b$ with $x$ as a random variable. $S$ is the set of
the observed locations and
$S=\\{\bm{s}_{1},\,\bm{s}_{2},\,\ldots,\,\bm{s}_{n}\\}\subset D$. $N$ is the
set of the number of adjacent locations $n_{i}$ at location $\bm{s}_{i}$,
where $N=\\{n_{1},\,n_{2},\,\ldots,\,n_{n}\\}$. $C$ is a set of pairs of
adjacent locations, and $(i,\,j)\in C$ indicates that two locations
$\bm{s}_{i},\,\bm{s}_{j}$ are spatially close to each other.
$\lambda_{k,1},\,\lambda_{k,2}$ are parameters that determine the strength of
the shrinkage. For example, the following is obtained.

$\displaystyle C=\\{(i,j)\mid w_{j}(\bm{s}_{i})>0\\}.$

The first term on the right-hand side of the equation (2) is a Laplace
distribution with mean $0$ for the coefficient at each location based on the
idea of Bayesian Lasso (Park and Casella, 2008). The role of this term is to
shrink the coefficient to $0$ regardless of location. As it tends to shrink to
$0$ of the estimated value of the coefficient at a location with a small size
adjacent set to $0$, the estimation is more numerically stable than GWR. The
second term on the right side of the equation (2) is based on the idea of
Bayesian Fused Lasso (Kyung et al., 2010), setting a Laplace distribution with
mean $0$ for the difference of coefficients at adjacent locations. This allows
us to estimate similar coefficients for adjacent locations.

As the coefficients $\bm{\beta}_{k}$ have a Laplace distribution for the prior
distribution, the full conditional distribution for Gibbs sampling cannot be
derived explicitly without modification. We apply parameter augmentation
(Andrews and Mallows, 1974). For calculating the posterior distribution, the
prior distribution of the coefficient $\bm{\beta}_{k}$ can be rewritten using
the parameter
$\tau_{k,\,i}^{2},\,\omega_{k,\,i,\,j}^{2}\,(i=1,\,2,\,\ldots,\,n;j=1,\,2,\,\ldots,\,n,\,i\neq
j;k=1,\,2,\,\ldots,\,p)$ as follows, where
$\tau_{k,i}^{2}>0,\,\omega_{k,\,i,\,j}^{2}>0$.

$\displaystyle\bm{\beta}_{k}|\sigma^{2},T_{k},\,\Omega_{k}$
$\displaystyle\sim{\rm{MN}}(\bm{0}_{n},\,\sigma^{2}\bm{\Sigma}_{k}^{-1}),\,$
$\displaystyle T_{k}|\lambda_{k,\,1}^{2}$
$\displaystyle\sim\prod_{i=1}^{n}\frac{\lambda_{k,\,1}^{2}}{2}{\rm{exp}}\left\\{-\frac{\lambda_{k,\,1}}{2}\tau_{k,i}^{2}\right\\},\,$
$\displaystyle\Omega_{k}|\lambda_{k,\,2}^{2}$
$\displaystyle\sim\prod_{(i,\,j)\in\bm{C}}\frac{\lambda_{k,\,2}^{2}}{2}{\rm{exp}}\left\\{-\frac{\lambda_{k,\,2}}{2}\omega_{k,i,j}^{2}\right\\}.$

Here,
$T_{k}=\\{\tau_{k,1}^{2},\,\tau_{k,2}^{2},\,\ldots,\tau_{k,n}^{2}\\},\,\Omega_{k}=\\{\omega_{k,1,2}^{2},\,\omega_{k,1,3}^{2},\,\ldots,\,\omega_{k,n-1,n}^{2}\\}$
and
$\omega_{k,\,i,\,j}^{2}=\omega_{k,\,j,\,i}^{2}\quad(k=1,\,2,\,\ldots,\,p)$.
Let $\sigma^{2}\bm{\Sigma}_{k}^{-1}$ be the covariance matrix of the
coefficients defined as

$\displaystyle(\bm{\Sigma}_{k}^{-1})_{ij}=\left\\{\begin{array}[]{ll}\frac{1}{\tau_{k,i}^{2}}+\sum_{(i,\ell)\in
C_{i}}\frac{1}{\omega_{k,\,i,\ell}^{2}}&(i=j)\\\
-\frac{1}{\omega_{k,i,j}^{2}}&(i\neq j,\,(i,\,j)\in C)\\\
0&({\rm{otherwise}}),\end{array}\right.$

where $C_{i}\subset C$ is a set comprising only pairs of locations adjacent to
$\bm{s}_{i}$.

Based on Kyung et al. (2010), the models for the other parameters are set as

$\displaystyle\lambda_{k,\,1}^{2}|r_{1},\,q_{1}\sim{\rm{Gamma}}(r_{1},\,q_{1}),\quad\lambda_{k,\,2}^{2}|r_{2},\,q_{2}\sim{\rm{Gamma}}(r_{2},\,q_{2}),$
$\displaystyle\sigma^{2}|r,\,q\sim{\rm{IGamma}}\left(\frac{r}{2},\,\frac{q}{2}\right),\quad
w_{j}(\bm{s}_{i})|h=f(d_{ij}|h),$

where ${\rm{Gamma}}(a,\,b),\,{\rm{IGamma}}(a,\,b)$ denote the gamma
distribution and Inverse Gamma distribution with shape parameter $a$ and
reciprocal scale parameter $b$. Meanwhile
$r,\,q,\,r_{1},\,q_{1},\,r_{2},\,q_{2}$ are hyperparameters determined by the
analyst. The function $f(d_{ij}|h)$ determines the weights of two locations
$\bm{s}_{i},\,\bm{s}_{j}$ with the distance $d_{ij}$ as an argument. For
example, the following equation is used.

$\displaystyle
f(d_{ij}|h)=\begin{cases}1-\left(\frac{d_{ij}}{h}\right)^{2}&\,\left(\left\|\boldsymbol{s}_{i}-\boldsymbol{s}_{j}\right\|<h\right)\\\
0&\text{ (otherwise) }\end{cases}.$

$h$ is a parameter that determines the degree of spatial correlation. It may
be determined by the analyst based on domain knowledge, or it can be
incorporated as a model. For example, based on Boehm Vock et al. (2015); Ma et
al. (2021), we have

$\displaystyle h|a\sim{\rm{U}}(0,\,a),$ (3)

where ${\rm{U}}(0,\,a)$ represents a Uniform distribution with lower limit $0$
and upper limit $a$, where $a$ is the hyperparameter determined by the
analyst. For example, the maximum distance between locations in a region is
set as $a$. Sampling methods of the posterior distribution of $h$ include the
Metropolis–Hastings Algorithm (MH method) (Ma et al., 2021). The method of
determining the optimal adjacency for each location is described in 2.3.

###  2.2 Parameter estimation

The probability density function of the simultaneous posterior distributions
of the parameters for BGWSR has a complex form. Therefore, Gibbs sampling and
MH method are used to estimate the posterior distribution of the parameters.
For $\bm{\beta}_{k},\,\sigma^{2},\,\lambda_{k,\,1}^{2},\,\lambda_{k,\,2}^{2}$,
Gibbs sampling is used because the probability density function of the
conditional distribution can be derived explicitly. For the parameter $h$, we
use the MH method based on Ma et al. (2021).

First, We show the conditional distribution for Gibbs sampling. The full
conditional distribution of the coefficients $\bm{\beta}_{k}$ is given as
follows,

$\displaystyle\begin{split}\bm{\beta}_{k}|&\bm{X}_{k},\,\bm{y}_{k},\,\bm{\Sigma}_{k}^{-1},\,W,\,\sigma^{2},\,T_{k},\,\Omega_{k}\\\
&\sim{\rm{MN}}\left(\left(\bm{X}_{k}\sum_{i=1}^{n}(\bm{W}(\bm{s}_{i}))^{2}\bm{X}_{k}+\bm{\Sigma}_{k}^{-1}\right)^{-1}\bm{X}_{k}\sum_{i=1}^{n}(\bm{W}(\bm{s}_{i}))^{2}\bm{y}_{k},\,\sigma^{2}\left(\bm{X}_{k}\sum_{i=1}^{n}\bm{W}(\bm{s}_{i})^{2}\bm{X}_{k}+\bm{\Sigma}_{k}^{-1}\right)\right),\end{split}$

where $W$ denotes the set of weight matrices for all locations, and
$W=\\{\bm{W}(\bm{s}_{1}),\,\bm{W}(\bm{s}_{2}),\,\ldots,\,\bm{W}(\bm{s}_{n})\\}$.
The conditional distribution of the parameters $T_{k},\,\Omega_{k}$ added by
parameter augmentation is as follows.

$\displaystyle\frac{1}{\tau_{k,i}^{2}}\Big{|}\beta_{k}(\bm{s}_{i}),\,\sigma^{2},\,\lambda_{k,\,1}$
$\displaystyle\sim{\rm{IGauss}}\left(\sqrt{\frac{\lambda_{k,\,1}^{2}\sigma^{2}}{\beta_{k}(\bm{s}_{i})^{2}}},\,\lambda_{k,\,1}^{2}\right),$
$\displaystyle\frac{1}{\omega_{k,i,j}^{2}}\Big{|}\beta_{k}(\bm{s}_{i}),\,\beta_{k}(\bm{s}_{j}),\,\sigma^{2},\,\lambda_{k,\,2}$
$\displaystyle\sim{\rm{IGauss}}\left(\sqrt{\frac{\lambda_{k,\,2}^{2}\sigma^{2}}{(\beta_{k}(\bm{s}_{i})-\beta_{k}(\bm{s}_{j}))^{2}}},\,\lambda_{k,\,2}^{2}\right),$

where ${\rm{IGauss}}(\mu,\,\lambda)$ denotes the Inverse Gaussian distribution
with mean $\mu$ and shape parameter $\lambda$.

The full conditional distribution of the parameters
$\sigma^{2},\,\lambda_{k,\,1},\,\lambda_{k,\,2}$ is given below. We discuss
the parameter $h$ in the next section.

$\displaystyle\sigma^{2}|\bm{W}(\bm{s}_{i})\bm{y}_{k},\,\bm{W}(\bm{s}_{i})\bm{X}_{k},\,\bm{\beta}_{k}\sim{\rm{IGamma}}\left(r_{*},\,q_{*}\right),$
$\displaystyle\quad\quad r_{*}=r+\frac{pn^{2}}{2},$ $\displaystyle\quad\quad
q_{*}=q+\frac{\sum_{i=1}^{n}\sum_{k=1}^{p}\|\bm{W}(\bm{s}_{i})\bm{y}_{k}-\bm{W}(\bm{s}_{i})\bm{X}_{k}\bm{\beta}_{k}\|^{2}}{2},$
$\displaystyle\lambda_{k,\,1}^{2}|\bm{y}_{k},\,\bm{X}_{k},\,\bm{\beta}_{k},\,\sigma^{2},\,T_{k},\,\Omega_{k}\sim{\rm{Gamma}}(r_{1}^{*},\,q_{1}^{*}),$
$\displaystyle\quad\quad r_{1}^{*}=r_{1}+n,$ $\displaystyle\quad\quad
q_{1}^{*}=q_{1}+\frac{1}{2}\sum_{i=1}^{n}\tau_{k,i}^{2},$
$\displaystyle\lambda_{k,\,2}^{2}|\bm{y}_{k},\,\bm{X}_{k},\,\bm{\beta}_{k},\,\sigma^{2},\,T_{k},\,\Omega_{k}\sim{\rm{Gamma}}(r_{2}^{*},\,q_{2}^{*}),$
$\displaystyle\quad\quad r_{2}^{*}=r_{2}+\frac{1}{2}\sum_{i=1}^{n}n_{i},$
$\displaystyle\quad\quad q_{2}^{*}=q_{2}+\frac{1}{2}\sum_{(i,\,j)\in
C}\omega_{k,\,i,\,j}^{2}.$

###  2.3 Adjacency estimation

In BGWSR, BGWR, and GWR, the parameter $h$, which determines the range of
locations to be adjacent, is the same for all locations. Therefore, these
methods implicitly assume that the observed locations are obtained uniformly
within the region. However, actual data may include dense areas with many
locations and sparse areas with few locations. In an area with sparse
observation locations, numerous adjacent locations are required to have
adequate data for the coefficient estimation. Nevertheless, in an area with
dense locations, it is necessary to reduce the range of neighboring locations
to capture local changes in the coefficients successfully. Therefore, we
consider the range of adjacent locations for each location. Specifically, the
MH method is used for sampling, allowing the parameters responsible for
determining the range of adjacency to vary from one location to another. The
sampling scheme of $h$ in Equation (3) is the special case of this sampling
scheme.

We set the prior distribution of the parameter $h(\bm{s}_{i})$ that determines
the adjacent location at location $\bm{s}_{i}$ as

$\displaystyle
h(\bm{s}_{i})|a\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}{\rm{U}}(0,\,a),$

where $h(\bm{s}_{i})$ is a parameter that determines the adjacency range at
location $\bm{s}_{i}$. By setting $h(\bm{s}_{i})$ as a different parameter for
each location, we can expect to estimate a smaller adjacency range in areas
with dense observation locations and large local changes in the coefficients.
Subsequently, in areas with sparse observation locations, a large adjacency
range is expected to be estimated.

We use the MH method to sample the parameter $h(\bm{s}_{i})$. The sampling
procedure for $h(\bm{s}_{i})$ at location $\bm{s}_{i}$ using the MH method is
as follows.

  1. 1.

We use the Normal distribution as the proposal distribution. A new
$h(\bm{s}_{i})$ candidate $h_{\rm{prop}}(\bm{s}_{i})$ is sampled using the
proposal distribution as follows.

$\displaystyle
h_{\rm{prop}}(\bm{s}_{i})|\sigma_{h}^{2}\sim{\rm{N}}(h_{\rm{cur}}(\bm{s}_{i}),\sigma_{h}^{2}).$

In the above, $h_{\rm{cur}}(\bm{s}_{i})$ is the sampled $h(\bm{s}_{i})$ before
the update, and $\sigma_{h}^{2}$ is the hyperparameter representing the
variance of the proposal distribution.

  2. 2.

Decide whether to accept $h_{\rm{prop}}(\bm{s}_{i})$ based on the posterior
probability. Specifically, calculate $\alpha$ as defined in the following
equation.

$\displaystyle\alpha$
$\displaystyle=\min\left(1,\,\frac{p(h_{\rm{prop}}(\bm{s}_{i})|{\rm{data}})}{p(h_{\rm{cur}}(\bm{s}_{i})|{\rm{data}})}\right)$
$\displaystyle=\min\left(1,\,\frac{\prod_{k=1}^{p}p(\bm{W}_{h_{\rm{prop}}(\bm{s}_{i})}(\bm{s}_{i})\bm{y}_{k}|\bm{W}_{h_{\rm{prop}}(\bm{s}_{i})}(\bm{s}_{i})\bm{X}_{k},\,\bm{\beta}_{k},\,\sigma^{2})}{\prod_{k=1}^{p}p(\bm{W}_{h_{\rm{cur}}(\bm{s}_{i})}(\bm{s}_{i})\bm{y}_{k}|\bm{W}_{h_{\rm{cur}}(\bm{s}_{i})}(\bm{s}_{i})\bm{X}_{k},\,\bm{\beta}_{k},\,\sigma^{2})}\right)$
$\displaystyle=\min\left(1,\,\frac{\exp\left\\{-\frac{1}{2\sigma^{2}}\|\sum_{k=1}^{p}(\bm{W}_{h_{\rm{prop}}(\bm{s}_{i})}(\bm{s}_{i})\bm{y}_{k}-\bm{W}_{h_{\rm{prop}}(\bm{s}_{i})}(\bm{s}_{i})\bm{X}_{k}\bm{\beta}_{k})\|^{2}\right\\}}{\exp\left\\{-\frac{1}{2\sigma^{2}}\|\sum_{k=1}^{p}(\bm{W}_{h_{\rm{cur}}(\bm{s}_{i})}(\bm{s}_{i})\bm{y}_{k}-\bm{W}_{h_{\rm{cur}}(\bm{s}_{i})}(\bm{s}_{i})\bm{X}_{k}\bm{\beta}_{k})\|^{2}\right\\}}\right).$

In the above, $\bm{W}_{h}(\bm{s}_{i})$ is the matrix $\bm{W}(\bm{s}_{i})$
representing the weights of the neighbors when $h=h$.

  3. 3.

Using $u\sim{\rm{U}}(0,\,1)$, update $h$ as follows.

$\displaystyle
h_{\rm{cur}}(\bm{s}_{i})=\left\\{\begin{array}[]{ll}h_{\rm{prop}}(\bm{s}_{i})\quad(u<\alpha)\\\
h_{\rm{cur}}(\bm{s}_{i})\quad({\rm{otherwise}})\end{array}\right..$

Using these, the sampling method for BGWSR parameter estimation is represented
by the following Algorithm 1. Here $\bm{y},\,\bm{X},\,\bm{B}$ denote the set
of objective variable vectors, the set of explanatory variable matrices, and
the set of coefficient vectors, respectively,
$y=\\{\bm{y}_{1},\,\bm{y}_{2},\,\ldots,\,\bm{y}_{p}\\},\,X=\\{\bm{X}_{1},\,\bm{X}_{2},\,\ldots,\,\bm{X}_{p}\\},\,B=\\{\bm{\beta}_{1},\,\bm{\beta}_{2},\,\ldots,\,\bm{\beta}_{p}\\}$.

0:
$\bm{s}_{i},\,X,\,y,\,r,\,q,\,r_{1},\,q_{1},\,r_{2},\,q_{2},\,\sigma_{h}^{2},\,t_{\rm{max}}\,\,(i=1,\,2,\,\ldots,\,n;k=1,\,2,\,\ldots,\,p)$

0: $\widehat{B}$

1: Initialize
$\bm{\beta}_{k}^{(0)},\,T_{k}^{(0)},\,\Omega_{k}^{(0)},\,{\sigma^{2}}^{(0)},\,\lambda_{k,\,1}^{(0)},\,\lambda_{k,\,2}^{(0)},\,h(\bm{s}_{i})^{(0)}\,(i=1,\,2,\,\ldots,\,n;k=1,\,2,\,\ldots,\,p)$.

2: Initialize $\bm{W}^{(0)}$ using $w_{j}(\bm{s}_{i})|h(\bm{s}_{i})^{(0)}$.

3: Initialize $B^{(0)}$ using
$\bm{\beta}_{1}^{(0)},\,\bm{\beta}_{2}^{(0)},\,\ldots,\,\bm{\beta}_{p}^{(0)}$.

4: Initialize ${\bm{\Sigma}_{k}^{-1}}^{(0)}$ using
$T_{k}^{(0)},\,\Omega_{k}^{(0)}$.

5: for $t=0,\,1,\,\ldots,\,t_{\rm{max}}$ do

6: for $(k=1,\,2,\,\ldots,\,p)$ do

7: Sampling
$\bm{\beta}_{k}^{(t+1)}\sim\pi(\bm{\beta}_{k}|\bm{X}_{k},\,\bm{y}_{k},\,W^{(t)},\,{\sigma^{2}}^{(t)},\,{\bm{\Sigma}_{k}^{-1}}^{(t)})$.

8: end for

9: Update $B^{(t+1)}$ using
$\bm{\beta}_{1}^{(t+1)},\bm{\beta}_{2}^{(t+1)},\,\ldots,\,\bm{\beta}_{p}^{(t+1)}$.

10: for k=1, 2, …, p do

11: Sampling
$T_{k}^{(t+1)}\sim\pi(T_{k}|\bm{\beta}_{k}^{(t+1)},\,{\sigma^{2}}^{(t)},\,\lambda_{k,\,1}^{(t)})$.

12: Sampling
$\Omega_{k}^{(t+1)}\sim\pi(\Omega_{k}|\bm{\beta}_{k}^{(t+1)},\,{\sigma^{2}}^{(t)},\,\lambda_{k,\,2}^{(t)})$.

13: Update ${\bm{\Sigma}_{k}^{-1}}^{(t)}$ using
$T_{k}^{(t)},\Omega_{k}^{(t)}$.

14: end for

15: Sampling
${\sigma^{2}}^{(t+1)}\sim\pi(\sigma^{2}|X,\,y,\,B^{(t+1)},\,W^{(t)},\,r,\,q)$.

16: for k=1, 2, …, p do

17: Sampling
${\lambda_{k,\,1}^{2}}^{(t+1)}\sim\pi(\lambda_{k,\,1}^{2}|T_{k}^{(t+1)},\,r_{1},\,q_{1})$.

18: Sampling
${\lambda_{k,\,2}^{2}}^{(t+1)}\sim\pi(\lambda_{k,\,2}^{2}|\Omega_{k}^{(t+1)},\,r_{2},\,q_{2})$.

19: Sampling
$h(\bm{s}_{i})^{(t+1)}\propto\pi(h(\bm{s}_{i})|y,\,X,\,B^{(t+1)},\,{\sigma^{2}}^{(t+1)},\,\sigma_{h}^{2})$
by MH method.

20: end for

21: Update $W^{(t+1)}$ using $w_{j}(\bm{s}_{i})|h(\bm{s}_{i})^{(t+1)}$.

22: end for

Algorithm 1 Algorithm for BGWSR parameter sampling

###  2.4 Objective variable Prediction

The weight parameter $w_{j}(\bm{s}_{i^{*}})$ at any prediction location
$\bm{s}_{i^{*}}\in D$ is given by

$\displaystyle w_{j}(\bm{s}_{i^{*}})=$
$\displaystyle\left\\{\begin{array}[]{ll}1-\left(\frac{d_{i^{*}j}}{h(\bm{s}_{j})}\right)^{2}&\left(d_{i^{*}j}<h(\bm{s}_{j})\right)\\\
0&\left({\rm{otherwise}}\right)\end{array}\right.,\,(j=1,\,2,\,\ldots,\,n).$
(6)

Using this, we predict the coefficient
$\bm{\beta}(\bm{s}_{i^{*}})=(\beta_{1}(\bm{s}_{i^{*}}),\,\beta_{2}(\bm{s}_{i^{*}}),\,\ldots,\,\beta_{p}(\bm{s}_{i^{*}}))^{\prime}$
at the prediction location $\bm{s}_{i^{*}}$ by the following formula.

$\displaystyle\widehat{\bm{\beta}}(\bm{s}_{i^{*}})=\sum_{j=1}^{n}w_{j}(\bm{s}_{i^{*}})\bm{\beta}(\bm{s}_{j}).$
(7)

Using $\widehat{\bm{\beta}}(\bm{s}_{i^{*}}),\,\bm{x}(\bm{s}_{i^{*}})$, we
predict the objective variable $y(\bm{s}_{i^{*}})$ at the prediction location
$\bm{s}_{i^{*}}$ with the following formula.

$\displaystyle\widehat{y}(\bm{s}_{i^{*}})=\bm{x}(\bm{s}_{i^{*}})^{\prime}\widehat{\bm{\beta}}(\bm{s}_{i^{*}}).$
(8)

Therefore, the algorithm for predicting the coefficients and the objective
variable at the prediction locations is given by algorithm 2, where $n^{*}$ is
the number of locations predicted.

0:
$\bm{s}_{i},\,\bm{s}_{i^{*}},\,\bm{X}_{k},\,\bm{X}_{k}^{*},\,\bm{y}_{k},\,r,\,q,\,r_{1},\,q_{1},\,r_{2},\,q_{2},\,\sigma_{h}^{2},\,t_{\rm{max}}\,(i=1,\,2,\,\ldots,\,n;i^{*}=1,\,2,\,\ldots,\,n^{*};k=1,\,2,\,\ldots,\,p)$

0:
$\widehat{\bm{\beta}}(\bm{s}_{i^{*}}),\,\widehat{y}(\bm{s}_{i^{*}})\,(i^{*}=1,\,2,\,\ldots,\,n^{*})$

1: for $t=0,\,2,\,\ldots,\,t_{\rm{max}}$ do

2: for $i^{*}=1,2,\ldots,n^{*}$ do

3: for $j=1,2,\ldots,n$ do

4: Caluculate $w_{j}(\bm{s}_{i^{*}})^{(t)}$ from equation (6).

5: end for

6: Sampling $\widehat{\bm{\beta}}(\bm{s}_{i^{*}})^{(t)}$ from equation (7).

7: Sampling $\widehat{y}(\bm{s}_{i^{*}})^{(t)}$ from equation (8).

8: end for

9: end for

Algorithm 2 Algorithm for predicting coefficients and objective variables at
the prediction location by BGWSR

##  3 Numerical studies

In this section, we describe the details and results of numerical studies
conducted to evaluate the prediction performance of the proposed method.

###  3.1 Settings

In data generation, we first generated the location, followed by the data
associated with the location. In all scenarios, we randomly generate $n=1000$
spatial locations $\bm{s}_{1},\,\bm{s}_{2},\,\ldots,\,\bm{s}_{n}$ in the
square domain $[-1,1]\times[0,2]$ as follows.

$\displaystyle\bm{s}_{i}=(s_{i1},\,s_{i2})^{\prime},\,s_{i1}\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}{\rm{U}}(-1,1),\,s_{i2}\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}{\rm{U}}(0,2).$

The objective variables and covariates are generated as follows, where $p=3$.

$\displaystyle
y(\bm{s}_{i})=\sum_{k=1}^{p}x_{k}(\bm{s}_{i})\beta_{k}(\bm{s}_{i})+\varepsilon(\bm{s}_{i}),$
$\displaystyle\quad i=1,\,2,\,\ldots,\,n,$ $\displaystyle
x_{k}(\bm{s}_{i})\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}{\rm{N}}(0,\,1),$
$\displaystyle\quad i=1,\,2,\,\ldots,\,n;\,k=1,\,2,\,\ldots,\,p,$
$\displaystyle\varepsilon(\bm{s}_{i})\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}{\rm{N}}(0,\,\sigma^{2}),$
$\displaystyle\quad i=1,\,2,\,\ldots,\,n.$

The true values of the coefficients consider the following three structures of
coefficients.

  1. (a1)

The coefficients differ smoothly.

  2. (a2)

The clusters are different in the left and right regions, and the coefficients
within the clusters are the same.

  3. (a3)

The clusters are different in the upper and lower regions, and the
coefficients within the clusters are the same.

In (a1), the values are different for each location. The correlation of
coefficients between adjacent locations is strong, and the coefficients have
similar values. This structure of the coefficients is assumed for GWR and
BGWR. In (a2) and (a3), we assume a cluster structure depending on location so
that spatially close locations belong to the same cluster and have the same
coefficient values. We generated the true value of the coefficient at (a1) by

$\displaystyle\bm{\beta}_{k}=(\beta_{k}(\bm{s}_{1}),\,\beta_{k}(\bm{s}_{2}),\,\ldots,\,\beta_{k}(\bm{s}_{n}))^{\prime}\sim{\rm{MN}}(\beta_{k}^{*}\bm{1}_{n},\,\sigma^{2}\bm{H}),\,\quad
k=1,\,2,\,\ldots,\,p.$

Here, we set the covariance matrix $\bm{H}\in\mathbb{R}^{n\times n}$ so that
there is a correlation between adjacent locations,

$\displaystyle\bm{H}=\begin{pmatrix}1&w_{2}(\bm{s}_{1})&\cdots&w_{n}(\bm{s}_{1})\\\
w_{1}(\bm{s}_{2})&1&\cdots&w_{n}(\bm{s}_{2})\\\ \vdots&\vdots&\ddots&\vdots\\\
w_{1}(\bm{s}_{n})&w_{2}(\bm{s}_{n})&\cdots&1\end{pmatrix},$

where $w_{j}(\bm{s}_{i})=f(d_{ij}|h)=\exp\\{-(d_{ij}^{2}/2h^{2})\\}$. We set
the mean vector of coefficients
$\bm{\beta}^{*}=(\beta_{1}^{*},\,\beta_{2}^{*},\,\beta_{3}^{*})^{\prime}$, the
parameter of error variance $\sigma^{2}$, and bandwidth $h$ as described
below.

$\displaystyle\bm{\beta}^{*}=(1,1,1)^{\prime},\quad\sigma^{2}=1^{2},\quad
h=0.5.$

We generated the true value of the coefficient at (a2) and (a3) by

$\displaystyle\bm{\beta}(\bm{s}_{i})=\bm{\beta}^{(g)},\quad\bm{s}_{i}\in
C_{g}.$

Let $C_{g}$ denote the set of locations belonging to the $g$-th cluster
($g=1,\,2$). We set the true values of the coefficients as follows.

$\displaystyle\bm{\beta}^{(1)}=(1,1,1)^{\prime},$
$\displaystyle\bm{\beta}^{(2)}=(2,2,2)^{\prime}.$

We set the clusters of locations based on

$\displaystyle\ref{a2}:\quad\bm{s}_{i}\in$
$\displaystyle\left\\{\begin{array}[]{ll}C_{1}\quad(s_{i1}\leq 0)\\\
C_{2}\quad(s_{i1}>0),\end{array}\right.$
$\displaystyle\ref{a3}:\quad\bm{s}_{i}\in$
$\displaystyle\left\\{\begin{array}[]{ll}C_{1}\quad(s_{i2}\leq 0)\\\
C_{2}\quad(s_{i2}>0).\end{array}\right.$

Figure 1 shows an example of coefficients generated in scenarios (a1), (a2)
and (a3).

(a) Scenario (a1)

(b) Scenario (a2)

(c) Scenario (a3)

Figure 1: Examples of coefficients for each scenario

Next, we describe the division rule between observed and prediction locations.
In the numerical studies, we extracted the predicted and observed locations
from the generated locations. We randomly select $50$ locations for the
prediction. We set the number of observation locations to $100$ and consider
the following two conditions for the selection of observation locations.

  1. (b1)

The observed location is uniformly obtained within the region.

  2. (b2)

The observed location is non-uniformly obtained within the region.

In (b1), we randomly select locations from the generated data. In (b2), we
extract locations from the generated data by setting the number of locations
to $1:4$ on the left and right sides of the entire region. (b1) is the
situation assumed by GWR and BGWR. (b2) is the situation assumed by the BGWSR,
where there are not many locations near a given point in the sparse left-hand
area. Therefore, in GWR and BGWSR, the estimated coefficients are expected to
differ among adjacent locations. Figure 2 shows an example of locations
generated in the scenarios of (b1) and (b2). The dashed lines indicate the
left and right sides of the target area, divided by the center $0$; in the
scenario of (b2), the observation locations on the left side are sparse.

(a) Scenario (b1)

(b) Scenario (b2)

Figure 2: Examples of observation locations for each scenario

Table 1 shows the scenarios for the numerical studies. The scenario combining
setting (a3) and (b2) is the same as the scenario with the rotated region in
setting 3 and is omitted in these numerical studies because it does not affect
the weights used to estimate coefficients. The situation assumed for GWR and
BGWR is Scenario 1; the situations assumed for BGWSR are Scenarios 2, 4, and
5.

Table 1: Scenarios | true values of coefficients |  | Left-right ratio of observation locations
---
(Whether there is sparse or dense)
Scenario1 | smoothly (a1) |  $1:1$ (b1)
Scenario2 | smoothly (a1) |  $1:4$ (b2)
Scenario3 | left and right clusters (a2) |  $1:1$ (b1)
Scenario4 | left and right clusters (a2) |  $1:4$ (b2)
Scenario5 | upper and lower clusters (a3) |  $1:4$ (b1)

###  3.2 Comparative methods and Hyperparameters settings

We use our proposed method, BGWSR, and the existing methods, GWR and BGWR, in
the numerical studies. We apply two BGWSR methods: BGWSR with Adjacency
Estimation (BGWSR-AE), which determines the adjacency range differently for
each location, and BGWSR, which uses a common adjacency range for all
locations. We set the hyperparameters in BGWSR and BGWSR-AE as follows:

$\displaystyle(r,\,q,\,r_{1},\,q_{1},\,r_{2},\,q_{2},\,a)^{\prime}=(0.1,\,0.1,\,0.1,\,0.1,\,0.1,\,0.1,\,3)^{\prime}.$

The parameters $r_{\rm{BGWR}},\,q_{\rm{BGWR}}$ for the prior distribution of
the error variance $\sigma^{2}$ in BGWR were similarly set by

$\displaystyle(r_{\rm{BGWR}},\,q_{\rm{BGWR}})^{\prime}=(0.1,\,0.1)^{\prime}.$

The method for determining the adjacent locations of BGWR is based on Ma et
al. (2021), sampling the bandwidth $h_{\rm{BGWR}}$ using the MH method. The
prior distribution of $h_{\rm{BGWR}}$ is the following Uniform distribution.

$\displaystyle h_{\rm{BGWR}}\sim{\rm{U}}(0,\,3).$

We selected $\\{0.10,\,0.15,\,\ldots,\,3.00\\}$ as the candidate bandwidth
$h_{\rm{GWR}}$ in GWR and determined it using 5-fold cross validation.

###  3.3 Performance measures

We evaluate in the numerical studies the estimation performance of the
coefficients and objective variable at the observed and prediction location.
We use the mean squared error (MSE) as the performance measure and define the
MSE at the observation location as

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left(\widehat{\beta}_{k}(\bm{s}_{i})-\beta_{k}(\bm{s}_{i})\right)^{2},$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left(\widehat{y}(\bm{s}_{i})-y(\bm{s}_{i})\right)^{2},$

where
$\widehat{\bm{\beta}}(\bm{s}_{i})=(\widehat{\beta}_{1},\,\widehat{\beta}_{2},\,\ldots,\,\widehat{\beta}_{p})^{\prime},\,\widehat{y}(\bm{s}_{i})$
are the estimated value of the coefficients and the objective variable at
observed locations; $\bm{\beta}(\bm{s}_{i}),\,y(\bm{s}_{i})$ are the true
value of the coefficient and the objective variables at observed locations.
The MSE at the prediction location is defined as

$\displaystyle\frac{1}{n^{*}}\sum_{i^{*}=1}^{n^{*}}\left(\widehat{\beta}_{k}(\bm{s}_{i^{*}})-\beta_{k}(\bm{s}_{i^{*}})\right)^{2},$
$\displaystyle\frac{1}{n^{*}}\sum_{i^{*}=1}^{n^{*}}\left(\widehat{y}(\bm{s}_{i^{*}})-y(\bm{s}_{i^{*}})\right)^{2},$

where $\bm{s}_{i^{*}}\in D$ is the prediction location, and
$\widehat{\bm{\beta}}(\bm{s}_{i^{*}}),\,\widehat{y}(\bm{s}_{i^{*}})$ are the
predicted value of the coefficients and the objective variables at prediction
location $\bm{s}_{i^{*}}$. The true value of the coefficient at location
$\bm{s}_{i^{*}}$ is $\bm{\beta}(\bm{s}_{i^{*}})$, and that of the objective
variable is $y(\bm{s}_{i^{*}})$.

###  3.4 Summary of Numerical Studies Results

In the numerical studies, we generated data for each scenario $10$ times each
and calculated the median of the resulting evaluation metrics for each method
applied. Table 2 shows the median MSE for each method.

In Scenario 1, with respect to the MSE for the observation locations, the
BGWSR performed the best for the objective variables and some of the
coefficients, and the BGWSR-AE performed equally well. Similarly, for the MSE
for the prediction location, the BGWSR performed best, and the BGWSR-AE
performed equally well for the objective variables and most of the
coefficients. Figure 3 shows an example of the true values of the coefficients
and the estimates generated by each method for the location of the observation
in Scenario 1. BGWSR-AE tended to produce estimates closer to the overall
average rather than abrupt changes in the coefficients. This is because the
application of Bayesian Fused Lasso fails to capture the rapid changes.

In Scenario 2, for the MSE for the observation locations, BGWSR-AE yielded the
best results for the objective variables and some coefficients. Based on the
results of the MSE for the prediction locations, BGWSR-AE and BGWSR had a high
performance for the objective variables and coefficients. In Scenario 2,
BGWSR-AE and BGWSR can guarantee the similarity of locations by Bayesian Fused
Lasso. Figure 4 shows an example of the true values of the coefficients and
the estimated values by each method for the two locations in Scenario 2.
BGWSR-AE and BGWSR estimate values close to the true values even in regions
where observed locations are sparse. BGWR and GWR have locations with
extremely large estimated values in sparse areas. This is because the BGWR and
GWR set the range of adjacent locations to be relatively small to match the
dense locations. Consequently, extreme values are estimated for locations
where the data size is small.

In Scenario 3, BGWSR-AE and BGWSR exhibited high performance for objective
variables and coefficients at the observed location. According to the MSE for
the prediction location, the BGWSR performed the best for the objective
variables and coefficients. Figure 5 shows the true values of the coefficients
for the observation locations for Scenario 3 and examples of the estimated
values for each method. Scenario 3 has a cluster structure for locations,
suggesting that BGWSR-AE and BGWSR can estimate coefficients more uniformly
within clusters than the other methods. However, there are many locations
where the coefficients are shrunk to $0$ as a whole compared to the true
values. Nevertheless, in the BGWR, the estimation is unstable near the edge of
the cluster. For the GWR, the estimated coefficients varied gently near the
edge of the cluster and within clusters. This suggests that the estimation
performance of BGWR and GWR deteriorates near cluster boundaries.

In Scenario 4, BGWSR-AE had high performance for the objective variables and
coefficients at the observation locations. For the MSE for the prediction
location, the BGWSR-AE performed the best for the objective variables and most
of the coefficients, and the BGWSR performed equally well. Figure 6 shows an
example of the true values of the coefficients and the estimates generated by
each method in Scenario 4. The results suggest that BGWSR-AE and BGWSR can
estimate similarity coefficients within clusters compared to the other
methods. Additionally, BGWSR-AE and BGWSR provide numerically stable
coefficients even in areas with sparse observation locations. This suggests
that the Bayesian Fused Lasso effectively considers the similarity of adjacent
locations. Nonetheless, for the BGWR and GWR, the estimation became unstable
near the edge of the cluster. In areas with sparse observation locations,
extremely large or small coefficients were predicted due to the small data
size.

In Scenario 5, relative to the MSE for the observation locations, BGWSR-AE,
and BGWSR had high performance on the objective variables and coefficients at
the observation locations. BGWSR performed the best for the prediction
locations. Meanwhile the GWR yielded the best results on the objective
variables, and the BGWSR performed similarly well. Figure 7 shows an example
of the true values of the coefficients and the estimated values produced by
each method in Scenario 5. The results suggest that BGWSR-AE and BGWSR can
estimate coefficients more similarly within clusters than BGWR and GWR. For
BGWSR-AE, BGWSR, and GWR, the prediction of coefficients in the upper clusters
yielded larger estimates compared to the true values across clusters due to
Bayesian Fused Lasso. For BGWR and GWR, the coefficient estimates were
unstable at the edge of the domain and in areas with sparse observation
locations.

Table 2: Median MSE for each method  senario | method | observation location |  | prediction location
---|---|---|---|---
$\bm{\beta}_{1}$ | $\bm{\beta}_{2}$ | $\bm{\beta}_{3}$ | $\bm{y}$ |  | $\bm{\beta}_{1}$ | $\bm{\beta}_{2}$ | $\bm{\beta}_{3}$ | $\bm{y}$
1 | BGWSR-AE | 0.025 | 0.052 | 0.027 | 0.037 |  | 0.037 | 0.057 | 0.042 | 0.189
BGWSR | 0.024 | 0.051 | 0.027 | 0.036 |  | 0.032 | 0.068 | 0.033 | 0.161
BGWR | 0.046 | 0.083 | 0.044 | 0.167 |  | 0.041 | 0.075 | 0.037 | 0.187
GWR | 0.033 | 0.027 | 0.035 | 0.110 |  | 0.046 | 0.035 | 0.046 | 0.164
2 | BGWSR-AE | 0.024 | 0.041 | 0.027 | 0.032 |  | 0.042 | 0.057 | 0.059 | 0.172
BGWSR | 0.032 | 0.040 | 0.031 | 0.032 |  | 0.046 | 0.062 | 0.056 | 0.165
BGWR | 0.046 | 0.054 | 0.034 | 0.156 |  | 0.048 | 0.068 | 0.065 | 0.180
GWR | 0.057 | 0.036 | 0.047 | 0.119 |  | 0.152 | 0.064 | 0.219 | 0.323
3 | BGWSR-AE | 0.045 | 0.031 | 0.048 | 0.054 |  | 0.070 | 0.083 | 0.080 | 0.223
BGWSR | 0.036 | 0.028 | 0.052 | 0.055 |  | 0.056 | 0.049 | 0.059 | 0.145
BGWR | 0.077 | 0.075 | 0.089 | 0.199 |  | 0.057 | 0.060 | 0.078 | 0.187
GWR | 0.046 | 0.044 | 0.065 | 0.167 |  | 0.081 | 0.068 | 0.068 | 0.259
4 | BGWSR-AE | 0.034 | 0.041 | 0.031 | 0.052 |  | 0.085 | 0.086 | 0.088 | 0.293
BGWSR | 0.053 | 0.054 | 0.048 | 0.053 |  | 0.088 | 0.095 | 0.089 | 0.306
BGWR | 0.090 | 0.062 | 0.084 | 0.224 |  | 0.109 | 0.094 | 0.078 | 0.345
GWR | 0.052 | 0.071 | 0.049 | 0.210 |  | 0.090 | 0.139 | 0.114 | 0.400
5 | BGWSR-AE | 0.045 | 0.040 | 0.053 | 0.039 |  | 0.078 | 0.086 | 0.104 | 0.278
BGWSR | 0.033 | 0.036 | 0.033 | 0.039 |  | 0.058 | 0.052 | 0.064 | 0.164
BGWR | 0.073 | 0.064 | 0.074 | 0.157 |  | 0.075 | 0.069 | 0.068 | 0.196
GWR | 0.049 | 0.047 | 0.078 | 0.142 |  | 0.067 | 0.074 | 0.088 | 0.160

Figure 3: True and estimated values of coefficients for observation locations
in Scenario 1

Figure 4: True and estimated values of coefficients for observation locations
in Scenario 2

Figure 5: True and estimated values of coefficients for observation locations
in Scenario 3

Figure 6: True and estimated values of coefficients for observation locations
in Scenario 4

Figure 7: True and estimated values of coefficients for observation locations
in Scenario 5

##  4 Application to land price data

In this section, we evaluate the prediction performance of the proposed method
by applying it to real data to demonstrate its usefulness.

###  4.1 Data Summary and Applicable Settings

The data used are the official land price data for Tokyo for the year $2021$,
published by the Ministry of Land, Infrastructure, Transport and Tourism
(Ministry of Land and Tourism, 2021). The data were obtained for $1219$
locations, and $114$ variables such as land use classification and type of
frontage road were obtained based on the location information of latitude and
longitude. Figure 8 plots the obtained data by location information and shows
the land prices in 2021 in colors. In this study, we randomly selected $130$
locations in the area around the city center, which is surrounded by the black
square in Figure 8. Subsequently, we randomly selected $100$ observation
locations and $50$ prediction locations. Figure 9 shows the observation
locations and prediction locations. The observation locations are mostly
available in the left half of the city center, while the data are less
available in the right half of the city center—the data are coarse-grained in
terms of observation locations. The variables used are shown in table 3. In
this application, latitude and longitude were used as location, price as the
objective variable, and the volume ratio, land area, and station distance as
explanatory variables.

Table 3: List of variables used in real data example | variables | details
---|---|---
location | latitude |  $Y$-coordinate obtained by World Geodetic System (sec)
longitude |  $X$-coordinate obtained by World Geodetic System (sec)
objective variable | price | land prices in 2021 (Yen/${\rm{m}}^{2}$, Yen/$10a$ for forest land)
covariates | floor area ratio | Ratio of total building area (total floor space) to site area (%)
land area | land area (${\rm{m}}^{2}$)
distance | Distance to the nearest station (m)
Figure 8: land price in Tokyo

(a) observed location

(b) prediction location

Figure 9: Locations used in the analysis

We set up the hyperparameters for each method and the method of determining
the neighboring locations in the same way as in the numerical studies of
refchap:numerical study.

###  4.2 Results

Table 4 shows the MSE of the objective variable prediction by each method at
the prediction location. The MSE of BGWSR was the smallest. It was followed by
BGWSR-AE with the smallest MSE and GWR with the largest MSE. Figure 10 shows
the true value of the price at the prediction location and the predicted value
by each method, colored according to the magnitude of the value and visualized
on the map. As shown in Figure 10, BGWSR-AE produced the closest prediction to
the true value. Compared to the other methods, BGWSR-AE did not produce
extremely small values as predictions, suggesting that the Bayesian Lasso
effect stabilizes the estimated values. However, BGWSR and GWR produced
predictions that were extremely large or small compared to the true value in
the vicinity of the city center and near the boundaries of the region.

Table 5 shows the estimated coefficients for each method at the observation
location, colored according to the magnitude of the value and visualized on a
map. The $\beta_{1},\,\beta_{2},\,\beta_{3}$ are the coefficients for the
floor area ratio, land area, and distance. First, focusing on the estimated
coefficients for the floor area ratio, the coefficients are large in the city
center and decrease gradually as one moves away from the city center, as is
the case for all the methods. For BGWSR-AE and BGWSR, the coefficients were
estimated to vary gradually around the city center. Nevertheless, for BGWSR
and GWR, the estimated coefficients differed between the city center and areas
far from it. Next, focusing on the estimated coefficients for land area, all
the methods had a negative value for the coefficient in urban centers.
However, the coefficients of BGWSR, BGWR, and GWR are close to $-1$ in the
city center, while those of BGWSR-AE are around $0.4$. Finally, in the BGWSR-
AE, the estimated coefficients for station distance were larger near Chiyoda-
ku and smaller away from there. Nonetheless, for BGWSR, BGWR, and GWR, the
coefficients were estimated to be large in the upper right of the region and
decreased toward the lower left.

Figure 11 shows a visualization of the predicted values of the objective
variable by BGWSR-AE, BGWSR, and BGWR, colored by the size of the range of the
$95\%$ credible interval. From Figure 11, BGWSR-AE and BGWSR have more
locations with smaller interval widths compared to BGWR. In particular, the
interval widths were smaller for BGWSR-AE and BGWSR even at the edge of the
region. The ratio of true values included in the $95\%$ confidence interval
was the same for all methods. This suggests that BGWSR-AE and BGWSR have less
uncertainty than BGWR. Furthermore, BGWSR-AE and BGWSR are useful because they
have smaller MSEs than BGWR and GWR.

Table 4: MSE of predicted and true values method | ${\rm{MSE}}_{\bm{y}}$
---|---
BGWSR-AE | 0.092
BGWSR | 0.082
BGWR | 0.135
GWR | 0.171

(a) TRUE

(b) BGWSR-AE

(c) BGWSR

(d) BGWR

(e) GWR

Figure 10: True and predicted price at the prediction location | $\beta_{1}$ | $\beta_{2}$ | $\beta_{3}$
---|---|---|---
BGWSR-AE |  |  |
BGWSR |  |  |
BGWR |  |  |
GWR |  |  |
Table 5: Estimated coefficients at observation locations

(a) BGWSR-AE

(b) BGWSR

(c) BGWR

Figure 11: the Range of Credible interval of price at the prediction location

##  5 Conclusion

This paper proposes a BGWSR by combining the idea of Bayesian Fused Lasso and
BGWR. The proposed method enables us to consider spatial correlations even
when the number of observation locations is small. Moreover, by introducing a
method of determining the range of adjacent locations for each location into
BGWSR, the prediction performance is effectively maintained even when the
observation locations are heterogeneous.

This paper shows the usefulness of the proposed method through numerical
studies. The results show that the proposed method has better prediction
performance for the coefficients and the objective variable compared to
estimated by GWR and BGWR, where observation locations are obtained
heterogeneously in space. The results suggest that the proposed method
effectively considers the similarity of the coefficients between adjacent
locations since the estimated values of the coefficients are close among the
adjacent locations. Furthermore, when the locations have a cluster structure,
the estimated values of the proposed method are the same coefficients within
the clusters. We also applied the existing and proposed methods to land price
data for Tokyo, where the density of observation locations was spatially
heterogeneous. The results showed that the proposed method performed better
than the existing methods in predicting the objective variables. Furthermore,
we evaluated the uncertainty of the predicted values of the objective
variables. Consequently, we confirmed that the proposed method has a smaller
$95\%$ credible interval than BGWR, even at the edge of the domain and in
areas with sparse observation locations, and that most locations include the
true value within the interval. The results suggest that the proposed method
has less uncertainty than BGWR.

The proposed method can be widely used for data obtained from heterogeneous
observation patterns of location, and when spatial autocorrelation is assumed
in the coefficients. For example, the soil data used in Koh et al. (2020) have
observation locations where there are sparse and dense densities of locations.
In addition, data for seismic intensity observation locations in earthquake
intensity prediction and infection status in epidemiology are often obtained
in a spatially heterogeneous pattern. The proposed method is thus useful for
such data.

BGWSR could be potentially improved in several aspects. First, it is necessary
to find a way to determine the function of the weights used to estimate the
adjacency locations. The numerical studies and real data examples presented in
this paper use Bi-square, but existing studies have proposed Boxcar, Gaussian,
and other weight functions (Fotheringham et al., 2003). It is necessary to
select the appropriate approach among these weighting function based on the
characteristics of the data to be applied. For example, when using Boxcar, the
predicted values of the coefficients tend to be similar across the whole set.
If the purpose of the estimation is to emphasize the mean structure of the
coefficients, the use of Boxcar therefore may be appropriate.

The second point of potential improvement is enhancing interpretability for
high-dimensional data. BGWSR uses a Laplace distribution for the prior
distribution of the coefficients. Hence, the estimates of the coefficients are
not exactly $0$ in the methods that use the posterior mean or posterior median
as the estimator of the coefficient. This may reduce interpretability. In the
case of multivariate data analysis, variable selection methods are often used
and known to improve interpretability (Tibshirani, 1996). For achieving this
issue, one of the directions is the use of Horseshoe prior for the prior
distribution (Carvalho et al., 2010), and Spike and Slab prior (Ročková and
George, 2018). This allows to drive the estimation of some coefficient
estimates to $0$ and is expected to improve interpretability. For BGWSR,
similarly changing the prior distribution is expected to improve
interpretability through variable selection when applied to high-dimensional
data.

## References

  * (1)
  * Ahmed and De Marsily (1987) Ahmed, S. and De Marsily, G.(1987). Comparison of geostatistical methods for estimating transmissivity using data on transmissivity and specific capacity, Water Resources Research, 23(9), 1717–1737.
  * Andrews and Mallows (1974) Andrews, D. F. and Mallows, C. L.(1974). Scale mixtures of normal distributions, Journal of the Royal Statistical Society: Series B (Methodological), 36(1), 99–102.
  * Bailey et al. (1995) Bailey, T. C., Gatrell, A. C. et al.(1995). Interactive Spatial Data Analysis, Longman Scientific & Technical Essex, England.
  * Boehm Vock et al. (2015) Boehm Vock, L. F., Reich, B. J., Fuentes, M., and Dominici, F.(2015). Spatial variable selection methods for investigating acute health effects of fine particulate matter components, Biometrics, 71(1), 167–177.
  * Brunsdon et al. (1996) Brunsdon, C., Fotheringham, A. S., and Charlton, M. E.(1996). Geographically weighted regression: A method for exploring spatial nonstationarity, Geographical Analysis, 28(4), 281–298.
  * Carvalho et al. (2010) Carvalho, C. M., Polson, N. G., and Scott, J. G.(2010). The horseshoe estimator for sparse signals, Biometrika, 97(2), 465–480.
  * Cline (1973) Cline, A. K.(1973). Curve fitting using splines under tension(Curve fitting by application of splines under tension, discussing polynomial interpolation drawbacks and linear system solution for unknown second derivatives), Atmospheric Technology, 3, 60–65.
  * Fischer and Getis (2010) Fischer, M. M. and Getis, A.(2010). Handbook of Applied Spatial Analysis: Software Tools, Methods and Applications, Springer, New York.
  * Fotheringham et al. (2003) Fotheringham, A. S., Brunsdon, C., and Charlton, M.(2003). Geographically Weighted Regression: the Analysis of Spatially Varying Relationships, John Wiley & Sons, New York.
  * Fotheringham et al. (2017) Fotheringham, A. S., Yang, W., and Kang, W.(2017). Multiscale geographically weighted regression (MGWR), Annals of the American Association of Geographers, 107(6), 1247–1265.
  * Jena et al. (2023) Jena, R., Pradhan, B., Gite, S., Alamri, A., and Park, H.-J.(2023). A new method to promptly evaluate spatial earthquake probability mapping using an explainable artificial intelligence (XAI) model, Gondwana Research, 123, 54–67.
  * Koh et al. (2020) Koh, E.-H., Lee, E., and Lee, K.-K.(2020). Application of geographically weighted regression models to predict spatial characteristics of nitrate contamination: Implications for an effective groundwater management strategy, Journal of Environmental Management, 268, 110646.
  * Krige (1951) Krige, D. G.(1951). A statistical approach to some basic mine valuation problems on the Witwatersrand, Journal of the Southern African Institute of Mining and Metallurgy, 52(6), 119–139.
  * Kyung et al. (2010) Kyung, M., Gill, J., Ghosh, M., and Casella, G.(2010). Penalized regression, standard errors, and bayesian lassos, Bayesian Analysis, 5(2), 369–412.
  * Ministry of Land and Tourism (2021) Ministry of Land, T., Infrastructure and Tourism(2021). National Land Price Information (Land Price Public Notice Data), . https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-L01-v3_0.html (acquired in 26/12/2023).
  * LeSage (2004) LeSage, J. P.(2004). A family of geographically weighted regression models, in Advances in Spatial Econometrics: Methodology, Tools and Spplications, Springer, Berlin Heidelberg, 241–264.
  * Li and Sang (2019) Li, F. and Sang, H.(2019). Spatial Homogeneity Pursuit of Regression Coefficients for Large Datasets, Journal of the American Statistical Association, 114(527), 1050–1062.
  * Ma et al. (2021) Ma, Z., Xue, Y., and Hu, G.(2021). Geographically weighted regression analysis for spatial economics data: A Bayesian recourse, International Regional Science Review, 44(5), 582–604.
  * Park and Casella (2008) Park, T. and Casella, G.(2008). The Bayesian lasso, Journal of the American Statistical Association, 103(482), 681–686.
  * Ročková and George (2018) Ročková, V. and George, E. I.(2018). The spike-and-slab lasso, Journal of the American Statistical Association, 113(521), 431–444.
  * Subedi et al. (2018) Subedi, N., Zhang, L., and Zhen, Z.(2018). Bayesian geographically weighted regression and its application for local modeling of relationships between tree variables, iForest-Biogeosciences and Forestry, 11(5), 542\.
  * Tibshirani (1996) Tibshirani, R.(1996). Regression shrinkage and selection via the lasso, Journal of the Royal Statistical Society: Series B (Methodological), 58(1), 267–288.
  * Tibshirani et al. (2005) Tibshirani, R., Saunders, M., Rosset, S., Zhu, J., and Knight, K.(2005). Sparsity and smoothness via the fused lasso, Journal of the Royal Statistical Society: Series B (Methodological), 67(1), 91–108.
  * Wang et al. (2023) Wang, X., Zhu, Z., and Zhang, H. H.(2023). Spatial heterogeneity automatic detection and estimation, Computational Statistics & Data Analysis, 180, 107667.
  * Yu et al. (2007) Yu, D., Wei, Y. D., and Wu, C.(2007). Modeling spatial dimensions of housing prices in Milwaukee, WI, Environment and Planning B: Planning and Design, 34(6), 1085–1102.
  * Zhong et al. (2023) Zhong, Y., Sang, H., Cook, S. J., and Kellstedt, P. M.(2023). Sparse spatially clustered coefficient model via adaptive regularization, Computational Statistics & Data Analysis, 177, 107581.


# Accelerating Electron Dynamics Simulations through Machine Learned Time
Propagators

Karan Shah  Attila Cangi

###### Abstract

Time-dependent density functional theory (TDDFT) is a widely used method to
investigate electron dynamics under various external perturbations such as
laser fields. In this work, we present a novel approach to accelerate real
time TDDFT based electron dynamics simulations using autoregressive neural
operators as time-propagators for the electron density. By leveraging physics-
informed constraints and high-resolution training data, our model achieves
superior accuracy and computational speed compared to traditional numerical
solvers. We demonstrate the effectiveness of our model on a class of one-
dimensional diatomic molecules. This method has potential in enabling real-
time, on-the-fly modeling of laser-irradiated molecules and materials with
varying experimental parameters.

Machine Learning, ICML

##  1 Introduction

Time-Dependent Density Functional Theory (TDDFT) (Runge & Gross, 1984) is a
widely used method to study the evolution of electronic structure under the
influence of time-dependent potentials. It is used to calculate various
excited state properties such as excitation energies (Adamo & Jacquemin,
2013), charge transfer (Maitra, 2017), stopping power (Yost et al., 2017),
optical absorption spectra (Jacquemin et al., 2011) and non-linear optical
properties (Goncharov, 2014). Due to its favorable balance between accuracy
and computational cost, TDDFT has been applied in various domains including
photocatalysis, biochemistry, nanoscale devices and the study of light-matter
interactions in general.

For weak perturbations, the linear response (LR) formalism of TDDFT is used to
calculate the excitation spectrum of a system. It is calculated using the
Casida equation (Casida, 1995) as an eigenvalue problem that describes the
first order response of the density. In contrast, the electron density is
directly propagated in time under the real time (RT) formalism. RTTDFT can be
used to calculate the nonlinear response of the density under strong
perturbations such as ultrafast electron dynamics with strong laser fields.

There are multiple components that must be decided for RTTDDFT calculations.
These include preparation of the initial density and orbitals, choice of the
exchange-correlation functional, the form of the time-dependent potential and
the choice of the time-propagation scheme. The time evolution is a significant
fraction of the computation required (Castro et al., 2004; Gómez Pueyo et al.,
2018).

Machine learning (ML) for accelerating scientific simulations is a rapidly
growing area of research (Carleo et al., 2019). A variety of unsupervised and
data-driven models have been developed to solve differential equations across
a wide range of domains (Karniadakis et al., 2021).

Neural operators (NOs) (Kovachki et al., 2023) are a class of models that map
function-to-function spaces, as opposed to the finite-dimensional vector space
mappings of neural networks. This is especially useful for partial
differential equation (PDE) problems where experimental or simulation data is
available. Fourier neural operators (FNOs) (Li et al., 2020) are a type of NO
which represent operator weights in Fourier space. The main advantages of FNOs
are that they generalize well across function spaces, and being resolution-
invariant, they can be used for inference on higher resolution grids than the
training set grids. FNOs have also been used for forward and inverse PDE
problems in various domains (Azizzadenesheli et al., 2024).

While many applications of ML have been developed for ground state DFT (Snyder
et al., 2012; Brockherde et al., 2017; Fiedler et al., 2022), there have been
relatively fewer efforts focused on TDDFT. Some applications of ML for TDDFT
include development of exchange-correlation potentials (Yang & Whitfield,
2023; Suzuki et al., 2020) and predicting properties such as spectra and
stopping power (Ward et al., 2023). In this work, we demonstrate the
effectiveness of FNOs in propagating the electron density in time under the
TDDFT framework. Instead of propagating orbitals in time as done
conventionally in terms of the time-dependent Kohn-Sham equations, we use the
FNO propagator to directly evolve the density. This has two advantages: the
computational cost does not scale with the number of orbitals and larger
propagation time steps can be used, thus using fewer iterations. Section 2
contains a brief description of TDDFT, FNOs and our proposed autoregressive
model. Section 3 contains results obtained with this model for the time
evolution of one-dimensional diatomic molecules under the influence of an
oscillating laser pulse. We show that the model can be generalized across
ionic configurations, is faster and more accurate than comparable numerical
simulations and can also be used for higher resolution grids. In Section 4, we
discuss the physical viability of the predicted densities, and then offer
conclusion in Section 5.

We use atomic units a.u. ($\hbar=m_{e}=e=1)$ unless specified otherwise.

##  2 Methods

###  2.1 Time-Dependent Density Functional Theory

The time-dependent many-body Schrödinger equation describes the dynamics of
systems composed of multiple interacting particles. A common problem is
interacting electrons in atoms and molecules, for which the Schrödinger
equation is

$i\frac{\partial}{\partial
t}\Psi(\underline{\bm{r}},t)=\hat{H}(t)\Psi(\underline{\bm{r}},t),$ (1)

where $\underline{\bm{r}}$ represents the collective coordinates of the $N$
electrons, $\Psi(\underline{\bm{r}},t)$ is the many-body wavefunction and
$\hat{H}$ is the Hamiltonian operator, representing the total energy of they
system. For this problem the form of $\hat{H}$ is
$\hat{H}(t)=\hat{T}+\hat{W}+\hat{V}(t)$, where
$\hat{T}=\sum_{j=1}^{N}-\nabla_{j}^{2}/2$ is the kinetic energy operator,
$\hat{W}=\frac{1}{2}\sum_{\begin{subarray}{c}j,k\\\ j\neq
k\end{subarray}}^{N}1/\left|\mathbf{r}_{j}-\mathbf{r}_{k}\right|$ represents
the electron-electron interaction, and
$\hat{V}(t)=v_{\mathrm{{\scriptscriptstyle ext}}}(\bm{r},t)$ is the time-
dependent external potential operator which includes the static external
potential due to the ions and an external time-dependent potential that drives
the system. Note that in writing this Schrödinger equation we are assuming the
Born-Oppenheimer approximation, i.e. we are treating the dynamics of the ions
as occurring on a much longer time scale than the dynamics of the electrons,
and can therefore treat the ions as classical point-like particles. The
initial state $\Psi(\underline{\bm{r}})=\Psi(\underline{\bm{r}},t_{0})$ is
obtained by solving the time-independent Schrödinger equation

$\hat{H}\Psi(\underline{\bm{r}})=E\Psi(\underline{\bm{r}})$ (2)

which is an eigenvalue problem where $E$ denotes the eigenvalues and
corresponds to the total energy of the system.

The many-body wavefunction $\Psi$ contains all the information about the
system and can be used to calculate its properties. However for 3 spatial
dimensions, we need to solve a system of $3N$ variables and the computational
cost scales exponentially with $N$, rendering it intractable for all but the
simplest systems. DFT and TDDFT allow us to solve Eq. 2 and Eq. 1 respectively
by reformulating the problem in terms of electron density rather than the
wavefunction, dramatically reducing the number of variables from $3N$ to $3$.

The Runge-Gross theorem (Runge & Gross, 1984) states that for a system with a
given ground state many-body wavefunction
$\Psi_{0}=\Psi_{0}\left(\bm{r},t_{0}\right)$, there exists a unique mapping
between the potential and the time-dependent density. The density can be
obtained by solving a system of fictitious non-interacting particles governed
by the time-dependent Kohn-Sham (TDKS) equations (van Leeuwen, 1998)

$\left[-\frac{1}{2}\nabla^{2}+v_{s}[n](\bm{r},t)\right]\phi_{j}(\bm{r},t)=i\frac{\partial\phi_{j}(\bm{r},t)}{\partial
t},j=1,...,N$ (3)

where the electron density

$\displaystyle n(\bm{r},t)=\sum_{j}\,|\phi_{j}(\bm{r},t)|^{2},j=1,...,N$ (4)

is the quantity of interest. The Kohn-Sham potential $v_{s}[n](\bm{r},t)$ is a
functional of the density $n$ which is a function of position $\bm{r}$ and
time $t$. The correspondence between density and potential is established
through the Kohn-Sham potential

$\displaystyle v_{\mathrm{{\scriptscriptstyle
S}}}[n](\bm{r},t)=v_{\mathrm{{\scriptscriptstyle
ext}}}(\bm{r},t)+v_{\mathrm{{\scriptscriptstyle
H}}}[n](\bm{r},t)+v_{\mathrm{{\scriptscriptstyle XC}}}[n](\bm{r},t).$ (5)

where $v_{\mathrm{{\scriptscriptstyle ext}}}(\bm{r},t)$ is the external
potential, $v_{\mathrm{{\scriptscriptstyle H}}}(\bm{r},t)$ is the Hartree
potential and $v_{\mathrm{{\scriptscriptstyle XC}}}(\bm{r},t)$ is the
exchange-correlation potential. The external potential is composed of the
potential of the ions and any time-dependent external perturbation such as a
laser field.

While in theory, $v_{\mathrm{{\scriptscriptstyle XC}}}(\bm{r},t)$ depends on
the density at all previous time steps, the adiabatic approximation is often
used in practice (Ullrich, 2012). We use the adiabatic local density
approximation (ALDA). Under ALDA, $v_{\mathrm{{\scriptscriptstyle
XC}}}(\bm{r},t)$ is approximated as:

$v_{\mathrm{{\scriptscriptstyle
XC}}}^{\text{ALDA}}[n,\Psi_{0},\Phi_{0}](\mathbf{r},t)=v_{\mathrm{{\scriptscriptstyle
XC}}}^{\text{LDA}}[n](\mathbf{r},t),$ (6)

where $v_{\mathrm{{\scriptscriptstyle XC}}}^{\text{LDA}}$ is the LDA exchange-
correlation potential in ground-state DFT.

These coupled equations are solved through iterative numerical algorithms
which represent a significant computational cost in the TDDFT workflow (Castro
et al., 2004).

####  2.1.1 Time Propagators

The general form of the time evolution operator $\hat{U}$ for a time domain
$T$ is given by:

$\phi_{i}(\mathbf{r},T)=\hat{U}(T,0)\phi_{i}(\mathbf{r},0).$ (7)

In practice, it is applied multiple times for shorter time periods over the
domain

$\hat{U}(T,0)=\prod_{i=0}^{N-1}\hat{U}\left(t_{i}+\Delta t_{i},t_{i}\right)$
(8)

with $t_{0}=0,t_{i+1}=t_{i}+\Delta t_{i}$, and $t_{N}=T$. Additionally, the
evolution operator is unitary $\hat{U}^{\dagger}(t+\Delta
t,t)=\hat{U}^{-1}(t+\Delta t,t)$, which is required for conserving the
electron density, and follows time-reversal symmetry $\hat{U}(t+\Delta
t,t)=\hat{U}^{-1}(t,t+\Delta t)$. A time propagation algorithm should account
for these properties. Multiple time propagation methods have been developed
for this evolution (Castro et al., 2004; Gómez Pueyo et al., 2018). For this
work, we use the Crank-Nicholson method (Crank & Nicolson, 1947) to generate
the reference data. We then train a FNO to directly propagate the density in
time as opposed to propagating single particle orbitals.

####  2.1.2 Model Systems

We consider an external potential

$\displaystyle
v_{\text{ext}}(\mathbf{r},t)=v_{\text{ion}}(\mathbf{r})+v_{\text{las}}(t)\,,$
(9)

where we simulate a class of one-dimensional diatomic molecules that serves as
a model system with two interacting electrons under the static ionic
potential:

$\displaystyle
v_{\text{ion}}(\mathbf{r})=-\frac{Z_{1}}{\sqrt{(\mathbf{r}-\frac{d}{2})^{2}+a^{2}}}-\frac{Z_{2}}{\sqrt{(\mathbf{r}+\frac{d}{2})^{2}+a^{2}}},$
(10)

where $Z_{1}$ and $Z_{2}$ denote the charge of the atomic wells, $d$ denotes
the bond length, and $a$ is a softening parameter for numerical stability.

The system is excited with a laser given by a time-dependent external
perturbing potential:

$\displaystyle v_{\text{las}}(t)=A\sin{\omega t}.$ (11)

With the dipole approximation, we can treat the laser as spatially constant
because the size of our molecule is much small than the wavelength of the
laser. Given the ground state prepared with $v_{\text{ion}}(\mathbf{r})$, we
propagate the density under the influence of the laser $v_{\text{las}}(t)$.

The domain is defined by [-L, L] for time [0, T] with discretization $\Delta
x$ and $\Delta t$. We use fixed boundary conditions

$\displaystyle\phi_{i}(-L,t)=\phi_{i}(L,t)=0,\quad\forall t\in[0,T].$ (12)

The ground state density, ground state potential terms and shape of the laser
pulse are shown in in Figure LABEL:fig:example_system. We use the LDA
functional due to its simplicity and accuracy for such simple one-dimensional
systems. The density evolution is calculated by solving the time-dependent
Kohn-Sham equations with the external potential $v_{\text{ext}}(r,t)$ using
the Crank-Nicholson scheme. The numerical simulation is performed using the
Octopus RT-TDDFT code (Tancogne-Dejean et al., 2020).

###  2.2 Fourier neural operators

Figure 3: Autoregressive FNO architecture for predicting the density at time
$t+1$ based on $k$ previous time steps.

Neural operators (NOs) (Lu et al., 2019; Anandkumar et al., 2020) extend
neural networks by mapping functions to functions instead of finite-
dimensional vectors. While a neural network maps $\mathbb{R}^{n}$ to
$\mathbb{R}^{m}$, a neural operator maps function spaces,
$G:\mathcal{A}\rightarrow\mathcal{U}$. Our goal is to approximate the non-
linear map $\mathcal{G}^{\dagger}:\mathcal{A}\rightarrow\mathcal{U}$ with a
neural operator $\mathcal{G}_{\theta}$, parameterized by
$\theta\in\mathbb{R}^{p}$.

Training involves observations $\\{(a_{i},u_{i})\\}_{i=1}^{N}$ where
$u_{i}=\mathcal{G}^{\dagger}(a_{i})$. The objective is to find parameters
$\theta^{*}$ minimizing the loss:

$\theta^{*}=\min_{\theta\in\mathbb{R}^{p}}\frac{1}{N}\sum_{i=1}^{N}\left\|u_{i}-\mathcal{G}_{\theta}(a_{i})\right\|_{\mathcal{U}}^{2}.$
(13)

A neural operator $\mathcal{G}_{\theta}(a)$ is defined by:

$\mathcal{G}_{\theta}(a)=Q(v_{L}(v_{L-1}(\dots v_{1}(P(a)))),$ (14)

with layers $v_{l+1}$ as:

$v_{l+1}(x)=\sigma_{(l+1)}\left(W_{l}v_{l}(x)+\left(\mathcal{K}_{l}(a;\lambda)v_{l}\right)(x)\right),$
(15)

where $\mathcal{K}_{l}$ is a non-local kernel integral operator:

$\left(\mathcal{K}_{l}(a;\lambda)v_{l}\right)(x)=\int_{\Omega_{l}}\kappa_{l}(x,y,a(x),a(y);\lambda)v_{l}(y)dy.$
(16)

The kernel function $\kappa_{l}(x,y,a(x),a(y))$ depends on $x,y,a(x)$, and
$a(y)$, and is parameterized by $\lambda$. $W_{l}$ is a learnable weights
matrix which corresponds to a linear transformation and $\sigma_{l+1}$ is a
component-wise non-linear activation function. $P(a)$ preprocesses $a$ into
higher dimensional space $v_{0}$, and $Q(v_{L})$ post-processes $v_{L}$ back
into $\mathcal{U}$.

FNOs (Li et al., 2020) use the Fourier transform for efficient computation.
With $\kappa(x-y)$ as a convolution operator, we use the convolution theorem:

$\left(\mathcal{K}_{l}(a;\lambda)v_{l}\right)(x)=\mathcal{F}^{-1}\left(\mathcal{F}\left(\kappa_{l}\right)\cdot\mathcal{F}\left(v_{l}\right)\right)(x),$
(17)

where $\mathcal{F}$ and $\mathcal{F}^{-1}$ are the Fourier transform and its
inverse, respectively. Parameterization in Fourier space is determined by the
number of modes $k_{max}$ and width of the convolutional layers, allowing
efficient kernel computation and effective capture of global patterns using
Fast Fourier Transform.

###  2.3 Modeling the time-dependent electron density with Fourier neural
operators

An autoregressive model is used to predict future values based on past
observations. The model takes a slice of $T_{\text{in}}$ time steps and
predicts the next time step. The input sequence is then updated to include the
new predicted value while discarding the oldest time step.

Let $\mathbf{n}_{t}$ be the density grid at time step $t$. The input to the
model at time step $t$ is a sequence of $T_{\text{in}}$ previous observations,
denoted as:

$\mathbf{N}_{t}=\left[\mathbf{n}_{t-T_{\text{in}}+1},\mathbf{n}_{t-T_{\text{in}}+2},\ldots,\mathbf{n}_{t}\right],$
(18)

with $\mathbf{N}_{0}$ consisting of the initial $T_{in}$ time steps beginning
with the ground state density $\mathbf{n}_{0}$.

The model predicts the next time step $\hat{\mathbf{n}}_{t+1}$ as:

$\hat{\mathbf{n}}_{t+1}=\mathcal{G}_{\theta}(\mathbf{N}_{t})$ (19)

where $\mathcal{G}_{\theta}(\cdot)$ represents the autoregressive FNO model
parameterized by $\theta$.

After predicting $\hat{\mathbf{n}}_{t+1}$, the input sequence is updated for
the next prediction. The new input sequence $\mathbf{N}_{t+1}$ is formed by
appending $\hat{\mathbf{n}}_{t+1}$ and removing the oldest time step
$\mathbf{n}_{t-T_{\text{in}}+1}$:

$\mathbf{N}_{t+1}=\left[\mathbf{n}_{t-T_{\text{in}}+2},\mathbf{n}_{t-T_{\text{in}}+3},\ldots,\mathbf{n}_{t},\hat{\mathbf{n}}_{t+1}\right]$
(20)

This architecture is shown in Figure 3. This process is repeated iteratively
to generate further predictions.

The loss function for this model is the mean squared error (MSE) between the
predicted and true densities, averaged over all time steps and all systems.
Let $\mathcal{D}$ denote the set of all systems, and let $\mathcal{T}$ denote
the set of all predicted time steps for a given system. The loss function
$\mathcal{L}(\theta)$ is defined as:

$\mathcal{L}(\theta)=\frac{1}{|\mathcal{D}|}\sum_{d\in\mathcal{D}}\frac{1}{|\mathcal{T}|}\sum_{t\in\mathcal{T}}\left\|\mathbf{n}_{t}^{(d)}-\hat{\mathbf{n}}_{t}^{(d)}\right\|_{2}^{2},$
(21)

where $\mathbf{n}_{t}^{(d)}$ is the true density and
$\hat{\mathbf{n}}_{t}^{(d)}$ is the predicted density at time step $t$ for
system $d$.

To ensure that the density is conserved through time, we calculate the norm of
density

$N=\int_{-\infty}^{\infty}n(x,t)dx\approx\sum_{i}\hat{n}_{t,i}^{(d)}\Delta x$
(22)

and add a loss term

$\mathcal{L}_{norm}(\theta)=\lambda\left(\sum_{i}\hat{n}_{t,i}^{(d)}\Delta
x-2\right)^{2},$ (23)

where $i$ denotes index over grid points $x_{i}$ and $\lambda$ is weight term.
N = 2 for the reference system.

To evaluate model performance, we use the mean absolute error (MAE) and mean
squared error (MSE) between the reference and predicted densities.

##  3 Results

As baseline, we generate a set of 729 systems with varying $Z_{1},Z_{2},d$
parameters with fixed softening parameter $a=1.0$. The domain parameters are
$L=6.0$ a.u. and $T=10$ a.u. (0.241 fs). The discretization used is $\Delta
x=0.05$ a.u. and $\Delta t=0.025$ a.u. (0.6 as). The laser parameters are
$A=0.75$ a.u. (corresponding to peak intensity $I=1.97\times 10^{16}\text{ W
cm}^{-2}$) and $\omega=2.0$ Ha (corresponding to $\lambda=22.78$ nm). Using
this reference set, we generate datasets with discretization $\Delta t=0.2$
a.u. (4.83 as), resulting in $T_{n}=51$ time slices, each on a grid of
$X_{n}=241$ spatial points. The training dataset consists of 600 systems, the
validation dataset consists of 10 systems and performance is evaluated on a
test dataset of 100 systems. The networks consist of 3 Fourier layers, trained
using the Adam optimizer with learning rate = 0.001, halving every 200 epochs.

###  3.1 Effect of varying input time slice

We first show that by incorporating more information in the input time slices,
the performance of the network increases. In this case, the network consists
of layers with width 64 and 16 modes. Prediction metrics are evaluated on
$T_{n}-T_{in}$ time slices. As the number of input time slices $T_{in}$
increases, the errors decrease at the expense of the initial data required. We
show the error metrics for four models with $T_{in}={5,10,15,20}$ in Table 1.
As a balance between the accuracy and data requirements, a larger model is
created with $T_{in}=10$, consisting of layers with 16 modes and width 128,
trained for 2000 epochs. This model serves as the baseline for further
comparisons. The reference and predicted densities along with errors for a
representative system are shown in Figure LABEL:fig:snapshots.

###  3.2 Comparison with coarse solver

We compare the accuracy and speed of the baseline model with a coarse solver.
The baseline model is trained on a grid with $\Delta t=0.2$ a.u., using data
generated with a finer grid $\Delta t=0.025$ a.u. Numerical results are
calculated on a comparable coarse grid $\Delta t=0.2$ a.u. in Octopus. Both
these results are compared with reference fine grid numerical calculations. As
show in Table 1, the baseline model is roughly $2\times$ faster than the
numerical solver with error reductions of about $4\times$ MAE and $10\times$
MSE.

###  3.3 Time offset

We show that the model generalizes well when used to predict the density
evolution on a time grid offset by $\frac{\Delta T}{2}=0.1$ a.u. The domain in
this case becomes [$0.1,T-0.1$]. As show in Table 1, the MAE=5.779$\times
10^{-3}$ and prediction time per step $t=1.75$ ms is comparable to the
performance of the baseline model, which was evaluated on the same time grid
as the training time grid.

###  3.4 Spatial super-resolution

One key advantage of FNOs is discretization invariance. We show the
performance of the model trained on a dataset with $\Delta x=0.05$ a.u.
($X_{n}=241$) on a higher resolution grid with $\Delta x=0.025$ a.u.
($X_{n}=481$). We get accurate predictions with a reasonable increase in
errors but with the same inference speed, as shown in Table 1. This is
especially useful because the time taken per step increases to 9 ms for the
numerical solver on this grid while the model inference time stays constant at
1.75 ms.

Model | MAE ($\times 10^{-3}$) | MSE ($\times 10^{-4}$) | Time (ms)
---|---|---|---
FNO $T_{in}$ = 05 | 8.595 | 4.18 | 1.58
FNO $T_{in}$ = 10 | 7.087 | 4.10 | 1.62
FNO $T_{in}$ = 15 | 6.101 | 3.16 | 1.60
FNO $T_{in}$ = 20 | 5.652 | 3.32 | 1.60
Baseline | 5.118 | 2.56 | 1.75
Octopus (Coarse) | 19.616 | 24.69 | 4.00
Time Offset | 5.779 | 2.78 | 1.75
Super-resolution | 7.335 | 6.56 | 1.75
Table 1: Performance metrics of various models. The time column indicates the
calculation/inference time to obtain the next density time slice.

##  4 Discussion

In this section, we summarize the numerical results and discuss the
performance of the FNO time propagator for calculating physically viable
densities. For the baseline case, the MAE in density is 5.118 $\times 10^{-3}$
with a prediction time of 1.75 ms per time step. This is in excellent
agreement to the reference data with a fraction of the computational cost.
This model can be used for offset time grids and for higher-resolution space
grids as well without additional training. To compare the nature of calculated
densities, we evaluate the model on the properties of time evolution operators
and use the density to calculate an observable.

###  4.1 Properties of time evolution operators

As discussed in 2.1.1, time propagators for the Kohn-Sham equations follow
certain properties.

  * •

Density conservation: The norm of the density
$N=\int_{-\infty}^{\infty}n(x,t)dx$ must be conserved in the model system. For
our system with two electrons, $N=2$. For the baseline system, the FNO
predicted densities have a mean norm of $1.99992\pm 0.00039$, averaged over
all predicted time steps over all systems. The deviation of the predicted norm
increases with time, as shown in Figure 7.

  * •

Time reversal symmetry: A stable time-propagator must follow time-reversal
symmetry. As a qualitative check, densities are predicted using the baseline
model in reverse time order. However, the MAE jumps to 0.0751, with even the
MAE on the training set jumping to 0.0741.

The FNO model follows density conservation, which is also constrained in the
loss function during training. However, time reversal symmetry is broken as
shown by the increase in MAE for reverse time order predictions.

###  4.2 Dipole Moment Calculation

$\mu(t)=\int_{-\infty}^{\infty}x\,n(x,t)\,dx\approx\sum_{i}x_{i}\,n_{t,i}\,\Delta
x$ (24)

To demonstrate that the model predicted densities can be used for calculating
observables, we also compare the dipole moments calculated using predicted and
reference densities. For the baseline, the dipole moment MSE is 0.02976 with a
mean average percentage error of 8.34 %. The error in dipole moment
accumulates across time, as shown in Figure 8.

Figure 7: Evolution of predicted density norm over time. The blue points
represent the norms of predicted density, averaged across all systems, at
varying time steps with the error bars representing the stand deviation. The
orange line represents the true norm. The subplot shows the absolute error
between reference and predicted norm in time. Figure 8: Absolute error between
the reference and predicted dipole moment, averaged across all systems, at
different time steps.

The performance of the FNO time propagator in calculating these physically
relevant quantities shows that while density predictions are accurate, the
model can be improved further by incorporating more physics-based constraints
such as time reversal symmetry and stronger density conservation. Since the
typical time scale for electron dynamics is in the range of a few hundred
femtoseconds at most, achieving small error accumulation on this time scale
seems feasible based on our results.

###  4.3 Conclusion

We show that machine learned time propagators have the potential to accelerate
TDDFT calculations. A promising direction would be the development of ML time
propagators that can predict the time evolution given just the ground state
density and the shape of a laser pulse. Such extensions would include having
the laser as a separate input for better generalization across laser
parameters and to encode time reversal symmetry in the model structure.
Extending this work to physics-informed three-dimensional models would enable
on-the-fly modeling of the electronic response properties of laser-excited
molecules and materials in various scattering experiments that are conducted
at photon sources around the globe. This would enable fast simulations that
generalize well over the input parameters of the experimental setup. Rapid
modeling would also enable the design of laser pulses to precisely control
quantum dynamics under quantum optimal control theory (Werschnik & Gross,
2007).

## Acknowledgments

This work was supported by the Center for Advanced Systems Understanding
(CASUS), which is financed by Germany’s Federal Ministry of Education and
Research (BMBF) and by the Saxon state government out of the State budget
approved by the Saxon State Parliament. We acknowledge funding from the
Helmholtz Association’s Initiative and Networking Fund through Helmholtz AI.

## References

  * Adamo & Jacquemin (2013) Adamo, C. and Jacquemin, D.  The calculations of excited-state properties with Time-Dependent Density Functional Theory.  _Chemical Society Reviews_ , 42(3):845–856, January 2013.  ISSN 1460-4744.  doi: 10.1039/C2CS35394F.  URL https://pubs.rsc.org/en/content/articlelanding/2013/cs/c2cs35394f.  Publisher: The Royal Society of Chemistry.
  * Anandkumar et al. (2020) Anandkumar, A., Azizzadenesheli, K., Bhattacharya, K., Kovachki, N., Li, Z., Liu, B., and Stuart, A.  Neural Operator: Graph Kernel Network for Partial Differential Equations.  In _ICLR 2020 Workshop on Integration of Deep Neural Models and Differential Equations_ , February 2020.
  * Azizzadenesheli et al. (2024) Azizzadenesheli, K., Kovachki, N., Li, Z., Liu-Schiaffini, M., Kossaifi, J., and Anandkumar, A.  Neural operators for accelerating scientific simulations and design.  _Nature Reviews Physics_ , 6(5):320–328, May 2024.  ISSN 2522-5820.  doi: 10.1038/s42254-024-00712-5.  URL https://www.nature.com/articles/s42254-024-00712-5.  Publisher: Nature Publishing Group.
  * Brockherde et al. (2017) Brockherde, F., Vogt, L., Li, L., Tuckerman, M. E., Burke, K., and Müller, K.-R.  Bypassing the Kohn-Sham equations with machine learning.  _Nature Communications_ , 8(1):872, October 2017.  ISSN 2041-1723.  doi: 10.1038/s41467-017-00839-3.  URL https://www.nature.com/articles/s41467-017-00839-3.  Publisher: Nature Publishing Group.
  * Carleo et al. (2019) Carleo, G., Cirac, I., Cranmer, K., Daudet, L., Schuld, M., Tishby, N., Vogt-Maranto, L., and Zdeborová, L.  Machine learning and the physical sciences.  _Reviews of Modern Physics_ , 91(4):045002, December 2019.  doi: 10.1103/RevModPhys.91.045002.  URL https://link.aps.org/doi/10.1103/RevModPhys.91.045002.  Publisher: American Physical Society.
  * Casida (1995) Casida, M. E.  Time-Dependent Density Functional Response Theory for Molecules.  In _Recent Advances in Computational Chemistry_ , volume 1, pp. 155–192. WORLD SCIENTIFIC, November 1995.  ISBN 978-981-02-2442-4 978-981-283-058-6.  doi: 10.1142/9789812830586˙0005.  URL http://www.worldscientific.com/doi/abs/10.1142/9789812830586_0005.
  * Castro et al. (2004) Castro, A., Marques, M. A. L., and Rubio, A.  Propagators for the time-dependent Kohn–Sham equations.  _The Journal of Chemical Physics_ , 121(8):3425–3433, August 2004.  ISSN 0021-9606.  doi: 10.1063/1.1774980.  URL https://doi.org/10.1063/1.1774980.
  * Crank & Nicolson (1947) Crank, J. and Nicolson, P.  A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type.  _Mathematical Proceedings of the Cambridge Philosophical Society_ , 43(1):50–67, 1947.  doi: 10.1017/S0305004100023197.
  * Fiedler et al. (2022) Fiedler, L., Shah, K., Bussmann, M., and Cangi, A.  Deep dive into machine learning density functional theory for materials science and chemistry.  _Physical Review Materials_ , 6(4):040301, April 2022.  doi: 10.1103/PhysRevMaterials.6.040301.  URL https://link.aps.org/doi/10.1103/PhysRevMaterials.6.040301.  Publisher: American Physical Society.
  * Goncharov (2014) Goncharov, V.  _Non-Linear Optical Response in Atoms, Molecules and Clusters: An Explicit Time Dependent Density Functional Approach_.  Springer, September 2014.  ISBN 978-3-319-08320-9.
  * Gómez Pueyo et al. (2018) Gómez Pueyo, A., Marques, M. A. L., Rubio, A., and Castro, A.  Propagators for the Time-Dependent Kohn–Sham Equations: Multistep, Runge–Kutta, Exponential Runge–Kutta, and Commutator Free Magnus Methods.  _Journal of Chemical Theory and Computation_ , 14(6):3040–3052, June 2018.  ISSN 1549-9618.  doi: 10.1021/acs.jctc.8b00197.  URL https://doi.org/10.1021/acs.jctc.8b00197.  Publisher: American Chemical Society.
  * Jacquemin et al. (2011) Jacquemin, D., Mennucci, B., and Adamo, C.  Excited-state calculations with TD-DFT: from benchmarks to simulations in complex environments.  _Physical Chemistry Chemical Physics_ , 13(38):16987–16998, 2011.  doi: 10.1039/C1CP22144B.  URL https://pubs.rsc.org/en/content/articlelanding/2011/cp/c1cp22144b.  Publisher: Royal Society of Chemistry.
  * Karniadakis et al. (2021) Karniadakis, G. E., Kevrekidis, I. G., Lu, L., Perdikaris, P., Wang, S., and Yang, L.  Physics-informed machine learning.  _Nature Reviews Physics_ , 3(6):422–440, June 2021.  ISSN 2522-5820.  doi: 10.1038/s42254-021-00314-5.
  * Kovachki et al. (2023) Kovachki, N., Li, Z., Liu, B., Azizzadenesheli, K., Bhattacharya, K., Stuart, A., and Anandkumar, A.  Neural Operator: Learning Maps Between Function Spaces With Applications to PDEs.  _Journal of Machine Learning Research_ , 24(89):1–97, 2023.  ISSN 1533-7928.  URL http://jmlr.org/papers/v24/21-1524.html.
  * Li et al. (2020) Li, Z., Kovachki, N. B., Azizzadenesheli, K., Liu, B., Bhattacharya, K., Stuart, A., and Anandkumar, A.  Fourier Neural Operator for Parametric Partial Differential Equations.  In _International Conference on Learning Representations_ , October 2020.
  * Lu et al. (2019) Lu, L., Jin, P., and Karniadakis, G. E.  DeepONet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators.  _CoRR_ , abs/1910.03193, 2019.
  * Maitra (2017) Maitra, N. T.  Charge transfer in time-dependent density functional theory.  _Journal of Physics: Condensed Matter_ , 29(42):423001, September 2017.  ISSN 0953-8984.  doi: 10.1088/1361-648X/aa836e.  URL https://dx.doi.org/10.1088/1361-648X/aa836e.  Publisher: IOP Publishing.
  * Runge & Gross (1984) Runge, E. and Gross, E. K. U.  Density-Functional Theory for Time-Dependent Systems.  _Physical Review Letters_ , 52(12):997–1000, March 1984.  doi: 10.1103/PhysRevLett.52.997.  URL https://link.aps.org/doi/10.1103/PhysRevLett.52.997.  Publisher: American Physical Society.
  * Snyder et al. (2012) Snyder, J. C., Rupp, M., Hansen, K., Müller, K.-R., and Burke, K.  Finding Density Functionals with Machine Learning.  _Physical Review Letters_ , 108(25):253002, June 2012.  doi: 10.1103/PhysRevLett.108.253002.  URL https://link.aps.org/doi/10.1103/PhysRevLett.108.253002.  Publisher: American Physical Society.
  * Suzuki et al. (2020) Suzuki, Y., Nagai, R., and Haruyama, J.  Machine learning exchange-correlation potential in time-dependent density-functional theory.  _Physical Review A_ , 101(5):050501, May 2020.  doi: 10.1103/PhysRevA.101.050501.  URL https://link.aps.org/doi/10.1103/PhysRevA.101.050501.  Publisher: American Physical Society.
  * Tancogne-Dejean et al. (2020) Tancogne-Dejean, N., Oliveira, M. J. T., Andrade, X., Appel, H., Borca, C. H., Le Breton, G., Buchholz, F., Castro, A., Corni, S., Correa, A. A., De Giovannini, U., Delgado, A., Eich, F. G., Flick, J., Gil, G., Gomez, A., Helbig, N., Hübener, H., Jestädt, R., Jornet-Somoza, J., Larsen, A. H., Lebedeva, I. V., Lüders, M., Marques, M. A. L., Ohlmann, S. T., Pipolo, S., Rampp, M., Rozzi, C. A., Strubbe, D. A., Sato, S. A., Schäfer, C., Theophilou, I., Welden, A., and Rubio, A.  Octopus, a computational framework for exploring light-driven phenomena and quantum dynamics in extended and finite systems.  _The Journal of Chemical Physics_ , 152(12):124119, March 2020.  ISSN 0021-9606.  doi: 10.1063/1.5142502.  URL https://doi.org/10.1063/1.5142502.
  * Ullrich (2012) Ullrich, C.  _Time-Dependent Density-Functional Theory: Concepts and Applications_.  Oxford Graduate Texts. Oxford University Press, Oxford, New York, February 2012.  ISBN 978-0-19-956302-9.
  * van Leeuwen (1998) van Leeuwen, R.  Causality and Symmetry in Time-Dependent Density-Functional Theory.  _Physical Review Letters_ , 80(6):1280–1283, February 1998.  doi: 10.1103/PhysRevLett.80.1280.  URL https://link.aps.org/doi/10.1103/PhysRevLett.80.1280.  Publisher: American Physical Society.
  * Ward et al. (2023) Ward, L., Blaiszik, B., Lee, C.-W., Martin, T., Foster, I., and Schleife, A.  Accelerating Electronic Stopping Power Predictions by 10 Million Times with a Combination of Time-Dependent Density Functional Theory and Machine Learning, November 2023.  URL http://arxiv.org/abs/2311.00787.  arXiv:2311.00787 [cond-mat].
  * Werschnik & Gross (2007) Werschnik, J. and Gross, E. K. U.  Quantum optimal control theory.  _Journal of Physics B: Atomic, Molecular and Optical Physics_ , 40(18):R175, September 2007.  ISSN 0953-4075.  doi: 10.1088/0953-4075/40/18/R01.  URL https://dx.doi.org/10.1088/0953-4075/40/18/R01.
  * Yang & Whitfield (2023) Yang, J. and Whitfield, J.  Machine-learning Kohn–Sham potential from dynamics in time-dependent Kohn–Sham systems.  _Machine Learning: Science and Technology_ , 4(3):035022, August 2023.  ISSN 2632-2153.  doi: 10.1088/2632-2153/ace8f0.  URL https://dx.doi.org/10.1088/2632-2153/ace8f0.  Publisher: IOP Publishing.
  * Yost et al. (2017) Yost, D. C., Yao, Y., and Kanai, Y.  Examining real-time time-dependent density functional theory nonequilibrium simulations for the calculation of electronic stopping power.  _Physical Review B_ , 96(11):115134, September 2017.  doi: 10.1103/PhysRevB.96.115134.  URL https://link.aps.org/doi/10.1103/PhysRevB.96.115134.  Publisher: American Physical Society.