# Standard Dual Quaternion Optimization and Its Applications in Hand-Eye
Calibration and SLAM

Liqun Qi111Department of Applied Mathematics, The Hong Kong Polytechnic
University, Hung Hom, Kowloon, Hong Kong; Department of Mathematics, School of
Science, Hangzhou Dianzi University, Hangzhou 310018 China
(maqilq@polyu.edu.hk).

###### Abstract

Several common dual quaternion functions, such as the power function, the
magnitude function, the $2$-norm function and the $k$th largest eigenvalue of
a dual quaternion Hermitian matrix, are standard dual quaternion functions,
i.e., the standard parts of their function values depend upon only the
standard parts of their dual quaternion variables. Furthermore, the sum,
product, minimum, maximum and composite functions of two standard dual
functions, the logarithm and the exponential of standard unit dual quaternion
functions, are still standard dual quaternion functions. On the other hand,
the dual quaternion optimization problem, where objective and constraint
function values are dual numbers but variables are dual quaternions, naturally
arises from applications. We show that to solve an equality constrained dual
quaternion optimization problem, we only need to solve two quaternion
optimization problems. If the involved dual quaternion functions are all
standard, the optimization problem is called a standard dual quaternion
optimization problem, and some better results hold. Then, we show that the
dual quaternion optimization problems arising from the hand-eye calibration
problem and the simultaneous localization and mapping (SLAM) problem are
equality constrained standard dual quaternion optimization problems.

Key words. Standard dual quaternion functions, dual quaternion optimization,
quaternion optimization, hand-eye calibration, simultaneous localization and
mapping.

##  1 Introduction

Dual quaternions have wide applications in robotics, 3D motion modelling and
control, and computer graphics [1, 2, 4, 9, 10, 13, 14, 15, 16, 17, 23].

According to [18], the magnitudes of dual quaternions and the $2$-norms of
dual quaternion vectors are dual numbers. A total order was defined for dual
numbers in [18]. Thus, dual quaternion optimization problems, where objective
and constraint function values are dual numbers but variables are dual
quaternions, naturally arise. In particular, the least squares approach
results in such optimization problems. In Section 5, we will present such dual
quaternion optimization problems arising from hand-eye calibration. However,
how to solve a dual quaternion optimization problem is still a problem. First,
there is no calculus concepts for dual quaternion functions. Even the limit
concept is not established for dual quaternion functions. Second, unlike
quaternion optimization problems, which can be transformed to one-level real
optimization problems, general dual quaternion optimization problems are
essentially bilevel optimization problems. These are difficulties.

Fortunately, some common dual quaternion functions are easy to handle. Several
common dual quaternion functions, such as the power function, the magnitude
function, the $2$-norm function and the $k$th largest eigenvalue of a dual
quaternion Hermitian matrix, are standard dual quaternion functions. The
standard parts of the function values of these functions depend upon only the
standard parts of their dual quaternion variables. Furthermore, the sum,
product, minimum, maximum and composite functions of two standard dual
functions, the logarithm and the exponential of a standard unit dual
quaternion functions, are still standard dual quaternion functions. Then we
show that to solve an equality constrained dual quaternion optimization
problem, we only need to solve two quaternion optimization problems. For a
dual quaternion optimization problem, if its objective and constraint
functions are all standard, then we call it a standard dual quaternion
optimization problem, and some better results hold there. Finally, we show
that the dual quaternion optimization problems arising from the hand-eye
calibration problem are equality constrained standard dual quaternion
optimization problems, and thus relatively easier to handle. This provides a
method to solve this problem.

In the next section, we present some preliminary knowledge on dual numbers,
dual quaternions and quaternion optimization. Standard dual quaternion
functions are studied in Section 3. A solution method for the equality
constrained dual quaternion optimization problem is proposed, and standard
dual quaternion optimization is discussed in Section 4. In Section 5, we show
that equality constrained standard dual quaternion optimization problem has
some better properties. Finally, in Sections 6 and 7, we show that the dual
quaternion optimization problems arising from the hand-eye calibration problem
and the simultaneous localization and mapping (SLAM) problem are equality
constrained standard dual quaternion optimization problems.

Scalars, vectors and matrices are denoted by small letters, bold small letters
and capital letters, respectively. Dual numbers and dual quaternions are
distinguished by a hat symbol.

##  2 Dual Numbers, Quaternions and Dual Quaternions

###  2.1 Dual Numbers

The set of the dual numbers is denoted as $\hat{\mathbb{R}}$. Following the
literature such as [23], we use the hat symbol to denote dual numbers and dual
quaternions. A dual number $\hat{q}$ has the form $\hat{q}=q+q_{d}\epsilon$,
where $q$ and $q_{d}$ are real numbers, and $\epsilon$ is the infinitesimal
unit, satisfying $\epsilon^{2}=0$. The quaternion $q$ is called the real part
or the standard part of $\hat{q}$, and the quaternion $q_{d}$ is called the
dual part or the infinitesimal part of $\hat{q}$. The infinitesimal unit
$\epsilon$ is commutative in multiplication with real numbers, complex numbers
and quaternion numbers. If $q\not=0$, $\hat{q}$ is said to be appreciable,
otherwise, $\hat{q}$ is said to be infinitesimal.

A total order was introduced in [18] for dual numbers. Given two dual numbers
$\hat{p},\hat{q}\in\hat{\mathbb{R}}$, $\hat{p}=p+p_{d}\epsilon$,
$\hat{q}=q+q_{d}\epsilon$, where $p$, $p_{d}$, $q$ and $q_{d}$ are real
numbers, we say that $\hat{p}\leq\hat{q}$, if either $p<q$, or $p=q$ and
$p_{d}\leq q_{d}$. In particular, we say that $\hat{p}$ is positive,
nonnegative, nonpositive or negative, if $\hat{p}>0$, $\hat{p}\geq 0$,
$\hat{p}\leq 0$ or $\hat{p}<0$, respectively.

###  2.2 Quaternions

The set of the quaternions is denoted by $\mathbb{Q}$. A quaternion $q$ has
the form $q=q_{0}+q_{1}\mathbf{i}+q_{2}\mathbf{j}+q_{3}\mathbf{k},$ where
$q_{0},q_{1},q_{2}$ and $q_{3}$ are real numbers, $\mathbf{i},\mathbf{j}$ and
$\mathbf{k}$ are three imaginary units of quaternions, satisfying

$\mathbf{i}^{2}=\mathbf{j}^{2}=\mathbf{k}^{2}=\mathbf{i}\mathbf{j}\mathbf{k}=-1,~{}~{}\mathbf{i}\mathbf{j}=-\mathbf{j}\mathbf{i}=\mathbf{k},~{}~{}\mathbf{j}\mathbf{k}=-\mathbf{k}\mathbf{j}=\mathbf{i},~{}~{}\mathbf{k}\mathbf{i}=-\mathbf{i}\mathbf{k}=\mathbf{j}.$

The real part of $q$ is Re$(q)=q_{0}$. The imaginary part of $q$ is
Im$(q)=q_{1}\mathbf{i}+q_{2}\mathbf{j}+q_{3}\mathbf{k}$. The multiplication of
quaternions satisfies the distribution law, but is noncommutative.

The conjugate of $q=q_{0}+q_{1}\mathbf{i}+q_{2}\mathbf{j}+q_{3}\mathbf{k}$ is
$q^{*}:=q_{0}-q_{1}\mathbf{i}-q_{2}\mathbf{j}-q_{3}\mathbf{k}.$ The magnitude
of $q$ is $|q|=\sqrt{q_{0}^{2}+q_{1}^{2}+q_{2}^{2}+q_{3}^{2}}.$ It follows
that the inverse of a nonzero quaternion $q$ is $q^{-1}={q^{*}/|q|^{2}}.$ For
any two quaternions $p$ and $q$, we have $(pq)^{*}=q^{*}p^{*}$.

A quaternion is said imaginary if its real part is zero. If $q$ is imaginary,
then $q^{*}=-q$. In the literature [23], it is called a vector quaternion.
Various 3D vectors, such as position vectors, displacement vectors, linear
velocity vectors, and angular velocity vectors, can be represented as
imaginary quaternions.

If $|q|=1$, then $q$ is called a unit quaternion, or a rotation quaternion. A
spatial rotation around a fixed point of $\theta$ radians about a unit axis
$(x_{1},x_{2},x_{3})$ that denotes the Euler axis is given by the unit
quaternion

$q=\cos(\theta/2)+x_{1}\sin(\theta/2)\mathbf{i}+x_{2}\sin(\theta/2)\mathbf{j}+x_{3}\sin(\theta/2)\mathbf{k}=e^{{\theta\over
2}x}.$ (1)

where the unit axis $x$ is an imaginary unit quaternion
$x=x_{1}\mathbf{i}+x_{2}\mathbf{j}+x_{3}\mathbf{k}$. We may also write

$\ln q={\theta\over 2}x.$ (2)

A unit quaternion $q$ is always invertible and $q^{-1}=q^{*}$.

The collection of $n$-dimensional quaternion vectors is denoted by
${\mathbb{Q}}^{n}$. For
$\mathbf{x}=(x_{1},x_{2},\cdots,x_{n})^{\top},\mathbf{y}=(y_{1},y_{2},\cdots,y_{n})^{\top}\in{\mathbb{Q}}^{n}$,
define $\mathbf{x}^{*}\mathbf{y}=\sum_{i=1}^{n}x_{i}^{*}y_{i}$, where
$\mathbf{x}^{*}=(x_{1}^{*},x_{2}^{*},\cdots,x_{n}^{*})$ is the conjugate
transpose of $\mathbf{x}$.

###  2.3 Dual Quaternions and Unit Dual Quaternions

The set of dual quaternions is denoted by $\hat{\mathbb{Q}}$. A dual
quaternion $\hat{q}\in\hat{\mathbb{Q}}$ has the form

$\hat{q}=q+q_{d}\epsilon,$ (3)

where $q,q_{d}\in\mathbb{Q}$ are the standard part and the dual part of
$\hat{q}$, respectively. If $q\not=0$, then we say that $\hat{q}$ is
appreciable. If $q$ and $q_{d}$ are imaginary quaternions, then $\hat{q}$ is
called an imaginary dual quaternion.

The conjugate of $\hat{q}$ is

$\hat{q}^{*}=q^{*}+q_{d}^{*}\epsilon.$ (4)

Thus, if $\hat{q}=\hat{q}^{*}$, then $\hat{q}$ is a dual number. If $\hat{q}$
is imaginary, then $\hat{q}^{*}=-\hat{q}$.

The magnitude of $\hat{q}$ was defined in [18] as

$|\hat{q}|:=\left\\{\begin{aligned} |q|+{(qq_{d}^{*}+q_{d}q^{*})\over
2|q|}\epsilon,&\ {\rm if}\ q\not=0,\\\ |q_{d}|\epsilon,&\ {\rm
otherwise},\end{aligned}\right.$ (5)

which is a dual number.

For two dual quaternions $\hat{p}=p+p_{d}\epsilon$ and
$\hat{q}=q+q_{d}\epsilon$, their addition and multiplications are defined as

$\hat{p}+\hat{q}=\left(p+q\right)+\left(p_{d}+q_{d}\right)\epsilon$

and

$\hat{p}\hat{q}=pq+\left(pq_{d}+p_{d}q\right)\epsilon.$

A dual number is always commutative with a dual quaternion or a dual
quaternion vector.

A dual quaternion $\hat{q}$ is called invertible if there exists a quaternion
$\hat{p}$ such that $\hat{p}\hat{q}=\hat{q}\hat{p}=1$. A dual quaternion
$\hat{q}$ is invertible if and only if $\hat{q}$ is appreciable. In this case,
we have

$\hat{q}^{-1}=q^{-1}-q^{-1}q_{d}q^{-1}\epsilon.$

If $|\hat{q}|=1$, then $\hat{q}$ is called a unit dual quaternion. A unit dual
quaternion $\hat{q}$ is always invertible and we have
${\hat{q}}^{-1}={\hat{q}}^{*}$. The 3D motion of a rigid body can be
represented by a unit dual quaternion. We have

$\hat{q}\hat{q}^{*}=(q+q_{d}\epsilon)(q^{*}+q_{d}^{*}\epsilon)=qq^{*}+(qq_{d}^{*}+q_{d}q^{*})\epsilon=\hat{q}^{*}\hat{q}.$

Thus, $\hat{q}$ is a unit dual quaternion if and only if $q$ is a unit
quaternion, and

$qq_{d}^{*}+q_{d}q^{*}=q^{*}q_{d}+q_{d}^{*}q=0.$ (6)

Suppose that there is a rotation $q\in{\mathbb{Q}}$ succeeded by a translation
$p^{b}\in{\mathbb{Q}}$, where $p^{b}$ is an imaginary quaternion. The
superscripts $b$ and $s$ represent the relation of the rigid body motion with
respect to the body frame attached to the rigid body and the spatial frame
which is relative to a fixed coordinate frame. Then the whole transformation
can be represented by unit dual quaternion $\hat{q}=q+q_{d}\epsilon$, where
$q_{d}={1\over 2}qp^{b}$. Note that we have

$qq_{d}^{*}+q_{d}q^{*}={1\over
2}\left[q(p^{b})^{*}q^{*}+qp^{b}q^{*}\right]={1\over
2}q\left[(p^{b})^{*}+p^{b}\right]q^{*}=0.$

Thus, a transformation of a rigid body can be represented by a unit dual
quaternion

$\hat{q}=q+{\epsilon\over 2}qp^{b},$ (7)

where $q$ is a unit quaternion to represent the rotation, and $p^{b}$ is the
imaginary quaternion to represent the translation or the position. In (7), $q$
is the attitude of the rigid body, while $\hat{q}$ represents the
transformation.

Combining (7) with (2), we have

$\ln\hat{q}={1\over 2}(\theta x+\epsilon p^{b}).$ (8)

A unit dual quaternion $\hat{q}$ serves as both a specification of the
configuration of a rigid body and a transformation taking the coordinates of a
point from one frame to another via rotation and translation. In (7), if
$\hat{q}$ is the configuration of the rigid body, then $q$ and $p^{b}$ are the
attitude of and position of the rigid body respectively.

Denote the collection of $n$-dimensional dual quaternion vectors by
${\hat{\mathbb{Q}}}^{n}$.

For $\hat{\mathbf{x}}=(\hat{x}_{1},\hat{x}_{2},\cdots,\hat{x}_{n})^{\top}$,
$\hat{\mathbf{y}}=(\hat{y}_{1},\hat{y}_{2},\cdots,\hat{y}_{n})^{\top}\in{\hat{\mathbb{Q}}}^{n}$,
define
$\hat{\mathbf{x}}^{*}\hat{\mathbf{y}}=\sum_{i=1}^{n}\hat{x}_{i}^{*}\hat{y}_{i}$,
where
$\hat{\mathbf{x}}^{*}=(\hat{x}_{1}^{*},\hat{x}_{2}^{*},\cdots,\hat{x}_{n}^{*})$
is the conjugate transpose of $\hat{\mathbf{x}}$. We say $\hat{\mathbf{x}}$ is
appreciable if at least one of its component is appreciable.

For $\hat{\mathbf{x}}=(\hat{x}_{1},\hat{x}_{2},\cdots,\hat{x}_{n})^{\top}$, by
[18], if not all of $\hat{x}_{i}$ are infinitesimal, its $2$-norm is defined
as

$\|\hat{\mathbf{x}}\|_{2}=\sqrt{\sum_{i=1}^{n}|\hat{x}_{i}|^{2}}.$ (9)

If all $\hat{x}_{i}$ are infinitesimal, we have
$\hat{x}_{i}=(x_{i})_{d}\epsilon$ for $i=1,2,\ldots,n$. Then we have

$\|\hat{\mathbf{x}}\|_{2}=\sqrt{\sum_{i=1}^{n}|(x_{i})_{d}|^{2}}\epsilon.$
(10)

##  3 Standard Dual Quaternion Functions

We call a function $\hat{f}:{\hat{\mathbb{Q}}}^{n}\to{\hat{\mathbb{R}}}$ a
dual quaternion function. Let

$\hat{f}(\hat{\mathbf{x}})=f(\hat{\mathbf{x}})+\epsilon
f_{d}(\hat{\mathbf{x}})$

and

$\hat{\mathbf{x}}=\mathbf{x}+\epsilon\mathbf{x}_{d},$

where $f(\hat{\mathbf{x}}),f_{d}(\hat{\mathbf{x}})\in\mathbb{Q}$,
$\mathbf{x},\mathbf{x}_{d}\in{\mathbb{Q}}^{n}$ are the standard parts and dual
parts of $\hat{f}(\hat{\mathbf{x}})$ and $\hat{\mathbf{x}}$, respectively. If
$f(\hat{\mathbf{x}})\equiv f(\mathbf{x})$ for all
$\hat{\mathbf{x}}\in{\hat{\mathbb{Q}}}^{n}$, then we say that $\hat{f}$ is a
standard dual quaternion function.

A simple example is the power function. Consider the power function
$\hat{f}(\hat{x})=(\hat{x})^{m}$ for a positive integer $m$. Let
$\hat{x}=x+\epsilon x_{d}$. Then
$\hat{f}(\hat{x})=(\hat{x})^{m}=x^{m}+m\epsilon x^{m-1}x_{d}$, i.e.,
$f(\hat{x})=x^{m}=\hat{f}(x)$. By (5), the magnitude function is also a
standard quaternion function. By (5), (9) and (10), the $2$-norm function is
also a standard quaternion function. A fourth example is
$\hat{f}(\hat{A})=\hat{\lambda}_{k}(\hat{A})$, where $\hat{A}$ is the an
$n\times n$ dual quaternion Hermitian matrix, $\hat{\lambda}_{k}(\hat{A})$ is
the $k$th largest eigenvalue of $A$, $1\leq k\leq n$. By Theorem 4.1 of [19],
this is also a standard dual quaternion function.

Furthermore, many operations preserve standard dual quaternion functions.

######  Theorem 3.1.

Suppose that $\hat{f},\hat{g}:{\hat{\mathbb{Q}}}^{n}\to{\hat{\mathbb{R}}}$ are
two standard dual quaternion functions. Then their sum, product, minimum and
maximum functions are still standard dual quaternion functions.

###### Proof.

Let
$\hat{h}(\hat{\mathbf{x}})=\hat{f}(\hat{\mathbf{x}})+\hat{g}(\hat{\mathbf{x}})$.
Then $\hat{h}=h+\epsilon h_{d}$, where

$h(\hat{\mathbf{x}})=f(\hat{\mathbf{x}})+g(\hat{\mathbf{x}})=f(\mathbf{x})+g(\mathbf{x})=h(\mathbf{x}),$

i.e., $\hat{h}$ is also a standard dual quaternion function. This proves the
first conclusion.

Let
$\hat{h}(\hat{\mathbf{x}})=\hat{f}(\hat{\mathbf{x}})\hat{g}(\hat{\mathbf{x}})$.
Then $\hat{h}=h+\epsilon h_{d}$, where

$h(\hat{\mathbf{x}})=f(\hat{\mathbf{x}})g(\hat{\mathbf{x}})=f(\mathbf{x})g(\mathbf{x})=h(\mathbf{x}),$

i.e., $\hat{h}$ is also a standard dual quaternion function. This proves the
second conclusion.

Let
$\hat{h}(\hat{\mathbf{x}})=\min\\{\hat{f}(\hat{\mathbf{x}}),\hat{g}(\hat{\mathbf{x}})\\}$.
Then $\hat{h}=h+\epsilon h_{d}$, where

$h(\hat{\mathbf{x}})=\min\\{f(\hat{\mathbf{x}}),g(\hat{\mathbf{x}})\\}=\min\\{f(\mathbf{x}),g(\mathbf{x})\\}=h(\mathbf{x}),$

i.e., $\hat{h}$ is also a standard dual quaternion function. This proves the
third conclusion.

The fourth conclusion can be proved similarly.

∎

######  Corollary 3.2.

Suppose that
$\hat{f}_{1},\cdots,\hat{f}_{m}:{\hat{\mathbb{Q}}}^{n}\to{\hat{\mathbb{R}}}$
are $m$ standard dual quaternion functions, where $m$ is a positive integer.
Then their sum, product, minimum and maximum functions are still standard dual
quaternion functions.

######  Corollary 3.3.

Suppose that $\hat{f}:{\hat{\mathbb{Q}}}^{n}\to{\hat{\mathbb{R}}}$ is a
standard dual quaternion functions, and $m$ is a positive integer. Then
$(\hat{f})^{m}$ is still a standard dual quaternion function.

######  Theorem 3.4.

Suppose that $\hat{f}:\hat{\mathbb{Q}}\to{\hat{\mathbb{R}}}$ is a standard
unit dual quaternion function. Then its logarithm and exponential functions
are still standard dual quaternion functions.

###### Proof.

We have

$\hat{f}(\hat{x})=f(\hat{x})+\epsilon f_{d}(\hat{x})=f(x)+\epsilon
f_{d}(\hat{x}),$

as $\hat{f}$ is a standard dual quaternion function. Since $f$ is a unit dual
quaternion function, we have $|f(x)|=1$. Thus, we may write

$\hat{f}(\hat{x})=f(x)+{\epsilon\over 2}f(x)p^{b}.$

Write

$f(x)=\cos(\theta/2)+y_{1}\sin(\theta/2)\mathbf{i}+y_{2}\sin(\theta/2)\mathbf{j}+y_{3}\sin(\theta/2)\mathbf{k}.$

Then $\theta/2$ and $y=y_{1}\mathbf{i}+y_{2}\mathbf{j}+y_{3}\mathbf{k}$ are
functions of $x$. By (8), we have

$\ln\hat{f}(\hat{x})={1\over 2}(\theta y+\epsilon p^{b}).$

Thus, $\ln\hat{f}$ is still a standard dual quaternion function. Similarly, we
may show that $e^{\hat{f}}$ is a standard dual quaternion function too. ∎

######  Theorem 3.5.

Suppose that $\hat{f}:\hat{\mathbb{Q}}\to{\hat{\mathbb{R}}}$ and
$\hat{g}:\hat{\mathbb{Q}}^{n}\to{\hat{\mathbb{R}}}$ are two standard unit dual
quaternion functions. Then their composite function $\hat{h}=\hat{f}\circ
g:\hat{\mathbb{Q}}^{n}\to{\hat{\mathbb{R}}}$ is also a standard dual
quaternion function.

###### Proof.

For $\hat{\mathbf{x}}\in{\mathbb{Q}}^{n}$,

$h(\hat{\mathbf{x}})=f(\hat{g}(\hat{\mathbf{x}}))=f(g(\hat{\mathbf{x}}))=f(g(\mathbf{x}))=h(\mathbf{x}),$

where the second and the third equalities hold because $\hat{f}$ and $\hat{g}$
are standard unit dual quaternion functions respectively. Thus, $\hat{h}$ is
also a standard dual quaternion function. ∎

##  4 Equality Constrained Dual Quaternion Optimization

In this section, we propose a solution method for the following equality
constrained dual quaternion optimization problem (EQDQO):

$\min\\{\hat{f}(\hat{\mathbf{x}}):\hat{h}_{j}(\hat{\mathbf{x}})=0,\
j=1,\cdots,m\\},$ (11)

where $\hat{\mathbf{x}}\in{\hat{\mathbb{Q}}}^{n}$, $\hat{f}$ and $\hat{h}_{j}$
for $j=1,\cdots,m$ are dual quaternion functions. If $\hat{f}$ and
$\hat{h}_{j}$ for $j=1,\cdots,m$ are all standard dual quaternion functions,
then this problem is called an equality constrained standard dual quaternion
optimization problem (EQSDQO). EQSDQO arises in many applications, as the most
dual quaternion application problems are unit dual quaternion application
problems, and their constraints may only have the unit requirement. In the
next section, we will see that the dual quaternion optimization problem
arising from the hand-eye calibration problem is an EQSDQO.

We have
$\hat{f}(\hat{\mathbf{x}})=f(\mathbf{x}+\mathbf{x}_{d}\epsilon)+f_{d}(\mathbf{x}+\mathbf{x}_{d}\epsilon)\epsilon$.
We may regard $f$ and $f_{d}$ as real valued functions $\bar{f}$ and
$\bar{f}_{d}$, with $2n$-dimensional quaternion vector variables
$(\mathbf{x},\mathbf{x}_{d})$, i.e.,

$\bar{f}(\mathbf{x},\mathbf{x}_{d})\equiv
f(\mathbf{x}+\mathbf{x}_{d}\epsilon),\
\bar{f}_{d}(\mathbf{x},\mathbf{x}_{d})\equiv
f_{d}(\mathbf{x}+\mathbf{x}_{d}\epsilon).$

We now abuse the notation, and simply write

$f(\mathbf{x},\mathbf{x}_{d})\equiv f(\mathbf{x}+\mathbf{x}_{d}\epsilon),\
f_{d}(\mathbf{x},\mathbf{x}_{d})\equiv
f_{d}(\mathbf{x}+\mathbf{x}_{d}\epsilon),$

i.e., we now simply regard $f$ and $f_{d}$ as real valued functions, with
$2n$-dimensional quaternion vector variables $(\mathbf{x},\mathbf{x}_{d})$. We
may treat $h_{j}$ and $(h_{j})_{d}$ similarly. Then we have the following
theorem.

######  Theorem 4.1.

Suppose that $L^{I}$ is the optimal function value of the equality constrained
quaternion optimization problem (EQQOPI):

$\min\\{f(\mathbf{x},\mathbf{x}_{d}):h_{j}(\mathbf{x},\mathbf{x}_{d})=0,(h_{j})_{d}(\mathbf{x},\mathbf{x}_{d})=0,\
j=1,\cdots,m\\},$ (12)

and $(\mathbf{x}^{opt},\mathbf{x}_{d}^{opt})$ is a global optimal solution of
the following equality constrained quaternion optimization problem (EQQOPII):

$\min\\{f_{d}(\mathbf{x},\mathbf{x}_{d}):f(\mathbf{x},\mathbf{x}_{d})=L^{I},\
h_{j}(\mathbf{x},\mathbf{x}_{d})=0,\
(h_{j})_{d}(\mathbf{x},\mathbf{x}_{d})=0,\ j=1,\cdots,m\\}.$ (13)

Then $\hat{\mathbf{x}}^{opt}=\mathbf{x}^{opt}+\epsilon\mathbf{x}_{d}^{opt}$ is
a global optimal solution of EQDQO (11).

###### Proof.

Since $\hat{\mathbf{x}}^{opt}=\mathbf{x}^{opt}+\epsilon\mathbf{x}_{d}^{opt}$
satisfies all the constraints of EQDQO, and the standard part and the dual
part of $f$ attain their global minimal values, the conclusion follows. ∎

Let
$\mathbf{x}=\mathbf{x}_{0}+\mathbf{x}_{1}\mathbf{i}+\mathbf{x}_{2}\mathbf{j}+\mathbf{x}_{3}\mathbf{k}$
and
$\mathbf{x}_{d}=\mathbf{x}_{4}+\mathbf{x}_{5}\mathbf{i}+\mathbf{x}_{6}\mathbf{j}+\mathbf{x}_{7}\mathbf{k}$.
Then we may further regard $f(\mathbf{x},\mathbf{x}_{d})$,
$f_{d}(\mathbf{x},\mathbf{x}_{d})$, $h_{j}(\mathbf{x},\mathbf{x}_{d})$ and
$(h_{j})_{d}(\mathbf{x},\mathbf{x}_{d})$ as real valued functions with real
vector variables $\mathbf{x}_{k}$ for $k=0,\cdots,7$. Then we only need to
apply ordinary nonlinear optimization techniques to analyze and to solve
(12-13). However, we may apply the approach in [20] to simplify the language.

In the following, as in [20], we say that the real valued functions
$f,f_{d},h_{j},(h_{j})_{d}$ are continuously differentiable, if they are
continuously differentiable with respect to real vector variables
$\mathbf{x}_{k}$ for $k=0,\cdots,7$. We now take $f$ as an example, and the
same argument applies to $f_{d},h_{j}$ and $(h_{j})_{d}$. Let
$\mathbb{H}={\mathbb{Q}}^{n}\times{\mathbb{Q}}^{n}$. Denote the partial
derivatives and gradient of $f$ at $(\mathbf{x},\mathbf{x}_{d})$ as

${\partial\over\partial\mathbf{x}}f(\mathbf{x},\mathbf{x}_{d})={\partial\over\partial\mathbf{x}_{0}}f(\mathbf{x},\mathbf{x}_{d})+{\partial\over\partial\mathbf{x}_{1}}f(\mathbf{x},\mathbf{x}_{d})\mathbf{i}+{\partial\over\partial\mathbf{x}_{2}}f(\mathbf{x},\mathbf{x}_{d})\mathbf{j}+{\partial\over\partial\mathbf{x}_{3}}f(\mathbf{x},\mathbf{x}_{d})\mathbf{k},$
${\partial\over\partial\mathbf{x}_{d}}f(\mathbf{x},\mathbf{x}_{d})={\partial\over\partial\mathbf{x}_{4}}f(\mathbf{x},\mathbf{x}_{d})+{\partial\over\partial\mathbf{x}_{5}}f(\mathbf{x},\mathbf{x}_{d})\mathbf{i}+{\partial\over\partial\mathbf{x}_{6}}f(\mathbf{x},\mathbf{x}_{d})\mathbf{j}+{\partial\over\partial\mathbf{x}_{7}}f(\mathbf{x},\mathbf{x}_{d})\mathbf{k},$
$\nabla
f(\mathbf{x},\mathbf{x}_{d})=\left({\partial\over\partial\mathbf{x}}f(\mathbf{x},\mathbf{x}_{d}),{\partial\over\partial\mathbf{x}_{d}}f(\mathbf{x},\mathbf{x}_{d})\right).$

Then $\nabla f(\mathbf{x},\mathbf{x}_{d})\in\mathbb{H}$. With the R-linear
independence concept and Theorem 4.3 of [20], we may have the first optimality
conditions for (12-13). This is essentially the linear independence constraint
qualification in real nonlinear optimization.

##  5 Equality Constrained Standard Dual Quaternion Optimization

If $\hat{f}$ and $\hat{h}_{j}$ for $j=1,\cdots,m$ are all standard dual
quaternion functions, then (EQQOPI) has the form

$\min\\{f(\mathbf{x}):h_{j}(\mathbf{x})=0,(h_{j})_{d}(\mathbf{x},\mathbf{x}_{d})=0,\
j=1,\cdots,m\\},$ (14)

and (EQQOPII) has the form

$\min\\{f_{d}(\mathbf{x},\mathbf{x}_{d}):f(\mathbf{x})=L^{I},\
h_{j}(\mathbf{x})=0,\ (h_{j})_{d}(\mathbf{x},\mathbf{x}_{d})=0,\
j=1,\cdots,m\\}.$ (15)

Then we may have some better results. See the theorem and its corollary below.
We call such a dual quaternion optimization problem a standard dual quaternion
optimization problem.

######  Theorem 5.1.

Suppose that $f$, $h_{j}$ for $j=1,\cdots,m$ are continuously differentiable.

Assume $(\mathbf{x}^{opt},\mathbf{x}_{d}^{opt})$ is an optimal solution of
(14). If $\left\\{\nabla h_{j}(\mathbf{x}^{opt}):j=1,\cdots,m\right\\}$ is
R-linearly independent in the sense of [20], and
$\left\\{{\partial\over\partial\mathbf{x}_{d}}(h_{j})_{d}(\mathbf{x}^{opt},\mathbf{x}_{d}^{opt}):j=1,\cdots,m\right\\}$
is R-linearly independent in the sense of [20], then there are real Lagrangian
multipliers $\lambda_{j}$ for $j=1,\cdots,m$, such that

$\nabla f(\mathbf{x}^{opt})+\sum_{j=1}^{m}\lambda_{j}\nabla
h_{j}(\mathbf{x}^{opt})=0,$ (16)

and

$h_{j}(\mathbf{x}^{opt})=0,\
(h_{j})_{d}(\mathbf{x}^{opt},\mathbf{x}_{d}^{opt})=0,\ {\rm for}\
j=1,\cdots,m.$ (17)

Assume $(\mathbf{x}^{opt},\mathbf{x}_{d}^{opt})$ is an optimal solution of
(15). If $\left\\{\nabla f(\mathbf{x}^{opt}),\nabla
h_{j}(\mathbf{x}^{opt}):j=1,\cdots,m\right\\}$ is R-linearly independent in
the sense of [20], and
$\left\\{{\partial\over\partial\mathbf{x}_{d}}(h_{j})_{d}(\mathbf{x}^{opt},\mathbf{x}_{d}^{opt}):j=1,\cdots,m\right\\}$
is R-linearly independent in the sense of [20], then there are real Lagrangian
multipliers $\sigma$ and $\lambda_{j}$ for $j=1,\cdots,m$, such that

$\nabla f_{d}(\mathbf{x}^{opt},\mathbf{x}_{d}^{opt})+\sigma\nabla
f(\mathbf{x}^{opt})+\sum_{j=1}^{m}\lambda_{j}\nabla
h_{j}(\mathbf{x}^{opt})=0,$ (18)

and

$f(\mathbf{x}^{opt})=L^{I},\ h_{j}(\mathbf{x}^{opt})=0,\
(h_{j})_{d}(\mathbf{x}^{opt},\mathbf{x}_{d}^{opt})=0,\ {\rm for}\
j=1,\cdots,m.$ (19)

###### Proof.

By Theorem 4.3 of [20] and the property of block triangular matrices, if
$\left\\{\nabla h_{j}(\mathbf{x}^{opt}):j=1,\cdots,m\right\\}$ is R-linearly
independent in the sense of [20], and
$\left\\{{\partial\over\partial\mathbf{x}_{d}}(h_{j})_{d}(\mathbf{x}^{opt},\mathbf{x}_{d}^{opt}):j=1,\cdots,m\right\\}$
is R-linearly independent in the sense of [20], then there are real Lagrangian
multipliers $\lambda_{j}$ and $\mu_{j}$ for $j=1,\cdots,m$, such that

$\nabla f(\mathbf{x}^{opt})+\sum_{j=1}^{m}\left[\lambda_{j}\nabla
h_{j}(\mathbf{x}^{opt})+\mu_{j}{\partial\over\partial\mathbf{x}}(h_{j})_{d}(\mathbf{x}^{opt},\mathbf{x}_{d}^{opt})\right]=0,$
(20)
$\sum_{j=1}^{m}\mu_{j}{\partial\over\partial\mathbf{x}_{d}}(h_{j})_{d}(\mathbf{x}^{opt},\mathbf{x}_{d}^{opt})=0,$
(21)

and

$h_{j}(\mathbf{x}^{opt})=0,\
(h_{j})_{d}(\mathbf{x}^{opt},\mathbf{x}_{d}^{opt})=0,\ {\rm for}\
j=1,\cdots,m.$ (22)

Since
$\left\\{{\partial\over\partial\mathbf{x}_{d}}(h_{j})_{d}(\mathbf{x}^{opt},\mathbf{x}_{d}^{opt}):j=1,\cdots,m\right\\}$
is R-linearly independent in the sense of [20], by (21), all $\mu_{j}=0$. We
have the first result of this theorem. The second result can be derived
similarly. ∎

Actually, Theorem 5.1 may be stated and proved in the Language of real
optimization. That way is somewhat tedious.

######  Corollary 5.2.

Let $m=1$ and denote $h\equiv h_{1}$. Suppose that $f$ and $h$ are
continuously differentiable.

Assume $(\mathbf{x}^{opt},\mathbf{x}_{d}^{opt})$ is an optimal solution of
(14). Assume that $\nabla h(\mathbf{x}^{opt})\not=\mathbf{0}$, and
${\partial\over\partial\mathbf{x}_{d}}h_{d}(\mathbf{x}^{opt},\mathbf{x}_{d}^{opt})\not=\mathbf{0}$.
Then there is a real Lagrangian multiplier $\lambda$, such that

$\nabla f(\mathbf{x}^{opt})+\lambda\nabla h(\mathbf{x}^{opt})=0,$ (23)

and

$h(\mathbf{x}^{opt})=0,\ h_{d}(\mathbf{x}^{opt},\mathbf{x}_{d}^{opt})=0.$ (24)

Assume $(\mathbf{x}^{opt},\mathbf{x}_{d}^{opt})$ is an optimal solution of
(15). If $\left\\{\nabla f(\mathbf{x}^{opt}),\nabla
h(\mathbf{x}^{opt})\right\\}$ is R-linearly independent in the sense of [20],
then there are real Lagrangian multipliers $\sigma$ and $\lambda$, such that

$\nabla f_{d}(\mathbf{x}^{opt},\mathbf{x}_{d}^{opt})+\sigma\nabla
f(\mathbf{x}^{opt})+\lambda\nabla h(\mathbf{x}^{opt})=0$ (25)

and

$f(\mathbf{x}^{opt})=L^{I},\ h(\mathbf{x}^{opt})=0,\
h_{d}(\mathbf{x}^{opt},\mathbf{x}_{d}^{opt})=0.$ (26)

In practice, we may only find an approximal value of $L^{I}$, an error
analysis is necessary in the application of this approach.

This approach cannot be extended to the inequality constraint case naturally,
the optimal index set there is dependent upon the concrete solution.

##  6 Hand-Eye Calibration

The hand-eye calibration problem (the sensor-actuator problem) is an important
application problem in robot research. They can be solved by using dual
quaternions [10, 15, 16]. In this section, we formulate it as an equality
constrained dual quaternion optimization problem, and show that this problem
is a standard dual quaternion optimization problem.

In 1989, Shiu and Ahmad [21] and Tsai and Lenz [22] formulated the hand-eye
calibration problem as a matrix equation

$AX=XB,$ (27)

where $X$ is the transformation matrix from the camera (eye) to the gripper
(hand), $A=A_{2}A_{1}^{-1}$ and $B=B_{2}^{-1}B_{1}$, $A_{i}$ is the
transformation matrix from the camera to the world coordinate system, and
$B_{i}$ is the transformation matrix from the robot base to the gripper [10].
They are all $4\times 4$ matrices, $A_{i}$ and $B_{i}$ can be measured, while
$X$ is the variable matrix. Assume that there are $n+1$ measurements (poses)
$A_{i}$ and $B_{i}$ for $i=1,\cdots,n+1$. Denote $A^{(i)}=A_{i+1}A_{i}^{-1}$
and $B^{(i)}=B_{i+1}^{-1}B_{i}$, for $i=1,\cdots,n$. Then the problem is to
find the best solution $X$ of

$A^{(i)}X=XB^{(i)},$ (28)

for $i=1,\cdots,n$. In 1999, Daniilidis [10] proposed to use dual quaternions
to solve this problem. See Figure 1 for the geometry of this model.

Figure 1: $AX=XB$ Hand-Eye Calibration Model.

In 1994, Zhuang, Roth and Sudhaker [26] generalized (27) to

$AX=YB,$ (29)

where $Y$ is the transformation matrix from the world coordinate system to the
robot base. Assume that we have $n$ measurements (poses) for $i=1,\cdots,n$.
Then the problem is to find the best solution $X$ and $Y$ of

$A_{i}X=YB_{i},$ (30)

for $i=1,\cdots,n$. In 2010, Li, Wang and Wu [15] proposed to use dual
quaternions to solve this problem. Also see [16]. See Figure 2 for the
geometry of this model.

Figure 2: $AX=YB$ Hand-Eye Calibration Model.

As the 3D movement of a rigid body may also be expressed by unit dual
quaternions, we may rewrite (28) and (30) as

$\hat{a}^{(i)}\hat{x}=\hat{x}\hat{b}^{(i)},$ (31)

for $i=1,\cdots,n$, and

$\hat{a}_{i}\hat{x}=\hat{y}\hat{b}_{i},$ (32)

for $i=1,\cdots,n$, where $\hat{a}^{(i)}$, $\hat{b}^{(i)}$, $\hat{a}_{i}$,
$\hat{b}_{i}$, $\hat{x}$ and $\hat{y}$ are corresponding unit dual
quaternions. We may solve (31) and (32) by the following minimization
problems:

$\min\left\\{\sum_{i=1}^{n}|\hat{a}^{(i)}\hat{x}-\hat{x}\hat{b}^{(i)}|:|\hat{x}|^{2}=1\right\\},$
(33)

and

$\min\left\\{\sum_{i=1}^{n}|\hat{a}_{i}\hat{x}-\hat{y}\hat{b}_{i}|:|\hat{x}|^{2}=1,|\hat{y}|^{2}=1\right\\}.$
(34)

Then (33) and (34) are equality constrained dual quaternion optimization
problems. By the analysis in Section 3, all the functions involved are
standard dual quaternion optimization problems. Thus, they are equality
constrained standard dual quaternion optimization problems.

Note that by (5), the magnitudes in the objective functions of (33) and (34)
cannot be replaced by their squares as for a dual quaternion number $\hat{q}$,
$|\hat{q}|=0$ and $|\hat{q}|^{2}=0$ are not equivalent.

Consider (33). We have
$\hat{f}(\hat{x})=\sum_{i=1}^{n}|\hat{a}^{(i)}\hat{x}-\hat{x}\hat{b}^{(i)}|$
and $\hat{h}(\hat{x})=|\hat{x}|^{2}-1$. Then EQSDQO has the form

$\min\\{\hat{f}(\hat{x}):\hat{h}(\hat{x})=0\\}.$ (35)

By (9) and (10),

$h(x,x_{d})=h(x)=|x|^{2}-1,$

and

$h_{d}(x,x_{d})=xx_{d}^{*}+x_{d}x^{*}.$

Similarly, we have

$f(x,x_{d})=f(x)=\sum_{i=1}^{n}|a^{(i)}x-xb^{(i)}|,$

and $f_{d}(x,x_{d})$ can be computed by (5). Then, (EQQOPI) has the form

$\min\\{f(x):h(x)=0,h_{d}(x,x_{d})=0,\\},$ (36)

and (EQQOPII) has the form

$\min\\{f_{d}(x,x_{d}):f(x)=L^{I},\ h(x)=0,\ h_{d}(x,x_{d})=0\\}.$ (37)

We may apply Theorems 4.1 and 5.1 as well as Corollary 5.2 now. In particular,
we have $\nabla h(x^{opt})=2x^{opt}\not=0$ as $|x^{opt}|^{2}=1$. The problem
(34) can be treated similarly.

##  7 Pose Graph Optimization

The simultaneous localization and mapping (SLAM) problem is a very hot topic
in robotic research [5]. The applications of SLAM include not only field
robots, but also underwater navigation [25], unmanned aerial vehicle (UAV)
[3], planetary exploration rover [24]. The graph-based approach poses the SLAM
problem to the pose graph optimization problem. In 2016, Cheng, Kim, Jiang and
Che [9] studied dual quaternion-based graph SLAM. In this section, we show
that the dual quaternion pose graph optimization problem is a standard dual
quaternion optimization problem.

Pose graph optimization estimates $n$ robot poses from $m$ relative pose
measurements. The dual quaternion pose graph problem can be visualized as a
directed graph $G=(V,E)$ [6], where each vertex $i\in V$ corresponds to a
robot pose $\hat{x}_{i}$ for $i=1,\cdots,n$, and each directed edge (arc)
$(i,j)\in E$ corresponds to a relative measurement $\hat{q}_{ij}$, a unit dual
quaternion. Thus, $\hat{x}_{i}$ for $i=1,\cdots,n$, are $n$ unknown unit dual
quaternion variables. The cardinality of the edge set $|E|=m$ as we have $m$
relative measurements. By the pose relation, the relative pose $\hat{y}_{ij}$
should be

$\hat{y}_{ij}=\hat{x}_{i}^{*}\hat{x}_{j}.$

Then the error at the edge $(i,j)$ is

$\hat{e}_{ij}=\hat{q}_{ij}-\hat{x}_{i}^{*}\hat{x}_{j}.$

We may regard $\hat{\bf e}=(\hat{e}_{ij}:(i,j)\in E)$ as an $m$-dimensional
dual quaternion vector, i.e., $\hat{\bf e}\in{\hat{\mathbb{Q}}}^{m}$. Then we
have the following minimization model for this pose graph optimization
problem:

$\min\left\\{\|\hat{\bf e}\|_{2}:|\hat{x}_{i}|^{2}=1,i=1,\cdots,n\right\\},$
(38)

where $\|\hat{\bf e}\|_{2}$ is defined by (9) and (10). By the discussion in
Section 3, this is an equality constrained standard dual quaternion
optimization problem.

We will study the algorithm aspect of applying standard dual quaternion
optimization approach to the hand-eye calibration problem and the SLAM problem
in further coming papers.

Acknowledgment I am thankful to Chen Ling, Ziyan Luo and Zhongming Chen for
the discussion on standard dual quaternion optimization, to Wei Li for the
discussion on hand-eye calibration, to Jiantong Cheng for the discussion on
SLAM, to Guyan Ni for introducing Jiantong Cheng to me, and to Chen Ouyang and
Jinjie Liu for Figures 1 and 2.

## References

  * [1] H. Abaunza, P. Castillo, A. Victorino and R. Lozano, “Dual quaternion modeling and control of a quad-rotor aerial manipulator”, Journal of Intelligent & Robotic Systems 88 (2017) 267-283.
  * [2] G. Brambley and J. Kim, “Unit dual quaternion-based pose optimization for visual runway observations”, IET Cyber-Systems and Robotics 2 (2020) 181-189.
  * [3] M. Bryson and S. Sukkarieh, “Building a robust implementation of bearing - only inertial SLAM for a UAV”, Journal of Field Robotics 24 (2007) 113-143.
  * [4] S. Bultmann, K. Li and U.D. Hanebeck, “Stereo visual SLAM based on unscented dual quaternion filtering”, 2019 22th International Conference on Information Fusion (FUSION) (2019) 1-8.
  * [5] C. Cadena, L. Carlone, H. Carrillo, Y. Latif, D. Scaramuzza, J. Neira, I. Reid and J.J. Leonard, “Past, present and future of simultaneous localization and mapping: Toward the robust-perception age”, IEEE Transactions on Robotics 32 (2016) 1309-1332.
  * [6] L. Carlone, R. Tron, K. Daniilidis and F. Dellaert, “Initialization techniques for 3D SLAM: A survey on Rotation and its use in pose graph optimization”, IEEE International Conference on Robotics and Automation (ICRA) (2015) 4597-4604.
  * [7] Y. Chen, L. Qi and X. Zhang, “Color image completion using a low-rank quaternion matrix approximation”, Pacific Journal of Optimization 18 (2022) 55-75.
  * [8] Y. Chen, X. Xiao and Y. Zhou, “Low-rank quaternion approximation for color image processing”, IEEE Transactions on Image Processing 29 (2019) 1057-7149.
  * [9] J. Cheng, J. Kim, Z. Jiang and W. Che, “Dual quaternion-based graph SLAM”, Robotics and Autonomous Systems 77 (2016) 15-24.
  * [10] K. Daniilidis, “Hand-eye calibration using dual quaternions”, The International Journal of Robotics Research 18 (1999) 286-298.
  * [11] J. Flamant, S. Miron and D. Brie, “A general framework for constrained convex quaternion optimization”, IEEE Transactions on Signal Processing, doi:10.1109/TSP.2021.3137746.
  * [12] Z. Jia, M.K. Ng and G. Song, “Robust quaternion matrix completion with applications to image inpaiting”, Numerical Linear Algebra with Applications 26 (2019) e2245.
  * [13] B. Kenright, “A beginners guide to dual-quaternions”, 20th International Conference in Central Europe on Computer Graphics, Visualization and Computer Vision, Plzen, Czech, 2012.
  * [14] G. Leclercq, Ph. Lefévre and G. Blohm, “3D kinematics using dual quaternions: Theory and applications in neuroscience”, Frontiers in Behavioral Neuroscience 7 (2013) Article 7, 1-25.
  * [15] A. Li, L. Wang and D. Wu, “Simultaneous robot-world and hand-eye calibration using dual-quaternions and Kronecker product”, International Journal of the Physical Sciences 5 (2010) 1530-1536.
  * [16] W. Li, N. Lv, M. Dong and X. Lou, “Simultaneous robot-world/hand-eye calibration using dual quaternion”, Robot 40 (in Chinese) (2018) 301-308.
  * [17] X. Peng, X. Shi and Y. Gong, “Dual-quaternion-based modelling and control for motion tracking of a tumbling target”, Journal of Systems Engineering and Electronic 30 (2019) 985-994.
  * [18] L. Qi, C. Ling and H. Yan, “Dual quaternions and dual quaternion vectors”, Communications on Applied Mathematics and Computation DOI: 10.1007/s42967-022-00189-y.
  * [19] L. Qi and Z. Luo, “Eigenvalues and singular values of dual quaternion matrices”, to appear in: Pacific Journal of Optimization.
  * [20] L. Qi, Z. Luo, Q. Wang and X. Zhang, “Quaternion matrix optimization: Motivation and analysis”, Journal of Optimization Theory and Applications, DOI: 10.1007/s10957-021-01906-y.
  * [21] Y. Shiu and S. Ahmad, “Calibration of wrist-mounted robotic sensors by solving homogeneous transform equation of the form $AX=XB$”, IEEE Transactions on Robotics and Automation 5 (1989) 16-27.
  * [22] R. Tsai and R. Lenz, “A new technique for fully autonomous and efficient 3D robotic hand/eye calibration”, IEEE Transactions on Robotics and Automation 5 (1989) 345-358.
  * [23] X. Wang, C. Yu and Z. Lin, “A dual quaternion solution to attitude and position control for rigid body coordination”, IEEE Transactions on Robotics 28 (2012) 1162-1170.
  * [24] E. Wei, S. Jin and Q. Zhang, “Autonomous navigation of mars probe using X-ray pulsars: Modeling and results”, Advances in Space Research 51 (2013) 849-857.
  * [25] S.B. Williams, P. Newman and G. Dissanayake, “Autonomous underwater simultaneous localization and map building”, IEEE International Conference on Robotics and Automation (ICRA), (2000) 1793-1798.
  * [26] H.Q. Zhuang, Z. Roth and R. Sudhakar, “Simultaneous robot-world and tool/flange calibration by solving homogeneous transformation of the form $AX=YB$”, IEEE Transactions on Robotics and Automation 10 (1994) 549-554.


# Understanding Police Force Resource Allocation using Adversarial Optimal
Transport with Incomplete Information

Yinan Hu Tandon School of Engineering, New York University, Brooklyn, NY,
11201 Juntao Chen Graduate School of Art and Sciences, Fordham University,
Bronx, NY, 10458 Quanyan Zhu the same as 1st author

Adversarial optimal transport has been proven useful as a mathematical
formulation to model resource allocation problems to maximize the efficiency
of transportation with an adversary, who modifies the data. It is often the
case, however, that only the adversary knows which nodes are malicious and
which are not. In this paper we formulate the problem of seeking adversarial
optimal transport into Bayesian games. We construct the concept of Bayesian
equilibrium and design a distributed algorithm that achieve those equilibria,
making our model applicable to large-scale networks.
Keywords: game theory, crime control, Markov games

JEL Codes: key1, key2, key3

##  1 Introduction

The problem of resource matching and allocation has a wide application in many
areas such as wireless networks Georgiadis et al. (2006), data centers Ferris
and Riveros (2013). The ultimate goal of resource allocation is to find a plan
to distribute resources to the target so as to maximize the resource
efficiency. Optimal transport Kantorovich (1942) is one of the centralized
schemes to find plans for resource matching and allocations Galichon (2018);
Villani (2008).

In a standard network transport scheme, a dispatcher designs the resource
allocation plan that maximizes the accumulated utility function for all
dispatchers Hughes and Chen (2021b) under the assumption that the data
received from the dispatchers are truthful. This assumption may not always
hold when some dispatchers present selfish Felegyhazi et al. (2006); Raya et
al. (2008) or even malicious behaviors when submitting their data Buttyan and
Hubaux (2007). Malicious nodes modify the data sent to the dispatcher to
disrupt the fairness or the efficiency of the resulting transport plan.
Previously there are worksHughes and Chen (2021a) that develop a game-
theoretic formulation to characterize the interactions between the adversary
and the dispatcher in resource allocation processes under the assumption that
all the information in the game is common knowledge. That is, the dispatcher
knows before the game all the possible modifications that the adversary may
adopt upon the data, while the adversary is aware of all possible transport
plans and they both know what their opponent know, and so on. However, a
general game between an adversary and a defender is one with incomplete
information (see chapter 9 of Maschler et al. (2013)) where the defender is
not aware of the adversary’s the benignity or malice. In a network transport
scenario, where there are multiple target nodes that receive resources, the
dispatcher may not know which target nodes are benign and truthfully report
their data and which ones are malicious with the intent to modify the data
along the links connected to them. A malicious target node may alter the data
along the link connected to one particular source node or along the ones
connected to several source nodes.

To this end, we propose a Bayesian game formulation to characterize the
asymmetric information and adopt the concept of ‘types’, raised in Harsanyi’s
model Harsanyi (1967, 1968) to capture the adversary’s attack mode. Each of
the adversary’s type represent the links on which the data are compromised.
The Bayesian game formulation allows the dispatcher to associate a common
distribution over all possible types, so his goal is to maximize the expected
utility on the distribution over the type space. The adversary aims at
minimizing the aggregated utility of all nodes over his possible modifications
given his type. The proposed scenario concerns a variety of attacks in
resilient resource allocations such as jamming attacks Garnaev and Trappe
(2014) and network topology attacks Shao et al. (2020).

The planner’s utility functions are usually monotonic with respect to the
‘weights’ along the links in the network Srikant and Ying (2013b). The amount
of resources transported along links with large weight coefficients tend to
dominate in the network and thus causes non-smooth transport plans. Inspired
by Cuturi (2013), we add an extra term, named a regularizer Ferradans et al.
(2014b), in the utility function with a coefficient that leverages between the
solving traditional optimal transport problem and preserving the smoothness of
the transport plan. In particular, we introduce the entropic regularizer
because the corresponding regularized optimal transport problems can be solved
in distributed algorithms and are thus scalable to large networks Genevay et
al. (2016).

A centralized planning of the transport within the network requires full
knowledge of the data on all links between source nodes and target nodes. As a
result, the computational complexity for an optimal transport with a
dispatcher grows exponentially with the number of participants. Our goal is to
design a decentralized algorithm to obtain a resilient transport plan so that
such algorithms can be applied in large-scale networks. We obtain the best
responses for source nodes through asynchronous dual pricing (ADP) algorithm
Huang et al. (2006) so that every source node only needs to solve its own
Bayesian game against target nodes it connects to.

The contributions of our paper are threefold: 1) we extend the zero-sum game-
theoretic framework raised by Hughes and Chen (2021a) into a Bayesian game
formulation to capture the adversary’s private information. We also extend the
previous static game into a multistage dynamic game formulation to capture the
adversary’s action in altering the belief updates; 2) We introduce an entropic
regularizer in the dispatchers’ utility functions, smoothing the resulting
transport plans; 3) we develop a distributed algorithm to solve our Bayesian
game model so that our formulation is applicable to large-scale networks,
especially those with ephemeral nodes and thus overcomes the curse of
dimensionality. 4) we apply our dynamic adversarial optimal transport scheme
to help developing an optimal plan of allocating police resources to combat
opportunistic criminals.

The rest of the paper is organized as follows. In section 2 we first consider
the problem of regularized optimal transport over networks without adversary,
after which we introduce the adversary that compromises the utility function
according to his type, which is unknown to the dispatcher. In section 4 we
come up with a distributed algorithm that solves the pairs of the best
response of strategies. We demonstrate our simple example in section 6.
Finally we conclude our paper in section 7.

###  1.1 Related work

Learning adaptive learning has been a useful tool to study the crime
neighborhood patrolling.

The topic of police resource allocation for patrolling has been attracting
attention for a long time Curtin et al. (2010). Many methods have been
implemented to characterize and study how to design optimal patrol plan to
counteract both opportunistic criminals and strategic adversarys. Different
optimal patrol plans have different goals. Authors in minimizing the response
time or maximizing the patrol coverage The proposed framework aims a

Authors in Zhang et al. (2016a) formulate the patrolling model with dynamical
Bayesian networks (DBN) and computes the optimal parameters by learning the
criminal’s behaviors. Author in Hespanha et al. (2000b) proposes a pursuit-
evasion game framework to capture the conflicting relationship between
patrolling police and criminals. However the goal of criminals in such a game
model is only to evade capture and does not include any reward function from
committing crimes. Authors in Jiang et al. (2013) adopts a Stackelberg
security games (SSG) with randomized strategies model to compute an optimal
patrolling schedule against passive criminals while incorporating uncertainty
in the execution of the patrolling plan. Instead of combatting opportunistic,
passive crimes, the proposed framework designs a patrolling plan against
strategic, organized criminals, who knows the police’s patrolling allocation
plan by methodical investigations and undermines the plan strategically to
counteract the police.

The formulation of optimal transport with adversaries has been serving as a
useful tool to be implemented to resource allocation problems..

Efficient methods in computing the problems over networks include…..Srikant
and Ying (2013a)

Asychronous dual pricing…

The rest of the work is organized as follows.

##  2 Formulation of source allocation without adversaries

Although different experts may have different orders, but in general there is
an agreement that the allocation of police forces serve the following
purposes: deterrence of crimes; apprehension of criminals; providing on-demand
non-crime police services to the neighborhoods; building trust and confidence
among neighborhoods and a sense of security.

The evaluation of the police work may include the change in crime rate;
reduction in police response time; an increase in arrests, etc.

In this section, we propose a formulation optimal transport with regularizers
over a bipartite network.

###  2.1 Discrete regularized optimal transport over networks

In this paper we study the resource allocation process over a bipartite
network $G$ that has $N$ source nodes and $M$ target nodes. We denote the sets
of source nodes and target nodes as
$J=\\{j_{1},j_{2},\dots,j_{N}\\},Q=\\{q_{1},q_{2},\dots,q_{M}\\}$,
respectively. Let $\mathcal{E}$ be the set of possible edges or links, each of
which connects a source node $j\in J$ and a target node $q\in Q$. We denote
$\\{j,q\\}$ as a generic link connecting a source node $j$ and a target node
$q$. Let $Q_{j}$ be the set of links that begin with the source node $j$.
Correspondingly, let $J_{q}$ be the set of links that end with the target node
$q$.

To characterize the connectivity within the bipartite network $G$, we
introduce a 0-1 matrix $B$ with its entries $(b_{jq})_{\begin{subarray}{c}j\in
J\\\ q\in Q\end{subarray}}$ defined as

$b_{jq}=\begin{cases}1&r=\\{j,q\\},\;\text{if}\;\\{j,q\\}\in\mathcal{E},\\\
0&\text{otherwise}.\end{cases}$ (1)

We characterize a generic dispatcher’s action $x\in
X=\mathbb{R}^{|J|\times|Q|}$ as a matrix representation of the rates of
transportation of resource from source nodes to target nodes, or the transport
plan. For convenience, we denote a generic entry, $x_{jq}$ to represent the
transport rate along the edge $\\{j,q\\}\in\mathcal{E}$. Notice the transport
rate should always be non-negative: $x\geqslant 0$ (here ‘$\geqslant$’ or
‘$\leqslant$’ stands for element-wise comparisons) For a typical source node
$j$, we associate its finite capacity $c_{j}>0$ to be limitation of the amount
of the resources distributed from the source node $j$.

####  2.1.1 The dispatcher’s utility function

We now formulate the planner’s utility function in a two-step way: we first
formulate a series of utility functions $\\{u_{jq}\\}_{\begin{subarray}{c}j\in
J\\\ q\in Q\end{subarray}}:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ for the
planner on the link $r$ (which connects the source node $j$ and the target
node $q$); then we extend $u_{jq}$ to a utility function
$U:\mathbb{R}^{|J|\times|Q|}\rightarrow\mathbb{R}_{+}$ defined on all
transport plans over the network. We have the following requirements for a
utility function in terms of the transporting rate $x_{jq}$.

######  Assumption 1 (Monoticity, convexity, and smoothness)

For all $j\in J,q\in Q,\;\\{j,q\\}\in\mathcal{E}$, the utility function
$u_{jq}:\mathbb{R}\rightarrow\mathbb{R}$ is monotone increasing, concave,
continuously differentiable in terms of $x_{jq}\in\mathbb{R}$.

We consider $u_{jq}$ as a sum of two terms: $u_{jq}=v_{jq}+\lambda h_{jq}$,
where $v_{jq}:\mathbb{R}\rightarrow\mathbb{R}$ could be interpreted as a
‘revenue function’ and $h_{jq}$ is named as a regularizer and $\lambda$ is a
trade-off parameter that leverages the portion of the two functions upon the
utility function. The first term $v_{jq}$ represents the ‘revenue’ that the
planner makes by transporting the resource at the rate $x_{jq}$. Possible
candidates for $v_{jq}$Zhang and Zhu (2019) include: linear functions, log-
concave functions, etc, depending on the speed of the increase of the revenue
with respect to the resources. We choose
$v_{jq}=m_{jq}x_{jq},\;\forall\\{j,q\\}\in\mathcal{E}$. The coefficient
$m_{jq}$ reflects the ‘perception’‘ of crime activities on the route
$\\{j,q\\}$ and larger coefficients lead to a hot spot upon the crime heat map
among the neighborhoods. The function
$h_{jq}:\mathbb{R}_{+}\rightarrow\mathbb{R}$ in (4) serves the role of a
regularizer Cuturi (2013) to prevent the planner to allocate all of the
resources to a specific target node. We impose the following assumption on the
regularizer function:

######  Assumption 2 (Fairness)

The fairness functions $h_{jq},\;\\{j,q\\}\in\mathcal{E}$ are monotonically
increasing, convex and differentiable.

Typical examples of regularizers include the Kullback–Leibler (KL)-divergence,
total variation (TV) functions Ferradans et al. (2014a), quadratic functions
Essid and Solomon (2018). In this paper, we adopt the entropic regularizer for
every link: for all $x_{jq}\in\mathbb{R}_{+}$,

${h_{jq}(x_{jq})}:=-\lambda{x_{jq}\log x_{jq}},\;\forall j\in J,\;q\in Q.$ (2)

Now we extend $u_{jq}$ to a utility function of transport rates over all
connections in the network, which is denoted as $U$. We assume $U$ is
separable and additive in terms of the transport rate in every node in the
network, that is, there is no explicit coupling between the transportation
rate on any other two links. We set the utility function in terms of transport
rate $x$ of the network as follows.

$U(x)=V(x)+\lambda H(x):=\sum_{\\{j,q\\}\in\mathcal{E}}{v_{jq}(x_{jq})+\lambda
h_{jq}(x_{jq})},$ (3)

where we assume $V,H$, are also both separable and additive in terms of the
links in the network
$H(x)=\sum_{\\{j,q\\}\in\mathcal{E}}{h_{jq}(x_{jq})},\;V(x)=\sum_{\\{j,q\\}\in\mathcal{E}}{v_{jq}(x_{jq})}$.
In an adversary-free scenario, the goal of the discrete optimal transport over
a network is to find an optimal transport plan $x^{*}$ that maximizes the
planner’s the utility function, with the constraint that the transportation of
resources does not overwhelm any target rate. We can state the problem of
seeking an optimal transport plan as an optimization problem Srikant and Ying
(2013a) as follows:

$\displaystyle\underset{x\geqslant 0}{\max}$
$\displaystyle\sum_{\\{j,q\\}\in\mathcal{E}}{u_{jq}(x_{jq})},$ (4) s.t.
$\displaystyle\;\;Bx\leqslant c.$

Referring to the form of objective function in (4), we write the formulation
of optimization problem that the dispatcher faces as follows:

$\displaystyle\underset{x\geqslant 0}{\max}$
$\displaystyle\;\sum_{\\{j,q\\}\in\mathcal{E}}{m_{jq}x_{jq}-\lambda x_{jq}\log
x_{jq}},$ (5) s.t. $\displaystyle\;\;B\;\text{vec}(x)\leqslant c,$

The resulting formulation is sometimes called regularized optimal transport
Cuturi (2013).

For the purpose of demonstration, we first depict the plot to illustrate the
role of ‘fairness’ term. Let the number of source nodes and target nodes to be
$|J|=2,|Q|=3$. Let the ‘perception’ matrix $m$ to be

$m=\begin{bmatrix}1&3&5\\\ 2&5&1\end{bmatrix}.$ (6)

Let the capacity (maximum source available for each node) to be
$c=\begin{bmatrix}4&3&5\end{bmatrix}^{T}$. Then as a bipartite network where
every source node is connected to every destination node, the connection
matrix $B\in\mathbb{R}^{|J|\times|Q|}$ can be written as

$B=\begin{bmatrix}1&1&1&0&0&0\\\ 0&0&0&1&1&1\end{bmatrix}.$ (7)

To solve the dispatcher’s optimization problem, we first list the Lagrangian

$\mathcal{L}_{j}(x,p)=\sum_{q}{m_{jq}x_{jq}-\lambda x_{jq}\log
x_{jq}}-p_{j}(\sum_{q}{x_{jq}}-c_{j}),$ (8)

where we notice that the variables $p_{j}$ are ‘dual pricing’ variables.
According to Srikant and Ying (2013a), we solve the optimal transport plan
through the following iterations

$\displaystyle x^{(k+1)}_{jq}$
$\displaystyle=\exp\left(\frac{m_{jq}-p^{(k)}_{j}}{\lambda}-1\right),\;\forall
j\in J,q\in Q,$ (9) $\displaystyle p^{(k+1)}_{j}$
$\displaystyle=p^{(k)}_{j}-\gamma(\sum_{q}{x_{jq}}-c_{j})^{+}_{p^{(k)}_{j}},j\in
J,$ (10)

where $\gamma$ refers to the step size and

$(x)^{+}_{y}=\begin{cases}x&y>0;\\\ max(x,0)&y=0\end{cases}$ (11)

##  3 Adversarial regularized optimal transport

We now introduce the concept of ‘adversarial organizer’, who manages the
allocation of opportunistic criminals among neighborhoods to jeopardize the
safety of those neighborhoods and thus the undermine the revenue generated by
police allocation. Different from Zhang et al. (2016a) and Curtin et al.
(2010) that study individual opportunistic criminals, we treat crime activity
as a group phenomenon since the structure of social networks of individuals
play an important role in characterizing his delinquent behaviors.

To characterize the asymmetric information structure and the unknown influence
on the data in seeking an optimal transport, we consider a Bayesian game
$\mathcal{G}$ between a dispatcher and an adversary.

##### The adversary’s information structure

The adversary has private information named his type Harsanyi (1967), which
can be constructed in the following way. We denote $\theta_{q}=\\{1,2\\}$ as
adversary’s type space on the destination node $j$ and a generic type on the
target node is denoted as $\theta_{q}$. We denote $\theta_{q}=1$ if the
criminals at the target node $j$ commits minor crimes and $\theta_{q}=2$ if
they commit major crimes. Minor crimes incurs smaller loss on dispatcher’s
revenue Then the overall adversary’s type $\theta$ can be obtained by
computing Cartesian product among all the ‘sub-types’ as follows:

$\displaystyle\theta=\underset{q\in
Q}{\bigotimes}\;\theta_{q},\;\Theta=\underset{q\in
Q}{\bigotimes}\;\Theta_{q},\;\theta\in\Theta,\;\theta_{q}\in\Theta_{q},\;\forall
j\in J.$ (12)

##### Beliefs/types

In the Bayesian game setting, the dispatcher does not know the type the target
nodes. He develops a common belief $(\mu(\theta))_{\theta\in\Theta}$ , over
the adversary’s type space (12) that reflects the conditions of the target
nodes. We assume the criminal actions taken at every target node $q$ has a
unique type $\theta_{q}\in\\{1,2\\}$, where $\theta_{q}=1$ means minor crime
offenders, while $\theta_{q}=2$ means major crime offenders. We therefore have

$\displaystyle\mu(\theta)=\prod_{q\in Q}{\mu(\theta_{q})}.$ (13)

##### Action spaces

Similar to the situation in section 2, a typical dispatcher’s action can be
represented by a transport plan $x\in\mathbb{R}^{|J|\times|Q|}$, with each row
$x_{j},\;j\in J$ being the transport rate vector from a specific source node
$j$ and $x_{jq}$ as the transport rate along the link $\\{j,q\\}$. We shall
assume that the transport rate is always non-negative and upper bounded
because of the limited capacity of resource that every node owes. Referring to
the constraints in (15), we construct the strategy space for the dispatcher as
$X$:

$\displaystyle X_{j}$
$\displaystyle:=\\{x_{j}\in\mathbb{R}^{M}\;\Big{|}x_{j}\geqslant
0,\;\sum_{q\in Q}{x_{jq}}\leqslant c_{j}\\},$ (14) $\displaystyle X$
$\displaystyle:=\underset{j\in J}{\bigotimes}X_{j}=\\{x\in\mathbb{R}^{N\times
M}\;\Big{|}\;x\geqslant 0,\;Bx\leqslant c\\},$ (15)

where $B$ is the connection matrix described in (1) and $c$ is the vector of
capacities of the source nodes.

A generic behavioral strategy of the adversary is denoted as
$\xi:\Theta\rightarrow\Xi(\theta),\;\xi(\theta)=(\xi_{q}(\theta_{q}))_{\begin{subarray}{c}j\in
J\\\ q\in Q\end{subarray}}$, whose entries represent the increase of crime
perception regarding the neighborhood $\\{j,q\\}$. Meanwhile, the magnitude of
a criminal’s influence upon the perception of the crime activity of a
neighborhood should be bounded.

$\xi_{q}(\theta_{q}):=\begin{cases}\left\\{\xi_{q}:\theta_{q}\rightarrow\mathbb{R}^{d}_{+}\Big{|}\;{\;\xi_{q}(\theta_{q})\leqslant\underline{n}_{q}},\;q\in
Q\right\\}&\theta_{q}=1,\\\
\left\\{\xi_{q}:\theta_{q}\rightarrow\mathbb{R}^{d}_{+}\Big{|}\;{\;\xi_{q}(\theta_{q})\leqslant\bar{n}_{q},\;q\in
Q}\right\\}&\theta_{q}=2,\end{cases}$ (16)

where $m_{jq}$ is a parameter in the utility functions representing
‘perception’ of crime rate upon every neighborhood. The adversary’s action
space from the network’s viewpoint can be written as a direct product of
action spaces of every target node. Thus we have

$\displaystyle\Xi(\theta)$
$\displaystyle=\bigotimes_{\\{j,q\\}\in\mathcal{E}}{\xi_{q}(\theta_{q})}=\bigotimes_{j\in
J}{\Xi_{j}(\theta_{q})}.$ (17)

There are ‘hot spots’ of crime regions, where gang bangers gather around and
organized crimes happen and there are also ‘green’ areas with very few crime
reports every year. Dispatching police resources into ‘hot spots’ incurs
higher revenue since ‘hotpot’ places require more police to help suppress
illegal activities.

##### Stage Utility functions

We now formulate the dispatcher and the adversary’s utitlity/cost functions
the game.

The dispatcher leverages two purposes: on one hand, he makes a ‘revenue’ by
allocating resources along different links; on the other hand, he also induces
‘fairness’ in the allocation process, making sure that the resources
transported along different links do not vary too much. He maximizes his
expectation of the sum of revenue and the fairness penalty over his belief
upon the types of the target nodes. We denote
${U}^{(p)}:X\times\Xi\rightarrow\mathbb{R}$ as the utility function for the
dispatcher given the strategy profile of $(x,\xi)\in X\times\Xi$ as

$U^{(p)}(x,\xi)=\sum_{\begin{subarray}{c}\\{j,q\\}\in\mathcal{E}_{j}\\\
\theta_{q}\in\\{1,2\\}\end{subarray}}{(m_{jq}+\theta_{q}\xi_{q}(\theta_{q}))\mu(\theta_{q})x_{jq}}-\lambda
x_{jq}\log x_{jq}.$ (18)

The adversary’s utility function
$U^{(a)}:X\times\Xi\times\Theta\rightarrow\mathbb{R}$ not only depends on the
actions of the adversary and of the dispatcher, but also depends on the
adversary’s private type, which is a Cartesian product of the ‘sub-types’‘ of
every target node. According to Becker (1968), from the economist’s viewpoint
the adversary’s utility function
$U^{(a)}:\Xi\times\Theta\rightarrow\mathbb{R}$ consists of also two terms: the
‘revenue’ minus the ‘expectation of punishment’. The ‘revenue’ stands for the
‘benefits’, either financial or social, that criminals directly obtained from
committing crimes. We hereby assumes the revenue is negation of the police’s
revenue obtained from allocation polices forces, that is, the first term in
(28).

The terms of ‘expected punishment’ is a product of the probability of getting
caught and the cost of punishment (either fine, probation, imprisonment or
even execution). We assumes criminals who commit major crimes will benefit a
lot more than the criminals who commit minors, but also felony criminals run a
higher risk of getting caught as well as a receiving serious punishment.
Inspired by Ehrlich (1973), the criminals’ probability of getting caught on
every neighborhood relies positively on the amount of police resources and
negatively on the intensity of committing the crime. So we can express the
criminal’s expected penalty from committing criminal action $\xi_{q}$ based on
the type $\theta_{q}\in\Theta_{q}$ as follows:

$\nu_{jq}(\xi_{q}(\theta_{q}))=c_{jq}\xi_{q}^{-\beta_{2}}(\theta_{q})x^{\beta_{1}}_{jq},$
(19)

where $\beta_{1},\beta_{2}\in[0,1]$ are coefficients. For every type of the
adversary $\theta\in\Theta$, we thus formulate the adversarial organizer’s
utility function $U^{(a)}$ as

$U^{(a)}(x,\xi(\theta))={\sum_{\\{j,q\\}\in\mathcal{E}_{j}}{c_{jq}\xi_{q}^{-\beta_{2}}(\theta_{q})x^{\beta_{1}}_{jq}+(m_{jq}+\theta_{q}\xi_{q}(\theta_{q}))x_{jq}}},$
(20)

To summarize and adopt the theories of potential games Monderer and Shapley
(1996) , the adversary faces the following optimization problem:

$\displaystyle\underset{\begin{subarray}{c}\xi_{1}\in\Xi(1)\\\
\xi_{2}\in\Xi(2)\end{subarray}}{\min}$
$\displaystyle\;{\sum_{\\{j,q\\}\in\mathcal{E}_{j}}{c_{jq}(\xi^{-\beta_{2}}_{jq}(2)+\xi^{-\beta_{2}}_{jq}(1))x^{\beta_{1}}_{jq}+(m_{jq}+(\xi_{q}(1)+2\xi_{q}(2))x_{jq}}},$
(21) $\displaystyle\text{s.t}.$
$\displaystyle\;\;{\xi_{q}(1)}\leqslant\underline{n}_{q},\;q\in Q,$
$\displaystyle\;{\xi_{q}(2)}\leqslant\bar{n}_{q},\;\;q\in Q$

We can now define the Bayesian game to characterize the adversarial optimal
transport as follows:

######  Definition 1 (Bayesian game of adversarial transport)

We define the Bayesian game of adversarial transport $\mathcal{G}$ as a tuple

$\mathcal{G}=\langle\Xi,X,U^{(a)},U^{(p)},\Theta,\mu\rangle,$ (22)

where $\Xi,X$ are action spaces of the adversary and the dispatcher
characterized in (15) and (17), and $U^{(a)},U^{(p)}$ are utility functions of
the dispatcher and the adversary defined in (28) and (20) respectively. The
adversary has a non-singleton type space denoted as $\Theta$ associated with a
common belief $\mu$.

##  4 Analysis of Static Game

In section 2 we come up with formulations of a static Bayesian game and a
dynamic Bayesian game to characterize the optimal transport under adversarial
setting. We now solve the Bayesian equilibrium concept raised in section.

###  4.1 Existence of a solution

####  4.1.1 Bayesian equilibrium

The dispatcher aims at maximizing his expected regularized utility function,
while the adversary attempts to maximize his utility function. Thus we come up
with the definition of the Bayesian equilibrium for regularized adversarial
optimal transport problem as follows.

######  Definition 2 (Bayesian equilibrium for adversarial OT)

A Bayesian equilibrium for the game $\mathcal{G}$ (these are parameters of the
game) of adversarial transport is a strategy profile
$(x^{*},\\{\xi^{*}(\theta)\\}_{\theta\in\Theta}),$ such that neither player
has incentive to unilaterally deviate from the strategy to increase their
expected payoff. That is,

$\displaystyle\forall x\in
X,\;U^{(p)}(x,\xi^{*})\geqslant{U^{(p)}(x^{*},\xi^{*})}$ (23)
$\displaystyle\forall\theta\in\Theta,\;\xi(\theta)\in\Xi,\;U^{(a)}(x^{*},\xi(\theta))\leqslant
U^{(a)}(x^{*},\xi^{*}(\theta)),$ (24)

Given the adversary’s behavioral strategy $\xi^{*}$ (but not his type, which
is private information for the adversary), the dispatcher faces the following
optimization problem:

$\displaystyle x^{*}\in\arg\underset{x\in X}{\max}\;{U}^{(p)}(x,\xi^{*}),$
(25)

while given $x^{*}$ the adversary faces the maximization problem for all
$\theta\in\Theta$ from (21)

$\xi^{*}(\theta)\in\arg\underset{\xi(\theta)\in\Xi}{\min}\;{U}^{(a)}(x^{*},\xi(\theta)).$
(26)

To obtain a Bayesian equilibrium concept of the game is equivalent to solving
(25) and (26) simultaneously. Since the utility functions for both players are
convex and their strategy spaces are compact, we can know that there exists a
Bayesian equilibrium.

We have the following existence proposition upon an Bayesian Nash equilibrium
of the adversarial optimal transport game $\mathcal{G}$.

###### Proposition 1

Let $\mathcal{G}$ be the Bayesian games mentioned in definition 22. Let
$(x^{*},\\{\xi^{*}(\theta)\\}_{\theta\in\Theta})$ be a strategy profile
meeting the requirements of Bayesian equilibrium in definition 2. Then such an
equilibrium exists in the game $\mathcal{G}$.

Proof. According to Harsanyi Harsanyi (1968), it suffices to prove the
existence of a Nash equilibrium, which is obtained before both players learn
their types. We observe that the best response function for both players (26)
and (25) are upper semi-continuous. Meanwhile, referring to (15) and (17), we
know that the strategy spaces of both players are compact and convex. So by
Kakutani’s fixed-point theorem Kakutani (1941) there exists a Nash equilibrium
for the game where the types of the players are not revealed. As a result,
there exists a Bayesian equilibrium for the game $\mathcal{G}$.

###  4.2 Dynamic adversarial optimal transport

The crime activities within one region also vary by time as depicted in Figure
We now introduce the dynamic version of the Bayesian game $\mathcal{G}$ as in
definition 22 to further study the temporal relations of police resource
allocation under the presence of criminal behaviors.

We can define the history of the dispatcher’s plans $x^{(t)}\in X^{\otimes t}$
and the history of the criminal’s actions $\xi^{t}\in\Xi^{\otimes t}$ as
follows:

######  Definition 3 (History of the criminal’s actions and the dispatcher’s
plans)

For $t\in T$, we denote the history of the criminal’s actions
$\xi^{(t)}(\theta)\in\Xi^{t}$ and the history of the dispatcher’s plans
$x^{(t)}\in X^{\otimes t}$ as follows:

$\displaystyle x^{(t)}$ $\displaystyle:=(x^{1},x^{2},\dots,x^{t}),$ (27)
$\displaystyle\xi^{(t)}(\theta)$
$\displaystyle:=(\xi^{1}(\theta),\xi^{2}(\theta),\dots,\xi^{t}(\theta)),\forall\theta\in\Theta$

The planner’s utilities can be written as

$U^{t(p)}(x^{t},\xi^{t})=\sum_{\begin{subarray}{c}\\{j,q\\}\in\mathcal{E}_{j}\\\
\theta_{q}\in\\{1,2\\}\end{subarray}}{(m_{jq}+\theta_{q}\xi^{t}_{q}(\theta_{q}))\mu^{t}(\theta_{q})x^{t}_{jq}-\lambda
x^{t}_{jq}\log x^{t}_{jq}}.$ (28)

The criminal’s utilities

$\displaystyle\underset{\begin{subarray}{c}\xi^{t}(1)\in\Xi\\\
\xi^{t}(2)\in\Xi\end{subarray}}{\min}$
$\displaystyle{\sum_{\\{j,q\\}\in\mathcal{E}_{j}}{c_{jq}(\phi[\xi^{t}_{q}(1);\xi^{t-1}_{q}(1)]^{-\beta_{2}}+\phi[\xi^{t}_{q}(2);\xi^{t-1}_{q}(2)]^{-\beta_{2}})x^{\beta_{1}}_{jq}}}$
(29)
$\displaystyle+(m_{jq}+(\phi[\xi^{t}_{q}(1);\xi^{t-1}_{q}(1)]+\phi[\xi^{t}_{q}(2);\xi^{t-1}_{q}(2)])x_{jq},$
$\displaystyle\text{s.t}.$ $\displaystyle\;\;\sum_{j\in
J}{\xi^{t}_{q}(1)}\leqslant\underline{n}_{q},\;q\in Q$
$\displaystyle\;\sum_{j\in J}{\xi^{t}_{q}(2)}\leqslant\bar{n}_{q},\;\;q\in Q,$

where $\phi:\mathbb{R}\rightarrow\mathbb{R}$ is some smooth thresholding
function with the threshold set to be $\xi^{t-1}$.

setting up the threshold function $\phi$ A rough sketch of $\phi$ should be as
Figure 1. We aim to characterize a ‘thresholding behavior’ of $\xi^{t}$ upon
$\xi^{t-1}$. One choice is as follows:

$\phi(\xi^{t};\xi^{t-1})=\xi^{t-1}\mathbf{1}_{\\{\xi^{t}<\xi^{t-1}+\tau\\}}+(\xi^{t}-\tau)\mathbf{1}_{\\{\xi^{t}>\xi^{t-1}+\tau\\}},$
(30)

Figure 1: Possible thresholding function $\phi(\xi^{t};\xi^{t-1})$. Where
$\tau=\xi^{t-1}$ as the criminal’s action in the last stage plays the role of
threshold.

###  4.3 Formulation of dynamic Bayesian game

###  4.4 Solutions

For a dynamic Bayesian game, many equilibrium concepts can be raised Fudenberg
and Tirole (1991) such as perfect Bayesian Nash equilibrium, sequential
equilibrium, and even Bayesian Nash equilibrium. We adopt the following
perfect Bayesian Nash equilibrium (PBNE):

######  Definition 4 (Perfect Bayesian Nash equilibrium (PBNE))

Let $\mathcal{G}^{T}$ be the dynamic Bayesian game defined in A perfect
Bayesian Nash equilibrium for the dynamic Bayesian game $\mathcal{G}_{T}$
consists of the history of strategy profile $(x^{(t)*},\Xi^{(t)*})$ and a
belief system $\mu^{*}$ such that the following requirements are met:

  1. 1.

(Sequential rationality for the dispatcher) The dispatcher’s optimal transport
plan $x^{t*}$ at every stage is a maximizer of the dispatcher’s utilities
given the previous strategy profiles and the current belief $\mu^{t}$:

$x^{t}\in\underset{x}{\arg\max}{\;U^{t(p)}(x,\xi^{t*}|x^{(t)*},\xi^{(t)*})}$
(31)

  2. 2.

(Sequential rationality for the criminal) The criminal’s optimal strategy at
every stage $\xi^{t*}(\theta)$ given its type $\theta\in\\{1,2\\}^{|Q|}$ (the
type does not change through stage) is a minimizer of its own cost function
$U^{t(a)}$ given the history of the strategy profiles $x^{(t)*},\xi^{(t)*}$ as
follows:

$\xi^{t}(\theta)\in\underset{\xi(\theta)\in\Xi}{\arg\max}\;U^{t(a)}(x^{t*},\xi(\theta)|x^{(t)*},\xi^{(t)*}),\;\forall\theta\in\Theta.$
(32)

  3. 3.

(Belief update) The belief system of the game is updated after the
dispatcher’s strategy and the criminal’s strategy is reveal at stage $t$:

$\mu^{t+1}(\theta|\xi^{t})=\frac{\mu^{t}(\theta)\xi^{t}(\theta)}{\sum_{\theta^{\prime}\in\Theta}{{\mu^{t}(\theta^{\prime})\xi^{t}(\theta^{\prime})}}},\;\forall
t\in T,\theta\in\Theta.$ (33)

##  5 Distributed Algorithms

|
---|---
Figure 2: The NYC crime heat map expressed in regions in different time
periods. Left: May 2022, Right: June 2022.

##  6 Numerical Experiments

##### Crime Neighborhood

The crime city patrol has been a variety of problems in building smart cities.
In particular, a fair, in-time patrol schedule that allocates police resources
among neighborhoods is crucial in detecting and deterring urban crimes. The
most naive way of arranging police forces for patrolling is by planning
manual, fixed patrol schedules ahead of time and let the police to commit to
the schedule. However fixed patrol schedules do not respond well to
emergencies and opportunistic crimes that are committed during the
intermittence of shifts. Another branch of models are security pursuit evasion
games (PEG)Hespanha et al. (2000a), which formulates the police as pursuers
and criminals as evaders. However, the main goal of criminals is not just to
evade from police, but to benefit from committing crimes. Another landscape
model based on game theories are Stackelberg security games (SSG)Tambe (2011),
in which the criminals are followers and police are leaders. Authors in Pita
et al. (2011) exploits and implements a Stackelberg security game formulation
upon the on-demand allocation of security forces among airports. To explore
temporal correlations of criminals, authors in Brown et al. (2014) , formulate
the patrolling of polices as a Markov decision process (MDP). Authors in Zhang
et al. (2016b) formulates the patrolling game as a dynamic Bayesian network to
better predict the criminals’ migration among different neighborhoods.

All of the models optimizes police force allocation regarding a particular
goal such as minimizing response time, lowering the crime rate, etc. In real
practice the allocation of police needs to leverage between different goals
and so are criminals. Hence we apply our dynamical Bayeisan optimal transport
model into formulating police dispatch. We assume there are two types of
criminals: minor criminal offenders and major criminal offenders.Such types
are private that the police can only develop a common prior distribution yet
do not know the exact realizations. We assume that at every stage, criminals
act simultaneously with police patrolling and can adjust their locations of
committing crimes at the next moment by observing the police force allocation
for the moment. The police forces do not have access to the criminals’ actions
at previous stages, but criminals do and they can refer to their history of
actions to help make decisions for the current stage. We aim to deploy our
formulation of adversarial dynamic optimal transport into crime activities to
study temporal and spatial correlation of crimes and how crimes are affect in
the presences of police force allocations.

As a concrete example, we present numerical experiments of illustrate the
equilibrium achieved both in static and in dynamic Bayesian games to solve
optimal transport under the adversarial settings. We consider first a network
that contains 2 source nodes and 3 target nodes whose topology is depicted in
Figure 3.

I

Figure 3: The scheme of resource allocation problem with $m=3$ police stations
(source nodes) and $n=4$ neighborhoods (target nodes).

##### Parameter Setting

We set the capacities for the source nodes as $c_{1}=4,c_{2}=3$ and the
coefficient matrix (meaning the cost of transportation along every link) is
chosen as follows:

$M=\begin{bmatrix}1&3&5\\\ 2&5&1\end{bmatrix}.$ (34)

When there is criminal involved (the Bayesian adversarial OT case), the upper
bound matrix $N$ is chosen as

$\underline{n}=\begin{bmatrix}6&4&4\end{bmatrix},\bar{n}=\begin{bmatrix}8&10&10\end{bmatrix}.$
(35)

The coefficients related to probability of getting caught are

$C=\begin{bmatrix}1&2&3\end{bmatrix}.$ (36)

the indices $\beta_{1}=0.5,\beta_{2}=0.5$.

Figure 4: The optimal transport plan of police resource $x^{*}$ without the
fairness term. The dispatcher’s utility function is chosen to be (28) with
$\lambda=0$. We choose the perception matrix $M$ to be in (34).  Figure 5: The
optimal plan of police resource allocation $x$ with fairness term taken into
account. The perception coefficient $M$ is set in (34). The dispatcher’s
utility is set in (28) with the parameter $\lambda=3$.  Figure 6: The
convergence of criminal actions into equilibrium of the static Bayesian game
$\mathcal{G}$ in definition 2. We set the perception coefficient matrix $M$ as
in (34), the upper bounds of criminal $\underline{n},\bar{n}$ actions under
type $\theta_{q}=1,2$ to be in (35). The coefficient matrix $C$ corresponding
to the punishment is set in (36).

Figure 7: The dispatcher’s equilibrium strategies $x^{*}$ for the static
Bayesian adversarial optimal transport game $\mathcal{G}$. The perception
matrix $M$ is given in (34). The criminal’s upper bound activity is given in
(35). The parameter $\lambda=3$.

u

Figure 8: The iterations of the dispatcher’s and the criminal’s utilities for
the Bayesian adversarial optimal transport game $\mathcal{G}$.

Since every source node $j$ has two edges connected to it and the data on both
edges could be truthful or modified, we could associate the type space of
$j$’s adversary as a set containing four elements:
$\theta_{q}=\\{\theta_{j_{1}q},\theta_{j_{1}q_{2}}\\}$, where
$\theta_{j_{1}q_{1}},\theta_{j_{1}q_{2}}$ could both be $0$ or $1$. For
simplicity we associate a uniform distribution as a common prior over the
adversary’s type space, that is, the probability that every type is adopted is
$\frac{1}{16}$.

1\. We first illustrate the impact of the regularization term upon smoothing
the transport plans. We set $\lambda=0$ and $\lambda=3$ in adversary-free
scenario in (5) and illustrate the corresponding transport plans in Figure 4
and Figure 5, respectively. We observe that without regularization, every
source node tend to distribute all of their resources along the most weighted
link they are connected to. With the entropic regularization term, however,
every source node begins to distribute their resources more smoothly among its
links.

2\. The next example illustrates the influence of the adversary with unknown
type upon the social utility of the planner formulated in static games with
incomplete information as in Figure 8. We compute the Bayesian equilibrium
defined in plot the expected utilities of the dispatcher with and without
attack in (a). We set $\lambda=3$ and demonstrate the converging transport
plans as well as adversary’s actions as Bayesian equilibrium of the game in
Figure 7 and in Figure 6, respectively. We observe that while the adversary
modifies and increase the perception $m$ on some link, he increases the
transport rate along that link, which is decided by the dispatcher.

3\. Next we demonstrate how the criminal behavior varies through time.

##  7 Conclusion

We formulate the problem of seeking an optimal transport plan over networks
under the adversarial setting as a Bayesian game where the adversary and the
dispatcher have asymmetric information. Specifically, the dispatcher does not
know which of the source nodes are malicious and has to develop optimal
transport plan based on his subjective probability distribution over the
conditions of the source nodes. The asymmetric information happens often in
most security setting than complete information. We have also developed a
distributed algorithm to solve the Bayesian games, suggesting that large-scale
implementation of our formulation is possible since every user can have their
own transport rates updated simultaneously without knowing the global topology
of the network.

Extensions of the current work can go along several directions. One example is
to consider correlated strategies of the dispatchers/the dispatcher and the
adversary, while maintaining the assumption that the game is of incomplete
information for the dispatchers. Another direction is to consider various
information structures for the adversary and the dispatcher. For instance, the
dispatcher Our model also has also the potentiality to be applied large
communication networks or urban science context.

## References

  * Becker [1968] Gary S Becker.  Crime and punishment: An economic approach.  In _The economic dimensions of crime_ , pages 13–68. Springer, 1968\.
  * Brown et al. [2014] Matthew Brown, Sandhya Saisubramanian, Pradeep Varakantham, and Milind Tambe.  Streets: game-theoretic traffic patrolling with exploration and exploitation.  2014\.
  * Buttyan and Hubaux [2007] Levente Buttyan and Jean-Pierre Hubaux.  _Security and cooperation in wireless networks: thwarting malicious and selfish behavior in the age of ubiquitous computing_.  Cambridge University Press, 2007.
  * Curtin et al. [2010] Kevin M Curtin, Karen Hayslett-McCall, and Fang Qiu.  Determining optimal police patrol areas with maximal covering and backup covering location models.  _Networks and spatial economics_ , 10(1):125–145, 2010.
  * Cuturi [2013] Marco Cuturi.  Sinkhorn distances: Lightspeed computation of optimal transport.  In _Advances in neural information processing systems_ , pages 2292–2300, 2013.
  * Ehrlich [1973] Isaac Ehrlich.  Participation in illegitimate activities: A theoretical and empirical investigation.  _Journal of political Economy_ , 81(3):521–565, 1973.
  * Essid and Solomon [2018] Montacer Essid and Justin Solomon.  Quadratically regularized optimal transport on graphs.  _SIAM Journal on Scientific Computing_ , 40(4):A1961–A1986, 2018.
  * Felegyhazi et al. [2006] Mark Felegyhazi, J-P Hubaux, and Levente Buttyan.  Nash equilibria of packet forwarding strategies in wireless ad hoc networks.  _IEEE Transactions on Mobile computing_ , 5(5):463–476, 2006.
  * Ferradans et al. [2014a] Sira Ferradans, Nicolas Papadakis, Gabriel Peyré, and Jean-François Aujol.  Regularized discrete optimal transport.  _SIAM Journal on Imaging Sciences_ , 7(3):1853–1882, 2014a.
  * Ferradans et al. [2014b] Sira Ferradans, Nicolas Papadakis, Gabriel Peyré, and Jean-François Aujol.  Regularized discrete optimal transport.  _SIAM Journal on Imaging Sciences_ , 7(3):1853–1882, 2014b.
  * Ferris and Riveros [2013] James Michael Ferris and Gerry Edward Riveros.  Methods and systems for matching resource requests with cloud computing environments, March 19 2013.  US Patent 8,402,139.
  * Fudenberg and Tirole [1991] Drew Fudenberg and Jean Tirole.  Game theory, 1991.  _Cambridge, Massachusetts_ , 393(12):80, 1991\.
  * Galichon [2018] Alfred Galichon.  _Optimal transport methods in economics_.  Princeton University Press, 2018.
  * Garnaev and Trappe [2014] Andrey Garnaev and Wade Trappe.  Fair resource allocation under an unknown jamming attack: a bayesian game.  In _2014 IEEE International Workshop on Information Forensics and Security (WIFS)_ , pages 227–232. IEEE, 2014.
  * Genevay et al. [2016] Aude Genevay, Marco Cuturi, Gabriel Peyré, and Francis Bach.  Stochastic optimization for large-scale optimal transport.  In _Advances in neural information processing systems_ , pages 3440–3448, 2016.
  * Georgiadis et al. [2006] Leonidas Georgiadis, Michael J Neely, and Leandros Tassiulas.  _Resource allocation and cross-layer control in wireless networks_.  Now Publishers Inc, 2006.
  * Harsanyi [1967] John C Harsanyi.  Games with incomplete information played by “bayesian” players, i–iii part i. the basic model.  _Management science_ , 14(3):159–182, 1967.
  * Harsanyi [1968] John C Harsanyi.  Games with incomplete information played by “bayesian” players part ii. bayesian equilibrium points.  _Management Science_ , 14(5):320–334, 1968.
  * Hespanha et al. [2000a] Joao P Hespanha, Maria Prandini, and Shankar Sastry.  Probabilistic pursuit-evasion games: A one-step nash approach.  In _Proceedings of the 39th IEEE Conference on Decision and Control (Cat. No. 00CH37187)_ , volume 3, pages 2272–2277. IEEE, 2000a.
  * Hespanha et al. [2000b] Joao P Hespanha, Maria Prandini, and Shankar Sastry.  Probabilistic pursuit-evasion games: A one-step nash approach.  In _Proceedings of the 39th IEEE Conference on Decision and Control (Cat. No. 00CH37187)_ , volume 3, pages 2272–2277. IEEE, 2000b.
  * Huang et al. [2006] Jianwei Huang, Randall A Berry, and Michael L Honig.  Distributed interference compensation for wireless networks.  _IEEE Journal on Selected Areas in Communications_ , 24(5):1074–1084, 2006.
  * Hughes and Chen [2021a] Jason Hughes and Juntao Chen.  Resilient and distributed discrete optimal transport with deceptive adversary: A game-theoretic approach.  _IEEE Control Systems Letters_ , 2021a.
  * Hughes and Chen [2021b] Jason Hughes and Juntao Chen.  Fair and distributed dynamic optimal transport for resource allocation over networks.  In _2021 55th Annual Conference on Information Sciences and Systems (CISS)_ , pages 1–6. IEEE, 2021b.
  * Jiang et al. [2013] Albert Xin Jiang, Zhengyu Yin, Chao Zhang, Milind Tambe, and Sarit Kraus.  Game-theoretic randomization for security patrolling with dynamic execution uncertainty.  In _AAMAS_ , pages 207–214, 2013.
  * Kakutani [1941] Shizuo Kakutani.  A generalization of brouwer’s fixed point theorem.  _Duke mathematical journal_ , 8(3):457–459, 1941\.
  * Kantorovich [1942] Leonid Vitalievich Kantorovich.  On the translocation of masses.  In _Dokl. Akad. Nauk. USSR (NS)_ , volume 37, pages 199–201, 1942\.
  * Maschler et al. [2013] Michael Maschler, Eilon Solan, and Shmuel Zamir.  _Game Theory_.  Cambridge University Press, 2013.  doi: 10.1017/CBO9780511794216.
  * Monderer and Shapley [1996] Dov Monderer and Lloyd S Shapley.  Potential games.  _Games and economic behavior_ , 14(1):124–143, 1996.
  * Pita et al. [2011] James Pita, Milind Tambe, Christopher Kiekintveld, Shane Cullen, and Erin Steigerwald.  Guards—innovative application of game theory for national airport security.  In _Twenty-Second International Joint Conference on Artificial Intelligence_ , 2011.
  * Raya et al. [2008] Maxim Raya, Mohammad Hossein Manshaei, Márk Félegyházi, and Jean-Pierre Hubaux.  Revocation games in ephemeral networks.  In _Proceedings of the 15th ACM conference on Computer and communications security_ , pages 199–210, 2008.
  * Shao et al. [2020] Guangru Shao, Rui Wang, Xue-Fang Wang, and Kun-Zhi Liu.  Distributed algorithm for resource allocation problems under persistent attacks.  _Journal of the Franklin Institute_ , 357(10):6241–6256, 2020.
  * Srikant and Ying [2013a] Rayadurgam Srikant and Lei Ying.  _Communication networks: an optimization, control, and stochastic networks perspective_.  Cambridge University Press, 2013a.
  * Srikant and Ying [2013b] Rayadurgam Srikant and Lei Ying.  _Communication networks: an optimization, control, and stochastic networks perspective_.  Cambridge University Press, 2013b.
  * Tambe [2011] Milind Tambe.  _Security and game theory: algorithms, deployed systems, lessons learned_.  Cambridge university press, 2011.
  * Villani [2008] Cédric Villani.  _Optimal transport: old and new_ , volume 338.  Springer Science & Business Media, 2008.
  * Zhang et al. [2016a] Chao Zhang, Shahrzad Gholami, Debarun Kar, Arunesh Sinha, Manish Jain, Ripple Goyal, and Milind Tambe.  Keeping pace with criminals: An extended study of designing patrol allocation against adaptive opportunistic criminals.  _Games_ , 7(3):15, 2016a.
  * Zhang et al. [2016b] Chao Zhang, Shahrzad Gholami, Debarun Kar, Arunesh Sinha, Manish Jain, Ripple Goyal, and Milind Tambe.  Keeping pace with criminals: An extended study of designing patrol allocation against adaptive opportunistic criminals.  _Games_ , 7(3):15, 2016b.
  * Zhang and Zhu [2019] Rui Zhang and Quanyan Zhu.  Consensus-based distributed discrete optimal transport for decentralized resource matching.  _IEEE Transactions on Signal and Information Processing over Networks_ , 5(3):511–524, 2019.

## Appendix A. Placeholder