# MemDA: Forecasting Urban Time Series with Memory-based Drift Adaptation

Zekun Cai<EMAIL_ADDRESS>0000-0002-5773-1395 The University of
TokyoTokyoJapan ,  Renhe Jiang<EMAIL_ADDRESS>The University of
TokyoTokyoJapan ,  Xinyu Yang<EMAIL_ADDRESS>Tencent
CorporationBeijingChina ,  Zhaonan Wang<EMAIL_ADDRESS>The
University of TokyoTokyoJapan ,  Diansheng Guo<EMAIL_ADDRESS>Tencent CorporationBeijingChina ,  Hill Hiroki Kobayashi
<EMAIL_ADDRESS>The University of TokyoTokyoJapan ,  Xuan Song
<EMAIL_ADDRESS>The University of TokyoTokyoJapan and  Ryosuke
Shibasaki<EMAIL_ADDRESS>The University of TokyoTokyoJapan

(2023)

###### Abstract.

Urban time series data forecasting featuring significant contributions to
sustainable development is widely studied as an essential task of the smart
city. However, with the dramatic and rapid changes in the world environment,
the assumption that data obey Independent Identically Distribution is
undermined by the subsequent changes in data distribution, known as concept
drift, leading to weak replicability and transferability of the model over
unseen data. To address the issue, previous approaches typically retrain the
model, forcing it to fit the most recent observed data. However, retraining is
problematic in that it leads to model lag, consumption of resources, and model
re-invalidation, causing the drift problem to be not well solved in realistic
scenarios. In this study, we propose a new urban time series prediction model
for the concept drift problem, which encodes the drift by considering the
periodicity in the data and makes on-the-fly adjustments to the model based on
the drift using a meta-dynamic network. Experiments on real-world datasets
show that our design significantly outperforms state-of-the-art methods and
can be well generalized to existing prediction backbones by reducing their
sensitivity to distribution changes.

time series, concept drift, domain adaptation, urban computing

††journalyear: 2023††copyright: acmlicensed††conference: Proceedings of the
32nd ACM International Conference on Information and Knowledge Management;
October 21–25, 2023; Birmingham, United Kingdom.††booktitle: Proceedings of
the 32nd ACM International Conference on Information and Knowledge Management
(CIKM ’23), October 21–25, 2023, Birmingham, United Kingdom††price:
15.00††isbn: 979-8-4007-0124-5/23/10††doi: 10.1145/3583780.3614962††ccs:
Computing methodologies Knowledge representation and reasoning††ccs: Computing
methodologies Artificial intelligence††ccs: Applied computing Figure 1.
Illustration of traffic speed concept drift in Beijing from May 8 to July 25,
2022. (a) Grayscale ribbon of road speeds, with highlighted data for three
sampled roads. (b) The ridgeline plot exhibits weekly changes in the road
speed distribution. (c) The ridgeline plot exhibits weekly changes in the
error distribution of the predictive model, and the subplot shows the trend of
the daily error.

##  1\. Introduction

Accurate and reliable forecasting of urban multivariate time series data, such
as traffic speed, travel demand, and electricity consumption, has been
intensively studied as an essential task since it contributes to reducing
traffic congestion, energy waste, and greenhouse gas emissions (Ruan et al.,
2020; Jia et al., 2017; Yao et al., 2018; Jiang et al., 2019, 2023b). The
advancement of machine learning, particularly deep learning (Vaswani et al.,
2017; Li et al., 2018), has facilitated both research and practical
applications in this field. However, the underlying assumption that these
algorithms rest on - independent and identical distributions (i.i.d.) - may
not always hold true in the fast-paced, dynamic urban environment teeming with
implicit event-induced uncertainty. As an illustrative example, consider the
changes in daily traffic speeds in Beijing as shown in Fig. 1. The grayscale
ribbon plot in Fig. 1(a) demonstrates that the road speed decreases over
Mid-2022, indicating changes in the road speed distribution. The evolution of
this distribution, as captured week-by-week in the ridgeline plots of Fig.
1(b), shows the shift from a concentration around 75 km/h to a broader
distribution around 10-75 km/h. The drift distribution imposes vulnerability
on models trained from earlier data as they are not expected to generalize
well to the new coming data, known as the Concept Drift problem (Widmer and
Kubat, 1996). As a result, Fig. 1(c) highlights an increase in the error of a
pre-trained model over time, with error distribution spreading from low to
medium and high. These changes reflect the realities of a city in flux, where
social, economic, and environmental factors can drastically affect urban
patterns, leading to concept drift. It is considered a primary reason for the
declining effectiveness of intelligent information systems, as the change in
data distribution weakens the induced correlation patterns between new and
past data (Lu et al., 2018). Addressing this challenge can enhance the equity,
safety, and sustainability of cities and communities, assisting urban
decision-makers and stakeholders in using data-driven approaches to tackle
urban problems amidst continuous social changes.

Learning from changing data has attracted considerable research interest. For
non-stationary time series with concept drift, studies are conducted to: (1)
detect the occurrence of concept drift by error rate-based (Gama et al., 2004;
Bifet and Gavalda, 2007) or data distribution-based (Kifer et al., 2004;
Bhatia et al., 2022) detectors; (2) train a new model, fine-tune the current
model, or ensemble the new and existing model (Li et al., 2022; Du et al.,
2021; You et al., 2021) to fit the latest data. However, such incremental
learning approaches can only adjust the model to narrow the gap among the
drifting distributions. Awaiting the accumulation of new data inevitably leads
to lagging. The latest adapted model could fail again if the drift continues.
And even many scenarios do not qualify for frequent retraining as periodic
updates are costly. Considering that urban models are often associated with
substantial economic benefits and sustainable social development, this prompts
us to explore an evolving approach beyond retraining that can be proactively
adjusted on-the-fly based on inputs.

Nevertheless, developing a self-configuring learner is non-trivial. A crucial
consideration is the computational efficiency. Addressing concept drift
implies the necessity for the model to process extensive historical data,
providing a comprehensive understanding of past data distribution to
accurately perceive ongoing drifts. However, the assimilation and retention of
large volumes of historical data impose a significant computational burden.
The model must therefore exhibit the capacity to manage historical data in an
intelligent manner, maintaining long-term historical data whilst picking the
most valuable information. Second, the navigation of the stability-plasticity
dilemma cannot be overlooked. An adaptive model with essential levels of
responsiveness to new trends must achieve an equilibrium between maintaining
prior learned knowledge (stability) and accommodating new data patterns
(plasticity). This poses high requirements in deep learning models as the
parameter space of the model can be extremely large. Searching this space for
the optimal configuration involves considerable degrees of freedom and can be
practically infeasible for large models.

To confront the identified challenges, we propose a novel drift adaptation
network designed around two key components: a dual memory module and a
strategically adjustable meta-dynamic network. Firstly, our model incorporates
a dual memory module to enhance computational efficiency while handling
extensive historical information. This module consists of a Replay Memory for
time-efficient drift embedding and a Pattern Memory for preserving significant
patterns over extended periods. Then a novel drift-adaptive meta-dynamic
network is utilized to generate adjustable parameters for bridging the encoder
and decoder, allowing the model to respond to changing data patterns with
minimal degrees of freedom. Our approach is both simple and general, as it
only requires tracking a few parameters to model the drift and can be highly
compatible with any modern urban deep learning models as a plug-and-play
module for handling concept drift. Compared to existing models, the new
network improves prediction under non-stationary time series in a robust and
self-generating manner. Our contribution is summarized in the following
sections.

  * •

We first propose to model the urban concept drift by finding learnable
components in non-stationarity, which facilitates a more nuanced understanding
of urban dynamics and introduces new perspectives for prediction in such
environments.

  * •

We propose MemDA, an on-the-fly adaptation architecture that implicitly
perceives drift without retraining and dynamically adjusts the corresponding
model parameters.

  * •

The proposed model is pluggable and backbone-agnostic. Extensive experiments
on four urban datasets demonstrate the superior performance of our
design111https://github.com/deepkashiwa20/Urban_Concept_Drift.

##  2\. Related Work

###  2.1. Concept Drift and Model Adaptation

Time-varying patterns in the data render obsolete the knowledge previously
gained from old and static data, thereby reducing the predictive capability of
the model for new data (Lu et al., 2018; Khamassi et al., 2018). As the gap
between the latest and upcoming data distribution is more likely narrow, the
fresh training model can be utilized as a tool to address the issue. There are
three types of approaches based on trade-offs for preserving historical
information: the straightforward way is to retrain a new model with the most
recent data. An explicit sliding window drift detector is commonly employed to
determine when to retrain (Bifet and Gavalda, 2007; Li et al., 2022).
Alternatively, in some cases, retaining and fine-tuning the old model can save
efforts in training the new model. New base predictors can be added to
existing collections of models, or data-adaptive learning models be created.
(You et al., 2021) designs multiple predictors using a deep ensemble approach
and trains a classifier to predict which predictor to use. (Zhao et al., 2020)
proposes a weighted method for reusing previous models, where each model is
associated with a metric that represents its reusability for the current
observation. (Du et al., 2021) uses a segmentation algorithm to divide the
time series and uses a transfer learning approach to integrate knowledge under
various concepts. However, training-based approaches are just attempting to
catch up with the evolution of the concept. They can only build models after
drift and inevitably suffer from model aging.

Normalizing the input into a normalized distribution is also used to deal with
non-stationarity. RevIN(Kim et al., 2021) to remove and rebuild the
discrepancy in the mean and standard deviation of the input sequences. The
Non-stationary Transformer(Liu et al., 2022) mitigates the over-
stationarization problem by calculating the normalized attention. Dish-TS(Fan
et al., 2023) learns sample-level normalization parameters. However, The
normalization-based approaches simply attempt to align the magnitudes of the
samples, with a weak capacity for modeling changes in data patterns and
conditional probability distributions.

###  2.2. Urban Time Series Forecasting

Building time series models to forecast future traffic speed, taxi demand,
electricity consumption, crowd in/outflow, etc., are widely studied as key
technologies for smart city applications. (Yao et al., 2018) and (Cui et al.,
2019) employ an LSTM structure to connect hidden embeddings for urban data
forecasting. (Lai et al., 2018), (Li et al., 2018), (Bai et al., 2020), and
(Zhao et al., 2019) use GRU for more efficient modeling of traffic volume and
speed. Meanwhile, convolution-based modules, such as 1D CNN and TCN, are shown
to be competent with temporal data as well. (Yu et al., 2018), (Shih et al.,
2019), (Wu et al., 2019), and (Wu et al., 2020) for predicting urban data
demonstrate the efficiency and effectiveness of such designs. Attention
mechanisms in natural language processing have also been adapted. (Yao et al.,
2019; Guo et al., 2019) design models with multiple inputs based on the
periodicity of the data (hours, days, or weeks) and subsequently use attention
to capture the dynamic temporal correlation and assign greater weight to
important moments. (Zheng et al., 2020) and (Xu et al., 2020) also propose
similar attention modules for better global/long-term temporal modeling.
However, mainstream efforts in urban data forecasting assume constant data
patterns. The pre-trained models are highly vulnerable to the drift in data.

Another line focuses on the prediction of urban data under events. (Song et
al., 2014) predicts human movements during earthquakes using a Hidden Markov
Model. (Jiang et al., 2019, 2023a) models the crowd dynamics during four big
events by mobility momentum; EAST-Net (Wang et al., 2022) enhances the
robustness of the deep model using heterogeneous information networks and
memory bank. However, prediction in the event case is closer to emergency
response, which is concerned with the effect of the model on localized
outliers, whereas the concept drift is focused on persistent changes in the
patterns underlying the data.

##  3\. Preliminaries

Concept Drift. Concept drift refers to the underlying distribution of
streaming data varies over time caused by changes in unobserved hidden
variables (Lu et al., 2018). It will lead to model performance degradation as
the statistical properties of the target domain do not match the history.
Formally, given a set of instances obtained at different timesteps
$\\{(x_{0},y_{0}),...,(x_{T},y_{T})\\}$, where $x_{t}$ is the feature and
$y_{t}$ is the target, concept drift at timestamp $t$ can be described as
$p_{0:t-1}(x,y)\neq p_{t:T}(x,y)$, where $p(\cdot)$ denotes the joint
probability distribution. The source of concept drift can be decomposed into
two parts: changes in the marginal probability distribution $p(x)$ and changes
in the conditional probability distribution $p(y|x)$.

Adaptive Forecasting. Given a specified granularity of urban time series data,
the tensor $x_{t}\in\mathbb{R}^{N\times C}$ can be represented as the
observation at timestamp $t$, where $N$ denotes the number of urban sensors
(or regions) and C denotes the feature channel (e.g., traffic speed, crowd).
Then given the historical observations $\\{x_{i}|i\in[0,t]\\}$, building
adaptive model for the future $\alpha$ steps of urban data
$\hat{Y}_{t}=\\{\hat{x_{i}}|i\in[t+1,t+\alpha]\\}$ is to obtain such
parameters $\theta_{t}$ that can minimize the objective function
$\mathcal{L(\cdot)}$ as follows:

(1)
$\theta_{t}=\mathop{\arg\\!min}_{\theta_{t}}{\mathcal{L}(\hat{Y}_{t},\\{x_{i}|i\in[t+1,t+\alpha]\\})}.$

where $x_{i}$ are drawn from a time-varying distribution. $\theta_{t}$ is
learned from the training data and will be continuously adapted to test data.

Table 1. Notations

Symbol | Description
---|---
$x_{t}$ | Urban data observation tensor at timestamp $t$
$N,C$ | Dimension of $x_{t}$
$T$ | Total number of training samples
$p(x,y)$ | Joint probability distribution of variables $x$ and $y$
$\alpha$ | Prediction leading time
$\mathit{X}_{t}$ | Observations of urban data over time $[t-\alpha,t]$
$\mathcal{X}_{t}$ | A collection of observations constructed by periodicity
$p$ | Data sampling volume in one day
$k$ | the k-th look-back day for drift embedding
$P_{k}$ | Look-back timestamp for drift embedding, $P_{k}=p\times k$
$K$ | Total look-back days
$\mathcal{F}$ | Neural network functions
$\mathit{Z}_{t}$ | Embedding result of $\mathit{X}_{t}$ by Encoder
$C_{e}$ | Dimension of embedding generated by encoder
$\mathit{E}_{t}$ | Embedding result of $\mathcal{X}_{t}$ by Encoder and Replay Memory
$\mathit{M}$ | Pattern Memory
$L,D$ | Dimension of Pattern Memory
$\mathit{V}_{t}$ | The query result of $\mathit{E}_{t}$ returned by Pattern Memory
$H_{t}$ | Drift Embedding result of $\mathcal{X}_{t}$ by Dual-Memory
$W_{t}$ | Drift adaptation parameters
$N_{s}$ | Dimension of similarity generated by NTN
$I(X;Y)$ | Mutual information between variables $X$ and $Y$

Figure 2. The framework of the memory-based drift adaptation model.

##  4\. Methodology

###  4.1. Rationale

Although the distribution of urban data may change due to a variety of
factors, the inherent periodicity of urban dynamics follows a consistent
temporal regularity(Gonzalez et al., 2008). This persistent feature enables
the model to maintain its predictive capabilities, as the accumulation and
integration of fresh patterns facilitate adaptive learning even when
distribution undergoes significant transformations. It serves as an anchor for
our model, guiding the model toward the interpretation and management of
concept drift. Moreover, we acknowledge the sustained features to be encoded
in the pre- and post-drift data retain significant commonalities. These
sustained features can be efficiently captured by the stable component of the
model, which means modifications to some key model parameters should be
prioritized instead of updating all parameters. This strategy strikes a
balance between stability and plasticity, ensuring the resilience of the model
to varying data patterns and distributions.

Motivated by these, we design a Memory-based Drift Adaptation (MemDA) network
illustrated in Fig. 2, where we first propose to use a Dual-Memory to encode
periodicity implied in the long sequential inputs to counteract the drift, and
later use meta dynamic networks to generate adaptation parameters governing
the integration of input information. The notations are listed in Table 1.

###  4.2. Dual-Memory for Drift Embedding

Periodic inputs are essential for coping with urban data dynamics, while long
sequential inputs are the informational foundation for perception to drift.
Therefore, as shown in Fig. 3(a), we accordingly construct Drift-Aware Inputs
$\mathcal{X}_{t}$ along the time axis. Assuming that the data are sampled at a
frequency of $p$ per day, the inputs are constructed by fetching up to $K$
days of data as follows:

$\displaystyle\mathcal{X}_{t}$
$\displaystyle=\\{\mathit{X}_{t}\\}\cup\\{\mathit{X}_{t-P_{k}},\mathit{X}_{t+\alpha-
P_{k}}|k\in[1,K]\\},$ where $P_{k}=p\times k$ and $\mathit{X}_{t}$,
$\mathit{X}_{t-P_{k}}$, $\mathit{X}_{t+\alpha-P_{k}}$ are:
$\displaystyle\mathit{X}_{t}$
$\displaystyle=\\{x_{t-\alpha},x_{t-\alpha+1},...,x_{t}\\},$
$\displaystyle\mathit{X}_{t-P_{k}}$ $\displaystyle=\\{x_{t-\alpha-
P_{k}},x_{t-\alpha+1-P_{k}},...,x_{t-P_{k}}\\},$ (2)
$\displaystyle\mathit{X}_{t+\alpha-P_{k}}$
$\displaystyle=\\{x_{t+1-P_{k}},x_{t+2-P_{k}},...,x_{t+\alpha-P_{k}}\\}.$

The segment $\mathit{X}_{t}$ consists of the closest $\alpha$ observations up
to time $t$. The segments $\mathit{X}_{t-P_{k}}$ and $\mathit{X}_{t+\alpha-
P_{k}}$ consist of the same number of observations as $\mathit{X}_{t}$ but
from before and after the same clock $t$ prior to $k$ days. Here, $P_{k}$
represents the time gap between the backtrack and the present. Thus,
$\mathcal{X}_{t}$ contains $2K+1$ segments and each segment has $\alpha$
observations.

Figure 3. Design of model inputs switches from Plain Blocks (a) to Replay
Memory enhanced (b) inputs.

However, for such plain long time series (i.e., $(2K+1)\times\alpha$ steps),
traditional approaches (Zhang et al., 2017; Guo et al., 2019; Jiang et al.,
2021a) that utilize multi-branch to encode the sequence will be no longer
applicable as the computational complexity exponentially increases as the
input grows. Therefore, we propose to use a Replay Memory (RM) for efficient
drift embedding. RM is essentially a queue that stores all the hidden
embedding of $\mathit{X}_{t}$ computed by the encoder, forming
$\mathit{Z}\in\mathbb{R}^{T\times N\times C_{e}}$, where $T$ is the number of
training samples and $C_{e}$ is the output dimension of the encoder. During
training, the RM fetches the embeddings $\mathit{Z}_{t-P_{k}}$ and
$\mathit{Z}_{t+\alpha-P_{k}}$ of the past segments $\mathit{X}_{t-P_{k}}$ and
$\mathit{X}_{t+\alpha-P_{k}}$ to avoid duplicate computing. Without loss of
generality, we use $\mathcal{F}_{e}$ to denote the encoder network with
parameter $\theta_{e}$. As shown in Fig. 3(b), the embedding result
$\mathit{E}_{t}$ for Drift-Aware Inputs $\mathcal{X}_{t}$ is converted to the
$\mathit{X}_{t}$ encoding and RM replaying as:

(3)
$\mathit{E}_{t}=\underbrace{\\{\mathcal{F}_{e}(\mathit{X}_{t},\theta_{e})\\}}_{\text{Encoding}}\cup\underbrace{\\{\mathit{Z}_{t-P_{k}},\mathit{Z}_{t+\alpha-
P_{k}}|k\in[1,K]\\}}_{\text{Replaying}}.$

The embedding of $\mathit{X}_{t}$ is subsequently updated to RM by
$\mathit{Z}_{t}=\mathcal{F}_{e}(\mathit{X}_{t},\theta_{e})$. Note that
replayed memories are collected from the previous epoch while $\mathit{Z}_{t}$
is obtained at the current epoch (i.e., one epoch lag) when training. However,
as the epoch proceeds, they will gradually converge to the same embedding
space.

While periodic long inputs ensure the possibility of tracking data drift, the
model cannot extend segments indefinitely, and distribution information may
not be limited to periodicity. Therefore, we further augment the drift
embedding $E_{t}$ with a Pattern Memory (PM). The PM explicitly records the
prototype of normal patterns during training and serves to perceive out-of-
distribution during inference. Specifically, the Pattern Memory
$M\in\mathbb{R}^{L\times D}$ is a parameter matrix with $L$ real-valued $D$
dimensional prototypes. For the $i$-th embedding in $\mathit{E}_{t}$, $M$ is
queried to find similar prototypes by attention. Formally, we have:

(4) $\left\\{\begin{aligned} Q_{i}&=E_{t}^{i}\mathit{W}_{Q}+b_{Q}\\\
\varphi_{j}&=\frac{e^{Q_{i}\ast M_{j}}}{\sum_{j=1}^{L}e^{Q_{i}\ast M_{j}}}\\\
\mathit{V}_{t}^{i}&=(\sum_{j=1}^{L}\varphi_{j}\cdot
M_{j})\mathit{W}_{V}+b_{V},\end{aligned}\right.$

where $\varphi_{j}$ is the attention score corresponding to $j$-th prototype;
$\mathit{W_{Q}}$, $b_{Q}$, $\mathit{W_{V}}$, $b_{V}$ are linear transformation
parameters; $\mathit{V}_{t}^{i}$ denotes the query result of
$\mathit{E}_{t}^{i}$.

Then, $\mathit{E}_{t}$ and $\mathit{V}_{t}$ are concatenated to form the final
drift embedding $\mathit{H}_{t}=[\mathit{E}_{t},\mathit{V}_{t}]$, which
encodes the changes in data distribution and implies the essential information
for prediction. In general, the weighted fusion of embedding can integrate all
the information and can be fed to the decoder for future predictions. One
mainstream approach (Zhang et al., 2017; Yao et al., 2019; Bai et al., 2020)
is to use the trainable parameter to estimate the influence weights of each
segment as:

(5)
$\hat{Y}_{t}=\mathcal{F}_{proj.}(\sum_{i=1}^{4K+2}W^{i}\odot\mathit{H}_{t}^{i})$

where $\odot$ is Hadamard product; $W\in\mathbb{R}^{4K+2}$ are learnable
parameters, and $\mathcal{F}_{proj.}$ projects the result to the prediction.
Next, we introduce a dynamic network to make the key parameters $W$ self-
configurable.

###  4.3. Meta-Dynamic for Drift Adaptation

We characterize model adaptation in terms of Information Bottlenecks (IB). The
IB problem was introduced in (Tishby et al., 2000) as an information-theoretic
framework for learning. It considers quantifying the amount of information of
hidden representation $H$ extracted from the observation $X$ contains about
another relevant target signal $Y$ by using Mutual Information (MI). The goal
is to find a $H$ that minimizes the mutual information $I(X;H)$ while keeping
$I(H;Y)$ above the threshold. IB theory suggests that deep neural networks
construct Markov chains of the hidden layer representations about the input
$X$, and DNNs trained using SGD are essentially solving the compression-
prediction trade-off problem (Goldfeld and Polyanskiy, 2020), which means
deeper layers correspond to smaller $I(X;H)$ but larger $I(H;Y)$.

When it comes to the concept drift problem, in the encoder-decoder network,
the encoder essentially acts as an information bottleneck that filters the
noise from the data and forces the network to extract the typical patterns of
high-dimensional data, of which outputs are then transmitted to the decoder.
This makes $I(X;H)$ and $I(H;Y)$ a function of the en/decoder parameters. When
drift occurs, the encoder is more resistant to drift as it is more focused on
information compression. This leads us to fix the encoder parameters
$\theta_{e}$ to reduce the degrees of freedom while preserving the ability to
extract features over different target domains. For the adaptor, since we want
to boost the MI between the hidden embedding $\mathit{H}_{t}$ and $Y_{t}$, the
fusion parameter $W$ should be actively and dynamically adjusted to different
magnitudes under different concepts, i.e., varying with time $t$. For example,
the weight of historical information is evenly distributed when the data is
stable, while the influence of history is reduced when the data is out of
distribution. To generate the appropriate $W_{t}$, it is necessary to decide
how much each segment embedding contributes to the future, which means it
should be proportional to the MI of $\mathit{H}_{t}$ and $Y_{t}$, i.e.,

(6) $\displaystyle W_{t}$
$\displaystyle=softmax(I(\mathit{H}_{t}^{1};Y_{t}),...,I(\mathit{H}_{t}^{4K+2};Y_{t}))$
$\displaystyle=softmax(\mathbb{E}_{Y_{t}}\\{D_{KL}(p(\mathit{H}_{t}^{1})\|p(\mathit{H}_{t}^{1}|Y_{t})),...,$
$\displaystyle
D_{KL}(p(\mathit{H}_{t}^{4K+2})\|p(\mathit{H}_{t}^{4K+2}|Y_{t}))\\}).$

However, since the computation of MI requires knowledge of the joint and
marginal distributions of two random variables, the distributions that can be
depicted for high-dimensional urban datasets are stretched compared to the
true continuous distribution, especially when the joint distribution changes
in drift scenarios. This makes it tough to compute MI directly. Nevertheless,
it is possible to approximate the MI between high-dimensional variables using
a neural network with gradient descent optimization(Belghazi et al., 2018).
Therefore, we propose to use a meta-dynamic network to approximate the
contribution of each segment by dynamically generating the corresponding
$W_{t}^{i}$ conditional on the input. Specifically, to perceive the drift, we
construct three lists of temporal alignment tuples,

(7) $\displaystyle Pair_{x}$
$\displaystyle=\\{(\mathit{Z}_{t},\mathit{Z}_{t-P_{1}})\\}\cup\\{(\mathit{Z}_{t-P_{k}},\mathit{Z}_{t-P_{k+1}})\\},$
$\displaystyle Pair_{y}$ $\displaystyle=\\{(\mathit{Z}_{t+\alpha-
P_{k}},\mathit{Z}_{t+\alpha-P_{k+1}})\\},$ $\displaystyle Pair_{xy}$
$\displaystyle=\\{([\mathit{Z}_{t-P_{k}},\mathit{Z}_{t+\alpha-
P_{k}}],[\mathit{Z}_{t-P_{k+1}},\mathit{Z}_{t+\alpha-P_{k+1}}])\\}.$

The constructed embedding pairs are used to measure the evolution of the
distribution over the marginal probability $p(x)$ and the conditional
probability $p(y|x)$. Each pair serves as an anchor to facilitate determining
the contribution of each segment embedding to the prediction. Each
$(\mathit{Z}_{i},\mathit{Z}_{j})\in Pair$ is fed into a Neural Tensor Network
(NTN) (Socher et al., 2013) to model the relation:

(8) $\displaystyle sim_{k}$
$\displaystyle=\sigma(\mathit{Z}_{i}^{T}W_{S}\mathit{Z}_{j}+V[\mathit{Z}_{i},\mathit{Z}_{j}]^{T}+b_{S}),$
$\displaystyle W_{t}$
$\displaystyle=\mathcal{F}_{meta}([sim_{1},sim_{2},...]),$

where $W_{S}\in\mathbb{R}^{C_{e}\times C_{e}\times N_{s}}$,
$V\in\mathbb{R}^{N_{s}\times 2C_{e}}$, and $b_{S}\in\mathbb{R}^{N_{s}}$ are
weight tensor. $N_{s}$ is a hyperparameter that controls the drift scores
generated by the model for each pair from different perspectives. $\sigma$ is
an activation function. $\mathcal{F}_{meta}$ denotes a meta layer that can be
implemented in various ways, which generates fusion parameters based on the
similarity matrix to adjust to the drift on-the-fly. The generated $W_{t}$ are
loaded into the fusion layer replacing the static $W$ in Eq. 5, and a
multilayer CNN is employed as a decoder to project the result to the
prediction. The groundtruth values are taken to calculate the loss with
prediction to optimize the model.

###  4.4. Training Details

During the training phase, we train the model by random sampling batch, and
the RM is updated between epochs. During the testing phase, we feed the
samples into the model in chronological order, and the RM is updated
sequentially. Note that in the training phase, we maintain a full number of
embeddings in RM as many as $T$, while in the testing phase and online
environment, we only retain the most recent $p\times K$ embeddings to prevent
unlimited growth. The PM is updated as the gradient descends and is fixed in
the test.

##  5\. Experiments

Table 2. Summary of Experimental Datasets

Dataset | Period | #$N$ | Temporal | Drift
---|---|---|---|---
PeMS | 2020/01/01$\sim$2020/07/31 | 325 | 5 minutes | Sudden
Beijing | 2022/05/12$\sim$2022/07/25 | 3126 | 5 minutes | Sudden
Electricity | 2012/01/01$\sim$2012/06/30 | 370 | 1 hour | Incremental
COVID-CHI | 2019/07/01$\sim$2020/12/31 | 112 | 2 hour | Incremental

###  5.1. Datasets

We collected data from four cities covering different data sources and concept
drift types, which are presented in Table 2. The PeMS and Beijing datasets are
collected from the traffic speeds of major roads in California and Beijing.
The Electricity dataset contains the electricity consumption. And the COVID-
CHI (Wang et al., 2022) dataset is the demand for shared bicycles collected
from Chicago. To better understand how the data distribution has changed, we
show sampled data in Fig. 4, where each dataset exhibits different types of
concept drift. These drifts are caused by known events (i.e., COVID,
Seasonality) and result in significant variations in the data distribution.

Figure 4. Two nodes from each of the four datasets are drawn for time series
plots. The black vertical lines indicate the cut-off positions of the training
and testing datasets. Table 3. Performance Comparisons with state-of-the-art
Urban Time Series Prediction Models

Model | PeMS | Beijing | Electricity | COVID-CHI
---|---|---|---|---
RMSE | MAE | MAPE | RMSE | MAE | MAPE | RMSE | MAE | MAPE | RMSE | MAE | MAPE
CLD | 4.022 | 1.784 | 3.140% | 11.028 | 5.225 | 15.029% | 81.325 | 41.088 | 15.422% | 15.879 | 6.789 | 86.680%
ARIMA | 3.135 | 1.215 | 2.301% | 9.184 | 4.640 | 16.400% | 94.334 | 64.808 | 18.222% | 16.721 | 6.903 | 80.334%
AGCRN | 2.572 | 1.191 | 2.094% | 8.715 | 4.238 | 14.376% | 95.242 | 50.422 | 15.451% | 16.279 | 7.087 | 85.224%
DCRNN | 2.500 | 1.158 | 2.030% | 8.280 | 3.811 | 13.838% | 125.335 | 125.335 | 25.754% | 17.453 | 7.175 | 87.617%
MTGNN | 2.426 | 1.148 | 2.007% | 10.976 | 5.120 | 19.928% | 96.475 | 51.562 | 15.410% | 18.721 | 7.434 | 89.097%
STGCN | 2.556 | 1.212 | 2.114% | 10.284 | 5.000 | 15.284% | 132.209 | 72.802 | 22.580% | 22.462 | 8.687 | 85.064%
ASTGCN | 2.402 | 1.157 | 2.078% | 9.156 | 4.250 | 17.130% | 74.074 | 39.388 | 13.584% | 15.242 | 6.830 |  77.223%
GMAN | 2.849 | 1.298 | 2.308% | 11.555 | 6.448 | 21.941% | 138.436 | 81.285 | 26.263% | 22.462 | 8.687 | 85.064%
StemGNN | 2.414 | 1.144 | 2.030% | 8.575 | 4.294 | 15.605% | 103.375 | 56.803 | 19.366% | 17.231 | 7.132 | 85.32%
GW-Net | 2.385 | 1.125 | 1.976% | 7.754 | 3.564 | 11.916% | 97.087 | 51.684 | 16.120% | 16.210 | 6.857 | 80.038%
STTN | 2.441 | 1.168 | 2.040% | 9.886 | 4.748 | 17.150% | 93.125 | 50.671 | 15.776% | 14.853 | 6.934 | 79.155%
EAST-Net | 2.399 | 1.118 | 1.988% | 11.671 | 5.970 | 21.119% | 100.789 | 55.158 | 18.714% | 17.543 | 7.832 | 82.321%
MemDA | 2.297 | 1.053 | 1.861% | 6.720 | 3.192 | 9.913% | 67.413 | 34.814 | 12.186% | 14.003 | 6.115 | 72.088%
$\Delta$% | 3.70% | 5.82% | 5.78% | 13.34% | 10.43% | 16.81% | 8.99% | 11.61% | 10.29% | 5.73% | 9.93% | 9.93%

Table 4. Performance Comparisons with state-of-the-art Adaptive Models

Model | PeMS | Beijing | Electricity | COVID-CHI
---|---|---|---|---
RMSE | MAE | MAPE | RMSE | MAE | MAPE | RMSE | MAE | MAPE | RMSE | MAE | MAPE
Backbone | 2.385 | 1.125 | 1.976% | 7.754 | 3.564 | 11.916% | 97.087 | 51.684 | 16.120% | 16.210 | 6.857 | 80.038%
RevIN | 2.541 | 1.129 | 1.956% | 8.116 | 3.707 | 11.200% | 85.352 | 45.434 | 15.847% | 15.472 | 6.789 | 80.049%
$\Delta$% | No Improvement | No Improvement | 12.1% | 12.1% | 1.7% | 4.6% | 1.0% | -
Dish-TS | 2.519 | 1.144 | 1.973% | 7.738 | 3.586 | 10.499% | 73.509 | 41.465 | 13.726% | 14.807 | 6.435 | 74.182%
$\Delta$% | No Improvement | No Improvement | 24.3% | 19.8% | 14.8% | 14.8% | 12.7% | 14.0%
MemDA | 2.297 | 1.053 |  1.861% | 6.720 | 3.192 | 9.913% | 67.413 | 34.814 | 12.186% | 14.003 | 6.115 |  72.088%
$\Delta$% | 3.7% | 6.4% | 5.8% | 13.3% | 10.4% | 16.8% | 30.6% | 32.6% | 24.4% | 19.4% | 17.1% | 16.4%

Table 5. Variants Performance Evaluation

Variant | PeMS | Beijing | Electricity | COVID-CHI
---|---|---|---|---
RMSE | MAE | MAPE | RMSE | MAE | MAPE | RMSE | MAE | MAPE | RMSE | MAE | MAPE
Backbone | 2.385 | 1.125 | 1.976% | 7.754 | 3.564 | 11.916% | 97.087 | 51.684 | 16.120% | 16.210 | 6.857 | 80.038%
RM | 2.359 | 1.114 | 1.975% | 7.747 | 3.607 | 12.422% | 69.945 | 36.093 | 12.411% | 17.578 | 7.362 | 84.249%
RM+PM | 2.319 | 1.060 | 1.856% | 8.748 | 4.048 | 14.627% | 71.184 | 36.652 | 12.875% | 14.470 | 6.296 | 74.447%
MemDA | 2.297 | 1.053 | 1.861% | 6.720 | 3.192 | 9.913% | 67.413 | 34.814 | 12.186% | 14.003 | 6.115 |  72.088%

###  5.2. Settings

We annotate the time of events causing drift on each dataset. Based on this we
divide the data before the drift into the training set and the rest as the
test set, and further selected 20% of the training set as the validation set.
The cut-off positions for training and testing are shown in Fig. 4. As for the
model, the predicted leading time $\alpha$ is set to 12. The hyperparameters
are chosen based on the performance of the validation set. The look-back days
$K$ for drift embedding is set to 2, with the encoding dimension $C_{e}$ being
256, the number of prototypes $L$ of Pattern Memory being 20, and the
dimension $D$ being 32. The GW-Net (Wu et al., 2019) is chosen as the encoder
backbone for the model as its proven superior performance (Jiang et al.,
2021b). We set the $N_{s}$ in the drift adaptation module as 5, and the meta-
layer is implemented by a linear layer. As for the training, MAE is chosen to
be optimized using Adam, with a batch size set to 64 and a learning rate of
0.001. The best model is selected by early stop with a maximum of 200 epochs.
The experiments were performed on a GPU server with four GeForce RTX 3090
graphic cards. Root Mean Square Error (RMSE), Mean Absolute Error (MAE), and
Mean Absolute Percentage Error (MAPE) are selected to evaluate the results.

###  5.3. Baselines

####  5.3.1. Non-adaptive Models

We implemented state-of-the-art methods on urban data prediction for
comparison, including:

  * •

CopyLastDay (CLD): Observations from the same moment of the previous day are
used as predicted values.

  * •

ARIMA: AutoRegressive Integrated Moving Average is a widely used statistical
model for time series prediction by identifying patterns in data.

  * •

AGCRN (Bai et al., 2020): A node parameter learning and graph generation model
for traffic forecasting tasks. It automatically captures the spatial and
temporal correlation of covariates in time series data by adaptive modules.

  * •

DCRNN (Li et al., 2018): A model for modeling multivariate time-series
prediction tasks using bidimensional graph diffusion and recurrent neural
networks.

  * •

MTGNN (Wu et al., 2020): An advanced multivariate time-series prediction model
using active graph learning and multi-kernel TCN modules.

  * •

STGCN (Yu et al., 2018): A spatio-temporal network that combines graph
convolution with 1D convolution.

  * •

ASTGCN (Guo et al., 2019): A novel attention-based spatial-temporal graph
convolution model. It features three temporally separate components to model
periodic dependencies.

  * •

GMAN (Zheng et al., 2020): A model using multi-attentive networks with gated
fusion to fit nonlinear temporal correlations and adaptively fuse information.

  * •

StemGNN (Cao et al., 2020): A predictive model for modeling the spatial and
temporal dependence of multivariate time series in the spectral domain.

  * •

GW-Net (Wu et al., 2019): A state-of-the-art multivariate modeling model using
parametric graph inputs and a WaveNet-like temporal dilated structure.

  * •

STTN (Xu et al., 2020): A model in which various relationships related to
different factors in covariates are modeled jointly by spatial and temporal
transformers.

  * •

EAST-Net (Wang et al., 2022): A novel model for modeling event-induced non-
stationary prediction problems in multimodal data via heterogeneous mobility
information networks.

####  5.3.2. Adaptive Models

The latest advanced adaptive models for non-stationary time series are also
chosen to compare the effect on the concept drift problem, including:

  * •

RevIN (Kim et al., 2021): An instance normalization method for against
distribution shift that suppresses non-stationary information in one input
layer and restores it in the output layer.

  * •

Dish-TS (Fan et al., 2023): An advanced version of RevIN, which designs a
coefficient net that maps input sequences into learnable distribution to
relieve distribution change.

For a fair comparison, we select adaptive models that can also function
pluggable for comparison. The Backbone model is chosen as GW-Net based on its
optimal performance on the four datasets.

###  5.4. Performance Evaluation

####  5.4.1. Overall Performance Comparison

Table 3 shows the results of MemDA and baselines under three evaluation
metrics. The proposed MemDA consistently outperforms the best baseline model
on all four datasets, achieving a reduction in MAE of approximately 6% to 12%.
These results demonstrate the effectiveness of the adaptive model under drift
scenarios. Furthermore, we observed that the model performance on test data
varied. Among the two speed datasets, the Beijing dataset presented a greater
challenge due to the impact of lifting lockdowns, which caused traffic speeds
to drift from uniformly smooth to fluctuating. The proposed model exhibited a
greater advantage over the baselines when dealing with such new patterns,
resulting in an average improvement of 13% for the three evaluation metrics.
The PeMS dataset, on the other hand, drifted from congested to uncongested,
which can be interpreted as an oversampling from the original distribution.
Therefore, the concept drift problem may not cause significant issues, but our
approach still achieved a 6% improvement in the MAE metric. For the
Electricity datasets, ASTGCN, which also utilized periodically constructed
inputs, achieved the best results among the baseline models. However, our
proposed framework still achieves an additional 11% improvement, demonstrating
the potential of self-configuration capabilities in fully exploiting long-term
historical information. Finally, all baseline models showed poor results on
the COVID-CHI dataset, failing to match the simplest rule-based model
CopyLastDay in terms of MAE. Nevertheless, MemDA still stands out and shows
strong robustness in the drift.

Table 4 shows the results of MemDA and SOTA adaptive models. The experimental
results presented clearly demonstrate that our proposed MemDA model
outperforms the other existing adaptive models. Moving onto the analysis of
these results, RevIN and Dish-TS dynamically normalize the input data during
neural network computation, which can only alleviate issues caused by changes
in the marginal probability distribution $p(x)$. However, they overlook
changes in the conditional probability distribution, which makes them
disadvantaged in comparison. In addition, the Beijing and PeMS datasets
exhibit significant variations in data scale (refer to Fig. 4). Both RevIN and
Dish-TS fail to show improvement in these scenarios, indicating a potential
inability to handle drastic distribution drift effectively. In contrast, MemDA
shows improved performance, demonstrating its robust prediction logic that can
be particularly effective for prediction in challenging environments
characterized by substantial concept drift.

Figure 5. Visualization of the variation of the adaptation parameter $W_{t}$
on PeMS dataset. Where (a) is from the RM+PM, (b) is from the Meta, and (c) is
from the proposed model.

####  5.4.2. Ablation Test

A series of variant experiments are conducted to verify the effectiveness of
each module. The variant with only the original encoder network (GW-Net) is
referred to as the Backbone, to which we add RM and PM modules. The Meta
variant removes the meta-dynamic layer and directly generates $W_{t}$ through
a linear network connected to $E_{t}$. Results are shown in Table 5, where we
can observe that (1) Simply augmenting the long-term input via memory modules
does not always improve performance, with negative effects for Beijing and
COVID-CHI datasets. (2) Simply generating $W_{t}$ through the embedding of
input, as is done in the Meta variant, does not lead to satisfactory results.
In contrast, the MemDA evaluates drift by embedding pairs and generates
$W_{t}$ accordingly, demonstrating a proven adaptation capability and the
effectiveness of the meta-dynamic layer.

To interpret how MemDA copes with drift, we visualize the fusion parameter
$W_{t}$ in Fig. 5. Specifically, we show the solidified $W$ (the variant of
RM+PM) in (a), $W_{t}$ generated by the Simple Meta module (the variant of
RM+PM+Meta) in (b), and $W_{t}$ of MemDA in (c). The visualization results
demonstrate that the MemDA module enables the adaptation model to be aware of
changes in data distribution and to make proactive adjustments to the model
parameters in response to different concepts at different times. In contrast,
the solidified model can only maintain the prediction logic that fits the
training set, and the Simple Meta variant, although making some adjustments in
response to different concepts, demonstrates a weak adaptation capability that
cannot effectively track drift.

Figure 6. Sensitivity Analysis Figure 7. Case Study. (a) Average speed change
of Beijing test set. (b) The daily error of Backbone and MemDA model. (c,d)
Density error gray ribbon of Backbone/MemDA. (e,f) The model prediction
results of two roads at different time periods in the test set.

####  5.4.3. Sensitivity Analysis

We experimented the effects of Pattern Memory size $L$ and dimension $D$ of
MemDA on the PeMS dataset. The results are presented in Fig. 6(a,b). The model
achieves optimal performance with a memory size of 32. This optimal point
represents a balance between storing enough prototype patterns to accurately
detect out-of-distribution instances and the ability to generalize by
discarding noise. Similarly, optimal performance can also only be achieved
with a suitable memory dimension.

We further tested the effect of look-back days $K$. As a comparison, we
constructed a Plain model, in which naive uses multiple branches with a linear
fusion layer to process the input segments from different days (same as Fig. 3
(a)). The results are shown in Fig. 6(c). MemDA demonstrates a continuous
improvement of the effect with increasing look-back days, while the Plain
model fluctuates and does not show significant trends. This demonstrates that
dealing with concept drift is not a simple matter of adding inputs. And even
sometimes the increase in look-back days isn’t beneficial (Plain model 1 to
2). As the distribution of data changes over time, a model that heavily relies
on past data could be misled by outdated information that no longer applies to
the current. This is a risk when increasing input without an adaptive
mechanism.

A comparison of the running speed of the Plain and MemDA is shown in Fig.
6(d). It demonstrates that the computational pressure of the Plain model shows
an exponential increase as the number of backtracking days increases, and it
takes more than 10 times more time to look back at a week. While MemDA
enhanced by Replay Memory presents an efficient ability of slow linear
increase in computation time, which is crucial for enhancing model effects in
many real-time scenarios.

###  5.5. Generalizability Evaluation

The properties of MemDA enable it to be easily integrated with existing urban
data forecasting models. We utilize the backbone of three widely used models,
MTGNN, GMAN, and GW-Net as the encoder to test the encoder-agnostic
performance of the proposed model. The results are presented in Table 6, which
demonstrates the generalizability and applicability of MemDA across various
models.

Table 6. Encoder-agnostic Performance Evaluation

Backbone | Beijing | Electricity
---|---|---
RMSE | MAE | MAPE | RMSE | MAE | MAPE
MTGNN | 10.976 | 5.120 | 19.928% | 96.475 | 51.562 | 15.410%
+MemDA | 9.411 | 4.379 | 16.923% | 69.381 | 35.416 | 11.517%
$\Delta$% | 14.26% | 14.48% | 15.08% | 28.08% | 31.31% | 25.26%
GMAN | 11.555 | 6.448 | 21.941% | 138.436 | 81.285 | 26.263%
+MemDA | 10.005 | 6.011 | 19.041% | 91.030 | 49.400 | 15.072%
$\Delta$% | 13.41% | 6.78% | 13.22% | 34.24% | 39.23% | 42.61%
GW-Net | 7.754 | 3.564 | 11.916% | 97.087 | 51.684 | 16.120%
+MemDA | 6.720 | 3.192 | 9.913% | 67.413 | 34.814 | 12.186%
$\Delta$% | 13.37% | 10.43% | 16.81% | 30.56% | 32.64% | 24.40%

##  6\. Case Study

We do case studies from both macro and micro perspectives. Fig. 7 (a,b,c,d)
shows the city-wide average speed of the Beijing test dataset, the daily error
time series of the Backbone (GW-Net) and MemDA model, the daily error density
gray ribbon distribution of the Backbone and MemDA model, respectively. As the
average road speed in Beijing dwindles, traffic conditions progressively
worsen. This downturn significantly disrupts the predictability of traffic
patterns, triggering a substantial decline in the performance of the Backbone
model. This is reflected in the model’s increasing mean prediction error, with
a considerable spread ranging from low to high errors. In stark contrast, the
incorporation of MemDA into the Backbone manages to keep the error
distribution tightly controlled and centered on the lower values.

A microscopic analysis of the performance was focused on two representative
sample roads selected from the initial, intermediate, and final phases of the
test set. The zoomed-in views are shown in Fig. 7(e,f). As can be seen, both
roads with different speed patterns undergo significant pattern shifts over
time (marked with yellow boxes), and the Backbone model just barely works near
the training periods. However, as the concept drifts, Backbone’s predictions
stray further from the actual values. In essence, the Backbone model appears
to be memorizing speed patterns rather than predicting them. The introduction
of MemDA to it, however, imparts a remarkable transformation to the Backbone,
enabling it to track the changes in data patterns and always maintain
robustness after drift.

##  7\. Conclusion

In this study, we target the problem of concept drift in urban time series
data prediction. A general adaptive urban model leverages a dual-memory
module, and a meta-dynamic network is designed to enable self-configuration.
Extensive experiments demonstrate the effectiveness of the proposed framework.
In the future, we aim to explore online approaches for constructing robust
models as well as additional strategies for addressing spatial drift.

###### Acknowledgements.

This work was supported by Japan Science and Technology Agency (JST) SPRING,
Grant Number JPMJSP2108, and the 2022 Tencent Rhino-Bird Research Elite
Program.

## References

  * (1)
  * Bai et al. (2020) L Bai, L Yao, C Li, X Wang, and C Wang. 2020.  Adaptive Graph Convolutional Recurrent Network for Traffic Forecasting. In _34th Conference on Neural Information Processing Systems_.
  * Belghazi et al. (2018) Mohamed Ishmael Belghazi, Aristide Baratin, Sai Rajeswar, Sherjil Ozair, Yoshua Bengio, Aaron Courville, and R Devon Hjelm. 2018\.  Mine: mutual information neural estimation.  _arXiv preprint arXiv:1801.04062_ (2018).
  * Bhatia et al. (2022) Siddharth Bhatia, Arjit Jain, Shivin Srivastava, Kenji Kawaguchi, and Bryan Hooi. 2022\.  MemStream: Memory-Based Streaming Anomaly Detection. In _Proceedings of the ACM Web Conference 2022_. 610–621.
  * Bifet and Gavalda (2007) Albert Bifet and Ricard Gavalda. 2007.  Learning from time-changing data with adaptive windowing. In _Proceedings of the 2007 SIAM international conference on data mining_. SIAM, 443–448.
  * Cao et al. (2020) Defu Cao, Yujing Wang, Juanyong Duan, Ce Zhang, Xia Zhu, Conguri Huang, Yunhai Tong, Bixiong Xu, Jing Bai, Jie Tong, et al. 2020\.  Spectral temporal graph neural network for multivariate time-series forecasting. In _Proceedings of the 34th International Conference on Neural Information Processing Systems_. 17766–17778.
  * Cui et al. (2019) Zhiyong Cui, Kristian Henrickson, Ruimin Ke, and Yinhai Wang. 2019.  Traffic graph convolutional recurrent neural network: A deep learning framework for network-scale traffic learning and forecasting.  _IEEE Transactions on Intelligent Transportation Systems_ 21, 11 (2019), 4883–4894.
  * Du et al. (2021) Yuntao Du, Jindong Wang, Wenjie Feng, Sinno Pan, Tao Qin, Renjun Xu, and Chongjun Wang. 2021.  Adarnn: Adaptive learning and forecasting of time series. In _Proceedings of the 30th ACM International Conference on Information & Knowledge Management_. 402–411.
  * Fan et al. (2023) Wei Fan, Pengyang Wang, Dongkun Wang, Dongjie Wang, Yuanchun Zhou, and Yanjie Fu. 2023\.  Dish-TS: A General Paradigm for Alleviating Distribution Shift in Time Series Forecasting. In _Proceedings of the AAAI Conference on Artificial Intelligence_ , Vol. 37. 7522–7529.
  * Gama et al. (2004) Joao Gama, Pedro Medas, Gladys Castillo, and Pedro Rodrigues. 2004\.  Learning with drift detection. In _Brazilian symposium on artificial intelligence_. Springer, 286–295.
  * Goldfeld and Polyanskiy (2020) Ziv Goldfeld and Yury Polyanskiy. 2020.  The information bottleneck problem and its applications in machine learning.  _IEEE Journal on Selected Areas in Information Theory_ 1, 1 (2020), 19–38.
  * Gonzalez et al. (2008) Marta C Gonzalez, Cesar A Hidalgo, and Albert-Laszlo Barabasi. 2008\.  Understanding individual human mobility patterns.  _nature_ 453, 7196 (2008), 779–782.
  * Guo et al. (2019) Shengnan Guo, Youfang Lin, Ning Feng, Chao Song, and Huaiyu Wan. 2019.  Attention based spatial-temporal graph convolutional networks for traffic flow forecasting. In _Proceedings of the AAAI conference on artificial intelligence_ , Vol. 33. 922–929.
  * Jia et al. (2017) Chuanmin Jia, Shiqi Wang, Xinfeng Zhang, Shanshe Wang, and Siwei Ma. 2017.  Spatial-temporal residue network based in-loop filter for video coding. In _2017 IEEE Visual Communications and Image Processing (VCIP)_. IEEE, 1–4.
  * Jiang et al. (2021a) Renhe Jiang, Zekun Cai, Zhaonan Wang, Chuang Yang, Zipei Fan, Quanjun Chen, Kota Tsubouchi, Xuan Song, and Ryosuke Shibasaki. 2021a.  DeepCrowd: A deep model for large-scale citywide crowd density and flow prediction.  _IEEE Transactions on Knowledge and Data Engineering_ (2021).
  * Jiang et al. (2019) Renhe Jiang, Xuan Song, Dou Huang, Xiaoya Song, Tianqi Xia, Zekun Cai, Zhaonan Wang, Kyoung-Sook Kim, and Ryosuke Shibasaki. 2019.  Deepurbanevent: A system for predicting citywide crowd dynamics at big events. In _Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining_. 2114–2122.
  * Jiang et al. (2023a) Renhe Jiang, Zhaonan Wang, Yudong Tao, Chuang Yang, Xuan Song, Ryosuke Shibasaki, Shu-Ching Chen, and Mei-Ling Shyu. 2023a.  Learning Social Meta-knowledge for Nowcasting Human Mobility in Disaster. In _Proceedings of the ACM Web Conference 2023_. 2655–2665.
  * Jiang et al. (2023b) Renhe Jiang, Zhaonan Wang, Jiawei Yong, Puneet Jeph, Quanjun Chen, Yasumasa Kobayashi, Xuan Song, Shintaro Fukushima, and Toyotaro Suzumura. 2023b.  Spatio-temporal meta-graph learning for traffic forecasting. In _Proceedings of the AAAI Conference on Artificial Intelligence_ , Vol. 37. 8078–8086.
  * Jiang et al. (2021b) Renhe Jiang, Du Yin, Zhaonan Wang, Yizhuo Wang, Jiewen Deng, Hangchen Liu, Zekun Cai, Jinliang Deng, Xuan Song, and Ryosuke Shibasaki. 2021b.  Dl-traff: Survey and benchmark of deep learning models for urban traffic prediction. In _Proceedings of the 30th ACM international conference on information & knowledge management_. 4515–4525.
  * Khamassi et al. (2018) Imen Khamassi, Moamar Sayed-Mouchaweh, Moez Hammami, and Khaled Ghédira. 2018.  Discussion and review on evolving data streams and concept drift adapting.  _Evolving systems_ 9, 1 (2018), 1–23.
  * Kifer et al. (2004) Daniel Kifer, Shai Ben-David, and Johannes Gehrke. 2004\.  Detecting change in data streams. In _Proceedings of the Thirtieth international conference on Very large data bases-Volume 30_. 180–191.
  * Kim et al. (2021) Taesung Kim, Jinhee Kim, Yunwon Tae, Cheonbok Park, Jang-Ho Choi, and Jaegul Choo. 2021\.  Reversible instance normalization for accurate time-series forecasting against distribution shift. In _International Conference on Learning Representations_.
  * Lai et al. (2018) Guokun Lai, Wei-Cheng Chang, Yiming Yang, and Hanxiao Liu. 2018\.  Modeling long-and short-term temporal patterns with deep neural networks. In _The 41st international ACM SIGIR conference on research & development in information retrieval_. 95–104.
  * Li et al. (2022) Wendi Li, Xiao Yang, Weiqing Liu, Yingce Xia, Jiang Bian, et al. 2022\.  DDG-DA: Data Distribution Generation for Predictable Concept Drift Adaptation. In _Proceedings of the AAAI Conference on Artificial Intelligence_.
  * Li et al. (2018) Yaguang Li, Rose Yu, Cyrus Shahabi, and Yan Liu. 2018\.  Diffusion Convolutional Recurrent Neural Network: Data-Driven Traffic Forecasting. In _International Conference on Learning Representations_.
  * Liu et al. (2022) Yong Liu, Haixu Wu, Jianmin Wang, and Mingsheng Long. 2022\.  Non-stationary Transformers: Exploring the Stationarity in Time Series Forecasting. In _Advances in Neural Information Processing Systems_.
  * Lu et al. (2018) Jie Lu, Anjin Liu, Fan Dong, Feng Gu, Joao Gama, and Guangquan Zhang. 2018\.  Learning under concept drift: A review.  _IEEE Transactions on Knowledge and Data Engineering_ 31, 12 (2018), 2346–2363.
  * Ruan et al. (2020) Sijie Ruan, Jie Bao, Yuxuan Liang, Ruiyuan Li, Tianfu He, Chuishi Meng, Yanhua Li, Yingcai Wu, and Yu Zheng. 2020.  Dynamic public resource allocation based on human mobility prediction.  _Proceedings of the ACM on interactive, mobile, wearable and ubiquitous technologies_ 4, 1 (2020), 1–22.
  * Shih et al. (2019) Shun-Yao Shih, Fan-Keng Sun, and Hung-yi Lee. 2019.  Temporal pattern attention for multivariate time series forecasting.  _Machine Learning_ 108, 8 (2019), 1421–1441.
  * Socher et al. (2013) Richard Socher, Danqi Chen, Christopher D Manning, and Andrew Ng. 2013.  Reasoning with neural tensor networks for knowledge base completion.  _Advances in neural information processing systems_ 26 (2013).
  * Song et al. (2014) Xuan Song, Quanshi Zhang, Yoshihide Sekimoto, and Ryosuke Shibasaki. 2014.  Prediction of human emergency behavior and their mobility following large-scale disaster. In _Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining_. 5–14.
  * Tishby et al. (2000) Naftali Tishby, Fernando C Pereira, and William Bialek. 2000\.  The information bottleneck method.  _arXiv preprint physics/0004057_ (2000).
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin. 2017.  Attention is all you need.  _Advances in neural information processing systems_ 30 (2017).
  * Wang et al. (2022) Zhaonan Wang, Renhe Jiang, Hao Xue, Flora D Salim, Xuan Song, and Ryosuke Shibasaki. 2022\.  Event-Aware Multimodal Mobility Nowcasting. In _Proceedings of the AAAI Conference on Artificial Intelligence_ , Vol. 36. 4228–4236.
  * Widmer and Kubat (1996) Gerhard Widmer and Miroslav Kubat. 1996.  Learning in the presence of concept drift and hidden contexts.  _Machine learning_ 23, 1 (1996), 69–101.
  * Wu et al. (2020) Zonghan Wu, Shirui Pan, Guodong Long, Jing Jiang, Xiaojun Chang, and Chengqi Zhang. 2020\.  Connecting the dots: Multivariate time series forecasting with graph neural networks. In _Proceedings of the 26th ACM SIGKDD international conference on knowledge discovery & data mining_. 753–763.
  * Wu et al. (2019) Zonghan Wu, Shirui Pan, Guodong Long, Jing Jiang, and Chengqi Zhang. 2019.  Graph wavenet for deep spatial-temporal graph modeling. In _Proceedings of the 28th International Joint Conference on Artificial Intelligence_. 1907–1913.
  * Xu et al. (2020) Mingxing Xu, Wenrui Dai, Chunmiao Liu, Xing Gao, Weiyao Lin, Guo-Jun Qi, and Hongkai Xiong. 2020.  Spatial-temporal transformer networks for traffic flow forecasting.  _arXiv preprint arXiv:2001.02908_ (2020).
  * Yao et al. (2019) Huaxiu Yao, Xianfeng Tang, Hua Wei, Guanjie Zheng, and Zhenhui Li. 2019.  Revisiting spatial-temporal similarity: A deep learning framework for traffic prediction. In _Proceedings of the AAAI conference on artificial intelligence_ , Vol. 33. 5668–5675.
  * Yao et al. (2018) Huaxiu Yao, Fei Wu, Jintao Ke, Xianfeng Tang, Yitian Jia, Siyu Lu, Pinghua Gong, Jieping Ye, and Zhenhui Li. 2018.  Deep multi-view spatial-temporal network for taxi demand prediction. In _Proceedings of the AAAI conference on artificial intelligence_ , Vol. 32.
  * You et al. (2021) Xiaoyu You, Mi Zhang, Daizong Ding, Fuli Feng, and Yuanmin Huang. 2021.  Learning to learn the future: Modeling concept drifts in time series prediction. In _Proceedings of the 30th ACM International Conference on Information & Knowledge Management_. 2434–2443.
  * Yu et al. (2018) Bing Yu, Haoteng Yin, and Zhanxing Zhu. 2018.  Spatio-temporal graph convolutional networks: a deep learning framework for traffic forecasting. In _Proceedings of the 27th International Joint Conference on Artificial Intelligence_. 3634–3640.
  * Zhang et al. (2017) Junbo Zhang, Yu Zheng, and Dekang Qi. 2017.  Deep spatio-temporal residual networks for citywide crowd flows prediction. In _Thirty-first AAAI conference on artificial intelligence_.
  * Zhao et al. (2019) Ling Zhao, Yujiao Song, Chao Zhang, Yu Liu, Pu Wang, Tao Lin, Min Deng, and Haifeng Li. 2019.  T-gcn: A temporal graph convolutional network for traffic prediction.  _IEEE Transactions on Intelligent Transportation Systems_ 21, 9 (2019), 3848–3858.
  * Zhao et al. (2020) Peng Zhao, Le-Wen Cai, and Zhi-Hua Zhou. 2020.  Handling concept drift via model reuse.  _Machine Learning_ 109, 3 (2020), 533–568.
  * Zheng et al. (2020) Chuanpan Zheng, Xiaoliang Fan, Cheng Wang, and Jianzhong Qi. 2020\.  Gman: A graph multi-attention network for traffic prediction. In _Proceedings of the AAAI conference on artificial intelligence_ , Vol. 34. 1234–1241.


K. Chen Y. Ji, M. Möller and C. Vuik

*Corresponding author.

# A short report on the preconditioned Anderson acceleration

Kewang Chen  Ye Ji  Matthias Möller  Cornelis Vuik  College of Mathematics and
Statistics, Nanjing University of Information Science and Technology, Nanjing,
210044, China  School of Mathematical Sciences, Dalian University of
Technology, Dalian, 116024, China  Delft Institute of Applied Mathematics,
Delft University of Technology, Delft, 2628CD, the Netherlands
<EMAIL_ADDRESS>
###### Abstract

[Summary] In this report, we present a versatile and efficient preconditioned
Anderson acceleration (PAA) method for fixed-point iterations. The proposed
framework offers flexibility in balancing convergence rates (linear, super-
linear, or quadratic) and computational costs related to the Jacobian matrix.
Our approach recovers various fixed-point iteration techniques, including
Picard, Newton, and quasi-Newton iterations. The PAA method can be interpreted
as employing Anderson acceleration (AA) as its own preconditioner or as an
accelerator for quasi-Newton methods when their convergence is insufficient.
Adaptable to a wide range of problems with differing degrees of nonlinearity
and complexity, the method achieves improved convergence rates and robustness
by incorporating suitable preconditioners. We test multiple preconditioning
strategies on various problems and investigate a delayed update strategy for
preconditioners to further reduce the computational costs.

###### keywords:

Anderson acceleration, preconditioning, fixed-point iteration, nonlinear
system

††articletype: Report

##  1 Introduction

In 1962, Anderson 1, 2 developed an extrapolation algorithm to accelerate the
convergence of nonlinear fixed-point iterations. Since then, this has achieved
remarkable success and widespread use in various fields, such as computational
chemistry (known as Pulay mixing) and electronic structure computations
(referred to as Anderson mixing). Within the applied mathematics community,
this technique is now recognized as Anderson acceleration (AA). Unlike Picard
iterations, which utilize only one previous iterate, the Anderson Acceleration
method $AA(m)$ linearly combines a set of $m$ previous iterates to approximate
the minimization of the linearized fixed-point residual. Anderson acceleration
methods are regarded as “essentially equivalent" to nonlinear GMRES methods 3,
4, 5 and the direct inversion in the iterative subspace method 6, 7. They also
share similarities with methods based on quasi-Newton updating 8, 9. Recently,
some authors also investigated the non-stationary AA with dynamic window sizes
and damping factors. For example, Evans et al. 10 proposed a heuristic
strategy for selecting damping factors based on the gain at each iteration,
while Pollock and Rebholz 11 developed a strategy to dynamically alternate
window sizes. More recently, inspired by the hybrid linear solver GMRES
Recursive (GMRESR) 12, 13, Chen and Vuik 14 introduced a composite Anderson
acceleration method incorporating two window sizes.

AA has been widely observed to accelerate convergence, convergence analysis,
however, has only been reported more recently. In 2015, Toth and Kelley 15
initially demonstrated that the stationary version of AA without damping is
locally $r$-linearly convergent if the fixed-point map is a contraction.
Subsequently, Evans et al. 10 extended this result to AA with damping factors.
Pollock et al. 16 recently applied AA to the Picard iteration for solving
steady incompressible Navier–Stokes equations and showed that acceleration
improved the convergence rate of the Picard iteration. Sterck and He 17
further extended this result to a more general fixed-point iteration of the
form $\mathbf{x}=\mathbf{g}(\mathbf{x})$, given knowledge of the spectrum of
$\mathbf{g}^{\prime}(\mathbf{x})$ at the fixed-point $\mathbf{x}^{*}$.
Additionally, Wang et al. 18 investigated the asymptotic linear convergence
speed of sAA when applied to the Alternating Direction Method of Multipliers
(ADMM) method. Global convergence and sharper local convergence results of AA
remain active research topics. For more related results about AA and its
applications, we refer interested readers to the papers 19, 20, 21, 22, 23,
24, 25, 26, 27 and the references therein.

Despite the extensive utilization and investigation of AA over several
decades, most authors have primarily concentrated on analyzing its classical
version applied to Picard iterations. Nonetheless, the potential of
preconditioned Anderson acceleration (PAA) and its applicability to other
iterative methods have remained largely unexplored, with only a few reported
results 28, 29, 30, 31. For instance, Herbst and Levitt 29 proposed using
effective preconditioners $\mathbf{P}$ for solving the nonlinear equations
$\mathbf{f(\mathbf{x})}=\mathbf{0}$ through preconditioned fixed-point maps
$\mathbf{g(x)}=\mathbf{x}-\mathbf{P}\mathbf{f(\mathbf{x})}$. In the specific
case of Newton’s method, where $\mathbf{P}$ approximates the inverse of
$\mathbf{f^{\prime}(\mathbf{x^{*}})}$, Pollock and Schwartz28 introduced and
investigated the Anderson accelerated Newton method, which employs Anderson
acceleration to expedite Newton’s iterations. More recently, Ji et al.31
applied the preconditioned Anderson acceleration framework to a PDE-based
elliptic parameterization technique in isogeometric analysis. However, the
relationship between preconditioned Anderson acceleration and other iterative
methods remains unclear, and the convergence analysis of PAA is currently
missing.

This study presents a comprehensive investigation and analysis of a versatile
preconditioned Anderson acceleration (PAA) framework. Initially, we
demonstrate that PAA can encompass various fixed-point iteration techniques,
including Picard, Newton, and quasi-Newton iterations. The approach of PAA can
be viewed as employing Anderson acceleration both as its own preconditioner
and as an accelerator for quasi-Newton methods when their convergence is
insufficient. Subsequently, under suitable assumptions regarding non-singular
preconditioning matrices, we establish a convergence result for the proposed
PAA method. The flexibility of the PAA framework allows for striking a balance
between the convergence rate (linear, super-linear, or quadratic) and the
computational cost associated with approximating the Jacobian matrix. This
adaptability renders the method suitable for a wide array of problems
characterized by different levels of nonlinearity. By incorporating an
appropriate preconditioner into the Anderson acceleration method, we achieve
superior convergence rates and increased robustness. Lastly, we conduct a
series of practical preconditioning strategies on various problems to enhance
numerical stability and convergence speed. Through extensive numerical
experiments, we demonstrate the effectiveness of the proposed update strategy.

##  2 Preconditioned Anderson acceleration

Recall that the Picard iteration:

$\mathbf{x}_{k+1}=\mathbf{g}(\mathbf{x}_{k})=\mathbf{x}_{k}-\mathbf{f}(\mathbf{x}_{k})$
(1)

can only converges linearly. However, when Anderson acceleration is used to
accelerates this Picard iteration, a super-linear convergence rate could be
achieved. Therefore, to further improve the numerical stability and the
convergence speed, Anderson acceleration can also be used to accelerate the
quasi-Newton method and even the Newton’s iteration method.

Specifically, we introduce non-singular matrices $\mathbf{M}_{k}$ as the
preconditioners at iteration $k$. With this, we define the following
preconditioned fixed-point iteration scheme, which is also known as the quasi-
Newton iteration scheme:

$\mathbf{x}_{k+1}=\mathbf{x}_{k}-\mathbf{M}_{k}^{-1}\mathbf{f}(\mathbf{x}_{k}).$
(2)

Next, we utilize the Anderson acceleration (AA) algorithm to expedite the
convergence of this preconditioned iteration (2). This approach can be seen as
using AA as a preconditioner for AA itself, or as a means to accelerate quasi-
Newton methods in cases where the quasi-Newton approach alone proves
insufficient.

Denoting the most recent iterates by
$\mathbf{x}_{k-m},\cdots,\mathbf{x}_{k}\in\mathbb{R}^{n}$ and the
corresponding outputs
$\mathbf{w}_{k-m}=-\mathbf{M}_{k}^{-1}\mathbf{f}(\mathbf{x}_{k-m}),\cdots,\mathbf{w}_{k}=-\mathbf{M}_{k}^{-1}\mathbf{f}(\mathbf{x}_{k})\in\mathbb{R}^{n}$.
Then the Anderson acceleration takes the latest $m+1$ steps into account:

$\mathbf{x}_{k+1}^{AA}=\bar{\mathbf{x}}_{k}+\beta\bar{\mathbf{w}}_{k},$ (3)

where the damping parameter $\beta\in(0,1]$ is a real number,
$\bar{\mathbf{x}}_{k}$ and $\bar{\mathbf{f}}_{k}$ denotes the weighted average
of the previous iterates and residuals, respectively. Specifically,

$\bar{\mathbf{x}}_{k}=\mathbf{x}_{k}-\sum_{j=1}^{m}\gamma_{j}^{k}(\mathbf{x}_{k-m+j}-\mathbf{x}_{k-m+j-1})=\mathbf{x}_{k}-\mathbf{E}_{k}\Gamma_{k},$
(4)

$\bar{\mathbf{w}}_{k}=-\mathbf{M}_{k}^{-1}\mathbf{f}_{k}-\sum_{j=1}^{m}\gamma_{j}^{k}(\mathbf{w}_{k-m+j}-\mathbf{w}_{k-m+j-1})=\mathbf{w}_{k}-\mathbf{W}_{k}\Gamma_{k},$
(5)

where $\mathbf{E}_{k}$ and $\mathbf{W}_{k}$ are given by

$\displaystyle\mathbf{E}_{k}$
$\displaystyle=\left[(\mathbf{x}_{k-m+1}-\mathbf{x}_{k-m})\ \
(\mathbf{x}_{k-m+2}-\mathbf{x}_{k-m+1})\ \ \cdots\ \
(\mathbf{x}_{k}-\mathbf{x}_{k-1})\right],$ (6) $\displaystyle\mathbf{W}_{k}$
$\displaystyle=\left[(\mathbf{w}_{k-m+1}-\mathbf{w}_{k-m})\ \
(\mathbf{w}_{k-m+2}-\mathbf{w}_{k-m+1})\ \ \cdots\ \
(\mathbf{w}_{k}-\mathbf{w}_{k-1})\right].$

The scalars $\Gamma_{k}=[\gamma_{1},\cdots,\gamma_{m}]^{T}$ are chosen so as
to minimize the $l_{2}$ norm of the residual

$\Gamma_{k}=\rm{argmin}_{\Gamma\in\mathcal{R}^{m}}\|\mathbf{w}_{k}-\mathbf{W}_{k}\Gamma\|^{2}_{2}.$
(7)

Then, we obtain

$\displaystyle\mathbf{x}_{k+1}^{AA}$
$\displaystyle=\bar{\mathbf{x}}_{k}+\beta\bar{\mathbf{w}}_{k}$ (8)
$\displaystyle=\mathbf{x}_{k}+\beta\mathbf{M}_{k}^{-1}\mathbf{f}(\mathbf{x}_{k})-\left(\mathbf{E}_{k}+\beta\mathbf{W}_{k}\right)\left[\mathbf{W}_{k}^{\mathrm{T}}\mathbf{W}_{k}\right]^{-1}\mathbf{W}_{k}^{\mathrm{T}}\mathbf{M}_{k}^{-1}\mathbf{f}(\mathbf{x}_{k}).$

If we set

$\mathbf{C}_{k}=\frac{1}{\beta}\left[\left(\mathbf{E}_{k}+\beta\mathbf{W}_{k}\right)\left[\mathbf{W}_{k}^{\mathrm{T}}\mathbf{W}_{k}\right]^{-1}\mathbf{W}_{k}^{\mathrm{T}}-\beta\mathbf{I}\right],$
(9)

the update of the preconditioned Anderson acceleration reads

$\mathbf{x}_{k+1}^{AA}=\mathbf{x}_{k}-\beta\mathbf{C}_{k}\mathbf{M}_{k}^{-1}\mathbf{f}(\mathbf{x}_{k}).$
(10)

Next, we give our motivations on how to choose those $\mathbf{M}_{k}$ and
introduce several practical preconditioners.

###  2.1 Motivation on choosing suitable preconditioner

To determine suitable preconditioners for our proposed preconditioned Anderson
acceleration method, we first examine the convergence rate of the
unaccelerated iteration scheme (2). The convergence results of iteration (2)
are summarized in the following theorem.

######  Theorem 2.1.

For the nonlinear system (LABEL:eq:4601), let $\mathbf{x}^{*}$ be a solution
of $\mathbf{f}(\mathbf{x})=\mathbf{0}$ and suppose that
$\mathbf{f}(\mathbf{x})$ is continuously differentiable and second-order
derivative exists, then the preconditioned fixed-point iteration scheme in (2)

$\mathbf{x}_{k+1}=\mathbf{x}_{k}-\mathbf{M}_{k}^{-1}\mathbf{f}(\mathbf{x}_{k})$

converges if

$\left\|\mathbf{I}-\mathbf{M}_{k}^{-1}\bm{\mathcal{J}}_{k}\right\|<1,$ (11)

where $\bm{\mathcal{J}}_{k}$ is the $n\times n$ Jacobian matrix of
$\mathbf{f}(\mathbf{x})$ at iteration $k$. Moreover, if we set the local error
$\mathbf{e}_{k}=\mathbf{x}^{*}-\mathbf{x}_{k}$, then the relationship between
the error in any two consecutive iterates can be shown to be

$\mathbf{e}_{k+1}=\left(\mathbf{I}-\mathbf{M}_{k}^{-1}\bm{\mathcal{J}}_{k}\right)\mathbf{e}_{k}-\frac{1}{2}\mathbf{M}_{k}^{-1}\mathbf{E}_{k}^{\mathrm{T}}\mathbf{H}_{k}\mathbf{E}_{k}+\mathcal{O}(\|\mathbf{e}_{k}\|^{3}).$
(12)

Here,
$\mathbf{E}_{k}=diag\\{\mathbf{e}_{k}^{\mathrm{T}},\ldots,\mathbf{e}_{k}^{\mathrm{T}}\\}$
is a diagonal matrix composed of the transposed residual vector
$\mathbf{e}_{k}^{\mathrm{T}}$, and

$\mathbf{H}_{k}=diag\\{\mathbf{H}^{(1)}(\mathbf{x}),\mathbf{H}^{(2)}(\mathbf{x}),\ldots,\mathbf{H}^{(n)}(\mathbf{x})\\}$

is a block-diagonal matrix composed of $\mathbf{H}^{(i)}(\mathbf{x})$ which is
the Hessian matrix of the $i$-th nonlinear equation $\mathbf{f}^{(i)}(x)$
($i=1,2,\ldots,n$) at iteration $k$.

######  Proof 2.2.

Let us assume that $\mathbf{x}^{*}$ is the solution to the nonlinear system
$\mathbf{f}(\mathbf{x})=\mathbf{0}$, and $\mathbf{x}_{k}$ is an estimate for
$\mathbf{x}^{*}$ at the $k$-th iteration, such that
$\|\mathbf{e}_{k}\|=\|\mathbf{x}^{*}-\mathbf{x}_{k}\|\ll 1$. Based on Taylor’s
expansion theorem, we have:

$\mathbf{0}=\mathbf{f}(\mathbf{x}^{*})=\mathbf{f}(\mathbf{x}_{k})+\bm{\mathcal{J}}_{k}\mathbf{e}_{k}+\frac{1}{2}\mathbf{E}_{k}^{\mathrm{T}}\mathbf{H}_{k}\mathbf{E}_{k}+\mathcal{O}(\|\mathbf{e}_{k}\|^{3}),$
(13)

Utilizing the preconditioned fixed-point iteration (2), we obtain:

$\mathbf{f}({\mathbf{x}_{k}})=\mathbf{M}_{k}\left(\mathbf{x}_{k}-\mathbf{x}_{k+1}\right)=\mathbf{M}_{k}\left((\mathbf{x}^{*}-\mathbf{x}_{k+1})-(\mathbf{x}^{*}-\mathbf{x}_{k})\right)=\mathbf{M}_{k}\left(\mathbf{e}_{k+1}-\mathbf{e}_{k}\right).$
(14)

Substituting the above formula (14) into (13) results in

$\displaystyle\mathbf{0}$
$\displaystyle=\mathbf{f}(\mathbf{x}_{k})+\bm{\mathcal{J}}_{k}\mathbf{e}_{k}+\frac{1}{2}\mathbf{E}_{k}^{\mathrm{T}}\mathbf{H}_{k}\mathbf{E}_{k}+\mathcal{O}(\|\mathbf{e}_{k}\|^{3})$
(15)
$\displaystyle=\mathbf{M}_{k}\left(\mathbf{e}_{k+1}-\mathbf{e}_{k}\right)+\bm{\mathcal{J}}_{k}\mathbf{e}_{k}+\frac{1}{2}\mathbf{E}_{k}^{\mathrm{T}}\mathbf{H}_{k}\mathbf{E}_{k}+\mathcal{O}(\|\mathbf{e}_{k}\|^{3}).$

Multiplying $\mathbf{M}_{k}^{-1}$ on both sides of equation (15), we have

$\displaystyle\mathbf{0}$
$\displaystyle=\left(\mathbf{e}_{k+1}-\mathbf{e}_{k}\right)+\mathbf{M}_{k}^{-1}\bm{\mathcal{J}}_{k}\mathbf{e}_{k}+\frac{1}{2}\mathbf{M}_{k}^{-1}\mathbf{E}_{k}^{\mathrm{T}}\mathbf{H}_{k}\mathbf{E}_{k}+\mathcal{O}(\|\mathbf{e}_{k}\|^{3})$
(16)
$\displaystyle=\mathbf{e}_{k+1}-\left(\mathbf{I}-\mathbf{M}_{k}^{-1}\bm{\mathcal{J}}_{k}\right)\mathbf{e}_{k}+\frac{1}{2}\mathbf{M}_{k}^{-1}\mathbf{E}_{k}^{\mathrm{T}}\mathbf{H}_{k}\mathbf{E}_{k}+\mathcal{O}(\|\mathbf{e}_{k}\|^{3}),$

where $\mathbf{I}$ is the $n\times n$ identity matrix.

Finally, we arrive at

$\mathbf{e}_{k+1}=\left(\mathbf{I}-\mathbf{M}_{k}^{-1}\bm{\mathcal{J}}_{k}\right)\mathbf{e}_{k}-\frac{1}{2}\mathbf{M}_{k}^{-1}\mathbf{E}_{k}^{\mathrm{T}}\mathbf{H}_{k}\mathbf{E}_{k}+\mathcal{O}(\|\mathbf{e}_{k}\|^{3}).$

In particular, convergence of the preconditioned fixed-point iteration
requires

$\left\|\mathbf{I}-\mathbf{M}_{k}^{-1}\bm{\bm{\mathcal{J}}}_{k}\right\|<1.$

######  Remark 2.3.

For the linear problem,

$\mathbf{f}(\mathbf{x})=\mathbf{Ax-b}=\mathbf{0},\
\mathbf{x}\in\mathbb{R}^{n},$ (17)

suppose the matrix $\mathbf{A}$ can be split as

$\mathbf{A=D+R},$ (18)

where all the off-diagonal components of the matrix $\mathbf{D}$ and the
diagonal components of $\mathbf{R}$ are identically zero. By using this
decomposition, the preconditioned fixed-point iteration (2) converges if

$\left\|\mathbf{I}-\mathbf{D}^{-1}\mathbf{A}\right\|<1,$ (19)

which recovers the results from Pratapa and Suryanarayana32.

######  Remark 2.4.

If we set the preconditioner

$\mathbf{M}_{k}=\bm{\mathcal{J}}_{k},$ (20)

then we obtain the Newton iteration. In this case, (12) boils down to

$\displaystyle\mathbf{e}_{k+1}$
$\displaystyle=\left(\mathbf{I}-\mathbf{M}_{k}^{-1}\bm{\mathcal{J}}_{k}\right)\mathbf{e}_{k}-\frac{1}{2}\mathbf{M}_{k}^{-1}\mathbf{E}_{k}^{\mathrm{T}}\mathbf{H}_{k}\mathbf{E}_{k}+\mathcal{O}(\|\mathbf{e}_{k}\|^{3})$
(21)
$\displaystyle=\left(\mathbf{I}-\bm{\mathcal{J}}_{k}^{-1}\bm{\mathcal{J}}_{k}\right)\mathbf{e}_{k}-\frac{1}{2}\mathbf{M}_{k}^{-1}\mathbf{E}_{k}^{\mathrm{T}}\mathbf{H}_{k}\mathbf{E}_{k}+\mathcal{O}(\|\mathbf{e}_{k}\|^{3})$
$\displaystyle=-\frac{1}{2}\bm{\mathcal{J}}_{k}^{-1}\mathbf{E}_{k}^{\mathrm{T}}\mathbf{H}_{k}\mathbf{E}_{k}+\mathcal{O}(\|\mathbf{e}_{k}\|^{3})$

therefore, we have

$\|\mathbf{e}_{k+1}\|\propto\|\mathbf{e}_{k}\|^{2}$ (22)

which means the Newton’s method achieves quadratic convergence as we expected.

This implies that, when commencing with a suitably accurate initial guess and
updating the preconditioner $\mathbf{M}_{k}$ at each iteration, our
preconditioned AA scheme (2) can attain a quadratic convergence rate, aligning
with the standard Newton iteration. Nevertheless, we encounter a similar
problem as faced by Newton-type solvers: the computationally intensive task of
frequently updating the full Jacobian, which contradicts our initial intention
of employing AA to reduce computational overhead. To address this concern, we
propose several more practical preconditioning strategies aimed at effectively
solving the involved nonlinear systems.

Building upon the above discussions, we summarize several typical
preconditioning strategies for $\mathbf{M}_{k}$ in Table 1.

Table 1: The selection of preconditioners $\mathbf{M}_{k}$ and its pros and cons Preconditioners $\mathbf{M}_{k}$ |  Pros and Cons
---|---
1\. Constant preconditioner, i.e., $\mathbf{M}_{k}=\alpha\mathbf{I}$, where $\alpha$ is a constant value. |  This preconditioning strategy is straightforward and computationally economical. However, it requires manual tuning of the $\alpha$ value. If $\alpha$ is set to 1, it reverts to the case of no preconditioning. It’s worth noting that as per (12), if the preconditioners $\mathbf{M}_{k}$ significantly deviate from the Jacobian, convergence may become limited to linear or even diverge in certain cases.
2\. Full Jacobian preconditioner, i.e., $\mathbf{M}_{k}=\bm{\mathcal{J}}_{k}$. |  Essentially, this preconditioning strategy uses AA to accelerate Newton’s method 28. However, it comes with computational expenses since we need to update this preconditioner at each iteration.
3\. Diagonal Jacobian or Block-diagonal Jacobian preconditioner, i.e., $\mathbf{M}_{k}=\rm{diag}(\bm{\mathcal{J}}_{k})\ or\ \mathbf{M}_{k}=\rm{diagBlock}(\bm{\mathcal{J}}_{k})$. |  To strike a balance between a “good preconditioner" and “computational overhead", the (block) diagonal Jacobian preconditioner can be a promising candidate due to its relatively inexpensive computation.
4\. Other approximations of Jacobian as preconditioners. |  For example, we can use rank one correction and rank two correction of the inverse of Jacobian matrix in L-BFGS 33 as a preconditioner. Besides, we could also use AA itself as a preconditioner for AA14. These preconditioners are computationally cheaper and they can approximate the inverse of Jacobian matrix directly.
5\. Updating the preconditioners $\mathbf{M}$ every $N_{update}$ steps instead of updating them each iteration. |  Unlike Newton-type solvers and gradient-based methods, the AA iteration scheme involves a linear combination of several previous iterates. This characteristic allows us to avoid frequent updates of preconditioners, which can be computationally expensive.

###  2.2 Overview of the proposed algorithm and its relations with other
methods

In this section, we present the details of the proposed algorithm and explore
its connections with other methods. The general preconditioned Anderson
acceleration scheme, incorporating a dynamic preconditioning strategy, is
summarized in the following Algorithm 1.

Input: $\mathbf{x}_{0}$: Initial guess;

$\mathbf{f}(\mathbf{x})$: Nonlinear system;

$m$: Window size used in AA;

$N_{max}$: Maximum iterations;

$tol$: Tolerance of residual norm for convergence test;

$N_{update}$: Update the preconditioner every $N_{update}$ steps.

Output: $\mathbf{x}_{k+1}$: Numerical solution

for  _$k=0\mathrm{;}$ $k\leq N_{max}\mathrm{;}$ $++k$_ do

Compute residual vector $\mathbf{f}_{k}=\mathbf{f}(\mathbf{x}_{k})$;

// Check the termination criteria

if  _ $\|\mathbf{f}_{k}\|<tol$_ then

break;

end if

// Update preconditioner

if  _$k$ is evenly divisible by $N_{update}$_ then

Compute preconditioning matrix $\mathbf{M_{k}}$;  // Table 1

end if

// Anderson acceleration

Compute $m_{k}=\min\\{m,k\\}$;

Compute preconditioned residual vector $\mathcal{F}_{k}$ by solving
$\mathbf{M_{k}}\mathcal{F}_{k}=-\mathbf{f}_{k}$;

Determine
$\bm{\alpha}^{(k)}=\left(\alpha_{0}^{(k)},\cdots,\alpha_{m_{k}}^{(k)}\right)^{\mathrm{T}}$
that solves
$\displaystyle\mathop{\arg\min}\limits_{\bm{\alpha}=(\alpha_{0},\cdots,\alpha_{m_{k}})^{\mathrm{T}}}$
$\displaystyle\left\|\sum_{i=0}^{m_{k}}\mathcal{F}_{k-m_{k}+i}\bm{\alpha}\right\|^{2}_{2}$
s.t. $\displaystyle\sum_{i=0}^{m_{k}}\alpha_{i}=1;$

Update
$\displaystyle\mathbf{x}_{k+1}=\sum_{i=0}^{m_{k}}\alpha_{i}^{(k)}\left(\mathbf{x}_{k-m_{k}+i}+\mathcal{F}_{k-m_{k}+i}\right)$;

end for

Algorithm 1 Preconditioned Anderson acceleration: PAA(m)

######  Remark 2.5.

The proposed preconditioned Anderson acceleration $(\mathbf{PAA}(m))$ method
can be viewed as a nonlinear acceleration method applied to a general
iteration scheme:

$\mathbf{x}_{k+1}=\mathbf{U}(\mathbf{x}_{k}),$ (23)

where $\mathbf{U}(\mathbf{x})$ is an update function. The $\mathbf{PAA}(m)$
method can accelerate various fixed-point iteration methods, including the
Picard iteration, the Newton iteration, and the quasi-Newton iteration.
Specifically, the $\mathbf{PAA}(m)$ method accelerates the Picard iteration
when $\mathbf{M}_{k}=\mathbf{I}$, i.e.,
$\mathbf{U}(\mathbf{x}_{k})=\mathbf{x}_{k}-\mathbf{f}(\mathbf{x}_{k})$ (see
Anderson1), accelerates the Newton iteration when
$\mathbf{M}_{k}=\bm{\mathcal{J}}_{k}$, i.e.,
$\mathbf{U}(\mathbf{x}_{k})=\mathbf{x}_{k}-[\bm{\mathcal{J}}_{k}]^{-1}\mathbf{f}(\mathbf{x}_{k})$
(see Pollock and Schwartz 28), and accelerates the quasi-Newton iteration when
the preconditioner matrix $\mathbf{M}_{k}$ lies between these two matrices and
$\mathbf{U}(\mathbf{x}_{k})=\mathbf{x}_{k}-\mathbf{M}_{k}^{-1}\mathbf{f}(\mathbf{x}_{k})$.
Moreover, the $\mathbf{PAA}(m)$ method can recover all the above unaccelerated
iteration methods by setting the window size $m=0$ and using different
$\mathbf{M}_{k}$ in $\mathbf{PAA}(m)$.

The key advantage of the $\mathbf{PAA}(m)$ method lies in its ability to
easily control and manipulate the trade-off between convergence rate (linear,
superlinear, or quadratic) and the cost of computing or approximating the
Jacobian matrix (identity matrix, diagonal Jacobian matrix, full Jacobian
matrix, or other approximation Jacobian matrix). This control can be achieved
by tuning the window size $m$, the preconditioners $\mathbf{M}_{k}$, and the
frequency of updating these preconditioners $\mathbf{N}_{update}$ in
$\mathbf{PAA}(m)$.

######  Remark 2.6.

Other aspects to further improve the performance of AA: Using damping factor
or dynamical damping factors $\beta$, see Evans et al.10 and Chen and Vuik26;
Using filtering for Anderson acceleration, see Haelterman et al.34 and Pollock
and Rebholz35; Using multiple window sizes, see Pollock and Rebholz11.

##  3 Numerical experiments

###  3.1 Benchmark results

The problems considered here are the systems of nonlinear equations from the
standard test set in literature. Throughout this section, each method was
terminated on the residual reduced below the tolerance
$\|\mathbf{f}(\mathbf{x_{k}})\|<10^{-10}$. All these experiments are done in
the MATLAB 2021b environment.

######  Example 3.1.

[2D nonlinear polynomial systems36] Apply preconditioned AA to solve the
following $\mathbf{f}(\mathbf{x})=\mathbf{0}$ with

$\left\\{\begin{aligned} f_{1}(x_{1},x_{2})&=(x_{1}-1)+(x_{2}-3)^{2},\\\
f_{2}(x_{1},x_{2})&=\epsilon(x_{2}-3)+\frac{3}{2}(x_{1}-1)(x_{2}-3)+(x_{2}-3)^{2}+(x_{2}-3)^{3}.\end{aligned}\right.$
(24)

In this simple example, we know from Kelley and Suresh36 that there are two
distinct roots $\mathbf{x}_{+}=(1,3)^{\mathrm{T}}$ and
$\mathbf{x}_{-}=(1-\eta^{2},3+\eta)^{\mathrm{T}}$ with
$\eta=1-\sqrt{1+2\epsilon}$ for $\epsilon>0$. When $\epsilon=0$, the two roots
collapse into the same solution $\mathbf{x}=(1,3)^{\mathrm{T}}$. We apply
preconditioned AA with different preconditioning strategies to solve this
problem for both $\epsilon=0$ and $\epsilon=10^{-6}$. For each method, we use
$5$ random initial guesses $(x_{0},y_{0})\in[-1,3]\times[1,5]$ in our
numerical experiments.

From Fig. 1, we observe that the original Anderson acceleration method
$AA(1)$, without preconditioning, exhibits poor performance and may even
diverge for certain initial guesses in both cases. A simple constant
preconditioning strategy, $\mathbf{M}_{k}=\alpha\mathbf{I}$ with $\alpha=0.1$,
improves the convergence rate slightly, but it remains slow. On the other
hand, when better preconditioners, such as the diagonal Jacobian and the full
Jacobian preconditioning strategies, are employed, the performance is
significantly improved. However, it is essential to consider the computational
costs of these preconditioners when comparing their performance for solving
large systems of nonlinear equations.

((a)) $\epsilon=0$

((b)) $\epsilon=10^{-6}$

Figure 1: Using Preconditioned AA to solve a 2D nonlinear polynomial system
with 5 random initial guesses.

######  Example 3.2.

[Trigonometric nonlinear equations, n=5,50,50037] Apply preconditioned AA to
solve the following larger nonlinear trigonometric systems. Let

$h_{i}(x)=n-\sum_{j=1}^{n}\cos(x_{j})+i[1-\cos(x_{i})]-\sin(x_{i}),\
i=1,2,\cdots,n,$ (25)

we create a manufacturing fixed-point
$\mathbf{x}^{*}=(\frac{\pi}{4},\cdots,\frac{\pi}{4})^{\mathrm{T}}$ of the
following system

$n-\sum_{j=1}^{n}\cos(x_{j})+i(1-\cos(x_{i}))-\sin(x_{i})=h_{i}(\frac{\pi}{4}),\
i=1,2,\cdots,n.$ (26)

Then, we define the nonlinear equations as:

$f_{i}(x)=\left(h_{i}(x)-h_{i}(\frac{\pi}{4})\right),\ i=1,2,\cdots,n.$ (27)

with real solution
$\mathbf{x}^{*}=(\frac{\pi}{4},\cdots,\frac{\pi}{4})^{\mathrm{T}}$.

In this numerical test, we adopt an example from Moré et al.37 such that a
manufacturing solution is created. The initial guesses
$\mathbf{x}_{0}\in\mathbb{R}^{n}$ used here are five random vectors close to
$\mathbf{x}^{*}$ such that
$\mathbf{x}_{0}\in(\pi/4-0.05,\pi/4+0.05)\times\cdots\times(\pi/4-0.05,\pi/4+0.05)$.
We apply preconditioned AA with different preconditioning strategies to solve
equations (27) for $n=5,50,500$. As the dimension $n$ increases, we use larger
window sizes for Anderson acceleration. The results of our experiments are
illustrated in Fig. 2.

When $n=5$, we see from Fig. 2(a) that the straightforward constant
preconditioning approach using $\mathbf{M}_{k}=\alpha\mathbf{I}$ with
$\alpha=0.1$ and $\alpha=0.01$, yields poorer results than the original
$AA(3)$ method. However, the diagonal Jacobian preconditioning and the full
Jacobian strategies performs much better than original $AA(3)$ as we expected.
This indicates that the simple constant preconditioner may not always works
and a good preconditioner is very important in solving nonlinear equations.

When $n=50$ and $n=500$, as we can see from Fig. 2(b) and Fig. 2(c) that the
original Anderson acceleration without preconditioning does not converge. The
simple constant preconditioning strategy performs slightly better but still
diverges or converges very slowly. Again, if we choose better preconditioners,
the diagonal Jacobian preconditioning and the full Jacobian strategies
converge very fast.

However, as $n$ grows larger, the cost of computing a full Jacobian
preconditioner appears to be progressively increasing, see Table. 2. For
example, when $n=500$, the Anderson acceleration with diagonal Jacobian
preconditioner only needs around 0.09 seconds to find the solution while the
full Jacobian preconditioning strategy (Anderson accelerates Newton’s
iteration) needs almost around 16 seconds to achieve the same accuracy.
Moreover, to further improve the time efficiency of preconditioned AA, we can
also update our preconditioners every $N_{update}$ steps instead calculating
it at each iteration. For instance, when employing Anderson acceleration with
the full Jacobian preconditioning strategy, setting the number of steps
between preconditioner updates to $N_{update}=2$ leads to superior performance
in comparison to computing preconditioners at every iteration, as demonstrated
by the results presented in Table 2.

((a)) n=5

((b)) n=50

((c)) n=500

Figure 2: Using preconditioned AA to solve trigonometric nonlinear equations. Table 2: Total average time (seconds) used for different preconditioned AA methods. Methods |  $n=5$ |  $n=50$ |  $n=500$
---|---|---|---
AA | 0.0040 | - | -
DiagJacobian preconditioner | 0.0030 | 0.0040 | 0.0909
FullJacobian preconditioner | 0.0038 | 0.0411 | 15.8453
FullJacobian preconditioner, $N_{update}=2$ | 0.0021 | 0.0322 | 11.7176

######  Example 3.3.

[The nonlinear Bratu problem5] The Bratu problem is a nonlinear PDE boundary
value problem as follows:

$\displaystyle\Delta u+\lambda\ e^{u}$ $\displaystyle=0,\ \ in\ \
\Omega=[0,1]\times[0,1],$ $\displaystyle u$ $\displaystyle=0,\ \ on\ \
\partial\Omega.$

This example is from Walker5, where a stationary Anderson acceleration with
window size $m=50$ is used to solve the Bratu problem. This problem has a long
history, we refer the reader to Glowinski et al.38 and the references in those
papers.

In this experiment, we employe a centered-difference discretization approach
to approximate the solution of the Bratu problem on two different grid sizes:
$32\times 32$ and $64\times 64$. The value of parameter $\lambda$ in the Bratu
problem was set to 6, and the initial guess for the solution was chosen as
$u_{0}=(1,1,\cdots,1)^{T}$ for all cases. It was observed that the Picard
iteration diverged when no preconditioning was applied. We also note here that
the preconditioning matrix of the linear part refers to the matrix $A$, where
$A$ is a matrix for the discrete Laplace operator. The numerical results are
shown in Fig. 3.

From Fig. 3, we observe that Anderson acceleration, utilizing a diagonal
Jacobian preconditioner (either for only the linear component or for the
complete nonlinear system), notably expedites the Picard iteration. Although
the variant utilizing the full Jacobian achieves even faster convergence,
computing the complete Jacobian incurs significant computational expense. It
necessitates more time to achieve comparable accuracy compared to the diagonal
preconditioning strategies. Consequently, deriving a "reasonably good and
affordable" approximation of the complete Jacobian is crucial when solving
nonlinear systems with high dimensions. The preconditioners utilized in
preconditioned Anderson acceleration (PAA) must be adjusted accordingly to
obtain optimal acceleration.

((a)) 32 $\times$ 32

((b)) 64 $\times$ 64

Figure 3: Using preconditioned AA to solve the nonlinear Bratu problem.

######  Example 3.4.

[The stationary nonlinear convection-diffusion problem14] Solve the following
2D nonlinear convection-diffusion equation in a square region:

$\epsilon(-u_{xx}-u_{yy})+(u_{x}+u_{y})+ku^{2}=f(x,y),\ \
(x,y)\in\Omega=[0,1]\times[0,1]$ (28)

with the source term

$f(x,y)=2\pi^{2}\sin(\pi x)\sin(\pi y)$ (29)

and zero boundary conditions: $u(x,y)=0$ on $\partial\Omega$.

In this numerical experiment, we use a centered-difference discretization on a
$32\times 32$ grid. We choose $\epsilon=0.1$ and $\epsilon=0.01$,
respectively. Those $\epsilon$ values indicate the competition between the
diffusion and convection effect. We take $k=3$ in the above problem and use
$u_{0}=(1,1,\cdots,1)^{T}$ as an initial approximate solution in all cases.
Again, the preconditioning matrix of the linear part refers to the matrix $A$,
where $A$ is a matrix for the discrete Laplace operator.

The results are shown in Fig. 4, which are consistent with that of solving the
Bratu problem. Particularly, when $\epsilon=0.01$, the linear part becomes
less important and the nonlinear term dominates. Fig.4(b) illustrates that
solely utilizing the linear portion of the diagonal Jacobian in a simple
Picard iteration leads to divergence. However, our preconditioned Anderson
acceleration method, utilizing the same preconditioner, achieves convergence.
Furthermore, implementing a superior preconditioner (i.e., the diagonal
Jacobian of the complete nonlinear equations) results in significantly faster
convergence, as anticipated.

((a)) $\epsilon=0.1$

((b)) $\epsilon=0.01$

Figure 4: Using preconditioned AA to solve the stationary nonlinear
convection-diffusion problem.

##  4 Conclusions and outlook

Our current study introduces a preconditioned Anderson acceleration scheme.
This method provides the flexibility to adjust the balance between convergence
rate (linear, super-linear, or quadratic) and the cost of computing or
approximating the Jacobian matrix. It can recover different fixed-point
iteration techniques such as the Picard, Newton, and quasi-Newton iterations.
The preconditioned Anderson acceleration method is versatile and can be
adapted to a wide range of problems with varying levels of nonlinearity. By
incorporating an appropriate preconditioner, improved convergence rates and
enhanced robustness can be achieved. We test several preconditioning
strategies, including a constant preconditioner, diagonal Jacobian or block-
diagonal Jacobian preconditioner, and full Jacobian preconditioner.
Additionally, we discuss a delayed update strategy of the preconditioner,
where it is updated every few steps to reduce the computational cost per
iteration.

Overall, our proposed framework is an effective tool for solving nonlinear
problems, enabling users to tailor the convergence behavior to their specific
needs while maintaining low computational costs.

## References

  * 1 Anderson DG. Iterative procedures for nonlinear integral equations. Journal of the ACM (JACM) 1965; 12(4): 547–560.
  * 2 Anderson DG. Comments on “Anderson acceleration, mixing and extrapolation”. Numerical Algorithms 2019; 80(1): 135–234.
  * 3 Carlson NN, Miller K. Design and application of a gradient-weighted moving finite element code I: in one dimension. SIAM Journal on Scientific Computing 1998; 19(3): 728–765.
  * 4 Oosterlee CW, Washio T. Krylov subspace acceleration of nonlinear multigrid with application to recirculating flows. SIAM Journal on Scientific Computing 2000; 21(5): 1670–1690.
  * 5 Walker HF, Ni P. Anderson acceleration for fixed-point iterations. SIAM Journal on Numerical Analysis 2011; 49(4): 1715–1735.
  * 6 Pulay P. Convergence acceleration of iterative sequences. The case of SCF iteration. Chemical Physics Letters 1980; 73(2): 393–398.
  * 7 Lin L, Yang C. Elliptic preconditioner for accelerating the self-consistent field iteration in Kohn–Sham density functional theory. SIAM Journal on Scientific Computing 2013; 35(5): S277–S298.
  * 8 Fang Hr, Saad Y. Two classes of multisecant methods for nonlinear acceleration. Numerical Linear Algebra with Applications 2009; 16(3): 197–221.
  * 9 Haelterman R, Degroote J, Van Heule D, Vierendeels J. On the similarities between the quasi-Newton inverse least squares method and GMRES. SIAM Journal on Numerical Analysis 2010; 47(6): 4660–4679.
  * 10 Evans C, Pollock S, Rebholz LG, Xiao M. A proof that Anderson acceleration improves the convergence rate in linearly converging fixed-point methods (but not in those converging quadratically). SIAM Journal on Numerical Analysis 2020; 58(1): 788–810.
  * 11 Pollock S, Rebholz LG. Anderson acceleration for contractive and noncontractive operators. IMA Journal of Numerical Analysis 2021; 41(4): 2841–2872.
  * 12 Vuik C. Solution of the discretized incompressible Navier-Stokes equations with the GMRES method. International Journal for Numerical Methods in Fluids 1993; 16(6): 507–523.
  * 13 Vorst V. dHA, Vuik C. GMRESR: a family of nested GMRES methods. Numerical Linear Algebra with Applications 1994; 1(4): 369–386.
  * 14 Chen K, Vuik C. Composite Anderson acceleration method with two window sizes and optimized damping. International Journal for Numerical Methods in Engineering 2022; 123(23): 5964–5985.
  * 15 Toth A, Kelley CT. Convergence analysis for Anderson acceleration. SIAM Journal on Numerical Analysis 2015; 53(2): 805–819.
  * 16 Pollock S, Rebholz LG, Xiao M. Anderson-accelerated convergence of Picard iterations for incompressible Navier–Stokes equations. SIAM Journal on Numerical Analysis 2019; 57(2): 615–637.
  * 17 Sterck HD, He Y. On the asymptotic linear convergence speed of Anderson acceleration, Nesterov acceleration, and nonlinear GMRES. SIAM Journal on Scientific Computing 2021; 43(5): S21–S46.
  * 18 Wang D, He Y, De Sterck H. On the asymptotic linear convergence speed of Anderson acceleration applied to ADMM. Journal of Scientific Computing 2021; 88(2): 1–35.
  * 19 Sterck HD. A nonlinear GMRES optimization algorithm for canonical tensor decomposition. SIAM Journal on Scientific Computing 2012; 34(3): A1351–A1379.
  * 20 Brune PR, Knepley MG, Smith BF, Tu X. Composing scalable nonlinear algebraic solvers. SIAM Review 2015; 57(4): 535–565.
  * 21 Toth A, Ellis JA, Evans T, et al. Local improvement results for Anderson acceleration with inaccurate function evaluations. SIAM Journal on Scientific Computing 2017; 39(5): S47–S65.
  * 22 Peng Y, Deng B, Zhang J, Geng F, Qin W, Liu L. Anderson acceleration for geometry optimization and physics simulation. ACM Transactions on Graphics (TOG) 2018; 37(4): 1–14.
  * 23 Zhang J, Peng Y, Ouyang W, Deng B. Accelerating ADMM for efficient simulation and optimization. ACM Transactions on Graphics (TOG) 2019; 38(6): 1–21.
  * 24 Zhang J, O’Donoghue B, Boyd S. Globally convergent type-I Anderson acceleration for nonsmooth fixed-point iterations. SIAM Journal on Optimization 2020; 30(4): 3170–3197.
  * 25 Bian W, Chen X, Kelley C. Anderson acceleration for a class of nonsmooth fixed-point problems. SIAM Journal on Scientific Computing 2021; 43(5): S1–S20.
  * 26 Chen K, Vuik C. Non-stationary Anderson acceleration with optimized damping. arXiv preprint arXiv:2202.05295 2022\.
  * 27 Khatiwala S. Fast Spin-Up of Geochemical Tracers in Ocean Circulation and Climate Models. Journal of Advances in Modeling Earth Systems 2023; 15(2): e2022MS003447.
  * 28 Pollock S, Schwartz H. Benchmarking results for the Newton–Anderson method. Results in Applied Mathematics 2020; 8: 100095.
  * 29 Herbst MF, Levitt A. Black-box inhomogeneous preconditioning for self-consistent field iterations in density functional theory. Journal of Physics: Condensed Matter 2020; 33(8): 085503.
  * 30 Li Z, Li J. A fast Anderson-Chebyshev acceleration for nonlinear optimization. In: PMLR. ; 2020: 1047–1057.
  * 31 Ji Y, Chen K, Möller M, Vuik C. On an improved PDE-based elliptic parameterization method for isogeometric analysis using preconditioned Anderson acceleration. Computer Aided Geometric Design 2023; 102: 102191.
  * 32 Pratapa PP, Suryanarayana P, Pask JE. Anderson acceleration of the Jacobi iterative method: An efficient alternative to Krylov methods for large, sparse linear systems. Journal of Computational Physics 2016; 306: 43–54.
  * 33 Liu DC, Nocedal J. On the limited memory BFGS method for large scale optimization. Mathematical programming 1989; 45(1-3): 503–528.
  * 34 Haelterman R, Bogaers AE, Scheufele K, Uekermann B, Mehl M. Improving the performance of the partitioned QN-ILS procedure for fluid–structure interaction problems: Filtering. Computers & Structures 2016; 171: 9–17.
  * 35 Pollock S, Rebholz LG. Filtering for Anderson acceleration. arXiv preprint arXiv:2211.12953 2022\.
  * 36 Kelley CT, Suresh R. A new acceleration method for Newton’s method at singular points. SIAM journal on numerical analysis 1983; 20(5): 1001–1009.
  * 37 Moré JJ, Garbow BS, Hillstrom KE. Testing unconstrained optimization software. ACM Transactions on Mathematical Software (TOMS) 1981; 7(1): 17–41.
  * 38 Glowinski R, Keller HB, Reinhart L. Continuation-conjugate gradient methods for the least squares solution of nonlinear boundary value problems. SIAM journal on scientific and statistical computing 1985; 6(4): 793–832.