#  $n$-valued Coset Groups and Dynamics

Mikhail Kornev

###### Abstract

We obtain asymptotic and exact formulae of growth functions for some families
of $n$-valued coset groups. We also describe connections between the theory of
$n$-valued groups and Symbolic Dynamics.

###### Contents

  1. 1 Introduction
  2. 2 Preliminaries
  3. 3 Growth of $n$-valued groups
  4. 4 Connection with Symbolic Dynamics

##  1 Introduction

In 1971, S. P. Novikov and V. M. Buchstaber gave the construction, which was
predicted by the theory of characteristic classes. This construction describes
a multiplication, with a product of any pair of elements being a non-ordered
multiset of $n$ points. It led to the notion of $n$-valued groups, which was
given axiomatically and developed by V. M. Buchstaber. Many authors are
studying the theory of $n$-valued groups with applications in various areas of
Mathematics and Mathematical Physics. Since 1996, V. M. Buchstaber and A. P.
Veselov have been studying the applications of the theory of $n$-valued groups
to discrete dynamical systems [2]. In 2010, V. Dragović showed that the
associativity equation for some 2-valued group explains the Kovalevskaya top
integrability mechanism [5]. Recently, M. Chirkov obtained results about the
growth function of some $n$-valued group on complex numbers
$\operatorname{\mathbb{C}}$ [8]. He also established connections between the
compositions of integer numbers and Venn diagrams. See references for further
relevant results.

In Section 2, we briefly introduce the main constructions from the theory of
$n$-valued groups and fix some notations following [4].

In Section 3, we get exact and asymptotic formulae for the growth functions of
the $n$-vaued coset groups

$\operatorname{\mathbb{G}}_{\varphi}(\langle a,b\ |\
a^{m}=b^{m}=1\rangle)\text{ and
}\operatorname{\mathbb{G}}_{\varphi}\left(\langle a_{1},...,a_{s}\ |\
a_{1}^{2}=...=a_{s}^{2}=1\rangle\right),$

where $\varphi$ is the automorphism permuting cyclically the generators of
group presentations in each case. Surprisingly, this problem is related to
$n$-bonacci numbers, and we use the result of [9] to succeed.

In Section 4, we consider more closely the case
$\operatorname{\mathbb{G}}_{\varphi}(\langle a,b\ |\ a^{3}=b^{3}=1\rangle)$,
constructing a fruitful tree which enables us to establish an interesting
combinatorial observation in Proposition 6, and connections with Symbolic
Dynamics.

## Acknowledgments

The author is grateful to Dr. Victor Buchstaber for the inspiring atmosphere
and fruitful discussions during the preparation of this paper.

##  2 Preliminaries

Let us give some definitions and examples following paper [4].

######  Definition 1.

An _$n$ -valued multiplication_ on a set $X$ is a map

$\nu:X\times
X\to\operatorname{\mathrm{sym}}^{n}(X)\,:\,\nu(x,y)=x*y=[z_{1},z_{2},\dots,z_{n}],\;z_{k}=(x*y)_{k},$

where $\operatorname{\mathrm{sym}}^{n}(X):=X^{n}/\Sigma_{n}$ is the
$n$-symmetric power of $X$, or, equivalently, the set consisting of all
$n$-multisets $[x_{1},...,x_{n}]$ (or, sets for short), such that the
following conditions are satisfied:

  * •

Associativity. The $n^{2}$-sets

$\displaystyle[x*(y*z)_{1},x*(y*z)_{2},\dots,x*(y*z)_{n}],$
$\displaystyle[(x*y)_{1}*z,(x*y)_{2}*z,\dots,(x*y)_{n}*z]$

are equal for all $x,y,z\in X$.

  * •

_Unit._ There is an element $e\in X$, for which $e*x=x*e=[x,x,\dots,x]$ for
all $x\in X$.

  * •

_Inverse map._ There is a map $\mathrm{inv}\colon X\to X$ with the property

$e\in\mathrm{inv}(x)*x\text{ and }\;e\in x*\mathrm{inv}(x)\text{ for all }x\in
X.$

If a set $X$ has an $n$-valued multiplication, then $X$ is called an
$n$-valued group. Sometimes, we denote $n$-valued groups via blackboard bold,
like $\operatorname{\mathbb{G}}$.

######  Example 1.

Any 1-valued group is an ordinary group.

######  Example 2.

Consider the semigroup of nonnegative integers $\mathbb{Z}_{+}$. Define the
multiplication
$\nu\colon\mathbb{Z}_{+}\times\mathbb{Z}_{+}\to\operatorname{\mathrm{sym}}^{2}\mathbb{Z}_{+}$
by the formula $x*y=[x+y,|x-y|]$. The unit $e=0$. The inverse map
$\mathrm{inv}(x)=x$. For the associativity, one has to verify that the two
$4$-sets

$[x+y+z,|x-y-z|,x+|y-z|,|x-|y-z||]$

and

$[x+y+z,|x+y-z|,|x-y|+z,||x-y|-z|]$

are equal for all nonnegative integers $x,y,z$ (an exercise for the reader).

As in the case of ordinary groups, one can define homomorphisms of $n$-valued
groups. Hence, the class of all $n$-valued groups forms a category
$\operatorname{\mathit{n}-\mathcal{G}\mathit{rp}}$.

######  Definition 2.

A map $f:X\to Y$ of $n$-valued groups is called _a homomorphism_ if

  * •

$f(e_{X})=e_{Y}$.

  * •

$f(\mathrm{inv}_{X}(x))=\mathrm{inv}_{Y}(f(x))$ for all $x\in X$.

  * •

$\nu_{Y}(f(x),f(y))=(f)^{n}\nu_{X}(x,y)$ for all $x,y\in X$.

In particular, the following diagram commutes:

$\textstyle{X\times
X\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{\nu_{X}}$$\scriptstyle{f\times
f}$$\textstyle{\operatorname{\mathrm{sym}}^{n}(X)\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{Y\times
Y\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{\operatorname{\mathrm{sym}}^{n}(Y)}$

There is a way to construct many $n$-valued groups [4]. Given an ordinary
($1$-valued) group $G$ with the mutliplication $\nu_{0}$, the unit $e_{G}$,
and the inverse element $\operatorname{\mathrm{inv}}_{G}(u)=u^{-1}$, and given
a finite subgroup $H$ of automorphism group $\operatorname{\mathrm{Aut}}(G)$.
Consider the quotient space $X:=G/\varphi(H)$ of $G$ by the action of the
group $\operatorname{\mathrm{im}}(\varphi)$ with the quotient map $\pi:G\to X$
and $n$-valued multiplication

$\nu:X\times X\to\operatorname{\mathrm{sym}}^{n}(X)$

defined by the formula

$\nu(x,y)=[\pi(\nu_{0}(u,h(v)))\ |\ h\in H]$ (1)

where $u\in\pi^{-1}(x)$ and $v\in\pi^{-1}(y)$.

######  Theorem 1 (V. Buchstaber, [4]).

The multiplication $\nu$ from (1) defines an $n$-valued group structure on the
orbit space $X=G/\varphi(H)$ with the unit $e_{x}=\pi(e_{G})$ and the inverse
map
$\operatorname{\mathrm{inv}}_{X}(x)=\pi(\operatorname{\mathrm{inv}}_{G}(u))$,
where $u\in\pi^{-1}(x)$.

######  Definition 3.

The $n$-valued group $\operatorname{\mathbb{G}}_{\varphi}(G)$ from Theorem 1
is called the coset group.

######  Example 3.

Consider
$G=\operatorname{\mathbb{Z}}/2\ast\operatorname{\mathbb{Z}}/2\cong\langle a,b\
|\ a^{2}=b^{2}=e\rangle$ and the finite cyclic group of automorphisms
generated by the element $\varphi:a\mapsto b$ of order $2$. It is easy to see
that the underlying set of $X=G/\langle\varphi\rangle$ consists of elements

$u_{2n}=[(ab)^{n},(ba)^{n}],$ $u_{2n+1}=[b(ab)^{n},a(ba)^{n}]$

for every $n\geqslant 0$, and the 2-valued mutliplication in the coset group
$\operatorname{\mathbb{G}}_{\varphi}(G)$ being

$u_{k}\ast u_{\ell}=[u_{k+\ell},u_{|k-\ell|}].$

Particularly, the 2-valued groups $\operatorname{\mathbb{G}}_{\varphi}(G)$ and
$\operatorname{\mathbb{Z}}_{+}$ from Example 2 are isomorphic. It means that
$\operatorname{\mathbb{Z}}_{+}$ is a $2$-valued coset group (but not every
$n$-valued group is a coset group).

It is time to discuss the "bridge" between the theory of $n$-valued groups and
dynamical systems.

######  Definition 4.

An $n$-valued dynamics $T$ on a space (not necessarily an $n$-valued group)
$X$ is a map $T:X\to\operatorname{\mathrm{sym}}^{n}(X)$.

Thinking about $X$ as a state space, an $n$-valued dynamics $T$ defines
possible states $T(x)=[x_{1},...,x_{n}]$ at the moment $t+1$ as a state
function of $x$ at the moment $t$.

######  Example 4.

Consider a polynomial with respect to the variable $y$:

$F(x,y)=b_{0}(x)y^{n}+b_{1}(x)y^{n-1}+...+b_{n}(x),\
x,y\in\operatorname{\mathbb{C}},$

with $b_{i}:\operatorname{\mathbb{C}}\to\operatorname{\mathbb{C}}$ being
arbitrary functions. The equation $F(x,y)=0$ defines an $n$-valued dynamics

$T:\operatorname{\mathbb{C}}\to\operatorname{\mathrm{sym}}^{n}(\operatorname{\mathbb{C}}),\
x\mapsto[y_{1},...,y_{n}],$

where $[y_{1},...,y_{n}]$ is the $n$-set of roots of $F(x,y)=0$.

For an $n$-valued dynamics $T:X\to\operatorname{\mathrm{sym}}^{n}(X)$, there
is a fruitful characteristic measuring the growth of $T$.

######  Definition 5.

For any $x\in X$ and an $n$-valued dynamics
$T:X\to\operatorname{\mathrm{sym}}^{n}(X)$, the growth function
$\xi_{x}:\operatorname{\mathbb{N}}\to\operatorname{\mathbb{N}}$ sends $k$ to
the number of different points in

$\bigcup\limits_{i=0}^{k}T^{i}(x)=\\{\ast\\}\cup T(x)\cup T^{2}(x)\cup...\cup
T^{k}(x),$

where $\ast$ denotes a marked point in the underlying set of $X$.

In this context, there is a general problem of characterizing such $n$-valued
dynamics that functions $\xi_{x}$ have a polynomial growth (see Figure 1).

Figure 1: Polynomial and exponential growth

$n$-valued dynamical systems come, for example, from an action of $n$-valued
groups on a space:

######  Definition 6.

An action of $n$-valued group $\operatorname{\mathbb{G}}$ on a space (not
necessarily an $n$-valued group) $X$ is defined to be the map

$f:\operatorname{\mathbb{G}}\times X\to\operatorname{\mathrm{sym}}^{n}(X),\
f(g,x)=g\cdot x=[x_{1},...,x_{n}],$

such that

  * •

for any elements $g_{1}$, $g_{2}\in\operatorname{\mathbb{G}}$ and any $x\in X$
the following $n^{2}$-multisets coincide:

$g_{1}(g_{2}x)=[g_{1}x_{1},...,g_{1}x_{n}]$

and

$(g_{1}g_{2})x=[h_{1}x,...,h_{n}x],$

where $g_{2}x=[x_{1},...,x_{n}]$ and $g_{1}g_{2}=[h_{1},...,h_{n}]$,

  * •

$ex=[x,...,x]$ for any $x\in X$.

######  Example 5 (The left action of $\operatorname{\mathbb{G}}$ on itself).

In the setting of Definition 6, take $X$ and $f$ to be the group
$\operatorname{\mathbb{G}}$ and the operation on $\operatorname{\mathbb{G}}$,
respectively.

The next class of $n$-valued groups will be in the focus of this paper.

######  Definition 7.

An $n$-valued group $\operatorname{\mathbb{G}}$ is called cyclic if there is
an element $g\in\operatorname{\mathbb{G}}$, such that any element $h$ of the
underlying set of $\operatorname{\mathbb{G}}$ belongs to a multiset $g^{k}$
for some $k>0$.

##  3 Growth of $n$-valued groups

In this section, we consider the class of such finitely presented groups
$G=\langle f_{1},...,f_{n}\ |\ R\rangle$ that admit an automorphism
$\varphi\in\operatorname{\mathrm{Aut}}(G)$ cyclically permuting $n$ generators
$f_{1},...,f_{n}$. Let $\mathbb{G}_{\varphi}(G)$ denote the $n$-valued group
arising from the group $G$ and the automorphism
$\varphi\in\operatorname{\mathrm{Aut}}(G)$. Note that
$\operatorname{\mathbb{G}}_{\varphi}(G)$ is cyclic in the sense of Definition
7. We study the growth function $\xi_{g}(k)$ associated with the element
$g=[f_{1},...,f_{n}]$. In this setting, we assume the marked point from
Definition 5 to be the empty word $\Lambda$, or equivalently, the identity
element of $G$. For short, we will use the notation $\xi_{k}$ instead of
$\xi_{g}(k)$.

######  Example 6.

In Example 3, we have the linear growth $\xi_{k}=k+1$.

######  Proposition 1.

For the group $\mathbb{Z}/3\ast\mathbb{Z}/3=\langle a,b\ |\
a^{3}=b^{3}=1\rangle$ and the automorphism $a\mapsto b$ the corresponding
$2$-valued coset group has the growth function

$\xi_{k}=F_{k+3}=\cfrac{1}{\sqrt{5}}\left(\left(\cfrac{1+\sqrt{5}}{2}\right)^{k+3}-\left(\cfrac{1-\sqrt{5}}{2}\right)^{k+3}\right)-1.$

In particular, the growth is exponential:
$\xi_{k}\sim\cfrac{\varphi^{k+3}}{\sqrt{5}}$ for $k\to\infty$, where
$\varphi=\cfrac{1+\sqrt{5}}{2}$.

###### Proof.

We will count the number of words of length $k$ at the step at which
$[a,b]^{k}$ appears. At step zero, we have an empty word corresponding to the
unit of the 2-valued group. The first step produces the element $[a,b]$. At
the second step, we have: $[a,b]^{2}=[[a^{2},b^{2}],[ab,ba]]$ — two elements.

At step $k$, all words in the normal form of the group
$\mathbb{Z}/3\ast\mathbb{Z}/3$ will appear. Any such word has the form (up to
automorphism)

$a^{k_{1}}b^{k_{2}}a^{k_{3}}...,$

where $k_{i}\in\\{1,2\\}$. It means there will be as many such words of length
$k$ as there are in how many ways the number $k$ can be divided into summands
from the set $\\{1,2\\}$, taking into account the order. Let us denote the
last quantity by $S_{k}$.

From step $k$, you can go to step $k-1$, taking $1$ as the first summand among
$k_{i}$’s — in this case, the remaining summands can be chosen in $S_{k-1}$
ways, or you can take $2$ as the first summand — in this case, the ones can be
chosen in $S_{k-2}$ ways. It means we have the recurrence relation

$S_{k}=S_{k-1}+S_{k-2},S_{0}=1,S_{1}=1.$

Recall that the recurrent sequence

$F_{n}=F_{n-1}+F_{n-2}$

with initial conditions

$F_{0}=0,F_{1}=1,$

is called the Fibonacci sequence. There is the well-known Binet formula

$F_{n}=\frac{1}{\sqrt{5}}\left(\left(\frac{1+\sqrt{5}}{2}\right)^{n}-\left(\frac{1-\sqrt{5}}{2}\right)^{n}\right).$

In our case,

$S_{k}=F_{k+1}.$

It is easy to prove by induction that

$F_{1}+...+F_{k}=F_{k+2}-1.$

Hence

$\xi_{k}=S_{0}+S_{1}+...+S_{k}=F_{k+3}-1=\cfrac{1}{\sqrt{5}}\left(\left(\cfrac{1+\sqrt{5}}{2}\right)^{k+3}-\left(\cfrac{1-\sqrt{5}}{2}\right)^{k+3}\right)-1.$

∎

######  Proposition 2.

For the group $\mathbb{Z}/2\ast\mathbb{Z}/2\ast\mathbb{Z}/2=\langle a,b,c\ |\
a^{2}=b^{2}=c^{2}=1\rangle$ with the cyclic automorphism $a\mapsto b\mapsto
c$, the corresponding $3$-valued group has the growth function
$\xi_{k}=2^{k}$. In particular, the growth is exponential.

###### Proof.

The normal form, up to automorphism, for the group from the conditions has the
form in which the first letter of the word is the letter $a$, and identical
letters are not adjacent. Let us prove by induction that the number $S_{k}$ of
distinct words in the normal form of length $k$ in the multiset $[a,b,c]^{k}$
is equal to $2^{k-1}$ for $k>0$, and $S_{0}=1$. When $k=0$, we have the empty
word. Suppose that at step $k-1$, we have $2^{k-2}$ different words of length
$k-1$. Then, we add one letter to the right of each word of length $k-1$ so
that the length of the normal form of the word increases by $1$. It can be
done in two ways for one fixed word. There are exactly $2^{k-2}$ words of
length $k-1$, so there will be $2\cdot 2^{k-2}=2^{k-1}$ words of length $k$.
The statement about $\xi_{k}$ now follows easily.

∎

######  Proposition 3.

For the group $\left(\mathbb{Z}/2\right)^{\ast s}=\langle a_{1},...,a_{s}\ |\
a_{1}^{2}=...=a_{s}^{2}=1\rangle$ with the cyclic automorphism
$\varphi:a_{i}\mapsto a_{i+1}$ $($the indices are taken modulo $s)$ the
corresponding s-valued coset group has the following growth function

$\xi_{k}=\begin{cases}\dfrac{\left(s-1\right)^{k}-1}{s-2}+1,&s\geqslant 3\\\
k+1,&s=2\end{cases}$

In particular, the growth is polynomial if and only if $s=2$.

###### Proof.

Similar to the proof of Proposition 2. ∎

Before moving on to the next example, let us define a sequence of $n$-bonacci
numbers.

######  Definition 8.

A sequence of integers $\\{F_{k}^{(n)}\\}$ is called $n$-bonacci if it
satisfies the recurrence relation

$F_{k}^{(n)}=F_{k-1}^{(n)}+...+F_{k-n}^{(n)}$ (2)

and initial conditions $F_{0}=...=F_{n-2}=0,\ F_{n-1}=1.$

In what follows, instead of the cumbersome notation $F_{n}^{(n)}$, we will
write $F_{n}$.

Using the standard recurrent sequence technique, we can obtain the following

######  Proposition 4 (see, for example, in [9]).

The explicit formula for the $k$-th term of the $n$-bonacci sequence
$($generalized Binet’s formula$)$ is as follows

$F_{k}=\sum\limits_{i=1}^{n}\frac{\lambda_{i}-1}{(n+1)\lambda_{i}-2n}\lambda_{i}^{k-n+1},$
(3)

where $\\{\lambda_{i}\ |\ 1\leqslant i\leqslant n\\}$ is the set of roots of
the characteristic polynomial $\chi(\lambda)=\lambda^{n}-\lambda^{n-1}-...-1$.

###### Proof.

From the theory of linear recurrences with constant coefficients, it is known
that the general solution of equation (2) has the form

$F_{k}=C_{1}\lambda_{1}^{k}+...+C_{n}\lambda_{n}^{k},$

where $\lambda_{i}$ are the same as in the conditions. It remains to find the
constants $C_{i}$ from the linear system of equations

$\begin{cases}C_{1}+...+C_{n}&=0,\\\
C_{1}\lambda_{1}+...+C_{n}\lambda_{n}&=0,\\\ ...\\\
C_{1}\lambda_{1}^{k-2}+...+C_{n}\lambda_{n}^{k-2}&=0,\\\
C_{1}\lambda_{1}^{k-1}+...+C_{n}\lambda_{n}^{k-1}&=1\end{cases}$

It can be easily done using Cramer’s rule and the properties of Vandermonde’s
determinant. The system has a unique solution since the characteristic
polynomial $\chi(\lambda)$ has no multiple roots. Due to the symmetry of the
system for cyclic permutations, we restrict ourselves to finding $C_{1}$. We
have

$C_{1}=\frac{W_{1}}{W}=(-1)^{n+1}\frac{\prod\limits_{1<i<j\leqslant
n}(\lambda_{j}-\lambda_{i})}{\prod\limits_{1\leqslant i<j\leqslant
n}(\lambda_{j}-\lambda_{i})}=\frac{(-1)^{n+1}}{\prod\limits_{1<j\leqslant
n}(\lambda_{j}-\lambda_{1})}=\frac{1}{\chi^{\prime}(\lambda_{1})},$

where $W$ is the Vandermonde determinant resulting from the system above,
$W_{1}$ is the determinant of the system with the first column replaced by a
column of free terms of the system (as usual, in Cramer’s rule).

It remains to find the value of $\chi^{\prime}(\lambda_{1})$. To do this,
consider the polynomial

$P(\lambda)=(\lambda-1)\chi(\lambda)=\lambda^{n+1}-2\lambda^{n}+1.$

It is clear that

$P^{\prime}(\lambda_{1})=(\lambda_{1}-1)\chi^{\prime}(\lambda_{1}).$

On the other hand,

$P^{\prime}(\lambda_{1})=(n+1)\lambda_{1}^{n}-2n\lambda_{1}^{n-1}.$

Hence, taking into account the fact that

$\lambda_{1}^{n-1}=\frac{1}{\lambda_{1}(2-\lambda_{1})}=\frac{1}{\lambda_{1}\cdot\lambda_{1}^{-n}}=\frac{1}{\lambda_{1}^{-n+1}},\text{
(since $2-\lambda_{1}=\lambda_{1}^{-n})$}$

it follows that

$C_{1}=\frac{1}{\chi^{\prime}(\lambda_{1})}=\frac{\lambda_{1}-1}{(n+1)\lambda_{1}-2n}\lambda_{1}^{-n+1}.$

∎

From Rouché’s theorem and Descartes’ rule of signs, it easily follows that the
polynomial

$\chi(\lambda)=\frac{\lambda^{n+1}-2\lambda^{n}+1}{\lambda-1}$

from Proposition 4 ($\lambda\neq 1$) has exactly one positive root greater
than $1$ and $n$ complex roots inside the unit circle centered at zero on the
complex plane. It means that among the terms of sum (3), there is precisely
one whose absolute value is greater than $1$. It turns out that the sum of all
other terms is always less than $1/2$ (see [12] for the connection of this
fact with random walks on the number line). Hence, the following holds:

######  Theorem 2 (Z. Du, G. P. Dresden’14, [9]).

For the $k$th term of the $n$-bonacci sequence we have

$F_{k}^{(n)}=\mathrm{rnd}\,\left(\frac{r-1}{(n+1)r-2n}r^{k-n+1}\right),$

where $r$ is the positive root of the polynomial
$\chi(\lambda)=\lambda^{n}-\lambda^{n-1}-...-1$, and $\mathrm{rnd}\,(\cdot)$
is the nearest integer.

Apply this result to the case of interest to us.

######  Proposition 5.

For the group $\mathbb{Z}/m\ast\mathbb{Z}/m=\langle a,b\ |\
a^{m}=b^{m}=1\rangle$ with the cyclic automorphism $\varphi:a\mapsto b$ of
order $2$, the growth function

$\xi_{k}\sim\frac{r^{k+1}}{mr-2(m-1)},$

for $m\geqslant 3$, $k\to\infty$, and $r$ is the positive root of the
polynomial $\chi(\lambda)=\lambda^{n}-\lambda^{n-1}-...-1$. In particular, the
growth of the corresponding $2$-valued coset group is polynomial if and only
if $m=2$.

###### Proof.

The normal form in this case looks like

$a^{k_{1}}b^{k_{2}}a^{k_{3}}...,$

where $k_{i}\in\\{1,2,...,m-1\\}$. Suppose that $m\geqslant 3$. Let $S_{k}$
denote the number of distinct words of length $k$ in the multiset $[a,b]^{k}$.
Assume that

$S_{-(m-2)}=...=S_{-1}=0\text{ and }S_{0}=1.$

We have the recurrence relation

$S_{k}=S_{k-1}+S_{k-2}+...+S_{k-m+1}.$

This recurrence relation, together with the initial conditions, determines the
$(m-1)$-bonacci sequence shifted by $m-2$

$S_{k}^{(m-1)}=F_{k+m-2}^{(m-1)},$

where $k\geqslant-(m-2).$

Thus, by Theorem 2

$S_{k}=\mathrm{rnd}\,\left(\frac{r-1}{mr-2(m-1)}r^{k}\right).$

Since $r>1$ (see the discussion after the proof of Proposition 4), we get:

$\xi_{k}=S_{1}+S_{2}+...+S_{k}\sim\sum\limits_{i=1}^{k}\frac{r-1}{mr-2(m-1)}r^{i}\sim\frac{r^{k+1}}{mr-2(m-1)}$

for $k\to\infty$.

∎

##  4 Connection with Symbolic Dynamics

In this section, we more closely study the 2-valued group
$\mathbb{G}_{\varphi}\left(\mathbb{Z}/3\ast\mathbb{Z}/3\right)$ considered
above in the context of combinatorics on words (see, e.g., [10, 13]).

For short, let us give the following

######  Definition 9.

We will call cubeless any word in the normal form of the group
$\mathbb{Z}/3\ast\mathbb{Z}/3=\langle a,b\ |\ a^{3}=b^{3}=1\rangle$.

From the normal form of the group $\mathbb{Z}/3\ast\mathbb{Z}/3$, it is clear
that at step $k$ all possible cubeless words in the alphabet $\\{a,b\\}$
appear.

Let us construct an oriented tree $\Gamma$ (see Figure 2), the vertices of
which will contain elements of the 2-valued group $\mathbb{G}$ as follows. At
step $0$, we start with the vertex corresponding to the empty word $\Lambda$ —
this is the root of our tree. At step $1$, we form a vertex $[a,b]$ adjacent
to the root. At step $2$, we draw two edges from the last vertex, each
corresponding to assigning the letter $a$ or $b$ to the right. You will get
two vertices with words of length $2$: $[a^{2},b^{2}]$ and $[ab,ba]$. At step
$k$, we start with all cubeless words of length $k-1$ and draw $1$ or $2$
edges from each vertex, guided by the following principle: if a word ends with
the first degree of a letter, then exactly $2$ edges will depart from this
vertex, corresponding multiplication by $a$ or $b$; if a word ends with the
square of a letter, then precisely one edge corresponding to this additional
letter will depart from this vertex.

Figure 2:

In this tree, at level $k$, there are $F_{k+1}$ vertices (as we saw in
Proposition 1).

Interestingly, one can draw two paths on this graph, corresponding to the two
well-known infinite words in Symbolic Dynamics. To define them, introduce the
notion of morphism.

######  Definition 10.

Let $A$ and $B$ be alphabets. A morphism is a map $\operatorname{\mathcal{F}}$
between the sets of words $A^{\ast}$ and $B^{\ast}$ in the corresponding
alphabets satisfying

$\operatorname{\mathcal{F}}(xy)=\operatorname{\mathcal{F}}(x)\operatorname{\mathcal{F}}(y)$

for all words $x,y\in A^{\ast}$, i.e., $\operatorname{\mathcal{F}}$ is a
homomorphism of monoids.

In some cases, the sequence of words $\\{\operatorname{\mathcal{F}}^{n}(a)\\}$
has a limit $\operatorname{\mathcal{F}}^{\infty}(a)$ (see [10] for details).

######  Example 7 ([10]).

Consider a morphism

$\operatorname{\mathcal{F}}:\\{a,b\\}^{\ast}\to\\{a,b\\}^{\ast},\ a\mapsto
ab,\ b\mapsto a.$

The infinite Fibonacci word is defined to be
$\Phi:=\operatorname{\mathcal{F}}^{\infty}(a)$

$\Phi=abaababaabaababaababaabaababaaba...$

######  Example 8 ([10]).

A morphsim

$\operatorname{\mathcal{F}}:\\{a,b\\}^{\ast}\to\\{a,b\\}^{\ast},\ a\mapsto
ab,\ b\mapsto ba$

defines the _Thue-Morse sequence_
$\mathcal{T}:=\operatorname{\mathcal{F}}^{\infty}(a)$

$\mathcal{T}=abbabaabbaababbabaababbaabbabaab...$

Since the paths in the tree $\Gamma$ correspond to all possible cubeless
words, the Fibonacci word $\Phi$ and the Thue-Morse sequence $\mathcal{T}$
define some two paths. In Figure 3, the path corresponding to the Fibonacci
word is shown in purple, and the path corresponding to the Thue-Morse sequence
is shown in red (with the overlap along the first two edges starting from the
root $\Lambda$).

Figure 3:

Note the following property of the constructed tree $\Gamma$.

######  Lemma 1.

At the $k$-th level of the tree $\Gamma$, all cubeless words of length $k$
starting with the letter $a$ are located from top to bottom in the standard
lexicographic order.

###### Proof.

Induction on $k$. When $k=0$, we have an empty word, and the statement is
obvious. Consider the set $L_{k}$ of all vertices of level $k$. Let
$v_{1},v_{2}\in L_{k}$, and $v_{1}<v_{2}$ in the lexicographic order. Then the
inequality

$\mathrm{Child}(v_{1})<\mathrm{Child}(v_{2})$

will be satisfied between the children of the vertices $v_{1}$ and $v_{2}$,
respectively, i.e. $u_{1}<u_{2}$ for any $u_{1}\in\mathrm{Child}(v_{1})$ and
$u_{2}\in\mathrm{Child}(v_{2})$. It is so because additions occur on the
right, and words are compared on the left in the lexicographic order. ∎

######  Lemma 2.

For an arbitrary vertex $v$ of the tree $\Gamma$, consider the subtree
$\Gamma_{v}$ starting from this vertex. Then the numbers of vertices at each
level of the subtree $\Gamma_{v}$ form the Fibonacci sequence starting either
from the Fibonacci number $F_{1}$ or $F_{2}$ depending on whether the word
$w(v)$ $($with the first letter $a)$ at the vertex $v$ ends with the first or
second degree of a letter.

###### Proof.

Let us look at the last letter of the word $w(v)$ and forget about the
remaining letters. It is either the first power of a letter, say $a$, or the
second power, say $a^{2}$ (without loss of generality). In the first case, a
subtree $\Gamma_{[a,b]}$ will grow from $v$, and in the second case, the
original tree $\Gamma_{\Lambda}=\Gamma$. But for the original tree $\Gamma$,
the numbers of vertices at each level form the Fibonacci sequence by
constructing $\Gamma$, as noted in Proposition 1. ∎

From this lemma, we get

######  Proposition 6.

Consider the sequence

$\\{\Theta_{k}\\}=\\{\Psi,\Psi a,\Psi aa,\Psi aab,...\\}$

of subwords of the infinite cubeless word

$\Psi aabaabaab...=:\Psi(aab),$

where $\Psi$ is any cubeless word whose the last letter is different from the
letter $a$. Then the number $Q_{k}$ of cubeless words which are not smaller
than the word $\Theta_{k}$ in the standard lexicographic order satisfies the
recurrence relation

$Q_{k}=Q_{k-1}+Q_{k-2}.$

###### Proof.

The word $\Psi aabaabaab...=\Psi(aab)$ corresponds to an infinite path
$\gamma$ in the tree $\Gamma$. Each $\Theta_{k}$ bijectively corresponds to
some vertex on this path. The set of words $\bigcup\limits_{k}Q_{k}$ we are
interested in corresponds to all vertices lying under $\gamma$. But the last
set is the union of all the vertices of the subtrees. Now, the statement
follows from the summation of all recurrence relations for each subtree
according to Lemma 2. ∎

## References

  * [1] V. G. Bardakov, A. Y. Vesnin, A Generalization of Fibonacci Groups., Algebra and Logic 42, 73–91 (2003)
  * [2] V. M. Buchstaber, A. P. Veselov, Integrable correspondences and algebraic representation of multivalued groups, FIM, ETH, Zurich, 1995
  * [3] V. M. Buchstaber, E. G. Rees, Multivalued groups, their representations and Hopf algebras, Transform. Groups, 2, no. 4, Birkhäuser, Boston, 1997, 325–349
  * [4] V. M. Buchstaber, $n$-valued groups: theory and applications, Mosc. Math. J., 6:1 (2006), 57–84
  * [5] V. M. Buchshtaber, V. Dragovic, Two-valued groups, Kummer varieties and integrable billiards , arXiv: 1011.2716, 2010
  * [6] V. M. Buchstaber, A. P. Veselov, Conway topograph, $\mathrm{PGL}_{2}(\operatorname{\mathbb{Z}})$-dynamics and two-valued groups, Russian Math. Surveys, 74:3 (2019), 387–430
  * [7] V. M. Buchstaber, A. P. Veselov, A. A. Gaifullin, Classification of involutive commutative two-valued groups, Russian Math. Surveys, 77:4 (2022), 651–727
  * [8] M. Chirkov, Growth function for an $n$-valued dynamics, arXiv: 2309.14251, 2023
  * [9] Z. Du, G. P. Dresden, A simplified binet formula for k-generalized fibonacci numbers, Journal of Integer Sequences, 17, 2014
  * [10] M. Lothaire, Algebraic Combinatorics on Words, Cambridge: Cambridge University Press; 2002, doi:10.1017/CBO9781107326019
  * [11] D. Johnson, J. Wamsley, D. Wright, The Fibonacci groups, Proc. London Math. Soc. (3) 29 (1974), 577–592
  * [12] User mathe, the forum discussion,
https://bbs.emath.ac.cn/thread-667-4-1.html#pid9145

  * [13] An. A. Muchnik, Yu. L. Pritykin, A. L. Semenov, Sequences close to periodic, Uspekhi Mat. Nauk, 64:5(389) (2009), 21–96; Russian Math. Surveys, 64:5 (2009), 805–871


# Improving the Efficiency of Payments Systems Using Quantum
Computing††thanks: The views expressed in the paper are solely those of the
authors and do not necessarily represent the views of the Bank of Canada and
other affiliations. The authors wish to thank Alex Koszegi, Tim Doyle, Peguy
Pierre-Louis, and Jonathan Garven from D-Wave Systems for their contributions
to the project. We also thank Maryam Haghighi for her inputs and Noorain
Noorani for his assistance to the project. We also thank Narayan Bulusu,
Rodney Garratt, Victoria Goliber, Isaiah Hull and Hossein Sadeghi for their
comments. †Corresponding authors<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
Christopher McMahon
GoodLabs Studio, University of Waterloo
&Donald McGillivray
GoodLabs Studio
&Ajit Desai${\ddagger}$
Bank of Canada
&Francisco Rivadeneyra
Bank of Canada
&Jean-Paul Lam${\ddagger}$
GoodLabs Studio, University of Waterloo
&Thomas Lo
GoodLabs Studio
&Danica Marsden
Bank of Canada
&Vladimir Skavysh
Bank of Canada

###### Abstract

High-value payment systems (HVPS) are typically liquidity-intensive as the
payment requests are indivisible and settled on a gross basis. Finding the
right order in which payments should be processed to maximize the liquidity
efficiency of these systems is an NP-hard combinatorial optimization problem,
which quantum algorithms may be able to tackle at meaningful scales. We
developed an algorithm and ran it on a hybrid quantum annealing solver to find
an ordering of payments that reduced the amount of system liquidity necessary
without substantially increasing payment delays. Despite the limitations in
size and speed of today’s quantum computers, our algorithm provided
quantifiable efficiency improvements when applied to the Canadian HVPS using a
30-day sample of transaction data. By reordering each batch of 70 payments as
they entered the queue, we achieved an average of C$240 million in daily
liquidity savings, with a settlement delay of approximately 90 seconds. For a
few days in the sample, the liquidity savings exceeded C$1 billion. This
algorithm could be incorporated as a centralized preprocessor into existing
HVPS without entailing a fundamental change to their risk management models.

_K_ eywords Quantum algorithm $\cdot$ Combinatorial optimization $\cdot$ NP-
hard problem $\cdot$ High-value payments system

JEL Codes: C61, D83, E42, E58

##  1 Introduction

High-value payments systems (HVPS), used to settle transactions between large
financial institutions, are part of the core financial infrastructure of every
country.111To get a sense of their importance, Canada’s HVPS processes payment
values equivalent to annual GDP every week. Central banks that operate or
oversee these systems are mandated to ensure their safety and efficiency. Most
HVPS around the world—including Canada’s Lynx, the US’s Fedwire, and the
Eurosystem’s TARGET2—are real-time gross settlement (RTGS) systems that settle
each payment request on an individual basis, i.e., without netting the
offsetting positions of participants. These systems are liquidity-intensive:
participants are required to have liquidity before payments are
processed.222In 2021, Canadian financial institutions participating in HVPS
allocated C$15-20 billion in daily liquidity to process around forty thousand
payments each day. Liquidity has an opportunity cost, incentivizing
participants to wait for incoming payments before making their own. This
behaviour reduces the efficiency of a system and can cause delays in payments
and even gridlock [1].

Given the systemic importance of HVPS, improving their efficiency is an
important area of research for central banks and academics alike [2, 3, 4].
Improving the liquidity efficiency of a HVPS can be done in two ways. The
first is to change the incentives of the participants for the timing of
payment submissions, so that payments are more coordinated. The second is,
once payments have been submitted to the system, change the order in which
payments will be settled into a sequence with lower liquidity demand.
Inherently, a re-ordering of payments requires some delay in payments.

In this paper we take the second approach. We propose a novel algorithm that
can reduce the liquidity requirements in HVPS without significantly increasing
the delay of payments. Our algorithm searches over the space of orderings of
payment requests. Reordering payments leads to an NP-hard combinatorial
optimization problem that grows to an enormous size with just a few
outstanding payments $(O(n!))$.333Our problem has _factorial_ time complexity
similar to that of the well-known traveling salesman problem, which is NP-hard
combinatorial optimization problem [5]. Quantum techniques, though not
expected to be able to solve such problems in polynomial time, are nonetheless
well suited to tackle this problem by providing speedup via heuristics leading
to near optimal solutions. Following recent technical advances, quantum
computing, especially when combined with classical computing in a hybrid
algorithm, is reaching a stage where it can provide superior solutions to some
of these problems when compared with classical-only approaches [6, 7, 8].

Economically, the proposed algorithm can be interpreted as a centralized
payments pre-processor. In this setup, participants submit their payment
requests to this algorithm before submitting them to the actual system.
Subsequently, after the optimization using quantum computing, our algorithm
returns to the participants a suggested ordering of payments that reduces the
total liquidity required to settle those payments relative to settling them in
the original order in which they were submitted. While participants could
effectively veto the proposal by not submitting their payments in the
suggested order, in practice, if the suggested re-orderings almost always
reduce aggregate liquidity demands, participants are unlikely to do so
strategically.444Participants chose their daily liquidity allocation using
risk management models attempting to cover a range of scenarios of liquidity
needs. In other words, participants’ choices are driven by expected daily
liquidity needs [9]. On an intraday basis, liquidity is rarely managed by
choosing specific payments to delay. Likewise, we think it is unlikely that
participants would evaluate the algorithm proposal batch by batch. The
proposed setup as pre-processing algorithm provides a simple way to
incorporate such algorithms into existing payment systems.

We formulated the reordering of payments in HVPS as a mixed binary
optimization (MBO) and translated it into a quadratic unconstrained binary
optimization (QUBO) format for processing on a generic quantum device.555MBO
problems arise in cases with discrete and continuous variables or in the cases
with inequality constraints; in our case, we have inequality constraints (see
Equation 17). Also, transformation to QUBO is necessary due to the binary
nature of quantum bits (or qubit) states and their physical coupling. See [8,
10, 11] for further details on QUBO formulation. Note the qubit is the
fundamental unit of quantum computing—the quantum version of the classic
binary bit. Due to the size of the problem (pertaining to the use of realistic
dataset) compared with the current quantum annealer topology (number and
connectivity of qubits), the final problem was posed as a constrained
quadratic model (CQM) and optimization performed on D-Wave Systems’
infrastructure using their hybrid quantum solver for CQMs [12]. The algorithm
was run on a sample of 30 days of non-urgent transactions from the Canadian
HVPS.666Note, we use randomly chosen days and settlement data from
LVTS—Canada’s HVPS before Lynx.

Our algorithm provided significant improvements in the payment system’s
efficiency by reducing the required liquidity to process payments with minimal
added settlement delay. The input for the optimizer is a given number of
queued payments (a "batch size"). While the set of potential re-orderings
grows with the batch size—potentially offering greater improvements in
liquidity efficiency—a larger batch size also requires a longer wait to
accumulate the payments into the batch and a longer processing time for the
algorithm. Among a variety of queue sizes, a batch size of 70 payments struck
a balance between optimizing liquidity savings and minimizing introduced
delays in settlements. With that batch size, we find liquidity savings in 26%
of them, providing an average daily savings of C$240 million. On average, the
added settlement delay was 90 seconds, mainly due to the time required to
accumulate the batch of 70 payments (the time required to run the algorithm,
in this case, is only 5 seconds). On half of the days in the sample, savings
exceeded C$107 million, and on three days, C$997 million. We also observed
that the total liquidity saved for each participant is proportional to the
participants’ total incoming and outgoing transaction values.

The batch size was constrained by the current size and speed limitations of
D-Wave’s CQM hybrid solver used in this analysis.777The CQM hybrid solver can
handle model with 100,000 constraints, 500,000 binary and integer variables
[12]. To test capabilities of different batch sizes, we ran the algorithm on
two typical days in our sample using batches of 140 payments. We achieved
significantly higher savings on those days (C$326 and C$94 million above the
savings of the 70-payment batch size) with an average delay of 3 minutes per-
batch (to accumulate the batch of 140 payments). This suggests that liquidity
savings could increase substantially as the capabilities of quantum computers
scale and that the more relevant limitation for this algorithm could be the
time required to accumulate larger batches of payments.888In our case, the
time required to accumulate payments grows linearly with batch size (see
Figure 3).

We proceed as follows. In Section 2, we discuss the related literature. We
discuss the methodology in Section 3 and present an overview of RTGS payment
systems and the schematic of the proposed algorithm as a centralized payment
preprocessor. In the same section, we present the payments queue optimization
algorithm and discuss the quantum annealing-based hybrid solver used to solve
this optimization problem. In Section 4 we provide a brief overview of the
sample of the data from Canada’s HVPS used to test the quantum algorithm. The
results and discussion follow this in Section 5 where we investigate the
performance of our algorithm. Finally, in Section 6 we conclude by discussing
the opportunities and challenges of the proposed algorithm and the scope for
future research. Appendix A provides further details on the formulation of the
objective function.

##  2 Literature Review

The problem of liquidity intensity of RTGS systems is well-known and
alternatives have been explored, in particular liquidity-saving mechanisms
(LSM) [2, 13, 14, 15]. The purpose of LSMs is to reduce the liquidity
requirements of a HVPS without the credit risks associated with delay net
settlement models999An alternative type of payment systems are use deferred
net settlement (DNS). In these systems, payments are submitted by participants
and accepted by the system but not settled in that moment. Instead, payment
exposures are accumulated and settled at given intervals after calculating the
net positions of the participants, therefore creating credit risk. See [16]
for an overview of LSMs.. LSMs temporarily accumulate payments requests in a
queue until potential offsetting positions are found and settlement occurs.
Theoretically, the efficiency improvements of LSMs depend on the liquidity
regime [17, 2, 18], and empirically, new evidence is emerging questioning
their effectiveness for improving efficiency of the HVPS [4].

In practice, LSMs are designed using heuristics and evaluated using simulation
approaches; therefore, they do not search the entire solution space to find
the optimal solution [19, 3, 20]. Compared to the LSM literature, our approach
is more direct and transparent: our algorithm attempts to evaluate the entire
space of re-orderings looking for liquidity savings.

A newer strand of literature is exploring an alternative approach by taking
seriously the incentives of the participants arising from the costs they face
for liquidity provision and the benefits of settling their payments. [21]
proposed a mechanism using the Shapley value cost allocation method by which
participants receive take-it-or-leave-it offers to contribute the liquidity
necessary to ensure that the welfare-maximizing netting proposals are always
accepted. The challenges of this approach are that it requires knowing the
valuation of the individual benefits of settling each payment and calculating
Shapley values becomes computationally intensive for large sets of payments.

Quantum computing technology is evolving rapidly, and there is growing
interest in its potential for solving certain classes of problems, such as
optimization, faster than classical computers [6]. Using the quantum
properties of superposition, entanglement, and tunneling, quantum computers
can search large input parameter spaces faster than any classical computer
could even theoretically achieve [22, 23]. Researchers have begun exploring
using quantum computing to solve various problems arising in
economics—including queuing problems that arise in financial markets [7, 8].
Closely related to our application, [24] extend an algorithm for a MBO problem
applied to transaction settlement. They optimize small batches of securities
trading transactions using gate-based quantum computing. However, their
optimization focuses on determining the number of transactions that can be
settled while allowing for netting, when the latter does not apply to RTGS
systems.

##  3 Methodology

In this section, we first provide an overview of a typical HVPS to understand
the incentives of participants and argue how a pre-processing algorithm could
be incorporated into the settlement process. Next, the rationale for choice of
computational resources is described. We then show how such an optimization
problem can be generically solved using quantum computing resources, alone or
in conjunction with a hybrid solver incorporating classical resources too. We
present the most general QUBO form of the objective function which could be
used by any quantum or classical solver, as well as the input format required
for the specific CQM hybrid annealing solver we opted to use due to current
physical limitations of quantum computing hardware. Lastly, we describe the
outputs of the CQM solver, and how this is converted to a final solution for
the pre-processing of payments in a given batch.

###  3.1 High-Value Payment Systems

In HVPS, the participating institutions seek to process the payment requests
received from their clients. To do so, they apportion collateral to the
central bank in exchange of the liquidity nececesary to settle those payment
requests. The system settles the payments in the order they arrive to the
system if they satisfy risk controls. The efficiency of these systems is
determined by the amount of liquidity that the participants choose to allocate
to the system and the timing of their payment requests. Participants can delay
submitting their payments to await for incoming payments to fund their own
payments reducing the need for their initial liquidity allocation. This gives
rise to the trade-off between liquidity and delay [1, 25]. Our algorithm seeks
to find better solutions to the liquidity management problem by optimizing the
order of payments without significantly increasing the settlement delay.

Figure 1 is a stylized schematic of a wholesale payment system. The figure
also shows a generic quantum optimizer as a central pre-processing mechanism
between participants and the payment system. A key challenge for new
algorithms seeking to improve the efficiency of existing payments systems is
how to incorporate them without incurring a fundamental change to the system.
As a pre-processor, quantum optimization algorithms could be used to propose
orderings that tangibly benefit the participants, incentivizing them to submit
the payments to the system in the order suggested by the algorithm.

Therefore, in our approach each participant submits its requests to a
centralized queue denoted in the ordered set $\\{P_{1},P_{2},\dots,P_{n}\\}$.
Each element in this set contains information about the value of the payment,
the two participants involved (payer and payee), and the time the request was
submitted. Once there are enough payments in the queue (a batch), the quantum
optimizer processes those payments and provides a proposal for an optimized
queue in the order $\\{Q_{1},Q_{2},\dots,Q_{n}\\}$ that minimizes liquidity
requirements and contains the new submission timing. If all banks accept the
proposal, those payments are submitted to the system in the new order proposed
by the quantum optimizer. As payments are settled throughout the day, each
participant’s liquidity position (balance) changes; this information—necessary
for the optimization process—is fed to the quantum optimizer by the payment
system.

Our approach is to have a fixed size $n$ of the batch of payments that will be
optimized by the quantum optimizer. Once the number of payments received
reaches that size the optimization is triggered. Since the flow of payment
requests varies though the day, the optimization occurred with varied
frequency (see Figure 3). To evaluate the algorithm we will report delay
statistics that consider the time needed to collect enough payments (wait-
time) for the batch and the time required by the combined compute resources -
hybrid of quantum and classical - to evaluate the optimized queue (processing
time). Given the wait time, we chose to test the algorithm on non-urgent
payments.101010Our approach with a fixed batch size is potentially not
suitable for urgent payments because it requires some wait time until a batch
of payments of a given size has accumulated. This can be alleviated to a
certain extent using flexible batch sizes. In the Canadian HVPS, however,
typically urgent payments represent only about 1.5% of the total volume of
payments, and they are typically very large in value and usually come at
specific times of the day.

Figure 1: Schematic of the quantum optimizer as a pre-processor to the payment
system. There are $N$ participating banks; they receive exogenous payments
demand from their clients’ throughout the day. Banks submit payments to the
central queue in the order $\\{P_{1},P_{2},\dots,P_{n}\\}$. The quantum
optimizer then processes those payments and proposes a new order
$\\{Q_{1},Q_{2},\dots,Q_{n}\\}$. If each bank accepts the proposal, those
payments will be submitted to the HVPS. Each participant’s initial and
intraday liquidity positions (updated after collateral apportionment at the
start of the day or after payment settlement throughout the day) are provided
to the quantum optimizer by the payment system.

###  3.2 Quantum Formulation of Queue Optimization Algorithm

The goal of the optimization is to minimize the aggregate liquidity needed in
the system to settle a given set of payments, ordered with index $i=1\ldots
n$. To measure the liquidity used in the system, we start by defining $N(i)$
as a participant’s net position before a payment $i$ is settled. A
participant’s net position is their liquidity balance accounting for all
incoming and outgoing payments up to that time in the day.

Next, we define $mNDP$ as the maximum net debit position experienced by a
participant at any time thus far in the day. If a participant begins the day
with $N=\$0$, then the end-of-day $mNDP$ after all payments have settled
measures the minimum liquidity a participant would have required to process
its payments without incurring in a negative position. We can therefore think
of $mNDP$ as equivalent to a participant’s minimum liquidity requirements.
Reducing the sum of $mNDP$ across all participants, holding payment values and
settlement delay fixed, represents an improvement in the liquidity efficiency
of the system. Further, let $b$ represent a participant’s liquidity
allocation.

The quantum optimizer pre-processor rearranges the queue of payment
instructions with an alternative index $t=1\ldots n$ that meets the goal of
minimizing total liquidity needs subject to the constraints that no
participant gets at any point into a negative liquidity position and that no
payment remains in the queue.

One can construct a Hamiltonian, in QUBO format, which aligns the minimum
energy of a quantum annealing system to the optimal queue order. While some
problems have a mapping already known in the literature, for example [10], our
problem did not. Thus we derived a formulation in Appendix A for a QUBO which
is the most general form of a problem which can be easily run on any type of
current quantum computer or even classical solvers111111We use $H$ to
represent the objective to be minimized over as in quantum annealing the
energy is being minimized. This is typically analogous to applying the
Hamiltonian operator to the state space. Note that this Hamiltonian is
slightly different from the one encountered in economic applications of
optimal control theory, which often contains a discount factor and hence need
transformation and the new operator is commonly refereed to as the _current-
value_ Hamiltonian [26, 27].

$\displaystyle H=$
$\displaystyle\sum_{\alpha}b_{\alpha}+\lambda_{1}\sum_{\alpha,t}\left(b_{\alpha}+N_{\alpha}(0)+{mNDP}_{\alpha}+\sum_{i}\sum_{\tau\leq
t}f(\alpha,i)\,x_{i,\tau}-s_{\alpha}(t)\right)^{2}$
$\displaystyle+\lambda_{2}\sum_{i}\left(1-\sum_{t}x_{i,t}\right)^{2}+\lambda_{2}\sum_{t}\left(1-\sum_{i}x_{i,t}\right)^{2},$
(1)

where:

$b_{\alpha}$ |  is the participant $\alpha$’s liquidity allocation. When minimized, it is the amount of additional liquidity that $\alpha$ would need to have to permit settlement of the payments in that batch (i.e., to avoid gridlock). Equivalently, it is the increase of $mNDP_{\alpha}$ after this batch is settled
---|---
$N_{\alpha}(0)$ |  is participant $\alpha$’s net position immediately before this batch is settled
$mNDP_{\alpha}$ |  is participant $\alpha$’s maximum net debit position incurred during previous batches earlier in the day
$x_{i,t}$ |  is a binary decision variable that indicates whether payment $i$ is settled ($x_{i,t}=1$) or not ($x_{i,t}=0$) at position $t$ in the final queue
$s_{\alpha}(t)$ |  is a non-negative real number acting as a slack variable that permits the enforcement of an inequality
$f(\alpha,i)$ |  is a function that returns $v$ (the value or amount of payment $i$) if $\alpha$ is the payee; $-v$ if $\alpha$ is the payer; and 0 otherwise
$\lambda_{1}$, $\lambda_{2}$ |  are Lagrange multipliers that enforce constraints as quadratic penalty terms

Minimizing $H$ over the set of variables
$\\{b_{\alpha},x_{i,t},s_{\alpha}(t)\\}$ $\forall\alpha,i,t$ finds the optimal
ordering $\\{x_{i,t}\\}$ that requires the least increase to the aggregate
$mNDP$. Note that the optimal ordering may be many-times degenerate.

###  3.3 Quantum Annealing and Hybrid Solvers

To solve an optimization problem, several tested commercial and open-source
optimization software packages represent the state-of-the-art. They implement
pre-solvers and heuristics to simplify the problem and find a solution (SCIP,
CPLEX, and Gurobi are some of the most popular, [28, 29, 30]). In addition,
they greatly simplify many technical issues in solving these optimizations,
especially the non-linear restrictions, and allow for greater focus on their
modelling. However, as stated earlier, for problems which scale with number of
variables as ours does, classical computing resources begin to fall short,
hence the desire to explore quantum methods.

A variety of different technologies are being explored as the best means to
build quantum hardware, but they are primarily classified into universal gate-
based quantum computers—that can process logic and run general algorithms—and
quantum annealers—used exclusively to solve optimization problems. While both
technologies can theoretically solve the optimization problems explored here,
currently available quantum annealers can solve much larger optimization
problems than gate-based quantum computers can, and their hybrid solver
counterparts larger problems still [12].

D-Wave Systems’ latest quantum annealer (aka QPU or quantum processing unit)
accessible via the cloud has over 5,000 qubits with 15-way qubit-to-qubit
connectivity (the “pegasus" topology). In addition to their quantum annealer,
D-Wave Systems offers access to 3 hybrid solvers which make use of both
classical and quantum compute resources –the Constrained Quadratic Model
(CQM), Binary Quadratic Model (BQM), and Discrete Quadratic Model (DQM). Each
of these solvers allows for different levels of flexibility in the
construction of the models, and numbers and types of variables, and numbers of
constraints.

In the case of the variable matrix mapping $i$ to $t$, the number of variables
to be considered here grows as the square of the batch size, quickly
overwhelming the size of the quantum annealer (QPU). Thus for our problem, as
the number of payments to be reordered increases, even 5,000 qubits with
15-way connectivity (which equates to around 180 fully connected qubits or
variables) are not enough to embed the problem in order to perform the anneal.
Since the CQM hybrid solver can handle more constraints (100,000) and binary
and integer variables (500,000) than any other available quantum hardware, we
opted to use it.121212The main downside our using CQM hybrid solvers’ that the
details of their inner workings are not known. On D-Wave Systems’ quantum
annealer, parameters such as annealing time, number of reads, and annealing
path can be controlled. For the hybrid solvers, these parameters are tuned
automatically, and the only input parameter (optional) is a time limit for the
calculation. Aside from the returned solutions, only the amount of QPU time
used per calculation can be retrieved.

The CQM solver takes objectives and constraints directly as inputs, thus the
problem was submitted specifically as:

Objective: $\displaystyle\text{min}\sum_{\alpha}b_{\alpha}$ (2) Subject to:
$\displaystyle b_{\alpha}+N_{\alpha}(0)+{mNDP}_{\alpha}+\sum_{i}\sum_{\tau\leq
t}f(\alpha,i)\,x_{i,\tau}\geq 0,\forall\alpha$ (3)
$\displaystyle\sum_{i}x_{i,t}=\lambda,\forall t$ (4)
$\displaystyle\sum_{t}x_{i,t}=\lambda,\forall i$ (5)

negating the need for slack variables and one Lagrange parameter on the user’s
end. The explanation of variables from Section 3.2 otherwise remains the same,
including that $x_{i,t}$ are binary decision variables.

The economic intuition of the problem is neatly summarized by the main
constraint of the problem (Equation 3). There, an increase in the first term
indicates a deterioration in the solution to the problem by requiring
additional liquidity in the system above the $mNDP$ observed up to the
processing of that batch. The last term in that constraint is the value of
settled payments, which is the product of the variable that indicates if a
payment has settled or not ($x_{i,t}$) and the function of the value of the
payment ($f$). Although this term changes discretely, it is useful to think in
marginal terms. A marginal increase in this term indicates that more payment
value is settled which if the constraint binds, can only be satisfied by
increase in liquidity, worsening the solution. Therefore, by searching over
different sets $\\{x_{i,t}\\}$, the optimizer explores the trade-off between
settling more payment value and increasing the liquidity necessary to do so.

Note that, using CQM solver each batch is processed as if it were happening in
real time, and only analyzed to yield aggregate results afterwards. That is,
we did not look ahead to use knowledge of what would come later in the day to
organize each earlier batch. Thus our method, while true to what could be
achieved in real world scenarios, optimizing batch by batch, doesn’t always
result in achieving the minimum liquidity requirements over the course of a
full day. This would require, among other things, predicting which payments
may appear in subsequent batches. As will be shown, on average, institutions
are nontheless still positively affected, proportional to the value of the
transactions they conduct.

###  3.4 From CQM Solver Results to the Optimal Queue Order

Finding a solution amounts to finding the matrix which represents the
reordering of the payments from $i$ to $t$, i.e. $\\{x_{i,t}\\}$. For example,
if there are seven payments in a queue, the CQM solver might return:

$\\{x_{i,t}\\}=\begin{pmatrix}0&1&0&0&0&0&0\\\ 0&0&0&0&1&0&0\\\
0&0&0&0&0&1&0\\\ 1&0&0&0&0&0&0\\\ 0&0&0&1&0&0&0\\\ 0&0&0&0&0&0&1\\\
0&0&1&0&0&0&0\end{pmatrix}$ (6)

which would mean that first the 4th payment in the original queue should be
processed, then the 1st, and so on. To find these solutions, our algorithm
takes the payments in the batch, generates the objective and constraints
imposed by each payment in that batch, and runs several times, returning a
reordering solution $\\{x_{i,t}\\}$ each time.

The reason to run multiple times is due to the nature of quantum computation.
Each run of the annealer potentially ends in a different state based on the
length of the anneal (analogous to “cooling time" in a physical anneal) and
the underlying probability distribution of the states (the lower the energy of
the final state, the higher its probability of occurring)131313The problem is
run on the QPU by first starting from the state where all qubits are in a
superposition of the basis states $\ket{0}$ and $\ket{1}$. This corresponds to
the lowest energy state of the initial, “tunneling” or “mixing” Hamiltonian.
This initial Hamiltonian is then evolved to a Hamiltonian encoding the problem
which is being solved (in this case given by Equation 1 \- enacted through
coupling parameters and biases on the qubits. If this change is slow enough,
the adiabatic theorem guarantees that the system will end up in the lowest
energy state of the final Hamiltonian, which corresponds to an optimal
solution of the problem, the number of which correspond to the degeneracy of
the ground state.. In practice, time constraints mean that we may reach the
near-optimal-but-not-best solution in a given run. However, more solutions
provide more confidence based on the statistics of the solution set that one
of the optimal reorderings was found (there may be multiple which are equally
good).

Our first step in post-processing the set of returned reordering solutions is
to check if any violated the constraints in the Hamiltonian, i.e., an
infeasible solution.141414For instance, if there are more than one non-zero
element in each row or column of the output matrix in Equation 6 then that is
an infeasible solution. This information is provided by the CQM solver. If so,
those infeasible solutions are discarded. The remaining feasible solutions are
then searched classically for the ones that provided the most savings,
yielding our final result and the payments processing order for that batch.

The optimality of the quantum solution is challenging to quantify for larger
problems since optimal solutions become too computationally expensive to
verify (and is the reason we turn to quantum computation in the first place).
However, an approximate measure of optimality can be inferred from the number
of feasible candidate solutions returned from the CQM solver. The CQM solver’s
default run time for batches of 70 payments, for example, provides about 50-70
solutions. However, as the problem gets larger with queue size, fewer feasible
solutions came back from CQM. This indicates that the quality of the solutions
went down. Thus if the CQM solver returned a histogram with many (few)
feasible results, it should correlate with a better (worse) solution. Indeed,
we observed this correlation (see Figure 5).

##  4 Data

To test the performance of the proposed quantum algorithm, we randomly sampled
30 days from Canada’s LVTS settlement data between Jan 2015 to Dec 2017. From
those 30 days, we used only non-urgent payment requests made between 8 am and
6 pm.151515In the LVTS, non-urgent payment requests are submitted to Tranche
2, and urgent payment requests are submitted to Tranche 1. The urgent payments
could have a high cost of delay; therefore, we exclude them from the sample.
The payment requests submitted after 6 pm and before 8 am are special
transactions; therefore, they are also excluded from the chosen sample (see
[31]). In our chosen sample, after excluding non-urgent payments, we have
approximately 23,000 transactions per business day. On some busy days,
however, the LVTS processed a higher volume of transactions.161616On average,
the LVTS processed approximately 37,000 transactions per business day.
However, on some busy days, such as the end of the month, quarter, or the next
day after a national holiday, the LVTS processed about 45,000 to 50,000
transactions. Therefore, to test the benefits of a quantum optimizer on such
days, we ensured that ten high volume days were included in our 30 days
sample.171717Before random sampling of 30 days, we exclude a few days from our
sample on which volume is minimal, for instance, on Canadian provincial
holidays or holidays in the US.

In our sample, the frequency at which payment requests are submitted varies
with the time of day, as shown in Figure 2 (left). We see a much higher volume
of payments settled between 8 to 9 am. The number of payment requests remains
steady from about 9 am until about 4 pm before they tapper down towards the
end of the day. On the other hand, the payment amount submitted to the system
varies almost uniformly throughout the day; However, the share of the number
of payments with high value is slightly higher towards the end of the day in
our sample, as shown in Figure 2 (right).

Figure 2: (Left) the number of payments and (right) payment amount submitted
to the system over time for each day in our sample.

##  5 Results and Discussion

This section presents the results of various simulation exercises performed on
the 30 days of sample from Canada’s LVTS settlement data using the D-Wave
system’s CQM hybrid solver. First, we present the results concerning the
computational time to accumulate and solve queue optimization problems for
different batch sizes. Next, we examine the solution quality of the quantum
optimizer. This is followed by the performance of the quantum algorithm on the
entire sample. Finally, we investigate the performance of the quantum
algorithm for a chosen day in our sample.

###  5.1 Computation Time

The reordering of a queue using the quantum optimizer consists of the
following steps:

  1. 1.

accumulating the batch of payments in the order they arrive,

  2. 2.

building the CQM object formulated in Section 3.3,

  3. 3.

sending the object to D-Wave Systems’ hybrid CQM solver to perform
optimization,

  4. 4.

post-processing the solver’s output to extract the solution.

For step 1, the wait time to collect a batch of payments grows linearly with
the batch size (Figure 3). However, for steps 2 and 3, the processing time and
memory required by the CQM object grow cubically with batch size. Figure 4
shows the time to compile and solve the CQM object versus queue size ranging
between 20 to 700 payments. The cubic growth is due to the scaling of the
number of biases which grows like $O(n^{3})$ in Equation 1. Finally, in step 4
the CQM solver’s output is processed to get the best solution, i.e., the order
in which submit the payments.

Figure 3: Average waiting time to fill a queue for a given queue size. The
orange star highlights the queue size of 70 payments. Figure 4: Time needed to
compile the CQM object (blue) and then find the solution using D-Wave Systems’
hybrid CQM solver (red), along with a best fit. Note that the same day’s
payments are used to compute the time, but truncated to the desired queue
size.

For real-world applications of such algorithms, solving the quantum reordering
should not take longer than collecting the batch of transactions in order to
avoid added delays. For a queue size ranging between 20 to 700 payments, the
average wait-time ranges from 1 to 15 minutes (Figure 3). However, this varies
with the time of day as the frequency at which payment requests are submitted
changes throughout the day (Figure 2, left). For a batch of 70 payments, the
average wait time is about 90 seconds and average compile and processing time
is 5 seconds each. The cubic growth in time shown in Figure 4 limits the
maximum queue size to around 140 transactions, which takes on average about
115 seconds to compile and solve the problem. However, there is room for
potential improvement by calculating the CQM on more powerful computers with
multiprocessing and reducing network latency by calculating the problem on
servers physically closer to the solvers. This could be explored in the
future.

####  5.1.1 Classical Comparison

To compare our optimization with a classical solver, we formulated our
algorithm to run on SCIP, a non-commercial solver for mixed-integer nonlinear
programming [28, 29]. It utilizes pre-solving and heuristics to divide the
problem into subproblems and tighten domain variables to solve recursively,
similar to CPLEX [30]. Unfortunately, depending on the set of transactions
submitted, the structure of the network graph, and the mNDP’s lower bound, it
generally takes over 2 hours to solve and occasionally fails after over 24
hours for a batch size of 70, run on a R5.4xlarge AWS instance (16 vCPUs at
3.1GHz, 128GB RAM). This suggests that the our problem is better suited to
being run on DWave’s hybrid quantum and classical solver.

###  5.2 Solution Quality

Using a representative set of 700 payments, Figure 5 shows the aggregate mNDP
versus batch size for single queues containing the first $n$ payments, with
$n$ ranging from 70 to 700, (black dots). The comparison to a FIFO queue
arrangement is also shown (blue line), as well as the number of feasible
solutions returned by the CQM solver (red bars). Up to and including a batch
size of 300, the solution from the CQM solver is better than FIFO; however,
the results become worse than FIFO as batch size continues to increase. This
corresponds to a drastic decrease in the number of feasible results returned.
For batch sizes above 250, the default time limit parameter for the
calculation (chosen algorithmically by the CQM solver based on multiple
properties of the CQM input) was not enough to obtain any feasible results and
had to be overwritten.

The data point corresponding to a batch size of 700 used 3 hours of hybrid
quantum computation time—set manually—to achieve a solution that is superior
to FIFO. These results suggest that we have found the maximum problem size
($n\lesssim 200$) that current hybrid classical-quantum annealing technology
can solve in a practically-reasonable amount of time. However, it also implies
that as the technology continues to develop and scale, future hybrid solvers
or quantum annealers with more entangled qubits which can search larger state
spaces will be able to handle larger queues.

Figure 5: The aggregate maximum net debit position (mNDP; left) and number of
feasible solutions (right) with respect to batch size. The mNDP of FIFO is
included for reference (blue line). Note that the hybrid CQM solver
computation time limit for batches of size 300, 400 and 500 was set to 5, 11
and 26 minutes, respectively (2 times the CQM default setting). For the batch
size of 700, a time limit of 3 hours was used (6 times the default).

###  5.3 Performance of Quantum Algorithm on the Entire Sample

As mentioned earlier, the optimization of settlement order was conducted by
dividing each day into batches and solving each batch on the CQM solver. The
batches were solved using each participant’s mNDP and net position from the
previous batch. The summary of the results for the entire sample is presented
in Table 1 and 2 for a batch size of $n=70$. Note that these tables represent
two different scenarios. For Table 1, we ran a horse race between the FIFO and
CQM. They both started with the same initial conditions (for instance, the
same mNDP at 8 am). Consequently, every batch throughout the day was the same
for both FIFO and CQM. At the end of the day (6 pm), we report the results
(mNDP). However, for Table 2, we ran the batch-by-batch horse race between
FIFO and CQM. In both FIFO and CQM, we started with the same initial
conditions for every batch and checked which one gave a better result, and
they were reported.

Table 1: Summary statistics of the sample data showing the total value of payments settled in Tranche 2 between 8am and 6pm. The days are sorted in descending order of total number of batches where each batch contains $n=70$ payments. A batch is considered optimizable if settlement by the original order in which they came, (i.e., FIFO), increases the aggregate mNDP. The improved batches column shows the number in which the quantum optimizer using the CQM solver was able to reduce the aggregate mNDP. The last column shows the total liquidity saved by the end of the day using the CQM solver. Date |  Value Settled ($Bn) |  Total Batches |  Improved Batches |  End-of-Day Savings ($Mn)
---|---|---|---|---
2017-06-30 |  158.00 |  492 |  89 |  70.76
2017-12-22 |  116.70 |  392 |  80 |  160.92
2017-06-01 |  136.88 |  390 |  70 |  -21.27
2015-03-31 |  163.12 |  387 |  102 |  30.44
2017-05-31 |  142.81 |  387 |  64 |  40.17
2016-01-29 |  150.80 |  386 |  88 |  997.39
2017-12-01 |  123.99 |  378 |  77 |  212.67
2016-07-15 |  138.24 |  371 |  78 |  1,264.70
2016-06-01 |  140.73 |  369 |  92 |  65.04
2015-01-30 |  139.19 |  360 |  55 |  406.47
2016-12-30 |  126.59 |  346 |  97 |  1,104.92
2017-06-15 |  130.75 |  336 |  87 |  57.54
2017-11-15 |  108.92 |  335 |  95 |  56.13
2015-05-15 |  128.94 |  317 |  97 |  24.32
2017-02-21 |  120.83 |  311 |  86 |  206.77
2015-04-02 |  111.10 |  306 |  91 |  52.29
2015-04-15 |  133.56 |  300 |  85 |  40.38
2015-12-16 |  153.42 |  300 |  85 |  14.40
2015-10-26 |  123.29 |  296 |  76 |  299.37
2017-11-28 |  91.92 |  290 |  83 |  371.20
2017-09-20 |  167.21 |  284 |  51 |  33.35
2017-11-16 |  107.90 |  282 |  62 |  56.09
2017-09-26 |  101.61 |  276 |  72 |  213.15
2017-10-12 |  92.62 |  275 |  84 |  97.37
2016-01-08 |  103.40 |  271 |  65 |  170.64
2015-04-20 |  117.34 |  265 |  73 |  241.88
2015-03-26 |  111.17 |  262 |  87 |  167.41
2015-03-23 |  112.44 |  257 |  102 |  60.53
2015-03-04 |  117.99 |  256 |  95 |  160.16
2015-02-23 |  101.36 |  240 |  84 |  542.75
Average |  125.76 |  323 |  81 |  239.93

As summarized in Table 1, the quantum solver can provide significant end-of-
day liquidity savings—averaged at $239.9M with a median of 128.8M across 30
days and following an exponential distribution. The most significant savings
was $1.26B on 2016-07-15. However, out of the 9,717 batches examined across
these days, 49% had no participant who moved their mNDP, and thus no
rearrangement could provide liquidity savings. 25.3% of batches had an
optimized rearrangement that provided liquidity savings. The remaining 25.7%
had participants that grew their mNDP but were unoptimizable for reasons such
as the participants only had outgoing payments, FIFO was already an optimal
solution, or there was a large payment that could not be offset.

As summarized in Table 2, on a typical day, an optimized batch saved $13.25M
on average, with a median of $0.79M and a maximum of $339.41M. Overall, the
maximum savings achieved by any one batch was $1.29B on July 15, 2016.
Occasionally, the order returned by the CQM hybrid solver resulted in less-
efficient liquidity usage than FIFO. This occurred only seven times out of all
the 9,717 batches and never occurred more than once on any day. The loss in
savings was less than $10k in 5 of those instances and reached as high as
$2.96M (Table 2). These instances can be attributed to stochastic noise in the
quantum annealer and could easily be removed in a practical implementation of
this system by immediately re-running the affected batches (if time allowed)
or defaulting to FIFO.

Table 2: Comparison of the results across the batches of $n=70$ payments on each day in the sample using the order from CQM solver against original FIFO order. Average and median savings are calculated using only the batch of payments that saw improvement (for example, using 89 improved batches on 2017-06-30 in Table 1). The number of worsened batches using CQM settlement compared to FIFO are due to stochastic noise in the quantum solver. In these 30 days, out of the 9,717 batches tested, only 7 came back worse than FIFO (0.072% error rate). These errors could be easily identified and corrected by re-running the batches. Date |  Average Savings ($Mn) |  Median Savings ($Mn) |  Maximum Savings ($Mn) |  Num. Worsened Batches |  Additional Liquidity Required ($Mn)
---|---|---|---|---|---
2017-06-30 |  8.24 |  0.51 |  195.91 |  \- - |  \- -
2017-12-22 |  17.76 |  1.02 |  536.12 |  \- - |  \- -
2017-06-01 |  13.37 |  0.57 |  446.99 |  1 |  -0.0060
2015-03-31 |  8.59 |  1.01 |  129.17 |  \- - |  \- -
2017-05-31 |  11.99 |  0.45 |  205.99 |  \- - |  \- -
2016-01-29 |  19.93 |  0.64 |  568.12 |  \- - |  \- -
2017-12-01 |  5.85 |  0.38 |  150.78 |  \- - |  \- -
2016-07-15 |  20.52 |  0.42 |  1,289.03 |  \- - |  \- -
2016-06-01 |  6.20 |  0.31 |  92.39 |  \- - |  \- -
2015-01-30 |  9.48 |  0.48 |  387.78 |  \- - |  \- -
2016-12-30 |  17.00 |  0.43 |  929.20 |  \- - |  \- -
2017-06-15 |  9.31 |  1.71 |  141.24 |  \- - |  \- -
2017-11-15 |  5.60 |  0.36 |  331.91 |  \- - |  \- -
2015-05-15 |  6.52 |  0.68 |  163.31 |  \- - |  \- -
2017-02-21 |  14.79 |  0.46 |  266.57 |  \- - |  \- -
2015-04-02 |  22.44 |  1.40 |  597.34 |  \- - |  \- -
2015-04-15 |  10.71 |  1.45 |  295.93 |  1 |  -0.0012
2015-12-16 |  15.12 |  0.44 |  318.21 |  \- - |  \- -
2015-10-26 |  11.44 |  0.87 |  234.21 |  1 |  -0.0084
2017-11-28 |  13.61 |  1.43 |  200.57 |  \- - |  \- -
2017-09-20 |  21.71 |  2.18 |  331.22 |  \- - |  \- -
2017-11-16 |  14.78 |  0.91 |  219.75 |  \- - |  \- -
2017-09-26 |  13.14 |  0.50 |  361.89 |  1 |  -2.96
2017-10-12 |  8.22 |  0.32 |  171.14 |  1 |  -0.0076
2016-01-08 |  12.68 |  0.62 |  457.88 |  1 |  -0.0013
2015-04-20 |  23.63 |  0.25 |  334.58 |  \- - |  \- -
2015-03-26 |  13.48 |  1.02 |  190.14 |  \- - |  \- -
2015-03-23 |  15.47 |  1.65 |  241.80 |  \- - |  \- -
2015-03-04 |  8.52 |  0.62 |  139.14 |  \- - |  \- -
2015-02-23 |  17.29 |  0.60 |  253.94 |  1 |  -0.21
Average |  13.25 |  0.79 |  339.41 |  0.16% |  -0.11

To examine how equally the savings were distributed across participants, in
Figure 6, we compare each participant’s liquidity savings (using a batch size
of 70) to their total transaction participation. The percentage of the total
amount saved for each participant is compared to their percentage of total
incoming and outgoing transactions. The Pearson’s correlation coefficient
between percent saved and percent of incoming and outgoing transactions is
r=0.961 and 0.967, respectively. This high correlation implies the savings are
somewhat fairly distributed based on how much liquidity each participant
moves. In other words, if a participant is involved in payments responsible
for half of the total settled value in a day, they will likely see about half
the savings achieved for that day.

Figure 6: Amount of liquidity saved per bank as a percentage of total savings
incurred (red) compared with the percentage of the total value of incoming
(gold) and outgoing (green) payments using a fixed batch size.

In most cases, including ours, the payment values follow a Pareto
distribution, maintained in the total settled values for each batch (Figure 7
\- left). In our sample, the average value settled for a batch size of 70 is
$0.4 billion, but some settle for $3.5 billion. The non-normal distribution of
payments makes it difficult to find bilateral or multilateral netting
opportunities. To further reduce the optimization capability, most batches
have more payees than payers. Figure 7 (right) shows the number of payees and
payers per batch using transparency to represent their frequency. On average,
there are eleven payers and eight payees. The presence of a lone payee in
batches acts as a liquidity sink removing it from the liquidity pool. We infer
that a strongly interconnected network of equal value payments in each batch
would provide an opportunity for larger liquidity savings.

Figure 7: (Left) Total value (in CAD) settled per batch in our sample. (Right)
The scatter plot of number of unique payers and payees involved in each batch
in our sample (darker blue dots show more repetitions and lighter blue dots
shows fewer repetitions.

###  5.4 Performance of Quantum Algorithm on a Chosen Day

Next, as an example, we show the evolution of the aggregate mNDP for a typical
day (Sept 20, 2017) in our sample in Figure 8 for a batch size of $n=70$. At
the start of the day, as participants send their payments, we see a sharp rise
in aggregated mNDP. As the day progresses, however, we see a slower increase
in aggregated mNDP. For many batches during the day, no participant had moved
their mNDP, and thus no rearrangement could provide liquidity savings. These
are seen as flat sections in the aggregate mNDP curve in Figure 8.

Sometimes, a large payment (or group of payments) in a single batch can cause
the aggregate mNDP to converge for the CQM solver and FIFO, erasing savings
that the CQM had achieved earlier in the day. Examples of this can be seen at
14:20 and more extremely at 16:05 in Figure 8. This cannot be avoided without
increasing the queue size or using a more complicated queue-building algorithm
that allows payments to move between queues to create larger graph cycles and
offset the required liquidity needed to settle the large payment(s). Clearly,
for a chosen batch size, there is significant heterogeneity across the
batches, and there are limitations on the improvements which can be achieved
with this simple pre-processing setup.

Figure 8: The aggregate max net debit position (mNDP) versus time of the day
for Sept 20, 2017. The blue line shows the original payment order (FIFO), the
orange and green lines respectively shows the orders returned from the quantum
annealer with the day divided into batch sizes of 70 and 140 payments.

To help visualize the queue reordering process in a typical batch of 70
payments (on Sept 20, 2017), each participant’s debit and credit balance
evolution is shown in Figure 9. The top chart uses the FIFO order, and the
bottom chart is for the order proposed by quantum optimization. Each line
represents each bank’s balance as the transactions are completed sequentially
in the prescribed order. Since we optimize the aggregate mNDP, we want to
limit the depth to which the participants drop below zero. In this example,
the mNDP for bank 14 originally was $17 million. The optimization rearranged
the transactions to prevent bank 14 (blue line) from taking any debit position
while not penalizing other participants. Doing so, and reordering other banks’
payments, saved the system 18 million in aggregate mNDP over these 70
transactions. However, it should be noted that not all transactions in this
batch can be optimized. Two examples are banks 13 (teal) and 16 (orange),
which only have outgoing payments; therefore, no scope exists for optimizing
their mNDP through queue reordering any of their transactions in this batch.

Figure 9: The change in individual banks credit or debit position for a given
queue of 70 transactions immediately after each payment settlement for Sept
20, 2017. The top is for FIFO order ($73 million Aggregate mNDP), and the
bottom is for the order proposed by the hybrid quantum CQM solver ($55 million
Aggregate mNDP).

To understand the network effects in terms of the number of payees and payers
per batch, in Figure 10 we show the participants that cause the aggregate mNDP
to increase (if at all) for each batch as a stacked bar plot with FIFO
settlement on Sept 20, 2017. The increase is usually the result of only one
participant (a notable exception is batch 98, for example). Therefore,
minimizing these participants’ mNDP would provide the most significant
liquidity savings. However, the queue rearrangement can do little if they have
a large negative change to their net position in a single batch; i.e., they do
not have incoming payments to offset their one or more large outgoing
payments. For this reason, FIFO and CQM with 70 and 140 queue size track each
other with varying offset and can exhibit the convergence that erases gains
discussed above (see Figure 8). Also evident are groupings throughout the day
where one participant moves the aggregate mNDP in succession. These are cycles
where participants have repeated outgoing transactions and a low balance.

Figure 10: The change in Aggregate max net debit position against the
accumulated batches of payments for Sept 20, 2017. The blue line represents
the Aggregated mNDP (on the left axis) and stacked bars represents magnitude
of the change in mNDP caused by the individual participants (on the right
axis).

Finally, the algorithm was run by doubling the batch size to $n=140$ for two
test days: Sept 20 and Oct 12, 2017. Similar to Table 1, we ran a horse race
between the FIFO and CQM, where they both started with the same initial
conditions and settled every batch throughout the day using FIFO and CQM
order, and at the end of the day, we report the mNDP. As detailed in Table 3,
this produced substantial improvements in savings, as expected. Most notably,
the number of batches that saw higher savings over FIFO increased by close to
15% for both days with respect to the number of batches that could potentially
be optimized. This was certainly due to the increase in both the amount of
liquidity available for recycling and the increased connectivity between
payers and payees in a single batch’s network graph.

Table 3: Comparison of the results for two days using the order proposed by
the CQM hybrid solver and the original order (FIFO) with batch size $n=140$.
The improvements over the same method using the the batch size of $n=70$ are
shown in bold.

Date |  Num. of Batches |  Optimizable Batches |  Num. Improved Batches |  Average Savings ($1M) |  Median Savings ($1M) |  Maximum Savings ($1M) |  End-of-Day Savings ($1M)
---|---|---|---|---|---|---|---
2017-09-20 |  142 |  87 |  50 |  31.46 |  5.93 |  508.41 |  359.56
Improvement Over n=70: |  +13.9% |  +9.75 |  +3.75 |  +177.19 |  +326.21
2017-10-12 |  138 |  110 |  73 |  16.84 |  1.16 |  171.15 |  191.84
Improvement Over n=70: |  +14.5% |  +8.62 |  +0.84 |  +0.01 |  +94.47

On Sept 20, 2017, we saw an impressive improvement in its average, median, and
maximum savings of approximately 45%, 172%, and 53%, respectively. On Oct 12,
2017, we saw an even greater improvement in its average and median savings of
approximately 105% and 263%, respectively, while its maximum savings saw a
negligible improvement of less than 1%. However, maximum savings as a
performance metric is the most dependent on how (and if) the day’s largest
payments are grouped into the queues, and it is expected to be the least
consistent of the three measures.

Arguably the most encouraging improvement was in the end-of-day savings, where
these two days saw an increase of approximately $326M (978%) and $192M (97%),
respectively. These dramatic changes can be attributed to both the overall
superior performance of the larger queue size, as well as a reduction in the
loss of mid-day savings. This is particularly evident for 20 Sept, 2017, also
depicted in Figure 8.

##  6 Conclusions and Future Work

We have shown that substantial liquidity savings can be achieved with the
addition of a pre-processing layer where payments are placed in a queue and
reordered using an algorithm making use of hybrid quantum-classical resources
before being submitted to the RTGS for settlement. Using conventional
computing hardware, optimizing the queue requires prohibitively long
computation times to achieve results for all but trivial queue sizes. However,
quantum computing via hybrid annealing reduces computation times such that
queues that are large enough to provide significant savings-potential can be
optimized. On a sample of 30 representative days worth of transactions from
Canada’s HVPS, we demonstrated this for queue sizes of 70 and 140, with
average aggregate daily liquidity savings of C$239.93M and C$275.70M,
respectively. Furthermore, these savings tended to be distributed
proportionately to the participants’ payments values.

Although much larger savings were achieved on an intraday basis (one batch
reached nearly C$1.3B), the limit in the potential savings is given by the
timing and structure of payment flows between participants, especially when
large payments appear in queues with few or no other payments involving that
payer. In such cases, minimal improvements are possible. This suggests avenues
for future improvements. For example, combining a by-pass component (used
today in traditional LSMs) that sets aside larger payments allowing other
payments to be optimized. Likewise, the progress of quantum technology will
permit optimizing larger queue sizes in the future which could help overcome
such limitations. Lastly, we may extend our algorithm to implement a flexible
batch size, based on a combination of queue time, the value to process, and
the balance of payers and payees, to improve the chances of offsetting
payments occurring in the same batch.

In the end, we found the answer that we sought: quantum technologies readily
available today _can_ be used to optimize the order of payments. While the
trade-off of liquidity and delay remains, the small delay required by our
quantum pre-processor can provide substantial liquidity savings in payment
systems. Ultimately, we demonstrated that a quantum method could reduce
liquidity usage without significantly increasing settlement delay, and this
paves the way for future research as these technologies mature.

## References

  * [1] Morten Bech and Rodney Garratt.  The intraday liquidity management game.  Journal of Economic Theory, 109(2):198–219, 2003.  doi: 10.1016/S0022-0531(03)00016-4.
  * [2] Antoine Martin and James McAndrews.  Liquidity-saving mechanisms.  Journal of Monetary Economics, 55(3):554–567, 2008.
  * [3] Francisco Rivadeneyra and Nellie Zhang.  Liquidity usage and payment delay estimates of the new canadian high value payments system.  Technical Report 9, Bank of Canada Staff Discussion Paper, 2020.
  * [4] Biliana Alexandrova-Kabadjova, Anton Badev, Saulo Benchimol Bastoso, Evangelos Benos, Freddy Cepeda-López, James Chapman, Martin Diehl, Ioana Duca-Radu, Rodney Garratt, Ronald Heijmans, Anneke Kosse, Antoine Martin, Thomas Nellen, Thomas Nilsson, Jan Paulick, Andrei Pustelnikov, Francisco Rivadeneyra, Mario Rubem do Coutto Bastoso, and Sara Testi.  Intraday liquidity around the world.  Mimeo, 2022.
  * [5] Michael Jünger, Gerhard Reinelt, and Giovanni Rinaldi.  The traveling salesman problem.  Handbooks in operations research and management science, 7:225–330, 1995.
  * [6] National Academies of Sciences, Engineering and Medicine.  Quantum Computing: Progress and Prospects.  The National Academies Press, Washington, DC, 2019.
  * [7] Daniel J. Egger, Claudio Gambella, Jakub Marecek, Scott McFaddin, Martin Mevissen, Rudy Raymond, Andrea Simonetto, Stefan Woerner, and Elena Yndurain.  Quantum computing for finance: State-of-the-art and future prospects.  IEEE Transactions on Quantum Engineering, 1:1–24, 2020.
  * [8] Isaiah Hull, Or Sattath, Eleni Diamanti, and Göran Wendin.  Quantum technology for economists.  Available at SSRN 3745608, 2020.  http://dx.doi.org/10.2139/ssrn.3745608.
  * [9] Bank of England.  Liquidity saving mechanism user guide.  Technical report, Bank of England, 2021.  https://www.bankofengland.co.uk/-/media/boe/files/payments/liquidity-saving-mechanism-user-guide.pdf.
  * [10] Andrew Lucas.  Ising formulations of many NP problems.  Frontiers in Physics, 2:5, 02 2014.  doi: 10.3389/fphy.2014.00005.
  * [11] Fred Glover, Gary Kochenberger, and Yu Du.  Quantum bridge analytics I: a tutorial on formulating and using QUBO models.  4OR, 17, 12 2019.  doi: 10.1007/s10288-019-00424-y.
  * [12] D-Wave Systems.  Hybrid solver for constrained quadratic models.  Technical report, D-Wave Systems White Paper, 2021.  https://www.dwavesys.com/media/rldh2ghw/14-1055a-a_hybrid_solver_for_constrained_quadratic_models.pdf.
  * [13] Martin Diehl and Uwe Schollmeyer.  Liquidity-saving mechanisms: quantifying the benefits in TARGET2.  Diagnostics for the financial markets–computational studies of payment system, 2009.  https://helda.helsinki.fi/bof/bitstream/handle/123456789/9381/170846.pdf?sequence=1#page=412.
  * [14] Nick Davey and Daniel Gray.  How has the liquidity saving mechanism reduced banks’ intraday liquidity costs in CHAPS?  Bank of England Quarterly Bulletin, page Q2, 2014.  https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2464096.
  * [15] Enghin Atalay, Antoine Martin, and James McAndrews.  Quantifying the benefits of a liquidity-saving mechanism.  FRB of New York Staff Report Number 447, 2010.  http://dx.doi.org/10.2139/ssrn.1607104.
  * [16] Ben Norman.  Liquidity saving in real-time gross settlement systems: An overview.  Journal of Payments Strategy & Systems, 4(3):261–276, 2010.
  * [17] Antoine Martin and James McAndrews.  An economic analysis of liquidity-saving mechanisms.  Economic Policy Review, pages 25–39, 01 2008.  doi: 10.2139/ssrn.1141361.
  * [18] Marius Jurgilas and Antoine Martin.  Liquidity-saving mechanisms in collateral-based RTGS payment systems.  Annals of Finance, 9:29–60, 2013.
  * [19] Marco Galbiati and Kimmo Soramaki.  Liquidity-saving mechanisms and bank behaviour.  Technical report, Bank of England Working Paper, 2010.  https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.173.1040&rep=rep1&type=pdf.
  * [20] Francisco Rivadeneyra and Nellie Zhang.  Payment coordination and liquidity efficiency in the new canadian wholesale payments system.  Bank of Canada Staff Discussion Paper 2022-3, Bank of Canada, 2022.
  * [21] Rodney J Garratt.  An application of Shapley value cost allocation to liquidity savings mechanisms.  Journal of Money, Credit and Banking, 2019.  doi: 10.1111/jmcb.12853.
  * [22] Andrew Steane.  Quantum computing.  Reports on Progress in Physics, 61(2):117, 1998.
  * [23] Jack D Hidary.  Quantum Computing: An Applied Approach.  Springer, 2019.  doi: 10.1007/978-3-030-23922-0.
  * [24] Lee Braine, Daniel Egger, Jennifer Glick, and Stefan Woerner.  Quantum algorithms for mixed binary optimization applied to transaction settlement.  IEEE Transactions on Quantum Engineering, PP:1–1, 03 2021.  doi: 10.1109/TQE.2021.3063635.
  * [25] Pablo S Castro, Ajit Desai, Han Du, Rodney Garratt, and Francisco Rivadeneyra.  Estimating policy functions in payment systems using reinforcement learning.  Available at SSRN 3743017, 2020.  https://www.bankofcanada.ca/wp-content/uploads/2021/02/swp2021-7.pdf.
  * [26] Alpha C. Chiang.  Elements of Dynamic Optimization.  Waveland Press,, 1992.
  * [27] David Cass and Karl Shell.  The Hamiltonian approach to dynamic economics.  Academic Press, 2014.
  * [28] Ksenia Bestuzheva, Mathieu Besançon, Wei-Kun Chen, Antonia Chmiela, Tim Donkiewicz, Jasper van Doornmalen, Leon Eifler, Oliver Gaul, Gerald Gamrath, Ambros Gleixner, Leona Gottwald, Christoph Graczyk, Katrin Halbig, Alexander Hoen, Christopher Hojny, Rolf van der Hulst, Thorsten Koch, Marco Lübbecke, Stephen J. Maher, Frederic Matter, Erik Mühmer, Benjamin Müller, Marc E. Pfetsch, Daniel Rehfeldt, Steffan Schlein, Franziska Schlösser, Felipe Serrano, Yuji Shinano, Boro Sofranac, Mark Turner, Stefan Vigerske, Fabian Wegscheider, Philipp Wellner, Dieter Weninger, and Jakob Witzig.  The SCIP Optimization Suite 8.0.  ZIB-Report 21-41, Zuse Institute Berlin, December 2021.
  * [29] Stefan Vigerske and Ambros Gleixner.  SCIP: global optimization of mixed-integer nonlinear programs in a branch-and-cut framework.  Optimization Methods and Software, 33:1–31, 06 2017.
  * [30] Rimmi Anand, Divya Aggarwal, and Vijay Kumar.  A comparative analysis of optimization solvers.  Journal of Statistics and Management Systems, 20(4):623–635, 2017\.
  * [31] Neville Arjani and Darcey McVanel.  A primer on canada’s large value transfer system.  Technical report, Bank of Canada Staff Working Paper, 2006.  https://www.bankofcanada.ca/wp-content/uploads/2010/05/lvts_neville.pdf.

##  Appendix A Formulation of the Objective Function as Hamiltonian

Queue optimization is done analytically through the minimization of an
objective function. A number of payments, $n$, is collected in a queue. The
payments are ordered with index $i=1\ldots n$, while the rearranged queue has
the alternative index $t=1\ldots n$. The optimization of the Hamiltonian
returns the mapping $i\rightarrow t$ that minimizes the amount of liquidity
needed in the system (the aggregate max net debit position). We define a
participant $\alpha$’s net position in the final queue at position $t$ as
$N_{\alpha}(t)$. Each participant’s maximum net debit position incurred from
previous batches is included as $mNDP_{\alpha}$ (zero for the initial batch).
Thus the collateral needed by each participant to ensure they do not go into a
negative liquidity position is

${mNDP}_{\alpha}=|\min\\{N_{\alpha}(\text{previous\ queues}),0\\}|.$ (7)

The participant can access liquidity from either $N(t)$ (if their position
happens to be a net credit) or this collateral. More precisely, their total
liquidity available to settle payment $t$ is effectively

$\ell_{\alpha}(t)=N_{\alpha}(t)+{mNDP}_{\alpha}.$ (8)

A participant’s instantaneous net position can be calculated using

$N_{\alpha}(t)=N_{\alpha}(0)+\sum_{i=1}^{n}\sum_{\tau=1}^{t}f(\alpha,i)\,x_{i,\tau},$
(9)

where $N_{\alpha}(t)\geq-{mNDP}_{\alpha}$, $N_{\alpha}(0)$ is participant
$\alpha$’s net position immediately before the settlement of this batch
(starting balance), $x_{i,\tau}$ is a binary decision variable that indicates
whether payment $i$ is settled ($x_{i,\tau}=1$) or not ($x_{i,\tau}=0$) at
position $\tau$ in the final queue,

$f(\alpha,i)=\begin{cases}\ v,&\text{if }\alpha\text{ is the payee of payment
}i\\\ -v,&\text{if }\alpha\text{ is the payer of payment }i\\\ \
0,&\text{otherwise}\end{cases},$ (10)

and $v$ is the value of payment $i$.

Ideally, we seek a mapping $i\rightarrow t$ such that $\ell_{\alpha}(t)\geq 0\
\forall\alpha,t$. However, such a mapping may not exist for a given batch. For
the general case, we define the positive real number $b_{\alpha}$ that is a
hypothetical amount of liquidity one would need to add to participant
$\alpha$’s account in order to keep their balance non-negative for all
payments $t$. We now have the inequality

$b_{\alpha}+N_{\alpha}(0)+{mNDP}_{\alpha}+\sum_{i=1}^{n}\sum_{\tau=1}^{t}f(\alpha,i)\,x_{i,\tau}\geq
0$ (11)

that can always be satisfied for the trivial case
$b_{\alpha}\rightarrow\infty$. The optimal mapping occurs when $b_{\alpha}$ is
collectively minimized across $\alpha$. Quite simply,

$\min_{b_{\alpha},x_{i,t}}\sum_{\alpha}b_{\alpha}.$ (12)

Note that the optimal $b_{\alpha}$ represents the required change that will
have to be made to $mNDP_{\alpha}$ before the next queue is optimized.
Equivalently, this is the additional collateral that the central bank will
demand from the participant $\alpha$ to provide the liquidity before it allows
settlement of the current queue to commence.

This problem is of the mixed binary optimization (MBO) type where the
optimization is over both continuous ($b_{\alpha}$) and binary ($x_{i,t}$)
variables. We can easily transform this into a quadratic unconstrained binary
optimization (QUBO) problem for processing on a quantum device [11, 10].
First, the inequality constraint given in Equation 11 must be written as an
equality constraint with the addition of a slack variable $s_{\alpha}(t)$:

$b_{\alpha}+N_{\alpha}(0)+\text{mNDP}_{\alpha}+\sum_{\tau}^{t}\sum_{i=1}^{\tau}f(\alpha,i)\,x_{i,\tau}-s_{\alpha}(t)=0$
(13)

where $s_{\alpha}(t)\geq 0$ is an additional continuous variable over which
optimization must occur. Second, this constraint must be enforced in the QUBO
as a quadratic penalty term with a Lagrange multiplier. Third, all continuous
variables ($b_{\alpha}$, $s_{\alpha}(t)$) must be discretized to the nearest
cent (or larger base value with the loss of some resolution) and then
represented with the so-called “log trick.” For example, $b_{\alpha}$ is
replaced with the sum of a geometric series of binary variables
$\beta_{\alpha}$,

$b_{\alpha}=\sum_{j}2^{j}\beta_{\alpha}^{(j)},$ (14)

where $\beta_{\alpha}^{(0)}$ represents the least-significant digit in the
binary representation of $b_{\alpha}$, $\beta_{\alpha}^{(1)}$ the second-
least-significant digit, etc.

Finally, two more constraints must be added to ensure that every payment is
settled exactly once; to guarantee the mapping $i\rightarrow t$ is one-to-one.
To ensure that every $i$ corresponds to exactly one $t$, we have

$1-\sum_{t}x_{i,t}=0\ \forall i,$ (15)

and similarly that every $t$ corresponds to exactly one $i$, we have

$1-\sum_{i}x_{i,t}=0\ \forall t.$ (16)

In the end, we are left with the equation

$\displaystyle H=$
$\displaystyle\sum_{\alpha}b_{\alpha}+\sum_{\alpha,t}\lambda_{1}\left(b_{\alpha}+N_{\alpha}(0)+\text{mNDP}_{\alpha}+\sum_{i}\sum_{\tau\leq
t}f(\alpha,i)\,x_{i,\tau}-s_{\alpha}(t)\right)^{2}$
$\displaystyle+\sum_{i}\lambda_{2}\left(1-\sum_{t}x_{i,t}\right)^{2}+\sum_{t}\lambda_{2}\left(1-\sum_{i}x_{i,t}\right)^{2},$
(17)

where the log trick has been omitted for brevity. Minimizing this Hamiltonian
over all variables

$\min_{b_{\alpha},\,x_{i,t},\,s_{\alpha}(t)}H$ (18)

yields the optimal reordering of the payment queue. Note that the same
Lagrange multiplier, $\lambda_{2}$, can be assumed for both one-to-one
constraints without loss of generality due to their symmetry.