# Extreme points of general transportation polytopes

Patrice Koehl
Department of Computer Science
University of California, Davis Davis, CA 95616
<EMAIL_ADDRESS>
###### Abstract

Transportation matrices are $m\times n$ non-negative matrices whose row sums
and row columns are equal to, or dominated above with given integral vectors
$R$ and $C$. Those matrices belong to a convex polytope whose extreme points
have been previously characterized. In this article, a more general set of
non-negative transportation matrices is considered, whose row sums are bounded
by two integral non-negative vectors $R_{min}$ and $R_{max}$ and column sums
are bounded by two integral non-negative vectors $C_{min}$ and $C_{max}$. It
is shown that this set is also a convex polytope whose extreme points are then
fully characterized.

_Keywords_ transportation polytopes $\cdot$ assignment problems $\cdot$
stochastic matrices

##  1 Introduction

Imagine that there are $n$ agents that are assigned to perform $m$ tasks, and
assume that each agent can perform an integer number of tasks and each task
can only be assigned to an integer number of agents. A manager in charge of
the assignment agent / task will consider the “cost” of assigning each agent
to each task and will design an “optimal assignment”, namely an assignment
that leads to a minimal overall cost. Finding a solution to this seemingly
simple practical problem has led to the development of a gem in the
mathematics and statistics communities, namely the optimal transport (OT)
problem.

A mathematical formulation of the problem is to consider two sets of points
$S_{1}$ of size $n$ (the agents) and $S_{2}$ of size $m$ (the tasks). Each
point $k$ in $S_{1}$ (resp $S_{2}$) is assigned a “mass” $m_{1}(k)$, the
amount of tasks it can perform (resp $m_{2}(k)$, the number of agents needed
to perform this task). The cost of an assignment between $S_{1}$ and $S_{2}$
is encoded as a non-negative matrix $T=[T(i,j)]$ with $i\in\\{1\ldots,n\\}$
and $j\in\\{1\ldots,m\\}$. The OT problem can then be formulated as finding a
matrix $G$ of correspondence between points in $S_{1}$ and points in $S_{2}$
that minimizes the total transport cost $U$ defined as

$U(G)=\sum_{i=1}^{n}\sum_{j=1}^{m}T(i,j)G(i,j).$

The minimum of $U$ is to be found for the values of $G(i,j)$ that satisfy the
following constraints

$\displaystyle\forall(i,j),\quad G(i,j)\geq 0,$ (1a) $\displaystyle\forall i,\
\ \quad\sum\limits_{j}G(i,j)=m_{1}(i),$ (1b) $\displaystyle\forall j,\ \
\quad\sum\limits_{i}G(i,j)=m_{2}(j).$ (1c)

What makes the OT problem so interesting is that its solution includes two
essential components. First, it defines a distance between the distributions
considered ($m_{1}$ and $m_{2}$). These distances have enabled statisticians
and mathematicians to derive a geometric structure on the space of probability
distributions (Villani [2008], Peyré and Cuturi [2018]). Second, it also
provides the optimal transportation plan $G$ between the distributions; this
optimal plan defines a registration, thereby enabling alignment between the
distributions. Applications of OT have exploded in the recent years, in
domains such as applied mathematics, machine learning, computer vision, and
linguistics (see Villani [2008], Cotar et al. [2013], Santambrogio [2015],
Peyré and Cuturi [2018], among others).

The formulation defined above corresponds to the Monge-Kantorovich version of
the optimal transport problem (Kantorovich [1942], Villani [2008]). A
simplified version of the OT problem is to set the masses to 1 (i.e. each
agent can only perform one task and each task is performed by one agent) and
to assume balance (same number of agents and tasks). The transportation matrix
$A$ then becomes integral (its values are integers) and satisfies the
following constraints:

$\displaystyle\forall(i,j),\quad A(i,j)\in\\{0,1\\},$ (2a)
$\displaystyle\forall i,\ \ \quad\sum\limits_{j}A(i,j)=1,$ (2b)
$\displaystyle\forall j,\ \ \quad\sum\limits_{i}A(i,j)=1.$ (2c)

This simplified problem that has its root in Monge’s original transport
problem Monge [1781] is itself a classical problem in combinatorial
optimization referred to as the assignment problem or alternatively, using the
language of graph theory, as the bipartite weighted matching problem (for a
comprehensive analysis of assignment problems, see for example Burkard et al.
[2009]).

Our focus in this paper is on the structures of the set of matrices that
either satisfy the constraints 1 or 2. Those sets of matrices belong to so-
called transportation polytopes that have been studied extensively in the
literature (see for example, Birkhoff [1946], Brualdi and Gibson [1976],
Brualdi [1980], Muchlis [1992], Gill and Linusson [2009], Ziegler [2012], De
Loera and Kim [2013], Cao and Chen [2019], Chen et al. [2021]). In the special
case $n=m$ and $m_{1}(i)=m_{2}(j)=1\quad\forall(i,j)$, the matrices $G$ are
doubly stochastic. The set of all those matrices form a convex polytope
$\Omega_{n}$. Under the same conditions, the set of matrices $A$ is the set of
permutations $P_{n}$. The set $P_{n}$ is the set of extreme points of
$\Omega_{n}$; as such, any matrix $G$ in $\Omega_{n}$ can be expressed as a
linear combination of permutation matrices in $P_{n}$:

######  Theorem 1.

(Birkhoff [1946], Von Neumann [1953]) An $n\times n$ matrix G is doubly
stochastic if and only if there is a finite set of permutations matrices
$P_{1},\ldots,P_{N}$ and corresponding non-negative real numbers
$\alpha_{1},\ldots,\alpha_{N}$ with $\alpha_{1}+\ldots+\alpha_{N}=1$ such that
$A=\alpha_{1}P_{1}+\ldots+\alpha_{N}P_{N}$.

There are many proofs available for this theorem, referred to as the Birkhoff
- von Neumann theorem; we only refer here to the original proofs Birkhoff
[1946], Von Neumann [1953]. This theorem has proved useful to establish
convergence for algorithms developed to solve the assignment problem (see for
example Kosowsky and Yuille [1994], Koehl and Orland [2021].

The equations 2 refer to a very specific balanced assignment problem. There
are, however many more types of assignment problems, usually referred to as
unbalanced problems. They consider a number of agents that differ from the
number of tasks (using the example and terminology from above). The most
common formulation assigns a single agent to one task, leaving some agents
and/or some task unmatched. It is referred to in the literature as the
$k$-cardinality assignment problem (Dell’Amico and Martello [1997]). Some
other formulations allow for multiple jobs to be assigned to the same agent to
compensate for the imbalance (when the number of tasks is bigger than the
number of agents), with possibly additional constraints such as each agent is
allocated at least one task. Transportation matrices for such problems belong
to a more general class of transportation polytopes than the doubly stochastic
matrices and doubly substochastic matrices. In this paper, we propose to
characterize those transportation polytopes.

We define these more general assignment problems as follows. The cost of
transport between $S_{1}$ and $S_{2}$ is still encoded as a non-negative
matrix $T=[T(i,j)]$ with $i\in\\{1\ldots,n\\}$ and $j\in\\{1\ldots,m\\}$. The
general assignment problem is formulated as finding a matrix $G$ of
correspondence between points in $S_{1}$ and points in $S_{2}$ that minimizes
the total transport cost $U=\sum_{i=1}^{n}\sum_{j=1}^{m}T(i,j)G(i,j)$ whose
values $G(i,j)$ satisfy the following constraints

$\displaystyle\forall(i,j),\ \ \quad G(i,j)\in\\{0,1\\},$ (3a)
$\displaystyle\forall i,\ \ \quad r_{i}\leq\sum\limits_{j}G(i,j)\leq R_{i},$
(3b) $\displaystyle\forall j,\ \ \quad c_{j}\leq\sum\limits_{i}G(i,j)\leq
C_{j},$ (3c)

where $r_{i}$ and $R_{i}$ are given positive integers satisfying $0\leq
r_{i}\leq R_{j}$ and similarly $c_{j}$ and $C_{j}$ are given positive integers
satisfying $0\leq c_{j}\leq C_{j}$. A relaxed version of these constraints is
to replace equation 3a with $0\leq G(i,j)\leq 1$, i.e. allowing $G(i,j)$ to be
fractional. Finally, a possible additional constraint is to preset the total
number of assignments between $S_{1}$ and $S_{2}$ to a given integer value
$k$, i.e. $\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{m}G(i,j)=k$. Possible
solutions of the relaxed general assignment problem belong to sets of matrices
that are defined below.

######  Definition 1.

Let $R_{min}=(r(1),\ldots,r(n))$, $R_{max}=(R(1),\ldots,R(n))$,
$C_{min}=(c(1),\ldots,c(m))$, and $C_{max}=(C(1),\ldots,C(m))$ be four non-
negative integral vectors satisfying

$\displaystyle\forall i,\ \ \quad 0\leq r(i)\leq R(i),$ $\displaystyle\forall
j,\ \ \quad 0\leq c(j)\leq C(j).$

Let $A$ be a non-negative matrix of size $n\times m$ and let us denote the row
sum vector of $A$ as $RA$ and the column sum vector of $A$ as $CA$. Let us
also define $\sigma(A)$ to be the sum of all elements of $A$, i.e.
$\sigma(A)=\sum_{i=1}^{n}\sum_{j=1}^{m}A(i,j)$. The transportation polytope
$\mathcal{U}(R_{min}^{max},C_{min}^{max})$ is the set of $n\times m$ matrices
$A$ that satisfy

$\displaystyle\forall(i,j),\ \ \quad 0\leq A(i,j)\leq 1,$
$\displaystyle\forall i,\ \ \quad r(i)\leq RA(i)\leq R(i),$
$\displaystyle\forall j,\ \ \quad c(j)\leq CA(j)\leq C(j).$

The transportation polytope $\mathcal{U}^{k}(R_{min}^{max},C_{min}^{max})$ is
the set of $n\times m$ matrices $A$ that satisfy

$\displaystyle\mathcal{U}^{k}(R_{min}^{max},C_{min}^{max})=\left\\{A\in\mathcal{U}(R_{min}^{max},C_{min}^{max})|\quad\sigma(A)=k\right\\}.$

Denote $\mathcal{P}(R_{min}^{max},C_{min}^{max})$ the set of all matrices in
$\mathcal{U}(R_{min}^{max},C_{min}^{max})$ whose entries are either 0 or 1
with a similar definition for $\mathcal{P}^{k}(R_{min}^{max},C_{min}^{max})$
with respect to $\mathcal{U}^{k}(R_{min}^{max},C_{min}^{max})$. We now state
our main result.

######  Theorem 2.

The transportation polytopes $\mathcal{U}(R_{min}^{max},C_{min}^{max})$ and
$\mathcal{U}^{k}(R_{min}^{max},C_{min}^{max})$ satisfy the following
properties:

  * a)

$\mathcal{U}(R_{min}^{max},C_{min}^{max})$ is the convex hull of all matrices
in $\mathcal{P}(R_{min}^{max},C_{min}^{max})$,

  * b)

$\mathcal{U}^{k}(R_{min}^{max},C_{min}^{max})$ is the convex hull of all
matrices in $\mathcal{P}^{k}(R_{min}^{max},C_{min}^{max})$.

### Remarks:

  * i)

If $R_{min}=R_{max}=1_{n}$, $C_{min}=C_{max}=1_{m}$ (where $1_{n}$ and $1_{m}$
are vectors of one of size $n$ and $m$, respectively), $n=m=k$,
$\mathcal{U}^{k}(R_{min}^{max},C_{min}^{max})=\Omega_{n}$ and
$\mathcal{P}^{k}(R_{min}^{max},C_{min}^{max})=P_{n}$ and theorem 2 is then
equivalent to the Birkhoff-Von Neuman theorem for doubly stochastic matrices,
theorem 1.

  * ii)

If $R_{min}=0_{n}$, $C_{min}=0_{m}$, $R_{max}=1_{n}$, $C_{max}=1_{m}$,
$\mathcal{U}(R_{min}^{max},C_{min}^{max})$ is the set of doubly substochastic
matrices, $\mathcal{P}(R_{min}^{max},C_{min}^{max})$ is the set of
subpermutation matrices; a specific version of theorem 2 was established (see
for example Mirsky [1959], Čihák [1970]).

  * iii)

Similar to case ii), if $R_{min}=0_{n}$, $C_{min}=0_{m}$, $R_{max}=1_{n}$,
$C_{max}=1_{m}$, $\mathcal{U}^{k}(R_{min}^{max},C_{min}^{max})$ is the set of
doubly substochastic matrices with total sum $k$, and
$\mathcal{P}(R_{min}^{max},C_{min}^{max})$ is the set of subpermutation
matrices of rank $k$; a specific version of theorem 2 was established by
Mendelsohn and Dulmage for square matrices (Mendelsohn and Dulmage [1958]),
and later by Brualdi and Lee for rectangular matrices (Brualdi and Lee
[1978]).

##  2 A simple proof of Birkhoff - von Neumann theorem

As mentioned above, there are many proofs available for the Birkhoff - von
Neumann theorem, some of which belong now to textbooks. Here we describe a
simple proof. It is not original, but will serve as the basis for elements of
the proof of theorem 2.

First, we note that it is straightforward to show that $\Omega_{n}$ is a non-
empty compact convex set in $\mathbb{R}^{n\times m}$ and that any matrix in
$P_{n}$ is an extreme point of $\Omega_{n}$. Since any non-empty compact
convex set is the convex hull of its extreme points, to finish the proof, we
only need to show that any matrix $A\in\Omega_{n}\smallsetminus P_{n}$ is not
an extreme point of $\Omega_{n}$.

We prove first the following claim.

######  Claim 1.

Let $A\in\Omega_{n}$. If a row or column of $A$ contains a fractional value,
then it contains at least two.

###### Proof.

Let $A(i,j)$ be a fractional value of $A$. Since

$\displaystyle\sum_{k}A(i,k)=1,$

there exists $j_{2}\in[1,n]\smallsetminus\\{j\\}$ such that $A(i,j_{2})$ is
fractional. Similarly, since

$\displaystyle\sum_{k}A(k,j)=1,$

there exists $i_{2}\in[1,n]\smallsetminus\\{i\\}$ such that $A(i_{2},j)$ is
fractional. ∎

Let now $A$ be a matrix in $\Omega_{n}\smallsetminus P_{n}$. There exists a
pair $(i_{1},j_{1})$ such that $A(i_{1},j_{1})$ is fractional. Based on claim
1, there exists $j_{2}\in[1,n]$ with $j_{2}\neq j_{1}$ such that
$A(i_{1},j_{2})$ is non integral. Similarly, we can find $i_{2}\in[1,n]$ with
$i_{2}\neq i_{1}$ such that $A(i_{2},j_{2})$ is fractional. We can continue in
this manner, leading to a path
$\left((i_{1},j_{1}),(i_{1}j_{2}),\ldots\right)$ with fractional values in
$A$. As $n$ and $m$ are finite, we will ultimately reach a pair that we have
already visited. This means that we have identified a loop $L$ among all edges
between $S_{1}$ and $S_{2}$; the cardinality of this cycle is even (bipartite
graph). We write this cycle as

$\displaystyle L=\\{(a_{1},b_{1}),(a_{2},b_{2}),\ldots,(a_{2M},b_{2M})\\},$

where $2M=|L|$. We define the matrix $N$

$\displaystyle\begin{cases}N(i,j)&=0\quad(i,j)\notin C\\\
N(a_{2k},b_{2k})&=1\quad k\in\\{1,\ldots,M\\}\\\ N(a_{2k-1},b_{2k-1})&=-1\quad
k\in\\{1,\ldots,M\\}.\end{cases}$

Let us now define

$\displaystyle\epsilon_{max}=\min\\{A(a_{1},b_{1}),\ldots,A(a_{2M},b_{2M}),1-A(a_{1},b_{1}),\ldots,1-A(a_{2M},b_{2M})\\}.$

As all elements in the loop $L$ are fractional, $0<\epsilon_{max}<1$. For
$\epsilon\in(0,\epsilon_{max}]$, we define $E_{1}=A+\epsilon N$ and
$E_{2}=A-\epsilon N$. As two consecutive pairs in $L$ lead to the addition and
subtraction of the same quantity $\epsilon$ on one row or one column of $A$,
it is easy to verify that $E_{1}$ and $E_{2}$ are doubly stochastic and
therefore belong to $\Omega_{n}$. Since $A=\frac{1}{2}(E_{1}+E_{2})$, $A$ is
not an extreme point of $\Omega_{n}$.

##  3 The polytopes $\mathcal{U}(R,C)$ and $\mathcal{U}^{k}(R,C)$

Let us start with the simple case for which $R_{min}=R_{max}=R$ and
$C_{min}=C_{max}=C$. We rewrite $\mathcal{U}(R_{min}^{max},C_{min}^{max})$ and
$\mathcal{U}^{k}(R_{min}^{max},C_{min}^{max})$ as simply $\mathcal{U}(R,C)$
and $\mathcal{U}^{k}(R,C)$. Those polytopes have been studied extensively (see
for example Jurkat and Ryser [1967], Brualdi [1980], Rothblum and Schneider
[1989], Brualdi and Hwang [2005], daFonseca and Mamede [2009], Brualdi [2006],
Cavenagh [2013], Chen et al. [2016]). We prove first the following lemma.

######  Lemma 1.

The polytope $\mathcal{U}(R,C)$ is non-empty if and only if
$\sum_{i=1}^{n}R(i)=\sum_{j=1}^{m}C(j)$. In addition, $\mathcal{U}^{k}(R,C)$
is non empty if and only if $k=\sum_{i=1}^{n}R(i)$, in which case
$\mathcal{U}^{k}(R,C)=\mathcal{U}(R,C)$.

###### Proof.

Let $A$ be a matrix in $\mathcal{U}(R,C)$. Then,

$\displaystyle\sigma(A)$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{m}A(i,j)=\sum_{i=1}^{n}\left(\sum_{j=1}^{m}A(i,j)\right)=\sum_{i=1}^{n}R(i)$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{m}\left(\sum_{i=1}^{n}A(i,j)\right)=\sum_{j=1}^{m}C(j).$

Therefore, $\sum_{i=1}^{n}R(i)=\sum_{j=1}^{m}C(j)$. Conversely, if
$S=\sum_{i=1}^{n}R(i)=\sum_{j=1}^{m}C(j)$, it is straightforward to build a
matrix $A$ that belongs to $\mathcal{U}(R,C)$ (for example
$A=\frac{1}{S}RC^{T}$).

Let $A$ be a matrix in $\mathcal{U}^{k}(R,C)$. By definition, $k=\sigma(A)$.
Since $\mathcal{U}^{k}(R,C)\subset\mathcal{U}(R,C)$, we get
$k=\sum_{i=1}^{n}R(i)$. In addition, let $A$ be a matrix in
$\mathcal{U}(R,C)$. Then $\sigma(A)=\sum_{i=1}^{n}R(i)=k$, therefore
$A\in\mathcal{U}^{k}(R,C)$. ∎

In the case considered here, theorem 2 is reduced to the following lemma:

######  Lemma 2.

The transportation polytopes $\mathcal{U}(R,C)$ and $\mathcal{U}^{k}(R,C)$
satisfy the following properties:

  * a)

$\mathcal{U}(R,C)$ is the convex hull of all matrices in $\mathcal{P}(R,C)$,

  * b)

$\mathcal{U}^{k}(R,C)$ is the convex hull of all matrices in
$\mathcal{P}^{k}(R,C)$.

The proof of property $a)$ of lemma 2 is similar to the simple proof provided
above for the Birkhoff - von Neumann theorem, once we replace claim 1 with the
following claim:

######  Claim 2.

Let $A\in\mathcal{U}(R,C)$. If a row or column of $A$ contains a fractional
value, then it contains at least two.

###### Proof.

Let $A(i,j)$ be a fractional value of $A$. Since

$\displaystyle\sum_{k}A(i,k)=R(i)$

and $R(i)$ is a non-negative integer, necessarily $R(i)\neq 0$ and there
exists $j_{2}\in[1,n]\smallsetminus\\{j\\}$ such that $A(i,j_{2})$ is
fractional. Similarly, since

$\displaystyle\sum_{k}A(k,j)=C(j)$

and $C(j)$ is a non-negative integer, necessarily $C(j)\neq 0$ there exists
$i_{2}\in[1,n]\smallsetminus\\{i\\}$ such that $A(i_{2},j)$ is fractional. ∎

The proof of part $b)$ requires in addition that the matrices $E_{1}$ and
$E_{2}$ built from the matrix $A$ in $\mathcal{U}^{k}(R,C)$ also belong to
$\mathcal{U}^{k}(R,C)$. This is a direct consequence of the fact that the loop
identified in the matrix $A$ has en even number of elements. As such the
matrix $N$ defined in equation Case 2) satisfies $\sigma(N)=0$, and therefore
$\sigma(E_{1})=\sigma(E_{2})=\sigma(A)=k$, i.e. $E_{1}$ and $E_{2}$ belong to
$\mathcal{U}^{k}(R,C)$.

##  4 The polytope $\mathcal{U}(R_{min}^{max},C_{min}^{max})$

We start with the following lemma:

######  Lemma 3.

$\mathcal{U}(R_{min}^{max},C_{min}^{max})$ is convex.

###### Proof.

The proof of lemma 3 is relatively straightforward. We provide it here for
sake of completeness. Let $A$ and $B$ be two matrices belonging to
$\mathcal{U}(R_{min}^{max},C_{min}^{max})$, $\alpha$ a real number in $[0,1]$,
and $E=\alpha A+(1-\alpha)B$. We have:

$\displaystyle E(i,j)=\alpha A(i,j)+(1-\alpha)B(i,j).$

As both $A(i,j)$ and $B(i,j)$ belong to $[0,1]$ and $[0,1]$ is convex,
$E(i,j)\in[0,1]$. Also, for $i\in[1,n]$,

$\displaystyle\sum_{j=1}^{m}E(i,j)=\alpha\sum_{j=1}^{m}A(i,j)+(1-\alpha)\sum_{j=1}^{n}B(i,j).$

As $\sum_{j=1}^{m}A(i,j)$ and $\sum_{j=1}^{m}B(i,j)$ belong to $[r(i),R(i)]$
and this interval is convex, $\sum_{j=1}^{n}E(i,j)\in[r(i),R(i)]$. This is
true for all $i\in[1,n]$.

Similarly,

$\displaystyle\sum_{i=1}^{n}E(i,j)=\alpha\sum_{i=1}^{n}A(i,j)+(1-\alpha)\sum_{1=1}^{n}B(i,j).$

As $\sum_{i=1}^{n}A(i,j)$ and $\sum_{i=1}^{n}B(i,j)$ belong to $[c(j),C(j)]$
and this interval is convex, $\sum_{i=1}^{n}E(i,j)\in[c(j),C(j)]$. This is
true for all $j\in[1,m]$.

Therefore, $\alpha A+(1-\alpha)B$ belongs to
$\mathcal{U}(R_{min}^{max},C_{min}^{max})$ and this set is convex. ∎

######  Lemma 4.

If $A$ belongs to $\mathcal{P}(R_{min}^{max},C_{min}^{max})$, then $A$ is an
extreme point of $\mathcal{U}(R_{min}^{max},C_{min}^{max})$.

###### Proof.

Again, the proof of lemma 4 is relatively straightforward. We provide it here
for sake of completeness. Let $A$ be a matrix of
$\mathcal{P}(R_{min}^{max},C_{min}^{max})$ and let us suppose that
$A=\frac{1}{2}(E+F)$ where $E$ and $F$ are two distinct matrices in
$\mathcal{U}(R_{min}^{max},C_{min}^{max})$.

For any entry $A(i,j)=0$, we have:

$\displaystyle\frac{1}{2}(E(i,j)+F(i,j))=0.$

Since $E(i,j)\geq 0$ and $F(i,j)\geq 0$, we have $E(i,j)=F(i,j)=0$.

For any entry $A(i,j)=1$, we have:

$\displaystyle\frac{1}{2}(E(i,j)+F(i,j))=1,$

i.e.,

$\displaystyle E(i,j)+F(i,j)=2.$

Since $0\leq E(i,j)\leq 1$ and $0\leq F(i,j)\leq 1$, we have
$E(i,j)=F(i,j)=1$.

Therefore, $A=E=F$, which is in contradiction with the hypothesis that $E$ and
$F$ are distinct. Therefore, $A$ is not the midpoint of a line segment whose
endpoints are in $\mathcal{U}(R_{min}^{max},C_{min}^{max})$, and so $A$ is an
extreme point of $\mathcal{U}(R_{min}^{max},C_{min}^{max})$. ∎

Finally, we prove the following lemma:

######  Lemma 5.

If a matrix $A$ is an extreme point of
$\mathcal{U}(R_{min}^{max},C_{min}^{max})$ then it belongs to
$\mathcal{P}(R_{min}^{max},C_{min}^{max})$.

###### Proof.

We use a proof by contrapositive, loosely inspired by the proof in Deng [2019]
for doubly substochastic matrices. Let $A$ be a matrix in
$\mathcal{U}(R_{min}^{max},C_{min}^{max})$.

We start with some definitions. An entry in $A$ being neither 0 nor 1 is
called a fractional entry, and a row (resp. column) containing at least one
fractional entry is called a fractional row (resp. column). A fractional line
is either a fractional row or a fractional column.

As
$A\in\mathcal{U}(R_{min}^{max},C_{min}^{max})\smallsetminus\mathcal{P}(R_{min}^{max},C_{min}^{max})$,
it has at least one fractional entry and therefore one fractional line. We
consider two cases.

  * Case 1)

All fractional lines include at least two fractional values.

The proof is then very similar to the simple proof provided above for the
Birkhoff - von Neumann theorem, with no need for claim 1. This means that we
can find a loop $L$ within $A$ whose cardinality is even (bipartite graph). As
before, we write this cycle as

$\displaystyle L=\\{(a_{1},b_{1}),(a_{2},b_{2}),\ldots,(a_{2M},b_{2M})\\},$

where $2M=|L|$. We then define the matrix $N$

$\displaystyle\begin{cases}N(i,j)&=0\quad(i,j)\notin C\\\
N(a_{2k},b_{2k})&=1\quad k\in\\{1,\ldots,M\\}\\\ N(a_{2k-1},b_{2k-1})&=-1\quad
k\in\\{1,\ldots,M\\}.\end{cases}$

We also define

$\displaystyle\epsilon_{max}=\min\\{A(a_{1},b_{1}),\ldots,A(a_{2M},b_{2M}),1-A(a_{1},b_{1}),\ldots,1-A(a_{2M},b_{2M})\\}.$

As all elements in the loop $L$ are fractional, $0<\epsilon_{max}<1$. For
$\epsilon\in(0,\epsilon_{max}]$, we define $E_{1}=A+\epsilon N$ and
$E_{2}=A-\epsilon N$. As two consecutive pairs in $L$ leads to the addition
and subtraction of the same quantity $\epsilon$ on one row or one column, it
is easy to verify that the row sums and column sums of $E_{1}$ and $E_{2}$ are
equal to the row sums and row columns $A$, and therefore $E_{1}$ and $E_{2}$
belong to $\mathcal{U}(R_{min}^{max},C_{min}^{max})$. Since
$A=\frac{1}{2}(E_{1}+E_{2})$, $A$ is not an extreme point of
$\mathcal{U}(R_{min}^{max},C_{min}^{max})$.

  * Case 2)

There exists at least one fractional line that includes a single fractional
value.

Let us assume that one such fractional line is a column $j_{1}$ (the proof
would be the same if it were a row) and let $A(i_{1},j_{1})$ be the only
fractional value on $j_{1}$. If row $i_{1}$ only contains one fractional
value, we terminate. Otherwise $i_{1}$ includes at least two fractional values
and therefore we can find $j_{2}\neq j_{1}$ such that $A(i_{1},j_{2})$ is
fractional. We continue in this manner, until either (a) we have reached a
pair that we have already visited, or (b) we have reached a line that only
contains one fractional value. If we are in case (a), we have identified a
loop of fractional values whose cardinality is even and the proof is then
similar to Case 1. If we are in case (b), we have identified an incomplete
loop, i.e. a cycle. If the cycle ends because the following column had a
single fractional entry, the cardinality of the cycle is even, and we can
again follow exactly the proof from Case 1. Otherwise, the cycle has an odd
number of elements; we write it as

$\displaystyle
C=\\{(i_{1},j_{1}),(i_{2},j_{2}),\ldots,(i_{2K-1},j_{2K-1})\\},$

where $2K-1=|C|$. We define the matrix $N$ as:

$\displaystyle\begin{cases}N(i,j)&=0\quad(i,j)\notin C\\\
N(i_{2k},j_{2k})&=1\quad k\in\\{1,\ldots,K-1\\}\\\
N(i_{2k-1},j_{2k-1})&=-1\quad k\in\\{1,\ldots,K\\}.\end{cases}$

Let
$\epsilon_{max}=\min\\{A(i_{1},j_{1}),\ldots,A(i_{2K-1},j_{2K-1}),1-A(i_{1},j_{1}),\ldots,1-A(i_{2K-1},j_{2K-1})\\}$.
Since all the values $A(i,j)$ when $(i,j)\in C$ are fractional,
$0<\epsilon_{max}<1$.

Let us define $E_{1}(\epsilon)=A+\epsilon N$ and $E_{2}(\epsilon)=A-\epsilon
N$. If $0<\epsilon<\epsilon_{max}$, all values in $E_{1}(\epsilon)$ and
$E_{2}(\epsilon)$ are in the interval $[0,1]$. In addition, by construction,
any line sum of $E_{1}(\epsilon)$ and $E_{2}(\epsilon)$ is equal to the sum of
the same line of $A$, with the exception of column $j_{1}$ and row $i_{2K-1}$
. Those two lines, however, include a single fractional value. As such, the
column sum $CA(j_{1})$ and the row sum $RA(i_{2K-1})$ in $A$ are fractional.
Since the bounds on those row / column are integer values,

$\displaystyle r(i_{2K-1})$ $\displaystyle<RA(i_{2K-1})<R(i_{2K-1}),$
$\displaystyle c(j_{1})$ $\displaystyle<CA(j_{1})<C(j_{1}),$

and it is possible to modify the values of $A(i_{1},j_{1})$ and
$A(i_{2K-1},j_{2K-1})$ without compromising the constraints on the
corresponding column and row, respectively. Let us define
$\epsilon_{2}=\min\\{CA(j_{1})-c(j_{1}),RA(i_{2K-1})-r(i_{2K-1}),C(j_{1})-CA(j_{1}),R(i_{2K-1})-RA(i_{2K-1})\\}$;
note the $\epsilon_{2}$ is strictly fractional, i.e. $0<\epsilon_{2}<1$. If we
set $\epsilon$ such that $0<\epsilon<\min(\epsilon_{max},\epsilon_{2})$, then
$E_{1}(\epsilon)$ and $E_{2}(\epsilon)$ belong to
$\mathcal{U}(R_{min}^{max},C_{min}^{max})$. Since
$A=\frac{1}{2}(E_{1}(\epsilon)+E_{2}(\epsilon))$, $A$ is not an extreme point
of $\mathcal{U}(R_{min}^{max},C_{min}^{max})$.

Note that there are no other cases, as any such case would imply that the
matrix $A$ does not contain any fractional line, and therefore would belong to
$\mathcal{P}(R_{min}^{max},C_{min}^{max})$. This concludes the proof of lemma
5. ∎

The proof of theorem 2 a) is then directly the consequence of lemma 3, 4, and
5.

##  5 The polytope $\mathcal{U}^{k}(R_{min}^{max},C_{min}^{max})$

The set of matrices $\mathcal{U}^{k}(R_{min}^{max},C_{min}^{max})$ is the
subset of the polytope $\mathcal{U}(R_{min}^{max},C_{min}^{max})$, whose
members $A$ satisfy the additional constraint that the total sum of their
elements, $\sigma(A)$ is set to a given positive integer value $k$. Note that
this set is non empty only if $0<k\leq\sum_{i=1}^{n}R(i)$ and
$k\leq\sum_{j=1}^{n}C(j)$. Here we show that
$\mathcal{U}^{k}(R_{min}^{max},C_{min}^{max})$ forms a convex polytope whose
extreme values are the matrices in
$\mathcal{P}^{k}(R_{min}^{max},C_{min}^{max})$ using the same strategy we used
for $\mathcal{U}(R_{min}^{max},C_{min}^{max})$.

######  Lemma 6.

$\mathcal{U}^{k}(R_{min}^{max},C_{min}^{max})$ is convex.

###### Proof.

The proof is straight forward. Briefly, let $A$ and $B$ be two matrices
belonging to $\mathcal{U}^{k}(R_{min}^{max},C_{min}^{max})$, $\alpha$ a real
number in $[0,1]$, and $E=\alpha A+(1-\alpha)B$. To prove that
$E\in\mathcal{U}^{k}(R_{min}^{max},C_{min}^{max})$, we need to show that:

  * 1)

$E(i,j)\in[0,1]$: this comes from the convexity of $[0,1]$ (see proof of lemma
3).

  * 2)

$r(i)\leq\sum_{j=1}^{m}E(i,j)\leq R(i)$ for all $i\in[1,n]$: this comes from
the convexity of $[r(i),R(i)]$ (see proof of lemma 3).

  * 3)

$c(j)\leq\sum_{i=1}^{n}E(i,j)\leq C(j)$ for all $j\in[1,m]$: this comes from
the convexity of $[c(j),C(j)]$ (see proof of lemma 3).

  * 4)

$\sigma(E)=\sum_{i=1}^{n}\sum_{j=1}^{m}E(i,j)=k$: this comes from the
linearity of the operator $\sigma$.

∎

######  Lemma 7.

If $A$ belongs to $\mathcal{P}^{k}(R_{min}^{max},C_{min}^{max})$, then $A$ is
an extreme point of $\mathcal{U}^{k}(R_{min}^{max},C_{min}^{max})$.

###### Proof.

Again, the proof of lemma 7 is straightforward and identical to the proof of
lemma 4. ∎

Finally, we prove:

######  Lemma 8.

If a matrix $A$ is an extreme point of
$\mathcal{U}^{k}(R_{min}^{max},C_{min}^{max})$ then it belongs to
$\mathcal{P}^{k}(R_{min}^{max},C_{min}^{max})$.

###### Proof.

We use a proof by contrapositive with similar ideas as those used for the
equivalent proof for lemma 5. Let
$A\in\mathcal{U}^{k}(R_{min}^{max},C_{min}^{max})\smallsetminus\mathcal{P}^{k}(R_{min}^{max},C_{min}^{max})$.
As such, it has at least one fractional entry and therefore one fractional
line. We consider two cases.

  * Case 1)

All fractional lines include at least two fractional values. The proof is then
equivalent to the proof of Case 1 for lemma 5.

  * Case 2)

There exists at least one fractional line that includes a single fractional
value.

Let us assume that one such fractional line is a column $j_{1}$ (the proof
would be the same if it were a row) and let $A(i_{1},j_{1})$ be the only
fractional value on $j_{1}$. Just like in Case 2 of lemma 5, we proceed from
this value to generate a cycle of fractional values. If the cardinality of
this cycle is even, we are done. A problem arises, however, if the cardinality
of the cycle is odd:

$\displaystyle
C=\\{(i_{1},j_{1}),(i_{2},j_{2}),\ldots,(i_{2K-1},j_{2K-1})\\},$

where $2K-1=|C|$. We could still define the matrix $N$ as:

$\displaystyle\begin{cases}N(i,j)&=0\quad(i,j)\notin C\\\
N(i_{2k},j_{2k})&=1\quad k\in\\{1,\ldots,K-1\\}\\\
N(i_{2k-1},j_{2k-1})&=-1\quad k\in\\{1,\ldots,K\\}.\end{cases}$

However, if we define $E_{1}(\epsilon)=A+\epsilon N$ and
$E_{2}(\epsilon)=A-\epsilon N$, even with an $\epsilon$ small enough, we have
$\sigma(E_{1}(\epsilon))=\sigma(A)-\epsilon$ and
$\sigma(E_{2}(\epsilon))=\sigma(A)+\epsilon$, i.e. $E_{1}(\epsilon)$ and
$E_{2}(\epsilon)$ do not belong to
$\mathcal{U}^{k}(R_{min}^{max},C_{min}^{max})$.

To address this issue, let us first define that a row $i$ (resp. a column $j$)
is mutable if and only if $r(i)<RA(i)<R(i)$ (resp. $c(j)<CA(j)<C(j)$) where
$RA(i)$ is the sum of the elements of row $i$ and $CA(j)$ is the sum of the
elements of column $j$. In other word, a row or column is mutable if it is
possible to add or remove a well chosen small non-zero value to one of its
elements without violating the constraints on its sum. For example, the column
$j_{1}$ of $A$ defined above is mutable: it contains a single fractional
value, hence $CA(j_{1})$ is fractional and therefore satisfies the definition
of mutable. We now make the following claim:

######  Claim 3.

If the matrix $A$ contains one mutable column, then it contains at least 2
mutable columns.

###### Proof.

Let $l$ be a mutable column of $A$. if all other columns $j$ of $A$ were non
mutable, then their column sums $CA(j)$ would all be integer values, and
$\sigma(A)=CA(l)+\sum_{j=1,j\neq l}^{m}CA(j)$ would then be fractional, which
is a contradiction with the fact that $\sigma(A)=k$. Therefore, there exists a
least another column of $A$ that is mutable. ∎

Based on claim 3, as $j_{1}$ is mutable, there exists another column
$j_{2}\neq j_{1}$ that is mutable. We consider then two cases:

    * a)

The cycle $C$ include $j_{2}$.

We truncate the cycle to the first vertex that belongs to the column $j_{2}$.
This truncated cycle, $C^{\prime}$, includes an even number of vertices. We
can then conclude as in _Case 1_ by considering this cycle.

    * b)

The cycle $C$ does not include $j_{2}$.

Let $A(i_{2},j_{2})$ be a fractional value on $j_{2}$. We build a second cycle
$D$ in $A$ starting from this value. If the cardinality of this cycle is even,
we can follow exactly the proof from Case 1. Otherwise, the cycle has an odd
number of elements; we write it as

$\displaystyle
D=\\{(i^{\prime}_{1},j^{\prime}_{1}),(i^{\prime}_{2},j^{\prime}_{2}),\ldots,(i^{\prime}_{2M-1},j^{\prime}_{2M-1})\\},$

where $2M-1=|D|$. We define the matrix $N^{\prime}$ as:

$\displaystyle\begin{cases}N^{\prime}(i,j)&=0\quad(i,j)\notin C\\\
N^{\prime}(i_{2k},j_{2k})&=-1\quad k\in\\{1,\ldots,M-1\\}\\\
N^{\prime}(i_{2k-1},j_{2k-1})&=1\quad k\in\\{1,\ldots,M\\}.\end{cases}$

For $\epsilon$ small enough (i.e. smaller than the characteristic $\epsilon$
for the two cycles), we define $E_{1}(\epsilon)=A+\epsilon N+\epsilon
N^{\prime}$ and $E_{2}(\epsilon)=A-\epsilon N-\epsilon N^{\prime}$ where the
matrix $N$ was previously defined in equation Case 2). With the right choice
of $\epsilon$, All values in $E_{1}(\epsilon)$ and $E_{2}(\epsilon)$ are in
the interval $[0,1]$. In addition, by construction, any row/column sum of
$E_{1}(\epsilon)$ and $E_{1}(\epsilon)$ are equal to to that of $A$, with the
exception of columns $j_{1}$, $j^{\prime}_{1}$ and rows $i_{2K-1}$, , and
$i^{\prime}_{2M-1}$. Those four lines, however, are mutable: as $j_{1}$,
$i_{2K-1}$, and $i^{\prime}_{2M-1}$ contain a single fractional value, there
are mutable, and $j^{\prime}_{1}=j_{2}$ was chosen as it is mutable. The
constraints on the sums of all rows and columns of $E_{1}(\epsilon)$ and
$E_{2}(\epsilon)$ are therefore satisfied. In addition,

$\displaystyle\sigma(E_{1}(\epsilon))$ $\displaystyle=$
$\displaystyle\sigma(A)+\epsilon\sigma{N}+\epsilon\sigma(N^{\prime})$
$\displaystyle=$ $\displaystyle k+\epsilon(1\times(K-1)+(-1)\times
K)+\epsilon(-1\times(M-1)+1\times M)$ $\displaystyle=$ $\displaystyle
k-\epsilon+\epsilon=k,$ $\displaystyle\sigma(E_{2}(\epsilon))$
$\displaystyle=$
$\displaystyle\sigma(A)-\epsilon\sigma{N}-\epsilon\sigma(N^{\prime})$
$\displaystyle=$ $\displaystyle k-\epsilon(1\times(K-1)+(-1)\times
K)-\epsilon(-1\times(M-1)+1\times M)$ $\displaystyle=$ $\displaystyle
k+\epsilon-\epsilon=k.$

Therefore, $E_{1}(\epsilon)$ and $E_{2}(\epsilon)$ belong to
$\mathcal{U}^{k}(R_{min}^{max},C_{min}^{max}$. Since
$A=\frac{1}{2}(E_{1}(\epsilon)+E_{2}(\epsilon))$, $A$ is not an extreme point
of $\mathcal{U}^{k}(R_{min}^{max},C_{min}^{max})$.

This concludes the proof of lemma 8. ∎

The proof of theorem 2 b) is then just the consequence of lemma 6, 7, and 8.

Acknowledgment. The work discussed here originated from a visit by P.K. at the
Institut de Physique Théorique, CEA Saclay, France, during the fall of 2022.
He thanks them for their hospitality and financial support.

## References

## References

  * Villani [2008] Cédric Villani.  _Optimal Transport: Old and New_.  Grundlehren der mathematischen Wissenschaften. Springer, 2008.
  * Peyré and Cuturi [2018] Gabriel Peyré and Marco Cuturi.  Computational optimal transport.  _arXiv e-prints_ , art. arXiv:1803.00567, 2018.
  * Cotar et al. [2013] Codina Cotar, Gero Friesecke, and Claudia Klüppelberg.  Density functional theory and optimal transportation with coulomb cost.  _Communications on Pure and Applied Mathematics_ , 66:548–599, 2013.
  * Santambrogio [2015] F. Santambrogio.  _Optimal Transport for applied mathematicians_.  progress in nonlinear differential equations and their applications. Springer, 2015.
  * Kantorovich [1942] L. Kantorovich.  On the transfer of masses.  _Dokl. Acad. Nauk. USSR_ , 37:7–8, 1942.
  * Monge [1781] G. Monge.  Mémoire sur la théorie des déblais et des remblais.  _Histoire de l’Académie Royale des Sciences et Mémoires de Mathématique et de Physique tirés des registres de cette Académie_ , 1784:666–704, 1781.
  * Burkard et al. [2009] R. Burkard, M. Dell’Amico, and S. Martello.  _Assignment problems_.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, United States, 2009.
  * Birkhoff [1946] Garrett Birkhoff.  Tres observaciones sobre el algebra lineal.  _Univ. Nac. Tucuman, Ser. A_ , pages 147–154, 1946.
  * Brualdi and Gibson [1976] Richard A Brualdi and Peter M Gibson.  Convex polyhedra of doubly stochastic matrices—iv.  _Linear Algebra and Its Applications_ , 15:153–172, 1976\.
  * Brualdi [1980] Richard A Brualdi.  Matrices of zeros and ones with fixed row and column sum vectors.  _Linear algebra and its applications_ , 33:159–231, 1980\.
  * Muchlis [1992] Ahmad Muchlis.  On the extremal points of a class of polytopes of matrices.  _Linear algebra and its applications_ , 162:369–383, 1992\.
  * Gill and Linusson [2009] Jonna Gill and Svante Linusson.  The k-assignment polytope.  _Discrete Optimization_ , 6:148–161, 2009.
  * Ziegler [2012] Günter M Ziegler.  _Lectures on polytopes_ , volume 152.  Springer Science & Business Media, 2012.
  * De Loera and Kim [2013] Jesús A De Loera and Edward D Kim.  Combinatorics and geometry of transportation polytopes: An update.  _Discrete geometry and algebraic combinatorics_ , 625:37–76, 2013.
  * Cao and Chen [2019] Lei Cao and Zhi Chen.  Partitions of the polytope of doubly substochastic matrices.  _Linear Algebra and its Applications_ , 563:98–122, 2019\.
  * Chen et al. [2021] Zhi Chen, Zelin Zhu, Jiawei Li, Lizhen Yang, and Lei Cao.  Extreme points of certain transportation polytopes with fixed total sums.  _The Electronic Journal of Linear Algebra_ , 37:256–271, 2021.
  * Von Neumann [1953] John Von Neumann.  A certain zero-sum two-person game equivalent to the optimal assignment problem.  _Contributions to the Theory of Games_ , 2:5–12, 1953.
  * Kosowsky and Yuille [1994] J. Kosowsky and A. Yuille.  The invisible hand algorithm: Solving the assignment problem with statistical physics.  _Neural Networks_ , 7:477–490, 1994.
  * Koehl and Orland [2021] P. Koehl and H. Orland.  Fast computation of exact solutions of generic and degenerate assignment problems.  _Phys. Rev. E_ , 103:042101, 2021.
  * Dell’Amico and Martello [1997] Mauro Dell’Amico and Silvano Martello.  The k-cardinality assignment problem.  _Discrete Appl. Math._ , 76:103–121, 1997.
  * Mirsky [1959] L. Mirsky.  On a convex set of matrices.  _Archiv der Mathematik_ , 10:88–92, 1959.
  * Čihák [1970] Pavel Čihák.  Remarks on doubly substochastic rectangular matrices.  _Commentationes Mathematicae Universitatis Carolinae_ , 11(2):385–391, 1970.
  * Mendelsohn and Dulmage [1958] Nathan S Mendelsohn and A Lloyd Dulmage.  The convex hull of sub-permutation matrices.  _Proceedings of the American Mathematical Society_ , 9:253–254, 1958.
  * Brualdi and Lee [1978] Richard A Brualdi and Geoffrey M Lee.  On the truncated assignment polytope.  _Linear Algebra and its Applications_ , 19:33–62, 1978\.
  * Jurkat and Ryser [1967] WB Jurkat and HJ Ryser.  Term ranks and permanents of nonnegative matrices.  _Journal of Algebra_ , 5:342–357, 1967.
  * Rothblum and Schneider [1989] Uriel G Rothblum and Hans Schneider.  Scalings of matrices which have prespecified row sums and column sums via optimization.  _Linear Algebra and its Applications_ , 114:737–764, 1989\.
  * Brualdi and Hwang [2005] Richard Brualdi and Suk-Geun Hwang.  A bruhat order for the class of (0, 1)-matrices with row sum vector r and column sum vector s.  _The Electronic Journal of Linear Algebra_ , 12:6–16, 2005\.
  * daFonseca and Mamede [2009] Carlos M daFonseca and Ricardo Mamede.  On (0, 1)-matrices with prescribed row and column sum vectors.  _Discrete mathematics_ , 309:2519–2527, 2009.
  * Brualdi [2006] Richard A Brualdi.  Algorithms for constructing (0, 1)-matrices with prescribed row and column sum vectors.  _Discrete Mathematics_ , 306:3054–3062, 2006.
  * Cavenagh [2013] Nicholas J Cavenagh.  Defining sets and critical sets in (0, 1)-matrices.  _Journal of Combinatorial Designs_ , 21:253–266, 2013.
  * Chen et al. [2016] Wei Chen, Yanfang Mo, Li Qiu, and Pravin Varaiya.  Constrained (0, 1)-matrix completion with a staircase of fixed zeros.  _Linear Algebra and its Applications_ , 510:171–185, 2016\.
  * Deng [2019] Lei Deng.  Convex set of doubly substochastic matrices.  In _2019 Ninth International Workshop on Signal Design and its Applications in Communications (IWSDA)_ , pages 1–4. IEEE, 2019.


# Quantum Bayes AI

Nick Polson111Nick Polson is Professor at Chicago Booth<EMAIL_ADDRESS>Vadim Sokolov is an Assistant Professor at Operations Research at George Mason
U<EMAIL_ADDRESS>Jianeng Xu is at Chicago Booth:
<EMAIL_ADDRESS>
Booth School of Business
University of Chicago
Vadim Sokolov
Department of Systems Engineering
and Operations Research
George Mason University
Jianeng Xu
Booth School of Business
University of Chicago

###### Abstract

Quantum Bayesian AI (Q-B) is an emerging field that levers the computational
gains available in Quantum computing. The promise is an exponential speed-up
in many Bayesian algorithms. Our goal is to apply these methods directly to
statistical and machine learning problems. We provide a duality between
classical and quantum probability for calculating of posterior quantities of
interest. Our framework unifies MCMC, Deep Learning and Quantum Learning
calculations from the viewpoint from von Neumann’s principle of quantum
measurement. Quantum embeddings and neural gates are also an important part of
data encoding and feature selection. There is a natural duality with well-
known kernel methods in statistical learning. We illustrate the behaviour of
quantum algorithms on two simple classification algorithms. Finally, we
conclude with directions for future research.

Key Words: Quantum Bayes, Quantum Deep Learning, Quantum Feature Extraction,
Machine Learning, Quantum Algorithms, Quantum Embedding, Quantum Learning,
Quantum Sensing, Quantum Entanglement, Quantum Encryption, Quantum MNIST.

##  1 Introduction

Quantum computers are devices that harness quantum mechanics to perform
computations in ways that classical computers cannot. Quantum computing
promises an exponential speed-up for Bayesian computation. Our goal is to
provide a review of existing quantum algorithms with an emphasis on their
applicability in Bayesian Learning. We view quantum embeddings as kernel
methods and as such as central to modern day machine learning problems.
Quantum state measurement is based on an idea of Von Neumann, (1932) where the
state of a physical system is used to measure quantities of interest,
Simulating the quantum system has counterparts in Bayesian learning such as
MCMC simulation and Deep Learning. For recent reviews of quantum computing and
its promise for big data analytics, see Wang, (2022); Hidary, (2021) and
Schuld and Killoran, (2019); Schuld, (2021).

Von Neumann’s principle of quantum measurement is central to the initial data
encoding problem. Feynman, (1986), in a seminal paper, discusses the
implications of quantum physics for mechanical computational methods. From a
statistical viewpoint, Smith, (1984) and Breiman, (2001) provide a framework
for Bayesian computation and black box neural networks for data science in the
21st century. Whilst the quantum mechanics underlying the construction of a
quantum computer is still in its infancy, there are many algorithms that are
available for speed-up of existing methods. In a recent experiment, Google AI
showed quantum supremacy in a stylized experiment.

Quantum embedding is an important part of the data encoding process and the
type of functions that are easy to evaluate with quantum simulation. First, we
create a density matrix of states of size $2^{b}\times 2^{b}$ from a $b$ Qubit
computer. This will form the equivalent of the kernel space that is
commonplace in learning methods. Machine learning methods essentially use
linear algebra algorithms on this augmented feature space and quantum
computers are directly applicable to such storage and calculations as they
typically provide an exponential speed up in such algorithms, e.g. principal
components, matrix inversion, Fourier transforms, inverse problem to a name a
few.

The rest of the paper is outlined as follows. Section 1.1 provides a
discussion of computability of functions (Church, Turing) and its
implementation in prediction and learning. Section 1.2 provides a review of
important quantum algorithms. Section 2 provides a review of quantum
probability ( Wang, (2011)) and simulation. Von Neumann’s approach to quantum
measurement is also described with particular emphasis on MCMC simulation,
Deep Learning and Quantum Simulation, Our goal is to provide a unified
overview of such methods as they pertain to quantum encoding. This provides a
unifying view for quantum implementation of machine learning algorithms and
provides an avenue for determining extensions. Section 3 describes quantum
Bayes AI methods including quantum embedding. Section 4 describes quantum
machine learning methods. Section 5 provides an example to quantum MNIST.
Finally, Section 6 concludes with directions for future research.

Much of the gains from quantum computing will be in regard to computability
and representation of multivariate functions. Deep learning achieves this via
a superposition of semi-affine functions whereas quantum computers use the
Bloch sphere and a variety of quantum embeddings.

###  1.1 Computability of Functions

The central question is simply the speed and computability of functions.
Church, (1936) defines the so-called $\lambda$-calculus where $\lambda$
denotes a function. Turing, (1937) famously wrote about computability of
functions (a.k.a. predictive rules). The promise of quantum computation and
physical measurement by simulation relies on the simple observation that it
only takes a few dozen Qubits to be able to calculate a large superposition of
functions. Compositions of functions (as opposes to additive structures) such
as deep learners are also naturally computable. Computer languages such as
Haskell are directly designed to perform such calculations. Polson and
Sokolov, (2020) provide a discussion of symbolic differentiation and
manipulation for deep learning. Gradients are also easily computed using
quantum simulation.

The fundamental problem is to reconstruct input-output representations. We
first represent them using quantum $0-1$ interaction terms. We represent it as
a mixture (data augmentation of $0-1$s) using quantum embedding, then simulate
the solution to Schrod̈inger’s equation to provide the von Neumann measurement
of the quantity of interest. Turing showed that you can represent an n dim
function (boolean) $\\{0,1\\}^{n}\rightarrow\\{0,1\\}$. In $\lambda$-calculus
and evaluate it as a composite function. Kolmogorov-Arnold theory discusses
the representation of multivariate functions. Qubit activation allows you to
represent Boolean functions and compute then quickly.

###  1.2 Quantum Algorithms

There are many quantum algorithms which include Grover’s algorithm, quantum
Bayesian Shor’s factorization algorithm which has implications in cryptography
and quantum versions of standard machine learning algorithms.

#### Quantum Cryptography: Shor’s Algorithm

Shor’s quantum factoring algorithm has implications in cryptography and the
strong Church-Turing thesis, see Wang and Liu, (2022) for further discussion.
For classical computer, finding all price factors of a given number is known
to be $O(\exp(n^{1/3}\log^{2/3}n))$. Shor (1994) showed theoretically that
quantum computer can factorize any number in $O(n^{2}\log n\log\log n)$
operations. This has implications for Quantum Encryption.

Harrow, (2020) shows how to use classical computers for data pre-processing
and quantum computers for sampling and optimisation required to build a
predictive model. This paper addresses the problem of loading a dataset onto a
quantum computer. Computationally, this task is very expensive, however a pre-
processing technique can be used to reduce the size of the training data set
required by a quantum computer. Practically, also shows how the Grover
algorithm can be used to solve several computational problems that arise in
machine learning using a reduced data set. The problem is formulated as the
following minimization problem

$\arg\min_{y\in Y}\sum_{x\in X}f(x,y)$

Here $Y$ indexes the set of candidate models, $X$ is the training data set and
$f$ is a loss function.

Grover’s Algorithm As input for Grover’s algorithm, suppose we have a function

$f:\\{0,1,\ldots,N-1\\}\to\\{0,1\\}.$

In the “unstructured database” analogy, the domain represent indices to a
database, and $f(x)=1$ if and only if the data that $x$ points satisfy the
search criterion. We additionally assume that only one index satisfies
$f(x^{\star})=1$, and we call this index $x^{\star}.$ Our goal is to identify
$x^{\star}$.

Harrow’s Algorithm. This applies Grover’s algorithm to the fundamental
estimator of Bayesian Machine learning, the MAP estimator, via the
optimisation problem,

${\rm argmin}_{y\in Y}\;\sum_{x\in X}f(x,y)+\lambda\phi(y)$

where $\phi(y)$ is a regularisation penalty and $f(\cdot,\cdot)$ is an
empirical loss function. See Peng et al., (2020) for further discussion.

Quantum Monte Carlo: Annealing and Tunneling This allows a speed up of
traditional Monte Carlo algorithms based on the notion of quantum tunneling.
As with all quantum measurement there is a preparation of states, denoted by
$|\psi\rangle$, then a quantum device is used for evolution of states allowing
the application of unitary transformation $U$. See Wang (2011, 2016, 2022) for
further discussion and algorithmic details. Quantum annealing utilizes the
physical process of a quantum system whose lowest energy, or equivalently, a
ground state of the dynamic system, gives a solution to the posed Monte Carlo
problem via the solution of the Schrödinger equation.

There are many applications of quantum algorithms. One area of application is
in quantum sensing.

#### Quantum Sensing

Degen et al., (2017) Quantum Sensing is an advanced sensor technology that
improves the accuracy of how we measure, navigate, study, explore, see, and
interact with the world around us by sensing changes in motion, and electric
and magnetic fields. The analyzed data is collected at the atomic level.

##  2 Classical Learning

Supervised statistical learning given data as a sequence of input-output
pairs. The goal is to determine a predictive rule for new input cases or to
describe with statistical uncertainty is described by predictive distribution
$p(y\mid x)$ of future observations.

How does one achieve such a goal given the empirical data of $N$ input-output
pairs $(y_{i},x_{i})_{i=1}^{N}$. We need to formulate an input-output map,
$y=f(x)$, the question is how to construct $f$?

The machine learning problem then is to find a good predictor rule from a
training dataset if input-output pairs $(y_{i},x_{i})_{i=1}^{N}$ of observed
data. The goal is to predict at a new high dimensional
$x_{i}=(x_{1i},\ldots,x_{pi})$. To achieve good generalisability we need to be
able to perform nonlinear dimension reduction and to find a suitable set of
features/factors. Hence, machine learning algorithms rely on interpolation.
Gaussian processes and deep learners provide two classes of interpolators that
have been shown to do well in many fields of application. We will provide
quantum speed-ups of Gaussian processes and a framework for quantum deep
learning. An area of future research, is whether quantum representations
provide other useful classes of interpolators.

Given a set of $(output,input)$ pairs $(y_{i},x_{i})_{i=1}^{N}$, the goal is
simply to find a mapping (a.k.a. data transformation), denoted by $\hat{f}$,
where

$y_{i}=f(x_{i})$

and then to find a generalised prediction rule, that holds $\forall
x\in\mathcal{X}$, denoted by

$\hat{y}(x)=\hat{f}(x)$

Here $\hat{f}$ is estimated from the ”training” dataset.

Bayesian predictive calculations start with a probabilistic model and map
(a.k.a. data generating process), denoted by $p(y|x)$. Under predictive mean
squared error (MSE), the optimal rule is then given by the predictive mean

$\hat{y}(x)=E(y\mid x)=\sum_{y\in Y}y\;p(y|x)$

The use of latent (hidden) variables, denoted by $z$, is also an important
feature of Bayesian thinking. They can be viewed as ”extending the art of the
conversation”, see Lindley, (1990). In quantum computing they can be viewed as
states of a physical system. With latent variables, $z$, we can write, the
optimal predictive rule as

$\hat{y}(x)=E\left[E(y\mid x,z)\right].$

This is simply the law of total probability. The key insight is that it is
easier to assess the inner expectations, the caveat being that we need to
evaluate large sums which are tailored made for quantum computers. This
includes the case of classification, where we need to evaluate classification
probabilities $p(y=1\mid x)$ via expectations of an indicator functions.

The predictive distribution is given by

$p(y|x)=\frac{p(x|y)p(y)}{p(x)}$

The introduction of latent states, $z$, allows the posterior to be broken into
simpler parts

$p(y|x)=\sum_{z\in Z}p(y|z,x)p(z|x)$

Again we require enumeration over the latent states $z$. These expectations
can be calculated via von Neumann’s quantum measurement simulation as a state
of a system. This is similar to the view of a particle filter where the
researcher simulates a stochastic process whose marginals over time as
realisations of the sequence of posterior distributions that one wishes to
calculate.

##  3 Quantum Bayes AI

We begin with a basic review of classical and quantum probability to contrast
the two approaches.

Classical Probability. In a simple coin tossing experiment, $y$ takes one of
two values $\\{0,1\\}$ with probability $p_{1},p_{0}$ where $p_{1}=1-p_{0}$.
Here $p_{0},p_{1}>0$ and $\lVert p\rVert_{1}=1$.

Quantum Probability. Quantum probability generalises this to the complex plane
with the $L^{2}$-norm for its complex coefficients. One can view this as a
data transformation to the Bloch sphere which is shown in Figure 1.

Figure 1: Bloch sphere

#### Qubit

Like a bit is a basic unit of information in traditional computing, The qubit
is a basic quantum unit and it encodes a two-state quantum system, such as
spin of the electron (up/down). We can measure a qubit (the simplest quantum
system) and obtain either $0$ with probability $\lVert a_{0}\rVert^{2}$ or $1$
with probability $\lVert a_{1}\rVert^{2}$. Using Dirac’s notation, besides
states $|0\rangle$ and $|1\rangle$, a qubit can also take states as their
super-positions, which are linear combinations of the form

$|\psi\rangle=a_{0}|0\rangle+a_{1}|1\rangle$

where $a_{0},a_{1}\in\mathbb{C}$ are called amplitudes such that $\lVert
a_{0}\rVert^{2}+\lVert a_{1}\rVert^{2}=1$. The notation
$|+\rangle=\frac{1}{2}|0\rangle+\frac{1}{2}|1\rangle$ is used for the halfway
state generated by superposition.

Qubits can evaluate boolean functions quickly particularly deep learning
representations based on superposition of affine functions. In general, a
system of $b$ qubits has $2^{b}$ computational basis states of the form
$|x_{1}x_{2}\ldots x_{b}\rangle$ where $x_{j}=0,1$ and $1\leq j\leq b$ that
generate a $2^{b}$-dim complex vector space and a superposition state that is
based on $2^{b}$-amplitudes. Hence, we only need a system of a few dozen
’qubits’ to provide a significant increase in computational speed as
$2^{52}\approx 10^{16}$. A system with $b=500$ qubits would have the ability
to model more than the number of atoms in the universe.

#### Quantum Entanglement.

Quantum entanglement is the physical phenomenon that occurs when a group of
particles are generated, interact, or share spatial proximity in a way such
that the quantum state of each particle of the group cannot be described
independently of the state of the others, including when the particles are
separated by a large distance. The topic of quantum entanglement is at the
heart of the disparity between classical and quantum physics: entanglement is
a primary feature of quantum mechanics lacking in classical mechanics.

Quantum Von Neumann Measurement. The idea is to ’measure’ the ensemble sum via
the simulation of the physical system. A quantum state is completely described
by density operators. A vector $|\psi\rangle\in\mathbb{H}$ corresponds to a
density operator $\rho=|\psi\rangle\langle\psi|$ the projection onto
$|\psi\rangle$. An ensemble state $\rho_{k}$ with probability $p_{k}$ of being
in vector $|\psi_{k}\rangle$ has density operator

$\rho_{k}=\sum_{k=1}^{K}p_{k}|\psi_{k}\rangle\langle\psi_{k}|$

Suppose that $\rho$ is generated by a $b$-qubit state, creating a $2^{b}\times
2^{b}$-Hermitian density matrix.

Quantum Features. From a data analytic and Bayesian learning perspective, one
can think of the mapping $\rho:x\rightarrow\rho(x)$ as a fixed feature map
providing a data encoding from $\mathcal{X}$ to density matrices. This quantum
embedding is the equivalent of feature selection in machine learning (see
Hoadley, 2000, Bhadra et al 2022, Polson, Sokolov, Xu (2022)) and is an
important part of which class of functions that can be represented by quantum
neural networks, The corresponding quantum kernel is

$k(x,x^{\prime})=tr(\rho(x)\rho(x^{\prime}))$

with the space of functions, for weight operator $W$ being of the form

$f_{W}(x)=tr(\rho(x)W)$

For example, Kernel ridge regression $y=f_{W}(x)+\epsilon$ will be
straightforward to calculate using quantum simulation and von Neumann
measurement.

Specifically, if $|\psi(t)\rangle$ is the state of the quantum system, then
$|\psi(t_{1})\rangle$ and $|\psi(t_{2})\rangle$ are related by

$|\psi(t_{2})\rangle=U(t_{1},t_{2})|\psi(t_{1})\rangle\;\;{\rm
with}\;\;U(t_{1},t_{2})=e^{iH(t_{2}-t_{1})}$

where $U$ is a unitary operator and $H$ is a Hamiltonian. In practice this
operator has to be discretized in simulation via Trotter’s (1959) formula. The
Trotter’s formula approximates $e^{-iH\delta}$ by $U_{\delta}$, which requires
only the evaluation of each $e^{-iH_{l}\delta}$, where

$U_{\delta}=[e^{-iH_{1}\delta/2},\ldots,e^{-iH_{L}\delta/2}][e^{-iH_{1}\delta/2},\ldots,e^{-iH_{L}\delta/2}],$

see also Wang, (2011).

With a quantum state prepared in a physical state, $\rho$, the observable $Y$
with values in $\\{y_{1},\ldots,y_{M}\\}$ has a probability distribution

$p_{\rho}(Y=y)=tr(\rho Q_{y})=\sum_{j=1}^{M}y_{j}p_{\rho}(Y=y_{j})$

The equivalence between everything can be thought of as an iterative map
(a.k.s. stochastic simulation)

###  3.1 General Paradigm

This allows us to construct a general paradigm for quantum measurement in
machine learning.

Quantum Simulation Here the initial complex state is subject to a unitary
transformation

$|\psi^{(t)}\rangle=U^{t}|\psi^{(0)}\rangle$

where $U=e^{-iH}$ is a unitary operator with $UU^{\dagger}=1$ and $H$ is the
Hamiltonian of the system.

In continuous-time we can argue as follows. The heat equation simulation is to
solve Schrödinger’s equation

$i\frac{\partial|\psi^{(t)}\rangle}{\partial t}=H|\psi^{(0)}\rangle$

starting from initial state $|\psi(0)\rangle$ with solution

$|\psi(t)\rangle=e^{-iHt}|\psi(0)\rangle$

The numerical evaluation of $e^{-iHt}$ is needed and this requires
discretization using Trotter’s formula (Trotter, 1959).

MCMC Simulation Here the initial probabilistic state $p^{(0)}$ is subject to a
stochastic transition matrix, $P$, that is time reversible and the evolution
of the system is given by

$p^{(t)}=P^{t}p^{(0)}$

In continuous-time $P$ can be generated by the infinitesimal generator (which
takes the place of the Hamiltonian in the quantum system)

Deep Learning Here the initial state $x^{(0)}$ is subject to an iterate map
(with $t$-layers)

$x^{(t)}=F^{t}x^{(0)}$

where the map $F$ is. a composition of ridge functions (see Polson, Sokolov
and Xu, 2022).

Hence, from the von Neumann measurement viewpoint, all there methods can be
though of as symbolic manipulation of linear operators leading to the physical
evaluation of the quantity of interest. The main advantage of quantum systems
is the exponential speed up as you can directly calculate the sums need for
Bayesian learning, Expectations and Predictive rules which we discuss in
Section 3. Deep learners have the advantage of representing multivariate
functions.

###  3.2 Quantum Data Encoding

As data grows in scale and complexity and algorithms such as deep learning
become more elaborate, computational techniques become ever more important.

#### Quantum Superposition

An intrinsic difference exists between bit representation in classical and
quantum systems. All classical algorithms prevent the simultaneous occurrence
of their states. The quantum system allows for simultaneous occurrence of
their states. The computational basis takes the form
$|x_{1},x_{2},\ldots,x_{b}\rangle$, where $x_{j}={0,1},~{}j=1,\ldots,b$. Any
quantum superposition state is a liner combination of $2^{b}$ possible base
states whose complex coefficients are called amplitudes. A $b$-qubit state
$|\psi\rangle$ may take u superposition state of the form

$|\psi\rangle=\sum_{x_{1},x_{2},\ldots,x_{b}}\alpha_{x_{1},x_{2},\ldots,x_{b}}|x_{1},x_{2},\ldots,x_{b}\rangle,\quad\mathrm{s.t.}\sum_{x_{1},x_{2},\ldots,x_{b}}|\alpha_{x_{1},x_{2},\ldots,x_{b}}|^{2}=1.$

#### Quantum Gates and Circuits.

There are many gates and circuits that can be used to transform quantum
probabilities on the Bloch sphere. For example, one of the simplest gate is a
quantum NOT gate maps $|0\rangle$ to $|1\rangle$ and $|1\rangle$ to
$|0\rangle$ and tranforms
$|\psi\rangle=\alpha_{0}|0\rangle+\alpha_{1}|1\rangle$ into
$|\psi\rangle=\alpha_{0}|1\rangle+\alpha_{1}|0\rangle$. The quantum NOT gate
can be represented by the Paule gate matrix

$\begin{pmatrix}1&0\\\ 0&1\end{pmatrix},$

which has the property of being a unitary operator in $\mathbb{C}^{2}$. Figure
2 shows other quantum gates.

Figure 2: Quantum Gates

Another useful data Encoding gate is the rotational gate. Here the encoding
gate $U_{\rho}(x)$ transforms data samples
$x=(x_{1},\ldots,x_{n})\in\mathbb{R}^{n}$ into a quantum state

$|x\rangle=U(x)|0\rangle=R_{x}(x_{1})\circledcirc
R_{x}(x_{2}),\ldots,R_{x}(x_{n})|0\rangle^{O\times n},$

where

$R_{x}(x)=e^{-i\phi\sigma_{x}/2}=\begin{pmatrix}\cos(\phi/2)&-i\sin(\phi/2)\\\
-i\sin(\phi/2)&\cos(\phi/2)\end{pmatrix}$

We will discuss the application of these later in the our discussion of
Quantum Neural Networks. Now we describe quantum entanglement.

##  4 Quantum Machine Learning

Schuld (2021) makes an important observation that, though quantum models maps
data into a high-dimensional feature space, they can be trained and selected
in a low-dimensional subspace. With the density matrices $\rho(x)$ given by
some embedding as feature vectors, quantum machine learning models are linear
in the feature space. And borrowing ideas from classical machine learning
kernel theory, similar results can be derived from the quantum kernel
$k(x,x^{\prime})=tr(\rho(x^{\prime})\rho(x))$. For example, the optimal models
that minimise the cost is linear expansion in terms of quantum kernel
functions:

$f(x)=\sum_{m=1}^{M}\alpha_{m}tr(\rho(x^{m})\rho(x))$

where $x^{m},m=1,...,M$ are training data. Hence, searching for the optimal
model based on a given dataset is essentially an $M$-dimensional optimisation
problem. Furthermore, by choosing a convex loss function, this help avoids the
barren plateaus in variational training and guarantees the optimality.

###  4.1 Quantum Kernels

Quantum kernel, the inner product of two feature vectors $\rho(x^{\prime})$
and $\rho(x)$, depends on the underlying data-encoding $\phi(x)$:

$k(x,x^{\prime})=tr(\rho(x^{\prime})\rho(x))=|\langle\phi(x^{\prime})|\phi(x)\rangle|^{2}.$

The simple amplitude encoding is defined as

$\phi(x)=|x\rangle\langle x|=\sum_{i,j=1}^{N}x_{i}x_{j}^{*}|i\rangle\langle
j|$

where $x$ is normalized such that $\|x\|^{2}=\sum_{i}x_{i}^{2}=1$ and $i,j$
are computational basis. The induced kernel is $k(x,x^{\prime})=|\langle
x^{\prime}|x\rangle|^{2}=|x^{\dagger}x^{\prime}|^{2}$.

Other interesting and common strategies including: basic encoding

$\phi(x)=|\langle i_{x}^{\prime}|i_{x}\rangle|^{2}$

which gives the Kronecker kernel $k(x,x^{\prime})=\delta_{x,x^{\prime}}$ where
the input binary string $x$ is represented by the integer
$i_{x}=\sum_{k=0}^{n-1}2^{k}x_{k}$, and coherent state encoding which gives a
Gaussian kernel. There is resemblance between quantum kernels and classical
machine learning kernels, and many of them have Fourier representation.

###  4.2 Quantum Models

A quantum model is defined as the expectation of the measurement $\mathcal{M}$
under the density matrix $\rho(x)$,

$f(x)=tr(\rho(x)\mathcal{M})$

and the measurement $\mathcal{M}$ can always be expressed as a linear
combination of $\rho(x^{k})$ where $x^{k}$ are from the data domain. More
specifically, there exist $\\{x^{k},k=1,2,...\\}$, such that

$tr(\rho(x)\mathcal{M})=tr(\rho(x)\mathcal{M}_{exp})$

for all $x$, where $\mathcal{M}_{exp}:=\sum_{k}\gamma_{k}\rho(x^{k})$. More
importantly, this indicates that any quantumn model can be represented as

$f(x)=tr\left(\rho(x)(\sum_{k}\gamma_{k}\rho(x^{k}))\right)=\sum_{k}\gamma_{k}k(x^{k},x)$

###  4.3 Optimizing Quantum Models

The representer theorem states that the optimal quantum models which minimize
the regularized empirical risk,

$\frac{1}{M}\sum_{m=1}^{M}L(y^{m},f(x^{m}))+\gamma\|f\|_{F}^{2}$

admits a representation of the form

$f(x)=\sum_{m=1}^{M}\alpha_{m}k(x^{m},x)$

Here $x^{m}$ are training samples and $F$ is the reproducing kernel Hilbert
space (RKHS) corresponding the encoding. Hence, the optimisation over the
regularized empirical risk is the same as finding the optimal measurement
$\mathcal{M}$. The vectorized $\mathcal{M}_{opt}$ is the familiar form

$|\mathcal{M}_{opt}\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}=\sum_{m}y^{m}\left(\sum_{m^{\prime}}|\rho(x^{m^{\prime}})\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}}\mathopen{\hbox{${\langle}$}\kern-1.94444pt\leavevmode\hbox{${\langle}$}}\rho(x^{m^{\prime}})|\right)^{-1}|\rho(x^{m})\mathclose{\hbox{${\rangle}$}\kern-1.94444pt\leavevmode\hbox{${\rangle}$}},$

This is analogous to the ordinary least square result for regression
coefficients as quantum models are linear in feature space.

###  4.4 Quantum FFT

The quantum Fourier transform is naturally calculated in a quantum measurement
system as it is already expressed in rotational form.

The QFT takes a “position” state $|x\rangle$ to the corresponding momentum
state $|p\rangle$ and is defined as follows

$\mathrm{QFT}_{q}(x)=\frac{1}{\sqrt{q}}\sum_{p=0}^{q-1}\exp(2\pi
ipx/q)|p\rangle,$

where $q$ is the dimension of the systems Hilbert space. See Weinstein et al.,
(2001) for further discussion. In terms of gates, the two but QFT corresponds
to the unitary operator

$\mathrm{QFT}_{4}=\frac{1}{2}\begin{pmatrix}1&-1&1&1\\\ 1&i&-1&-i\\\
1&-1&1&-1\\\ 1&-i&-1&i\end{pmatrix}$

Then we can simulate system forward in an iterative manner as described in our
general paradigm.

###  4.5 Quantum Linear Regression

In a high dimensional regression setting, we need to calculate
$\hat{\beta}=X^{\dagger}y$ where $X^{\dagger}$ is the Moore-Penrose pseudo-
inverse of $X$. The singular decomposition can be computed as
$X^{\dagger}=V\Sigma^{-1}U^{\dagger}$ and hence represented as a quantum
state, see Schuld et al., (2016). Then we need to calculate

$\bar{\beta}=\sum_{k=1}^{K}\sigma_{k}^{-1}v_{k}u_{k}^{T}y$

The optimal prediction rule under MSE, is given by

$\hat{y}(x)=\sum_{k=1}^{K}\sigma_{k}^{-1}x^{T}v_{k}u_{k}^{T}y$

We need to calculate a set of inner products with is easily computable using a
quantum computer which are naturally designed to calculate inner products. See
Harrow et al., (2009) for further discussion.

###  4.6 Quantum Gaussian Process Q-GP

A quantum algorithm for quantum algorithm for Gaussian process regression was
introduced in Zhao et al., (2019). The GP predictor is

$\hat{f}_{*}=k^{T}_{*}(K+\sigma^{2}I_{n})^{-1}y,$

depending on whether the goal is to compute the mean or variance, we choose
$|b\rangle=|y\rangle$ or $|b\rangle=|k_{*}\rangle$, where
$k_{*}=k(x_{*},x_{*})$ is the covariance of the target point with itself, with
the corresponding formula for the predictive variance. The quantum GP
algorithm simulates $K+\sigma^{2}I_{n}$ as a Hamiltonian acting on an input
state $|b\rangle$. Then performs quantum phase estimation to extract estimates
of the eigenvalues of $K_{\sigma}^{2}I_{n}$ and stores them as a weighted
superposition and then creates the desired prediction. Schuld, (2021) show how
to provide estimates of

$k_{*}^{T}\left(K+\sigma^{2}I_{n}\right)^{-1}y\;\;{\rm
and}\;\;k_{*}^{T}\left(K+\sigma^{2}I_{n}\right)^{-1}k_{*},$

written as two inner products.

###  4.7 Quantum Stochastic Gradient Descent (Q-SGD)

Gradients are much easier to calculate in quantum systems. This is due to the
following property of unitary transformations. To implement the gradient
descent, we approximate the derivative of the loss function by taking the
symmetric difference,

$\frac{df}{dx}(x)=\big{(}f(x+\epsilon)-f(x-\epsilon)\big{)}/(2\epsilon)+O(\epsilon^{2}).$
(1)

However the gradient can be analytically calculated when we use the following
unitary operators

$\exp\,(i\,\theta\,\Sigma)$ (2)

where $\Sigma$ is a generalized Pauli acting on a few qubits, that is,
$\Sigma$ is a tensor product of operators from the set
$\left\\{\sigma_{x},\sigma_{y},\sigma_{z}\right\\}$ acting on a few qubits.
The derivative with respect to $\theta$ gives an operator whose norm is
bounded by 1. Therefore, the gradient of the loss function with respect to
$\vec{\theta}$ is bounded by $L$, the number of parameters.

For this function the gradient is the generalised Pauli operator $\Sigma_{k}$

$\frac{d\text{loss}(\vec{\theta},z)}{d\theta_{k}}=2\operatorname{Im}\big{(}\left\langle
z,1\right|U_{1}^{\dagger}...U_{L}^{\dagger}Y_{n+1}U_{L}...U_{k+1}\Sigma_{k}U_{k}...U_{1}\left|z,1\right\rangle\big{)}$
(3)

Note that $Y_{n+1}$ and $\Sigma_{k}$ are both unitary operators. Define the
unitary operator

$\mathcal{U}(\vec{\theta})=U_{1}^{\dagger}...U_{L}^{\dagger}Y_{n+1}U_{L}...U_{k+1}\Sigma_{k}U_{k}...U_{1}$
(4)

No we can exress the gradient as

$\frac{d\text{loss}(\vec{\theta},z)}{d\theta_{k}}=2\operatorname{Im}\big{(}\left\langle
z,1\right|\mathcal{U}\left|z,1\right\rangle\big{)}.$ (5)

Here, we can think of $\mathcal{U}(\vec{\theta})$ as a quantum circuit
composed of $2L+2$ unitaries each of which depends on only a few qubits. The
quantum device allows $\mathcal{U}(\vec{\theta})$ to act on
$\left|z,1\right\rangle$. Using an auxiliary qubit we can measure the right-
hand side of (5). We start with

$\left|z,1\right\rangle\frac{1}{\sqrt{2}}\big{(}\left|0\right\rangle+\left|1\right\rangle\big{)}$
(6)

and act with $i\mathcal{U}(\vec{\theta})$ conditioned on the auxiliary qubit
being 1. This leads to

$\frac{1}{\sqrt{2}}\big{(}\left|z,1\right\rangle\left|0\right\rangle+i\mathcal{U}(\vec{\theta})\left|z,1\right\rangle\left|1\right\rangle\big{)}$
(7)

Performing a Hadamard on the auxiliary qubit gives

$\frac{1}{2}\big{(}\left|z,1\right\rangle+i\mathcal{U}(\vec{\theta})\left|z,1\right\rangle\left|0\right\rangle\big{)}+\frac{1}{2}\big{(}\left|z,1\right\rangle-i\mathcal{U}(\vec{\theta})\left|z,1\right\rangle\left|1\right\rangle\big{)}.$
(8)

Now measure the auxiliary qubit. The probability to get $0$ is

$\frac{1}{2}-\frac{1}{2}\operatorname{Im}\big{(}\left\langle
z,1\right|\mathcal{U}(\vec{\theta})\left|z,1\right\rangle\big{)}$ (9)

Thus, by making repeated measurements we can get a good estimate of the
imaginary part becomes an estimate of the $k$ﬂth component of the gradient.

##  5 Application: Quantum TensorFlow

We demonstrate an application of a quantum neural network to the problem of
image classification. It is assumed that input is an image that is represented
as a sequence of $n$ elements $z=(z_{1},\ldots,z_{n})$, with each $z_{i}$
being +1 or -1. The output $l(z)$ is also binary and takes values +1 and -1.
The basic building block of a quantum neural network is a unitary operator

$U_{a}(\theta).$

We assume that this operator acts on a subset of the qbits and is defined by
parameter $\theta$. Then the neural network is a composition of those unitary
operators.

$U(\theta)=U_{L}(\theta_{L})U_{L-1}(\theta_{L-1})\ldots U_{1}(\theta_{1}).$

Assume our output (readout) bit is 1, then the corresponding input $z$, we
construct a computational basis

$|z,1\rangle=|z_{1},\ldots,z_{n},1\rangle.$

Applying our neural network to the input $z$ gives the state

$U(\theta)|z,1\rangle.$

Then on the readout qubit, we measure a Pauli operator $\sigma$ which give us
+1 or -1. The goal is that this outcome matches the true label of the data
(e.g. image). Given that the outcome is uncertain, we can use multiple copies
of the outputs and then average them out.

#### Quantum Empirical Loss.

Then to train a quantum model, we use the following loss function

$\text{loss}(\vec{\theta},z)=1-l(z)\left\langle
z,1\right|U^{\dagger}(\vec{\,\theta}\,)Y_{n+1}U(\vec{\,\theta}\,)\left|z,1\right\rangle.$
(10)

Quantum TensorFlow: MNIST Now we illustrate how to use Tensorflow-Quantum
(Broughton et al.,, 2021) to demonstrate the quantum computing. We will
demonstrate how quantum computing works for the problem of binary
classification. Our framework described above has a binary readout bit, and we
train a model to discriminate between digits 3 and 6 in the MNIST dataset.

To represent the image as a seqeunce of +1 and -1 elements, we downscale MNIST
data set so that each image is 4$\times$4 and fits into a quantum computer.

    x_train_small = tf.image.resize(x_train, (4,4)).numpy()
    x_test_small = tf.image.resize(x_test, (4,4)).numpy()

Since, we introduce each pixel with a qbit, we make the data set binary. We do
it by simply “thresholding” the pixel values.

    THRESHOLD = 0.5
    x_train_bin = np.array(x_train_nocon > THRESHOLD, dtype=np.float32)
    x_test_bin = np.array(x_test_small > THRESHOLD, dtype=np.float32)

Figure 3 shows the original image of digit 6 and compares it with the
downscale image as well as binary image.

|  |
---|---|---
(a) Original Image | (b) Downscaled 4x4 Image | (c) Binary
Figure 3: Example of original sample and downscaled sample

Then we “flatten” the image into a vector and rotate each +1 qbit using an $X$
gate. We use the `cirq` library that allows to write, manipulate, and optimize
quantum circuits, and then running them on a simulated quantum computer.

    values = np.ndarray.flatten(image)
    qubits = cirq.GridQubit.rect(4, 4)
    circuit = cirq.Circuit()
    for i, value in enumerate(values):
            if value:
                circuit.append(cirq.X(qubits[i]))

Then we convert each image to TensorFlow Quantum circuits

    Ψx_train_tfcirc = tfq.convert_to_tensor(x_train_circ)Ψ

The figure below shows an example of a single quantum circuit layer.

Figure 4: Example of a Layer of a Quantum Neural Network

The resulting model predicts the label of an image rather well. In fact, the
accuracy is comparable to the accuracy of a classical neural network.

Figure 5: Comparison of prediction accuracy between classical and quantum
neural network

##  6 Discussion

Whilst quantum hardware in its infancy, the theoretical and quantum algorithms
are well-developed. Our goal is to show how to apply these to machine learning
and AI. In many cases, quantum von Neumann measurement provide an exponential
speedup in computations. We outline a unified framework for quantum machine
learning algorithms. We embed regression and classification methods, Bayesian
learning, MCMC simulations Hammersley and Handscomb, (1964) and deep learning
into quantum measurement systems. As with ML algorithms, quantum ML uses
kernel methods and feature extraction. Again storage and calculations of inner
products and matrix inversions provide exponential increases in speed.

Harnessing the power of quantum computing is one of the challenges in the
$21$st century.

## References

  * Breiman, (2001) Breiman, L. (2001).  Statistical Modeling: The Two Cultures (with comments and a rejoinder by the author).  Statistical Science, 16(3):199 – 231.
  * Broughton et al., (2021) Broughton, M., Verdon, G., McCourt, T., Martinez, A. J., Yoo, J. H., Isakov, S. V., Massey, P., Halavati, R., Niu, M. Y., Zlokapa, A., Peters, E., Lockwood, O., Skolik, A., Jerbi, S., Dunjko, V., Leib, M., Streif, M., Dollen, D. V., Chen, H., Cao, S., Wiersema, R., Huang, H.-Y., McClean, J. R., Babbush, R., Boixo, S., Bacon, D., Ho, A. K., Neven, H., and Mohseni, M. (2021).  Tensorflow quantum: A software framework for quantum machine learning.
  * Church, (1936) Church, A. (1936).  An unsolvable problem of elementary number theory.  American Journal of Mathematics, 58(2):345–363.
  * Degen et al., (2017) Degen, C. L., Reinhard, F., and Cappellaro, P. (2017).  Quantum sensing.  Rev. Mod. Phys., 89:035002.
  * Feynman, (1986) Feynman, R. P. (1986).  Quantum mechanical computers.  Found. Phys., 16(6):507–532.
  * Hammersley and Handscomb, (1964) Hammersley, J. and Handscomb, D. C. (1964).  Monte Carlo Methods.  Springer Science & Business Media.
  * Harrow, (2020) Harrow, A. W. (2020).  Small quantum computers and large classical data sets.  arXiv:2004.00026 [quant-ph].
  * Harrow et al., (2009) Harrow, A. W., Hassidim, A., and Lloyd, S. (2009).  Quantum Algorithm for Linear Systems of Equations.  Physical Review Letters, 103(15):150502.
  * Hidary, (2021) Hidary, J. D. (2021).  Quantum computing: an applied approach, volume 1.  Springer.
  * Lindley, (1990) Lindley, D. V. (1990).  The 1988 Wald Memorial Lectures: The Present Position in Bayesian Statistics.  Statistical Science, 5(1):44 – 65.
  * Peng et al., (2020) Peng, T., Harrow, A., Ozols, M., and Wu, X. (2020).  Simulating Large Quantum Circuits on a Small Quantum Computer.  Physical Review Letters, 125(15):150504.
  * Polson and Sokolov, (2020) Polson, N. and Sokolov, V. (2020).  Deep learning: Computational aspects.  Wiley Interdisciplinary Reviews: Computational Statistics, 12(5):e1500.
  * Schuld, (2021) Schuld, M. (2021).  Supervised quantum machine learning models are kernel methods.  arXiv:2101.11020 [quant-ph, stat].
  * Schuld and Killoran, (2019) Schuld, M. and Killoran, N. (2019).  Quantum machine learning in feature Hilbert spaces.  Physical Review Letters, 122(4):040504.
  * Schuld et al., (2016) Schuld, M., Sinayskiy, I., and Petruccione, F. (2016).  Prediction by linear regression on a quantum computer.  Physical Review A, 94(2):022342.
  * Smith, (1984) Smith, A. F. M. (1984).  Present position and potential developments: Some personal views: Bayesian statistics.  Journal of the Royal Statistical Society. Series A (General), 147(2):245–259.
  * Turing, (1937) Turing, A. M. (1937).  Computability and $\alpha$-definability.  Journal of Symbolic Logic, 2(4):153–163.
  * Von Neumann, (1932) Von Neumann, J. (1932).  Mathematische Grundlagen der Quantenmechanik.  Berlin, J. Springer.
  * Wang, (2011) Wang, Y. (2011).  Quantum monte carlo simulation.  The Annals of applied statistics, 5(2A):669–683.
  * Wang, (2022) Wang, Y. (2022).  When Quantum Computation Meets Data Science: Making Data Science Quantum.  Harvard Data Science Review, 4(1).  https://hdsr.mitpress.mit.edu/pub/kpn45eyx.
  * Wang and Liu, (2022) Wang, Y. and Liu, H. (2022).  Quantum computing in a statistical context.  Annual Review of Statistics and Its Application (to appear).
  * Weinstein et al., (2001) Weinstein, Y. S., Pravia, M., Fortunato, E., Lloyd, S., and Cory, D. G. (2001).  Implementation of the quantum fourier transform.  Physical review letters, 86(9):1889.
  * Zhao et al., (2019) Zhao, Z., Pozas-Kerstjens, A., Rebentrost, P., and Wittek, P. (2019).  Bayesian deep learning on a quantum computer.  Quantum Machine Intelligence, 1(1):41–51.