capbtabboxtable[][] heightadjust=object

# Posterior Estimation for Dynamic PET imaging using Conditional Variational
Inference

Xiaofeng Liu  Thibault Marin  Amal Tiss  Jonghye Woo  Georges El Fakhri  and
Jinsong Ouyang X. Liu, T. Marin, A. Tiss, J. Woo, G. El Fakhri, and J. Ouyang∗
are with the Gordon Center for Medical Imaging, Massachusetts General Hospital
and Harvard Medical School, Boston, MA 02114 USA.(∗corresponding author:
<EMAIL_ADDRESS>
###### Abstract

This work aims efficiently estimating the posterior distribution of kinetic
parameters for dynamic positron emission tomography (PET) imaging given a
measurement of time of activity curve. Considering the inherent information
loss from parametric imaging to measurement space with the forward kinetic
model, the inverse mapping is ambiguous. The conventional (but expensive)
solution can be the Markov Chain Monte Carlo (MCMC) sampling, which is known
to produce unbiased asymptotical estimation. We propose a deep-learning-based
framework for efficient posterior estimation. Specifically, we counteract the
information loss in the forward process by introducing latent variables. Then,
we use a conditional variational autoencoder (CVAE) and optimize its evidence
lower bound. The well-trained decoder is able to infer the posterior with a
given measurement and the sampled latent variables following a simple
multivariate Gaussian distribution. We validate our CVAE-based method using
unbiased MCMC as the reference for low-dimensional data (a single brain
region) with the simplified reference tissue model.

##  I Introduction

For many physical imaging systems, a common problem is to determine hidden
system parameters $\bm{x}$ from a measurement $\bm{y}$. For example, in
kinetic modeling of tau PET studies, our goal can be inferring the posterior
distribution of $p(\bm{x}|\bm{y})$ w.r.t. distribution volume ratio (DVR) for
a given measurement of time-activity curve (TAC) in a target brain region.

Usually, the physical model describing how measurable quantities $\bm{y}$
arise from the hidden parameters $\bm{x}$ is well-defined by a forward
process. For example (used in this work), the forward process can be,
$\bm{y}=f(\bm{x})+\varepsilon$, where $\bm{x}=\\{DVR,k_{2},R_{1}\\}$ in a
target brain region ($k_{2}$ is the rate constant from free to plasma
compartment, $R_{1}$ is the ratio of rate constants for transform from plasma
to free compartment), $\bm{y}$ is the measured TAC in the region, $f$
represents the kinetic modeling using simplified reference tissue model (SRTM)
[5] in PET, and noise $\varepsilon$ follows a Gaussian distribution. As a
result, the likelihood, $p(\bm{y}|\bm{x})$, is defined. If we assume a prior,
$p(\bm{x})$, based on our knowledge before the measurement. The posterior
distribution is determined as $p(\bm{x}|\bm{y})\propto
p(\bm{y}|\bm{x})p(\bm{(}x))$. The conventional approach to sample the
posterior distribution is to follow a rejection sampling scheme with Markov
Chain Monte Carlo (MCMC) [2]. MCMC is known to produce unbiased asymptotical
estimation [2], while it is impractical for high-dimensional data due to its
high computational cost. In this work, we propose a deep-learning-based
approach to sample posterior distribution.

##  II Methodology

In this work, we resort to the deep variational Bayes for the estimation of
the posterior distribution. Specifically, we propose to use a conditional
variational autoencoder (CVAE) framework and optimize its evidence lower bound
(ELBO) for accurate posterior estimation. To counteract the inherent
information loss of the forward process, we introduce latent variables
$\bm{z}$, which capture the loss of information in the forward process. Thus,
our approach explicitly learns to associate hidden parameters $\bm{x}$ with
the unique pairs [$\bm{y}$, $\bm{z}$] of measurements and latent variables.

Figure 1: The posterior estimation framework using CVAE. The network is
optimized with loss function $\mathcal{L}=\mathcal{L}_{MSE}+\mathcal{L}_{KL}$.

With the paired data in $\\{\bm{x}_{i},\bm{y}_{i}\\}_{i=1}^{M}$, we
reconstruct $\bm{x}$ to ensure $\bm{y}$ and $\bm{z}$ have complementary
information to recover $\bm{x}$. Additionally, we make sure that the density
$p(\bm{z})$ of the latent variables is shaped as a multivariate Gaussian
distribution. Thus, the CVAE represents the desired posterior
$p(\bm{x}|\bm{y})$ by a deterministic trained decoder, i.e., $\bm{x}=$
Decoder$(\bm{y},\bm{z})$, that transforms the known distribution $p(\bm{z})$
to $\bm{x}$-space, conditional on $\bm{y}$. The general idea is illustrated in
Fig. 1.

Based on the well-defined forward process, our CVAE framework aims to
efficiently infer the posterior $p(\bm{x}|\bm{y})$ for a given measurement
$\bm{y}$ and prior $p(\bm{x})$. The encoder
$[\bm{x},\bm{y}]\rightarrow{\bm{z}}$ encodes the latent variation $\bm{z}$
with $\bm{x}$ and $\bm{y}$. We would expect $\bm{z}$ to capture the
information lost in the forward process. This is implicitly achieved by
enforcing a information bottleneck [1] of the small size $\bm{z}$ to have
sufficient content to recover $\bm{x}$. The decoder
$[\bm{z},\bm{y}]\rightarrow\tilde{\bm{x}}$ reconstructs the parameters
$\bm{x}$ with $\bm{z}$ and $\bm{y}$. We empirically use three fully connected
layers for both encoder and decoder. Similar to original VAE [4], the encoder
has two output vectors, $i.e.,$ $\mu$ and $\sigma$. We then utilize the
reparametric trick $\bm{z}=\mu+\sigma\odot\epsilon$, where
$\epsilon\in\mathcal{N}(0,I)$. The posterior distribution of $\bm{z}$ is
$q(\bm{z}|\bm{x},\bm{y})\sim\mathcal{N}(\bm{z};\mu,\sigma^{2})$.

Figure 2: Comparison of the MCMC and our deep learning (DL) method. We show
the smoothed curves of 25 bins histograms. We set $DVR=1.25$, $k_{2}=0.03$ and
$R_{1}=0.51$ to calculate $\bm{y}^{*}$ with SRTM and sample 90,000 $\bm{z}$
following a Gaussian prior.  Figure 3: Detailed CVAE-based posterior
estimation framework for a single voxel. Only decoder is used for testing.

In training, the ELBO [4] of our CVAE includes the Kullback-Leibler (KL)
divergence and the reconstruction loss. In practice, the KL-divergence can be
computed as

$\displaystyle\mathcal{L}_{KL}(\bm{z};\mu,\sigma)=\frac{1}{2}\sum^{M_{z}}_{j=1}(1+{\rm
log}(\sigma_{j}^{2})-\mu_{j}^{2}-\sigma_{j}^{2}),\vspace{-10pt}$ (1)

where $M_{z}=10$ is the dimension of the latent code $\bm{z}$. For the
reconstruction error, we simply adopt the pixel-wise mean square error (MSE).
Let $\tilde{\bm{x}}$ be the reconstructed $\bm{x}$, their MSE loss can be
formulated as:

$\displaystyle\mathcal{L}_{MSE}(\bm{x},\tilde{\bm{x}})=\frac{1}{2}\sum^{M_{x}}_{j=1}||\bm{x}_{j}-\tilde{\bm{x}}_{j}||^{2}_{2},\vspace{-10pt}$
(2)

where ${M_{x}}=3$ indicates the dimension of $\bm{x}$ or $\tilde{\bm{x}}$.

In testing, only the decoder is used. Given an observation $\bm{y}^{*}$, we
sample $\bm{z}\sim\mathcal{N}(0,I)$ and concatenate each of them with
$\bm{y}^{*}$. If ${\bm{x}}\sim p(\bm{x}|\bm{y}^{*})$ and $\bm{y}=\bm{y}^{*}$,
we have following two claims: 1) $\tilde{\bm{x}}\sim p(\bm{x}|\bm{y}^{*})$
because $\mathcal{L}_{MSE}=0$; 2) the corresponding conditional probability
$q(\bm{z}|\bm{x},\bm{y}^{*})$ is equivalent to $\mathcal{N}(0,I)$ because
$\mathcal{L}_{KL}=0$. Since the distribution of $\tilde{\bm{x}}$ is solely
determined by $\bm{z}$ distribution and $\bm{y}$ input, the above two claims
guarantee that the decoder with inputs $\bm{z}\sim\mathcal{N}(0,I)$ and
$\bm{y}=\bm{y}^{*}$ outputs posterior distribution $p(\bm{x}|\bm{y}^{*})$.

##  III Experiments and results

The approach described above can be applied to many medial imaging problems.
In this work, we estimate the posterior distribution of DVR in a single brain
region given a measurement of TAC in the region using both our CVAE-based
approach and MCMC. The unbiased distribution estimated by MCMC is the
reference used to validate our approach. The forward process is defined by a
kinetic model and a noise model. The kinetic model is SRTM written as:

$\displaystyle\vspace{-10pt}\frac{dC_{T}(t)}{dt}=R_{1}\frac{dC_{R}(t)}{dt}+k_{2}C_{R}(t)-\frac{k_{2}}{DVR}C_{T}(t),\vspace{-10pt}$
(3)

where $C_{T}(t)$ and $C_{R}(t)$ are the activity concentrations in the target
and reference regions, respectively. The noise model is defined as:
$y_{n}=\int_{t_{n-1}}^{t_{n}}C_{T}(t)\mathrm{d}t+\epsilon_{n}$, where
$\frac{\epsilon_{n}}{\sigma\sqrt{\Delta t_{n}/T}}\sim\mathcal{N}(0,1)$ and
$T=\sum_{n=1}^{N}\Delta t_{n}$, $\Delta t_{n}={t_{n}}-{t_{n-1}}$, $\sigma$ is
the standard deviation.

With enough samples, the MCMC algorithm would produce the unbiased posterior
samples asymptotically for Bayesian inference [3]. Without loss of generality,
we simply set the number of time frame $N=40$, $\Delta t_{n}=0.5$ min in the
first ten frames and $\Delta t_{n}=1$ min in the later 30 frames, and
$\sigma\sim$ Gamma(1,1) in our forward model. We generate $M=10,000$ samples
for training. Specifically, we sample $DVR\sim\mathcal{N}(1,1)$,
$k_{2}\sim\mathcal{N}(0.05,0.01)$, $R_{1}\sim\mathcal{N}(0.5,0.1)$, and only
select the positive values.

Metropolis-Hastings sampler (MHS) [3], a variant of MCMC, was used to sample
from the ground truth posterior. In testing, we performed 100,000 iterations
of random walk Metropolis sampling for a generated measurement using SRTM with
10,000 burn-in steps. With the SRTM in a single voxel, PyMC111https://pymc-
devs.github.io/pymc/ takes about 10 mins to infer 90,000 samples. In testing,
our learning-based CVAE can infer 90,000 samples of a given $\bm{y}^{*}$ and
the sampled $\bm{z}$ within 30s. We note that MCMC takes a similar process
time for every $\bm{y}^{*}$. In contrast, the trained CVAE can make an
efficient posterior estimation for different $\bm{y}^{*}$.

The detailed network implementation using our approach is shown in Fig. 3. We
construct the paired dataset $\\{\bm{x}_{i},\bm{y}_{i}\\}_{i=1}^{M}$ by
randomly sampling $\bm{x}_{i}$ and $\sigma$ following a Gaussian and Gamma
distribution, respectively. Then, we generate its corresponding $\bm{y}_{i}$
with the forward kinetic model.

The comparison of the posterior estimation results is shown in Fig. 2. The
posterior distributions using MCMC and our proposed approach are in good
agreement.

##  IV Conclusions

We have proposed a CVAE framework for efficient posterior estimation of
dynamic PET imaging. The preliminary results for sampling DVR posterior
distribution show good agreement between unbiased MCMC and our proposed
approach. Our approach has much faster inference than MCMC. In the future, we
will apply it to high-dimensional imaging data.

## References

  * [1] A. A. Alemi, I. Fischer, J. V. Dillon, and K. Murphy.  Deep variational information bottleneck.  arXiv preprint arXiv:1612.00410, 2016.
  * [2] C. Andrieu, N. De Freitas, A. Doucet, and M. I. Jordan.  An introduction to mcmc for machine learning.  Machine learning, 50(1):5–43, 2003.
  * [3] S. Chib and E. Greenberg.  Understanding the metropolis-hastings algorithm.  The american statistician, 49(4):327–335, 1995.
  * [4] D. P. Kingma and M. Welling.  Auto-encoding variational bayes.  arXiv preprint arXiv:1312.6114, 2013.
  * [5] A. A. Lammertsma and S. P. Hume.  Simplified reference tissue model for pet receptor studies.  Neuroimage, 4(3):153–158, 1996.


# Accelerating Simulation of Stiff Nonlinear Systems using Continuous-Time
Echo State Networks

Ranjan Anantharaman
Massachusetts Institute of Technology
<EMAIL_ADDRESS>
&Yingbo Ma
Julia Computing Inc.
Shashi Gowda
Massachusetts Institute of Technology
&Chris Laughman
Mitsubishi Electric Research Laboratories
&Viral Shah
Julia Computing Inc.
&Alan Edelman
Massachusetts Institute of Technology
&Chris Rackauckas
Massachusetts Institute of Technology,
University of Maryland, Baltimore

###### Abstract

Modern design, control, and optimization often requires simulation of highly
nonlinear models, leading to prohibitive computational costs. These costs can
be amortized by evaluating a cheap surrogate of the full model. Here we
present a general data-driven method, the continuous-time echo state network
(CTESN), for generating surrogates of nonlinear ordinary differential
equations with dynamics at widely separated timescales. We empirically
demonstrate near-constant time performance using our CTESNs on a physically
motivated scalable model of a heating system whose full execution time
increases exponentially, while maintaining relative error of within 0.2 %. We
also show that our model captures fast transients as well as slow dynamics
effectively, while other techniques such as physics informed neural networks
have difficulties trying to train and predict the highly nonlinear behavior of
these models.

##  1 Introduction

Stiff nonlinear systems of ordinary differential equations, i.e., systems that
are ill-conditioned in a computational sense, are widely prevalent throughout
science and engineering [30, 26] and are characterized by dynamics with widely
separated time scales. These systems require highly stable numerical methods
to use large step-sizes reliably [8], and also tend to be computationally
expensive to solve. Even with state-of-the-art simulation techniques, design,
control, and optimisation of these systems remains intractable in many
realistic engineering applications [1]. To address these challenges, we
develop a data-driven technique, the continuous-time echo state network
(CTESN), to obtain an approximation to the system (called a "surrogate") whose
forward simulation time is relatively inexpensive while maintaining reasonable
accuracy on these highly ill-conditioned models.

A popular class of methods to produce surrogates is projection based model
order reduction, such as the proper orthogonal decomposition (POD) [1]. This
method computes "snapshots" of the trajectory and uses the singular value
decomposition of the linearization in order to construct a basis of a subspace
of the snapshot space, and the model is remade with a change of basis.
However, if the system is very nonlinear, the computational complexity of this
linearization-based reduced model can be almost as high as for the original
model. One way to overcome this difficulty is through empirical interpolation
methods [18]. Another recent method seeks to perform variable transformations
to change the dynamics of the system before performing model order reduction
on the transformed system [21]. However, such transformations could be system
specific and it is not known whether such transformations exist for every
system. Other methods to produce surrogates generally utilize the structural
phonomena known about highly regular systems like partial differential
equation discretizations [7].

Modern large-scale engineering models generated by symbolic systems such as
Modelica [17, 31] and
ModelingToolkit111https://github.com/SciML/ModelingToolkit.jl/ can be quite
heterogeneous and nonlinear without much repeated structure. Without the
structure inherent in problems like partial differential equations, these
problem-specific automated variable transformations and changes of bases could
be hard to implement.

Numerical stiffness [28] and multiscale dynamics represent an additional
challenge. While slow processes are relatively easy to capture, fast
transients are not. Some model-order reduction techniques simply eliminate the
fast modes and capture the slow modes accurately, such as the computational
singular perturbation (CSP) method [10], which results in a non-stiff system.
However, many systems (e.g., chemical reaction systems) require transient
activations to properly arrive at the overarching dynamics. In this case,
accurate representation of the fast dynamics can be required for appropriately
reconstructing the constituent behavior [11]. We will demonstrate how a new
surrogate technique, the CTESN, can overcome these difficulties and accurately
capture highly nonlinear heterogeneous stiff systems in this paper.

Figure 1: Training performance on the the Robertson equations. We see that our
network is able to reproduce the dynamics, capturing the slow (second row) and
fast transients (third row) with low error

##  2 Continuous-Time Echo State Networks

Figure 2: TEST performance on scalable heating system with 10 rooms. When
tested with parameters it has not seen in training, our surrogate is able to
reproduce the behaviour of the system to within 0.01% error.

Echo State Networks (ESNs) are a framework for computation derived from
recurrent neural network theory that maps signals into higher dimensional
computational spaces through the dynamics of a fixed, non-linear system called
a "reservoir". They are usually defined as high-dimensional discrete dynamical
systems with an output readout layer which is then fit to an input signal
[19]. While ESNs have traditionally been used as a stand-in for recurrent
neural networks [14, 16], recent work has applied ESNs to learning chaotic
systems [3, 4], nonlinear systems identification [13], bio-signal processing
[15], and robot control [20],

Here we describe a new variant of ESNs which we call continuous-time echo
state networks (CTESNs). Let $N$ be the dimension of our model, and let $P$ be
a Cartesian space of parameters under which the model is expected to operate.
The CTESN of dimension $N_{R}$ is defined as

$\displaystyle r^{\prime}$ $\displaystyle=\text{tanh}(Ar+W_{in}x(p^{*},t))$
(1) $\displaystyle x(t)$ $\displaystyle=W_{out}r(t),$ (2)

where $A$ is a fixed sparse random graph of dimension $N_{R}\times N_{R}$ and
$W_{in}$ is a fixed random dense matrix of dimensions $N_{R}\times N$. Given
these fixed values, the readout matrix $W_{out}$ is the only trained portion
of this network and is obtained through a least squares fit against the
original timeseries.

Figure 3: Reliability of surrogate through parameter space. We sampled our
grid at over 500 grid points and plotted a heatmap of test error through our
parameter space. We find our surrogate performs reliably even at the border of
our space with error within 0.1%

To obtain a surrogate that predicts the dynamics at new physical parameters,
the reservoir projection $W_{out}$ is fit against many solutions at
$p_{1},\dots,p_{n}$, where $n$ is the number of training data points sampled
from the parameter space. Using these fits, an interpolating function
$W_{out}(p)$ between the matrices can be trained. For the examples in this
paper, we have sampled the high-dimensional spaces using a Sobol low-
discrepancy sequence [27] and interpolated the $W_{out}$ using a radial basis
function provided by Surrogates.jl 222https://github.com/SciML/Surrogates.jl.
A prediction $\hat{x}(t)$ for a parameter vector $\hat{p}$ that is not in the
training set would simply be:

$\hat{x}(t)=W_{out}(\hat{p})r(t)$ (3)

A strong advantage of our method is its ease of implementation and ease of
training. Global $L_{2}$ fitting via stabilized methods like SVD are robust to
ill-conditioning, alleviating many of the issues encountered when attempting
to build neural surrogates of such equations. Also note that in this
particular case, the readout matrix is fit against the same reservoir time
series. Prediction does not need to simulate the reservoir, providing an extra
acceleration.

##  3 Case Studies

###  3.1 Robertson Equations and High Stiffness

We first consider Robertson’s chemical reaction equations:

$\displaystyle\dot{y_{1}}$ $\displaystyle=-0.04y_{1}+10^{4}y_{2}\cdot y_{3}$
(4) $\displaystyle\dot{y_{2}}$ $\displaystyle=0.04y_{1}-10^{4}y_{2}\cdot
y_{3}-3\cdot 10^{7}y_{2}^{2}$ (5) $\displaystyle\dot{y_{3}}$
$\displaystyle=3\cdot 10^{7}y_{2}^{2}$ (6)

where $y_{1},y_{2},y_{3}$ are the concentration of three chemical reaction
species. This system has widely separated time constants ($0.04,10^{4},3\cdot
10^{7}$), is well known to be very stiff [9, 25, 24] and is used as an example
to evaluate integrators of stiff ODEs [12]. Finding an accurate surrogate for
this system is difficult because it needs to capture both the stable slow
reacting system and the fast transients. Additionally, the surrogate needs to
be consistent with this system’s implicit physical constraints, such as the
conservation of matter ($y_{1}+y_{2}+y_{3}=1$) and positivity of the variables
($y_{i}>0$), in order to provide a stable solution. Figure 1 shows that
training our CTESN yields good performance for both the larger and smaller
quantities while accurately capturing the fast transients and respecting
physical constraints. Thus, for any set of parameters, if we know the
approximate reading out matrix for the reservoir, we can achieve accurate
predictions with $<1\%$ error.

To further demonstrate the difficulty of this equation, we used the
NeuralPDE.jl 333https://github.com/SciML/NeuralPDE.jl physics-informed neural
network (PINN) training framework to attempt to solve these equations using
PINNs. PINNs are a recent data-driven method to obtain surrogates for partial
differential equations by directly training the neural network to be the
solution of the differential equation [23]. This software has been verified to
automatically accurately solve many of the equations expounded upon in [23],
including nonlinear partial equations like the Burger’s equation,
demonstrating its accuracy and the robustness of the implementation. However,
when applied to Equation 4, the PINN training infrastructure was unable to
adequately train neural networks to the solution the equation, never lowering
below a relative error of 100% due to difficulties in capturing the short
transient. The reason for the difficulty can be attributed to recently
identified results in gradient pathologies in the training arising from
stiffness [29]. With a highly ill-conditioned Hessian in the training process
due to the stiffness of the equation, it would take an intractably long time
to train a PINN surrogate on this problem with current methods. We note that
in general, stiff systems of this form may be hard to capture by neural
networks directly as they show a bias towards low frequency functions [22].

Figure 4: Scaling performance of surrogate on heating system. Our surrogate’s
runtime grows in near-constant time while the benchmark problem shows
exponential growth in runtime.

###  3.2 Stiffly-Aware Surrogates of HVAC Systems

Our second test problem is a simplified, lumped-parameter model of a heating
system with a central heater supplying heat to several rooms through a
distribution network. The resulting system of equations is very stiff and has
localized activity. This provides a useful test for our method, exhibiting
stiff dynamics, and is a useful scalable benchmark used in the engineering
community [2].

The size of the heating system is scaled by a parameter $N$ which refers to
the number of users/rooms. This produces a system with $2N+1$ equations. This
“scalability” lets us test how our method scales. To train the surrogate, we
define a parameter space $P$ under which we expect it to operate. First, we
assume set point temperature of each room to be between $17\degree C$ to
$23\degree C$, which is what we expect users to set their room temperature to.
Each room is warmed by a heat conducting fluid, whose set point is likely
between $65\degree C$ and $75\degree C$. Thus the parameter space over which
we expect our surrogate to work is the rectangular space denoted by
$[17\degree C,23\degree C]\times[65\degree C,75\degree C]$. Figure 2
demonstrates that the training technique is accurately able to find matrices
$W_{out}$ which capture the stiff system within $0.01\%$ error. Figure 3
demonstrates that the interpolated $W_{out}(p)$ is able to adequately capture
the dynamics throughout the trained parameter space. Lastly, Figure 4
demonstrates the almost $O(1)$ cost of the surrogate evaluation, where at the
high end the surrogate accelerates a 201 dimensional stiff ODE system by over
15x.

##  4 Discussion, Conclusion & Future Work

We present CTESNs, a data-driven method for generating surrogates of nonlinear
ordinary differential equations with dynamics at widely separated timescales.
Our method maintains accuracy for different parameters in a chosen parameter
space, and shows favourable scaling with system size.

Known limitations: Our method utilizes the continuous nature of differential
equation solutions. Hybrid dynamical systems, such as those with event
handling [5], can introduce discontinuities into the system which will not be
adequately captured by the current method. Further extensions to the method
will be required to handle both derivative discontinuities and events present
in Filippov dynamical systems [6].

Reliability: In this work we have not proven any error bounds for our method,
and while the empirical results are promising, some formal intuition on its
reliability is a possible area of further study. In addition, more testing
needs to be done on CTESNs under controls and as components within larger
systems.

Performance: While we show promising scaling results, we have not shown tight
bounds on execution time. For the surrogate on the heating system, we do not
need solve the reservoir ODE every time we predict a time series because the
surrogate has been trained for various parameters against the same reservoir.
In other words, the reservoir ODE itself has no dependence on the parameter
space (only the readout matrix does). Hence, we can simply return the solution
to the reservoir ODE after training and just use it directly for prediction.
However, this approach is memory-intensive: the reservoir time series of
dimensions and the dense readout matrices of dimensions $N\times N_{R}$ are
both large. In addition, for CTESN formulations that take in input forcing
functions or have some dependence on parameters, the reservoir would need to
be simulated every single time a prediction is made. Future work will comment
on other such formulations. We do note in passing that simulating the
reservoir is quite fast in practice as it is non-stiff, and thus techniques
which regenerate reservoirs on demand will likely not incur a major run time
performance cost.

Opportunities: Since CTESNs are largely data driven (the physics dependence is
largely driven by using a candidate solution as "input" to the reservoir), an
opportunity exists to train surrogates of black-box systems. Further
opportunities could explore utilizing more structure within equations for
building a more robust CTESN or decrease the necessary size of the reservoir.

## Broader Impact

This section does not apply to our work.

## Acknowledgement

The information, data, or work presented herein was funded in part by the
Advanced Research Projects Agency-Energy (ARPA-E), U.S. Department of Energy,
under Award Number DE-AR0001222. The views and opinions of authors expressed
herein do not necessarily state or reflect those of the United States
Government or any agency thereof.

## References

  * [1] Peter Benner, Serkan Gugercin and Karen Willcox  “A survey of projection-based model reduction methods for parametric dynamical systems”  In _SIAM review_ 57.4 SIAM, 2015, pp. 483–531
  * [2] Francesco Casella  “Simulation of large-scale models in Modelica: State of the art and future perspectives”  In _Proceedings of the 11th International Modelica Conference_ , 2015, pp. 459–468
  * [3] Ashesh Chattopadhyay, Pedram Hassanzadeh, Devika Subramanian and Krishna Palem  “Data-driven prediction of a multi-scale Lorenz 96 chaotic system using a hierarchy of deep learning methods: Reservoir computing, ANN, and RNN-LSTM”  EarthArXiv, 2019
  * [4] Nguyen Anh Khoa Doan, Wolfgang Polifke and Luca Magri  “Physics-informed echo state networks for chaotic systems forecasting”  In _International Conference on Computational Science_ , 2019, pp. 192–198  Springer
  * [5] D Ellison  “Efficient automatic integration of ordinary differential equations with discontinuities”  In _Mathematics and Computers in Simulation_ 23.1 Elsevier, 1981, pp. 12–20
  * [6] Aleksei Fedorovich Filippov  “Differential equations with discontinuous righthand sides: control systems”  Springer Science & Business Media, 2013
  * [7] Michalis Frangos, Youssef Marzouk, Karen Willcox and B Bloemen Waanders  “Surrogate and reduced-order modeling: a comparison of approaches for large-scale statistical inverse problems [Chapter 7]”  John Wiley & Sons, 2010
  * [8] C William Gear  “Numerical initial value problems in ordinary differential equations”  In _nivp_ , 1971
  * [9] Matthias K Gobbert  “Robertson’s example for stiff differential equations”  In _Arizona State University, Technical report_ Citeseer, 1996
  * [10] M Hadjinicolaou and Dimitris A Goussis  “Asymptotic solution of stiff PDEs with the CSP method: the reaction diffusion equation”  In _SIAM Journal on Scientific Computing_ 20.3 SIAM, 1998, pp. 781–810
  * [11] Adrien Henry and Olivier C Martin  “Short relaxation times but long transient times in both simple and complex reaction networks”  In _Journal of the Royal Society Interface_ 13.120 The Royal Society, 2016, pp. 20160388
  * [12] ME Hosea and LF Shampine  “Analysis and implementation of TR-BDF2”  In _Applied Numerical Mathematics_ 20.1-2 Elsevier, 1996, pp. 21–37
  * [13] Herbert Jaeger  “Adaptive nonlinear system identification with echo state networks”  In _Advances in neural information processing systems_ , 2003, pp. 609–616
  * [14] Herbert Jaeger, Mantas Lukoševičius, Dan Popovici and Udo Siewert  “Optimization and applications of echo state networks with leaky-integrator neurons”  In _Neural networks_ 20.3 Elsevier, 2007, pp. 335–352
  * [15] Dhireesha Kudithipudi et al.  “Design and analysis of a neuromemristive reservoir computing architecture for biosignal processing”  In _Frontiers in neuroscience_ 9 Frontiers, 2016, pp. 502
  * [16] Mantas Lukoševičius and Herbert Jaeger  “Reservoir computing approaches to recurrent neural network training”  In _Computer Science Review_ 3.3 Elsevier, 2009, pp. 127–149
  * [17] Sven Erik Mattsson, Hilding Elmqvist and Martin Otter  “Physical system modeling with Modelica”  In _Control Engineering Practice_ 6.4 Elsevier, 1998, pp. 501–510
  * [18] VB Nguyen, M Buffoni, K Willcox and BC Khoo  “Model reduction for reacting flow applications”  In _International Journal of Computational Fluid Dynamics_ 28.3-4 Taylor & Francis, 2014, pp. 91–105
  * [19] Mustafa C Ozturk, Dongming Xu and José C Prıncipe  “Analysis and design of echo state networks”  In _Neural computation_ 19.1 MIT Press, 2007, pp. 111–138
  * [20] Athanasios Polydoros, Lazaros Nalpantidis and Volker Krüger  “Advantages and limitations of reservoir computing on model learning for robot control”  In _IROS Workshop on Machine Learning in Planning and Control of Robot Motion, Hamburg, Germany_ , 2015
  * [21] Elizabeth Qian, Boris Kramer, Benjamin Peherstorfer and Karen Willcox  “Lift & Learn: Physics-informed machine learning for large-scale nonlinear dynamical systems”  In _Physica D: Nonlinear Phenomena_ 406 Elsevier, 2020, pp. 132401
  * [22] Nasim Rahaman et al.  “On the spectral bias of neural networks”  In _International Conference on Machine Learning_ , 2019, pp. 5301–5310  PMLR
  * [23] Maziar Raissi, Paris Perdikaris and George E Karniadakis  “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations”  In _Journal of Computational Physics_ 378 Elsevier, 2019, pp. 686–707
  * [24] HH Robertson  “Numerical integration of systems of stiff ordinary differential equations with special structure”  In _IMA Journal of Applied Mathematics_ 18.2 Oxford University Press, 1976, pp. 249–263
  * [25] HH Robertson and J Williams  “Some properties of algorithms for stiff differential equations”  In _IMA Journal of Applied Mathematics_ 16.1 Oxford University Press, 1975, pp. 23–34
  * [26] Lawrence F Shampine and Charles William Gear  “A user’s view of solving stiff ordinary differential equations”  In _SIAM review_ 21.1 SIAM, 1979, pp. 1–17
  * [27] Ilya M Sobol’, Danil Asotsky, Alexander Kreinin and Sergei Kucherenko  “Construction and comparison of high-dimensional Sobol’generators”  In _Wilmott_ 2011.56 Wiley Online Library, 2011, pp. 64–79
  * [28] Gustaf Söderlind, Laurent Jay and Manuel Calvo  “Stiffness 1952–2012: Sixty years in search of a definition”  In _BIT Numerical Mathematics_ 55.2 Springer, 2015, pp. 531–558
  * [29] Sifan Wang, Yujun Teng and Paris Perdikaris  “Understanding and mitigating gradient pathologies in physics-informed neural networks”  In _arXiv preprint arXiv:2001.04536_ , 2020
  * [30] Gerhard Wanner and Ernst Hairer  “Solving ordinary differential equations II”  Springer Berlin Heidelberg, 1996
  * [31] Michael Wetter, Wangda Zuo, Thierry S Nouidui and Xiufeng Pang  “Modelica buildings library”  In _Journal of Building Performance Simulation_ 7.4 Taylor & Francis, 2014, pp. 253–270


# Two Classes of Power Mappings with Boomerang Uniformity 2

Zhen Li<EMAIL_ADDRESS>Haode Yan<EMAIL_ADDRESS>School of
Mathematics, Southwest Jiaotong University, Chengdu, 610031, China

###### Abstract

Let $q$ be an odd prime power. Let $F_{1}(x)=x^{d_{1}}$ and
$F_{2}(x)=x^{d_{2}}$ be power mappings over ${\mathrm{GF}}(q^{2})$, where
$d_{1}=q-1$ and $d_{2}=d_{1}+\frac{q^{2}-1}{2}=\frac{(q-1)(q+3)}{2}$. In this
paper, we study the the boomerang uniformity of $F_{1}$ and $F_{2}$ via their
differential properties. It is shown that, the boomerang uniformity of $F_{i}$
($i=1,2$) is 2 with some conditions on $q$.

###### keywords:

Power mapping , Differential uniformity , Boomerang uniformity

###### MSC:

94A60 , 11T06

##  1 Introduction

Differential cryptanalysis [1] is one of the most fundamental analysis method
to attack block cipher, which has attracted extensive researches for 30 years.
In order to measure the ability of a given function to resist differential
attack, Nyberg introduced the concept of differential uniformity in [25]. Let
${\mathrm{GF}}(q)$ denote the finite field containing $q$ elements and
${\mathrm{GF}}(q)^{*}={\mathrm{GF}}(q)\setminus\\{0\\}$. For a function $F$
from ${\mathrm{GF}}(q)$ to itself, the differential uniformity of $F$ is

$\Delta_{F}=\max_{a\in{\mathrm{GF}}(q)^{*}}\max_{b\in{\mathrm{GF}}(q)}\delta_{F}(a,b),$

where $\delta_{F}(a,b)=\\#\\{x\in{\mathrm{GF}}(q)~{}|~{}F(x+a)-F(x)=b\\}$. The
lower the quantity of $\Delta_{F}$ is, the stronger the ability of function
$F$ resisting differential attack. If $\Delta_{F}=1$, then $F$ is called
perfect nonlinear (PN). Known results on PN functions were presented in [6,
10, 11, 34]. If $\Delta_{F}=2$, then $F$ is called almost perfect nonlinear
(APN). For the known results on APN functions, the readers are referred to [5,
6, 12, 13, 15, 17, 18, 24, 25, 35].

Power mappings (i.e., monomials) with low differential uniformity serve as
good candidates for the design of S-boxes not only because of their strong
resistance to differential attacks but also for the usually low implementation
cost in hardware. For any power mappings $F(x)=x^{d}$, we have
$\delta_{F}(a,b)=\delta_{F}(1,\frac{b}{a^{d}})$ for arbitrary
$(a,b)\in{\mathrm{GF}}(q)^{*}\times{\mathrm{GF}}(q)$. Hence the differential
properties of the power function $F$ are completely determined by
$\delta_{F}(1,b)$ when $b$ varies through ${\mathrm{GF}}(q)$. In [2],
Blondeau, Canteaut and Charpin defined the differential spectrum of a power
function. Let
$\omega_{i}=\\#\\{b\in{\mathrm{GF}}(q)~{}|~{}\delta_{F}(1,b)=i\\}$ for $0\leq
i\leq\Delta_{F}$, where $\Delta_{F}$ is the differential uniformity of $F$.
The differential spectrum of $F$ is defined as the multiset

$\mathbb{S}=\\{\omega_{i}~{}|~{}0\leq i\leq\Delta_{F},\omega_{i}>0\\}.$

The differential spectrum is an important concept of cryptographic functions.
The differential spectrum of any PN function is
$\mathbb{S}=\\{\omega_{1}=q\\}$, while the differential spectrum of any APN
function over even characteristic finite field ${\mathrm{GF}}(2^{n})$ is
$\mathbb{S}=\\{\omega_{0}=\omega_{2}=2^{n-1}\\}$. For the results of power
mappings with known differential spectra, the readers are referred to [2, 3,
8, 14, 21, 22, 27, 28, 29, 30, 31, 32] and their references for more
information. In 1999, Wagner introduced a new cryptanalysis method against
block ciphers, namely, the boomerang attack [26]. It can be regard as a
generalization of the differential attack, and it allows new avenues of attack
for many ciphers previously deemed safe from differential cryptanalysis.
Nowadays, the boomerang distinguisher is one of the most powerful
distinguishers for identifying weaknesses in block ciphers. The boomerang
distinguisher has been recently deemed another highly important distinguisher
of block ciphers. In order to study this attack method in detailed, Cid et al.
introduced the Boomerang Connectivity Table (BCT) in EUROCRYPT 2018 [9]. To
quantify the resistance of a function against the boomerang attack, Boura and
Canteaut introduced the concept of boomerang uniformity [4]. Later in [20],
Li, Qu, Sun and Li generalized the definition of the boomerang uniformity. Let
$F$ be a function defined on ${\mathrm{GF}}(p^{n})$. For any
$a,b\in{\mathrm{GF}}(p^{n})^{*}$, denote by $\beta_{F}(a,b)$ the number of the
solutions $(x,y)$ in ${\mathrm{GF}}(p^{n})\times{\mathrm{GF}}(p^{n})$ of the
equation system

$\displaystyle\left\\{\begin{array}[]{ll}F(x)-F(y)&=b,\\\
F(x+a)-F(y+a)&=b.\end{array}\right.$

The boomerang uniformity of $F$ is defined as

$\beta_{F}=\max_{a,b\in{\mathrm{GF}}(p^{n})^{*}}\\{\beta_{F}(a,b)\\}.$

The lower the quantity of $\beta_{F}$ is, the stronger the ability of function
$F$ resisting boomerang attack. Recent progress on cryptographic functions
with known boomerang uniformity can refer to the survey written by Mesnager,
Mandal and Msahli [23].

The boomerang uniformity of power mappings attracts lots of attention. When
$F(x)=x^{d}$ is a power function, we have
$\beta_{F}=\max\limits_{b\in{\mathrm{GF}}(p^{n})^{*}}\\{\beta_{F}(1,b)\\}$
similarly. The known results are introduced as follows. For power functions
over even characteristic finite fields ${\mathrm{GF}}(2^{n})$, the boomerang
uniformity of the inverse function $x^{2^{n}-2}$ and the Gold function
$x^{2^{k}+1}$ were presented in [4]. Zha and Hu studied the boomerang
uniformity of power permutations $x^{2^{k}-1}$ in [33]. Calderini and Villa
studied the Bracken-Leader function $x^{2^{2k}+2^{k}+1}$, where $n=4k$ and $k$
is odd. They showed that the boomerang uniformity of the Bracken-Leader
function is upper bounded by 24 [7]. Recently, Hasan, Pal and Stanica studied
the power mapping $x^{q-1}$ over ${\mathrm{GF}}(q^{2})$, where $q=2^{m}$ is a
power of $2$. It was shown that the boomerang uniformity is $2$ if $m$ is
even, and is $4$ if $m$ is odd [16]. For odd characteristic finite fields,
there are few results. In [19], Jiang et al. determined the boomerang
uniformity of $x^{p^{n}-2}$ over ${\mathrm{GF}}(p^{n})$ and
$x^{\frac{3^{n}+3}{2}}$ over ${\mathrm{GF}}(3^{n})$ with odd $n$.

We focus on two power mappings over odd characteristic finite fields. Let
$F_{1}(x)=x^{q-1}$ and $F_{2}(x)=x^{\frac{(q-1)(q+3)}{2}}$ be power mappings
over ${\mathrm{GF}}(q^{2})$. In this paper, we present that the boomerang
uniformity of $F_{i}$ ($i=1,2$) is $2$. The rest of this paper is organized as
follows. In Section 2, we introduce some notation and the differential
properties of $F_{1}$ and $F_{2}$. In section 3, we investigate the boomerang
uniformity of $F_{1}$ and $F_{2}$ via their differential properties. Section 4
concludes this paper.

##  2 Preliminaries

We begin this section by fixing some notation which will be used throughout
this paper unless otherwise stated.

  * 1.

$q$ is an odd prime power.

  * 2.

Let $F_{1}(x)=x^{d_{1}}$ and $F_{2}(x)=x^{d_{2}}$ be power mappings over
${\mathrm{GF}}(q^{2})$, where $d_{1}=q-1$ and
$d_{2}=d_{1}+\frac{q^{2}-1}{2}=\frac{(q-1)(q+3)}{2}$.

  * 3.

$\Delta_{i}(x)=F_{i}(x+1)-F_{i}(x)=(x+1)^{d_{i}}-x^{d_{i}}$, where $i=1,2$.

  * 4.

For any $b\in{\mathrm{GF}}(q^{2})$, let
$\Delta^{-1}_{i}(b)=\\{x~{}|~{}\Delta_{i}(x)=b\\}$ and
$\delta_{i}(b)=\\#\Delta^{-1}_{i}(b)$, where $i=1,2$.

In the following, we give the differential properties of $F_{1}$ and $F_{2}$
with some conditions on $q$, which will be used in the sequel.

######  Lemma 2.1.

Let $q\not\equiv 2\pmod{3}$. We have,

  * 1.

$\Delta^{-1}_{1}(0)={\mathrm{GF}}(q)\setminus\\{0,-1\\}$, $\delta_{1}(0)=q-2$;

  * 2.

$\Delta^{-1}_{1}(1)=\\{0\\},\Delta^{-1}_{1}(-1)=\\{-1\\},\delta_{1}(1)=\delta_{1}(-1)=1$;

  * 3.

$\delta_{1}(b)\leq 2$ for all $b\neq 0,\pm 1$.

###### Proof.

We mainly study the derivative equation

$\Delta_{1}(x)=(x+1)^{q-1}-x^{q-1}=b$ (1)

for arbitrary but fixed $b\in{\mathrm{GF}}(q^{2})$. Clearly, $\Delta_{1}(0)=1$
and $\Delta_{1}(-1)=-1$. We assume that $b=0$ first, then $x\neq 0,-1$. By
(1), we obtain

$(1+\frac{1}{x})^{q-1}=1,$

then $x\in{\mathrm{GF}}(q)$. Moreover, any
$x\in{\mathrm{GF}}(q)\setminus\\{0,-1\\}$ gives $\Delta_{1}(x)=0$. We conclude
that $\Delta^{-1}_{1}(0)={\mathrm{GF}}(q)\setminus\\{0,-1\\}$ and
$\delta_{1}(0)=q-2$.

Next we assume that $b\neq 0$. Multiplying $x(x+1)$ on each side of (1) gives

$x-bx(x+1)=x^{q}.$ (2)

Then $x^{q}-b^{q}x^{q}(x^{q}+1)=x$. This with (2) leads to

$x^{2}+(1-\frac{2}{b})x+\frac{b^{q-1}+1-b^{q}}{b^{q+1}}=0,$ (3)

which is a quadratic equation of $x$. Therefore, $\delta(b)\leq 2$ for all
$b\neq 0,\pm 1$.

When $b=1$, it is obvious that $x=0$ is a solution of (1). When $x\neq 0,-1$,
the equation (2) becomes

$x^{q-2}=-1.$

Note that $\gcd(2(q-2),q^{2}-1)=2$ since $q\not\equiv 2\pmod{3}$. We obtain
$x=-1$, which is not considered. Then $\Delta^{-1}(1)=\\{0\\}$ and
$\delta_{1}(1)=1$.

Since $q-1$ is even, $x\in\Delta^{-1}_{1}(b)$ if and only if
$-x-1\in\Delta^{-1}_{1}(-b)$, then $\delta_{1}(b)=\delta_{1}(-b)$ for any
$b\in{\mathrm{GF}}(q^{2})$. The desired results follows. ∎

######  Lemma 2.2.

Let $q\not\equiv 2\pmod{3}$ and $q\equiv 3\pmod{4}$. We have

  * 1.

$\Delta^{-1}_{2}(0)={\mathrm{GF}}(q)\setminus\\{0,-1\\}$, $\delta_{2}(0)=q-2$;

  * 2.

$\Delta^{-1}_{2}(1)=\\{0\\},\Delta^{-1}_{2}(-1)=\\{-1\\},\delta_{2}(1)=\delta_{2}(-1)=1$;

  * 3.

$\delta_{2}(b)\leq 2$ for all $b\neq 0,\pm 1$.

###### Proof.

We mainly study the derivative equation

$\Delta_{2}(x)=(x+1)^{d_{2}}-x^{d_{2}}=b$ (4)

for arbitrary but fixed $b\in{\mathrm{GF}}(q^{2})$, where
$d_{2}=\frac{(q-1)(q+3)}{2}$ is even. Clearly, $\Delta_{2}(0)=1$ and
$\Delta_{2}(-1)=-1$. When $b=0$, $x\neq 0,-1$. By (4), we obtain

$(1+\frac{1}{x})^{d_{2}}=1,$

then $x\in{\mathrm{GF}}(q)$ since $\gcd(d_{2},q^{2}-1)=q-1$. We conclude that
$\Delta^{-1}_{2}(0)={\mathrm{GF}}(q)\setminus\\{0,-1\\}$ and
$\delta_{2}(0)=q-2$.

Next we assume that $b\neq 0$. Let $\chi(x)=x^{\frac{q^{2}-1}{2}}$ be the
quadratic multiplicative character of ${\mathrm{GF}}(q^{2})^{*}$. For $x\neq
0,-1$, we distinguish the following two cases.

Case I. $\chi(x+1)=\chi(x)$. Then (4) becomes

$(x+1)^{q-1}-x^{q-1}=b\chi(x).$ (5)

Moreover,

$b\chi(x)x(x+1)=x-x^{q}.$ (6)

Note that $(x-x^{q})^{q}=-(x-x^{q})$, $(x-x^{q})^{q-1}=-1$, then
$\chi(x-x^{q})=(x-x^{q})^{\frac{(q-1)(q+1)}{2}}=1$ since $q\equiv 3\pmod{4}$.
From (6) we conclude that if this case contributes solutions, then
$\chi(b)=1$.

By Lemma 2.1, (5) has at most two solutions for each value of $\chi(x)$. Note
that $x=x_{0}$ is a solution of $(x+1)^{q-1}-x^{q-1}=b$ if and only if
$x=-x_{0}-1$ is a solution of $(x+1)^{q-1}-x^{q-1}=-b$. However,
$\chi(-x_{0}-1)=\chi(x_{0}+1)=\chi(x_{0})$. This implies that (5) cannot have
solution for both two values of $\chi(x)$ simultaneously. This case
contributes at most two solutions.

Case II. $\chi(x+1)=-\chi(x)$. Then (4) becomes

$(x+1)^{q-1}+x^{q-1}=-b\chi(x).$ (7)

Multiplying $x(x+1)$ on both sides of (7), we have

$2x^{q+1}+x^{q}+x=-b\chi(x)x(x+1).$ (8)

Note that $-b\chi(x)x(x+1)\in{\mathrm{GF}}(q)$ and $\chi(x(x+1))=-1$, we
conclude that, if this case contributes solutions, then $\chi(b)=-1$.
Moreover, we obtain

$x^{q}=\frac{-b\chi(x)x(x+1)-x}{1+2x}$ (9)

from (8). From (9), we obtain that

$x^{2}+x+\frac{b^{q+1}+\chi(x)(b^{q}+b)}{-\chi(x)b(b^{q+1}-4)}=0.$ (10)

We mention that $b^{q+1}\neq 4$. Otherwise, $b^{q+1}=4$, then
$b^{\frac{q^{2}-1}{2}}=4^{\frac{q-1}{2}}=2^{q-1}=1$, which contradicts
$\chi(b)=-1$. Note that both $b^{q+1}+\chi(x)(b^{q}+b)$ and $b^{q+1}-4$ are in
${\mathrm{GF}}(q)$. For a given $\chi(x)$, the equation (10) is a quadratic
equation of $x$, the discriminant is
$\frac{(b+2\chi(x))^{2}b^{q-1}}{b^{q+1}-4}$, which is a square element in
${\mathrm{GF}}(q^{2})$. Let $x_{1}$ and $x_{2}$ be solutions of (10), then
$\chi(x_{1}x_{2})=\chi(\frac{b^{q+1}+\chi(x)(b^{q}+b)}{-\chi(x)b(b^{q+1}-4)})=-1$
since $\chi(b)=-1$. Only one of $x_{1}$ and $x_{2}$ is a square element, and
the other one is a nonsquare element. For each value of $\chi(x)$, there is at
most one solution. This case contributes at most two solutions.

For $b\neq 0,\pm 1$, the solutions are all in the above two cases. According
to the value of $\chi(b)$, we know that Cases I and II cannot have solution
simultaneously. Thus we obtain $\delta_{2}(b)\leq 2$ for $b\neq 0,\pm 1$.

At last, we consider $b=\pm 1$. $\Delta_{2}(0)=1$ and $\Delta_{2}(-1)=-1$ are
obvious. It can be checked that (4) has no solution in the two cases. Then
$\Delta^{-1}_{2}(1)=\\{0\\},\Delta^{-1}_{2}(-1)=\\{-1\\}$, and
$\delta_{2}(1)=\delta_{2}(-1)=1$. We complete the proof.

∎

##  3 The Boomerang Uniformity of $F_{1}$ and $F_{2}$

In this section, we determine the boomerang uniformity of $F_{1}(x)=x^{q-1}$
and $F_{2}(x)=x^{\frac{(q-1)(q+3)}{2}}$ over ${\mathrm{GF}}(q^{2})$. We have
the following theorems.

######  Theorem 3.1.

Let $q\not\equiv 2\pmod{3}$. The boomerang uniformity of $F_{1}(x)=x^{q-1}$
over ${\mathrm{GF}}(q^{2})$ is 2.

###### Proof.

To determine the boomerang uniformity of $F_{1}$, we mainly study the equation
system

$\displaystyle\left\\{\begin{array}[]{ll}x^{q-1}-y^{q-1}&=b\\\
(x+1)^{q-1}-(y+1)^{q-1}&=b\end{array}\right.$ (13)

for arbitrary but fixed $b\in{\mathrm{GF}}(q^{2})^{*}$. It is obvious that the
equation system (13) is equivalent to the following one,

$\displaystyle\left\\{\begin{array}[]{ll}x^{q-1}-y^{q-1}&=b\\\
(x+1)^{q-1}-(y+1)^{q-1}&=b\\\ \Delta_{1}(x)&=\Delta_{1}(y),\end{array}\right.$
(17)

where $\Delta_{1}(z)=(z+1)^{q-1}-z^{q-1}$ we defined before. Let
$\Delta_{1}(x)=\Delta_{1}(y)=c$ for some $c\in{\mathrm{GF}}(q^{2})$. By Lemma
2.1, we know that the equation $\Delta_{1}(x)=c$ (respectively,
$\Delta_{1}(y)=c$) has $0,1,2$ or $q-2$ roots. More precisely, define

$\Omega_{i}=\\{c\in{\mathrm{GF}}(q^{2})~{}|~{}\delta_{1}(c)=i~{}\\}$

for $i=0,1,2,q-2$. We mention that $\Omega_{q-2}=\\{0\\}$ and
$\Omega_{0}\cup\Omega_{1}\cup\Omega_{2}={\mathrm{GF}}(q^{2})^{*}$. When $c$
runs through ${\mathrm{GF}}(q^{2})$, we discuss in the following four cases.

Case 1. $c\in\Omega_{0}$. By the definition of $\Omega_{0}$, $\Delta_{1}(x)=c$
(respectively, $\Delta_{1}(y)=c$) has no solution in ${\mathrm{GF}}(q^{2})$.
Hence this case cannot contribute solution to (17).

Case 2. $c\in\Omega_{1}$. Let $\Delta_{1}^{-1}(c)=\\{z_{0}\\}$. In this case,
$\Delta_{1}(x)=\Delta_{1}(y)=c$ has one solution $(x,y)=(z_{0},z_{0})$. This
case cannot contribute solution to (17) since $b\neq 0$.

Case 3. $c\in\Omega_{q-2}$, i.e., $c=0$. By Lemma 2.1,
$\Delta_{1}^{-1}(0)={\mathrm{GF}}(q)\setminus\\{0,-1\\}$. When
$\Delta_{1}(x)=\Delta_{1}(y)=0$, i.e.,
$x,y\in{\mathrm{GF}}(q)\setminus\\{0,-1\\}$, we have $b=x^{q-1}-y^{q-1}=0$,
which is a contradiction. Hence this case cannot contribute solution to (17).

Case 4. $c\in\Omega_{2}$. We first consider the equation system (17). We know
that $x,y\neq 0,-1$. From the first two equations of (17), we obtain

$x^{q}y-xy^{q}=bxy$

and

$(x+1)^{q}(y+1)-(x+1)(y+1)^{q}=b(x+1)(y+1).$

Then

$b(x+y+1)=b(x+1)(y+1)-bxy=x^{q}+y-x-y^{q}.$

We have

$b(x+y+1)+(b(x+y+1))^{q}=0.$ (18)

Let $z_{1}$ and $z_{2}$ be the two distinct solutions of $\Delta_{1}(z)=c$.
Then the possible solutions of (17) are $(x,y)=(z_{1},z_{2})$ and
$(x,y)=(z_{2},z_{1})$ since $b\neq 0$. Note that
$z^{q-1}_{1}-z^{q-1}_{2}=-(z^{q-1}_{2}-z^{q-1}_{1})$, $(z_{1},z_{2})$ and
$(z_{2},z_{1})$ cannot be solutions of (17) simultaneously. Without loss of
generalization, we assume that $(z_{1},z_{2})$ is a solution of (17) for such
$c$ and the fixed $b$. By the proof of Lemma 2.1, equation (2) shows that
$z^{q-1}_{i}=1-c(z_{i}+1)$ for $i=1,2$. We have

$b=z^{q-1}_{1}-z^{q-1}_{2}=c(z_{1}-z_{2}).$ (19)

The equation (3) shows that $z_{1}+z_{2}=-1+\frac{2}{c}$ and
$z_{1}z_{2}=\frac{c^{q-1}+1-c^{q}}{c^{q+1}}$. Combining with (18) and (19), we
obtain the following relationship between $b$ and $c$.

$\displaystyle\left\\{\begin{array}[]{ll}(\frac{b}{c})^{q-1}&=-1\\\
(\frac{b}{c})^{2}&=1-\frac{4}{c^{q+1}}.\end{array}\right.$ (22)

For $c^{\prime}\in\Omega_{2}$, $c^{\prime}\neq c$, if
$\Delta_{1}(x)=\Delta_{1}(y)=c^{\prime}$ contributes another solution of (17),
we know that $c^{\prime}$ and $b$ also satisfy (22). Let $c^{\prime}=\alpha c$
for some $1\neq\alpha\in{\mathrm{GF}}(q^{2})$. Then we obtain $\alpha^{q-1}=1$
from the first equation of (22). Moreover, from the second equation of (22),
we have

$c^{2}(1-\frac{4}{c^{q+1}})=(\alpha c)^{2}(1-\frac{4}{(\alpha c)^{q+1}}).$

This with $\alpha^{q-1}=1$ leads to $\alpha^{2}=1$, then $\alpha=-1$, i.e.,
$c^{\prime}=-c$. Only $-c\in\Omega_{2}$ can contribute solution of (17) for
the given $b$. More precisely, when $z_{1}$ and $z_{2}$ are distinct solutions
of $\Delta_{1}(z)=c$, the solutions of $\Delta_{1}(z)=-c$ are $-z_{1}-1$ and
$-z_{2}-1$. One can easily check that, when $(z_{1},z_{2})$ is a solution of
(17), so is $(-z_{2}-1,-z_{1}-1)$. We conclude that, for all
$b\in{\mathrm{GF}}(q^{2})$, the maximum of the number of solutions of (13) is
2. We complete the proof. ∎

######  Theorem 3.2.

Let $q\not\equiv 2\pmod{3}$ and $q\equiv 3\pmod{4}$. The boomerang uniformity
of $F_{2}(x)=x^{d_{2}}$ over ${\mathrm{GF}}(q^{2})$ is 2, where
$d_{2}=\frac{(q-1)(q+3)}{2}$.

###### Proof.

We consider the equivalent equation system

$\displaystyle\left\\{\begin{array}[]{ll}x^{d_{2}}-y^{d_{2}}&=b\\\
(x+1)^{d_{2}}-(y+1)^{d_{2}}&=b\\\ \Delta_{2}(x)&=\Delta_{2}(y)\\\
\end{array}\right.$ (26)

for arbitrary but fixed $b\in{\mathrm{GF}}(q^{2})^{*}$, where
$\Delta_{2}(z)=(z+1)^{d_{2}}-z^{d_{2}}$ we defined before. Similarly, we
define

$\Omega_{i}=\\{c\in{\mathrm{GF}}(q^{2})~{}|~{}\delta_{2}(c)=i\\}$

for $i=0,1,2,q-2$. We mention that $\Omega_{q-2}=\\{0\\}$ and
$\Omega_{0}\cup\Omega_{1}\cup\Omega_{2}={\mathrm{GF}}(q^{2})^{*}$. Let
$\Delta_{2}(x)=\Delta_{2}(y)=c$, we discuss the equation system (26) when $c$
runs through ${\mathrm{GF}}(q^{2})$. Similar to the proof of Theorem 3.1, we
can prove that when $c\in\Omega_{0}\cup\Omega_{1}\cup\Omega_{q-2}$, the
corresponding pair $(x,y)$ which comes from $\Delta_{2}(x)=\Delta_{2}(y)=c$ is
not a solution of (26). We only consider the case that $c\in\Omega_{2}$. Let
$z_{1}$ and $z_{2}$ be the solutions of $\Delta_{2}(z)=c$, similar to the
proof of Theorem 2.2, possible solutions are $(z_{1},z_{2})$ and
$(z_{2},z_{1})$. We always have,

$b^{2}=(z_{1}^{d_{2}}-z_{2}^{d_{2}})^{2}=(\chi(z_{1})z^{q-1}_{1}-\chi(z_{2})z^{q-1}_{2})^{2}.$
(27)

According to Lemma 2.2, we discuss the relationship of $b$ and $c$ in the
following two cases.

Case 1. $z_{1}$ and $z_{2}$ belong to Case I of Lemma (2.2), i.e.,
$\chi(z_{1}+1)=\chi(z_{1})=\chi(z_{2}+1)=\chi(z_{2})$. In this case, $z_{1}$
and $z_{2}$ are the two solutions of

$(z+1)^{q-1}-z^{q-1}=c\chi(z).$

For a given value of $\chi(z_{1})$, by the proof of Lemma 2.1, we obtain
$z^{q-1}_{i}=1-c\chi(z_{1})(z_{i}+1)$ ($i=1,2$) from (2), and
$z_{1}+z_{2}=-1+\frac{2\chi(z_{1})}{c}$,
$z_{1}z_{2}=\frac{c^{q-1}+1-c^{q}\chi(z_{1})}{c^{q+1}}$ from (3). From (27) we
have

$b^{2}=(\chi(z_{1})z^{q-1}_{1}-\chi(z_{2})z^{q-1}_{2})^{2}=c^{2}(z_{1}-z_{2})^{2}=c^{2}-\frac{4}{c^{q-1}}.$

Case 2. $z_{1}$ and $z_{2}$ belong to Case II of Lemma (2.2), i.e.,
$\chi(z_{1}+1)=\chi(z_{2})=1,\chi(z_{1})=\chi(z_{2}+1)=-1$ (do not consider
the order of $z_{1}$ and $z_{2}$). Based on (10), $z_{1}$ and $z_{2}$ satisfy
of $z_{1}^{2}+z_{1}+\frac{c^{q+1}-c^{q}-c}{c(c^{q+1}-4)}=0$ and
$z_{2}^{2}+z_{2}+\frac{c^{q+1}+c^{q}+c}{-c(c^{q+1}-4)}=0$, respectively. We
have

$z_{1}=\frac{1}{2}(-1+\frac{(c-2)c^{\frac{q-1}{2}}}{\sqrt{c^{q+1}-4}})~{}~{}~{}\mathrm{or}~{}~{}~{}\frac{1}{2}(-1-\frac{(c-2)c^{\frac{q-1}{2}}}{\sqrt{c^{q+1}-4}}),$
$z_{2}=\frac{1}{2}(-1+\frac{(c+2)c^{\frac{q-1}{2}}}{\sqrt{c^{q+1}-4}})~{}~{}~{}\mathrm{or}~{}~{}~{}\frac{1}{2}(-1-\frac{(c+2)c^{\frac{q-1}{2}}}{\sqrt{c^{q+1}-4}}).$

There are four possible pairs of $(z_{1},z_{2})$. Moreover, by (9) we have

$z_{1}^{q-1}=\frac{c(z_{1}+1)-1}{1+2z_{1}}$

and

$z_{2}^{q-1}=\frac{-c(z_{2}+1)-1}{1+2z_{2}},$

then

$b=\chi(z_{1})z^{q-1}_{1}-\chi(z_{2})z^{q-1}_{2}=\frac{c(z_{2}-z_{1})+2(1+z_{1}+z_{2})}{(1+2z_{1})(1+2z_{2})}.$
(28)

Plugging the possible pairs $(z_{1},z_{2})$ into (28), we obtain

$b^{2}=c^{2}-\frac{4}{c^{q-1}},$

when
$(z_{1},z_{2})=(\frac{1}{2}(-1+\frac{(c-2)c^{\frac{q-1}{2}}}{\sqrt{c^{q+1}-4}}),\frac{1}{2}(-1-\frac{(c+2)c^{\frac{q-1}{2}}}{\sqrt{c^{q+1}-4}}))$
and
$(z_{1},z_{2})=(\frac{1}{2}(-1-\frac{(c-2)c^{\frac{q-1}{2}}}{\sqrt{c^{q+1}-4}}),\frac{1}{2}(-1+\frac{(c+2)c^{\frac{q-1}{2}}}{\sqrt{c^{q+1}-4}}))$.
The other two pairs give $b=0$ and we discard them.

In the above two cases, we always have $b^{2}=c^{2}-\frac{4}{c^{q-1}}$. Then
$(\frac{b}{c})^{2}=1-\frac{4}{c^{q+1}}\in{\mathrm{GF}}(q)$. We obtain
$(\frac{b}{c})^{q-1}=1$ or $-1$. In Case 1, we claim that
$(\frac{b}{c})^{q-1}=-1$. The proof is the same as the determination of the
first equation in (22). In what follows, we prove that
$(\frac{b}{c})^{q-1}=-1$ in Case 2. In Case 2, we have

$c=(z_{1}+1)^{d_{2}}-z^{d_{2}}_{1}=(z_{1}+1)^{q-1}+z^{q-1}_{1}=\frac{2z^{q+1}_{1}+z^{q}_{1}+z_{1}}{z_{1}(z_{1}+1)},$

and

$c^{q+1}=(\frac{2z^{q+1}_{1}+z^{q}_{1}+z_{1}}{z_{1}(z_{1}+1)})^{q+1}=\frac{(2z^{q+1}_{1}+z^{q}_{1}+z_{1})^{q+1}}{z^{q+1}_{1}(z_{1}+1)^{q+1}}.$

Then

$(\frac{b}{c})^{2}=1-\frac{4}{c^{q+1}}=1-\frac{4z^{q+1}_{1}(z_{1}+1)^{q+1}}{(2z^{q+1}_{1}+z^{q}_{1}+z_{1})^{q+1}}=(\frac{z^{q}_{1}-z_{1}}{2z^{q+1}_{1}+z^{q}_{1}+z_{1}})^{2}.$

Note that $2z^{q+1}_{1}+z^{q}_{1}+z_{1}\in{\mathrm{GF}}(q)$ and
$z^{q}_{1}-z_{1}\notin{\mathrm{GF}}(q)$, we assert that
$(\frac{b}{c})^{q-1}=-1$. By the discussion above, we have

$\displaystyle\left\\{\begin{array}[]{ll}(\frac{b}{c})^{q-1}&=-1\\\
(\frac{b}{c})^{2}&=1-\frac{4}{c^{q+1}},\end{array}\right.$

which is the same as the relationship shown in the proof of Theorem 3.1. Then
we can similarly prove that, the maximum number of the solutions of (26) for
all $b\in{\mathrm{GF}}(q^{2})$ is 2. The desired result follows. ∎

##  4 Concluding and Remarks

In this paper, two classes of power mappings with boomerang uniformity 2 were
presented. We used their differential properties to study their boomerang
uniformity. Our approach can be used to study the boomerang uniformity of
other power mappings. It would be interesting to find more infinite classes of
monomials with low boomerang uniformity.

## References

  * [1] E. Biham and A. Shamir,  Differential cryptanalysis of DES-like cryptosystems,  J. Cryptology, 4 (1991), 3–72.
  * [2] C. Blondeau, A. Canteaut and P. Charpin,  Differential properties of power functions,  Int. J. Inf. Coding Theory, 1 (2010), 149–170.
  * [3] C. Blondeau, A. Canteaut and P. Charpin,  Differential properties of ${x\mapsto x^{2^{t}-1}}$,  IEEE Trans. Inform. Theory, 57 (2011), 8127–8137.
  * [4] C. Boura and A. Canteaut,  On the boomerang uniformity of cryptographic sboxes,  IACR Trans. Symmetric Cryptology, 2018 (2018), 290–310.
  * [5] K. Browning, J. Dillon, M. McQuistan, and J. Wolfe,  An APN permutation in dimension six,  Post Proceedings of the 9th International Conference on Finite Fields and their Applications, Contemp. Math. Amer. Math. Soc., 518 (2010), 33–42.
  * [6] L. Budaghyan and T. Helleseth,  New perfect nonlinear multinomials over $F_{p^{2k}}$ for any odd prime $p$,  Sequences and their applications—SETA 2008, Lecture Notes in Comput. Sci., 5203 (2008), 403–414.
  * [7] M. Calderini and I. Villa,  On the boomerang uniformity of some permutation polynomials,  Cryptogr. Commun., 12 (2020), 1161–1178.
  * [8] S. Choi, S. Hong, J. No, and H. Chung,  Differential spectrum of some power functions in odd prime characteristic,  Finite Fields Appl., 21 (2013), 11–29.
  * [9] C. Cid, T. Huang, T. Peyrin, Y. Sasaki, and L. Song,  Boomerang connectivity table: a new cryptanalysis tool,  Advances in cryptology—EUROCRYPT 2018. Part II, Lecture Notes in Comput. Sci., 10821 (2018), 683–714.
  * [10] R. Coulter and R. Matthews,  Planar functions and planes of Lenz-Barlotti class II,  Des. Codes Cryptogr., 10 (1997), 167–184.
  * [11] C. Ding and J. Yuan,  A family of skew Hadamard difference sets,  J. Combin. Theory Ser. A, 113(2006), 1526–1535.
  * [12] H. Dobbertin,  Almost perfect nonlinear power functions on GF($2^{n}$): the Welch case,  IEEE Trans. Inform. Theory, 45(1999), 1271–1275.
  * [13] H. Dobbertin,  Almost perfect nonlinear power functions on $\rm{GF}(2^{n})$: the Niho case,  Inform. and Comput., 151(1999), 57–72.
  * [14] H. Dobbertin, T. Helleseth, P. Kumar, and H. Martinsen,  Ternary m-sequences with three-valued cross-correlation function : new decimations of Welch and Niho type,  IEEE Trans. Inform. Theory, 47(2001), 1473–1481.
  * [15] H. Dobbertin, D. Mills, E. Müller, A. Pott, and W. Willems,  APN functions in odd characteristic,  Discrete Math., 267(2003), 95–112.
  * [16] S. Hasan, M. Pal, and P. Stănică,  Boomerang uniformity of a class of power maps,  Des. Codes Cryptogr., 89(2021), 2627–2636.
  * [17] T. Helleseth, C. Rong, and D. Sandberg,  New families of almost perfect nonlinear power mappings,  IEEE Trans. Inform. Theory, 45(1999), 474–485.
  * [18] T. Helleseth and D. Sandberg,  Some power mappings with low differential uniformity,  Appl. Algebra Engrg. Comm. Comput., 8(1997), 363–370.
  * [19] S. Jiang, K. Li, Y. Li, and L. Qu,  Differential and boomerang spectrums of some power permutations,  Cryptogr. Commun., (2021).
  * [20] K. Li, L. Qu, B. Sun, and C. Li,  New results about the boomerang uniformity of permutation polynomials,  IEEE Trans. Inform. Theory, 65(2019), 7542–7553.
  * [21] N. Li, Y. Wu, X. Zeng, and X. Tang,  On the differential spectrum of a class of power functions over finite fields,  (2020).
  * [22] Z. Li and H. Yan,  Differential spectra of a class of power permutations with niho exponents,  Adv. Math. Commun., (2021).
  * [23] S. Mesnager, B. Mandal, and M. Msahli,  Survey on recent trends towards generalized differential and boomerang uniformities,  Cryptogr. Commun., (2021).
  * [24] G. Ness and T. Helleseth,  A new family of ternary almost perfect nonlinear mappings,  IEEE Trans. Inform. Theory, 53(2007), 2581–2586.
  * [25] K. Nyberg,  Differentially uniform mappings for cryptography,  Advances in cryptology—EUROCRYPT ’93 (Lofthus, 1993), Lecture Notes in Comput. Sci., 765(1994), 55–64.
  * [26] D. Wagner,  The boomerang attack,  International Workshop on Fast Software Encryption, Lecture Notes in Comput. Sci., 1636(1999), 156–170.
  * [27] Y. Xia, X. Zhang, C. Li, and T. Helleseth,  The differential spectrum of a ternary power mapping,  Finite Fields Appl., 64(2020).
  * [28] M. Xiong and H. Yan,  A note on the differential spectrum of a 4-uniform power function,  Finite Fields Appl., 48(2017), 117–125.
  * [29] M. Xiong, H. Yan, and P. Yuan,  On a conjecture of differentially 8-uniform power function,  Des. Codes Cryptogr., 86(2018), 1601–1621.
  * [30] H. Yan and C. Li.  Differential spectra of a class of power permutations with characteristic 5,  Des. Codes Cryptogr., 89(2021), 1181–1191.
  * [31] H. Yan, Y. Xia, C. Li, T. Helleseth, M. Xiong, and J. Luo,  The differential spectrum of the power mapping $x^{p^{n}-3}$,  (2021).
  * [32] H. Yan, Z. Zhou, J. Weng, J. Wen, T. Helleseth, and Q. Wang,  Differential spectrum of Kasami power permutations over odd characteristic finite fields,  IEEE Trans. Inform. Theory, 65(2019), 6819–6826.
  * [33] Z. Zha and L. Hu,  The boomerang uniformity of power permutations $x^{2^{k}-1}$ over $F_{2^{n}}$,  Ninth International Workshop on Signal Design and its Applications in Communications (IWSDA), (2019).
  * [34] Z. Zha, G. Kyureghyan, and X. Wang,  Perfect nonlinear binomials and their semifields,  Finite Fields Appl., 15(2009), 125–133.
  * [35] Z. Zha and X. Wang,  Almost perfect nonlinear power functions in odd characteristic,  IEEE Trans. Inform. Theory, 57(2011), 4826–4832.