11institutetext: Qing Zhou 22institutetext: School of Cyber Science and
Engineering, Huazhong University of Science and Technology, Wuhan 430074,
China
33institutetext: Songfeng Lu 44institutetext: School of Cyber Science and
Engineering, Huazhong University of Science and Technology, Wuhan 430074,
China
Corresponding author, 44email<EMAIL_ADDRESS>

# Hash function based on controlled alternate quantum walks with memory

Qing Zhou  Songfeng Lu

(Received: date / Accepted: date)

###### Abstract

We propose a new hash function QHFM based on controlled alternate quantum
walks with memory on cycles, where the $j$th message bit decides whether to
run quantum walk with one-step memory or to run quantum walk with two-step
memory at the $j$th time step, and the hash value is calculated from the
resulting probability distribution of the walker. Numerical simulation shows
that the proposed hash function has near-ideal statistical performance and is
at least on a par with the state-of-the-art hash functions based on quantum
walks in terms of sensitivity of hash value to message, diffusion and
confusion properties, uniform distribution property, and collision resistance
property; and theoretical analysis indicates that the time and space
complexity of the new scheme are not greater than those of its peers. The good
performance of QHFM suggests that quantum walks that differ not only in coin
operators but also in memory lengths can be combined to build good hash
functions, which, in turn, enriches the construction of controlled alternate
quantum walks.

###### Keywords:

Hash function Quantum walks with memory Controlled alternate quantum walks
Statistical properties time and space complexity

††journal: Journal name to insert

##  1 Introduction

As one of the principal tools of information security, cryptographic hash
functions not only act as essential components of identification, message
authentication, digital signatures, and random number generation, but also
play an important part in privacy amplification process of quantum key
distribution Bennett1995PrivAmplif . Classical hash functions based on hard
computational problems are, however, subject to an inherent security
limitation: the existence of one-way functions is still an open conjecture
that cannot be proved (a proof, with no assumptions, of existence would
establish P$\neq$NP Menezes1996Book ). As a result, they only satisfy
computational security and are challenged by cryptanalysis equipped with
quantum algorithms. Such a fact stimulates researchers to develop hash
functions with a higher level of security, such as hash functions based on (or
inspired by) quantum computing Ablayev2016qHash ; Ablayev2016Balanced ;
Vasiliev2016qHash ; Ziatdinov2016qHash ; Ablayev2013qHash ; Ablayev2020uniHash
; Yang2021qwHash ; Yang2019qwHash ; Yang2018qwHash264 ; Yang18qwHash221 ;
Li2018qwHash ; Cao2018qwHash ; Yang2016qHash ; Li2013qwHash , whose preimage
resistance property is ensured by quantum mechanics rather than hardness
assumptions.

There are two kinds of quantum-computing-based hash functions: classical-
quantum hash functions based on quantum one-way functions Ablayev2016qHash ;
Ablayev2016Balanced ; Vasiliev2016qHash ; Ziatdinov2016qHash ;
Ablayev2013qHash ; Ablayev2020uniHash  (hereafter, simply QOWF-based hash
functions) and classical-classical hash functions based on discrete quantum
walks on cycles Yang2021qwHash ; Yang2019qwHash ; Yang2018qwHash264 ;
Yang18qwHash221 ; Li2018qwHash ; Cao2018qwHash ; Yang2016qHash ; Li2013qwHash
(hereafter, simply QW-based hash functions). The former have balanced one-way
resistance property and collision resistance property that are well-defined
and strictly proved in the quantum setting, the latter take advantage of the
chaotic characteristics of quantum walks and belong to dedicated hash
functions, whose capabilities of collision resistance are difficult to prove
and are mainly assessed by means of statistical analysis. On the other hand,
of QOWF-based hash functions, the output length and the number of hashing
parameters are both positively correlated with the input size, while QW-based
hash functions map messages of arbitrary finite length to digests of fixed
length. In addition, the output length of QW-based hash functions can be
easily extended (to withstand brute-force attacks) by increasing the number of
nodes of the cycle or the number of hash bits “contributed” by each node, and
the hash result can be calculated classically. Thus, QW-based hash functions
are of greater practical utility before large-scale quantum computers are
built, for they can currently be used to improve the security of hash-
function-based schemes.

The essence of the design of QW-based hash functions is combining two or more
different quantum walk procedures governed by evolution operators
$\\{U_{0},U_{1},\dots\\}$ to construct a controlled alternate quantum walk
(CAQW) model, where the choice among $\\{U_{0},U_{1},\dots\\}$ at the $j$th
time step is determined by the $j$th bit of a binary string. Theoretically, a
valid CAQW model could be constructed if the walker can ”switch” freely among
$\\{U_{0},U_{1},\dots\\}$, and evolution operators that only differ in coin
transform naturally satisfy this requirement. Therefore, various quantum
walks, such as one-dimensional broken-line quantum walks Yang2021qwHash , one-
dimensional one-particle quantum walks Yang2019qwHash ; Yang2018qwHash264 ;
Yang18qwHash221 , two-dimensional one-particle quantum walks Li2018qwHash ,
quantum walks on Johnson graphs Cao2018qwHash , one-dimensional two-particle
(interacting) quantum walks Yang2016qHash ; Li2013qwHash , and one-dimensional
quantum walks with memory (QWM) Li2020QWM ; Dai2020QWM ; Zhou2019QW2M ;
Li2016QWM ; Gettrick2014QWM ; Gettrick2010QW1M  can all be used to construct
valid (but may not good) hash functions as long as they are modified to
utilize coin operators controlled by input messages. Among these walks, the
evolution of QWM is governed by three (rather than two) stages: flipping a
coin, determining the next direction according to the coin state and the
previous direction(s), and moving a step according to the new direction. Here,
an extra operator—the direction-determine transform—could be taken into
account when designing hash functions based on QWM: the alternately performed
evolution operators can differ in coin transform or direction-determine
transform, or both.

To examine the feasibility and utility of this idea, we combine two quantum
walks with different memory lengths, i.e., QW1M Gettrick2010QW1M  and QW2M
Zhou2019QW2M , to achieve a valid CAQW process, for different quantum walks
with unequal memory lengths typically have different direction-determine
transforms, and they can also use different coin operators. Based on this
walking process, we construct a new hash function (named QHFM) and then assess
its performance. Simulation results show that the statistical properties of
the proposed hash function are as good as those of the existing QW-based hash
functions, and theoretical analysis indicates that the space and time
complexity of QHFM are not greater than those of theirs peers.

##  2 One-dimensional Controlled Quantum Walk with One- and Two-step Memory

A one-dimensional controlled quantum walk with one- and two-step memory (CQWM)
takes place in the Hilbert space
$\mathcal{H}_{p}\otimes\mathcal{H}_{dr_{2}}\otimes\mathcal{H}_{dr_{1}}\otimes\mathcal{H}_{c}$
spanned by vectors $\Ket{x,dr_{2},dr_{1},c}$, where $c$ (with $c\in\\{0,1\\}$,
or $c\in\mathbb{Z}_{2}$) is the coin state, $dr_{1}$ (with
$dr_{1}\in\mathbb{Z}_{2}$) is the direction of the most recent step (0 stands
for left and 1 stands for right), $dr_{2}$ (with $dr_{2}\in\mathbb{Z}_{2}$) is
the direction of the penultimate step, and $x$ is the current position. If the
walker moves on a line, then $x\in\mathbb{Z}$ (all integers); and if the
walker moves on a cycle with $n$ nodes, then $x\in\\{0,1,2,\ldots,n-1\\}$ (or
$x\in\mathbb{Z}_{n}$).

Formally, the evolution of CQWM controlled by a $t$-bit string $msg=$
$(m_{1},m_{2},\dots,m_{t})\in\\{0,1\\}^{t}$ is the product of $t$ unitary
transforms

$U_{msg}=U^{(m_{t})}U^{(m_{t-1})}\cdots U^{(m_{2})}U^{(m_{1})},$ (1)

Here $U^{(m_{j})}$ ($1\leq j\leq t$) is the one-step transform controlled by
the $j$th bit of $msg$, and it is defined as

$U^{(m_{j})}=S\cdot\left(I_{n}\otimes
D^{(m_{j})}\right)\cdot\left(I_{4n}\otimes C^{(m_{j})}\right),$ (2)

where $C^{(m_{j})}$ is a $2\times 2$ coin operator controlled by $m_{j}$,
$I_{k}$ ($k=4n$ or $n$) is a $k\times k$ identity operator, $D^{(m_{j})}$ is
an $8\times 8$ direction-determine operator controlled by $m_{j}$, and $S$ is
the conditional shift operator controlled by the next direction. If $m_{j}=0$,
then $C^{(m_{j})}$ is parameterized by an angle $\theta_{0}$, i.e.,

$C^{(0)}=\begin{pmatrix}\text{cos}(\theta_{0})&\text{sin}(\theta_{0})\\\
\text{sin}(\theta_{0})&-\text{cos}(\theta_{0})\end{pmatrix},$ (3)

and $D^{(m_{j})}$ (becomes $D^{(0)}$) describes the direction-determine
process of QW1M; if $m_{j}=1$, then $C^{(m_{j})}$ is parameterized by another
angle $\theta_{1}$, i.e.,

$C^{(1)}=\begin{pmatrix}\text{cos}(\theta_{1})&\text{sin}(\theta_{1})\\\
\text{sin}(\theta_{1})&-\text{cos}(\theta_{1})\end{pmatrix},$ (4)

and $D^{(1)}$ describes the direction-determine process of QW2M.

The direction-determine transforms of QW1M Gettrick2010QW1M  and QW2M
Zhou2019QW2M  can be respectively written as
$DT_{0}:\Ket{dr_{1},c}\to\Ket{dr_{1}\oplus\bar{c},c}$ and
$DT_{1}:\Ket{dr_{2},dr_{1},c}\to\Ket{dr_{1},dr_{2}\oplus\bar{c},c}$, where
$\bar{c}\equiv 1\oplus c$, $dr_{1}\oplus\bar{c}$ in the first expression
specifies the next direction of the walker performing QW1M, and
$dr_{2}\oplus\bar{c}$ in the second expression specifies the next direction of
the walker performing QW2M. To enable QW1M and QW2M to be performed
alternately, one may add a redundant state $\Ket{dr_{2}}$ into QW1M and let
$D^{(0)}:\Ket{dr_{2},dr_{1},c}\to\Ket{dr_{2},dr_{1}\oplus\bar{c},c}$
determines the next direction when the controlling bit equals 0; otherwise,
the next direction is determined by $D^{(1)}=DT_{1}$.

Following Ref. Zhou2019QW2M , any 4-term basis state $\Ket{x,dr_{2},dr_{1},c}$
in
$\mathcal{H}_{p}\otimes\mathcal{H}_{dr_{2}}\otimes\mathcal{H}_{dr_{1}}\otimes\mathcal{H}_{c}$
can be rewritten as a 2-term basis state
$\Ket{x,j}=\Ket{x,2^{2}dr_{1}+2^{1}dr_{2}+2^{0}c}$ in
$\mathcal{H}_{p}\otimes\mathcal{H}^{8}$, where $\mathcal{H}^{8}$ is the
8-dimensional Hilbert space. conversely, from any 2-term basis state
$\Ket{x,j}$ ($j\in\mathbb{Z}_{8}$), one can deduce the coin value and the most
recent two directions as follows:

$\begin{cases}c=j\bmod 2,\\\ dr_{2}=(j\bmod 4-j\bmod 2)/2,\\\ dr_{1}=(j-j\bmod
4)/4.\end{cases}$ (5)

According to this correspondence, $D^{(0)}$ can be reformulated to

$\Ket{2^{2}dr_{1}+2^{1}dr_{2}+2^{0}c}\to\Ket{2^{2}(dr_{1}\oplus\bar{c})+2^{1}dr_{2}+2^{0}c},$

or, under the 2-term states

$D^{(0)}:\Ket{j}\to\Ket{4\left[(j-j\bmod 4)/4\oplus(j\bmod 2)\oplus
1\right]+j\bmod 4}.$ (6)

Analogously, $D^{(1)}$ can be expressed as

$D^{(1)}:\Ket{j}\to\Ket{4\left[(j\bmod 4-j\bmod 2)/2\oplus(j\bmod 2)\oplus
1\right]+(j-j\bmod 4)/2+j\bmod 2}.$ (7)

With formulas (6) and (7), one can verify that $D^{(0)}$ and $D^{(1)}$ are
both unitary.

Once the next direction, the new $dr_{1}$, is determined, the walker then
moves according to the shift operator controlled by $dr_{1}$. If the walk
takes place on a line, then the action of $S$ is expressed as
$\Ket{x,dr_{2},dr_{1},c}\to\Ket{x+2dr_{1}-1,dr_{2},dr_{1},c}$; if the walk
takes place on a cycle with $n$ nodes, then $S$ becomes
$\Ket{x,dr_{2},dr_{1},c}\to\Ket{x+2dr_{1}-1\pmod{n},dr_{2},dr_{1},c}$, which
can be reformulated (in 2-term states) to

$S:\Ket{x,j}\to\Ket{x+(j-j\bmod 4)/2-1\pmod{n},j}.$ (8)

In formula (8), the next position is calculated using modular arithmetic under
modulus $n$.

##  3 Hash Function Using Quantum Walks with One- and Two-step Memory on
Cycles

The proposed hash function is constructed by running CQWM on a circle with $n$
nodes under the control of the input message $msg$, where each node
contributes $m$ bits to the hash result $H(msg)$. The process of CQWM-based
hash function is described as follows:

  1. (1)

Select the values of parameters $(n,m,l,\theta_{0},\theta_{1},\alpha)$
satisfying the following constraints: $n$ is odd; $n\times m$ equals the bit
length of the hash value; $10^{l}\gg 2^{m}$; and
$\theta_{0},\theta_{1},\alpha\in(0,\pi/2)$.

  2. (2)

Initialize the walker in the state
$\Ket{\psi_{0}}=\text{cos}\alpha\Ket{0,1,0,0}+\text{sin}\alpha\Ket{0,1,0,1}$
(or, in the 2-term state
$\Ket{\psi_{0}}=\text{cos}\alpha\Ket{0,2}+\text{sin}\alpha\Ket{0,3}$).

  3. (3)

Apply $U_{msg}$ to $\Ket{\psi_{0}}$ and generate the resulting probability
distribution $prob=$
$(p_{0},p_{1},\dots,p_{n-1})$, where $p_{x}$ ($x\in\mathbb{Z}_{n}$) is the
probability that the particle locates at node $x$ when the walk is finished.

  4. (4)

The hash value of $msg$ is a sequence of $n$ blocks $H(msg)=B_{0}\lVert
B_{1}\lVert\dots\lVert B_{n-1}$, where each block $B_{x}$ is the $m$-bit
binary representation of $\lfloor p_{x}\cdot 10^{l}\rfloor\bmod 2^{m}$
($\lfloor\cdot\rfloor$ denotes the floor of a number), and $B_{x}\lVert
B_{x+1}$ denotes the concatenation of $B_{x}$ and $B_{x+1}$.

##  4 Statistical Performance Analysis

QHFM, like other QW-based hash functions Yang2021qwHash ; Yang2019qwHash ;
Yang2018qwHash264 ; Yang18qwHash221 ; Li2018qwHash ; Cao2018qwHash ;
Yang2016qHash ; Li2013qwHash , belongs to dedicated hash functions, whose
performances are mainly evaluated through statistical analysis. To make our
statistical tests reusable and usable by anyone else, we perform these tests
on a collection of items (i.e., input messages) randomly drawn from an open
dataset, named “arXiv Dataset”, of about 1.8 million records and upload the
complete MATLAB code for hash tests to “GitHub”. See the dataset at
https://www.kaggle.com/Cornell-University/arxiv and the test code at
https://github.com/Chloe-Zhouqing/Hash-functions-based-on-quantum-walks.

To make comparisons between the proposed scheme and the existing ones with
(detailed) experimental results Yang2021qwHash ; Yang2019qwHash ;
Yang2018qwHash264 ; Yang18qwHash221 ; Li2018qwHash ; Cao2018qwHash ;
Yang2016qHash  in a fair and informative manner, we consider seven “instances”
QHFM-$L$ ($L\in\\{296,264,221,200,195,136,120\\}$) of QHFM, where QHFM-$L$
produces $L$-bit hash values and will be compared with the existing QW-based
hash functions with $L$\- or close-to-$L$-bit output length (QHFM-136 and
QHFM-120 will be compared with the 128-bit scheme in Ref. Yang2016qHash ).
Different instances of QHFM share the same $l$ values, same $\theta_{0}$
values, same $\theta_{1}$ values, and the same $\alpha$ values, which are
taken to be 8, $\pi/4$, $\pi/3$, and $\pi/4$, respectively. Distinction
between QHFM-$L$ and QHFM-$L^{\prime}$ ($L\neq L^{\prime}$) lies in the values
of $n$ and $m$, which are listed in Table 1.

Table 1: Values of parameters chosen for the seven instances of the proposed hash scheme Hash Instances | $n$ | $m$
---|---|---
QHFM-296 | 37 | 8
QHFM-264 | 33 | 8
QHFM-221 | 17 | 13
QHFM-200 | 25 | 8
QHFM-195 | 15 | 13
QHFM-136 | 17 | 8
QHFM-120 | 15 | 8

###  4.1 Sensitivity of Hash Value to Message

Let $msg_{0}$ be an original message and $msg_{j}$ ($j\in\\{1,2,3\\}$) the
slightly modified result of $msg_{0}$, which are obtained under the following
four conditions:

  * Condition 1: Randomly choose an original message $msg_{0}$;

  * Condition 2: Flip a bit of $msg_{0}$ at a random position and then obtain the modified message $msg_{1}$;

  * Condition 3: Insert a random bit into $msg_{0}$ at a random position and then obtain $msg_{2}$;

  * Condition 4: Delete a bit from $msg_{0}$ at a random position and then obtain $msg_{3}$.

The sensitivity of hash value to message is assessed by comparing the hash
values $H(msg_{j})$ of the modified messages with the hash value $H(msg_{0})$
of the original one. In our sensitivity test, a record is randomly picked out
from the arXiv Dataset, then the article abstract within this record serves as
$msg_{0}$.

Corresponding to the conditions above, four hash values in hexadecimal format
produced by QHFM-195 are obtained as follows:

  * Condition 1: $H(msg_{0})=$“3 5A 2B 76 96 74 1C F7 51 09 2E AB 1F CB 6A C0 33 77 46 61 E5 D1 E4 38 EC”;

  * Condition 2: $H(msg_{1})=$“4 BC EC C7 0E A9 2B 5C 5C 93 34 30 69 E9 3A EC 1B D3 D3 95 7B 0F DF 5A 31”;

  * Condition 3: $H(msg_{2})=$“5 00 5C 40 AB AB 2F 26 9B AB D7 AF B5 23 4F 16 20 5C 63 A0 30 6D 5E 0C 15”;

  * Condition 4: $H(msg_{3})=$“0 5D 14 81 F1 29 CB E7 BE CB 01 F6 53 48 E8 90 D4 CD 35 C3 C7 55 DB 80 E8”.

Notice that the the first hexadecimal digit of the hash value under condition
$j$ only represents the first three (rather than four) bits of $H(msg_{j})$,
since the output length of QHFM-195 is not a multiple of four.

The plots of hash values $H(msg_{0})$, $H(msg_{1})$, $H(msg_{2})$, and
$H(msg_{1})$ in binary format are shown in Fig. 1, which indicates that a tiny
modification to the message could cause a significant change in the hash
value, and the positions of those changed bits are evenly distributed over the
entire interval $[1,195]$ of position numbers. A similar result can be
obtained using any other instance of QHFM; thus, the output digest of the
proposed hash scheme is highly sensitive to its input message.

Figure 1: (Color online) Plots of the hash values produced by QHFM-195 under
the four conditions, where $\text{C}j$ stands for Condition $j$
($j\in\\{1,2,3,4\\}$). In the $j$th subgraph with $j>1$, each asterisk (*)
marks a different bit between $H(msg_{j-1})$ and $H(msg_{0})$.

###  4.2 Diffusion and Confusion Properties

The test data for the diffusion and confusion properties of QHFM-$L$ is
collected by making $N$ random draws (with replacement) from the arXiv
Dataset. On each draw, an original message $msg_{0}$ is selected, then a
slightly modified result $msg_{1}$ of this message is obtained by inverting a
bit of $msg_{0}$ at a random position. Let $B_{i}$ be the Hamming distance
between the hash values of the original and modified messages obtained on the
$i$th draw and $N$ the number of draws, the diffusion and confusion properties
(reflecting the avalanche effect) of the proposed hash instances are assessed
based on the following four indicators:

  * $\bullet$

mean changed bit number
$\overline{B}=\begin{matrix}\sum_{i=1}^{N}B_{i}/N\end{matrix}$;

  * $\bullet$

mean changed probability $P=\left[\overline{B}/(n\times m)\right]\times
100\%$;

  * $\bullet$

standard deviation of the changed bit number
$\Delta
B=\sqrt{\left[1/(N-1)\right]\begin{matrix}\sum_{i=1}^{N}\left(B_{i}-\overline{B}\,\right)^{2}\end{matrix}}$;

  * $\bullet$

standard deviation of the changed probability
$\Delta
P=\sqrt{\left[1/(N-1)\right]\begin{matrix}\sum_{i=1}^{N}\left[B_{i}/(n\times
m)-P\right]^{2}\end{matrix}}\times 100\%$.

The ideal values of $\overline{B}$ and $P$ are $(n\times m)/2$ and $50\%$,
respectively; and smaller standard deviations ($\Delta B$ and $\Delta P$) are
more desirable. For a specific hash function with fixed output length,
$\overline{B}$ and $\Delta B$ are directly proportional to $P$ and $\Delta P$,
respectively; thus, only $P$ and $\Delta P$, or a combination of them, e.g.,
$I_{\text{DC}}=(\Delta P+|P-50\%|)/2\times 100\%$, would suffice to assess the
confusion and diffusion properties of this hash function: the smaller
$I_{\text{DC}}$, the better the avalanche effect achieved. The diffusion and
confusion test on QHFM-$L$ is performed with $N=10000$, and the simulation
results are presented in Table 2. For comparison, the reported results (with
$N\geq 10000)$ of the corresponding variables for the existing QW-based hash
schemes Yang2021qwHash ; Yang2019qwHash ; Yang2018qwHash264 ; Yang18qwHash221
; Li2018qwHash ; Cao2018qwHash ; Yang2016qHash  are also listed in the same
table, where the values for Yang21-296 and Yang18-221 (which have multiple
instances) are the test results for the representative instances of them.
Here, the representative instance of a hash scheme $X$ is the one that has the
$P$ result closest to $50\%$ (over all instances of $X$).

Table 2: Diffusion-and-Confusion-Test Results for The Proposed and Existing QW-Based Hash Functions Hash instances (Schemes) | $\overline{B}$ | $P(\%)$ | $\Delta B$ | $\Delta P(\%)$ | $I_{\text{DC}}(\%)$
---|---|---|---|---|---
QHFM-296 | 147.9101 | 49.9696 | 8.5997 | 2.9053 | 1.4679
QHFM-264 | 131.8667 | 49.9495 | 8.1378 | 3.0825 | 1.5665
QHFM-221 | 110.5313 | 50.0142 | 7.4455 | 3.3690 | 1.6916
QHFM-200 | 100.0205 | 50.0103 | 7.1654 | 3.5827 | 1.7965
QHFM-195 | 97.5591 | 50.0303 | 6.9844 | 3.5817 | 1.8060
QHFM-136 | 68.0530 | 50.0390 | 5.8782 | 4.3222 | 2.1806
QHFM-120 | 60.0914 | 50.0762 | 5.4699 | 4.5582 | 2.3172
Yang21-296 Yang2021qwHash  | 147.8640 | 49.9541 | 8.6141 | 2.9102 | 1.4781
Yang19-264 Yang2019qwHash  | 131.6803 | 49.8789 | 8.8877 | 3.3666 | 1.7439
Yang18-264 Yang2018qwHash264  | 132.1108 | 50.0420 | 8.0405 | 3.0457 | 1.5439
Yang18-221 Yang18qwHash221  | 112.7791 | 51.0313 | 8.2029 | 3.7117 | 2.3715
Li18-200 Li2018qwHash  | 99.9010 | 49.9505 | 7.1133 | 3.5567 | 1.8031
Cao18-195 Cao2018qwHash  | 124.7000 | 63.9600 | 6.4300 | 6.3000 | 10.1300
Yang16-128 Yang2016qHash  | 64.2894 | 50.2261 | 5.6686 | 4.4286 | 2.3274

The values of $I_{\text{DC}}$ suggest that the test results for QHFM-264 is
better than that for Yang19-264 but slightly poorer than that for Yang18-264,
and the results for other instances of the proposed hash scheme are better
than those for their peers (QHFM-296 vs. Yang21-296; QHFM-221 vs. Yang18-221;
QHFM-200 vs. Li18-200; QHFM-195 vs. Cao18-195; QHFM-136 and QHFM-120 vs.
Yang16-128). Thus, the diffusion and confusion properties of the proposed hash
function outperform or are at least on a par with the existing QW-based hash
schemes.

###  4.3 Uniform Distribution Analysis

Similar to the case of diffusion and confusion properties, the uniform
distribution property (reflecting the strict avalanche effect) could also be
assessed based on four indicators:

  * $\bullet$

mean number of draws with flipped hash bit (over $n\times m$ bit positions)
$\overline{T}=\begin{matrix}\sum_{j=1}^{n\times m}T_{j}/(n\times
m)\end{matrix}$;

  * $\bullet$

mean percentage of draws with flipped hash bit $Q=\overline{T}/N\times 100\%$;

  * $\bullet$

standard deviation of the number of draws with flipped hash bit
$\Delta T=\sqrt{\left[1/(n\times m-1)\right]\begin{matrix}\sum_{j=1}^{n\times
m}(T_{j}-\overline{T})^{2}\end{matrix}}$;

  * $\bullet$

standard deviation of the percentage of draws with flipped hash bit
$\Delta Q=\sqrt{\left[1/(n\times m-1)\right]\begin{matrix}\sum_{j=1}^{n\times
m}\left[T_{j}/N-Q\right]^{2}\end{matrix}}\times 100\%$;

where $T_{j}$ ($j\in\\{1,2,3,\dots,n\times m\\}$) is the number of draws on
which a bit-flip occurs in the hash value at the $j$th bit position after a
random message bit is inverted. The theoretical values of $\overline{T}$ and
$Q$ are $(n\times m)/2$ and $50\%$, respectively.

Since $\overline{T}$ and $\Delta T$ are directly proportional to $Q$ and
$\Delta Q$, respectively, the uniform distribution property of a hash function
could be evaluated using $|Q-50\%|$ and $\Delta Q$: the smaller they are, the
better the strict avalanche effect achieved. Additionally, the experimental
value of $Q$ is always equivalent to the value of $P$ if the test data (i.e.,
$N$ pairs of original and modified messages) used in the diffusion and
confusion test is re-used in the uniform distribution test. Such a result can
also be obtained through a simple reasoning: in
$P=\left(\begin{matrix}\sum_{i=1}^{N}B_{i}\end{matrix}\right)/(n\times m\times
N)\times 100\%$ and $Q=\left(\begin{matrix}\sum_{j=1}^{n\times
m}T_{j}\end{matrix}\right)/(n\times m\times N)\times 100\%$, both
$\begin{matrix}\sum_{i=1}^{N}B_{i}\end{matrix}$ and
$\begin{matrix}\sum_{j=1}^{n\times m}T_{j}\end{matrix}$ count the total number
of hash bits that are flipped over $N$ draws. Thus, $T$ or $Q$ alone is
insufficient for assessing the uniform distribution property of a hash
function, it should be considered along with $\Delta Q$.

The uniform distribution test on QHFM-$L$ is conducted as follows:

  1. (1)

Set $T_{j}=0$ for every bit position $j$ in the hash value.

  2. (2)

Randomly draw an article record from arXiv Dataset, take the abstract of this
article as the original message $msg_{0}$.

  3. (3)

Randomly flip a bit of $msg_{0}$ and then generate the modified message
$msg_{1}$.

  4. (4)

Compute the hash values of the two messages and get the digest pair
($H(msg_{0})$, $H(msg_{1})$); compare $H(msg_{0})$ with $H(msg_{1})$ bit by
bit, if $H(msg_{0})$ differs from $H(msg_{1})$ at the $j$th bit position, then
the value of $T_{j}$ is incremented by one.

  5. (5)

Repeat steps (2) to (4) $N$ times.

  6. (6)

Calculate $\overline{T}$, $Q$, $\Delta T$, and $\Delta Q$ from the obtained
data.

The data collected in step (2) is re-used for different instances of the
proposed hash scheme as well as for different hash properties (i.e., the
diffusion and confusion properties, the uniform distribution property, and the
collision resistant property). As a result, the experimental values of $P$ and
$Q$ for each instances are equal, which gives $N\times P=\overline{T}$ and
$|P-50\%|=|Q-50\%|$. On the other hand, for the existing schemes
Yang2021qwHash ; Yang2019qwHash ; Yang2018qwHash264 ; Yang18qwHash221 ;
Li2018qwHash ; Cao2018qwHash ; Yang2016qHash , the reported results of
$\overline{T}$ are not equivalent to the corresponding outcomes of $N\times
P$, this is probably because their input messages used in the uniform
distribution test are not the same as that used in the diffusion and confusion
test. Nevertheless, the reported values of $\overline{T}$ are generally close
to the corresponding results of $N\times P$.

Since the test results of $\Delta Q$ (or $\Delta T$) for the existing schemes
are unavailable for comparison, we collect reported data related to the
uniform distribution property from Refs. Yang2021qwHash ; Yang2019qwHash ;
Yang2018qwHash264 ; Yang18qwHash221 ; Li2018qwHash ; Cao2018qwHash ;
Yang2016qHash  as much as possible and list the corresponding results for the
proposed and existing schemes in Table 3, where “****.**, same” denotes a pair
of identical values, and “N/A” means “not available”. The values of $|Q-50\%|$
(in $\%$) for the existing schemes are deduced from the reported results of
$\overline{T}$: $|Q-50\%|=|\overline{T}/N\times 100\%-50\%|$, where $N=16383$
for Cao18-195 and $N=10000$ for others. Similar to Table 2, the values
presented in the 8th and 11th rows are results for the representative
instances of Yang21-296 and Yang18-221, respectively.

Table 3: Uniform-Distribution-Test Results for the Proposed and Existing QW-Based Hash Functions Instances (schemes) | $N\times P,\overline{T}\quad\;$ | $\Delta T$ | $\Delta Q(\%)$ | $|P-50\%|,|Q-50\%|(\%)$
---|---|---|---|---
QHFM-296 | 4996.96, same | 48.4334 | 0.4843 | 0.0304, same
QHFM-264 | 4994.95, same | 48.9253 | 0.4893 | 0.0505, same
QHFM-221 | 5001.42, same | 51.6083 | 0.5161 | 0.0142, same
QHFM-200 | 5001.03, same | 51.6897 | 0.5169 | 0.0103, same
QHFM-195 | 5003.03, same | 50.5134 | 0.5051 | 0.0303, same
QHFM-136 | 5003.90, same | 46.6002 | 0.4660 | 0.0390, same
QHFM-120 | 5007.62, same | 48.6068 | 0.4861 | 0.0762, same
Yang21-296 Yang2021qwHash  | 4995.41, 4998.1 | N/A | N/A | 0.0459, 0.019
Yang19-264 Yang2019qwHash  | 4987.89, 4996.6 | N/A | N/A | 0.1211, 0.034
Yang18-264 Yang2018qwHash264  | 5004.20, 5003.9 | N/A | N/A | 0.0420, 0.039
Yang18-221 Yang18qwHash221  | 5103.13, N/A | N/A | N/A | 1.0313, N/A
Li18-200 Li2018qwHash  | 4995.05, 4998.2 | N/A | N/A | 0.0495, 0.018
Cao18-195 Cao2018qwHash  | 10478.57, 6495.0 | N/A | N/A | 13.9600, 10.355
Yang16-128 Yang2016qHash  | 5022.61, 4973.5 | N/A | N/A | 0.2261, 0.265

It can be seen from the last column of Table 3 that the experimental values of
$P$ and $Q$ for QHFM-221, QHFM-195, and QHFM-136 together with QHFM-120 are
closer to their theoretical values than those for Yang18-221, Cao18-195, and
Yang16-128, respectively; and the values of $|Q-50\%|$ for the remaining
instances of QHFM are on a par with those for their peers. In addition, the
results of $\Delta Q$ for all instances of QHFM are very small, indicating
that the proposed hash scheme has a very good uniform distribution property.

To provide an intuitive description for this property of our scheme, we plot
the number of draws with flipped hash bit on every bit position of QHFM-195 in
Fig. 2, which suggests that the proposed scheme has a good resistance to
statistical attacks.

Figure 2: (Color online) Histogram of the 195-bit hash space, where $N=10000$.

###  4.4 Collision Resistance

The test data for the diffusion and confusion properties or the uniform
distribution property can also be used to analyzing the collision resistance
property, which is generally assessed in terms of two groups of indicators:
(1) the number of draws $W_{N}^{\text{\,e}}(\omega)$ (out of $N$ random
selections) on which the hash values of the original and modified messages
(i.e., $H(msg_{0})$ and $H(msg_{1})$) contain $\omega$ bytes with the same
value at the same location (here $\omega$ is also called the number of hits,
which satisfies $0\leq\omega\leq g\equiv\lceil(n\times m)/8\rceil$, and
$W_{N}^{\text{\,e}}(0)+W_{N}^{\text{\,e}}(1)+\dots+W_{N}^{\text{\,e}}(g)=N$);
and (2) the mean of the absolute difference per byte
$\bar{d}_{\text{byte}}^{\text{\,e}}$ between $H(msg_{0})$ and $H(msg_{1})$
over $N$ draws. If the results of $W_{N}^{\text{\,e}}(\omega)$ and
$\bar{d}_{\text{byte}}^{\text{\,e}}$ (the experimental values) are very close
to their theoretical values, then the related hash function could be regarded
as having a good property of collision resistance.

The number of hits on each draw can be obtained as follows: first, divide both
$H(msg_{0})$ and $H(msg_{1})$ into $g$ bytes (if $n\times m$ is not divisible
by 8, then add a prefix of $c=8-(n\times m)\bmod 8$ zeros to the hash values),
so that the two hash values can be expressed as $h=e_{1}\lVert e_{2}\lVert
e_{3}\lVert\dots\lVert e_{g}$ and $h^{\prime}=e^{\prime}_{1}\lVert
e^{\prime}_{2}\lVert e^{\prime}_{3}\lVert\dots\lVert e^{\prime}_{g}$,
respectively ($e_{j}$ and $e^{\prime}_{j}$ respectively represents the $j$th
byte of $h$ and $h^{\prime}$); second, compare $h$ and $h^{\prime}$ byte by
byte and compute $\omega$ according to

$\omega=\sum_{j=1}^{g}\delta\left[t(e_{j}),t(e^{\prime}_{j})\right],$ (9)

where $t(e_{j})$ is the decimal value of $e_{j}$ and $\delta[\cdot]$ is the
Kronecker delta function.

The theoretical value (denoted by $W_{N}^{\text{\,t}}(\omega)$) of
$W_{N}^{\text{\,e}}(\omega)$ is given by the product of $N$ and the
theoretical probability $P^{\text{\,t}}(\omega)$ that $\omega$ hits occur in
$(h,h^{\prime})$. Specifically, $P^{\text{\,t}}(\omega)$ is given by the
binomial distribution formula

$P^{\text{\,t}}(\omega)=\frac{g!}{\omega!(g-\omega)!}\left(\frac{1}{2^{8}}\right)^{\omega}\left(1\\!-\\!\frac{1}{2^{8}}\right)^{g-\omega},$
(10)

and the theoretical number of draws with $\omega$ hits is obtained by

$W_{N}^{\text{\,t}}(\omega)=\text{int}\left[N\times
P^{\text{\,t}}(\omega)\right],$ (11)

where $\text{int}[\cdot]$ denotes rounding a real number to its nearest
integer.

Since $P^{\text{\,t}}\equiv\\{P^{\text{\,t}}(\omega)\,|\,\omega=0,1,2,\dots
g\\}$ and
$P^{\text{\,e}}\equiv\\{W_{N}^{\text{\,e}}(\omega)/N\,|\,\omega=0,1,2,\dots
g\\}$ respectively describe the theoretical and experimental distributions of
$\omega$ (or simply hit distributions), the similarity or difference between
$\left\\{W_{N}^{\text{\,e}}(\omega)\,|\,\omega=0,1,\dots g\right\\}$
(hereafter, simply $\\{W_{N}^{\text{\,e}}(\omega)\\}$) and
$\left\\{W_{N}^{\text{\,t}}(\omega)\,|\,\omega=0,1,\dots g\right\\}$
(hereafter, simply $\\{W_{N}^{\text{\,t}}(\omega)\\}$) could be measured by
Kullback-Leibler divergence between $P^{\text{\,t}}$ and $P^{\text{\,e}}$,
i.e.,

$\displaystyle D_{\text{KL}}\left(P^{\text{\,e}}\lVert P^{\text{\,t}}\right)$
$\displaystyle=\sum_{\omega=0}^{g}P^{\text{\,e}}(\omega)\text{log}_{2}\left(\frac{P^{\text{\,e}}(\omega)}{P^{\text{\,t}}(\omega)}\right)$
(12)
$\displaystyle=\sum_{\omega=0}^{g}\frac{W_{N}^{\text{\,e}}(\omega)}{N}\text{log}_{2}\left(\frac{W_{N}^{\text{\,e}}(\omega)/N}{P^{\text{\,t}}(\omega)}\right);$

a smaller $D_{\text{KL}}\left(P^{\text{\,e}}\lVert P^{\text{\,t}}\right)$
indicates a closer similarity between $\\{W_{N}^{\text{\,e}}(\omega)\\}$ and
$\\{W_{N}^{\text{\,t}}(\omega)\\}$.

The absolute difference per byte between $h$ and $h^{\prime}$ is calculated by

$d_{\text{byte}}=\frac{1}{g}\sum_{j=1}^{g}\left|t(e_{j})-t(e^{\prime}_{j})\right|,$
(13)

and the theoretical value (denoted by $\bar{d}_{\text{byte}}^{\text{\,t}}$ )
of the mean of $d_{\text{byte}}$ (denoted by
$\bar{d}_{\text{byte}}^{\text{\,e}}$ ) over $N$ draws is
$\bar{d}_{\text{byte}}^{\text{\,t}}=85.33$ Yang2021qwHash .

The collision resistant test on QHFM-$L$ is performed with $N=10000$, and the
simulation results are shown in Table 4, where $W_{N}^{\text{\,e}}(4+)$
denotes the number of draws on which more than three hits occur in the hash
values of the original and modified messages, that is,
$W_{N}^{\text{\,e}}(4+)=N-[W_{N}^{\text{\,e}}(0)+W_{N}^{\text{\,e}}(1)+W_{N}^{\text{\,e}}(2)+W_{N}^{\text{\,e}}(3)]$.
One may notice that the sums of $W_{N}^{\text{\,t}}(\omega)$ (over all
$\omega$) for 264- and 128-bit hash schemes (or instances) are not equivalent
to $N$, this is due to the rounding operations performed on $N\times
P^{\text{\,t}}(\omega)$. The values of
$D_{\text{KL}}\left(P^{\text{\,e}}\lVert P^{\text{\,t}}\right)$ and
$\Delta\bar{d}_{\text{byte}}\equiv|\bar{d}_{\text{byte}}^{\text{\,e}}-\bar{d}_{\text{byte}}^{\text{\,t}}|$
for the existing schemes are deduced from the reported results of
$\\{W_{N}^{\text{\,e}}(\omega)\\}$ and $\bar{d}_{\text{byte}}^{\text{\,e}}$
(or the mean of $d_{\text{byte}}\times g$), respectively.

Table 4: Collision-resistance-test results for the proposed and existing QW-based hash functions Instances (schemes) | $\qquad\\{W_{N}^{\text{\,e}}(\omega)\\}$ | $\qquad\\{W_{N}^{\text{\,t}}(\omega)\\}$ | $D_{\text{KL}}\left(P^{\text{\,e}}\lVert P^{\text{\,t}}\right)$ | $\Delta\bar{d}_{\text{byte}}$
---|---|---|---|---
QHFM-296 | $8605,1312,81,2,0$ | $8652,1255,89,4,0$ | 0.000361 | 0.03
QHFM-264 | $8762,1159,74,5,0$ | $8788,1137,71,3,0$ | 0.000146 | 0.06
QHFM-221 | $8674,1230,93,3,0$ | $8962,984,52,2,0$ | 0.006711 | 2.48
QHFM-200 | $9071,895,34,0,0$ | $9068,889,42,1,0$ | 0.000302 | 0.03
QHFM-195 | $8066,1796,130,8,0$ | $9068,889,42,1,0$ | 0.069364 | 3.30
QHFM-136 | $9352,626,21,1,0$ | $9356,624,20,0,0$ | 0.000058 | 0.01
QHFM-120 | $9416,570,13,1,0$ | $9430,555,15,0,0$ | 0.000145 | 0.08
Yang21-296 Yang2021qwHash  | $8321,1547,110,22,0$ | $8652,1255,89,4,0$ | 0.008616 | 0.11
Yang19-264 Yang2019qwHash  | $9019,923,52,2,4$ | $8788,1137,71,3,0$ | 0.005647 | 4.43
Yang18-264 Yang2018qwHash264  | $8904,1026,68,2,0$ | $8788,1137,71,3,0$ | 0.000969 | 1.69
Yang18-221 Yang18qwHash221  | $9854,71,0,0,75$ | $8962,984,52,2,0$ | 0.188620 | N/A
Li18-200 Li2018qwHash  | $8982,989,25,4,0$ | $9068,889,42,1,0$ | 0.001689 | N/A
Cao18-195 Cao2018qwHash  | $16063,314,6,0,0$ | $14856,1456,69,2,0$ | 0.066791 | N/A
Yang16-128 Yang2016qHash  | $9367,617,16,0,0$ | $9393,589,17,0,0$ | 0.000151 | 1.88

The values of $D_{\text{KL}}\left(P^{\text{\,e}}\lVert P^{\text{\,t}}\right)$
indicate that the experimental result of hit distribution for QHFM-$L$ with
$L\geq 200$ has closer similarity to the theoretical distribution of $\omega$
than those for the existing ones with $L$-bit output length, and the Kullback-
Leibler divergence between $P^{\text{\,e}}$ and $P^{\text{\,t}}$ for QHFM-$L$
with $L<200$ is on a par with that for its peer (Cao18-195 or Yang16-128). As
for the average difference per byte in two hash values, the test results of
$\bar{d}_{\text{byte}}$ for QHFM-296, QHFM-264, and QHFM-136 (together with
QHFM-128) are closer to the theoretical value 85.33 than those for Yang21-296,
Yang19-264, and Yang18-128, and the differences between
$\bar{d}_{\text{byte}}^{\text{\,e}}$ and $\bar{d}_{\text{byte}}^{\text{\,t}}$
for the remaining instances of QHFM are very small. Therefore, the proposed
hash scheme has a good capability of collision resistance.

###  4.5 Resistance to Birthday Attacks

Since the proposed hash function has variable digest length, one can easily
obtain a QHFM instance that withstands birthday attacks by assigning
appropriate values (large enough) to the parameters $m$ and $n$ according to
the (cryptanalytic) hardware and software capabilities considered.

##  5 Time and Space Complexity Analysis

The hash value of an input message sent to a QW-based hash function can be
calculated by cascading three stages: (1) initializing the state of the
walker; (2) performing the underlying CAQW on a cycle according to the bit
values of the message; and (3) calculating the hash value from the resulting
probability distribution of the walker. The time and space complexity of the
proposed scheme or an existing one can thus be obtained by analyzing the
number of arithmetic operations taken by each stage of the related hashing
process. Since the hash value is computed classically (quantum transforms are
simulated by matrix multiplications), this section will concentrate on
classical complexity.

###  5.1 Time and space complexity of the proposed scheme

The quantum state of the walker after $t$ steps of CQWM ($t\geq 0$) can be
expressed as

$\Ket{\psi_{t}}=\sum_{x,j}A_{t}^{x,j}\Ket{x,j},$ (14)

where $A_{t}^{x,j}(x\in\mathbb{Z}_{n},j\in\mathbb{Z}_{8})$ is the amplitude of
the 2-term basis state $\ket{x,j}$ at time $t$, and the correspondence between
2-term and 4-term basis states is described by equation group (5). Before
($t=0$) and during ($t>0$) the walk, the state of the particle is identified
with these $8n$ amplitudes.

When $t=0$, the particle is in the state
$\Ket{\psi_{0}}=\text{cos}\alpha\Ket{0,2}+\text{sin}\alpha\Ket{0,3}$, which
gives $A_{0}^{0,2}=\text{cos}\alpha$, $A_{0}^{0,3}=\text{sin}\alpha$,
$A_{0}^{0,j}=0$ for $j\neq 2$ and $j\neq 3$, and $A_{0}^{x,j}=0$ for $x\neq
0$. Thus, the classical representation of the initial state $\Ket{\psi_{0}}$
can be specified using $8n$ assignments.

When $t>0$, if the $t$th message bit equals 0, the values of
$\\{A_{t}^{x,j}|x\in\mathbb{Z}_{n},j\in\mathbb{Z}_{8}\\}$ (hereafter, simply
$\\{A_{t}^{x,j}\\}$) are determined by $U^{(0)}=S(I_{n}\otimes
D^{(0)})(I_{4n}\otimes C^{(0)})$ and
$\\{A_{t-1}^{x,j}|x\in\mathbb{Z}_{n},j\in\mathbb{Z}_{8}\\}$ (hereafter, simply
$\\{A_{t-1}^{x,j}\\}$). For the sake of simplicity of notation, we denote
$C^{(0)}$ by $\bigl{(}\begin{smallmatrix}a_{0}&b_{0}\\\
c_{0}&d_{0}\end{smallmatrix}\bigr{)}$, then the action of this coin operator
on $\Ket{x,dr_{2},dr_{1},c}$ can be formulated as

$\displaystyle C^{(0)}\\!:\Ket{x,dr_{2},dr_{1},c}$
$\displaystyle\\!\to\\!\bar{c}(a_{0}\Ket{x,dr_{2},dr_{1},0}\\!+\\!c_{0}\Ket{x,dr_{2},dr_{1},1})+$
(15) $\displaystyle
c(b_{0}\Ket{x,dr_{2},dr_{1},0}\\!+\\!d_{0}\Ket{x,dr_{2},dr_{1},1}).$

Converting the 4-term states in expression (15) into 2-term states gives

$\displaystyle C^{(0)}:\Ket{x,j}\to$ $\displaystyle\left[(j\bmod 2)\oplus
1\right](a_{0}\Ket{x,j}+c_{0}\Ket{x,j+1})+$ (16) $\displaystyle(j\bmod
2)(b_{0}\Ket{x,j-1}+d_{0}\Ket{x,j}),$

where $j+1$ and $j-1$ are both calculated using modular arithmetic under
modulus 8. Combing expressions (16), (6), and (8), one can obtain the action
of $U^{(0)}$ on each 2-term basis state as well as on $\Ket{\psi_{t-1}}$ and
then deduce the relation between $\\{A_{t}^{x,j}\\}$ and
$\\{A_{t-1}^{x,j}\\}$. Specifically, the actions of $U^{(0)}$ on the
components $A_{t-1}^{x,j}\Ket{x,j}$ of $\Ket{\psi_{t-1}}$ are

$\displaystyle A_{t-1}^{x,0}\Ket{x,0}$ $\displaystyle\to
a_{0}A_{t-1}^{x,0}\Ket{x+1,4}+c_{0}A_{t-1}^{x,0}\Ket{x-1,1},$ (17)
$\displaystyle A_{t-1}^{x,1}\Ket{x,1}$ $\displaystyle\to
b_{0}A_{t-1}^{x,1}\Ket{x+1,4}+d_{0}A_{t-1}^{x,1}\Ket{x-1,1},$ $\displaystyle
A_{t-1}^{x,2}\Ket{x,2}$ $\displaystyle\to
a_{0}A_{t-1}^{x,2}\Ket{x+1,6}+c_{0}A_{t-1}^{x,2}\Ket{x-1,3},$ $\displaystyle
A_{t-1}^{x,3}\Ket{x,3}$ $\displaystyle\to
b_{0}A_{t-1}^{x,3}\Ket{x+1,6}+d_{0}A_{t-1}^{x,3}\Ket{x-1,3},$ $\displaystyle
A_{t-1}^{x,4}\Ket{x,4}$ $\displaystyle\to
a_{0}A_{t-1}^{x,4}\Ket{x-1,0}+c_{0}A_{t-1}^{x,4}\Ket{x+1,5},$ $\displaystyle
A_{t-1}^{x,5}\Ket{x,5}$ $\displaystyle\to
b_{0}A_{t-1}^{x,5}\Ket{x-1,0}+d_{0}A_{t-1}^{x,5}\Ket{x+1,5},$ $\displaystyle
A_{t-1}^{x,6}\Ket{x,6}$ $\displaystyle\to
a_{0}A_{t-1}^{x,6}\Ket{x-1,2}+c_{0}A_{t-1}^{x,6}\Ket{x+1,7},$ $\displaystyle
A_{t-1}^{x,7}\Ket{x,7}$ $\displaystyle\to
b_{0}A_{t-1}^{x,7}\Ket{x-1,2}+d_{0}A_{t-1}^{x,7}\Ket{x+1,7},$

where $x\pm 1$ are calculated using modular arithmetic under modulus $n$.
Summing up the transformed results on the right side, one can observe that the
amplitudes of the walker being at position $x$ at time $t-1$ contribute
$a_{0}A_{t-1}^{x,0}+b_{0}A_{t-1}^{x,1}$ to $A_{t}^{x+1,4}$,
$c_{0}A_{t-1}^{x,0}+d_{0}A_{t-1}^{x,1}$ to $A_{t}^{x-1,1}$, and
$a_{0}A_{t-1}^{x,2}+b_{0}A_{t-1}^{x,3}$ to $A_{t}^{x+1,6}$, etc.; here
$A_{t}^{x\pm 1,j}$ is the amplitude of $\Ket{x\pm 1\pmod{n},j}$ at time $t$.
Moreover, the amplitudes of being at position $x$ at time $t-1$ only
contribute to $A_{t}^{x-1,j}$ with $j\leq 3$ and to $A_{t}^{x+1,j}$ with
$j\geq 4$; conversely, the former 4 amplitudes (with $j\leq 3$) at an
arbitrary position $x$ are contributed by the amplitudes at position
$x+1\pmod{n}$, while the latter 4 amplitudes (with $j\geq 4$) at position $x$
are contributed by those at position $x-1\pmod{n}$. As a result, each
amplitude of being at position $x$ at time $t$, denoted by $A_{t}^{x,j}$, is
only contributed by the amplitudes of being at a single position
($x+1\pmod{n}$ or $x-1\pmod{n}$) at time $t-1$. Thus, the relation between
$\\{A_{t}^{x,j}\\}$ and $\\{A_{t-1}^{x,j}\\}$ after a step of QW1M can be
expressed as follows:

$\displaystyle
A_{t}^{x,0}=a_{0}A_{t-1}^{x+1,4}+b_{0}A_{t-1}^{x+1,5},A_{t}^{x,1}=c_{0}A_{t-1}^{x+1,0}+d_{0}A_{t-1}^{x+1,1},$
(18) $\displaystyle
A_{t}^{x,2}=a_{0}A_{t-1}^{x+1,6}+b_{0}A_{t-1}^{x+1,7},A_{t}^{x,3}=c_{0}A_{t-1}^{x+1,2}+d_{0}A_{t-1}^{x+1,3},$
$\displaystyle
A_{t}^{x,4}=a_{0}A_{t-1}^{x-1,0}+b_{0}A_{t-1}^{x-1,1},A_{t}^{x,5}=c_{0}A_{t-1}^{x-1,4}+d_{0}A_{t-1}^{x-1,5},$
$\displaystyle
A_{t}^{x,6}=a_{0}A_{t-1}^{x-1,2}+b_{0}A_{t-1}^{x-1,3},A_{t}^{x,7}=c_{0}A_{t-1}^{x-1,6}+d_{0}A_{t-1}^{x-1,7}.$

Analogously, if the $t$th message bit equals 1, the values of
$\\{A_{t}^{x,j}\\}$ are determined by $U^{(1)}=S(I_{n}\otimes
D^{(1)})(I_{4n}\otimes C^{(1)})$ and $\\{A_{t-1}^{x,j}\\}$. We denote
$C^{(1)}$ by $\bigl{(}\begin{smallmatrix}a_{1}&b_{1}\\\
c_{1}&d_{1}\end{smallmatrix}\bigr{)}$, then the action of $C^{(1)}$ on
$\Ket{x,j}$ is

$\displaystyle C^{(1)}:\Ket{x,j}\to$ $\displaystyle\left[(j\bmod 2)\oplus
1\right](a_{1}\Ket{x,j}+c_{1}\Ket{x,j+1})+$ (19) $\displaystyle(j\bmod
2)(b_{1}\Ket{x,j-1}+d_{1}\Ket{x,j}).$

A Combination of expressions (19), (7), and (8) gives the actions of $U^{(1)}$
on the components of $\Ket{\psi_{t-1}}$:

$\displaystyle A_{t-1}^{x,0}\Ket{x,0}$ $\displaystyle\to
a_{1}A_{t-1}^{x,0}\Ket{x+1,4}+c_{1}A_{t-1}^{x,0}\Ket{x-1,1},$ (20)
$\displaystyle A_{t-1}^{x,1}\Ket{x,1}$ $\displaystyle\to
b_{1}A_{t-1}^{x,1}\Ket{x+1,4}+d_{1}A_{t-1}^{x,1}\Ket{x-1,1},$ $\displaystyle
A_{t-1}^{x,2}\Ket{x,2}$ $\displaystyle\to
a_{1}A_{t-1}^{x,2}\Ket{x-1,0}+c_{1}A_{t-1}^{x,2}\Ket{x+1,5},$ $\displaystyle
A_{t-1}^{x,3}\Ket{x,3}$ $\displaystyle\to
b_{1}A_{t-1}^{x,3}\Ket{x-1,0}+d_{1}A_{t-1}^{x,3}\Ket{x+1,5},$ $\displaystyle
A_{t-1}^{x,4}\Ket{x,4}$ $\displaystyle\to
a_{1}A_{t-1}^{x,4}\Ket{x+1,6}+c_{1}A_{t-1}^{x,4}\Ket{x-1,3},$ $\displaystyle
A_{t-1}^{x,5}\Ket{x,5}$ $\displaystyle\to
b_{1}A_{t-1}^{x,5}\Ket{x+1,6}+d_{1}A_{t-1}^{x,5}\Ket{x-1,3},$ $\displaystyle
A_{t-1}^{x,6}\Ket{x,6}$ $\displaystyle\to
a_{1}A_{t-1}^{x,6}\Ket{x-1,2}+c_{1}A_{t-1}^{x,6}\Ket{x+1,7},$ $\displaystyle
A_{t-1}^{x,7}\Ket{x,7}$ $\displaystyle\to
b_{1}A_{t-1}^{x,7}\Ket{x-1,2}+d_{1}A_{t-1}^{x,7}\Ket{x+1,7}.$

Thus, the relation between $\\{A_{t}^{x,j}\\}$ and $\\{A_{t-1}^{x,j}\\}$ after
a step of QW2M can be expressed as

$\displaystyle
A_{t}^{x,0}=a_{1}A_{t-1}^{x+1,2}+b_{1}A_{t-1}^{x+1,3},A_{t}^{x,1}=c_{1}A_{t-1}^{x+1,0}+d_{1}A_{t-1}^{x+1,1},$
(21) $\displaystyle
A_{t}^{x,2}=a_{1}A_{t-1}^{x+1,6}+b_{1}A_{t-1}^{x+1,7},A_{t}^{x,3}=c_{1}A_{t-1}^{x+1,4}+d_{1}A_{t-1}^{x+1,5},$
$\displaystyle
A_{t}^{x,4}=a_{1}A_{t-1}^{x-1,0}+b_{1}A_{t-1}^{x-1,1},A_{t}^{x,5}=c_{1}A_{t-1}^{x-1,2}+d_{1}A_{t-1}^{x-1,3},$
$\displaystyle
A_{t}^{x,6}=a_{1}A_{t-1}^{x-1,4}+b_{1}A_{t-1}^{x-1,5},A_{t}^{x,7}=c_{1}A_{t-1}^{x-1,6}+d_{1}A_{t-1}^{x-1,7}.$

Relations (18) and (21) show that, given the amplitudes
$\\{A_{t-1}^{x_{0},j}|j\in\mathbb{Z}_{8}\\}$ of being at a fixed position
$x_{0}$ at time $t-1$, the values of the amplitudes
$\\{A_{t}^{x_{0},j}|j\in\mathbb{Z}_{8}\\}$ of being at $x_{0}$ at time $t$ can
be calculated using 16 multiplications and 8 additions, which means all
amplitudes at each time step of CQWM on a cycle with $n$ nodes can be obtained
using $O(n)$ basic arithmetic operations. To perform these operations, one
needs to store the old (or the initial) $8n$ amplitudes and their 8 possible
coefficients ($a_{1},a_{2},\dots,d_{1},d_{2}$) to calculate the new $8n$
amplitudes, and the values of both old and new amplitudes are refreshed at
each time step. If the input message $msg$ is a binary string of $M$ bits,
then the values of $\\{A_{M}^{x,j}\\}$ can be obtained using $O(Mn)$ basic
operations with $O(n)$ memory space. Finally, the hash value is computed from
$\\{A_{M}^{x,j}\\}$ using $O(n)$ multiplications and $O(n)$ modulo operations
with $O(n)$ space. Thus, the time and space complexity of QHFM with input
length $M$ are $O(Mn)$ and $O(n)$, respectively.

In particular, if one wants to obtain an $L$-bit hash value ($L$ is a multiple
of $m$) of $msg$ using QHFM-$L$, then the cycle utilized by QHFM-$L$ has
$L/m=O(L)$ nodes (here $m$ is constant with respect to the input length $M$);
in this case, the hash value is produced with $O(ML)$ time and $O(L)$ space.

###  5.2 Time and space complexity comparison of QW-based hash schemes

In a similar way, one can deduce the time and space complexity of the existing
QW-based hash functions Yang2021qwHash ; Yang2019qwHash ; Yang2018qwHash264 ;
Yang18qwHash221 ; Li2018qwHash ; Cao2018qwHash ; Yang2016qHash ; Li2013qwHash
with respect to the input and output length. To facilitate discussion, we
divide the existing schemes into four groups: (1) the hash functions based on
one-dimensional one-particle quantum walks Yang2021qwHash ; Yang2019qwHash ;
Yang2018qwHash264 ; Yang18qwHash221 , (2) the hash function based on two-
dimensional one-particle quantum walks Li2018qwHash , (3) the hash function
based on quantum walks on Johnson graphs Cao2018qwHash , and (4) the hash
functions based on one-dimensional two-particle quantum walks Yang2016qHash ;
Li2013qwHash .

Again, suppose all schemes produce hash values of length $L$. In this case,
the schemes in group (1) utilize a cycle with $O(L)$ nodes, and the amplitudes
of the particle being at each node at time $t$ can be calculated from the
amplitudes of being at the two neighbors of this node at time $t-1$ (possibly
calculated from the amplitudes of being at a single neighbor or remain
unchanged during broken-line quantum walks Yang2021qwHash ) using constant
number of (basic arithmetic) operations. In group (2), Li18-200 utilizes
cycles of length $O(\sqrt{L})$ in two-dimensional space, which lead to $O(L)$
positions for the walker, and the amplitudes of being at position $(x,y)$ at
time $t$ can be calculated from the amplitudes of being at $(x\pm 1,y\pm 1)$
at time $t-1$ using constant number of operations. In group (3), Cao18-195
utilizes a Johnson graph $J(n,1)$ with $n=O(L)$ nodes, and the amplitudes of
the particle being at each node at time $t$ can be calculated from the
amplitudes of being at the remaining $n-1$ nodes at time $t-1$ using $O(L)$
operations. Similar to group (2), schemes in group (4) also utilize a cycle
with $O(\sqrt{N})$ nodes, which leads to $O(L)$ position pairs for the two
particles, and the amplitudes of the first and second particles being
respectively at nodes $x$ and $y$ at time $t$ can be calculated from the
amplitudes of the two particles being respectively at $x\pm 1$ and $y\pm 1$ at
time $t-1$ using constant number of operations.

Thus, except for Cao18-195 Cao2018qwHash , which performs $O(L^{2})$
operations at each time step, the existing QW-based schemes Yang2021qwHash ;
Yang2019qwHash ; Yang2018qwHash264 ; Yang18qwHash221 ; Li2018qwHash ;
Yang2016qHash ; Li2013qwHash  take $O(L)$ time to calculate the amplitudes at
time $t$ from the amplitudes at time $t-1$. If the input message is of bit-
length $M$, then the resulting amplitudes at time $M$ can be obtained with
$O(ML^{2})$ and $O(ML)$ operations in Cao18-195 and the remaining schemes,
respectively. After that, for all these schemes Yang2021qwHash ;
Yang2019qwHash ; Yang2018qwHash264 ; Yang18qwHash221 ; Li2018qwHash ;
Cao2018qwHash ; Yang2016qHash ; Li2013qwHash , the hash value is computed from
the resulting $O(L)$ amplitudes with $O(L)$ operations. Therefore, the time
complexities of Cao18-195 and the other QW-based hash schemes are $O(ML^{2})$
and $O(ML)$, respectively. Since each amplitude at time $t$ is a linear
combination of the amplitudes of $O(1)$ or $O(L)$ positions at time $t-1$ in
all schemes, and the resulting probability distribution takes $O(L)$ space as
well, the space complexities of the existing schemes all equals $O(L)$.

As a result, the proposed scheme has the same time and space complexity as the
existing QW-based hash schemes except Cao18-195, whose time complexity is
slightly greater than that of the other schemes, including the proposed one.

##  6 Conclusion

In this paper, a new hash function QHFM based on quantum walks with one- and
two-step memory on circles is constructed, whose statistical properties as
well as time and space complexity are evaluated and compared with the existing
QW-based hash functions.

Unlike the existing analyses of hash schemes based on quantum walks without
memory, where a single indicator $\overline{T}$ is used to evaluate the
uniform distribution property, we adopted an additional indicator $\Delta Q$
to assess this property, since $\overline{T}$ alone is closely related to $P$,
implying that it also suggests the diffusion and confusion properties. In the
collision resistance analysis, we use Kullback-Leibler divergence to evaluate
the similarity between the experimental and theoretical distributions of
$\omega$, so that the difference between $\\{W_{N}^{\text{\,e}}(\omega)\\}$
and $\\{W_{N}^{\text{\,t}}(\omega)\\}$ can be indicated by a single number.

The analysis results show that QHFM has near-ideal statistical performance and
takes no more time and space than its peers, and they also suggest that
alternately running two quantum walks differing in more than one respects,
including coin operator and memory length, can also yield good hash functions.
Thus, it is unnecessary to restrict the component parts of a controlled
alternate quantum walk to a single kind of walk (equipped with controlled
coins). In future work, we will explore the possibility of combining two
quantum walks with more differences and investigate the effect of those
differences on the performance of the resulting hash function.

###### Acknowledgements.

The authors gratefully acknowledge the financial support from the China
Postdoctoral Science Foundation under Grant No. 2021M691148, the Hubei
Provincial Science and Technology Major Project of China under Grant No.
2020AEA011, and the Key Research $\&$ Development Plan of Hubei Province of
China under Grant No. 2020BAB100.

## References

  * (1) Bennett, C.H., Brassard, G., Crépeau, C.: Generalized privacy amplification. IEEE Trans. Inf. Theory 41(6), 1915-1923 (1995).
  * (2) Menezes, A.J., van Oorschot, P.C., Vanstone, S.A.: Handbook of Applied Cryptography, 1st ed. CRC press, Boca Raton, FL (1996)
  * (3) Ablayev, F., Ablayev, M., Vasiliev, A., Ziatdinov, M.: Quantum fingerprinting and quantum hashing. Computational and cryptographical aspects. Balt. J. Mod. Comput. 4(4), 860 (2016).
  * (4) Ablayev, F., Ablayev, M., Vasiliev, A.: On the balanced quantum hashing. J. Phys.: Conf. Ser. 681(1), 012019 (2016).
  * (5) Vasiliev, A.: Quantum hashing for finite abelian groups. Lobachevskii Journal of Mathematics 37(6), 753-757 (2016).
  * (6) Ziatdinov, M.: From graphs to keyed quantum hash functions. Lobachevskii Journal of Mathematics 37(6), 705-712 (2016).
  * (7) Ablayev, F.M., Vasiliev, A.V.: Cryptographic quantum hashing. Laser Phys. Lett. 11(2), 025202 (2013).
  * (8) Ablayev, F.M., Ziatdinov, M.T.: Universal Hash Functions from Quantum Procedures. Uchenye Zapiski Kazanskogo Universiteta. Seriya Fiziko-Matematicheskie Nauki 162(3), 259-268 (2020).
  * (9) Yang, Y.G., Dong, J.R., Yang, Y.L., Zhou, Y.H., Shi, W.M.: Usefulness of decoherence in quantum-walk-based hash function. Int. J. Theor. Phys., 1-13 (2021).
  * (10) Yang, Y.G., Bi, J.L., Li, D., Zhou, Y.H., Shi, W.M.: Hash function based on quantum walks. Int. J. Theor. Phys. 58(6), 1861-1873 (2019).
  * (11) Yang, Y.G., Bi, J.L., Chen, X.B., Yuan, Z., Zhou, Y.H., Shi, W.M.: Simple hash function using discrete-time quantum walks. Quantum Inf. Process. 17(8), 189 (2018).
  * (12) Yang, Y.G., Zhang, Y.C., Xu, G., Chen, X.B., Zhou, Y.H., Shi, W.M.: Improving the efficiency of quantum Hash function by dense coding of coin operators in discrete-time quantum walk. Sci. China Phys. Mech. 61(3), 030312 (2018).
  * (13) Li, D., Yang, Y.G., Bi, J.L., Yuan, J.B., Xu, J.: Controlled alternate quantum walks based quantum hash function. Sci. Rep. 8(1), 1-7 (2018).
  * (14) Cao, W.F., Zhang, Y.C., Yang, Y.G., Li, D., Zhou, Y.H., Shi, W.M.: Constructing quantum Hash functions based on quantum walks on Johnson graphs. Quantum Inf. Process. 17(7), 156 (2018).
  * (15) Yang, Y.G., Xu, P., Yang, R., Zhou, Y.H., Shi, W.M.: Quantum Hash function and its application to privacy amplification in quantum key distribution, pseudo-random number generation and image encryption. Sci. Rep. 6, 19788 (2016).
  * (16) Li, D., Zhang, J., Guo, F.Z., Huang, W., Wen, Q.Y., Chen, H.: Discrete-time interacting quantum walks and quantum Hash schemes. Quantum Inf. Process. 12(3), 1501-1513 (2013).
  * (17) Li, D., Mc Gettrick, M., Yang, Y.G., Xu, J., Wang, Y.: Quantum walks with memory provided by parity of memory. Int. J. Theor. Phys. 59(6), 1934-1943 (2020).
  * (18) Dai, W., Yuan, J., Li, D.: Discrete-Time Quantum walk with memory on the cayley graph of the dihedral group. Int. J. Theor. Phys. 59(1), 10-28 (2020).
  * (19) Zhou, Q., Lu, S.F.: One-dimensional quantum walks with two-step memory. Quantum Inf. Process. 18(12), 359 (2019).
  * (20) Li, D., Mc Gettrick, M., Gao, F., Xu, J., Wen, Q.Y.: Generic Quantum Walks with Memory on Regular Graphs. Phys. Rev. A 93(4), 042323 (2016).
  * (21) Mc Gettrick, M., Miszczak, J.A.: Quantum walks with memory on cycles. Physica A 399, 163-170 (2014).
  * (22) Mc Gettrick, M.: One dimensional quantum walks with memory. Quantum Inf. Comput. 10(5), 509 524 (2010).


# DeepPropNet - A Recursive Deep Propagator Neural Network for Learning
Evolution PDE Operators

Lizuo Liu, Wei Cai Department of Mathematics, Southern Methodist University,
Dallas, TX 75252. email<EMAIL_ADDRESS>Date: February 27, 2022.

###### Abstract

In this paper, we propose a deep neural network approximation to the evolution
operator for time dependent PDE systems over long time period by recursively
using one single neural network propagator, in the form of POD-DeepONet with
built-in causality feature, for a small time interval. The trained DeepPropNet
of moderate size is shown to give accurate prediction of wave solutions over
the whole time interval.

##  1 The Problem of Learning PDE Evolution Operator

Consider a second order evolution system

$\displaystyle u_{tt}$ $\displaystyle=\mathcal{L}u+f(x,t),$ (1) $\displaystyle
u(x,0)$ $\displaystyle=u_{0}(x),\quad x\in R^{1},$ $\displaystyle u_{t}(x,0)$
$\displaystyle=v_{0}(x).$ (2)

where $u$ could be a scalar or vector. For a model problem, we will consider
the inhomogeneous scalar wave (d’Alembert) equation

$\frac{\partial^{2}u}{\partial
t^{2}}-c^{2}\left(x,t\right)\frac{\partial^{2}u}{\partial x^{2}}=f(x,t),\quad
x\in\mathbb{R},\quad t\geq 0,$ (3)

with the source term $f$ that is compactly supported on a bounded space-time
domain $Q=S\times\left[0,T\right]$ where $S\subset\mathbb{R}$. This means that
the source term $f(x,t)$ differs from zero only on $S$ and operates for the
limited time interval $\left[0,T\right]$. When $c\left(x,t\right)=\text{const,
t}$he solution to the Cauchy problem $\left(\ref{1dwave}\right)$ is given by:

$u(x,t)=\frac{1}{2}(u_{0}(x-ct)+u_{0}(x+ct))+\frac{1}{2c}\int_{x-ct}^{x+ct}v_{0}(\xi)d\xi+\frac{1}{2c}\int_{0}^{t}\int_{x-c(t-\tau)}^{x+c(t+\tau)}f(\xi,\tau)d\xi
d\tau.$ (4)

The solution $u(x,t)$ can be viewed through an evolution operator
$P(u_{0},v_{0,}f(x,s),0\leq s\leq t)$, which maps the initial conditions and
the source term into the solution. Learning such a map between functions has
been actively studied recently with various types of operator learning
methods, including DeepONet [4] and Fourier Neural Operator [3]. The focus of
this paper is to find an efficient way to learn this evolution operator with
moderate size neural network for large time $t$. For lack of a precise term,
borrowing the term from quantum mechanics for the Green’s function propagator
[5], we shall name the operator $P$ as the propagator for the evolution
system. From (4), it is clear that if we like to train a neural network
operator for large time t, the size of the network will grow for increasing
time t. And, the amount of information to be input into a network will
increase dramatically as t grows as well. We will propose a recursive
propagator formulation for the evolution operator network.

First, the solution time interval $[0,T]$ will be divided into $N$ smaller
subintervals

$t_{0}=0<t_{1}<\cdots<t_{i}<\cdots t_{N},t_{i}=i\Delta t,\quad\Delta t=T/N,$
(5)

and for $t_{i}\leq t\leq t_{i+1},$ the solution is given by the propagator
with initial condition of the solution $u_{i}$ and its velocity
$v_{i}=\overset{.}{u}_{i}$, i.e.,

$u(x,t)=P(u_{i},v_{i},f(x,s),t_{i}\leq s\leq t),$ (6)

where the initial condition $u_{i}$ would have been given by the propagator
for the time block $t_{i-1}\leq t\leq t_{i}$.

The propagators in (6) will be approximated by a single neural network in the
form of DeepONet structure [2] with modification for time causality, namely,

$P(u_{i},v_{i},f(x,s),t_{i}\leq s\leq t)\sim
P_{\theta}(u_{i},v_{i},f(x,s),t_{i}\leq s\leq t),0\leq i\leq N-1.$ (7)

Therefore, the propagator $P_{\theta}(u_{0}(x),v_{0}(x),f(x,t))$ will be
trained to map the initial and force function data into the solution
$u(x,t),t_{0}\leq t\leq t_{1}$. Moreover, this same propagator $P_{\theta}$
will be trained to approximate the solution for time periodic $t_{1}\leq t\leq
t_{2}$ where the initial condition at $t_{1}$ can be computed with the
propagator for the previous time interval $[t_{0},t_{1}]$. This procedure will
be used recursively until we have trained the same propagator for the last
time interval $[t_{N-1},t_{N}]$. Taken all together, we arrive at a propagator
neural network for the whole time interval $[0,T]$ where the building block is
the single propagator of moderate size
$P_{\theta}(u_{0}(x),v_{0}(x),f(x,t)),t_{0}\leq t\leq t_{1})$. By controlling
the size of the $\Delta t$, the size of this propagator can be easily
controlled for efficiency as well as accuracy. As the evolution PDE system has
to observe the causality of the physical system, the DeepONet framework will
be modified to include the causality, a previously proposed causality DeepONet
in the study of dynamics system for modeling building response to seismic
waves [4] will be used for this purpose.

The rest of the paper is organized as follows. In section 2, we will review
the DeepONet [4] with time causality and extension with proper orthogonal
decomposition (POD) approach for efficient treatment of spatial dependence of
the solution. Section 3 will give the algorthmic outline of the DeepPropNet
and numerical results of the DeepPropNet will be presented in Section 4.
Finally, Section 5 gives a conclusion and some future work.

##  2 DeepONet with time causality and spatial POD

A causality DeepONet was proposed in [4] to handle the time casuality in
dynamics system and was shown to be very effective to predict the seismic
response of building. Here, we will just present the final form and for
details please refer to [4].

Causality-DeepONet: A DNN representation of an operator $\mathcal{G}(f)(t)$
for any continuous function $f(t)$ with retarded response for $t\in
K_{2}=[0,T]\subset$ $\mathbb{R}$ is given as

$\mathcal{G}(f)(t)\sim\sum_{k=1}^{N}\sum_{i=1}^{M}c_{i}^{k}\sigma_{b}\left(\sum_{j=1}^{\left\lceil\frac{t}{h}\right\rceil}\xi_{i,m-\left\lceil\frac{t}{h}\right\rceil+j}^{k}f\left(s_{j}\right)\right)\cdot\sigma_{trk}\left(\omega_{k}\cdot
t+\zeta_{k}\right),$ (8)

where $\\{s_{j}\\}_{j=1}^{m}\subset K_{1}=[0,T]\subset\mathcal{X},$coefficents
$c_{i}^{k},\xi_{ij}^{k},\omega_{k},\zeta_{k}$-all independent of continuous
functions $u$ $\in V\subset C(K_{1})$ and $t$.

To handle the spatial dependence of solution $u(x,t)$ for the evolution
system, we will adopt the idea of the POD-DeepONet[6, 7] which assumes there
is a set of global basis for the targeted output which could be found by SVD,
and the trunk net will be replaced by these basis, then the neural network is
learning the mapping between values at sensors to the singular values if we
consider using SVD as an example. Thus it may lose the ability to approximate
arbitrary values of the operator in the computational domain given the values
at sensors. The philosophy of POD-DeepONet is to find a set of basis globally,
thus compress the memory and computation needed, which is the crux for
learning the operator of the high dimensional problem$\left(d\gg
O\left(1\right)\right)$.

Lu et al.[7] discussed the POD-DeepONets and proposed a time marching scheme
for Fourier Neural Operator specifically to reduce the dimension of output,
and they also proposed the modified DeepONet with feature expansion by feeding
historical states of signal to trunk net as features. Bhattacharya et al.[6]
proposed that they reduce the dimension of input and output by PCA and then
learning the mapping between the reduced space by a neural network, and they
also discussed the approximation quality needed to get good mappings. Meuris
et al.[1] presented a procedure that learning the data-driven basis functions
harnessing the DeepONet machinery and the learned basis function will be used
by classical methods as custom basis to achieving high accurary for
computation of arbitrary complex domain.

As for a problem with causality, The POD-DeepONet’s philosophy guides us
either do the SVD for the whole outputs regardless the difference of temporal
variables and spatial variables, or do specific SVD timestep by timestep to
keep the causality for which we need to keep the basis for all time steps, or
consider solution at each time step for each case as an independent target and
find a common basis w.r.t spatial dimensions for all time steps & all cases.
These ideas either burns the high memory cost, or destroy the causality in a
brutal-force way. Therefore, a modification of POD-DeepONet for problem with
causality is crucial. We follow the idea of POD-DeepONet but only construct
the basis of spatial domain explicitly and utilize the causality DeepOnet to
handle the temporal-dependent coefficients of each spatial basis.

In the following, we consider the 1d wave equation (3 and assume the right
hand side has form

$f\left(x,t\right)=a_{0}\left(t\right)+\sum_{n=1}^{N}a_{n}\left(t\right)\cos\left(2n\pi
x\right)+b_{n}\left(t\right)\sin\left(2n\pi x\right),$ (9)

and the solution has a similar form

$u\left(x,t\right)=\psi_{0}\left(t\right)+\sum_{m=1}^{M}\psi_{m}\left(t\right)\cos\left(2m\pi
x\right)+\phi_{m}\left(t\right)\sin\left(2m\pi x\right),$ (10)

Next, we modify the Causality-DeepONet $\left(\ref{operator_appro}\right)$ by

$\mathcal{G}(\overrightarrow{f})(x,t)\sim\overrightarrow{\sigma}_{br}\left(\overrightarrow{f}\right)\odot\overrightarrow{\sigma}_{trk}\left(t\right)\cdot\overrightarrow{\sigma}_{basis}\left(x\right),$
(11)

where $\odot$ is the elementwise multiplication, $\cdot$ is the inner product,

$\overrightarrow{f}\left(s\right)=\left[a_{0}\left(s\right),a_{1}\left(s\right),b_{1}\left(s\right),\ldots,a_{N}\left(s\right),b_{N}\left(s\right)\right],$
$\displaystyle\begin{split}\overrightarrow{\sigma}_{br,i}\left(\overrightarrow{f}\right)={}&\sigma_{b}\left(\sum_{j=1}^{\left\lceil\frac{t}{h}\right\rceil}\xi_{i,m-\left\lceil\frac{t}{h}\right\rceil+j}^{0}a_{0}\left(s_{j}\right)\right.\\\
&\left.+\sum_{n=1}^{N}\sum_{j=1}^{\left\lceil\frac{t}{h}\right\rceil}\left(\xi_{i,m-\left\lceil\frac{t}{h}\right\rceil+j}^{n,a}a_{n}\left(s_{j}\right)+\xi_{i,m-\left\lceil\frac{t}{h}\right\rceil+j}^{n,b}b_{n}\left(s_{j}\right)\right)\right),\end{split}$
(12)
$\displaystyle\begin{split}\overrightarrow{\sigma}_{trk,i}\left(t\right)={}&\sigma_{trk}\left(\omega_{i}\cdot
t+\zeta_{i}\right),\end{split}$ (13)
$\displaystyle\begin{split}\overrightarrow{\sigma}_{basis}\left(x\right)={}&\\{1,\cos(2\pi
x),\sin(2\pi x),\cdots,\cos(2\pi Nx),\sin(2\pi Nx)\\}.\end{split}$ (14)

The activation functions of $\sigma_{b}$ and $\sigma_{trk}$ are ReLU. Note the
modes number $M$ in the CPOD-DeepONet(Causality POD-DeepONet)
$\left(\ref{1dsol}\right)$ need to be greater than or equal to the number of
given modes $N$ of right hand side in equation $\left(\ref{1drhs}\right)$. The
algorithm of the CPOD-DeepONet is shown in algorithm 1. To contain the memory
efficiency and following the fact that equation $\left(\ref{1dsol}\right)$ has
separation of variables, we will do outer product rather than computing point
by points as shown in the last line of algorithm. This is one of the crux to
handle higher dimensional problem.

Input : $t$, ${x}$,
$\left[a_{0}\left(s\right),a_{1}\left(s\right),b_{1}\left(s\right),\ldots,a_{N}\left(s\right),b_{N}\left(s\right)\right]$

Output : $u\left({x},t\right)$

$y_{t}\leftarrow\text{CDeepONet}\left(t,\left[a_{0}\left(s\right),a_{1}\left(s\right),b_{1}\left(s\right),\ldots,a_{N}\left(s\right),b_{N}\left(s\right)\right]\right)$
;

/* CDeepONet is the causality DeepONet. The output $y_{t}$ has size of
$[\text{batchsize}\times\text{Nt}\times\text{NBx}]$ */

$y_{x}\leftarrow\left[1,\cos\left(2\pi{x}\right),\sin\left(2\pi{x}\right),\ldots\cos\left(2M\pi{x}\right),\sin\left(2M{x}\right)\right]$;

/* The output $y_{x}$ has size of
$[\text{batchsize}\times\text{Nx}\times\text{NBx}]$ */

$u\left(x,t\right)\leftarrow\text{Einsum}\left(\text{'bxn,btn}\to\text{
bxt'},y_{x},y_{t}\right)$;

/* The Einsum is the Einstein summation convention,
$\left(\text{'bxn,btn}\to\text{ bxt'}\right)$ means the index change in the
Einstein summation convention. */

Algorithm 1 The algorithm of Causality POD-DeepONet

Loss function Given batch size $\mathcal{N}$ for training process and the
total number of the test records $N$, the loss function is defined as

$\mathcal{L}oss\left(\theta\right)=\frac{1}{\mathcal{N}}\sum_{i=1}^{\mathcal{N}}\sum_{j=1}^{N_{t}}\sum_{k=1}^{N_{{x}}}\left(\mathcal{G}\left(\overrightarrow{f}\right)\left({x}_{k},t_{j}\right)-y_{ijk}\right)^{2},$
(15)

where $N_{t}$ is the number of time step and $N_{x}$ is the number of points
on ${x}$ direction.

To evaluate the training process, the mean of the relative L2 error is
considered

$\mathcal{L}oss_{2}\left(\theta\right)=\frac{1}{\mathcal{N}}\sum_{i=1}^{\mathcal{N}}\frac{\left\|\mathcal{G}\left(\overrightarrow{f_{i}}\right)\left({x},t\right)-u_{i}\left({x},t\right)\right\|_{2}}{\left\|u_{i}\left({x},t\right)\right\|_{2}}.$
(16)

The relative L2 error in a complete epoch is defined as

$\mathcal{L}^{\mathcal{R}}_{\text{train}}=\frac{1}{B}\sum_{k=1}^{B}\frac{1}{\mathcal{N}}\sum_{i=1}^{\mathcal{N}}\frac{\left\|\mathcal{G}\left(\overrightarrow{f_{i}}\right)\left({x},t,\theta^{\left(k\right)}\right)-u_{ik}\left({x},t\right)\right\|_{2}}{\left\|u_{i}\left({x},t\right)\right\|_{2}},$
(17)

where $B$ is the number of batches, $\theta^{\left(k\right)}$ means the
parameters of neural network at $k$-th batch.

Similarly, we define the relative L2 error for the testing dataset

$\mathcal{L}_{\text{test}}^{\mathcal{R}}=\frac{1}{{N}}\sum_{i=1}^{N}\frac{\left\|\mathcal{G}\left(\overrightarrow{f_{i}}\right)\left({x},t\right)-u_{i}\left({x},t\right)\right\|_{2}}{\left\|u_{i}\left({x},t\right)\right\|_{2}}.$
(18)

Note the $N$ is the total number of test cases.

$\bullet$ Numerical Performance of Causality DeepONet with POD

Case 1: Constant wave-speed case In this test, we assume the wave speed is
$2$, and the right hand side $f\left(x\right)$ satisfies

$f\left(x,t\right)=\sum_{i=0}^{2}c_{i}t^{i}+\sum_{i=1}^{10}2(c^{2}-1)i^{2}\pi^{2}\left(a_{i}\cos(2i\pi
t)\cos\left(2i\pi x\right)+b_{i}\sin(2i\pi t)\sin(2i\pi x)\right),$ (19)

where $a_{i},b_{i},c_{i}$ are randomly sampled from $\left[0,1\right]$. The
exact solutions are from integrating the formula based on Duhamel’s method.
There are 21 basis based on 19, thus, the input size of the causality DeepONet
should be $21\times N_{t}$, where $N_{t}$ is the number of time steps. We set
$N_{t}=400$ for all other cases in this paper. We build the branch part and
trunk part of the causality DeepONet as $4$-layer neural network with $128$
hidden neurons each layer and $4$-layer neural network with $100$ hidden
neurons each layer respectively, both of whose outputs dimension is $500$ and
shrinking to dimension $N_{x}\times N_{t}$ by inner product with the
$x$-basis. The learning rates during training including the following cases
are $10^{-4}$. The batch size we choose for all the testing case is 20. The
results trained 500 epochs with 1000 cases are shown in Figure 1.

Figure 1: Random testing case after 500 epochs training for the simple case.
The maximum relative error is $5\%$.

Case 2: Variable wave-speed case

In this case we consider the case the wave speed is a function of $x,t$

$c\left(x,t\right)=\sqrt{\cos\left(2n\pi\left(t+x\right)\right)+1},$ (20)

and the exact solution is given by

$u\left(x,t\right)=\sum_{m=1}^{M}c_{m}\cos\left(2m\pi\left(t+x\right),\right)$
(21)

thus, the corresponding right hand side for (3) will be

$\displaystyle f\left(x,t\right)=\sum_{m=1}^{M}c_{m}\left(2m\pi\right)^{2}$
$\displaystyle\left[\cos\left(2\pi\left(m+n\right)t\right)\cos\left(2\pi\left(m+n\right)x\right)+\right.$
(22)
$\displaystyle-\sin\left(2\pi\left(m+n\right)t\right)\sin\left(2\pi\left(m+n\right)x\right)$
$\displaystyle+\cos\left(2\pi\left(n-m\right)t\right)\cos\left(2\pi\left(n-m\right)x\right)$
$\displaystyle\left.-\sin\left(2\pi\left(n-m\right)t\right)\sin\left(2\pi\left(n-m\right)x\right)\right].$

The inputs for the CPOD-DeepONet are the coefficients of basis funciton of $x$
in $\left(\ref{rhs_1dv}\right)$, since there are $4M$ basis, the input size of
the causality DeepONet should be $4M\times N_{t}$. The causality DeepONet is
as the same shape as in Case 1 but with different input dimensions. In this
case we assume $M=10$ and $n=10$ with 400 equal spaced time steps, and there
are 400 random sampled points on x-direction for each time step. We have
similar results as shown in Figure 2. It should be mentioned that the trained
CPOD-DeepONet can in fact be used to predict solution outside the trained
spatial domains as shown in Fig. 3.

Figure 2: Random testing case after 500 epochs training for the variable wave
speed case. The maximum relative error is $4.2\pm 2.3\%$.

Figure 3: Random testing case after 500 epochs training for the variable wave
speed case with random large domain on $x\in\left[-10,10\right]$. The maximum
relative error is $4.0\pm 1.4\%$.

##  3 A Recursive DeepPropNet for Learning Evolution PDE Operators

Following the semi-group formulation of evolution PDEs, the DeepPropNet
computes block by block along the time direction recursively. The initial
conditions for the PDEs will be included together with the force functions for
the time block as an input for DeepPropNet to assure the well-posedness of the
problem. And, the DeepPropNet is then constructed recursively with a Causality
DeepONet with POD to learn the PDE evolution operator over long time interval.

To simplify notation, we denote the initial conditions on the solution $u$ and
its velocity $v$ as

$u_{d}\left({x}\right)=\left[u_{0}\left({x}\right),v_{0}({x}),u_{1}\left({x}\right),v_{1}({x}),\ldots,u_{n}\left({x}\right),v_{n}({x})\right],$
(23)

and the DeepPropNet $P_{\theta}$. Following the notations of Causality
DeepONet with POD in equation $\left(\ref{cpod}\right)$, we define the
DeepPropNet as

$\mathcal{P}_{\theta}(\overrightarrow{f})(x,t)\sim\left\\{\overrightarrow{\sigma}_{br,v}\left(\overrightarrow{u}_{d}\right)\odot\overrightarrow{\sigma}_{trk,v}\left(t\right)+\overrightarrow{\sigma}_{br,c}\left(\overrightarrow{f}\right)\odot\overrightarrow{\sigma}_{trk,c}\left(t\right)\right\\}\cdot\overrightarrow{\sigma}_{basis}\left(x\right),$
(24)

where $\overrightarrow{\sigma}_{br,c}\left(\overrightarrow{f}\right)$ is the
branch net $\left(\ref{branch_cpod}\right)$,
$\overrightarrow{\sigma}_{trk,c}\left(t\right)$ is the corresponding trunk net
as shown in equation $\left(\ref{trunk_cpod}\right)$. Likewise,
$\overrightarrow{\sigma}_{br,v}\left(\overrightarrow{u}_{d}\right)\text{ and
}\overrightarrow{\sigma}_{trk,v}\left(t\right)$ are the trunk net and branch
net of the vanilla DeepONet,

$\displaystyle\overrightarrow{\sigma}_{br,v,i}\left(\overrightarrow{u}_{d}\right)=\sigma_{b}$
$\displaystyle\left(\sum_{n=0}^{N}\sum_{j=1}^{N}\xi_{i,j}^{n}u_{n}\left(s_{j}\right)+b_{i}\right),$
(25)

the corresponding $i$-th output of trunk net is

$\overrightarrow{\sigma}_{trk,v,i}\left(t\right)=\sigma_{trk}\left(\omega_{i}\cdot
t+\zeta_{i}\right).$ (26)

The notation is abused to some extent in above equations, but the
$\sigma_{trk}$’s and $\sigma_{b}$’s on the right hand side means different
neural networks. The definition of $\overrightarrow{\sigma}_{basis}$ is as the
same as $\left(\ref{basis}\right)$. Note the summation of
$\overrightarrow{\sigma}_{br,v}\left(\overrightarrow{u}_{d}\right)\odot\overrightarrow{\sigma}_{trk,v}\left(t\right)$
and
$\overrightarrow{\sigma}_{br,c}\left(\overrightarrow{f}\right)\odot\overrightarrow{\sigma}_{trk,c}\left(t\right)$
in the equation $\left(\ref{DP_CPOD}\right)$ is inspired from the
superposition of linear waves from different sources and other format could be
explored for nonlinear problems.

$\bullet$ Recursive Formulation of DeepPropNet: To predict the wave field in
the time block $\left[t_{0},t_{1}\right]$ at beginning, the Deep Propagator
DeepPropNet is to learn the mapping

$P_{\theta}:\left[u_{d}\left({x}_{j}\right),f\left({x}_{\xi},s\right),t_{0}\leq
s\leq t_{1}\right]\mapsto u\left({x},t\right)\quad t_{0}\leq t\leq
t_{1},{x}\in\mathbb{R}^{n}.$ (27)

Once the Deep Propagator is learned, the wave field in the next time block
$\left[t_{1},t_{2}\right]$ could be predicted by the Deep Propagator by using
the initial propagator as follows

$u\left(x,t\right)=P_{\theta}\left[P_{\theta}\left({x}_{j},t_{1}\right),{\dot{P}}_{\theta}\left({x}_{j},t_{1}\right),f\left({x}_{\xi},s\right),t_{1}\leq
s\leq t_{2}\right]\mapsto u\left({x},t\right)\quad t_{1}\leq t\leq
t_{2},{x}\in\mathbb{R}^{n},$ (28)

where the initial condition is replaced by the prediction of DeepPropNet at
time $t_{1}$.

The resulting propagator for the time block $[t_{1},t_{2}]$ then will be again
to be used to provide the initial condition for $t=t_{2}$ and the same initial
propagator network for the time period $[t_{0},t_{1}]$. This procedure can be
carried on recursively until the whole time period $[0,T]$ is covered,a global
propagator network is thus obtained. Since it is like the initial propagator
solver tracks the waves and propagates with the solutions along time direction
and the propagator itself is a deep neural network, we call it Deep
Propagator.

The schematics of the Deep Propagator DeepPropNet are shown in Algorithm 2.
This is another crux to solve high dimensional problem. As illustrated in
section 1, the input size could explode since the global dependence of source
term, by solving it block by block recursively combining with the memory-
efficient trick from separation of variables, learning high dimensional
evolutioin operators is manageable now.

Loss function Similar to the loss function in section 2, we will define the
loss function of DeepPropNet block by block in time. Given a batch size
$\mathcal{N}$ for training process, number of blocks $N_{b}$ and the total
number of the test records $N$, the loss function is defined as

$\mathcal{L}oss\left(\theta\right)=\frac{1}{\mathcal{N}}\sum_{i=1}^{\mathcal{N}}\sum_{n=1}^{N_{b}}\sum_{j=N_{t,n}}^{N_{t,n+1}}\sum_{k=1}^{N_{{x}}}\left(\mathcal{G}\left(\overrightarrow{f}\right)\left({x}_{k},t_{j}\right)-y_{ijk}\right)^{2},$
(29)

where the $n$-th block starts from $t_{N_{t,n}}$ but ends in $t_{N_{t,n+1}}$
and $N_{x}$ is the number of points on ${x}$ direction.

Likewise, the mean of the relative L2 error is considered as the relatvie L2
error in the total time-spatial domain

$\mathcal{L}oss_{2}\left(\theta\right)=\frac{1}{\mathcal{N}}\sum_{i=1}^{\mathcal{N}}\frac{\sum_{i=1}^{N_{b}}\left\|\mathcal{G}\left(\overrightarrow{f_{i}}\right)\left({x},t\right)-u_{i}\left({x},t\right)\right\|_{2,\Omega_{i}}}{\sum_{i=1}^{N_{b}}a_{i}z^{i}\left\|u_{i}\left({x},t\right)\right\|_{2,\Omega_{i}}},$
(30)

where $\Omega_{i}$ is the i-th block of time and spatial domain.

Input : $\left\\{\begin{aligned} &t\text{: Time step}\\\ &{x}\text{: Spatial
Coordinates}\\\ &u_{d}\left({x}\right)\text{: Initial conditions}\\\
&\left[a_{0}\left(s\right),a_{1}\left(s\right),b_{1}\left(s\right),\ldots,a_{N}\left(s\right),b_{N}\left(s\right)\right]\text{:
Coefficients of right hand side}\\\ &N\text{: Number of blocks at time
domain}.\end{aligned}\right.$

Output : $u\left({x},t\right)$

$u_{0}\left({x}\right)\leftarrow u_{d}\left({x}\right)$ for (  _$i=0;i
<N;i++$_ ) $\\{$

$u\left({x},t\right)\leftarrow\text{DeepPropaNet}\left(t,u_{0}\left({x}_{j}\right),\left[a_{0}\left(s\right),a_{1}\left(s\right),b_{1}\left(s\right),\ldots,a_{N}\left(s\right),b_{N}\left(s\right)\right]\right)$;

$u_{0}\left({x}\right)\leftarrow u\left({x},t_{i+1}\right)$;

$\\}$

Algorithm 2 The schematics of Deep Propagator

Input : $\left\\{\begin{aligned} &t\text{: Time step}\\\ &{x}\text{: Spatial
Coordinates}\\\ &u_{d}\left({x}\right)\text{: Initial conditions}\\\
&\left[a_{0}\left(s\right),a_{1}\left(s\right),b_{1}\left(s\right),\ldots,a_{N}\left(s\right),b_{N}\left(s\right)\right]\text{:
Coefficients of right hand side}\end{aligned}\right.$

Output : $u\left({x},t\right)$

$y_{t}\leftarrow\text{CDeepONet}\left(t,\left[a_{0}\left(s\right),a_{1}\left(s\right),b_{1}\left(s\right),\ldots,a_{N}\left(s\right),b_{N}\left(s\right)\right]\right)$
;

/* The output $y_{t}$ has size of
$[\text{batchsize}\times\text{Nt}\times\text{NBx}]$ */

$y_{i}\leftarrow\text{DeepONet}\left(t,u_{d}\left({x}_{j}\right)\right)$ ;

/* DeepONet is the vanilla DeepONet. The output $y_{i}$ has size of
$[\text{batchsize}\times\text{Nt}\times\text{NBx}]$ */

$y_{x}\leftarrow\left[1,\cos\left(2\pi{x}\right),\sin\left(2\pi{x}\right),\ldots\cos\left(2M\pi{x}\right),\sin\left(2M{x}\right)\right]$;

/* The output $y_{x}$ has size of
$[\text{batchsize}\times\text{Nx}\times\text{NBx}]$ */

$u\left(x,t\right)\leftarrow\text{Einsum}\left(\text{'bxn,btn}\to\text{
bxt'},y_{x},y_{t}+y_{i}\right)$

Algorithm 3 The algorithm of DeepPropaNet

##  4 Numerical Results of DeepPropNet

The equation and right hand sides are as the same as the non-constant wave
speed case $\left(\ref{wavespeed}\right)-\left(\ref{rhs_1dv}\right)$, but we
split 5 equal sized blocks in time domain and for each block, there are 80
sampled time locations for training. The causality DeepONet is as the same
shape as in Case 2 but with different input dimensions since the dependence of
time steps shrinks. The vanilla DeepONet is a composition of a $4$-layer fully
connected neural network with $128$ hidden neurons each layer as branch net
and a $4$-layer fully connected neural network with $100$ hidden neurons each
layer as trunk net, both of whose output dimensions are also $500$. The
initial data for each block will be given by true data. Figure 5 shows the
prediction of the exact solution up to time $t=1$ with a maximum relative
error of $3.4\pm 2.0\%$.

Training with DeepPropNet Prediction as Initial Condition at $t_{i}$. There
are two strategy to provide initial conditions for each block, either provide
the exact solution $u\left(x,t_{i}\right)$ or the prediction
$P_{\theta}\left[u_{0}\left(x_{j}\right),v_{0}\left(x_{j}\right),f\right]\left(x,t_{1}\right)$.
We compare the relative L2 error evolution of these $2$ different process and
as shown in Figure 6, they have similar loss convergence. The results are
shown in Figure 4 and Figure 5 respectively. This is a natural test about the
capability of the interpolation in training domain of the neural network, but
it offers a convenience that there is no further need to prepare exact data on
grids for initial conditions of each time blocks.

Figure 4: One of the random testing cases after 500 epochs training for the
variable wave speed case with DeepPropNet whose initial condition is given by
the prediction in $x\in\left[-0.5,1\right]$. The maximum relative error in the
whole domain is $3.2\pm 1.9\%$. (Left) Exact, (Middle) Prediction, (Right)
Relative Error.

Figure 5: One of the random testing cases after 500 epochs training for the
variable wave speed case with DeepPropNet whose initial condition is given by
true data in $x\in\left[-0.5,1\right]$. The maximum relative error in the
whole domain is $3.4\pm 2.0\%$. (Left) Exact, (Middle) Prediction, (Right)
Relative Error.

(a) The case using prediction
$P_{\theta}\left[u_{0}\left(x_{j}\right),v_{0}\left(x_{j}\right),f\right]\left(x,t_{1}\right)$
as initial condition for each time block.

(b) The case using exact solution $u\left({x},t_{i}\right)$ for initial
condition for each time block.

Figure 6: Relative L2 error evolution during training

##  5 Conclusion and Future works

In this paper, we proposed the recursive way to construct DeepPropNet - a DNN
propagator for evolution system over large time period by using a single
building block propagator over a small time period, thus reducing the overall
complexity and size of the neural network required. For the design of the
DeepPropNet, we also extended the Causality DeepONet with POD with specific
basis to alleviate the memory burden for large spatial variables. By
seperately handling the spacial and temporal domain, we gain not only the
memory efficiency but also the training boost. The preliminary numerical
results have shown the feasibility of this recursive DeepPropNet in predicting
the time evolution of wave propagations.

The proposed DeepPropNet here is based on a supervised learning approach where
the data can be generated by a separate numerical methods or observation data
or analytical solution when available. In theory, we could also use a
unsupervised learning procedure to train the DeepPropNet by using the
following residue of the PDEs as the loss function.

$\mathcal{L}\text{oss}=\int_{\Omega}\mid\mid\mathcal{L}_{t}P_{\theta}\left[u_{0},v_{0},f\right]-\mathcal{L}_{x}P_{\theta}\left[u_{0},v_{0},f\right]-f\mid\mid^{2}dxdt.$
(31)

It is also natural to extend this framework to learn high frequency problem,
since the spatial oscillating terms will be handled explicitly by the basis,
but the temporal oscillating terms could be handled by the Causality DeepONet,
based on our experience in the work [4].

Future work will be conducted on more complex evolution systems, the
unsupervised training, training over partial time domain, highly oscillating
problems and higher dimensional problems as well as initial boundary value
problems.

## References

  * [1] Brek Meuris, Saad Qadeer, and Panos Stinis.  Machine-learning custom-made basis functions for partial differential equations.  arXiv:2111.05307 [cs, math], November 2021.
  * [2] Lu L, Jin P, Pang G, Zhang Z, Karniadakis GE. Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators. Nature Machine Intelligence. 2021 Mar;3(3):218-29.
  * [3] Li Z, Kovachki NB, Azizzadenesheli K, Bhattacharya K, Stuart A, Anandkumar A. Fourier Neural Operator for Parametric Partial Differential Equations. International Conference on Learning Representations 2020 Sep 28.
  * [4] L.Z. Liu, K. Nath, W, Cai, Causality DeepONet for Nonlinear Operators for Building Seismic Wave Responses, manuscript, December, 2021.
  * [5] R. Shankar, Principles of Quantum Mechanics, Kluwer Academic/Plenum Publisher,1994.
  * [6] Kaushik Bhattacharya, Bamdad Hosseini, Nikola B. Kovachki, and Andrew M. Stuart.  Model Reduction And Neural Networks For Parametric PDEs.  The SMAI journal of computational mathematics, 7:121–157, July 2021.
  * [7] Lu Lu, Xuhui Meng, Shengze Cai, Zhiping Mao, Somdatta Goswami, Zhongqiang Zhang, and George Em Karniadakis.  A comprehensive and fair comparison of two neural operators (with practical extensions) based on FAIR data.  arXiv:2111.05512 [physics], November 2021.