# Transformers as Policies for Variable Action Environments

Niklas Zwingenberger

###### Abstract

In this project we demonstrate the effectiveness of the transformer encoder as
a viable architecture for policies in variable action environments. Using it,
we train an agent using Proximal Policy Optimisation (PPO) on multiple maps
against scripted opponents in the _Gym- $\mu$RTS_ environment. The final agent
is able to achieve a higher return using half the computational resources of
the next-best RL agent (Huang et al. 2021) which used the _GridNet_
architecture (Han et al. 2019).

##  1 Introduction

In the past years significant progress has been made in variable action
environments, particularly in Real Time Strategy games (RTS). For example the
agent AlphaStar(Vinyals et al. 2019) has achieved a grand-master level in
Starcraft II, consistently defeating human professional players. However, such
state-of-the-art agents require months of training on powerful machines
equipped with hundreds of CPUs and TPUs, which makes them impractical for
research. In response to this, Huang _et al._ have developed a more accessible
environment known as _Gym- $\mu$RTS_(Huang et al. 2021), suitable for training
on a single commercial machine.

Huang _et al._ also demonstrated 2 state-of-the-art methods to train in the
environment _Gym- $\mu$RTS_. These are _Unit Acton Simulation_ (UAS) (Huang et
al. 2021) and _GridNet_ (Han et al. 2019) respectively. However, we found that
both methods suffer from limitations, namely that UAS requires access to the
environment’s model and that _GridNet_ does not computationally scale well
with the size of the environment. Our new approach using a feature map and
transformer architecture mitigates these issues.

In this report, we begin by discussing the relevant background in section 2,
with the description of the architecture in section 3 and relevant results in
section 4. The agent’s implementation, trained model and recordings are made
available in Appendix A.

##  2 Background

###  2.1 Variable Action Environments

We can define a variable action environment as a tuple
$(S,A,P,\rho_{0},\gamma,r,T)$ where $S$ is the state-space, $A$ is a discrete
action space, $P:S\times A\times S\rightarrow[0,1]$ is the transition
probability, $\rho_{0}:S\rightarrow[0,1]$ is the initial state distribution,
$\gamma\in[0,1]$ is the discount factor, $r:S\times A\rightarrow\mathbb{R}$ is
the reward function and $T$ is the maximum episode length. The variable aspect
can be introduced by letting our action space $A=U_{1}\times U_{2}...\times
U_{k}$ be expressed as a Cartesian product of sub-actions $U_{i},i\in{1...k}$
where $k:S\rightarrow\mathbb{N}^{+}$ is the number of sub-actions that needs
to be taken, dependent on the current state.

It should noted that while $k$ can in theory be arbitrarily large, it in most
practical cases is bounded by some maximum number of sub-actions $N$. For
example, in Starcraft II, if we consider each unit $i$ to take some sub-action
from the space $U_{i}$, then $N$ must be equal to the maximum number of units
the game can support (accounts suggest that this ranges from 1700-6400,
depending on machine specs)(sta 2007). For our project, we can show that _Gym-
$\mu$RTS_ also has a clear upper bound $N$, which in practice is never reached
and our architecture takes advantage of this fact.

###  2.2 Transformers

The transformer architecture is a neural network that uses self-attention,
originally introduced for the task of machine translation (Vaswani et al.
2017). Since then it has seen a wide range of applications in NLP tasks
(Devlin et al. 2018), vision tasks (Park and Kim 2022), reinforcement learning
tasks (Bramlage and Cortese 2022) and multi-modal tasks (Lu et al. 2019). In
reinforcement learning in particular, decision transformers (Chen et al. 2021)
have been used for their sequence modelling properties to predict actions
across entire trajectories, which is useful for Non-Markovian environments.

In the context of this project it is sufficient to understand the transformer
layer and multi-head self-attention as presented by (Vaswani et al. 2017):

$\begin{split}\text{Attention}(\mathbf{Q},\mathbf{K},\mathbf{V})=\text{softmax}\left(\frac{\mathbf{Q}\mathbf{K}^{T}}{\sqrt{d_{k}}}\mathbf{V}\right),\\\
\text{MultiHead}(\mathbf{Q},\mathbf{K},\mathbf{V})=\text{Concat}\left(\text{head}_{1},...,\text{head}_{h}\right)\mathbf{W}^{O},\\\
\text{head}_{i}=\text{Attention}(\mathbf{Q}\mathbf{W}_{i}^{Q},\mathbf{K}\mathbf{W}_{i}^{K},\mathbf{V}\mathbf{W}_{i}^{V})\end{split}$

###  2.3 Proximal Policy Optimisation

Policy gradient methods are a popular choice when training agents with large
action spaces as they optimise the policy $\pi_{\theta}$ directly, which
usually allows for a more targeted exploration of promising actions. The
optimisation works by performing a gradient ascent on an objective function
$J(\theta)$, where $\theta$ are the model parameters that we optimise. We
commonly define $J$ to be the expected discounted return of a trajectory:

$\begin{split}J(\theta)=\mathbb{E}_{\tau\sim\pi_{\theta},\rho_{0}}[G_{0}]=\mathbb{E}_{\tau\sim\pi_{\theta},\rho_{0}}\left[\sum_{t=0}^{T-1}\gamma^{t}r_{t}\right],\\\
\tau=(s_{0},a_{0},r_{0},...s_{T-1},a_{T-1},r_{T-1})\end{split}$

Where $\tau$ is the trajectory obtained from sampling the initial distribution
$\rho_{0}$ and parameterised policy $\pi_{\theta}$. We maximise the
expectation by taking the policy gradient $\nabla_{\theta}J(\theta)$ and
updating the parameters using $\theta:=\theta+\alpha\nabla_{\theta}J(\theta)$,
where $\alpha\in(0,1]$ is the learning rate. The gradient itself is given by
(Sutton and Barto 2020):

$\nabla_{\theta}J(\theta)=\mathbb{E}_{\tau\sim\pi_{\theta},\rho_{0}}\left[\sum_{t=0}^{T-1}\nabla_{\theta}\log{\pi_{\theta}(a_{t}|s_{t})G_{t}}\right]$

As we are sampling entire trajectories, this gradient is subject to high
variance(Sutton and Barto 2020) and over the past years, many efforts have
been made to reduce it (Schulman et al. 2015). A particular enhancement of the
above is Proximal Policy Optimisation (PPO)(Schulman et al. 2017), with the
following objective function:

$\begin{split}J(\theta)=\mathbb{E}\left[\min(r_{t}(\theta)\hat{A}_{t},\text{clip}(r_{t}(\theta),1-\epsilon,1+\epsilon)\hat{A}_{t}\right],\\\
r_{t}(\theta)=\frac{\pi_{\theta}(a|s)}{\pi_{\theta_{old}}(a|s)},\\\
\hat{A}_{t}=\sum_{l=0}^{T-1}(\gamma\lambda)^{l}(r_{t}+\gamma
V_{\theta}(s_{t+l+1})-\gamma V_{\theta}(s_{t+l}))\end{split}$

PPO is an actor-critic algorithm, meaning it seeks to optimise both a
parameterised policy $\pi_{\theta}(a|s)$ and a value estimate $V_{\theta}(s)$.
The objective function limits the scope of gradient updates based on the
advantage $\hat{A}_{t}$. In this case, it is the Generalised Advantage
Estimate (GAE) (Schulman et al. 2018).

We have re-used Huang _et al._ ’s implementation of PPO (PPO 2022)(Engstrom et
al. 2020)(Andrychowicz et al. 2020) as we would like to draw comparisons
between our results, focusing on the neural network implementations underlying
$\pi_{\theta}(a|s)$ and $V_{\theta}(s)$.

###  2.4 $\mu$RTS Environment

The environment is a grid-world that supports 1-4 players, where each one
controls units which must gather resources and build structures and more units
to defeat each other in combat. An episode begins with each player’s units
starting in a mirrored section of a map and ends when only one player’s units
remain or some step limit is reached. This environment can be configured to
use different map sizes, layouts, partial observability, stochasticity and
more. In our project we will be using 1 agent to interact with scripted AIs in
the basesWorkers8x8.xml and basesWorkers16x16.xml maps. Both maps are 2 player
and fully observable, differing only in size ($8\times 8$ vs $16\times 16$)
and available resources (50 vs. 100). Please see Figure 1 for an example state
of the $16\times 16$ map.

Figure 1: Example of elements of a $\mu$RTS environment. In the above, there
is a blue player based in the top-left corner of the map and a red player
based in the bottom right. Both players are using workers to harvest resources
from mines and are currently fighting using combat units. Every unit can only
move to free horizontally or vertically adjacent cells and every action
(including moving) takes several time steps to complete.

#### Observation Space

During every step $t$ in an episode we receive an observation tensor $s_{t}$
of shape $(height\times width\times features)$. The first two dimensions
describe a particular cell in the grid-world and the features dimension
describes its state. For example, using the $8\times 8$ map, it has the shape
$(8\times 8\times 27)$, where the 27 features describe 5 pieces of
information: hit points (HP), resources, the owning player, unit type and
current action. Each piece of information is expressed using a one-hot
encoding, where for example we can express a cell having 2 HP as $[0,0,1,0,0]$
or having $\geq 4$ resources as $[0,0,0,0,1]$. For a full description see the
upper section of Table 2 with further examples in Appendix B.

Figure 2: Description observation and action space in $\mu$RTS environment
(Huang et al. 2021). $a_{r}$ is the maximum attack range and is usually set to
7.

#### Action Space

From the above we know that the observation restricts each cell to have at
most 1 unit on it, which limits how many units can be on the map at any time.
Each cell $(x,y)$ with a unit can take some sub-action $U_{x,y}$, allowing at
most $N=height\times width$ sub-actions to be taken per time step. The
components and possible values of an action are shown in the lower section of
Table 2. Using them we can show that for some cell $(x,y)$ there are $6\times
4\times 4\times 4\times 4\times 7\times 7^{2}=526848$ possible discrete values
that a sub-action can assume. This is quite large number and having a policy
network generate half a million logits for every grid cell is too expensive,
especially when taking the gradient afterwards.

Instead we make an independence assumption between each action component
(Huang et al. 2021) and define the probability of our policy
$\pi_{\theta,x,y}(u_{t}|s_{t})$ choosing a sub-action $u_{t}\in U_{x,y}$ given
some $s_{t}\in S$ as:

$\begin{split}\pi_{\theta}(a_{t}|s_{t})=\prod_{x\in height}\prod_{y\in
width}\pi_{\theta,x,y}(u_{t}|s_{t}),\\\
\pi_{\theta,x,y}(u_{t}|s_{t})=\prod_{u^{d}_{t}\in
D}\pi_{\theta}(u^{d}_{t}|s_{t}),\\\ D=\\{u_{t}^{\text{Source
Unit}},u_{t}^{\text{Action Type}},u_{t}^{\text{Move
Parameter}},u_{t}^{\text{Harvest Parameter}},\\\ u_{t}^{\text{Return
Parameter}},u_{t}^{\text{Produce Direction Parameter}},u_{t}^{\text{Produce
Type Parameter}},\\\ u_{t}^{\text{Relative Attack Position}}\\}\end{split}$

where $\pi_{\theta}$ is the joint policy of each individual policy
$\pi_{\theta,x,y}$ per grid cell and $D$ is the set of action components. This
way only $6+4+4+4+4+7+7^{2}=78$ logits are needed per sub-action.

Additionally, the _Gym- $\mu$RTS_ also provides built-in invalid action
masking (Huang and Ontañón 2022), which replaces the logits for nonsensical
actions with high negative values. Then when we take the softmax over a
distribution parameterised by these logits, the probability of such actions is
0. As an example, light, heavy and ranged combat units cannot produce any
buildings or units, so the masking sets the probability of the "Action Type
Parameter" action component assuming value "produce" to 0.

###  2.5 GridNet

The _GridNet_ approach is a policy that $\pi_{\theta,x,y}(a_{t}|s_{t})$ yields
a sub-action $u_{i}$ for every grid cell $i\in\\{1,2,...,height\times
width\\}$. This approach is convenient as it is quite easy to implement a
policy that takes a fixed-size input (e.g the grid world) and yields a fixed-
size output (action logits for every cell). The _GridNet_ itself uses a
convolutional neural net (CNN) for this. A problem is that in practice for
_Gym- $\mu$RTS_ and most other RTS games, controllable units do not exist on
every cell, so a large chunk of the sub-actions are computed only to be thrown
away. This computational burden also hinders the approach from being easily
scalable to larger grid sizes.

###  2.6 Unit Action Simulation

UAS works by iterating over unit $i\in\\{1,2,...N\\}$ controlled by the agent,
obtain a sub-action $u_{i}$, simulate how $u_{i}$ will play out in the
environment and get a new simulated state $s_{t}^{(i)}$ that is then used to
get the action for unit $i+1$. The key here is that only a policy
$\pi_{\theta,x,y}(u_{t}|s_{t}^{(i)})$ for sub-action space is needed rather
than the action space. This results in no wasted logits being generated for
empty grid cells as opposed to _GridNet_. The problem is that to acquire the
simulated state we need access to the environment’s transition probability
function $P$ to give us this next state. This effectively makes this approach
model-based and it also only reasonable works for deterministic environments,
as in a stochastic one, there would be several possible simulated states from
each sub-action.

##  3 Method

Given an observation $s_{t}$ of shape $(height\times width\times 27)$ we have
developed a scalable policy network that can convert it into $k_{t}$ logits of
shape $(k_{t}\times 78)$ where $k_{t}$ is the number units the agent controls.
Additionally, we yield a corresponding value estimate
$\hat{V}_{\theta}(s_{t})$. A visual overview of the architecture is provided
in Appendix C and with an explanation of each component being in the following
sections. Further, all used hyper-parameters are referenced in Appendix E.

###  3.1 Feature Map

The first step is to use a feature map $\phi(s_{t})$ that converts $s_{t}$
into $e_{t}=k_{t}+l_{t}+m_{t}$ individual representations, where $e_{t}$ is
the number of entities, $k_{t}$ is the number of agent-controlled units,
$l_{t}$ is the number of enemy-controlled units and $m_{t}$ is the number of
neutral units (e.g resource mines). Our implementation of $\phi(s_{t})$ is
shown in Algorithm 1 and expresses each representation as a one-hot vector of
the position (length $hw$) concatenated with the feature vector of the entity
(length $27$).

Algorithm 1 Feature Map $\phi(s_{t})$

Input: $s_{t}$ with shape $(h\times w\times 27)$
Output: $x_{t}$ with shape $(e_{t}\times(hw+27))$

1: Let $v_{t}=\text{reshape}(s_{t})$ with shape $(hw\times 27)$

2: Let $x_{t}\leftarrow\mathbf{0}$

3: for each entity index $i\in\\{1,2,...,e_{t}\\}$ do

4: $p=\text{position}(v_{t},i)$ {gets index position on grid}

5: $x_{t}[i,\mathbin{:}hw]=\text{onehot}(p)$

6: $x_{t}[i,hw\mathbin{:}]=v_{t}[p,:]$

7: end for

8: return $x_{t}$

For example, this means that for the $8\times 8$ map, $\phi(s_{t})$ yields an
$(e_{t}\times(64+27))=(e_{t}\times 91)$ matrix. Similarly, for the $16\times
16$ map, we obtain a $(e_{t}\times(256+27)=(e_{t}\times 283)$ matrix. We can
observe that for the larger map, the size of our representation grows
significantly. This does not scale well, but fortunately a one-hot encoding
can be represented much more densely by using a trainable embedding:

$\begin{split}\boldsymbol{p}_{embed}=\text{onehot}(p)\mathbf{W}_{embed}\\\
\boldsymbol{p}_{embed}\in\mathbb{R}^{d},\quad\mathbf{W}_{embed}\in\mathbb{R}^{(hw\times
d)}\end{split}$

This can mitigate the size increase caused by larger maps and in our
implementation, we set our embedding weights to be
$\mathbf{W}_{embed}\in\mathbb{R}^{(256\times 64)}$. As a result, the
$\phi(s_{t})$ generates the same size output of $(e_{t}\times 91)$ for both
the $8\times 8$ and $16\times 16$ map. We should note here, that this size
reduction is obviously not infinitely possible and that in particular the
representation size should be as least as large as the number of sub-action
logits (e.g $91>78$). Otherwise it becomes more difficult to obtain meaningful
logits.

###  3.2 Transformer Net

With an $x_{t}=\phi(s_{t})$ we can now input this into a transformer encoder.
The encoder is composed of 5 stacked, identical transformer layers each of
which has 7 attention heads and uses 512 neurons in the dense layer with a
ReLU activation (Agarap 2018) (other hyper-parameters detailed in Appendix).
In our implementation the entity representations in $x_{t}$ are ordered as
follows by player: player 1 (agent), player 2 (opponent) and neutral units. As
we do not use a positional embedding, this does not affect the encoder, but it
makes the encoder’s output $y_{t}\in\mathbb{R}^{e_{t}\times 91}$ easier to
process by the actor and critic.

###  3.3 Actor & Critic

Both the actor and critic networks share this encoder and use $y_{t}$. The
actor takes in all $k_{t}$ outputs related to the agent’s units and then
applies a shared single feed-forward layer (with $91\times 78+78=7176$
parameters) to each of them. This way the actor is invariant to the number of
inputs it receives and manages to yield variable sub-action logits
$z_{1},...z_{k_{t}}$.

The critic is a little more complex. As before, we can apply another single
feed-forward layer (with $91\times 1+1=92$ parameters), but this leaves us
with $e_{t}$ different value estimates
$\hat{V}_{\theta,1}(s_{t})...\hat{V}_{\theta,e_{t}}(s_{t})$ that must somehow
be combined into one value $\hat{V}_{\theta}(s_{t})$. A sum naturally comes to
mind:

$\hat{V}_{\theta}(s_{t})=\sum_{i=1}^{e_{t}}\hat{V}_{\theta,i}(s_{t})$

However, this creates some variance issues. We know that
$\hat{V}_{\theta,1}(s_{t})...\hat{V}_{\theta,e_{t}}$ are the values of player
units, enemy units and neutral resources. Intuitively, we would expect that a
large number of player units means the agent is doing well ($k_{t}>l_{t}$),
yielding a higher value estimate, whereas a large number of enemy units would
yield a lower estimate ($k_{t}<l_{t}$). Higher and lower is difficult to model
without knowing which value estimate belongs to which player. We could rely on
the transformer to provide this information for our critic, but this is an
unnecessary burden. Instead, we should try to obtain some more meaningful
metrics that the critic can learn to weight on its own:

$\begin{split}\hat{V}_{\theta}(s_{t})=\sum_{i\in\\{p_{1},p_{2},p_{3}\\}}w_{i,1}\Sigma_{i}+b_{i,1}+w_{i,2}\mu_{i}+b_{i,2}\\\
\Sigma_{i}=\sum_{\text{unit }j\text{ of player
}i}\hat{V}_{\theta,j}(s_{t}),\\\ \mu_{i}=\frac{\Sigma_{i}}{\text{unit count of
player }i}\end{split}$

Here we calculate 6 metrics: the sum $\Sigma_{i}$ and mean $\mu_{i}$ of each
player’s ($\in\\{p_{1},p_{2},p_{3}\\}$) value estimates (here $p_{3}$ refers
to the neutral units). We then apply another feed-forward layer (with weights
$w_{{p_{1}},1}...w_{{p_{3}},2}$ and biases $b_{{p_{1}},1}...b_{{p_{3}},2}$) so
the critic can choose how to weight each. This allows for the critic to make
the value estimates invariant of unit counts and weight units differently
based on their owner.

###  3.4 Training

In order to train this model the agent conducts exploratory episodes of up to
2000 steps in 24 parallel instances. The opposing player in these instances is
one of 4 scripted AIs:

  * •

randomBiasedAI \- chooses random moves for each unit, but is 5 times likelier
to harvest resources or attack.

  * •

workerRushAI \- implements a simple rush strategy where it constantly trains
workers to attack the closest enemy unit.

  * •

lightRushAI \- implements a simple rush strategy where it builds a barracks
and constantly sends light military units to attack.

  * •

coacAI \- a mixed AI that won the 2020 _$\mu$ RTS_ contest.

These AIs and their rankings are publicly available(Santiontanon 2018)(mic
2020) and there are 9 more available in the _Gym- $\mu$RTS_ environment, which
will be used for evaluation in Section 4.

As an episode can be quite long, we opted to use shaped rewards in order to
ensure a more stable and reproducible training (Schrittwieser et al. 2019).
The agent receives a reward whenever it completes any of the following:

  * •

$-10$/$0$/$+10$ for losing/drawing/winning, where drawing occurs when the
episode step limit is reached with both players still having units.

  * •

$+1$ for harvesting a resource from a mine.

  * •

$+1$ for attacking an enemy unit.

  * •

$+0.2$ for constructing a building.

  * •

$+1$ for creating a worker.

  * •

$+4$ for producing a combat unit.

After collecting ($24\times 256$) steps, we then run the PPO gradient ascent
and loss calculation for 4 epochs. Each particular update is conducted in
mini-batches of 4. Generally, the opponent selection, reward structure and
major training hyper-parameters are kept as similar as possible to (Huang et
al. 2021), so a more meaningful comparison can be made in the results.

##  4 Results

We trained 2 models for evaluation: one for the $8\times 8$ and one for
$16\times 16$ map. Both models used the same training hyper-parameters with
the exception of the $16\times 16$ model additionally using the embedding.
Accordingly, their trainable parameter counts are very similar: 645470 and
661854 respectively. They were each trained on 2 CPUs and 1 GPU (Tesla
V100-SXM2-16GB) with 16GB of RAM for 100 million steps.

The models were each evaluated for 100 games against each of the 13 scripted
AIs available in the Gym-$\mu$RTS environment, with a summary of their
performance in Table 1 and a more detailed breakdown in Appendix D.
Additionally, we benchmark against Huang _et al._ ’s results, comparing shaped
rewards and also training duration given that similar computational resources
were used.

Analysing the results on the $8\times 8$ map, the Transformer Net has achieved
a high win-rate of 91% within a reasonable training time, although not much
lower than any of the models used for the $16\times 16$ map. Its shaped reward
is significantly lower, but this is to be expected as the $8\times 8$ only has
half of the harvestable resources of the $16\times 16$ map. The proximity to
the opponent’s base also does not encourage the slow process of building
military buildings followed by training high reward military units. Instead,
the agent adopts a more refined "worker rush" strategy where it attacks early
and quickly only using workers.

More curious are the results of Transformer Net w. Embedding on the $16\times
16$ map. It appears to have achieved the highest shaped return of 189.7, but
the lowest win-rate of 76%. The likely cause here is an over-fitting of the
shaped rewards. A win/loss reward difference of 20 is simply not significant
enough (see Figure 4). An inspection of the evaluation recordings of the agent
supports this theory. The agent appears to have learned to cripple the
opposing player by defeating all units save for a single survivor. Then it
proceeds to artificially prolong the episode by continuing to harvest its own
resources and even starts harvesting the opponent’s resources. Eventually it
finishes off the remaining unit and ends the episode. Unfortunately, this
strategy does not always succeed and does not correlate strongly with
consistently winning episodes. However, we would argue that this is ultimately
an issue of reward design, not the model. Achieving the highest return in a
comparatively short training time shows a lot of promise.

Figure 3: Shaped return obtained during training of the transformer net on the
$8\times 8$ map. The chart on the left shows the total shaped return, whereas
the one on the right shows a breakdown by reward category. We can observe here
that the agent derived the most reward from harvesting resources and
attacking. Building barracks and training combat units were not rewarding for
it all.

Figure 4: Shaped return obtained during training of the transformer net on the $16\times 16$ map. The chart on the left shows the total shaped return, whereas the one on the right shows a breakdown by reward category. Here we can see that the agent obtained almost half of all rewards from harvesting resources, followed by producing military units. The win/loss reward was much less significant by comparison. Model | Map Size | Win-rate | Shaped Return | Training Time (hours)
---|---|---|---|---
GridNet Encoder-Decoder | $16\times 16$ | 0.89 | 180.8 | 117.03
UAS Impala CNN | $16\times 16$ | 0.91 | 137.6 | 63.67
Transformer Net | $8\times 8$ | 0.91 | 84.9 | 55.30
Transformer Net w. Embedding | $16\times 16$ | 0.76 | 189.7 | 58.87
Table 1: Evaluation win-rates, shaped return and training times for currently
known best performing RL agents in _Gym- $\mu$RTS_. The results for the
GridNet Encoder-Decoder and UAS Impala CNN are sourced from (Huang et al.
2021). It should be noted that most $8\times 8$ results are not really
comparable to the $16\times 16$ as the environments are too different. They
are kept in the same table solely for convenience.

##  5 Evaluation & Future Work

Overall, we believe this transformer-based policy to be a viable approach for
variable action environments such as _Gym- $\mu$RTS_. It has managed to
demonstrate a better reward maximisation than UAS and faster training time
than GridNet.

There are of course, a number of open problems. It is not yet known how this
transformer policy will work with sparse rewards nor how well it will
generalise from self-play. Further, it has yet to be determined how far this
can be scaled. Recall that a forward pass through a self-attention layer has a
computational complexity of $O(n^{2}d)$ (Vaswani et al. 2017), where $n$ is
the number of inputs and $d$, their length. In _Gym- $\mu$RTS_ we recorded a
maximum of $n=40$ entities being present on a $16\times 16$ map at any one
time (see Appendix F). Making $40^{2}$ comparisons is quite reasonable, but is
this feasible for $n=1000$ or more? This remains to be seen.

##  6 Acknowledgements

I would like to thank Shengyi Huang for his excellent paper(Huang et al. 2021)
and well-documented repository (Vwxyzjn 2022) that I could use as a starting
point for this work. Further, I would also like to thank him for taking the
time to answer my questions in a friendly and patient manner.

## References

  * sta (2007) 2007\.  Map Max.  https://starcraft.fandom.com/wiki/Map_max#:~:text=References%2C%20and%20hitting%20this%20limit.%22.
  * mic (2020) 2020\.  2020 Micro-RTS AI competition results.  https://sites.google.com/site/micrortsaicompetition/competition-results/2020-cog-results?authuser=0.
  * PPO (2022) 2022\.  The 37 implementation details of proximal policy optimization.  https://iclr-blog-track.github.io/2022/03/25/ppo-implementation-details/.
  * Agarap (2018) Agarap, A. F. 2018.  Deep learning using rectified linear units (relu).  _arXiv preprint arXiv:1803.08375_.
  * Andrychowicz et al. (2020) Andrychowicz, M.; Raichuk, A.; Stanczyk, P.; Orsini, M.; Girgin, S.; Marinier, R.; Hussenot, L.; Geist, M.; Pietquin, O.; Michalski, M.; Gelly, S.; and Bachem, O. 2020.  What Matters In On-Policy Reinforcement Learning? A Large-Scale Empirical Study.  _CoRR_ , abs/2006.05990.
  * Bramlage and Cortese (2022) Bramlage, L.; and Cortese, A. 2022.  Generalized attention-weighted reinforcement learning.  _Neural Networks_ , 145: 10–21.
  * Chen et al. (2021) Chen, L.; Lu, K.; Rajeswaran, A.; Lee, K.; Grover, A.; Laskin, M.; Abbeel, P.; Srinivas, A.; and Mordatch, I. 2021.  Decision Transformer: Reinforcement Learning via Sequence Modeling.  _CoRR_ , abs/2106.01345.
  * Devlin et al. (2018) Devlin, J.; Chang, M.; Lee, K.; and Toutanova, K. 2018.  BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding.  _CoRR_ , abs/1810.04805.
  * Engstrom et al. (2020) Engstrom, L.; Ilyas, A.; Santurkar, S.; Tsipras, D.; Janoos, F.; Rudolph, L.; and Madry, A. 2020.  Implementation Matters in Deep Policy Gradients: A Case Study on PPO and TRPO.  _CoRR_ , abs/2005.12729.
  * Han et al. (2019) Han, L.; Sun, P.; Du, Y.; Xiong, J.; Wang, Q.; Sun, X.; Liu, H.; and Zhang, T. 2019\.  Grid-Wise Control for Multi-Agent Reinforcement Learning in Video Game AI.  In _ICML_.
  * Huang et al. (2021) Huang, S.; Ontañón, S.; Bamford, C.; and Grela, L. 2021.  Gym-$\mu$RTS: Toward Affordable Full Game Real-time Strategy Games Research with Deep Reinforcement Learning.  _CoRR_ , abs/2105.13807.
  * Huang and Ontañón (2022) Huang, S.; and Ontañón, S. 2022.  A Closer Look at Invalid Action Masking in Policy Gradient Algorithms.  volume 35.
  * Lu et al. (2019) Lu, J.; Batra, D.; Parikh, D.; and Lee, S. 2019.  ViLBERT: Pretraining Task-Agnostic Visiolinguistic Representations for Vision-and-Language Tasks.  _CoRR_ , abs/1908.02265.
  * Park and Kim (2022) Park, N.; and Kim, S. 2022.  How Do Vision Transformers Work?
  * Santiontanon (2018) Santiontanon. 2018.  Artificial Intelligence - Microrts Wiki.  https://github.com/santiontanon/microrts/wiki/Artificial-Intelligence.
  * Schrittwieser et al. (2019) Schrittwieser, J.; Antonoglou, I.; Hubert, T.; Simonyan, K.; Sifre, L.; Schmitt, S.; Guez, A.; Lockhart, E.; Hassabis, D.; Graepel, T.; Lillicrap, T. P.; and Silver, D. 2019.  Mastering Atari, Go, Chess and Shogi by Planning with a Learned Model.  _CoRR_ , abs/1911.08265.
  * Schulman et al. (2015) Schulman, J.; Levine, S.; Moritz, P.; Jordan, M. I.; and Abbeel, P. 2015.  Trust Region Policy Optimization.  _CoRR_ , abs/1502.05477.
  * Schulman et al. (2018) Schulman, J.; Moritz, P.; Levine, S.; Jordan, M.; and Abbeel, P. 2018.  High-Dimensional Continuous Control Using Generalized Advantage Estimation.  _arXiv:1506.02438 [cs]_.  ArXiv: 1506.02438.
  * Schulman et al. (2017) Schulman, J.; Wolski, F.; Dhariwal, P.; Radford, A.; and Klimov, O. 2017.  Proximal Policy Optimization Algorithms.
  * Sutton and Barto (2020) Sutton, R. S.; and Barto, A. G. 2020.  _Reinforcement learning: An introduction_.  The MIT Press.
  * Vaswani et al. (2017) Vaswani, A.; Shazeer, N.; Parmar, N.; Uszkoreit, J.; Jones, L.; Gomez, A. N.; Kaiser, L.; and Polosukhin, I. 2017.  Attention Is All You Need.  _CoRR_ , abs/1706.03762.
  * Vinyals et al. (2019) Vinyals, O.; Babuschkin, I.; Chung, J.; Mathieu, M.; Jaderberg, M.; Czarnecki, W.; Dudzik, A.; Huang, A.; Georgiev, P.; Powell, R.; Ewalds, T.; Horgan, D.; Kroiss, M.; Danihelka, I.; Agapiou, J.; Oh, J.; Dalibard, V.; Choi, D.; Sifre, L.; Sulsky, Y.; Vezhnevets, S.; Molloy, J.; Cai, T.; Budden, D.; Paine, T.; Gulcehre, C.; Wang, Z.; Pfaff, T.; Pohlen, T.; Yogatama, D.; Cohen, J.; McKinney, K.; Smith, O.; Schaul, T.; Lillicrap, T.; Apps, C.; Kavukcuoglu, K.; Hassabis, D.; and Silver, D. 2019.  AlphaStar: Mastering the Real-Time Strategy Game StarCraft II.  https://deepmind.com/blog/alphastar-mastering-real-time-strategy-game-starcraft-ii/.
  * Vwxyzjn (2022) Vwxyzjn. 2022.  Vwxyzjn/gym-microrts-paper: The source code for the gym-microrts paper.  https://github.com/vwxyzjn/gym-microrts-paper.

##  Appendix A Source Code, Trained Models and Recordings

We provide additional resources under the following URLs:

  * •

source code and trained models: https://github.com/NiklasZ/transformers-for-
variable-action-envs

  * •

evaluation game recordings:
https://drive.google.com/drive/folders/1I63L0WyEyN0v1rvd7H5SN-27VbdRkvOq?usp=sharing

  * •

training run data and charts:
https://wandb.ai/niklasz/public_var_action_transformers

##  Appendix B Observation Examples

Consider the example board in Figure 5. We can express the top-left resource
cell using the following encoding:

$[1,0,0,0,0],[0,0,0,0,1],[0,1,0],[0,1,0,0,0,0,0,0],[1,0,0,0,0,0]$

Where the first vector tells us it has $0$ HP, the second tells us it has
$\geq 4$ resources, the 3rd tells us it is owned by the neutral player, the
4th tells us it is a resource and the last one tells us it is presently doing
nothing. Similarly, we can express the worker to the right of it using:

$[0,1,0,0,0],[0,1,0,0,0],[1,0,0],[0,0,0,0,1,0,0,0],[0,1,0,0,0,0]$

In this example we see the worker has $1$ HP and is carrying $1$ resource. It
belongs to player 1 and is a worker unit. Its current action (indicated by the
line pointing out of it) is to move.

Figure 5: The above depicts 2 tiles of resources in green, and 1 base and 2
workers for each player. Each base is currently constructing a worker. The
blue workers have a "1" on their sprite as they are each carrying a resource.

##  Appendix C Agent Neural Network Architecture

Figure 6: Neural Network Architecture of the Actor and Critic. We first apply
the feature map $\phi(s_{t})$ and then order the resulting vectors by unit
ownership: player 1, player 2, player 3 (neutral entities like resources).
Next we feed it through the transformer encoder layers. The actor then takes
the outputs from player 1 and creates action component logits
$z_{1}...z_{k_{t}}$, whereas the critic considers the outputs from all
players. From them, it calculates respective mean and sum values which it then
weights into a final value prediction $\hat{V}_{\theta}(s_{t})$. Lastly, note
that in the actual implementation we also include a batch dimension $B$,
meaning the real input has shape $B\times h\times w\times 27$ and the actor
output will be $B$ lists of logits $z_{1}...z_{k^{b}_{t}}$ where $k^{b}_{t}$
will differ for each list $b\in B$.

##  Appendix D Agent Win-Rates

Figure 7: Win-rate of the $8\times 8$ agent during training, smoothed via
exponentially moving average. It took the longest to train against the
workerRushAI and coacAI. Figure 8: Match stats of the $8\times 8$ agent during
evaluation. Ties occur when neither AI manages to completely eliminate the
opponent within the episode step limit $T=2000$. Figure 9: Win-rate of the
$16\times 16$ agent during training, smoothed via exponentially moving
average. We can observe that the agent managed to win most of the time against
its training opponents, but never achieved a decisive win-rate. Towards the
end, the win-rate even degraded, likely because the agent was not sufficiently
rewarded/punished for winning/losing. Figure 10: Match stats of the $16\times
16$ agent during evaluation. Ties occur when neither AI manages to completely
eliminate the opponent within the episode step limit $T=2000$.

##  Appendix E Agent Hyper-parameters

Hyper-parameter Name |  Setting |  Notes
---|---|---
Transformer Layers |  5 |  Chosen to be the best efficacy/efficiency trade-off from ranges of 3 to 8.
Transformer feed-forward neurons |  512 |  -
Transformer attention heads |  7 |  Chosen mostly as it fits the input size $91|7$
Transformer activation function |  ReLU |  -
Transformer Dropout |  0.1 |  -
Weight Initialisation |  $\mathcal{N}(0,2)$ |  -
Bias Initialisation |  0 |  -
Optimiser |  Adam with: $\alpha=2.5\mathrm{e}{-4}$ $\epsilon=1\mathrm{e}{-5}$ $\beta_{1}=0.9$ $\beta_{2}=0.999$ with a linear decay to 0. |  -
Max Training Steps |  100 Million |  This refers to the number of episode steps that have been used for gradient calculation.
Number of exploration steps |  256 |  How many steps to traverse per environment before running PPO.
Parallel bot environments |  24 |  Number of simultaneous environments for exploration.
Minibatch size |  4 |  -
Return discount factor |  $\gamma=0.99$ |  -
Generalised Advantage Estimate |  $\lambda=0.95$ |  -
Entropy coefficient |  $c=0.01$ |  Controls weighting of entropy loss.
Value Function coefficient |  $v_{coef}=0.5$ |  Controls weighting of value function loss.
Gradient clipping coefficient |  $\epsilon=0.1$ |  PPO clipping coefficient
Table 2: Training hyper-parameters

##  Appendix F Entity Distribution Counts

Figure 11: Distribution of entities counted in every state during evaluation
of the $8\times 8$ transformer net. On average around 12 entities are present,
corresponding to roughly $19\%$ of the cells in the grid. Figure 12:
Distribution of entities counted in every state during evaluation of the
$16\times 16$ transformer net. On average around 20 entities are present,
corresponding to roughly $8\%$ of the cells in the grid.


∎

11institutetext:  1The Rotman Institute of Philosophy and
The Ontario Research Center for Computer Algebra and
The School of Mathematical and Statistical Sciences
The University of Western Ontario
2School of Information Technology and Mathematical Sciences
University of South Australia

# Optimal Residuals and the Dahlquist Test Problem

Robert M. Corless1 C. Yalçın Kaya2 and Robert H. C. Moir1

(Received: date / Accepted: date)

###### Abstract

We show how to compute the optimal relative backward error for the numerical
solution of the Dahlquist test problem by one-step methods. This is an example
of a general approach that uses results from optimal control theory to compute
optimal residuals, but elementary methods can also be used here because the
problem is so simple. This analysis produces new insight into the numerical
solution of stiff problems.

###### Keywords:

stiff IVPbackward error residual optimal control

###### MSC:

65L0465L05 65L20 49M05

††journal: Numerical Algorithms

##  1 Introduction

The study of stiff differential equations and their efficient numerical
solution is by now a mature field. There are several, perhaps many, efficient
practical methods with freely available high quality implementations. The
literature on the theory of such methods is extensive. Surprisingly, it is not
yet complete: for instance, see the survey soderlind2014stiffness , which has
the intriguing title Stiffness 1952-2012: Sixty years in search of a
definition. That paper re-examines the fundamentals and thoroughly surveys the
literature, and proposes a new stiffness indicator that they claim is useful
both a priori for indicating stiffness and a posteriori for indicating varying
regions of the solution that stiffness is important.

This paper takes a different approach, that of optimal residuals, i.e.,
backward error, and uses it on the Dahlquist test problem to generate some new
observations about this, the simplest of all stiff problems. Indeed,
soderlind2014stiffness calls this problem “simplistic” and with good reason,
but surprisingly it still has things to teach us.

Trying to study stiff problems from the point of view of backward error
analysis is itself not new. For instance, there is the PhD thesis of W.L.
Seward and the paper enright1989achieving . But there is an intrinsic
dissonance: a stiff problem has the feature that errors are (often sharply)
damped as the integration proceeds forward in time, and thus it is not obvious
why one might prefer to look at backward error when, if the right method is
used, small forward error happens more or less automatically.

We claim that backward error, in particular optimal residuals, which we
explain below, really can be useful. One way to see the usefulness of this
type of analysis is in contrast to the classical stability analysis on the
Dahlquist test problem, $\dot{y}=\lambda y$, $\lambda\in\mathbb{C}$. In the
classical stability analysis there is an emphasis on A-stability, with an
understanding that, at least locally where linear stability assessments are
valid, a basic criterion for stability of a method is that the numerical
solution decays where the actual solution does. This condition is of
particular concern for stiff problems, since it entails that the method will
not be subject to stability restrictions on account of eigenvalues with large
negative real part.

The classical stability analysis really does explain most of the behaviour of
implicit methods on stiff problems, and much insight has been gained thereby.
This present analysis is a refinement only, that offers the possibility of
explaining some extra considerations which are “well-known” to practitioners,
namely that, under certain circumstances, a stiff method may well be stable
but not as accurate as one might wish. With the optimal residual approach a
rather different set of criteria emerge, which imply a de-emphasis of the
criterion of A-stability, being enhanced with a consideration of those regions
of the problem space that a given method can reproduce the dynamics of the
reference problem accurately; such regions often extend comfortably into the
right-half of the complex plane, even for explicit methods.

There are also deeper reasons why backward error analysis is useful on stiff
problems, as a result of the fact that a small backward error entails a small
perturbation of the dynamics of the problem. In this regard, one place where
optimal residuals are especially useful is in the solution of systems which
have nontrivial attracting sets: the problems can be stiff because decay to
the attracting set can be very strong, but small backward error is important
too in order to get the dynamics right on the attracting set, which might in
fact be actually chaotic. Of course, for a chaotic problem, small forward
error is not possible at all (at least, not for very long), but as explained
for instance in corless1994good small backward error can be perfectly
satisfactory as an explanation of the success of a numerical method on a
chaotic problem.

The question of what good a numerical method is, even one providing a small
backward error, for computation of an attracting set for a chaotic system,
deserves a repetition of the answer given in that reference. Every numerical
analyst knows that small forward error requires exponentially accurate initial
conditions and exponentially accurate integration. One way to explain a
successful computation is if some form of shadowing is invoked. That is, the
forward accuracy of the computed trajectory is explained by being“shadowed” by
an exact solution of the reference problem, typically from some nearby initial
condition. This is a form of backward error analysis moir2010reconsidering .
Computationally verifying that shadowing has occurred is expensive, however,
and while shadowing is generic, there are no a priori guarantees that
shadowing will occur or has occurred.

A cheaper and simpler method to explain the success of a numerical method on a
chaotic problem is to verify that the residual (also called the defect,
deviation, or slope error) is small. Also take note of the utility of
interpreting ordinary numerical errors as a modelling error in many
circumstances. If the model equations are written

$\dot{y}=f(t,y)$

after having been derived for some physical context, universally by neglecting
minor effects, and our numerical solution gives us $z(t)$ exactly satisfying

$\dot{z}(t)=f(t,z(t))+\Delta(t).$

or perhaps

$\dot{z}(t)=(I+\delta(t))f(t,z(t)),$

and $\Delta$ or $\delta$ are small compared to the neglected terms, then as
Nick Higham puts it “for all we know, we already have the exact solution,” for
error has been introduced into the dynamics by modelling assumptions.

We also emphasize that there must be some feature of this system that is
robust under perturbations, or else even the reference solution of
$\dot{y}=f(t,y)$ would be useless in the face of real-life perturbations. The
existence of such a feature for a given system was termed “well-enough
conditioning” in corless1994error . Nearly all models used have this property,
even chaotic ones. For instance, the dimension of the Lorenz attractor is
quite robust under forcings of this type: the attracting sets of

$\displaystyle\dot{y_{1}}$ $\displaystyle=$
$\displaystyle\sigma(y_{2}-y_{1})+\delta_{1}(y),$ $\displaystyle\dot{y_{2}}$
$\displaystyle=$ $\displaystyle y_{1}(\rho-y_{3})-y_{2}+\delta_{2}(y),$
$\displaystyle\dot{y_{3}}$ $\displaystyle=$ $\displaystyle y_{1}y_{2}-\beta
y_{3}+\delta_{3}(y),$

are remarkably close to one another even for quite large
$\Delta(t)=(\delta_{1}(t),\delta_{2}(t),\delta_{3}(t))^{\text{T}}$.
corless1994error

So backward error may be important in explaining the success or failure of
numerical methods for chaotic systems, which can be stiff. More than this,
however, we believe our approach to stability analysis on the Dahlquist test
problem may lead to a refinement of the classical explanation of the practical
success of various numerical methods.

We now consider backward error and the general context ofoptimal residuals
briefly, before moving on to an elementary analysis of the Dahlquist test
problem.

##  2 Optimal Residuals and Backward Error

###  2.1 Backward Error, Residuals and Interpolants

Most codes supply interpolants: for graphical output, for event location, or
for handling delays. Given an interpolant, which we will call $z(t)$, to the
numerical skeleton $(t_{k},y_{k})$ of our computed solution to the ODE
$\dot{y}=f(t,y)$, we define the residual $\Delta(t)$ as

$\Delta(t):=\dot{z}(t)-f(t,z(t))\>.$ (1)

As previously stated, this is sometimes called the defect. In one sense, this
is a kind of backward error, because the computed $z(t)$ can be interpreted as
the exact solution of the perturbed equation

$\dot{y}(t)=f(t,y(t))+\Delta(t)\>.$ (2)

Note that the residual as defined here is dependent on just which interpolant
$z(t)$ is used.

To be compatible with the accuracy of the numerical solutions they
interpolate, these interpolants should be $O(h^{p})$ accurate, but sometimes
are not. For example, in Matlab, ode45 uses a fifth-order Runge-Kutta Fehlberg
formula, but has only a fourth-order interpolant: so $z^{\prime}(t)$ will only
be third-order accurate. A consequence is that the residual, or defect,
$\Delta(t):=\dot{z}(t)-f(t,z)$ will sometimes be overestimated.

The work of enright1989new ; enright2000continuous , enright1991parallel ,
enright2007robust , muir2003pmirkdc , cao2004posteriori has shown that one can
interpolate the skeleton $\\{t_{k},y_{k}\\}$ of a numerical solution in a
practical way that gives the correct asymptotic size of the residual
$\|\Delta\|=O(h^{p})$ as the mean time step $h\rightarrow 0$, for a method of
order $p$. This gives a practical method that gives tolerance proportionality
and robust reliability for tight tolerances.

If instead one interpolates the skeleton badly, one will obviously get a large
residual. In one of their examples, Hubbard & West wonder if a smaller
residual (they term it “slope error”) might be achieved with some other
interpolant HubbardWest:1991 . This raises the question of how to find the
“best” interpolant, which gives the “smallest” residual. This smallest
residual is the one that will most accurately measure how good a job the
underlying method did in generating the skeleton.

This question is answered in general elsewhere (Corless, Kaya & Moir, in
preparation) where optimization methods are used to find interpolants from
$(t_{k},y_{k})$ to $(t_{k+1},y_{k+1})$ that minimize $\|\Delta(t)\|$ or
$\|\delta(t)\|$ (either the $2$-norm or $\infty$-norm are handled).

###  2.2 Interpolants via optimal control theory

For ease of exposition in this section we work entirely in $\mathbb{R}^{n}$.

Define the control variable vector $u(t):=\delta(t)$, where
$u:[t_{i},t_{i+1}]\to\mathbb{R}^{n}$, as a piecewise continuous vector
function. The problem of finding interpolants
$z:[t_{i},t_{i+1}]\to\mathbb{R}^{n}$ from $(t_{i},y_{i})$ to
$(t_{i+1},y_{i+1})$ that minimize the $L^{\infty}$-norm of the relative error
$\|\delta\|_{L^{\infty}}$ can be stated as an optimal control problem:

$\mbox{P1: }\left\\{\begin{array}[]{rl}\displaystyle\mbox{minimize}&\ \
\displaystyle\,\max_{t_{i}\leq t\leq t_{i+1}}\|u(t)\|_{\infty}\\\\[8.53581pt]
\mbox{subject to}&\ \ \dot{z}(t)=f(t,z(t))\,(1+u(t))\,,\mbox{ a.e. }t_{i}\leq
t\leq t_{i+1}\,,\\\\[5.69054pt] &\ \ z(t_{i})=y_{i},\ \
z(t_{i+1})=y_{i+1}\,,\end{array}\right.$

where $\|\cdot\|_{\infty}$ is the $\ell_{\infty}$-norm in $\mathbb{R}^{n}$.
Problem (P1), where $z(t)$ assumes the role of a state variable vector, can be
transformed into a more standard form to apply a maximum principle, as in
KayNoa2013 , as follows. Let $\alpha$ be a new parameter. Then one can rewrite
Problem (P1) equivalently as

$\mbox{P2: }\left\\{\begin{array}[]{rl}\displaystyle\mbox{minimize}&\ \
\alpha\\\\[8.53581pt] \mbox{subject to}&\ \
\dot{z}(t)=f(t,z(t))\,(1+u(t))\,,\mbox{ a.e. }t_{i}\leq t\leq
t_{i+1}\,,\\\\[5.69054pt] &\ \ z(t_{i})=y_{i},\ \
z(t_{i+1})=y_{i+1}\,,\\\\[5.69054pt] &\ \ |u_{j}(t)|\leq\alpha\,,\ \mbox{ a.e.
}t_{i}\leq t\leq t_{i+1}\,,\ j=1,\ldots,n\,.\end{array}\right.$

To apply the maximum principle, we define the Hamiltonian function for Problem
(P2):

$H(z,u,\psi,t):=\sum_{j=1}^{n}\,\psi_{j}\,f_{j}(t,z)\,(1+u_{j})\,,$

where $\psi:[t_{i},t_{i+1}]\to\mathbb{R}^{n}$ is referred to as the adjoint
(or costate) variable vector, such that

$\dot{\psi}_{j}=-\frac{\partial H}{\partial
z_{j}}=-\sum_{j=1}^{n}\,\psi_{j}\,\frac{\partial f_{j}}{\partial
z_{j}}(t,z)\,(1+u_{j})\,,\ j=1,\ldots,n\,.$ (3)

The maximum principle asserts that if $u$ is optimal then there exists a
nontrivial adjoint variable vector such that $u$ minimizes the Hamiltonian;
namely

$u(t)\in\operatorname*{argmin}_{|v|\leq\alpha}H(z(t),v,\psi(t),t)\,,\
\mbox{a.e. }t\in[t_{i},t_{i+1}]\,,$

i.e.,

$u_{j}(t)=\left\\{\begin{array}[]{ll}\ \ \alpha\,,&\mbox{if}\
\psi_{j}(t)\,f_{j}(t,z(t))<0\,,\\\\[2.84526pt] -\alpha\,,&\mbox{if}\
\psi_{j}(t)\,f_{j}(t,z(t))>0\,,\\\\[2.84526pt] \mbox{\small
undetermined}\,,&\mbox{if}\
\psi_{j}(t)\,f_{j}(t,z(t))=0\,.\end{array}\right.\quad j=1,\ldots,n\,.$ (4)

This implies that the components of the optimal control, or the optimal
backward error, will either be bang–bang (i.e., switching between the
constants $\alpha$ and $-\alpha$) or singular (when
$\psi_{j}(t)\,f_{j}(t,z(t))=0$ over some nontrivial subinterval). It is not
uncommon to encounter the error components to be of bang–bang or singular
types as preliminary studies show. However, this is a topic of further
research.

In this paper, we focus on the Dahlquist test equation, $\dot{y}=\lambda y$,
for which the optimality conditions derived above simplify drastically. For
Dahlquist with real $\lambda$, the adjoint equation (3) reduces to the scalar
ODE

$\dot{\psi}(t)=-\lambda\,(1+u(t))\,\psi(t)\,.$ (5)

Singular control can be ruled out for this case as follows. If the optimal
control is singular, then from (4) either $\psi(t)=0$ or $f(t,z(t))=z(t)=0$
over a subinterval $[s_{1},s_{2}]\subset[t_{i},t_{i+1}]$. Suppose that
$\psi(t)=0$ over $[s_{1},s_{2}]$. Then
$\dot{\psi}(t)=-\lambda\,(1+u(t))\,\psi(t)=0$ over $[s_{1},s_{2}]$, for any
$u$. However, this further means that $\psi(t)=0$ for all
$t\in[t_{i},t_{i+1}]$, which, being trivial, is not allowed by the maximum
principle. Now, suppose that $z(t)=0$ over $[s_{1},s_{2}]$. Then similarly
$z(t)=0$ for all $t\in[t_{i},t_{i+1}]$, which we can leave out as a very
special case with $y_{i}=0$.

Therefore the only possible optimal relative error is of bang–bang type. In
other words, $|\delta(t)|=\alpha$, i.e., the optimal relative error is of
constant magnitude.

Suppose that the control $u(t)=\alpha$ or $-\alpha$ on
$[t_{i},t_{i}+\varepsilon]$, with a small $\varepsilon>0$. Then (5) implies
that $\psi(t)$ does not change sign for any $t\in[t_{i},t_{i+1}]$. Thus,
either $u(t)=\alpha$ or $u(t)=-\alpha$ for all $t\in[t_{i},t_{i+1}]$. As a
result, the optimal interpolant is given by

$\dot{z}(t)=(1\pm\alpha)\,\lambda\,z(t)\,,$

whose solution can simply be obtained as

$z(t)=y_{i}\,e^{(1\pm\alpha)\,\lambda\,(t-t_{i})}\,.$

It follows that

$y_{i+1}=y_{i}\,e^{(1\pm\alpha)\,\lambda\,h}\,,$ (6)

or111Going from equation (6) to equation (7) is correct over the positive
reals but not quite right over $\mathbb{C}$: $\ln z$ means the principal
branch of the logarithm, with argument in $(-\pi,\pi]$. We will see how to
choose the correct complex branch in equation (9), in a later section.

$1\pm\alpha=\frac{1}{\lambda\,h}\,\ln\frac{y_{i+1}}{y_{i}}\,.$ (7)

Rearranging further yields the minimum value of the $L^{\infty}$-norm of the
relative error as

$\alpha=\left|1-\frac{1}{\lambda\,h}\,\ln\frac{y_{i+1}}{y_{i}}\right|.$

One may also look for a control defined over more than one piece of the
skeleton; this would allow the analysis to be used for multistep methods. In
this paper, however, we restrict attention to one-step methods.

Indeed for some problems the minimal $\delta$ is surprisingly easy to find:
see also corless2016variations for an investigation of Torricelli’s law using
this idea.

We will use this approach to re-examine the Dahlquist test problem,
$\dot{y}=\lambda y$. We will be able to find the interpolant giving the
optimal $\delta$ without difficulty using only elementary methods. Our
analysis will be valid for any one-step method. The results are quite
surprising.

##  3 Optimal Backward Error for $\dot{y}=\lambda y$

Without loss of generality suppose that our one-step numerical method gives,
for $\dot{y}=\lambda y$, $y(t_{n})=y_{n}$, and time-step $h$

$y_{n+1}=y_{n}+h\Phi(y_{n})=R(\mu)y_{n}\>,$

where $\mu=\lambda h\in\mathbb{C}$ and the time-step $h>0$. For instance, the
Explicit Euler method gives

$y_{n+1}=y_{n}+h\cdot\lambda y_{n}=(1+\mu)y_{n},$

while Implicit Euler gives

$y_{n+1}=y_{n}+h\cdot\lambda y_{n+1}$

or

$y_{n+1}=\frac{1}{1-\mu}y_{n}.$

A small table of $R(\mu)$ corresponding to various methods is given in Table
1. Note that Taylor series methods of order $p$ have

$R(\mu)=\sum_{k=0}^{p}\frac{\mu^{k}}{k!}.$

###  3.1 Classical stability regions

The classical stability regions are defined given the behaviour of numerical
methods on the Dahlquist test problem. The absolute stability region is
defined as the region in the complex plane where the value of $\mu=\lambda h$
guarantees that the Dahlquist test problem is uniformly bounded for all
forward time steps. Again the rationale for focusing on such a simple problem
is that it characterizes the local behaviour of a method on a scalar component
of a linearized version of any nonlinear problem.

For Runge-Kutta methods, the absolute stability region may be characterized in
terms of the quantity $R(\mu)$ as the region satisfying $|R(\mu)|\leq 1$.
Having this condition satisfied strictly ($<$) for those problems where the
exact solution decays (eigenvalues with negative real part) underlies the
criterion of A-stability. An A-stable method is any method that contains the
right-half plane in its stability region. This is an important criterion for
stiff problems, in the sense of problems that have eigenvalues with real parts
widely separated and negative, because an A-stable method will avoid stability
restrictions that force an extremely small time step when small time steps are
not required for an accurate solution. An exactly A-stable method is one that
has precisely the left-half plane, with the imaginary axis as boundary, as its
absolute stability region. The stability regions of the Euler method, the
implicit Euler method and the implicit midpoint rule are shown in figure 1.
The implicit midpoint rule is an example of an exactly A-stable method.

(a) Euler method

(b) Implicit midpoint rule

(c) Implicit Euler method

Figure 1: The classical stability regions for the one-stage $\theta$ methods
(a) $\theta=0$, (b) $\theta=\tfrac{1}{2}$, (c) $\theta=1$. The stability
region of (b) is the same for any exactly A-stable method.

Higher order methods have more complex and more interesting absolute stability
regions. The Runge-Kutta-Fehlberg (or RKF45) method includes a pair of Runge-
Kutta methods consisting of a fourth order method together with a fifth order
error estimator. The absolute stability regions for these methods are shown in
figure 2.

(a) Fourth order component of RKF45

(b) Fifth order component of RKF45

Figure 2: The classical stability regions for the components of the Fehlberg
method.

The shift to consider optimal residual stability regions entails a very
different perspective on numerical methods for stiff problems.

As is well-known, consideration of the relative forward error
$|R(\mu)e^{-\mu}-1|$ gives rise to the beautiful theory of Order Stars, which
considers the regions
$A_{+}:|R(\mu)e^{-\mu}|>1,A_{0}:|R(\mu)e^{-\mu}|=1,A_{-}:|R(\mu)e^{-\mu}|<1$.
For example, the three regions for the fifth-order method of RKF45 are shown
in figure 3. For a detailed account see iserles1991order . In hindsight, it is
not a surprise that a similar theory will arise from a consideration of
relative backward error.

method | $R(\mu)$
---|---
Euler | $1+\mu$
Backward Euler | $\frac{1}{1-\mu}$
Implicit Midpoint Rule | $\frac{1+\mu}{1-\mu}$
Taylor Series order $p$ | $\sum_{k=0}^{p}\frac{\mu^{k}}{k!}$
Table 1: Approximations to $e^{\mu}$ for some numerical methods Figure 3:
Order star for the fifth order component of RKF45. The $A_{-}$ region is
shaded, the $A_{0}$ boundary is dotted, and the $A_{+}$ region is the unshaded
part of the plane.

###  3.2 Elementary analysis

We now consider using only elementary tools to achieve the same result for the
Dahlquist test problem, with a one-step method that produces
$y_{n+1}=R(\mu)y_{n}$. Also without loss of generality consider only $t_{n}=0$
and $t_{n+1}=h$. We search for an optimal interpolant, which we call $z(t)$,
satisfying $z(0)=y_{n}$, $z(h)=y_{n+1}$ and

$\dot{z}(t)=\lambda(1+\delta(t))z(t),$ (8)

with $\delta(t)$ as small as possible. By constructing it explicitly we will
demonstrate its existence and its minimality for $\delta(t)$. We use the
$\infty$-norm to measure the size of $\delta$:

$\|\delta(t)\|_{\infty}=\max_{0\leq t\leq h}|\delta(t)|.$

Rearranging, and assuming $z(s)\neq 0$ anywhere in $0\leq s\leq h$,

$\frac{\dot{z}(t)}{z(t)}-\lambda=\lambda\delta(t),$

so

$\int_{0}^{h}\frac{\dot{z}(s)}{z(s)}ds-\lambda\int_{0}^{h}ds=\lambda\int_{0}^{h}\delta(s)ds,$

or

$\ln z(h)-\ln z(0)-\lambda h=\lambda\int_{0}^{h}\delta(s)\,ds-2\pi ik\>,$ (9)

where the integer $k$ is determined by the (as yet unknown) number of times
the path $z(s)$ winds around $z=0$. We will see later that $k=0$ often,
particularly for small step sizes, but there are important cases where $k\neq
0$, especially for higher-order methods. By adjusting $k$ if necessary we get

$\ln_{k}\left(\frac{z(h)}{z(0)}\right)-\mu=\lambda\int_{0}^{h}\delta(s)\,ds,$

where we have used David Jeffrey’s compact notation $\ln_{k}z$ for $\ln z+2\pi
ik$.

Since $\frac{z(h)}{z(0)}=\frac{y_{n+1}}{y_{n}}=R(\mu)$ we have

$\ln_{k}R(\mu)-\mu=\lambda\int_{0}^{h}\delta(s)ds.$

Taking absolute values and using the triangle inequality,

$\left|\ln_{k}R(\mu)-\mu\right|\leq|\lambda|\,\left|\int_{0}^{h}ds\right|\,\|\delta(t)\|_{\infty}$

or

$|\ln_{k}R(\mu)-\mu|\leq|\mu|\,\|\delta\|_{\infty}$

or

$\|\delta\|_{\infty}\geq\left|\frac{\ln_{k}R(\mu)}{\mu}-1\right|,$

where we ignore the uninteresting case of $\mu=0$.

This fundamental inequality gives a lower bound on any backward error
$\delta(t)$ capable of taking $y_{n}$ to $y_{n+1}$. We now show that this
lower bound is achieved, if we choose $\delta(t)$ to be constant. Suppose

$\delta=\frac{\ln_{k}R(\mu)}{\mu}-1\>.$ (10)

Then $\dot{z}=\lambda(1+\delta)z$ is well-defined and indeed

$z=e^{\lambda(1+\delta)t}z(0)=e^{\lambda\frac{\ln_{k}R(\mu)}{\mu}t}y_{n}.$

When $t=h$,

$\displaystyle z(h)$ $\displaystyle=$ $\displaystyle
e^{\lambda\frac{\ln_{k}R(\mu)}{\mu}h}y_{n}$ (11) $\displaystyle=$
$\displaystyle e^{\mu\frac{\ln_{k}R(\mu)}{\mu}}y_{n}$ $\displaystyle=$
$\displaystyle e^{\ln_{k}R(\mu)}y_{n}$ $\displaystyle=$ $\displaystyle
R(\mu)y_{n}$ $\displaystyle=$ $\displaystyle y_{n+1}$

as desired, for any choice of $k$.

Since the set
$\left\\{\left|\frac{\ln_{k}R(\mu)}{\mu}-1\right|\,:\,k\in\mathbb{Z}\right\\}$
is countable and nonnegative, it has a least member. We will use the $k$ that
picks out the least member of this set.

Lemma: To minimize $\left|\frac{\ln_{k}R(\mu)}{\mu}-1\right|$ over choices of
branch cut, we must choose

$k=\overline{\mathcal{K}}:=\left[\frac{\text{Im}(\mu-\ln
R(\mu))}{2\pi}\right],$

where $[a]$ is the nearest integer to $a$; in case of a tie, either integer
above or below will do.

###### Proof

We seek a choice of $k$ that minimizes
$\left|\frac{\ln_{k}R(\mu)}{\mu}-1\right|$, i.e., that minimizes
$|\varepsilon|$ for

$\varepsilon:=\frac{\ln_{k}R(\mu)}{\mu}-1=\frac{\ln R(\mu)+2\pi ik}{\mu}-1.$

Putting this over a common denominator, we get

$\varepsilon=\frac{\ln R(\mu)-\mu+2\pi ik}{\mu}\>.$ (12)

For the magnitude of $\varepsilon$ to be as small as possible, we choose $k$.
But this alters only the imaginary part of the numerator, because $k$ is an
integer and therefore $2\pi ik$ is purely imaginary. In order to make this
imaginary term as small as possible, it follows that we must have that the
integer $k$ cancels as much as possible of the imaginary part from equation
(12), or

$k=\overline{\mathcal{K}}:=\left[\frac{\text{Im}(\mu-\ln
R(\mu))}{2\pi}\right]\>.$ (13)

This formula is very reminiscent of the unwinding number from
jeffrey1996unwinding , but is different in detail.

###### Remark 1

We will be able to examine the optimal backward error of a very large class of
methods, essentially all one-step methods, with just formula (10) with the
integer $k$ specified by formula (13).

###### Remark 2

If ever an interpolant hits $z(s)=0$, from $\dot{z}=\lambda(1+\delta)z$ we
would have $z\equiv 0$ from then on, so unless $y_{n+1}=0$ we could not reach
it. Moreover, if $y_{n+1}=0$ then $\|\delta\|_{\infty}=\infty$ is necessary (a
small absolute $\Delta(t)$ is possible in that case, but the relative error
must go to infinity). As we saw in the general optimal control approach, this
corresponds to the case when the control is singular.

##  4 Results: Optimal backward error for some methods

We now examine properties of the optimal backward error along with contour
plots of the optimal $\infty$-norm optimal residuals. Note $|\delta|=1$
corresponds to $100$% error; $|\delta|$ larger than this means we are solving
a totally different equation. In the absence of systematic structure
preservation, we may regard any likeness of the solution to what we want as
coincidence.

Taking $|\delta|=0.05$ is analogous to the “$95\%$ confidence limit”. If
$|\delta|\leq 0.05$, we have the exact solution to a problem within $5$% of
the one we wanted to solve. We will also be concerned with the asymptotic
limit $|\delta|=\varepsilon$ as $\varepsilon\rightarrow 0$.

(a) Euler method

(b) Implicit midpoint rule

(c) Implicit Euler method

Figure 4: The $\infty$-norm optimal residual stability regions from equation
(10) for the one-stage $\theta$ methods (a) $\theta=0$, (b)
$\theta=\tfrac{1}{2}$, (c) $\theta=1$. In each image, the classical stability
region is shaded and the $\infty$-norm optimal residual stability region is
coloured, with contours at 5% intervals. Regions where $|\delta|>1$ are
coloured white. Note the significantly larger region where $\delta=0.05$ for
the second-order method with $\theta=\tfrac{1}{2}$. In essence this is why
second-order methods are more efficient than first-order methods, for a given
accuracy. Perhaps most important, the regions where $|\delta|>1$, i.e. those
regions where the problem solved is more than $100\%$ different to the problem
we wanted to solve, has a nontrivial intersection with the classical stability
region.

###  4.1 Theta methods

The rational approximation that arises from the theta-method

$y_{n+1}=y_{n}+\mu\left((1-\theta)y_{n}+\theta y_{n+1}\right)$ (14)

is

$R(\mu)=\frac{1+(1-\theta)\mu}{1-\theta\mu}\>.$ (15)

For $\theta=0$ we have explicit Euler, for $\theta=1$ we have implicit Euler,
and for $\theta=1/2$ we have the implicit midpoint rule. The contours of
$\delta$ for each of these three $\theta$ are plotted in figure 4. We see that
for explicit Euler there is a substantial portion of the classical stability
region $|\mu+1|<1$ where the computed solution will have $|\delta|>1$. This
means that even though the solution will decay, it must be the solution of a
problem that is more than $100\%$ different to the original problem. Likewise,
a substantial portion of the right-half plane will have $|\delta|>1$ for
implicit Euler; this is why such numerical methods can suppress (actual) chaos
corless1991numerical .

In contrast, the implicit midpoint rule (with $\theta=1/2$) doesn’t have such
large $100\%$ error zones (although it does, near the portions of the real
axis where $|\mu|>2$). Instead it has an interestingly larger region where
$|\delta|=0.05$ than either explicit or implicit Euler does. This is because
the method is second-order, and the substantially larger region “of $95\%$
confidence,” if you like, is a reflection of the value of a second-order
method. Higher order methods really can be more efficient, and this picture
shows why.

(a) Fourth order component of RKF45

(b) Fifth order component of RKF45

Figure 5: The $\infty$-norm optimal residual stability regions from equation
(10) for the components of the RK-Fehlberg method. In each image, the
classical stability region is shaded and the $\infty$-norm optimal residual
stability region is coloured, with contours at 5% intervals.

###  4.2 RKF45

In figure 5 we find the optimal stability regions for each member of the RKF45
pair. We see substantial regions in the left half plane where the optimal
residual must be larger than $1$ in magnitude; that is, where the backward
error must be larger than $100\%$. The regions where $|\delta|\leq 0.05$ are
comparable for each member of the pair, showing that if one method is
accurate, the other is likely to be as well. This information is complementary
to that of the classical stability regions for the pair in figure 2.

(a) Larger $g$ SDIRK method

(b) Smaller $g$ SDIRK method

Figure 6: The $\infty$-norm optimal residual stability regions from equation
(10) for the two singly-diagonally-implicit Runge-Kutta (SDIRK) methods of
order $3$. In each image, the classical stability region is shaded and the
$\infty$-norm optimal residual stability region is coloured, with contours at
5% intervals. Classical stability analysis suggests that the larger $\gamma$
method (left) is more stable. This analysis suggests that the smaller $\gamma$
method (right) is more accurate.

###  4.3 A third order SDIRK method

Singly-Diagonally Implicit Runge-Kutta (SDIRK) methods are important and
efficient implicit solvers for stiff problems. In wanner1991solving we find a
detailed derivation of a third-order SDIRK method that contains the diagonal
parameter $\gamma$. Two values of $\gamma$ make the method $3$rd order. It is
argued there that one value of $\gamma$ is to be preferred over the other, as
making the method have a larger (much larger) classical stability region.
However, the optimal backward error diagram in figure 6 shows that it is the
other value of $\gamma$ that is potentially more accurate.

Of course, real life is more complicated than either of these images would
show. See the massive review kennedy2016diagonally for a more nuanced picture.

Figure 7: Lanczos $\tau$-method Figure 8: The $\infty$-norm optimal residual
stability regions from equation (10) for the Lanczos Tau method with $n=1$.
The $\infty$-norm optimal residual stability region is shaded in each image,
with contours at 5% intervals.

###  4.4 The Lanczos $\tau$-method

We have included one figure (figure 8) with the results from solving
$\dot{y}=\lambda y$, $y(0)=y_{0}$, on $0\leq t\leq h$, by using the Lanczos
$\tau$-method (see for instance lanczos1988applied or ortiz1969tau ). For ease
of reference, we include a brief description of the method here. One starts
with a Chebyshev expansion for $\dot{y}$, namely

$\dot{y}(t)=\sum_{k=0}^{n}c_{k}T_{k}(\theta)\>,$ (16)

with $\theta=-1+2t/h$. Integrating $T_{0}$ with respect to $\theta$ gives
$T_{1}(\theta)$, integrating $T_{1}(\theta)$ gives
$T_{2}(\theta)/4+T_{0}(\theta)/4$, and thereafter (ignoring constants of
integration)

$\int
T_{k}(\theta)\,d\theta=\frac{1}{2k+2}T_{k+1}(\theta)-\frac{1}{2k-2}T_{k-1}(\theta)\>.$
(17)

Lanczos chose to expand the derivative of the unknown, and then integrate
using these simple formulas, because this was simpler to do by hand. The
details of the process are now of course automatable in a variety of ways.
Then integration of equation (16) gives an expression for $y(t)$, in terms of
Chebyshev polynomials of degree $n+1$ and below. One then uses the initial
condition $y(0)=y_{0}$ to identify the constant introduced on integration. One
subtracts $\lambda$ times this expression from the expression for $\dot{y}$,
and sets the coefficients of $T_{k}(\theta)$ to zero for $k=0$ up to $k=n$.
This leaves here a residual term containing $T_{n+1}(\theta)$. If we do this
using $n=1$ we find that our Chebyshev interpolant is

$z(t)=y_{0}\left(\frac{16}{(4-\mu)^{2}}T_{0}(\theta)+\frac{8\mu}{(4-\mu)^{2}}T_{1}(\theta)+\frac{\mu^{2}}{(4-\mu)^{2}}T_{2}(\theta)\right)\>,$
(18)

which can be simplified at $t=h$ to

$y_{1}=\frac{1+\mu/2+\mu^{2}/16}{1-\mu/2+\mu^{2}/16}\,y_{0}\>.$ (19)

This identifies $R(\mu)=(1+\mu/4)^{2}/(1-\mu/4)^{2}$ and from there we can
identify the optimal $\delta$, which is plotted in figure 8. Similar graphs
can be produced for larger $n$.

###### Remark 3

The Lanczos $\tau$-method is quite close in spirit to finding the minimal
$\infty$-norm residual, because Chebyshev expansions are, as is well-known,
near-minimal for real $\mu$ in this sense. One important difference here is
that because $\lambda\in\mathbb{C}$ we are concerned with the region in the
complex plane where $|\delta|$ is small (less than $1$ certainly and less
than, say, $0.05$ in practice). Comparison of the optimal backward error on
the real interval $(0,h)$ to the size of the Chebyshev residual from the
$\tau$ method (not shown here) shows that the Chebyshev method is less than a
factor of five worse. Reluctantly, however, we do not pursue the Lanczos
$\tau$-method further here.

(a) Order 2 Taylor method

(b) Order 4 Taylor method

(c) Order 8 Taylor method

(d) Order 16 Taylor method

Figure 9: The $\infty$-norm optimal residual stability regions from equation
(10) for the higher order Taylor methods of order (a) $2$, (b) $4$, (c) $8$
and (d) $16,$. The $\infty$-norm optimal residual stability region is shaded
in each image, with contours at 5% intervals.

(a) 2,2 Padé method

(b) 4,4 Padé method

(c) 8,8 Padé method

(d) 16,16 Padé method

Figure 10: The $\infty$-norm optimal residual stability regions from equation
(10) for the higher order Padé methods of signature (a) $2,2$, (b) $4,4$, (c)
$8,8$ and (d) $16,16$. The $\infty$-norm optimal residual stability region is
shaded in each image, with contours at 5% intervals.

###  4.5 Taylor series methods and Padé methods

Taylor series methods, including implicit Taylor series methods, and their
generalization to Hermite-Obreshkoff methods, remain of interest for practical
problems. See for instance neher2007taylor . For the Dahlquist test problem,
all of these methods lead to $R(\mu)$ being a Padé approximant to $\exp(\mu)$.
Figures 10 and 9 were generated using Maple’s built-in function for
computation of Padé approximants geddes1979symbolic . Notice that as the order
of the method increases, the size of the area enclosed by the $\delta=0.05$
contour increases. We remark that the unwinding number from formula (13) was
necessary to get this large central region correct, for high-order methods.
Without the correct branch of logarithm chosen, the central region only had a
vertical width of $2\pi$, corresponding to the range of the principal branch
of logarithm.

###  4.6 Asymptotic results

Formula 10 is not an asymptotic result. It does not rely on the time-step $h$
being small. The optimal $\delta$ is the optimal $\delta$, with almost no
reservations or caveats. It exists for the Dahlquist test problem for $O(1)$
intervals around $h=0$, regardless of the one-step numerical method used, so
long as the optimal interpolant does not go through $0$. Nonetheless there are
some interesting connections to asymptotic results. First, if the underlying
method has forward error $O(h^{p})$ as $h\to 0$, that is, it is a $p$th order
method, then the residual $\delta$ will also go to zero like $h^{p}$. That is,
$\delta=O(h^{p})$. This is in contrast to the distinction between so-called
local error (which is $O(h^{p+1})$) and forward error (also called global
error). For instance, if

$R(\mu)=\frac{(1+\mu/4)^{2}}{(1-\mu/4)^{2}}\>,$ (20)

which it is for the Lanczos $\tau$-method with $n=1$, then

$\delta=\frac{\ln_{k}(R(\mu))}{\mu}-1={\frac{1}{48}}{\mu}^{2}+{\frac{1}{1280}}{\mu}^{4}+O\left({\mu}^{6}\right)$
(21)

showing that this particular instance of the Lanczos $\tau$-method acts on the
Dahlquist test problem as if it was a second-order method (in fact it is in
general a second order method, but all this expansion does is illustrate that
it’s second-order on this problem). Note that for small $h$,
$R(\mu)=1+O(\mu^{2})$ and hence the unwinding number is $k=0$.

Because the residual is connected to the forward error by the Gröbner-Alexeev
nonlinear variation-of-constants formula hairer1993solving , the forward error
is also $O(h^{2})$ for this method.

Similar expansions for other methods confirm their numerical order.

##  5 Discussion

The classical stability analysis using the Dahlquist test problem is
informative in that it gives the fundamental reason for loss of stability on
taking too large a step size for a stiff problem. Yet the classical analysis
neglects the case of eigenvalues with positive real part (which are indeed
important in applications) and furthermore says little about accuracy: decay
is a qualitative feature, not a quantitative one, and one can have decay and
yet be $100\%$ wrong. As we see in the figures for some common methods,
computing the optimal backward error gives some complementary information,
namely the relative size of the figure enclosed by the $0.05$ contour level:
the larger the area, the larger the time-step that can be taken. We see indeed
that this region (and indeed even the $100\%$ region) can be much smaller than
the classical stability region. This indicates that even though the method may
be stable for time steps inside the classical stability region but outside the
$100\%$ contour, it certainly is not accurate, solving a problem that is more
than $100\%$ different to the one that was intended to be solved.

This observation is in some sense not new. Experienced numerical analysts knew
this, and knew that one had to take small time steps for accuracy anyway (when
it was needed). Still, this quantitative assessment of optimal relative
backward accuracy is indeed new. It seems possible that one might choose a
different method to solve some problems, using this criterion, than one might
choose using the classical stability criterion.

The analysis here applies to any one-step method.

###### Remark 4

Since the optimal $\delta$ is constant (for any one-step method), we have the
curious observation that constant-stepsize solution using a one-step method
gives the exact solution of $\dot{y}=\Lambda y$ where
$\Lambda=\lambda(1+\delta)=\lambda(\ln_{k}R(\mu)/\mu)$; that is, the optimal
interpolant is the exact solution of not only a nearby problem, but a nearby
problem of the same kind. That is, we have not only an optimal backward error,
but an optimal structured backward error.

That constant-stepsize one-step methods solved $\dot{y}=\Lambda y$ was known;
what was not known was that this solution has the optimality property derived
in this paper.

###### Remark 5

If the optimal backward error is large, then the numerical method has
necessarily solved a very different problem to the one intended. This is a
sure and certain indication that the underlying numerical method, that
generated the skeleton, has failed. In particular, when a numerical method
introduces spurious chaos into a nonchaotic system, the optimal backward error
must be too large. A more frequent failure detected (for nonchaotic systems)
will be when the automatically-chosen time step sizes are too large; because
the analysis of this paper and the more general paper in preparation do not
rely on the asymptotic limit as the mean stepsize goes to zero, detection of
failure is sure and certain.

##  6 Concluding Remarks

John C. Butcher was correct when, critiquing an earlier incarnation of this
idea, he said that the classical stability theory does indeed explain many
features of the behaviour of many numerical methods on stiff problems. It is a
highly-successful theory. The refinement of using Order Stars is also highly
successful, in that new proofs of classical order barriers were obtained using
it (for instance).

However, this optimal backward error idea for $y_{n+1}=R(\mu)y_{n}$ with
$\mu=h\lambda$, namely formula (10) with the unwinding number $k$ from formula
(13),

$\delta=\frac{\ln_{k}R(\mu)}{\mu}-1\>,$ (22)

with

$k=\overline{\mathcal{K}}:=\left[\frac{\text{Im}(\mu-\ln
R(\mu))}{2\pi}\right]\>,$ (23)

also explains, quantitatively, some features of the numerical solution of
stiff problems by one-step methods. In particular, it explains clearly why one
might sometimes wish to take smaller timesteps than strictly necessary for
stability reasons. That is, any adaptive step-size control that looks at
accuracy must be affected by these regions. We note that asymptotically as
$h\to 0$, the optimal residual approaches the local error per unit step, as
explained by a theorem of Stetter (this is cited and extended in
corless2013graduate ). There are other puzzles that this approach might help
to explain. We look forward, for instance, to finding an example problem for
which the smaller $\gamma$ SDIRK method performs better than the larger
$\gamma$ SDIRK method, as is predicted to exist by figure 6.

### Maple Worksheet

The Maple worksheet written by R.H.C. Moir used to produce the graphs in this
paper is at
publish.uwo.ca/~rcorless/Dahlquist/ResidualAnalysis-3D-Contours.mw.

### Acknowledgements

Part of this work was supported by the Natural Sciences and Engineering
Research Council of Canada. Thanks are also due to the Rotman Institute of
Philosophy, the Fields Institute for Research in the Mathematical Sciences,
the Ontario Research Center for Computer Algebra, and the University of South
Australia which supported a visit of the first author to the second. Special
thanks are owed to John C. Butcher for his hosting the beautiful ANODE
conferences over the years, and his long-term support of research in the
Runge-Kutta world.

## References

  * (1) Cao, Y., Petzold, L.: A posteriori error estimation and global error control for ordinary differential equations by the adjoint method.  SIAM Journal on Scientific Computing 26(2), 359–374 (2004)
  * (2) Corless, R.M.: Error backward.  Contemporary Mathematics 172, 31–31 (1994)
  * (3) Corless, R.M.: What good are numerical simulations of chaotic dynamical systems?  Computers & Mathematics with Applications 28(10-12), 107–121 (1994)
  * (4) Corless, R.M., Essex, C., Nerenberg, M.: Numerical methods can suppress chaos.  Physics Letters A 157(1), 27–36 (1991)
  * (5) Corless, R.M., Fillion, N.: A graduate introduction to numerical methods.  Springer (2013)
  * (6) Corless, R.M., Jankowski, J.E.: Variations on a theme of Euler.  SIAM Review 58(4), 775–792 (2016)
  * (7) Enright, W.: Continuous numerical methods for ODEs with defect control.  Journal of Computational and Applied Mathematics 125(1), 159–170 (2000)
  * (8) Enright, W., Seward, W.: Achieving tolerance proportionality in software for stiff initial-value problems.  Computing 42(4), 341–352 (1989)
  * (9) Enright, W.H.: A new error-control for initial value solvers.  Applied Mathematics and Computation 31, 288–301 (1989)
  * (10) Enright, W.H., Hayes, W.B.: Robust and reliable defect control for Runge–Kutta methods.  ACM Transactions on Mathematical Software (TOMS) 33(1), 1 (2007)
  * (11) Enright, W.H., Higham, D.J.: Parallel defect control.  BIT Numerical Mathematics 31(4), 647–663 (1991)
  * (12) Geddes, K.O.: Symbolic computation of Padé approximants.  ACM Transactions on Mathematical Software (TOMS) 5(2), 218–233 (1979)
  * (13) Hairer, E., Nørsett, S.P., Wanner, G.: Solving ordinary differential equations. i (1993)
  * (14) Hubbard, J., West, B.H.: Differential equations: a dynamical systems approach, vol. 5.;5, 18;18.;.  Springer-Verlag, New York (1991)
  * (15) Iserles, A., Norsett, S.P.: Order Stars: Theory and Applications, vol. 2.  CRC Press (1991)
  * (16) Jeffrey, D.J., Hare, D.E., Corless, R.M.: Unwinding the branches of the Lambert $W$ function.  Mathematical Scientist 21(1), 1–7 (1996)
  * (17) Kaya, C.Y., Noakes, J.L.: Finding interpolating curves minimizing ${L}^{\infty}$ acceleration in the Euclidean space via optimal control theory.  SIAM J. Control Optim. 51, 442–464 (2013)
  * (18) Kennedy, C.A., Carpenter, M.H.: Diagonally implicit Runge–Kutta methods for ordinary differential equations. a review (2016).  URL https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20160005923.pdf
  * (19) Lanczos, C.: Applied analysis.  Dover (1988)
  * (20) Moir, R.H.: Reconsidering backward error analysis for ordinary differential equations.  Master’s thesis, The University of Western Ontario. (2010)
  * (21) Muir, P., Pancer, R., Jackson, K.R.: PMIRKDC: a parallel mono-implicit Runge–Kutta code with defect control for boundary value ODEs.  Parallel Computing 29(6), 711–741 (2003)
  * (22) Neher, M., Jackson, K.R., Nedialkov, N.S.: On Taylor model based integration of ODEs.  SIAM Journal on Numerical Analysis 45(1), 236–262 (2007)
  * (23) Ortiz, E.L.: The tau method.  SIAM Journal on Numerical Analysis 6(3), 480–492 (1969)
  * (24) Söderlind, G., Jay, L., Calvo, M.: Stiffness 1952–2012: Sixty years in search of a definition.  BIT Numerical Mathematics 55(2), 531–558 (2014)
  * (25) Wanner, G., Hairer, E.: Solving ordinary differential equations ii.  Stiff and Differential-Algebraic Problems (1991)

##