#  Differentiating by prime numbers

Jack Jeffries  The author is an assistant professor of mathematics at the
University of Nebraska. His email address is<EMAIL_ADDRESS>His work
is supported by NSF CAREER Award DMS-2044833.

##

It is likely a fair assumption that you, the reader, are not only familiar
with but even quite adept at differentiating by $x$. What about
differentiating by 13? That certainly didn’t come up in my calculus class!
From a calculus perspective, this is ridiculous: are we supposed to take a
limit as 13 changes?

One notion of differentiating by 13, or any other prime number, is the notion
of _$p$ -derivation_ discovered independently by Joyal [MR789309] and Buium
[MR1387233]. $p$-derivations have been put to use in a range of applications
in algebra, number theory, and arithmetic geometry. Despite the wide range of
sophisticated applications, and the fundamentally counterintuitive nature of
the idea of differentiating by a number, $p$-derivations are elementary to
define and inviting for exploration.

In this article, we will introduce $p$-derivations and give a few basic ways
in which they really do act like derivatives by numbers; our hope is that you
will be inspired and consider adding $p$-derivations to your own toolkit!

##  $p$-derivations on $\mathbb{Z}$

First we want to discuss differentiating one number $n$, by another, $p$;
i.e., what we will call $p$-derivations on $\mathbb{Z}$. Before we succeed, we
need to abandon the notion of derivative as a limit as the input varies by a
small amount: the thing $p$ that we are differentiating by does not vary, and
the thing $n$ that we are differentiating does not even have an input!
Instead, we take a little inspiration from elementary number theory.

Let $p$ be a prime number. By Fermat’s little theorem, for any integer $n$, we
have

$n\equiv n^{p}\mod\ p,$

so we can divide the difference $n-n^{p}$ by $p$. The starting point of our
journey is that not only _can_ we divide by $p$ here, but we _should_. The
$p$-derivation on $\mathbb{Z}$ is the result of this process. Namely:

###### Definition 1

For a prime number $p$, the _$p$ -derivation_ on $\mathbb{Z}$ is defined as
the function ${\delta_{p}:\mathbb{Z}\to\mathbb{Z}}$ given by the formula

$\delta_{p}(n)=\frac{n-n^{p}}{p}.$

So, in particular, there is the $2$-derivation on $\mathbb{Z}$ and the
$13$-derivation on $\mathbb{Z}$ given respectively by

$\delta_{2}(n)=\frac{n-n^{2}}{2}\ \ \text{and}\ \
\delta_{13}(n)=\frac{n-n^{13}}{13}.$

Let’s plug in a few values:

$\begin{array}[]{ c | c c c }n&\delta_{2}(n)&\delta_{3}(n)&\delta_{5}(n)\\\ \hline\cr\vdots&\vdots&\vdots&\vdots\\\ -4&-10&20&204\\\ -3&-6&8&48\\\ -2&-3&2&6\\\ -1&-1&0&0\\\ 0&0&0&0\\\ 1&0&0&0\\\ 2&-1&-2&-6\\\ 3&-3&-8&-48\\\ 4&-6&-20&-204\\\ 5&-10&-40&-624\\\ 6&-15&-70&-1554\\\ \vdots&\vdots&\vdots&\vdots\end{array}$

A quick look at this table suggests a few observations, easily verified from
the definition:

  * •

Numbers are no longer “constants” in the sense of having derivative zero, but
at least $0$ and $1$ are.

  * •

These functions are neither additive nor multiplicative, e.g.:

$\delta_{p}(1)+\delta_{p}(1)\neq\delta_{p}(2),$
$\delta_{p}(1)\delta_{p}(2)\neq\delta_{p}(2).$

  * •

$\delta_{p}$ is an odd function, at least for $p\neq 2$.

  * •

The outputs of $\delta_{2}$ are just the negatives of the triangular numbers.

We might also note that the outputs are very large in absolute value, and
think that this operation is simply making a mess of our numbers. However,
something more informative occurs if we think about largeness of the outputs
from the point of view of $p$, namely, the _$p$ -adic order_ of $n$—the number
of copies of $p$ in its prime factorization. Writing $n=p^{a}m$ with
$\gcd(m,p)=1$, if $a>0$, we get

$\delta_{p}(p^{a}m)=\frac{p^{a}m-(p^{a}m)^{p}}{p}=p^{a-1}m(1-p^{ap-a}m^{p-1}).$

Since $p\geq 2$ and $a\geq 1$, we must have $ap-a\geq 1$, so $p$ does not
divide $1-p^{ap-a}m^{p-1}$. In particular, the $p$-derivation decreases the
$p$-adic order of a multiple of $p$ by exactly one. This leads to our first
comparison with old-fashioned $\frac{d}{dx}$:

######  Comparison 1 (Order-decreasing property)

  * •

If $f\in\mathbb{R}[x]$ is a polynomial and $x=r$ is a root of $f$ of
multiplicity $a>0$, then $x=r$ is a root of the polynomial
$\frac{d}{dx}(f(x))$ of multiplicity $a-1$.

  * •

If $n$ is an integer and $p$ is a prime factor of $n$ of multiplicity $a>0$,
then $p$ is a prime factor of the integer $\delta_{p}(n)$ of multiplicity
$a-1$.

In particular, if $r$ is a simple root or $p$ is a simple factor, then it is
no longer a root or factor of $\frac{d}{dx}(f(x))$ or $\delta_{p}(n)$
respectively.

Let’s check this against our table: the numbers $-2,2$, and $6$ that were
divisible by $2$ but not $4$ result in odd numbers when we apply $\delta_{2}$,
whereas $\pm 4$ returned even numbers no longer divisible by $4$. Note that
this order-decreasing property says nothing about what happens when you apply
$\delta_{2}$ to an odd number, and indeed, based on the table we observe that
even and odd numbers can result. You can convince yourself that

$\delta_{2}(n)\ \text{is}\ \begin{cases}\text{even}\ \text{if}\ n\equiv
0,1\mod 4\\\ \text{odd}\ \text{if}\ n\equiv 2,3\mod 4.\end{cases}$

We’ve observed already that these $p$-derivations on $\mathbb{Z}$ are not
additive. This can be a bit unsettling for those of us (like myself) who are
usually accustomed to the luxury of additive operators. However, any function
satisfying the order-decreasing property of $\delta_{p}$ above must not be
additive, since an additive function has to take multiples of $p$ to multiples
of $p$. However, the error term can be made concrete:

$\displaystyle\delta_{p}(m+n)-(\delta_{p}(m)+\delta_{p}(n))$
$\displaystyle=\frac{m^{p}+n^{p}-(m+n)^{p}}{p}$
$\displaystyle=-\sum_{i=1}^{p-1}\frac{\binom{p}{i}}{p}m^{i}n^{p-i}.$

All of the binomial coefficients $\binom{p}{i}$ appearing above are multiples
of $p$, so this expression is, given a particular value of $p$, a particular
polynomial in $m$ and $n$ with integer coefficients; let’s call it
$C_{p}(m,n)$ for convenience. This gives us the following “sum rule” for
$\delta_{p}$:

$\delta_{p}(m+n)=\delta_{p}(m)+\delta_{p}(n)+C_{p}(m,n).$ ($+$)

Products satisfy a rule with a similar flavor:

$\delta_{p}(mn)=m^{p}\delta_{p}(n)+n^{p}\delta_{p}(m)+p\delta_{p}(m)\delta_{p}(n).$
($\times$)

The fact that we have rules to break things down into sums and products gives
the basis for another comparison with old-fashioned $\frac{d}{dx}$:

######  Comparison 2 (Sum and product rules)

  * •

For polynomials $f(x),g(x)$, one can compute each of $\frac{d}{dx}(f+g)$ and
$\frac{d}{dx}(fg)$ as a fixed polynomial expression in the inputs
$f,g,\frac{d}{dx}(f),\frac{d}{dx}(g)$, namely
$\frac{d}{dx}(f+g)=\frac{d}{dx}(f)+\frac{d}{dx}(g)$ and
$\frac{d}{dx}(fg)=f\frac{d}{dx}(g)+g\frac{d}{dx}(f)$.

  * •

For integers $m,n$, one can compute each of ${\delta_{p}(m+n)}$ and
$\delta_{p}(mn)$ as a fixed polynomial expression in the inputs
$m,n,\delta_{p}(m),\delta_{p}(n)$, namely ($+$ ‣ $p$-derivations on
$\mathbb{Z}$) and ($\times$ ‣ $p$-derivations on $\mathbb{Z}$).

We might pause to ask whether we could have hoped for a simpler way to
differentiate by 13. If we want Comparsion 2 to hold, then the following
theorem of Buium provides a definitive answer.

######  Theorem 1 (Buium [MR1482984])

Any function ${\delta:\mathbb{Z}\to\mathbb{Z}}$ that satisfies

  * •

a sum rule $\delta(m+n)=S(m,n,\delta(m),\delta(n))$ for some polynomial $S$
with integer coefficients

  * •

a product rule $\delta(mn)=P(m,n,\delta(m),\delta(n))$ for some polynomial $P$
with integer coefficients

is of the form

$\delta(n)=\pm\frac{n-n^{p^{e}}}{p}+f(n)$

for some prime integer $p$, positive integer $e$, and polynomial $f$ with
integer coefficients.

That is, any function satisfying a sum rule and a product rule is a mild
variation on a $p$-derivation.

With the properties of $p$-derivations we have so far, we can recreate
analogues of some familiar aspects of calculus. For example, from the product
rule ($\times$ ‣ $p$-derivations on $\mathbb{Z}$) and a straightforward
induction, we obtain a power rule:

$\delta_{p}(n^{a})=\sum_{i=1}^{a}\binom{a}{i}p^{i-1}\delta_{p}(n)^{i}n^{(a-i)p}.$

Note that the $i=1$ term in the sum above, $an^{(a-1)p}\delta_{p}(n)$, looks a
bit like the power rule for usual derivatives. If we allow ourselves to extend
$\delta_{p}$ to a map on $\mathbb{Q}$, then we get an analogue of the quotient
rule:

$\delta_{p}\left(\frac{m}{n}\right)=\frac{n^{p}\delta_{p}(m)-m^{p}\delta_{p}(n)}{n^{2p}+pn^{p}\delta_{p}(n)}.$

Of the main cast of characters in a first class on derivatives, perhaps the
most conspicuous one missing at this point is the chain rule. Since there is
no way to compose a number with a number, we will need a notion of
$p$-derivations for functions to state a sensible analogue of the chain rule.

##  $p$-derivations for general commutative rings

One can define $p$-derivations for commutative rings with $1$.

###### Definition 2

Let $R$ be a commutative ring with $1$ and $p$ a prime integer. A _$p$
-derivation_ on $R$ is a function $\delta:R\to R$ such that
$\delta(0)=\delta(1)=0$ and $\delta$ satisfies the sum rule $(+)$ and the
product rule $(\times)$ above; i.e., for all $r,s\in R$,

$\delta_{p}(r+s)=\delta_{p}(r)+\delta_{p}(s)-\sum_{i=1}^{p-1}\frac{\binom{p}{i}}{p}r^{i}s^{p-i}$
($+$)

and

$\delta_{p}(rs)=r^{p}\delta_{p}(s)+r^{p}\delta_{p}(s)+p\delta_{p}(r)\delta_{p}(s).$
($\times$)

Evidently, the functions $\delta_{p}$ we defined on $\mathbb{Z}$ above are
$p$-derivations. In fact, for a fixed $p$, a simple induction and the sum rule
show that for any $p$-derivation $\delta$ on a ring $R$ and any $n$ in the
prime subring (image of $\mathbb{Z}$) of $R$, $\delta(n)=\delta_{p}(n)$.

The other basic example is as follows. Take the ring of polynomials in $n$
variables with integer coefficients, $R=\mathbb{Z}[x_{1},\dots,x_{n}]$. For
any polynomial $f(x_{1},\dots,x_{n})$, we can consider its $p$th power
$f(x_{1},\dots,x_{n})^{p}$, or we can plug in $p$th powers of the variables as
inputs to get $f(x_{1}^{p},\dots,x_{n}^{p})$. These are different, but they
agree modulo $p$ as a consequence of the “Freshman’s Dream”. Namely, in the
quotient ring $R/pR\cong\mathbb{Z}/p\mathbb{Z}[x_{1},\dots,x_{n}]$,

$(f+g)^{p}=f^{p}+\sum_{i=1}^{p-1}\binom{p}{i}f^{i}g^{p-i}+g^{p}=f^{p}+g^{p},$

since each $\binom{p}{i}$ is a multiple of $p$, and

$(fg)^{p}=f^{p}g^{p}$

as a consequence of commutativity, so the map ${f\mapsto f^{p}}$ is a ring
homomorphism in $R/pR$, called the _Frobenius map_. Thus, in $R/pR$, taking
$p$th powers before doing polynomial operations is just as good as after. So,
back in $R$ we _can_ divide the difference by $p$, and we _will_! Namely, we
can define the function

$\delta(f(x_{1},\dots,x_{n}))=\frac{f(x_{1}^{p},\dots,x_{n}^{p})-f(x_{1},\dots,x_{n})^{p}}{p},$

and this function is a $p$-derivation. Just so we can refer to this function
later, let’s call this the _standard_ $p$-derivation on
$\mathbb{Z}[x_{1},\dots,x_{n}]$ and denote it by $\delta_{\mathrm{st},p}$
(though this notation is not at all standard).

For example,

$\displaystyle\delta_{\mathrm{st},2}(x^{3}+5x)$
$\displaystyle=\frac{(x^{2})^{3}+5(x^{2})-(x^{3}+5x)^{2}}{2}$
$\displaystyle=-5x^{4}-10x^{2}.$

As this operator $\delta_{\mathrm{st},p}$ measures the failure of the
Freshman’s Dream, one might think of this as a Freshman’s Nightmare. In fact,
in large generality, $p$-derivations all arise from some freshman’s nightmare.
Let’s make this precise. Given a ring $R$, we say that a map $\Phi:R\to R$ is
a _lift of Frobenius_ if it is a ring homomorphism and the induced map from
$R/pR\to R/pR$ is just the Frobenius map, i.e., $\Phi(r)\equiv r^{p}\mod pR$
for all $R$. Given a $p$-derivation $\delta:R\to R$, the map $\Phi:R\to R$
given by

$\Phi(r)=r^{p}+p\delta(r)$

is a lift of Frobenius. Indeed, the congruence condition is automatic, and the
sum rule and product rule on $\delta$ translate exactly to the conditions that
$\Phi$ respects addition and multiplication. Conversely, if $p$ is a
nonzerodivisor on $R$, and $\Phi$ is a lift of Frobenius, then the map
$\delta(r)=\frac{\Phi(r)-r^{p}}{p}$ is a $p$-derivation: the freshman’s
nightmare associated to the lift of Frobenius $\Phi$.

It is worth noting that not every ring admits a $p$-derivation. For a quick
example, no ring $R$ of characteristic $p$ admits a $p$-derivation, since we
would have

$0=\delta(0)=\delta(p)=\delta_{p}(p)=1-p^{p-1}=1$

in $R$. Much more subtle obstructions exist, and it is an interesting question
to determine which rings admit $p$-derivations; see [MR4269423] for some
recent work on related questions.

Note that the power rule from before follows for any $p$-derivation on any
ring, since we just used the product rule to see it. The order-decreasing
property holds in general, too, at least if $p$ is a nonzerodivisor on
$R$—this follows from writing $s=p^{a}r$ with $p\nmid r$ and applying the
product rule:

$\displaystyle\delta$
$\displaystyle(p^{a}r)=p^{ap}\delta(r)+r^{p}\delta(p^{a})+p\delta(p^{a})\delta(r)$
$\displaystyle=p^{ap}\delta(r)+(p^{a-1}-p^{ap-1})r^{p}+p(p^{a-1}-p^{ap-1})\delta(r)$
$\displaystyle=p^{a-1}r^{p}+p^{a}(\delta(r)-p^{a(p-1)-1}r^{p}-p^{a(p-1)}\delta(r)).$

Let’s wrap up our cliffhanger from the previous section. Now that we have
$p$-derivations of polynomials, we have the ingredients needed for a chain
rule: given a polynomial $f(x)$ and a number $n$, we will think of the number
$f(n)$ as the composition of the function $f$ and the number $n$, and we can
try to compare $\delta_{p}(f(n))$ with $\delta_{\mathrm{st}}(f)$ and
$\delta_{p}(n)$. Here’s the chain rule:

$\delta_{p}(f(n))=\delta_{\mathrm{st},p}(f)(n)+\sum_{j=1}^{\deg(f)}{p^{j-1}}\frac{d^{j}f}{j!\,dx^{j}}(n^{p})\delta_{p}(n)^{j}.$

This is a bit more complicated than the original, but let’s notice in passing
that the $j=1$ term in the sum, $\frac{df}{dx}(n^{p})\delta_{p}(n)$, looks
pretty close to the classic chain rule, besides the $p$th power on $n$. The
curious reader is encouraged111For a hint, consider the lift of Frobenius on
$\mathbb{Z}[x]$ that sends $x\mapsto x^{p}+p\delta_{p}(n)$, and use Taylor
expansion to rewrite the associated $p$-derivation in terms of the standard
$p$-derivation and derivatives of $f$. to prove the formula above.

We have collected a decent set of analogues for the basics of differential
calculus for $p$-derivations. One can ask how far this story goes, and the
short answer is very far. Buium has developed extensive theories of arithmetic
differential equations and arithmetic differential geometry, building
analogues of the classical (nonarithmetic) versions of these respective
theories with $p$-derivations playing the role of usual derivatives. The
reader is encouraged to check out [MR2166202, MR3643159] to learn more about
these beautiful theories, though our story now diverges from these. Instead,
we will turn our attention towards using $p$-derivations to give some
algebraic results with geometric flavors.

## A Jacobian criterion

One natural geometric consideration is whether, and where, a shape has
_singularities_ : points that locally fail to look flat, due to some sort of
crossing or crinkle (or some harder to imagine higher-dimensional analogue or
a crossing or crinkle). For example, the double cone cut out by the equation
$z^{2}-x^{2}-y^{2}=0$ has a singularity at the origin where the two cones
meet, but any other point on the cone is not a singularity, see Figure 1.

Figure 1:  The cone of solutions of $z^{2}-x^{2}-y^{2}=0$. The origin is the
unique singular point.

We are going to consider shapes like this that are cut out by polynomial
equations, though to state the classical Jacobian criterion, we will consider
their solution sets over the complex numbers.

Since it is difficult to envision higher dimensional shapes (and impossible to
envision what we’re doing next!), it will be useful to give a somewhat more
algebraic heuristic definition of singularity. We will say that a point $x$ is
a _nonsingular_ point in $X$ if within $X$ one can locally cut out $x$ by
exactly $d=\dim(X)$-many equations without taking roots, and _singular_
otherwise. For example, the point $(1,0,1)$ in the cone is nonsingular, and I
claim that the two equations $y=0$, $z-1=0$ “work” for our definition: with
these two equations and the equation for $X$, we get

$y=z-1=z^{2}-x^{2}-y^{2}=0.$

Substituting in, we get $0=x^{2}-1=(x-1)(x+1)$, and “near $(1,0,1)$”, $x+1$ is
nonzero, so we can divide out and get $x-1=0$, so $x=1,y=0,z=1$. On the other
hand, $(0,0,0)$ is singular, and the two equations $y=0$, $z=0$ don’t “work”
for our definition: we have

$y=z=z^{2}-x^{2}-y^{2}=0,$

so $x^{2}=0$, but we need to take a root to get $x=0$.

The classical Jacobian criterion gives a recipe for the locus of all
singularities of a shape cut out by complex polynomials.

######  Theorem 2 (Jacobian criterion)

Let $X\subseteq\mathbb{C}^{n}$ be the solution set of the system of polynomial
equations

$f_{1}=\cdots=f_{m}=0.$

If the dimension of $X$ is $d=n-h$, and $f_{1},\dots,f_{m}$ generate a prime
ideal222We recall that an ideal $I$ is prime if it is proper and $gh\in I$
implies $g\in I$ or $h\in I$. Experts will recognize this condition as
overkill, but something needs to be done to avoid examples like $(x-y)^{2}=0$
whose solution set $X$ is the (complex) line $\\{(a,a)\\}$, which is
nonsingular, but $\frac{\partial f}{\partial x}=\frac{\partial f}{\partial
y}=0$ for every point in $X$, or $x^{2}-x=xyz=0$ whose solution set $Y$ is
two-dimensional but has a crossing singularity at $(0,0,1)$ that is not a
solution of the $1\times 1$ minors of the Jacobian matrix. in the polynomial
ring $\mathbb{C}[x_{1},\dots,x_{n}]$, then the set of singular points is the
solution set within $X$ of the system of polynomial equations

$\text{all $h\times h$ minors of}\ \begin{bmatrix}\frac{\partial
f_{1}}{\partial x_{1}}&\cdots&\frac{\partial f_{1}}{\partial x_{n}}\\\
\vdots&\ddots&\vdots\\\ \frac{\partial f_{m}}{\partial
x_{1}}&\cdots&\frac{\partial f_{m}}{\partial x_{n}}\end{bmatrix}=0.$

In particular, if $f$ is irreducible, the set of singular points of the
solution set $f=0$ is the solution set of

$f=\frac{\partial f}{\partial x_{1}}=\cdots=\frac{\partial f}{\partial
x_{n}}=0.$

For example, for the Whitney umbrella cut out by the polynomial
$f=x^{2}-y^{2}z$, the singular locus is cut out by the system

$x^{2}-y^{2}z=2x=2yz=y^{2}=0,$

which simplifies to $x=y=0$; the $z$-axis is where the shape crosses itself.

Figure 2:  The Whitney umbrella of solutions of ${x^{2}-y^{2}z=0}$. The
singular locus consists of the line where it crosses itself.

The notion of (non)singularity in geometry is generalized in algebra by the
notion of regular ring. For starters, prime ideals in algebra play the role of
points in geometry: this is motivated by Hilbert’s Nullstellensatz, which says
that for a quotient ring of the form333For a set of elements
$a_{1},\dots,a_{t}$ in a ring $R$, we use the notation
$(a_{1},\dots,a_{t})\colonequals\\{r_{1}a_{1}+\cdots+r_{t}a_{t}\ |\ r_{i}\in
R\\}$ for the _ideal generated by $\\{a_{1},\dots,a_{t}\\}$_.

$\frac{\mathbb{C}[x_{1},\dots,x_{n}]}{(f_{1},\dots,f_{m})},$

every maximal ideal is of the form

$\mathfrak{m}_{a}=(x_{1}-a_{1},\dots,x_{n}-a_{n})$

for some $a=(a_{1},\dots,a_{n})$ solution to

$f_{1}(a)=\cdots=f_{m}(a)=0;$

including all prime ideals leads to a better theory for general rings. Then we
say that a prime ideal in a ring $R$ is _nonsingular_ or _regular_ at a prime
ideal $\mathfrak{q}$ if $\mathfrak{q}$ can be generated “locally” by $h$
equations, where $h$ is the codimension of $\mathfrak{q}$ (how much
$\mathfrak{q}$ cuts down the dimension of $R$), and “locally”
means444Precisely, “locally” means that we work in the localization
$R_{\mathfrak{q}}$, and codimension refers to the height of the ideal
$\mathfrak{q}$. that one can divide by elements outside of $\mathfrak{q}$. In
the motivating geometric situation where $X$ is the solution set of
$f_{1}=\cdots=f_{m}=0$ over $\mathbb{C}$, the point $a\in X$ is nonsingular if
and only if $R=\mathbb{C}[x_{1},\dots,x_{n}]/(f_{1},\dots,f_{m})$ is
nonsingular at the maximal ideal $\mathfrak{m}_{a}$. A prime ideal is
_singular_ if it is not nonsingular.

Intuitively, when working over $\mathbb{Z}$ rather than over $\mathbb{C}$ or a
field, in addition to the geometric dimensions, there is an arithmetic
dimension that corresponds to the prime integers $p$ in $\mathbb{Z}$. To
detect singularity, it suffices to include $p$-derivations as a substitute for
derivatives in the $p$-direction! The following is a special case of a result
independently obtained by Saito [MR4412577] and Hochster and the author [HJ].

######  Theorem 3 (Saito, Hochster-Jeffries)

Let
${R=\displaystyle\frac{\mathbb{Z}[x_{1},\dots,x_{n}]}{(f_{1},\dots,f_{m})}}$
for some prime ideal $(f_{1},\dots,f_{m})$ of $\mathbb{Z}[x_{1},\dots,x_{n}]$.
Then the set of singular prime ideals $\mathfrak{q}$ of $R$ containing a fixed
prime integer $p$ is exactly the set of prime ideals containing $p$ and

$\begin{subarray}{c}\text{all $h\times h$}\\\ \text{ minors of}\end{subarray}\
\begin{bmatrix}\delta_{\mathrm{st},p}(f_{1})&(\frac{\partial f_{1}}{\partial
x_{1}})^{p}&\cdots&(\frac{\partial f_{1}}{\partial x_{n}})^{p}\\\
\vdots&&\ddots&\vdots\\\ \delta_{\mathrm{st},p}(f_{m})&(\frac{\partial
f_{m}}{\partial x_{1}})^{p}&\cdots&(\frac{\partial f_{m}}{\partial
x_{n}})^{p}\end{bmatrix},$

where $h=n+1-\dim(R)$. In particular, if $f$ is irreducible, the set of
singular prime ideals of $R=\frac{\mathbb{Z}[x_{1},\dots,x_{n}]}{(f)}$
containing $p$ is exactly the set of primes containing $p$ and

$\delta_{\mathrm{st},p}(f)\ ,\ \frac{\partial f}{\partial x_{1}}\ ,\ \ldots\
,\ \frac{\partial f}{\partial x_{n}}.$

###### Example 1

Let $n$ be a squarefree integer (excluding $0$ and $1$), $q$ a prime number,
and consider the ring $R=\mathbb{Z}[\sqrt[q]{n}]$. We claim that this admits a
singular prime ideal if and only if $\delta_{q}(-n)$ is a multiple of $q$.
Think555This is where excluding $0$ and $1$ is necessary. of $R$ as
$\mathbb{Z}[x]/(x^{q}-n)$. For a prime number $p$, the singular prime ideals
$\mathfrak{a}$ containing $p$ are those that contain

$p,\delta_{\mathrm{st},p}(x^{q}-n),\frac{d}{dx}(x^{q}-n).$

Using the sum rule for $\delta_{\mathrm{st},p}$ and the defining equation for
$R$, we have

$\delta_{\mathrm{st},p}(x^{q}-n)=\begin{cases}\delta_{p}(-n)&\text{ if $p\neq
2$}\\\ \delta_{p}(-n)-nx^{q}&\text{ if $p=2$}\end{cases}$

in $R$.

For $p\neq q$, from $\frac{d}{dx}(x^{q}-n)=qx^{q-1}$, we get that
$\mathfrak{a}$ must contain $x$, and from the defining equation, $n$ as well.
But if the integers $p,n,\delta_{p}(n)$ are in a proper ideal $\mathfrak{a}$,
since $p$ is a prime number, we must have $p\mid n$ and $p\mid\delta_{p}(n)$,
since $1$ would be a linear combination of these numbers otherwise. By the
Order Decreasing Property, $p^{2}\mid n$, contradicting that $n$ is
squarefree. So, there are no singular prime ideals containing ${p\neq q}$.

For $p=q\neq 2$ since $\frac{d}{dx}(x^{q}-n)$ is a multiple of $q$, using the
simplification above, a prime ideal containing $q$ is singular if and only if
it contains

$q,\delta_{q}(-n).$

But if $q$ and $\delta_{q}(-n)$ are in $\mathfrak{a}$, then
$q\mid\delta_{q}(-n)$; a singular prime ideal then occurs if and only if this
happens. The analysis for $p=q=2$ is similar (cf. [HJ]).

This has a consequence for a familiar object in elementary number theory. By
some standard results666Namely, $\mathbb{Z}[\sqrt[q]{n}]$ is a ring of
integers if and only if it is integrally closed in its fraction field. Since
this ring is an integral extension of $\mathbb{Z}$ generated by one element as
an algebra, it is a one dimensional Noetherian domain. Such a ring is
integrally closed in its fraction field if and only if it is nonsingular. in
commutative algebra, the ring $\mathbb{Z}[\sqrt[q]{n}]$ is a ring of algebraic
integers (for its fraction field $\mathbb{Q}(\sqrt[q]{n})$) if and only if it
has no singular prime ideals. Thus, we conclude that $\mathbb{Z}[\sqrt[q]{n}]$
is a ring of integers if and only if $q\nmid\delta_{q}(-n)$. In particular,
for $q=2$, using our earlier observation on when $\delta_{2}(n)$ is odd or
even, we recover the fact that for $n$ squarefree, $\mathbb{Z}[\sqrt{n}]$ is a
ring of integers if and only if $n\not\equiv 3\mod 4$.

## A Zariski-Nagata theorem for symbolic powers

Let’s recall another classical theorem relating algebra and geometry. To state
it, we need the notion of symbolic power of a prime ideal. Over a polynomial
ring over a field, or more generally, over a commutative Noetherian ring, any
ideal can be written as an intersection

$I=\mathfrak{q}_{1}\cap\cdots\cap\mathfrak{q}_{t}$ (3)

where the $\mathfrak{q}_{i}$ are _primary_ ideals: ideals with the property
that $xy\in\mathfrak{q}\Rightarrow x\in\mathfrak{q}\text{ or
}y^{n}\in\mathfrak{q}$ for some $n$. Such an expression is called a _primary
decomposition_. The existence of primary decomposition is a famous result of
Lasker in the polynomial case and Noether in the Noetherian case. This can be
thought of as a generalization of the Fundamental Theorem of Arithmetic: in
$\mathbb{Z}$, the primary ideals are just the ideals generated by powers of
primes, and a prime factorization

$n=p_{1}^{e_{1}}\cdots p_{n}^{e_{n}}$

corresponds to writing $(n)$ as an intersection of primary ideals:

$(n)=(p_{1}^{e_{1}})\cap\cdots\cap(p_{n}^{e_{n}}).$

There are two important differences with the Fundamental Theorem of
Arithmetic, though:

  * •

primary ideals are not powers of prime ideals in general, nor are powers of
prime ideals always primary, and

  * •

the collection of primary ideals $\mathfrak{q}_{i}$ appearing in the
decomposition (3) is not unique, but if the decomposition satisfies a simple
irredundancy hypothesis, then the components whose radical does not contain
any other component’s radical are uniquely determined.

In particular, if $\mathfrak{q}$ is a prime ideal and $n>1$, then according to
the first point above, $\mathfrak{q}^{n}$ may admit an interesting primary
decomposition, and as a consequence of the second point, the component with
radical $\mathfrak{q}$ is the same in any (irredundant) primary decomposition.
This is called the $n$th _symbolic power_ of $\mathfrak{q}$, denoted
$\mathfrak{q}^{(n)}$.

Symbolic powers arise in various contexts in algebra and geometry. For
example, they arise naturally in interpolation questions, they play a key role
in the proofs of various classical theorems such as Krull’s Principal Ideal
Theorem, and they have enjoyed a resurgence of interest in combinatorics in
connection with the Packing Problem of Conforti and Cornuéjols. The interest
reader is recommended to read the survey [MR3779569].

A classical pair of theorems of Zariski and Nagata gives a geometric
description of the symbolic power of an ideal in a polynomial ring over
$\mathbb{C}$. The result has various statements; we will give a differential
statement.

######  Theorem 4 (Zariski-Nagata Theorem)

Let ${X\subseteq\mathbb{C}^{n}}$ be the solution set of the system of
polynomial equations

$f_{1}=\cdots=f_{m}=0.$

Suppose that $f_{1},\dots,f_{m}$ generate a prime ideal $\mathfrak{q}$. Then
$\mathfrak{q}^{(r)}$ is exactly the set of polynomials
${f\in\mathbb{C}[x_{1},\dots,x_{n}]}$ such that

$\frac{\partial^{a_{1}+\cdots+a_{n}}f}{\partial x_{1}^{a_{1}}\cdots\partial
x_{n}^{a_{n}}}\Big{|}_{X}\equiv 0\text{ for all }a_{1}+\cdots+a_{n}<r.$

The same characterization is doomed to fail in $\mathbb{Z}[x]$: for example,
$\mathfrak{q}=(2)$ is a prime ideal with $\mathfrak{q}^{(2)}=(4)$; in
particular, $2\notin(4)$. But $2$ satisfies the derivative condition
corresponding to the right hand side above: taking $a_{1}=0$, we have
$2\in(2)$, and taking $a_{1}=1$, we have $\frac{\partial 2}{\partial
x}=0\in(2)$.

If you’ve been paying attention so far, you should be able to name the missing
ingredient. Indeed, $\delta_{\mathrm{st},2}(2)=-1\notin(2)$, so allowing
partial derivatives and a $2$-derivation is enough to take this element $2$
that isn’t in $\mathfrak{q}^{(2)}$ out of $\mathfrak{q}$, suggesting a way to
characterize $f\in\mathfrak{q}^{(2)}$ in terms of derivatives (including our
“derivative by 2”).

In fact, this works in general. The following analogue of the Zariski-Nagata
Theorem is a special case of a result of De Stefani, Grifo, and the author
[MR4080246].

######  Theorem 5 (De Stefani-Grifo-Jeffries)

Consider the ring $\mathbb{Z}[x_{1},\dots,x_{n}]$ and let
$\mathfrak{q}=(f_{1},\dots,f_{m})$ be a prime ideal. Suppose that
$\mathfrak{q}$ contains the prime integer $p$. Then $\mathfrak{q}^{(r)}$ is
exactly the set of polynomials $f\in\mathbb{Z}[x_{1},\dots,x_{n}]$ such that

$\delta_{\mathrm{st},p}^{a_{0}}\left(\frac{\partial^{a_{1}+\cdots+a_{n}}f}{\partial
x_{1}^{a_{1}}\cdots\partial x_{n}^{a_{n}}}\right)\in\mathfrak{q}\text{ for all
}a_{0}+a_{1}+\cdots+a_{n}<r.$

## Other applications

$p$-derivations have appeared in a range of sophisticated applications to
number theory and arithmetic geometry. We briefly list a few of these, and
encourage the reader to explore further.

### Effective bounds on rational points and $p$-jet spaces

The motivation for Buium’s original work on $p$-derivations was to give bound
the number of points on rational points on curves. For a complex curve $X$
defined over $\overline{\mathbb{Q}}$, there is a natural map from $X$ to an
algebraic group $A$, the Jacobian of $X$; the main result of Buium [MR1387233]
gives an effective bound, depending only on the genus of $X$ and the smallest
prime of a good reduction for $X$, on the number of points in $X$ that map to
torsion elements of $A$.

To establish these results, Buium constructs _$p$ -jet spaces_. In
differential geometry, the _jet_ of a function $f$ of order $k$ at a point is
the data of all of the values of the derivatives of $f$ up to order $k$ at
that point; the _jet space_ of order $k$ of a manifold $X$ is a manifold whose
points correspond to jets on $X$. Buium’s $p$-jet spaces are analogues of jet
spaces obtained by replacing usual derivatives with $p$-derivations. The
result mentioned above is then obtained by intersection theory on $p$-jet
spaces of curves.

As mentioned earlier, Buium has developed an extensive theory of arithmetic
differential geometry in analogy with classical differential geometry, for
which $p$-jet spaces form the starting point. We refer the reader to
[MR2166202] to learn more.

### Relationship with Witt vectors

The Witt vectors are a construction in number theory that generalizes the
relationship between the prime field $\mathbb{F}_{p}$ and the corresponding
ring of $p$-adic integers $\mathbb{Z}_{p}$. To be precise, there is a functor
$W$ from the category of rings to itself, called the functor of ($p$-typical)
Witt vectors, that maps $\mathbb{F}_{p}$ to $\mathbb{Z}_{p}$. More generally,
for any perfect field $F$ of characteristic $p$, $W(F)$ is a local ring with
maximal ideal $(p)$, and $W(F)/pW(F)\cong F$; in this way, one can think of
$W$ as a generalization of a construction of $\mathbb{Z}_{p}$ from
$\mathbb{F}_{p}$.

$p$-derivations have many interesting connections with Witt vectors; indeed,
they first arise in Joyal’s work to study Witt vectors. Namely, Joyal shows
that the forgetful functor from the category of rings with a $p$-derivation to
the category of rings is left adjoint to the Witt vector functor.

### Philosophy of the field with one element

Various formulas for enumerating basic objects over a finite field
$\mathbb{F}_{q}$ limit to combinatorially meaningful quantities as $q\to 1$;
for example,

$\lim_{q\to
1}\frac{|\mathrm{GL}_{n}(\mathbb{F}_{q})|}{|(\mathbb{F}_{q}^{\times})^{n}|}=n!,$

where $\mathrm{GL}_{n}(\mathbb{F}_{q})$ is the collection of linear
automorphisms of $\mathbb{F}_{q}^{n}$, the vector space of $n$-tuples over
$\mathbb{F}_{q}$, and $n!$ counts the number of permutations of $n$ elements.
Partially motivated by this phenomenon, and partially motivated by
transferring results over finite fields to other settings, there is a program
of inventing a notion of algebraic geometry over a “field with one element”:
the mythical field with one element is not literally a field with one element,
which would contradict the definition of field, but something with a different
structure than a field that admits some sort of geometry analogous to what one
would expect over a field with one element for quantitative or various other
reasons. On the algebraic side, the field with one element can be thought of
as a “deeper base ring” than $\mathbb{Z}$ (though not necessarily itself a
ring!).

There are many approaches in the literature to implementing the philosophy of
the field with one element. Most relevant for this article is the theory
established by Borger [Bor] as well as the closely related approach of Buium.
Roughly speaking, Borger proposes that a model for the field with one element
should be $\mathbb{Z}$ equipped with the collection of all its $p$-derivations
$\delta_{p}$. In particular, in Borger’s model, $\mathrm{GL}_{n}$ “over the
field with one element” is the symmetric group on $n$ letters, aligning with
the numerical coincidence noted in the previous paragraph.

### Unifying cohomology theories

Recent work of Bhatt and Scholze [MR4502597] has employed $p$-derivations to
relate various $p$-adic cohomology theories (étale, de Rham, crystalline). The
key ingredient is the notion of a _prism_ , which is a ring equipped with a
$p$-derivation subject to some conditions. We refer the reader to [MR4502597]
to learn more.

## Acknowledgments

The author is grateful to Alessandro De Stefani, Eloísa Grifo, Claudia Miller,
Steven Sam, and the referee for many helpful comments on this article.

## References


We study vehicle dispatching in autonomous mobility on demand (AMoD) systems, where a central operator assigns vehicles to customer requests or rejects these with the aim of maximizing its total profit. Recent approaches use multi-agent deep reinforcement learning (MADRL) to realize scalable yet performant algorithms, but train agents based on local rewards, which distorts the reward signal with respect to the system-wide profit, leading to lower performance. We therefore propose a novel global-rewards-based MADRL algorithm for vehicle dispatching in AMoD systems, which resolves so far existing goal conflicts between the trained agents and the operator by assigning rewards to agents leveraging a counterfactual baseline. Our algorithm shows statistically significant improvements across various settings on real-world data compared to state-of-the-art MADRL algorithms with local rewards. We further provide a structural analysis which shows that the utilization of global rewards can improve implicit vehicle balancing and demand forecasting abilities. Our code is available at <https://github.com/tumBAIS/GR-MADRL-AMoD>.

multi-agent learning, credit assignment, deep reinforcement learning, autonomous mobility on demand

§ INTRODUCTION

Within the coming years, in autonomous mobility on demand (AMoD) promises to significantly improve urban passenger transportation. An AMoD system enables its operator to exercise full control over the vehicle fleet, thus improving vehicle dispatching performance, provided that the operator uses an effective control algorithm. This ultimately yields increased vehicle utilization while preserving the comfort of individual mobility. By deciding which requests to accept and assigning vehicles to these requests, the operator faces a contextual multi-stage stochastic control problem. Various approaches to solve this problem exist, ranging from model predictive control <cit.> to deep reinforcement learning (DRL) <cit.>. The latter shows better performance due to its model-free nature. In this context, multi-agent approaches allow to achieve scalability to large instances. To train the agents' policy, existing work uses local, egoistic, per-agent rewards, which can lead to sub-optimal behavior from the perspective of the central operator interested in maximizing the system-wide profit.

In this work, we improve upon the state-of-the-art local-rewards algorithm (LRA) of <cit.> by proposing a novel way to train agents with global rewards. The core of our methodology is a new advantage function, which estimates the individual agents' contribution to the global reward, as such aligning their goal with the central operator, ultimately following a system-optimal policy. Our algorithm outperforms LRA by up to 2% on average across test dates and up to 6% on individual dates, which is a substantial performance improvement in AMoD. Additionally, our algorithm is as scalable as LRA, accordingly setting a new state-of-the-art for vehicle dispatching in AMoD systems.

§.§ Related work

Our work relates to two research streams: from an application perspective to AMoD fleet control, and methodologically to multi-agent deep reinforcement learning (MADRL). We review related literature from both streams in the following.

Algorithms for (A)MoD fleet control range from greedy heuristics <cit.>, to (stochastic) model predictive control <cit.>, combinations of optimization and supervised learning <cit.>, and DRL. Existing DRL-based algorithms differ from our approach, as they often use value-based learning <cit.>, which can decrease learning efficiency and speed. Others focus on explicit rebalancing instead of request assignment <cit.>. Furthermore, most previous work considers non-autonomous MoD, aiming to maximize revenues <cit.>, while AMoD should focus on profit maximization <cit.>. Finally, the  work of <cit.> uses policy-based hybrid MADRL for a profit-maximizing AMoD operator, but trains agents on local rewards, possibly leading to sub-optimal agent behavior. To the best of our knowledge, no work exists that addresses the problem of profit-maximizing vehicle dispatching in AMoD without the aforementioned shortcomings. We address this research gap by extending the work of <cit.> to incorporate global rewards, thus creating goal congruence between the agents and the system operator.

A crucial challenge in our setting is deriving per-agent contributions to global success from a shared reward signal, i.e., a credit assignment problem <cit.>. Solution approaches like inverse reinforcement learning <cit.> or value decomposition <cit.> are not applicable to our setting, because we cannot observe the behavior of an optimal agent and do not use Q-learning. An alternative approach is reward marginalization, based on difference rewards <cit.>, which uses functions (e.g., advantage functions) to estimate the contribution of individual agents to global rewards <cit.>. Reward marginalization approaches often use actor-critic algorithms. However, reward marginalization approaches often suffer from high learning variance and poor sample efficiency, as they are typically built on top of basic DRL algorithms. To the best of our knowledge, no work exists that combines reward marginalization approaches with low-variance actor-critic algorithms so far. We address this research gap by embedding reward marginalization into a Soft Actor-Critic (SAC) framework for discrete actions.

§.§ Contribution

To close the research gap outlined above, we develop the first scalable MADRL-based algorithm for vehicle dispatching in profit-maximizing AMoD that trains agents via global rewards. In this context, we combine SAC for discrete actions with credit assignment based on a counterfactual baseline to resolve goal conflicts between the trained agents and the operator's global objective. Our algorithm combines the benefits of low learning variance and sample efficiency of SAC with the benefits of credit assignment via a counterfactual baseline. Since this algorithm based on purely global rewards scales only to medium-sized problem instances, we additionally develop a scheduled algorithm that combines local rewards and global rewards with our counterfactual baseline. Thus, we obtain a powerful new MADRL algorithm with possible applications beyond AMoD. We evaluate our algorithm on real-world taxi data from New York City and show that it outperforms LRA of <cit.> by up to 2% on average across test dates and up to 6% on individual dates. This constitutes a substantial performance improvement for the AMoD application case, where a single percent improvement yields significant daily monetary gains. We further provide a structural analysis which shows that the utilization of global rewards can improve implicit vehicle balancing and demand forecasting abilities.

§ PROBLEM SETTING

Exemplary vehicle dispatching process.

We focus on the contextual multi-stage stochastic control problem from <cit.>, illustrated in Figure <ref>, to ensure comparability. In this problem setting, a central operator controls an AMoD fleet and dispatches vehicles to serve requests. We consider a discrete time horizon. During each time step, customers submit a variable number of new requests for point-to-point transportation, of which the operator has no prior knowledge. At the beginning of each time step, the operator makes a decision for the batch of requests that were placed during the previous time step. For each request, this decision is either to reject the request or to assign it to a vehicle. The operator can base its decision on fully observable state information, including the requests' origins and destinations as well as the vehicles' positions and already assigned but not yet completed requests. Each request has to be decided upon immediately and customers have a fixed maximum waiting time until they must be picked up. Based on the operator's decision, the system transitions to the next state: vehicles move towards their next destination and potentially pick up or drop off customers. Customers place new requests, while rejected requests leave the system. We simulate the requests by replaying historical data. Finally, the operator receives the system-wide profit as a reward. The profit is calculated as the sum of all revenues from accepted requests minus the costs for all vehicle movements. Revenues and costs are linear functions of the distance travelled. Since fixed costs do not depend on the dispatching problem, we only include operational costs (e.g., for fuel and maintenance). For further details, including a formal Markov decision process, see <cit.>.

§ METHODOLOGY

In the following, we first motivate our baseline algorithm <cit.> and extend it with a naive global reward allocation scheme (Section <ref>). We then focus on reward marginalization and show how to modify the existing Counterfactual Multi-Agent Policy Gradient (COMA) paradigm to SAC architectures (Sections <ref> & <ref>). We finally discuss enhancements to scale our algorithm to large-scale instances (Section <ref>).

Outline of base algorithm. Black parts are used during training and testing, gray parts only during training.

§.§ Base algorithm

Our problem setting has two peculiarities influencing the construction of an effective vehicle dispatching algorithm: the global action space is too large to use single-agent DRL and the number of requests varies between time steps. When addressing the first peculiarity with multi-agent learning, the need for agent coordination arises. The algorithm of <cit.> therefore uses a two-stage architecture, see Figure <ref>. The first stage consists of DRL-trained agents, the second one of an optimization-based central matching. In the first stage, each agent is the combination of a vehicle and a request. An agent can either accept the assignment of the respective request to the respective vehicle or reject this assignment, resulting in a per-agent action space of size two. As the number of requests varies, so does the number of agents. Each agent observes the global state, using an encoding and an attention mechanism to cope with the varying number of requests. Using an actor neural network with a softmax function, the agent obtains probabilities for acceptance and rejection. The per-agent post-processing selects the action by sampling during training or taking the action with the highest probability during testing. It then sets the agent's score to the acceptance probability in case of an acceptance and to zero in case of a rejection. The scores of all agents are submitted to the global matching, which coordinates the agents. This bipartite matching selects the global action by maximizing the submitted scores under the constraints that every request can be assigned at most once and that every vehicle can take at most one new request.

We train the agents using SAC for discrete actions <cit.>. As the use of multiple agents only serves to handle the large action space, we train a single actor and a single critic network and use them for all agents. The rewards we use for training are the profits obtained by the algorithm: when a request is assigned to a vehicle and served within the maximum waiting time, the respective agent immediately receives the profit as a reward.

To transition from this algorithm to  our contribution, we first follow a naive approach to include global rewards in training by replacing the per-agent rewards with global rewards when training the critic, resulting in a basic global rewards algorithm (GRA). We obtain global rewards by summing the profits of all agents at one time step. As this number can be substantially larger than per-agent rewards, we divide it by the average number of non-zero rewards per observation in the replay buffer to stabilize learning. Straightforwardly using this approach leads to a credit assignment problem, as the reward given to agents now depends on the actions of all agents, which generally complicates an agent's learning task. To mitigate this, we explore a credit assignment paradigm in the following.

§.§ Naive COMA

A suitable credit assignment paradigm for our setting is COMA, proposed by <cit.>. COMA fits our setting best, because it has a similar structure as SAC and is especially suitable for small per-agent action spaces. Following the main rationale of COMA, agents should maximize their contribution to the global reward instead of maximizing the global reward directly. As obtaining global rewards for several actions is computationally infeasible, COMA trains the critic on global rewards to approximate global state-action-values. The contribution of an agent to this global value is defined as the value of taking an action in contrast to the value of taking a default action. This default action is calculated as the policy-weighted average value of all possible actions the agent can take (counterfactual baseline). The advantage function of COMA therefore is

\begin{equation}
A_i(a_i|s,i)=Q_\theta(a_i|s,\bar a_{-i})-\sum_{a'_i}\pi_\phi(a'_i|s,i) \:Q_\theta(a'_i|s,\bar a_{-i}).
\end{equation}

In this and all following equations, $s$ denotes the global state, $a'_i$ an action agent $i$ can take as a reject/accept decision before the global matching, $a_i$ the action the agent actually takes, and $\bar a_{-i}$ the actions of all agents except agent $i$ after the global matching. Let $\pi_\phi(a_i|s,i) \in [0,1]$ be the probability that agent $i$ takes action $a_i$ in state $s$, following policy $\pi_\phi$ parameterized by network $\phi$. $Q_\theta(a_i|s,\bar a_{-i}) \in \mathbb{R}^2$ is the global state-action-value (Q-value) of action $a_i$ taken by agent $i$ in state $s$ estimated using network $\theta$, given other agents' actions $\bar a_{-i}$.

<cit.> use a sampling-based approach to estimate the actor's loss function and base its computation only on the action taken by the agent. The loss function thus reads $J_\pi(\phi)=\mathbb{E}_{s\sim D}\bigl[ \sum_i A_i(a_i|s,i)\bigr]$, with $D$ denoting the replay buffer. In contrast to that, <cit.> use SAC with discrete actions, which considers all possible actions of an agent in the loss function

\begin{equation}\label{eqn:ego_loss}
J_\pi(\phi)=\mathbb{E}_{s\sim D}\left[ \sum_i \sum_{a_i} \pi_\phi(a_i|s,i) \left(\alpha\log\pi_\phi(a_i|s,i)-\min_{j \in \{1,2\}} \left\{ Q_\theta^j(a_i|s,\bar a_{-i}) \right\} \right) \right],
\end{equation}

with $\alpha\in\mathbb{R}$ being the entropy coefficient. To use credit assignment via COMA in combination with the proven reliability and low variance of SAC, we need to integrate the baseline of COMA into the loss function of SAC. In the following, we use $\pi(a_i):=\pi_\phi(a_i|s,i)$ and $Q(a_i):=\min_{j \in \{1,2\}} \left\{ Q_\theta^j(a_i|s,\bar a_{-i}) \right\}$ for conciseness. Then, considering one instance of a batch and one agent, the loss function of SAC extended by the advantage of COMA reads

\begin{equation}\label{eqn:loss}
J_\pi(\phi|s,i)=\sum_{a_i}\pi(a_i) \left(\alpha\log\pi(a_i)-Q(a_i)+\sum_{a'_i}\pi(a'_i)\: Q(a'_i) \right).
\end{equation}

The loss function $J_\pi(\phi|s,i)$ as defined in Equation (<ref>) is equivalent to the entropy $J_\pi(\phi|s,i)=\sum_{a_i}\pi(a_i)\: \alpha\log\pi(a_i)$ of a plain SAC architecture.

For a proof of Proposition <ref>, we refer to Appendix <ref>.
Accordingly, using the loss function as derived in Equation (<ref>) does not allow to learn a meaningful policy such that we cannot straightforwardly apply the COMA paradigm to our SAC framework. This observation motivates us to study a novel approach to combine SAC for discrete actions with the COMA paradigm to learn a good policy with global rewards.

§.§ Adjusted COMA for SAC architectures

To solve the convergence problem outlined above, we have to adjust the loss function in Equation (<ref>). Using only the action taken by the agent for the loss function as in <cit.> does not lead to convergence even in small experimental instances, as it increases the loss function's variance. We therefore adjust $\pi(a'_i)$, changing the weighting of the default action in the baseline. This is possible from a theoretical perspective, as the exact specification of a default action is not derived from the idea of difference rewards <cit.>, but left to the discretion of the user.

Straightforwardly, we can define the default action by using an equally-weighted average instead of a policy-weighted average, resulting in the advantage function $A^\text{equ}_i(a_i)=Q(a_i)-\sum_{a'_i}\frac{1}{n_{a_i}}\: Q(a'_i)$, with $n_{a_i}$ being the number of actions per agent. We call this algorithm $\text{COMA}^\text{equ}$, which resolves the convergence problem, but has a disadvantage: when the actor network estimates  different probabilities for single actions, weighting all actions equally is not a reasonable default action. This problem is especially pronounced during late training, when the actor network is better at estimating action probabilities. We solve this issue by defining a second default action with the use of a target actor network $\bar\phi$. This network has the same structure and initialization as the actor network $\phi$ and is updated using exponential averages of the actor network's parameters, similar to target networks in Q-learning. The advantage function of this algorithm ($\text{COMA}^\text{tgt}$) is $A^\text{tgt}_i(a_i)=Q(a_i)-\sum_{a'_i}\pi_{\bar\phi}(a'_i)\: Q(a'_i)$. Since $\bar\phi$ differs from $\phi$, this algorithm solves the convergence problem as well. During early training, $\text{COMA}^\text{tgt}$ is not as suitable as $\text{COMA}^\text{equ}$, since the sub-optimal action probabilities of an untrained target actor network are a disadvantage compared to equally-weighted actions. Later in training, $\bar\phi$ can estimate better action probabilities, making $\text{COMA}^\text{tgt}$ superior to $\text{COMA}^\text{equ}$.

Since we now have one algorithm especially suitable for early learning and one especially suitable for later learning, we combine these two and obtain $\text{COMA}^\text{adj}$, based on a dynamic combination of the two newly introduced advantage functions. The advantage function of $\text{COMA}^\text{adj}$ thus reads

\begin{equation} \label{eqn:COMA_adj}
\begin{split}
A^\text{adj}_i(a_i)&=(1-\beta)A^\text{equ}_i(a_i)+\beta A^\text{tgt}_i(a_i) \\
&=Q(a_i)-(1-\beta)\sum_{a'_i}\frac{1}{n_{a_i}}\: Q(a'_i)-\beta\sum_{a'_i}\pi_{\bar\phi}(a'_i)\: Q(a'_i).
\end{split}
\end{equation}

Here, the hyperparameter $\beta\in[0,1]$ is the weight of the $\text{COMA}^\text{tgt}$ baseline and follows a schedule, starting at zero and ending at one, see Appendix <ref>. Then, the loss function of $\text{COMA}^\text{adj}$ reads

\begin{equation}
J^\text{adj}_\pi(\phi|s,i)=\sum_{a_i}\pi(a_i) \left(\alpha\log\pi(a_i)-A^\text{adj}_i(a_i)\right).
\end{equation}

With this loss function, $\text{COMA}^\text{adj}$ solves the credit assignment problem. In our experiments, $\text{COMA}^\text{adj}$ performs better than LRA, but has a scalability problem: when the number of agents increases beyond medium-sized problem instances, $\text{COMA}^\text{adj}$ fails to converge. Reasons for this are the diminishing influence of a single agent on global rewards and the overlap of many agents' actions when the number of agents increases, making learning per-agent Q-values difficult <cit.>. We therefore investigate how to scale $\text{COMA}^\text{adj}$.

§.§ Reward scheduling

Usually, one could resolve the scalability problem of $\text{COMA}^\text{adj}$ straightforwardly by adjusting the critic to accommodate value factorization <cit.>, but this approach is infeasible in our setting as the number of agents is variable. Similarly, learning the critic on a static mix of local and global rewards in a local-global-rewards algorithm (LGRA) does not solve the scalability problem, since any non-negligible share of global rewards distorts learning when increasing the number of agents. In addition, reward marginalization with a counterfactual baseline is problematic for partially local rewards.

Instead, we can train a single actor network using a weighted average policy loss function, consisting of the loss function for LRA, $J^\text{loc}_\pi (\phi|s,i)$, by <cit.> and $\text{COMA}^\text{adj}$. The loss function thus reads

\begin{equation}
J^\text{scd}_\pi(\phi|s,i)=(1-\kappa)\:J^\text{loc}_\pi (\phi|s,i)+\kappa\:J^\text{adj}_\pi(\phi|s,i),
\end{equation}

with $\kappa\in[0,1]$ being the weight of the loss function of $\text{COMA}^\text{adj}$. Again, $\kappa$ follows a schedule increasing linearly, following a polynomial pattern or jumping from zero to one at a specified point, see Appendix <ref>. This leads to a new algorithm, we call it $\text{COMA}^\text{scd}$. $\text{COMA}^\text{scd}$ solves the scalability problem, as it enables the learning of global Q-values when increasing the number of agents. The reason this works is the utilization of experience collected following a mixed policy: this way, more diverse experience is available than if using solely own experience, thus improving learning without the destabilizing influence of increasing the entropy. We can therefore train four critic networks from the beginning on, two for local and two for global rewards. Two networks each are necessary for SAC, where we always use the minimum of the two state-action-values to avoid value overestimation. Due to the influence of both local and global rewards, $\text{COMA}^\text{scd}$ can sometimes have a lower performance than algorithms purely based on global rewards, but makes up for this by being as scalable as LRA, thus sacrificing a portion of its performance for scalability.

§ NUMERICAL STUDIES

We benchmark our algorithms using the experimental design of <cit.>, which bases on New York Taxi data <cit.> in a hexagonal grid of Manhattan, with 38 large, 11 small, or 5 small zones. In this realm, we study five instances: two edge cases with high (5 zones, 15 vehicles) and low (11 zones, 6 vehicles) acceptance rates, two typical test cases (11 zones, 18 and 24 vehicles), and a comparatively large instance (38 zones, 100 vehicles). For details on the experimental design, we refer to Appendix <ref> and report hyperparameters for all models in Appendix <ref>.

Firstly, we test $\text{COMA}^\text{scd}$ on all five instances and benchmark it against LRA of <cit.> and a greedy algorithm, which considers requests in their order of submission, accepting profitable ones and rejecting all others. Secondly, we present an ablation study to show the superiority of $\text{COMA}^\text{scd}$ over our alternative algorithms that use global rewards. Thirdly, we discuss why $\text{COMA}^\text{scd}$ outperforms LRA.

§.§ Performance of $\text{COMA}^\text{scd}$

Relative test performance $\Delta\ [\%]$ of $\text{COMA}^\text{scd}$ vs. greedy and LRA for multiple test dates.

5 zones, 15 veh.     11 zones, 6 veh.     11 zones, 18 veh.     11 zones, 24 veh.     38 zones, 100 veh.
vs. LRA     0.5%     -0.2%     0.2%     1.9%     0.6%

p-value     0.05     0.52     0.22     0.00     0.00

vs. greedy     2.5%     2.8%     2.8%     3.5%     2.4%

p-value     0.00     0.00     0.01     0.00     0.00

Mean test performance improvement of $\text{COMA}^\text{scd}$ vs. LRA and greedy, including the respective Wilcoxon p-values.

We present results of our tests of $\text{COMA}^\text{scd}$ in Figure <ref> and Table <ref>. In all instances except the one with 11 zones and 6 vehicles, $\text{COMA}^\text{scd}$ outperforms LRA and the greedy algorithm on average, the former by up to 1.9% and the latter by up to 3.5%. On single test dates, $\text{COMA}^\text{scd}$ can outperform LRA by up to 6%. This improvement is significant, as the Wilcoxon p-values are at most 5% for the respective instances. While these relative improvements appear to be small, they are substantial in AMoD and of a similar magnitude as the improvements of previous state-of-the-art algorithms over their respective benchmarks <cit.>.

In the instance with a high acceptance rate (5 zones, 15 vehicles), the significant performance improvement of $\text{COMA}^\text{scd}$ compared to LRA is an especially positive result, as a vehicle is usually available for each request in this instance, limiting the improvement potential for DRL. In the instance with a low acceptance rate (11 zones, 6 vehicles), the improvement potential is similarly limited, as vehicles are rarely idle. Consequently, the performance of $\text{COMA}^\text{scd}$ is most similar to LRA in this instance. In contrast, the instance of 11 zones and 24 vehicles has a balanced ratio between the number of vehicles and requests. Here, the performance improvement of $\text{COMA}^\text{scd}$ is the largest of all instances, outperforming LRA by on average 1.9% and greedy by on average 3.5%. In the large instance (38 zones, 100 vehicles), $\text{COMA}^\text{scd}$ significantly improves performance by on average 0.6% compared to LRA, proving that the algorithm is applicable to large-scale environments. The lower performance improvement can be explained by the weight of $\text{COMA}^\text{adj}$ being required to increase more slowly in the loss function of $\text{COMA}^\text{scd}$ when the number of agents increases.

§.§ Ablation study

Relative test performance $\Delta\ [\%]$ of all algorithms vs. $\text{COMA}^\text{scd}$.

5 zones     11 zones     11 zones     11 zones     38 zones

15 veh.     6 veh.     18 veh.     24 veh.     100 veh.
GRA             $\ovoid$     $\ovoid$     –

$\text{COMA}^\text{equ}$             $\ovoid$     $\ovoid$     –

$\text{COMA}^\text{tgt}$             $\ovoid$     $\ovoid$     –

$\text{COMA}^\text{adj}$             $\ovoid$     $\ovoid$     –

LGRA                     –

Convergence of algorithms.   denotes stable convergence, $\ovoid$ unstable convergence (across random seeds),  –  no convergence.

2c11 zones

18 vehicles     24 vehicles
GRA     1.2% (0.02)     -2.7% (0.01)

$\text{COMA}^\text{equ}$     0.7% (0.03)     -3.5% (0.04)

$\text{COMA}^\text{tgt}$     2.3% (0.00)     -1.5% (0.23)

$\text{COMA}^\text{adj}$     0.3% (0.42)     -1.2% (0.41)

LGRA     2.3% (0.00)     -1.3% (0.55)

Test performance of algorithms vs. $\text{COMA}^\text{scd}$ (Wilcoxon p-value).

In the following, we discuss the performance of all proposed algorithms with respect to numerical stability and scalability across random seeds (Table <ref>) as well as computational performance (Table <ref> & Figure <ref>). As can be seen in Table <ref>, all algorithms show stable convergence for the small instances, while all but LRA, LGRA, and $\text{COMA}^\text{scd}$ exhibit stability issues already for medium-sized instances, failing to converge for one-third of the seeds and requiring about ten times as many training steps to converge compared to LRA for the remaining seeds. In contrast, LGRA and $\text{COMA}^\text{scd}$ converge on all seeds and require comparable to at maximum twice as many training steps compared to LRA. For the large instance, only $\text{COMA}^\text{scd}$ and the LRA baseline converge, with $\text{COMA}^\text{scd}$ again requiring similar to twice as many training steps. Figure <ref> and Table <ref> show the relative performance of all algorithms compared to $\text{COMA}^\text{scd}$ for the medium-sized instances over seeds for which all algorithms converged. As can be seen, the results are mixed: while pure global-rewards-based algorithms outperform $\text{COMA}^\text{scd}$ on average on the 18 vehicles instance, $\text{COMA}^\text{scd}$ outperforms all other algorithms on the 24 vehicles instance. To understand this ambiguous effect, we need to detail the algorithms' convergence behavior: increasing the instance from 18 to 24 vehicles technically requires to train 570 instead of 430 agents, which significantly challenges all purely global-rewards-based algorithms. In fact, if one looks at the variance of the validation reward of each algorithm (see Appendix <ref>), we observe converging but less stable learning behavior for all purely global-rewards-based algorithms, which explains the respective performance drop. While this observation manifests the robustness of $\text{COMA}^\text{scd}$ at a performance that improves  upon local-rewards-based algorithms, it also points at a promising direction for future research: if one manages to stabilize and scale the purely global-rewards-based algorithms, one will most likely obtain even better performing algorithms.

§.§ Structural analysis

Finally, we aim to understand the performance difference between the studied algorithms by analyzing the respective policy characteristics in Table <ref>, which details request rejection rates for our LRA baseline, the pure global-rewards-based algorithm $\text{COMA}^\text{adj}$, and $\text{COMA}^\text{scd}$. As can be seen, both COMA-based algorithms have a lower rejection rate compared to LRA, which explains their improved performance. This finding, as well as the relation between $\text{COMA}^\text{scd}$ and $\text{COMA}^\text{adj}$, are in line with the performance shown for the 24 vehicles instance in Figure <ref>.
To understand operational intricacies, we further analyze the difference between average rejection rates of empty destination zones and zones that contain more than two vehicles upon a request's arrival, which is 3.8% for LRA, 6.4% for $\text{COMA}^\text{adj}$ and 4.9% for $\text{COMA}^\text{scd}$. This indicates a stronger focus of $\text{COMA}^\text{scd}$ and $\text{COMA}^\text{adj}$ on implicit vehicle balancing, as these algorithms are more reluctant to send vehicles to already crowded zones. Such a focus on vehicle balancing stems from global reward structures and partially explains performance improvements: if vehicles are unbalanced, less overall requests can be served; individual vehicles might still obtain high local rewards, but the global reward decreases.

Beyond implicit balancing, we analyze the algorithms' anticipative performance, i.e., their capability to foresee future demand and consider it during decision-making. To do so, we analyze an algorithm's overperformance ratio (see Table <ref>), which compares the summed theoretical profits of future requests in the same zone following acceptance or rejection of an initially profitable request. We calculate this ratio by dividing the total theoretical profit after rejections by that after acceptances (see Appendix <ref> for details). As can be seen, $\text{COMA}^\text{adj}$ has the highest overperformance ratio, followed by LRA and $\text{COMA}^\text{scd}$. From the higher ratio, we conclude that $\text{COMA}^\text{adj}$ has better forecasting abilities, as requests after rejections are more profitable than requests after acceptances under its policy. In contrast, $\text{COMA}^\text{scd}$ appears to suffer from the mixture of local and global rewards, which can explain some of the performance gaps compared to algorithms with purely global rewards. A possible reason for the better forecasting abilities of $\text{COMA}^\text{adj}$ is that global Q-values incorporate more prescriptive information: since global rewards are less dependent on the actions of individual agents, it is easier to infer information about demand from them, such that agents trained using global Q-values might have access to better demand predictions.

measure     LRA     $\text{COMA}^\text{adj}$     $\text{COMA}^\text{scd}$

rejection rate     17.6$\%$     16.6$\%$     15.5$\%$

$\to$ rejection rate for destination zones without vehicles     16.4$\%$     12.9$\%$     13.6$\%$

$\to$ rejection rate for destination zones with $>$2 vehicles     20.2$\%$     19.3$\%$     18.5$\%$

ratio of overperformance rejections / acceptances     1.75     1.87     1.27

Rejection rates of generally profitable requests on the instance with 24 vehicles.

§ CONCLUSION

We study vehicle dispatching for AMoD systems, where a profit-maximizing central operator assigns vehicles to requests or rejects these. We propose a novel MADRL algorithm with global rewards and credit assignment based on a counterfactual baseline. Our algorithm combines the benefits of low learning variance and sample efficiency of SAC with the benefits of credit assignment of COMA. We show that a naive combination of SAC and COMA does not converge, and develop a stable and scalable algorithm that uses reward scheduling with a new advantage function based on an equally-weighted baseline and a target actor network baseline. We show that our algorithm improves upon the current state-of-the-art by up to 2% on average across test dates and up to 6% on individual dates from real-world data sets. This constitutes a substantial performance improvement for the AMoD application case where a single percent improvement yields significant daily monetary gains. We further provide a structural analysis which shows that the use of global rewards can improve implicit vehicle balancing and demand forecasting abilities. Our proposed algorithm is applicable beyond the area of AMoD, as it can be useful in any application where stable multi-agent deep reinforcement learning with credit assignment is required. Future work may apply our algorithm to other areas or focus on stabilizing purely global-rewards-based algorithms.

§ PROOF OF PROPOSITION <REF>

The combined loss function of SAC and COMA for one instance of a batch and one agent is

\begin{equation*}
J_\pi(\phi|s,i)=\sum_{a_i}\pi(a_i) \left(\alpha\log\pi(a_i)-Q(a_i)+\sum_{a'_i}\pi(a'_i)\: Q(a'_i) \right).
\end{equation*}

Considering the last summand separately, we get

\begin{align*}
&\sum_{a_i}\pi(a_i)\sum_{a'_i}\pi(a'_i)\: Q(a'_i) \\
&=\pi_1\:(\pi_1\: Q_1+\pi_2\: Q_2+...+\pi_n\: Q_n)
+...+\pi_n\:(\pi_1\: Q_1+\pi_2\: Q_2+...+\pi_n\: Q_n) \\&=\pi^2_1\:Q_1+\pi_1\:\pi_2\:Q_2+...+\pi_1\:\pi_n\:Q_n+\pi_2\:\pi_1\:Q_1+...+\pi_n\:\pi_1\:Q_1+...+\pi^2_n\:Q_n \\
&=\pi_1\:Q_1\:(\pi_1+\pi_2+...+\pi_n)+...+\pi_n\:Q_n\:(\pi_1+\pi_2+...+\pi_n) \\
&=\sum_{a_i}\pi(a_i)\: Q(a_i).
\end{align*}

Inserting that into the loss function, we obtain

\begin{equation*}
J_\pi(\phi|s,i)=\sum_{a_i}\pi(a_i)\: \alpha\log\pi(a_i)-\sum_{a_i}\pi(a_i)\: Q(a_i)+\sum_{a_i}\pi(a_i)\: Q(a_i)=\sum_{a_i}\pi(a_i)\: \alpha\log\pi(a_i),
\end{equation*}

which is just the entropy. If the loss function equals the entropy, the actor is not trained to maximize the probabilities of actions associated with the highest Q-values and thus cannot learn a meaningful policy. Note that this problem occurs for an arbitrary number of $n$ actions per agent (as displayed), including our setting with $n=2$ actions.

§ EXPERIMENTS

§.§ Experimental setup

Data: We use the same experimental setup as <cit.>: our data consists of yellow taxi trip records from 2015, excluding weekends and holidays. We define the time of request placement as the pickup time in the data set, and consider a one hour time interval per day as one episode, using a time step size of one minute. We discretize the island of Manhattan using small or large hexagonal zones. Neighboring zones' centers have a distance of 459 or 917 meters for the small and large zones, respectively. In our instances, 38 large zones roughly cover the southern two thirds of Manhattan, 11 small zones cover roughly three square kilometers in central Manhattan and 5 small zones denote the southern five of these zones, covering close to half of that area. For each instance, we delete all requests originating in or leading to zones outside the area under consideration and all requests that start and end within the same zone. Vehicles travel along the shortest route between zones. We assume a travel time of two time steps between neighboring small zones and five time steps between neighboring large zones. The maximum waiting time of customers between request submission and pickup is five minutes. We set the revenue to 5.00 USD per km and the costs to 4.50 USD per km, resulting in a profit margin of 10% in ideal cases. The maximum number of requests per time step never exceeds 7 for the instances with 5 small zones and 24 for the instances of 11 small zones. For the instances of 38 large zones, we use only every  $\text{20}^{\text{th}}$ request to make the problem size tractable for our available hardware. This results in a maximum of 20 requests per time step in these instances.

Hardware: We train the models on NVIDIA Tesla P100, V100, or A100 GPUs. Depending on the problem size and number of time steps for training, this training takes between 2 and 66 hours.

Benchmark: The greedy algorithm considers all requests in the order of their arrival. If a vehicle is available to serve the request with a positive profit and within the maximum waiting time, it is assigned to serve the request. If there are multiple available vehicles, the greedy algorithm assigns the vehicle that can serve the request with the highest profit. If no vehicle is available, the greedy algorithm rejects the request. In the instance with 5 zones and 15 vehicles, the greedy algorithm accepts 78% of requests. This acceptance rate is 22% in the instance with 11 zones and 6 vehicles, 47% in the instance with 11 zones and 18 vehicles, 55% in the instance with 11 zones and 24 vehicles, and 50% in the instance with 38 zones and 100 vehicles.

§.§ Hyperparameters

We mostly use the same hyperparameters as <cit.> for our algorithm, only adjusting the entropy coefficient $\alpha$, the learning rate, and the total number of training steps. A different number of random steps in the beginning, size of the replay buffer, batch size, or discount factor did not improve our results.

We train models until their validation performance is stable and does not increase further. Depending on the instance and algorithm, this results in 200,000 to 2,000,000 training steps. The number of required steps increases with the instance size, i.e., with the number of zones and vehicles. Across instances, LRA, $\text{COMA}^\text{scd}$ and LGRA require 200,000 to 400,000 steps, while GRA, $\text{COMA}^\text{equ}$, $\text{COMA}^\text{tgt}$ and $\text{COMA}^\text{adj}$ require 200,000 to 2,000,000 steps. We separate the learning rates for actor and critic, as global rewards and the inclusion of a baseline can change the models' convergence properties. Tuning the learning rates for each instance, we set both learning rates to $3\cdot10^{-4}$ for LRA and to values between $3\cdot10^{-4}$ and $2.4\cdot10^{-3}$ for GRA. For the algorithms based on COMA ($\text{COMA}^\text{equ}$, $\text{COMA}^\text{tgt}$, $\text{COMA}^\text{adj}$, $\text{COMA}^\text{scd}$), we set the learning rates to values between $6\cdot10^{-4}$ and $1.8\cdot10^{-3}$, with the actor having higher learning rates than the critic. We further tune the entropy coefficient $\alpha$ for each instance, with values ranging between 0.37 and 1.5.

At the end of each training run, we test the model with the best validation performance on the test data. The score of a single test date is the sum of all individual scores on that date. We repeat each training run with three different random seeds and average test scores over these random seeds. If a run does not converge to a reasonable performance, we exclude it from the average. This only applies to some of the algorithms in the ablation study in Section <ref>. All test scores in this paper are obtained in that way.

§.§ Training schedules

The policy loss function of $\text{COMA}^\text{adj}$ is a dynamic weighted average of the loss functions of $\text{COMA}^\text{equ}$ and $\text{COMA}^\text{tgt}$, with $\beta$ being the weight of the loss function of $\text{COMA}^\text{tgt}$. Over the course of the training, $\beta$ always starts at zero, ends at one, and can increase linearly or following a power function. Since $\text{COMA}^\text{tgt}$ shows better experimental performance than $\text{COMA}^\text{equ}$ (see Section <ref>), we test a power function with an exponent of 0.5 to increase the share of $\text{COMA}^\text{tgt}$ quicker than when using a linear schedule. In the instance of 11 zones and 18 vehicles, the maximum average performance of a model with a power function is 4.6% lower than the maximum average performance of a model with a linear function. In the instance of 11 zones and 24 vehicles, this number is 2.8%. We conclude that a linearly increasing $\beta$ works best.

The policy loss function of $\text{COMA}^\text{scd}$ is a dynamic weighted average of the loss functions of LRA and $\text{COMA}^\text{adj}$, with $\kappa$ being the weight of the loss function of $\text{COMA}^\text{adj}$. Over the course of the training, $\kappa$ always starts at zero and ends at one. It can increase linearly, following a power function, or it can jump from zero to one at any specified point during training. We present results of tests for different patterns of increasing $\kappa$ in Table <ref>. We draw four conclusions from these tests: firstly, power functions generally lead to the best performance, which is most likely the result of their smoother transition from local to global rewards in comparison to sudden jumps. Secondly, quickly increasing power functions with exponents between 0.01 and 0.5 work best in all instances. As a result of these quickly increasing functions, the share of $\text{COMA}^\text{adj}$ in the loss function rises faster in the best models than it would when following a linear schedule. This shows that a quick transition to global rewards enables utilizing their benefits for a larger part of the training. Thirdly, the best exponents vary between instances, i.e., a per-instance tuning of the schedule to increase $\kappa$ can lead to visible performance increases. Finally, the best performing exponent increases with the instance size. We thus conclude that large instances require a slower transition from local to global rewards to achieve best performance, in line with our results concerning the scalability problems of global-rewards-based algorithms.

instance     schedule type     exponent/jump     vs. best
4*5 zones, 15 veh.     power     0.01     0.0%

power     0.05     -1.5%

power     0.25     -1.6%

jump     0.01     -1.9%
5*11 zones, 6 veh.     power     0.25     0.0%

jump     0.01     -0.7%

power     0.125     -1.6%

jump     0.02     -2.1%

power     0.01     -4.4%
7*11 zones, 18 veh.     power     0.25     0.0%

jump     0.1     -0.4%

jump     0.25     -0.4%

power     2.00     -0.6%

power     1.00     -1.2%

power     0.50     -2.1%

jump     0.50     -3.1%
4*11 zones, 24 veh.     power     0.25     0.0%

jump     0.125     -1.0%

power     0.125     -1.6%

power     1.00     -1.7%
5*38 zones, 100 veh.     power     0.5     0.0%

power     1.00     -0.2%

power     0.25     -0.3%

jump     0.25     -0.3%

Performance of training schemes of $\text{COMA}^\text{scd}$, relative to the best performance per instance. The specifications used in Sections <ref> to <ref> are displayed in bold writing. The column exponent/jump denotes the exponent in case of a power function or the jump point relative to the number of training steps in case of a jump schedule. For example, if the jump point is 0.25 and the model is trained for 400,000 steps, $\kappa$ is zero for the first 100,000 steps and one for the remainder of the training.

§.§ Validation rewards

We show the validation rewards of GRA, $\text{COMA}^\text{equ}$, $\text{COMA}^\text{tgt}$, $\text{COMA}^\text{adj}$, LGRA, LRA, and $\text{COMA}^\text{scd}$ for the instance with 11 zones and 18 vehicles and the instance with 11 zones and 24 vehicles over the course of training in Figure <ref>. We observe three patterns in the validation rewards: firstly, in both instances, the purely global-rewards-based algorithms GRA, $\text{COMA}^\text{equ}$, $\text{COMA}^\text{tgt}$, and $\text{COMA}^\text{adj}$ need about ten times as many training steps to converge as LRA, while LGRA and $\text{COMA}^\text{scd}$ need about the same to double the number of steps until convergence. Secondly, most purely global-rewards-based algorithms display far larger differences between their best and worst validation performance in both instances, indicating unstable convergence behavior. Finally, the convergence speed decreases for the purely global-rewards-based algorithms except $\text{COMA}^\text{adj}$ in the instance with 24 vehicles compared to the one with 18 vehicles. For $\text{COMA}^\text{adj}$, the maximum and minimum performance diverges more in the larger instance. The algorithm LGRA also faces challenges when increasing the number of agents: to remain stable and converge quickly, we have to decrease the share of global rewards from 60% to 30%, thereby also lowering the positive influence of global rewards. In contrast, LRA and $\text{COMA}^\text{scd}$ have about the same stability and convergence speed in both instances.

This final observation provides evidence for the lower performance of purely global-rewards-based algorithms in the instance with 24 vehicles. As the learning is less stable, the learned policies are less reliable. With even validation performance curves of converging models being less stable, trained models are more likely to converge to a sub-optimal policy. While we allow all algorithms enough steps to converge, the slower learning of purely global-rewards-based algorithms can be an additional problem in practice. Overall, these results confirm our conclusion that learning using purely global rewards increases vehicle dispatching performance, but leads to problems when increasing the number of agents.

Validation rewards of algorithms over training steps. The black horizontal line indicates zero. The main line denotes the average validation reward over the three random seeds, the shaded area the maximum and minimum rewards at each step.

§.§ Calculation of overperformance ratio

We calculate the overperformance ratio as follows: firstly, we store the maximum profit of a request if it is positive together with the origin zone of this request. This is the profit that would be obtained if the request was served by the closest vehicle. Secondly, we store the theoretical profits that would be obtained from subsequent requests originating in the same zone if a vehicle was available at the same position as the vehicle closest to the original request. We store the theoretical profits for the ten time steps after the original request appears. Thirdly, we obtain the overprofits by subtracting the profit of the original request from the theoretical profits. We only count positive overprofits subsequently, as we only consider requests that would be more profitable than the original request. We sum these overprofits conditional on whether the original request is accepted or rejected. Finally, we divide the overprofit after a rejection by the overprofit after an acceptance.

A numerical illustration is as follows: consider a request, which could be served by a vehicle in the same zone with an initial profit of 10 USD, but the operator rejects this request. From the same zone, three further requests appear within the ten subsequent time steps. Using a vehicle within the same zone, they could be served with a profit of 5 USD, 12 USD, and 14 USD, respectively. These numbers are the theoretical profits. The associated overprofits are 0 USD, 2 USD, and 4 USD. Therefore, the total overprofit compared to the original request is 6 USD. Now, consider a second request, which is accepted by the operator, with an initial profit of 10 USD and subsequent theoretical profits of 8 USD, 11 USD, and 12 USD. The overperformance for this request is 3 USD. If these are the only two profitable requests available, the overperformance ratio is 6 divided by 3, which equals 2. As the overperformance ratio is larger than one, the subsequent requests after the rejection are more profitable than those after the acceptance. Accordingly, whenever the overperformance rato is larger than one, we can conclude that the analyzed algorithm is better at taking anticipative decisions than a greedy algorithm.

We thank the CORAIL research group at Polytechnique Montréal for valuable comments and discussions. The work of Heiko Hoppe was supported by a scholarship from the Max Weber-Programm Bayern.