# Active Localization using Bernstein Distribution Functions††thanks: This
research was supported by the Office of Naval Research, grant N000142112091

Camilla Tabasso, Venanzio Cichella C. Tabasso and V. Cichella are with the
University of Iowa, Department of Mechanical Engineering, 52240, Iowa City,
IA, email: {camilla-tabasso<EMAIL_ADDRESS>
(May 2022)

###### Abstract

In this work, we present a framework that enables a vehicle to autonomously
localize a target based on noisy range measurements computed from RSSI data.
To achieve the mission objectives, we develop a control scheme composed of two
main parts: an estimator and a motion planner. At each time step, new
estimates of the target’s position are computed and used to generate and
update distribution functions using Bernstein polynomials. A metric of the
efficiency of the estimator is derived based on the Fisher Information Matrix.
Finally, the motion planning problem is formulated to react in real time to
new information about the target and improve the estimator’s performance.

##  I Introduction

The use of autonomous vehicles to localize targets has become increasingly
popular in recent years as a more efficient, safer and cheaper alternative to
localization performed manually by humans. A few notorious applications of
this problem include search and rescue missions, where aerial vehicles are
used to locate victims [1, 2], surveillance and monitoring for both military
and civilian applications [3], and localization of objects tagged with RFID
tags by robotic agents for resource management in industrial applications [4,
5].

In many of these scenarios it is imperative to perform the localization in the
least amount of time. This raises the question of how can the efficiency of
this process be improved. One approach, which has been greatly analyzed,
relies on computing the optimal sensor placement to observe the target.
Specifically, it has been shown that the Cramer Rao Lower Bound (CRLB), which
expresses a lower bound on the variance of a given estimator, can be used to
derive a closed-form solution of the optimal sensor-target geometry [6, 7, 8].
In the context of using autonomous vehicles to perform the estimation, this
translates to the problem of generating trajectories such that this optimal
geometry is achieved. However, due to the fact that the position of the target
is unknown, offline trajectory planning is not suitable for such missions.
This leads to the need of designing a motion planning algorithm which is able
to react in real time as new information becomes available, and to adjust the
vehicle’s path to obtain an effective sensor-target geometry.

In the literature, several solutions to this problem have been proposed. In
[9], a method based on maximizing the observability of a maneuvering target
given relative bearing measurements was analyzed. In [10], the authors
formulate an optimal trajectory generation algorithm that maximizes the Fisher
Information Matrix (FIM) for cooperative target localization. In [11], the
localization problem was solved for multiple UAVs by generating non-causal
trajectories that maximized the FIM. In [12], a similar problem was solved by
generating optimal trajectories that minimize the weighted sum of CRLB
diagonal entries. The main limitation of these methods is that the generated
trajectories only optimize the target-sensor geometry. However, many
applications that rely on target localization may require the optimization of
several objectives to be successful. In [13], the authors propose a multi-
objective motion planning framework for multiple targets localization and
pursuit, based on maximizing the predicted Bayesian Fisher Information Matrix.
The problem is then solved by designing an MPC-based closed-loop controller.
However, because of the inherent complexity of the problem we consider, which
leads to nonlinearities in both the objective function and constraints,
closed-loop solutions cannot be computed. Therefore, we formulate the planning
problem for target estimation as a general nonlinear optimal planning problem,
which is then addressed by finding an approximate solution using a direct
method. We employ the direct method based on Bernstein approximants that was
initially proposed in [14]. Bernstein polynomials have been shown to posses
nice geometric and numerical properties that are particularly useful for the
computation and enforcement of feasibility and safety constraints. Here, we
exploit these properties and algorithms to address the problem of estimating a
target’s position. For a more in depth discussion on the benefits of Bernstein
polynomial-based approaches as compared to other direct methods for trajectory
generation, the reader is referred to [15, 16, 17, 18, 19].

Bernstein polynomials have also attracted a lot of attention as a mean to
estimate distribution functions since they are known to yield smooth estimates
and do not suffer from the Gibbs phenomenon when estimating discontinuous
functions, see [20, 21, 22, 23] and references therein. In [22], the
asymptotic properties of the Bernstein estimator were explored theoretically,
while in [23], the authors investigate the estimator’s boundary properties.
Furthermore, they provide a comparison of this method with respect to
empirical distribution estimators and kernel estimators in terms of mean-
integrated squared error. In the context of our framework, the use of
Bernstein polynomials for the estimation ensures that the algorithm can be
readily integrated with the motion planner, and allows for an overall
reduction in the computational load of the framework.

Given these motivations, we build upon previous work [24] and develop a
localization framework that leverages the use of Bernstein polynomials for
both estimation and motion planning. This approach allows us to consider
multi-objective scenarios where the vehicle must efficiently estimate the
position of a target while minimizing mission execution time, actuation effort
to extend battery life, and guarantee feasibility for the whole duration of
the mission. To ensure robustness to estimator faults and guarantee
convergence to the target, we propose a time-triggered approach to replan
trajectories for the whole mission duration (rather than for a shorter
receding horizon) using up-to-date estimates. Lastly, the framework proposed
can be easily scaled to plan trajectories for multiple vehicles.

This paper is organized as follows: preliminary results on the properties of
Bernstein polynomials are presented in Section II. The Bernstein estimator is
introduced in Section III, while the formulation of the motion planning
algorithm is given in Section IV. Finally, numerical results are shown in
Section V, while conclusions are presented in Section VI.

##  II Preliminaries

###  II-A Framework

We consider a static target located at $p^{\star}\in\mathbb{R}^{2}$ and
equipped with a beacon that generates a signal with known characteristics.
Furthermore, we consider an autonomous agent, equipped with a receiver able to
detect the emitted signal, moving along a trajectory, namely
$p\in\mathbb{R}^{2}$. Examples of such signals include, but are not limited
to, electromagnetic waves, radio frequency signals and mechanical vibrations.

The objective of the target localization framework is to locate a target as
quickly as possible. The framework proposed can be divided in two main
algorithms: an estimator and a motion planner. The estimator provides an
estimate of the target’s position based on range information, as well as an
estimation of the cumulative distribution function (CDF) and probability
density function (PDF) of the measurements. On the other hand, the motion
planner is tasked to generate feasible trajectories designed to increase the
efficiency of the estimator. A schematics of the proposed framework is shown
in Figure 1.

Figure 1: Localization framework.

###  II-B Bernstein Polynomials

As they will be widely used throughout this paper, we formally introduce
Bernstein polynomials and their main properties below.

Consider the following $m$-th order Bernstein polynomial:

$\begin{split}p_{m}(t)&=\sum_{j=0}^{m}\bar{c}_{j,m}b_{j,m}(t),\qquad
t\in[t_{0},t_{f}],\\\ \end{split}$ (1)

where $\bar{c}_{0,m},\ldots,\bar{c}_{m,m}\in\mathbb{R}^{d}$ are Bernstein
polynomial coefficients,

$b_{j,m}(t)=\binom{m}{j}\frac{(t-t_{0})^{j}(t_{f}-t)^{m-j}}{(t_{f}-t_{0})^{m}}\,,\qquad
t\in[t_{0},t_{f}],$ (2)

for $j=0,\ldots,m$, is the $m$-th order Bernstein basis, and
$\binom{m}{j}=\frac{m!}{j!(m-j)!}$ is the binomial coefficient. The following
properties of Bernstein polynomials are used in this paper.

######  Property 1 (Differentiation and integration)

The derivatives of the Bernstein polynomial $p_{m}(t)$ can be easily computed
as

$\begin{split}\dot{p}_{m}(t)=\sum_{i=0}^{m}\sum_{j=0}^{m}&\bar{c}_{j,m}{D}_{j,i}b_{i,m}(t),\\\
\ddot{p}_{m}(t)=\sum_{i=0}^{m}\sum_{j=0}^{m}&\bar{c}_{j,m}{D}_{j,i}^{2}b_{i,m}(t),\\\
\end{split}$ (3)

where ${D}_{j,i}$ is the $(j,i)$-th entry of a square differentiation matrix
[25]. The definite integral of $p_{m}(t)$ is computed as

$\int_{t_{0}}^{t_{f}}p_{m}(t)dt=w\sum_{j=0}^{m}\bar{c}_{j,m}\,,\qquad
w=\frac{t_{f}-t_{0}}{m+1}\,.$ (4)

######  Property 2 (Arithmetic operations)

The sum (difference) of two $m$-th order Bernstein polynomials can be computed
by adding (subtracting) their Bernstein coefficients. The product between two
Bernstein polynomials of orders $m$ and $n$, with coefficients
$\bar{f}_{m}=[\bar{f}_{0,m},\dots,\bar{f}_{m,m}]$ and
$\bar{g}_{n}=[\bar{g}_{0,n},\dots,\bar{g}_{n,n}]$, respectively, is a
Bernstein polynomial of order $m+n$ with coefficients given as [19, Chapter 5]

$\bar{c}_{k,m+n}=\sum_{j=\max(0,k-n)}^{\min(m,k)}\frac{\Big{(}\begin{array}[]{c}{m}\\\
{j}\end{array}\Big{)}\Big{(}\begin{array}[]{c}{n}\\\
{k-j}\end{array}\Big{)}}{\Big{(}\begin{array}[]{c}{m+n}\\\
{k}\end{array}\Big{)}}\bar{f}_{j,m}\hskip 2.84526pt\bar{g}_{k-j,n}.$

Using Properties 4 and 2, the following functions can be expressed as
Bernstein polynomials:

$\begin{split}||\dot{p}_{m}(t)||^{2}=\sum_{j=0}^{2m}\bar{v}_{j,m}b_{j,m}(t),\\\
||\ddot{p}_{m}(t)||^{2}=\sum_{j=0}^{2m}\bar{a}_{j,m}b_{j,m}(t),\\\
||{p}_{m}(t)-p||^{2}=\sum_{j=0}^{2m}\bar{d}_{j,m}b_{j,m}(t).\end{split}$ (5)

The equations above represent the squared norms of the first and second time
derivatives and the squared distance from an arbitrary point $p$,
respectively. The coefficients $\bar{v}_{j,m}$, $\bar{a}_{j,m}$ and
$\bar{d}_{j,m}$, $\forall j\in\\{0,\ldots,2m\\}$ can be obtained from
algebraic manipulation of the Bernstein coefficients of ${p}_{m}(t)$.

######  Property 3 (End point values)

The initial and final values of a Bernstein polynomial are equal to its first
and last Bernstein coefficients, e.g., ${p}_{m}(t_{0})=\bar{c}_{0,m}$ and
${p}_{m}(t_{f})=\bar{c}_{m,m}$.

Additional properties and algorithms for Bernstein polynomials can be found in
[16, 17] and the open-source implementation is available at [26].

##  III Estimator

The objective of the estimator is to find the target’s location based on noisy
measurements and to provide an estimate of their distribution. In this work,
we assume that the beacon provides Received Signal Strength Indicator (RSSI)
measurements, i.e., the relative signal strength measured by the receiver.
Based on these readings, it is possible to obtain range measurements between
the beacon and the receiver using the following relationship:
$r=10^{(P-RSSI)/10n}$, where $P$ and $n$ are constant parameters [27]. We
assume that the data is degraded by noise whose intensity is inversely
proportional to the squared distance between the beacon and the receiver. For
example, this behavior is observed in electromagnetic interference affecting
the data gathered by ARVA sensors, which are commonly used for search and
rescue missions of avalanche victims [24, 28]. This data is then used to
compute an estimate of the target’s position and update the distribution
function at every time step.

###  III-A Target Localization

Given a set of $n$ range measurements, the estimate of the transmitter’s
position is obtained by solving the following unconstrained optimization
problem.

###### Problem 1

Let $r_{k}$ be the range measurement obtained at time $t_{k}:=k\Delta_{s}$,
where $\Delta_{s}$ is the sampling time and $k\in\\{0,1,\ldots n\\}$. Find the
estimate of the transmitter’s location at time $t_{n}$ i.e., $\hat{p}_{n}$,
such that the following is minimized

$J_{n}=\sum_{k=1}^{n}||p(t_{k})-\hat{p}_{n}||-r_{k}.$ (6)

Because of its simplicity, and therefore short computational time, the
localization problem presented above can be solved continuously as new
readings are obtained.

###### Remark 1

We note that because the objective of our framework is not to design an
estimator, but rather to develop a motion planner that can actively improve
its performance, different types of estimator can be used for this step such
as Extended Kalman filter, least squares estimators etc.

###  III-B Smooth Estimation of Density Functions

Let $\\{\hat{p}_{1},\cdots\hat{p}_{n}\\}$ be the set of estimates obtained by
solving Problem 1 at time steps $\\{t_{1},\ldots,t_{n}\\}$. The empirical
distribution function of the target’s location, $F_{n}(\zeta)$,
$\zeta\triangleq[\zeta_{x},\zeta_{y}]\in Z\subseteq\mathbb{R}^{2}$,
$\zeta_{x},\zeta_{y}\in[\zeta_{\min},\zeta_{\max}]$, is given by

$F_{n}(\zeta)=\begin{bmatrix}F_{n}^{(x)}(\zeta_{x})\\\
F_{n}^{(y)}(\zeta_{y})\end{bmatrix}=\begin{bmatrix}\frac{1}{n}\sum_{i=1}^{n}I(\hat{x}_{n}\leq\zeta_{x})\\\\[3.0pt]
\frac{1}{n}\sum_{i=1}^{n}I(\hat{y}_{n}\leq\zeta_{y})\end{bmatrix},$ (7)

where $I(a\leq b)$ yields $1$ when $a\leq b$, and $0$ otherwise. Then, the
$m$-th order Bernstein estimator of the CDF, $F(\zeta)$, is defined as follows
[22, 23]:

$\hat{F}_{m,n}(\zeta)={\small\begin{bmatrix}\sum_{j=0}^{m}F_{n}^{(x)}\left(\frac{j}{m}(\zeta_{\max}-\zeta_{\min})+\zeta_{\min}\right)b_{j,m}(\zeta_{x})\\\
\sum_{j=0}^{m}F_{n}^{(y)}\left(\frac{j}{m}(\zeta_{\max}-\zeta_{\min})+\zeta_{\min}\right)b_{j,m}(\zeta_{y})\end{bmatrix}}.\vspace{1mm}$
(8)

Clearly, the choice of the order of polynomial $m$ is dependent on $n$ and
plays a crucial role in the estimation of $F(\zeta)$. In [22], the authors
show that for

$n^{2/3}\leq m\leq\left(\frac{n}{log(n)}\right)^{2},$

$\hat{F}_{m,n}(\zeta)$ is strongly consistent and as $m,n\to\infty$,

$||\hat{F}_{m,n}-F_{n}||=O\left((n^{-1}\log n)^{1/2}(m^{-1}\log
m))^{1/4}\right).$

In other words, $\hat{F}_{m,n}(\zeta)$ is used to smooth $F_{n}(\zeta)$, and
therefore provides a continuous approximation of $F(\zeta)$. Finally, the
density function, $f(\zeta)$ can be approximated by differentiating (8) using
Property 1, i.e.

$f(\zeta)\approx\hat{f}_{m,n}(\zeta)=\frac{d}{d\zeta_{x}}\hat{F}_{m,n}(\zeta)+\frac{d}{d\zeta_{y}}\hat{F}_{m,n}(\zeta)$
(9)

##  IV Optimal Motion Planner

At discrete times $t_{i}$ with $i\in\\{0,1,\ldots\\}$, the most up-to-date
estimate, $\hat{p}_{i}$, is sent to the planner. Then, the objective of the
motion planning algorithm is to compute desired trajectories, which maximize
the efficiency of the estimator until the target is successfully localized.

Recall that the CRLB defines a lower bound on the variance of an estimator
[29]. Based on the eigenvalues of the CRLB, it is possible to calculate the
length of the axis of the uncertainty ellipse [30]. Furthermore, the
reciprocal of the Cramer Rao Lower Bound, i.e., the Fisher Information Matrix
(FIM), represents the amount of information contained in a set of measurements
[31]. In turn, this can be used as a metric of the performance of the
estimator: maximizing the FIM decreases the area of the error covariance
ellipse, which ensures more accurate estimates [31, 32].

Given a measurement model and assuming that the target is stationary, the
Fisher Information Matrix for $t\in[t_{i},t_{f}]$ can be defined as follows:

$\mathcal{I}(p,\hat{p}_{i},t_{i},t_{f})=\frac{1}{\sigma^{2}}\bigintss_{t_{i}}^{t_{f}}\begin{bmatrix}\frac{\Delta
x(\tau)^{2}}{r(\tau)^{2}}&\frac{\Delta y(\tau)\Delta
x(\tau)}{r(\tau)^{2}}\vspace{3mm}\\\ \frac{\Delta y(\tau)\Delta
x(\tau)}{r(\tau)^{2}}&\frac{\Delta
y(\tau)^{2}}{r(\tau)^{2}}\end{bmatrix}d\tau,$ (10)

where

$\displaystyle\Delta x(t)$ $\displaystyle=x(t)-\hat{x}_{i},$
$\displaystyle\Delta y(t)$ $\displaystyle=y(t)-\hat{y}_{i},$ $\displaystyle
r(t)$ $\displaystyle=\sqrt{\Delta x(t)^{2}+\Delta y(t)^{2}}.$

and $\sigma^{2}$ is the variance of the measurements. It can be noticed that
$\mathcal{I}$ only depends on the geometry of the vehicle’s path. We can
leverage this insight by designing a cost function that maximizes the Fisher
Information Matrix, and by doing so improving the performance of the
estimator.

###  IV-A Motion planner

Recall that the intensity of the measurement’s noise is inversely proportional
to the squared distance between the beacon and the receiver. This motivates
the need to steer the vehicle on a trajectory toward the source of the signal,
which is unknown. To do so, we devise a cost that directs the vehicle to a
particular location in the search area in a manner that is proportional to the
probability that it is the true location of the target. This way, the use of
the density function, which is given by the Bernstein estimator, effectively
acts as a weight. In other words, at the beginning of the mission, when the
location of the transmitter is still quite uncertain, the maximization of
$\mathcal{I}$ is prioritized. Instead, when an estimate is found to have high
probability, the vehicle is steered more aggressively toward that position to
precisely pinpoint the target’s exact location. Lastly, we include a term in
the cost function which minimizes the actuation effort in order to preserve
battery life. With this in mind, the motion planning problem can be formulated
as follows.

###### Problem 2

Let $w_{1},w_{2},w_{3},w_{4}>0$, and
$f(\zeta)=[f_{x}(\zeta)\,f_{y}(\zeta)]^{\top}$ be the probability density
function of the target’s position estimates. At time $t_{i}$ compute a
trajectory $p:[t_{i},\,t_{f}]\to\mathbb{R}^{2}$ that minimizes the following:

$\begin{split}J=&w_{1}(t_{f}-t_{i})+w_{2}\displaystyle{\int_{t_{i}}^{t_{f}}}||\ddot{p}(\tau)||^{2}\,d\tau\\\
+&w_{3}\displaystyle\iint_{Z}[p(t_{f})-\zeta]^{T}\text{diag}\left(f(\zeta)\right)\hskip
1.42262pt[p(t_{f})-\zeta]\hskip 2.84526ptd\zeta\\\\[5.0pt]
-&w_{4}\log(\det\mathcal{I})\end{split}$ (11)

subject to

$\displaystyle p(t_{i})=p_{t_{i}},\quad\dot{p}(t_{i})=\dot{p}_{t_{i}},$ (12)
$\displaystyle h_{c}({p}(t),\dot{p}(t),\ddot{p}(t))\leq 0,\quad\forall
t\in[t_{i},t_{f}].$ (13)

###### Remark 2

In Problem 2, equation (12) ensures continuity with respect to the previous
trajectory at time $t_{i}$, while equation (13) ensures that the trajectory
satisfies the feasibility constraints of the vehicle, e.g., maximum and
minimum velocity, minimum distance with respect to obstacles, etc.

###### Remark 3

We opted to use the determinant of $\mathcal{I}$ as a measure of the
information captured by the estimator. However, alternative objective
functions can be used instead. Some of the most common choices include
maximizing the eigenvalues of $\mathcal{I}$, and maximizing the trace of
$\mathcal{I}$. For a more complete discussion on the criterion for
maximization of the Fisher Information Matrix, we refer the reader to [31].

Due to its complexity, Problem 2 cannot be solved analytically. Here, we use
numerical methods to approximate open-loop optimal solutions to Problem 2. In
particular, we employ the direct method that uses Bernstein polynomial
approximation to transcribe optimal control problems as finite-dimensional
optimization problems [25]. The motivation is that, as previously mentioned,
Bernstein polynomials offer favorable geometric properties that are
particularly useful when dealing with trajectory generation for autonomous
vehicles [14, 15].

Let the trajectory be approximated by a $d$-th order Bernstein polynomial.
Then, Problem 2 can be re-stated as follows.

###### Problem 3

Let $w_{1},w_{2},w_{3},w_{4}>0$ and $\hat{f}_{m,n}(\zeta)$ be the density
estimator shown in (9). At time $t_{i}>0$ compute a vector of Bernstein
polynomial coefficients $\bar{c}_{0,d},\,\cdots,\,\bar{c}_{d,d}$ such that
$p_{d}$ defined in (1), minimizes

$\begin{split}J=w_{1}&(t_{f}-t_{i})+w_{2}\displaystyle{\int_{t_{i}}^{t_{f}}}||\ddot{p}_{d}(\tau)||^{2}\,d\tau\\\
+&w_{3}\displaystyle\iint_{Z}[p_{d}(t_{f})-\zeta]^{T}\text{diag}\left(\hat{f}_{m,n}(\zeta)\right)\hskip
1.42262pt[p_{d}(t_{f})-\zeta]\hskip 2.84526ptd\zeta\\\\[5.0pt]
-&w_{4}\log(\det\mathcal{I})\end{split}$ (14)

subject to

$\displaystyle
p_{d}(t_{i})=p_{t_{i}},\quad\dot{p}_{d}(t_{i})=\dot{p}_{t_{i}},$ (15)
$\displaystyle h_{c}({p}_{d}(t),\dot{p}_{d}(t),\ddot{p}_{d}(t))\leq
0,\,\,\forall t\in[t_{i},\,t_{f}].$ (16)

By virtue of Properties 4-3 presented in Section II-B, the above problem
results into a nonlinear programming problem, which can be solved using off-
the-shelf nonlinear optimization solvers in near real time.

###  IV-B Replanning and Termination Criteria

Because the estimate of the target’s position is likely to change as the
mission unfolds, we employ a replanning strategy to exploit the most up-to-
date information. Specifically, as mentioned previously, the trajectories are
generated at constant time intervals until termination criteria are met.

As termination criteria, we again make use of probability distributions.
Specifically, in [22] the density estimator is proved to be asymptotically
normal as $m,n\to\infty$ and $m/n\to 0$. By virtue of this property, we can
easily obtain the standard deviation of the distribution. We then terminate
the algorithm when it reaches a $95\%$ confidence, or $2\sigma$, that the true
location of the target is within a given radius from the estimated one i.e.,
$2\sigma\leq r_{t}$, or the maximum final time has been reached, i.e. $t\geq
t_{f,max}$

##  V Numerical Results

In this section, simulation results conducted on MATLAB are given to prove the
efficacy of our proposed framework. The target localization problem presented
in Problem 1 is solved using MATLAB $\mathtt{fminunc}$, while the motion
planning problem given in Problem 2 is solved using MATLAB $\mathtt{fmincon}$.
To ensure feasibility of the trajectory, we impose a constraint on the maximum
velocity of the vehicle, i.e., $||\dot{p}(t)||\leq v_{max}$, with
$v_{max}=$1\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$$. New range measurements
become available at a frequency of $2\text{\,}\mathrm{Hz}$, and are degraded
by Gaussian noise with $\sigma=0.1$. The Bernstein estimator of the CDF is
generated with order $m=n^{3/4}+2$. The trajectories are generated with a
replanning interval $\Delta T=$5\text{\,}\mathrm{s}$$, and termination
criteria are set to $r_{t}=$2\text{\,}\mathrm{m}$$ and
$t_{f,max}=$350\text{\,}\mathrm{s}$$.

To better demonstrate the results, we performed two simulations: the first one
uses the motion planning problem given in Problem 3, while the second
simulation does not include the last term in the cost function, which maximize
the Fisher Information Matrix, $\mathcal{I}$. The results of the two
simulations are shown in blue and cyan, respectively.

The execution of both simulations is shown in Figure 2. It can be noticed that
in the first simulation, the vehicle follows a longer path around the target,
and by doing so is able to collect more meaningful measurements, which lead to
better estimates. This can be seen more clearly in Figure 3, which shows the
estimation error over time. It can be noticed that throughout the entire
simulation, this path results in a much smaller estimation error. On the other
hand, in the second simulation, the vehicle quickly moves toward the estimated
location of the target and subsequently circles around a small area. By doing
so, the agent collects range measurements close to each other, which do not
provide enough information to quickly converge toward the true location of the
target.

Finally, Figure 4 shows the CDF and PDF for both $x$ and $y$ coordinates at
the time of the last trajectory replanning. The true location of the target is
represented by the vertical red line. It can be seen that for the second
simulation, the probability curves are slightly skewed, and the density
functions have a lower peak, which indicates less confidence in the estimates.

Figure 2: Path traveled by the vehicle during the mission. Figure 3:
Estimation error over time. Figure 4: Cumulative density function (top row)
for $x$ and $y$ coordinates, respectively, and probability density function
(bottom row) for $x$ and $y$ coordinates, respectively. The true position of
the target is represented by the red vertical line.

##  VI Conclusions

In this article, we propose a target localization framework based on two main
algorithms: an estimator and a motion planner. The objective of the estimator
is to localize the target based on noisy range measurements obtained by
converting a beacon’s RSSI signal. Moreover, Bernstein polynomials are used to
smooth the empirical distribution function of the measurement, and thus
provide a continuous estimate of the density function. The Fisher Information
Matrix, which describes the optimal sensor-target geometry was derived for the
scenario considered. Then, the motion planning problem was formulated as a
multi-objective nonlinear optimal control problem to increase the
effectiveness of the estimator. Finally, simulation results were provided to
demonstrate the efficacy of the proposed algorithm. As part of future work, we
will extend our framework to simultaneously localize multiple targets. Lastly,
experimental results will be conducted to reflect real-life scenarios.

## References

  * [1] A. Albanese, V. Sciancalepore, and X. Costa-Pérez, “Sardo: An automated search-and-rescue drone-based solution for victims localization,” IEEE Transactions on Mobile Computing, vol. 21, no. 9, pp. 3312–3325, 2022.
  * [2] C. Sampedro, A. Rodriguez-Ramos, H. Bavle, A. Carrio, P. de la Puente, and P. Campoy, “A fully-autonomous aerial robot for search and rescue applications in indoor environments using learning-based techniques,” Journal of Intelligent & Robotic Systems, vol. 95, no. 2, pp. 601–627, 2019\.
  * [3] D. J. Pack, P. DeLima, G. J. Toussaint, and G. York, “Cooperative control of uavs for localization of intermittently emitting mobile targets,” IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics), vol. 39, no. 4, pp. 959–970, 2009.
  * [4] X. Liu, J. Zhang, S. Jiang, Y. Yang, K. Li, J. Cao, and J. Liu, “Accurate localization of tagged objects using mobile rfid-augmented robots,” IEEE Transactions on Mobile Computing, vol. 20, no. 4, pp. 1273–1284, 2021.
  * [5] A. Buffi and B. Tellini, “Measuring uhf-rfid tag position via unmanned aerial vehicle in outdoor scenario,” in 2018 IEEE 4th International Forum on Research and Technology for Society and Industry (RTSI), pp. 1–6, 2018.
  * [6] C. Yang, L. Kaplan, and E. Blasch, “Performance measures of covariance and information matrices in resource management for target state estimation,” IEEE Transactions on Aerospace and Electronic Systems, vol. 48, no. 3, pp. 2594–2613, 2012.
  * [7] M. Sadeghi, F. Behnia, and R. Amiri, “Optimal sensor placement for 2-d range-only target localization in constrained sensor geometry,” IEEE Transactions on Signal Processing, vol. 68, pp. 2316–2327, 2020.
  * [8] S. Xu, Y. Ou, and W. Zheng, “Optimal sensor-target geometries for 3-d static target localization using received-signal-strength measurements,” IEEE Signal Processing Letters, vol. 26, no. 7, pp. 966–970, 2019.
  * [9] P. Anjaly and A. Ratnoo, “Observability enhancement of maneuvering target with bearings-only information,” Journal of Guidance, Control, and Dynamics, vol. 41, no. 1, pp. 184–198, 2018.
  * [10] N. Crasta, D. Moreno-Salinas, A. M. Pascoal, and J. Aranda, “Multiple autonomous surface vehicle motion planning for cooperative range-based underwater target localization,” Annual Reviews in Control, vol. 46, pp. 326–342, 2018.
  * [11] S. Uluskan, “Noncausal trajectory optimization for real-time range-only target localization by multiple uavs,” Aerospace Science and Technology, vol. 99, p. 105558, 2020.
  * [12] H. Roh, M.-H. Cho, and M.-J. Tahk, “Trajectory optimization using cramér-rao lower bound for bearings-only target tracking,” in 2018 AIAA Guidance, Navigation, and Control Conference, p. 1591, 2018.
  * [13] N. T. Hung, N. Crasta, D. Moreno-Salinas, A. M. Pascoal, and T. A. Johansen, “Range-based target localization and pursuit with autonomous vehicles: An approach using posterior crlb and model predictive control,” Robotics and Autonomous Systems, vol. 132, p. 103608, 2020.
  * [14] V. Cichella, I. Kaminer, C. Walton, and N. Hovakimyan, “Optimal motion planning for differentially flat systems using bernstein approximation,” IEEE Control Systems Letters, vol. 2, no. 1, pp. 181–186, 2017.
  * [15] V. Cichella, I. Kaminer, C. Walton, N. Hovakimyan, and A. M. Pascoal, “Optimal multi-vehicle motion planning using bernstein approximants,” IEEE Transactions on Automatic Control, 2020.
  * [16] C. Kielas-Jensen and V. Cichella, “Bebot: Bernstein polynomial toolkit for trajectory generation,” in 2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 3288–3293, IEEE, 2019.
  * [17] C. Kielas-Jensen and V. Cichella, “Bernstein polynomial-based transcription method for solving optimal trajectory generation problems,” arXiv preprint arXiv:2010.09992, 2020.
  * [18] C. Kielas-Jensen, V. Cichella, T. Berry, I. Kaminer, C. Walton, and A. Pascoal, “Bernstein polynomial-based method for solving optimal trajectory generation problems,” Sensors, vol. 22, no. 5, p. 1869, 2022.
  * [19] R. T. Farouki, “The bernstein polynomial basis: A centennial retrospective,” Computer Aided Geometric Design, vol. 29, no. 6, pp. 379–419, 2012.
  * [20] R. A. Vitale, “A bernstein polynomial approach to density function estimation,” in Statistical inference and related topics, pp. 87–99, Elsevier, 1975.
  * [21] A. Tenbusch, “Two-dimensional bernstein polynomial density estimators,” Metrika, vol. 41, no. 1, pp. 233–253, 1994.
  * [22] G. J. Babu, A. J. Canty, and Y. P. Chaubey, “Application of bernstein polynomials for smooth estimation of a distribution and density function,” Journal of Statistical Planning and Inference, vol. 105, no. 2, pp. 377–392, 2002.
  * [23] A. Leblanc, “On estimating distribution functions using bernstein polynomials,” Annals of the Institute of Statistical Mathematics, vol. 64, no. 5, pp. 919–943, 2012.
  * [24] C. Tabasso, N. Mimmo, V. Cichella, and L. Marconi, “Optimal motion planning for localization of avalanche victims by multiple uavs,” IEEE Control Systems Letters, vol. 5, no. 6, pp. 2054–2059, 2021.
  * [25] V. Cichella, I. Kaminer, C. Walton, N. Hovakimyan, and A. Pascoal, “Bernstein approximation of optimal control problems,” arXiv preprint arXiv:1812.06132, 2018.
  * [26] “Bernstein/Bézier trajectory toolkit.” https://github.com/caslabuiowa/OptimalBezierTrajectoryGeneration.  Accessed: 2022-09-13.
  * [27] T. S. Rappaport et al., Wireless communications: principles and practice, vol. 2.  prentice hall PTR New Jersey, 1996.
  * [28] N. Mimmo, P. Bernard, and L. Marconi, “Avalanche victim search via robust observers,” IEEE Transactions on Control Systems Technology, pp. 1–12, 2020.
  * [29] R. J. Larsen and M. L. Marx, An introduction to mathematical statistics.  Prentice Hall, 2005.
  * [30] Y. Oshman and P. Davidson, “Optimization of observer trajectories for bearings-only target localization,” IEEE Transactions on Aerospace and Electronic Systems, vol. 35, no. 3, pp. 892–902, 1999.
  * [31] D. Ucinski, Optimal measurement methods for distributed parameter system identification.  CRC press, 2004.
  * [32] S. Ponda, R. Kolacinski, and E. Frazzoli, “Trajectory optimization for target localization using small unmanned aerial vehicles,” in AIAA guidance, navigation, and control conference, p. 6015, 2009.


11institutetext: Dipartimento di Fisica, Sapienza UniversitÃ di Roma, Piazzale
Aldo Moro 5, I-00185 Rome, Italy 22institutetext: Departamento de FÄsica
Teorica & CIAFF, Facultad de Ciencias, Universidad Autonoma de Madrid, Modulo
8, E-28049 Madrid, Spain 33institutetext: Institute for Astronomy, University
of Edinburgh, Edinburgh EH9 3HJ, UK

# Impact of filaments on galaxy cluster properties in The Three Hundred
simulation

S. Santoni 11<EMAIL_ADDRESS>M. De Petris 11 A. Ferragamo 11 G.
Yepes 22 W. Cui 2233

###### Abstract

Galaxy clusters and their filamentary outskirts reveal useful laboratories to
test cosmological models and investigate Universe composition and evolution.
Their environment, in particular the filaments of the Cosmic Web to which they
are connected, plays an important role in shaping the properties of galaxy
clusters. In this project, we analyse the gas filamentary structures present
in 324 regions of The Three Hundred hydrodynamical simulation extracted with
the DisPerSE filament finder. We estimate the number of gas filaments globally
connected to several galaxy clusters, i.e. the connectivity k, with a mass
range of $10^{13}\leq M_{200}\,h^{-1}\,M_{\odot}\leq 10^{15}$ at redshift
$z=0$. We study the positive correlation between the connectivity and mass of
galaxy clusters. Moreover, we explore the impact of filaments on the dynamical
state of clusters, quantified by the degree of relaxation parameter $\chi$.

##  1 Introduction

Clusters of galaxies, the largest gravitationally bound systems in the
Universe, reside at the nodes of the Cosmic Web Bond:1996  and are connected
by a multitude of filamentary structures. In the outskirts of the clusters,
matter and galaxies are funneled towards the centre through filaments. A
comprehensive knowledge of the filaments connected to galaxy clusters is
essential to understand the influence of the environment on galaxy cluster
properties and evolution. One way to quantify the filamentary skeleton around
the cluster is through the so-called connectivity $k$ Codis:2018 , defined as
the number of filaments globally connected to the galaxy cluster, estimated at
a specific aperture. This project aims to investigate the impact of filaments
on the main galaxy cluster properties, in particular in this work we focus on
their masses and dynamical state.

##  2 The Three Hundred project

In this work, we analyse the multiple zoom-in regions of The Three Hundred
hydrodynamical simulation Cui:2018 . The Three Hundred
project111https://the300-project.org aims to model 324 regions whose volumes
have cubic side lengths of $30\,h^{-1}\,Mpc$, centered on massive galaxy
clusters, with a mass $M_{200}>6.42\times 10^{14}\,h^{-1}\,M_{\odot}$. The
Three Hundred regions are re-simulated with higher resolution from the 324
most massive galaxy clusters at $z=0$ of the $1\,h^{-1}\,Gpc$ Dark Matter-only
MDPL2 MultiDark Klypin:2016  simulation. The cosmological parameters used in
the MDPL2 and The Three Hundred simulations are those measured by the Planck
mission Planck Collaboration:2016 . The regions have been re-simulated with
different baryonic models: Gadget-MUSIC Sembolini:2013 , Gadget-X Rasia:2015 ,
which is used in this work, and more recently Gizmo-Simba Cui:2022 . For each
region, 128 snapshots are available from redshift $z=17$ to $z=0$. The
simulated regions were analysed using the AHF halo finder Knollmann:2011
which self-consistently includes both gas and stars in the halo finding
process. The halo finder extracts haloes and estimates their properties, such
as the radius $R_{\Delta}$222The subscript $\Delta$ indicates the overdensity,
i.e. value of the ratio between the density of the cluster at that radius and
the critical density of the Universe $\rho_{c}=3H^{2}/(8\pi G)$ at the
cluster’s redshift. , mass $M_{\Delta}$ and density profile.

##  3 Methods

###  3.1 Cosmic Web extraction

The gas particle distribution of The Three Hundred regions are analysed with
DisPerSE Sousbie:2011 , a topological structure finder, designed to extract
the structures of the Cosmic Web. The finder identifies topologically
significant features in the input density field, which is obtained through a
Delaunay tesselation in the case of a 2D or 3D discrete distribution. The
noise introduced by the finite sampling of the distribution is quantified and
reduced with the persistence and topological simplification theories. The
persistence parameter, which quantifies the robustness of a topological pair,
is defined as the difference of the values of the critical points in the pair
and is used to filter low significant filaments. As an output, DisPerSE
provides the positions of the extreme points found in the distribution:
maxima, minima, saddle points and bifurcation points, where a filament splits
in two. The filaments are given as a set of segments connecting a maximum and
a saddle point.
For our analysis, we first binned the gas particle distribution of each region
at redshift $z=0$ in a three-dimensional grid of $30\,h^{-1}\,Mpc$ per side,
and within each region we give a pixel resolution of $150\,h^{-1}\,kpc$. Then,
to avoid sharp variations from one pixel to another, we applied a Gaussian
smoothing with a $\sigma$ of 4 pixels. Finally we applied an absolute
persistence cut of $0.2$, to focus on significant filaments connecting
clusters and haloes. The node distribution extracted from DisPerSE was
compared with the AHF halo catalogue of each region, to match the DisPerSE
maximum points to the simulated haloes and clusters. To avoid a possible low-
resolution contamination near the borders, we consider only the haloes inside
a sphere of $13\,h^{-1}\,Mpc$ radius from the centre of each region. The final
data set includes $3\times 10^{3}$ haloes and clusters with a mass range from
$10^{13}\leq M_{200}\,h^{-1}\,M_{\odot}\leq 5\times 10^{15}$.
The gas skeleton extracted from The Three Hundred simulated regions is a good
tracer of the overall matter distribution and accretion to galaxy clusters. In
particular, throughout the 324 simulated regions there is a good spatial
agreement between the gas filaments and both the Dark Matter and mock galaxy
filaments, both from 3D and 2D extractions Kuchner:2020 ; Kuchner:2021 .

###  3.2 Connectivity measurements

For each halo we estimated the connectivity $k$, which is defined in
Codis:2018  as the number of filaments globally connected to a cluster.
Different definitions are used in literature to estimate this parameter. In
this work, we compute the connectivity as the number of filaments crossing a
specific spherical surface at a radius $R_{\Delta}$ from the centre of the
halo. In particular, we estimate the connectivity at $R_{200}$ and at
$R_{500}$, defined respectively as $k_{200}$ and $k_{500}$. With this
definition, we are taking into consideration also the filaments that are
coming from substructures and bifurcation points that lie within the sphere,
which contribute to the cluster’s properties and therefore also to its
connectivity.

##  4 Results

In this work we investigate if the number of filaments connected to a cluster,
quantified by the connectivity, is correlated with the main properties of the
cluster itself, mainly its mass and dynamical state.

###  4.1 Connectivity and galaxy cluster mass

The number of filaments connected to a cluster is expected to correlate with
the mass of the cluster itself, as many studies show Aragon-Calvo:2010 ;
Sarron:2019 ; Darragh Ford:2019 ; Malavasi:2020 . The expected trend is for
the connectivity to increase with the mass of the cluster.
Taking the advantage to investigate a cluster’s sample with the largest mass
range, we analyse the connectivity of a set of haloes and clusters extracted
from The Three Hundred simulation at redshift $z=0$. Figure 1 shows the values
of the connectivity $k_{200}$ as a function of the clusters mass $M_{200}$. We
measured the mean and the standard deviation for each mass bin, chosen by
taking into account the overall mass distribution in The Three Hundred
simulation at $z=0$, shown in the bottom panel of the figure. We performed a
linear fitting, whose parameters are shown in Table 1, as the following: $\log
k_{200}=A\cdot\log M_{200}+B$.

Table 1: The fitted parameters for the $\log k_{200}-\log M_{200}$ relation | $\log k_{200}-\log M_{200}$
---|---
A | 0.298 $\pm$ 0.016
B | -3.78 $\pm$ 0.23

We compare our estimates of the connectivity to similar measurements from
literature studies, as shown in Figure 2. We compare our sample both to
simulated, as Aragon-Calvo:2010 , and observed data, such as Sarron:2019 ;
Darragh Ford:2019 ; Malavasi:2020 . The trend of the connectivity with the
mass is well in agreement, within the errors, with the ones from literature,
despite the differences presented in Table 2, and over a larger mass range.
More specifically, the connectivity is estimated at a fixed aperture in
Aragon-Calvo:2010 ; Sarron:2019 , while the connectivity from our work and
Darragh Ford:2019 ; Malavasi:2020  is estimated at an overdensity radius.
Other differences can arise from the Cosmic Web component chosen to extract
the filamentary skeletons, such as gas, Dark Matter or galaxies particles and
from the different filament finder used in the analyses. We refer to Aragon-
Calvo:2010 ; Sarron:2019 ; Darragh Ford:2019 ; Malavasi:2020  for a more
detailed description of the datasets used.

Figure 1: The connectivity of haloes and clusters $\log k_{200}$ plotted as a function of the mass $\log\,(M_{200}/h^{-1}M_{\odot})$ (grey points). In the top panel the mean and standard deviation values are also plotted. The bottom panel shows the mass distribution of haloes and clusters of The Three Hundred hydrodynamical simulation at redshift $z=0$ analysed in this work. Figure 2: Connectivity of The Three Hundred clusters, compared with literature values. For each work, we plot the mean values of the connectivity. The error bars refer to the standard deviation values for this work, Aragon-Calvo:2010 ; Sarron:2019 ; Malavasi:2020 , while for Darragh Ford:2019  they represent the errors on the mean values.  Table 2: A summary of the parameters of previous studies compared to this work. | Data | $k_{R}$ | M | CW extraction
---|---|---|---|---
The300 | Hydro simulation | $R_{200}$ | $M_{200}$ | 3D gas particles
AC+10 Aragon-Calvo:2010  | DM simulation | 3 Mpc | $M_{Vir}$ | 3D DM particles
S+19 Sarron:2019  | Observations | 1.5 Mpc | $M_{200}$ | 2D galaxies
DF+19 Darragh Ford:2019  | Observations | 1.5 $R_{Vir}$ | $M_{200}$ | 2D galaxies
M+20 Malavasi:2020  | Observations (Coma Cl.) | $R_{Vir}$ | $M_{200}$ | 3D z-slice galaxies

###  4.2 Connectivity and galaxy cluster dynamical state

In this subsection we investigate the correlation between the connectivity and
the dynamical state of The Three Hundred clusters. We quantify the dynamical
state with the degree of relaxation $\chi$, as defined by Haggar:2020 :

$\chi_{\Delta}=\left(\frac{(\frac{f_{s}}{0.1})^{2}+(\frac{\Delta_{r}}{0.04})^{2}+(\frac{|1-\eta|}{0.15})^{2}}{3}\right)^{-1/2}$

where $f_{s}$ is the sub-halo mass fraction, $\Delta_{r}$ is the centre-of-
mass offset and $\eta$ is the virial ratio. The threshold values for these
parameters were chosen following Cui:2017 . A cluster is considered as
dynamically relaxed when $\chi\geq 1$. To study the effect of connectivity on
the dynamical state of clusters, independently on their masses, we divide the
data set in three different connectivity and mass sub-samples. Respectively,
we consider weakly connected ($k_{200}<4$), medium connected ($k_{200}=4$) and
highly connected ($k_{200}>4$) clusters along with low mass ($M_{200}<7\times
10^{13}\,h^{-1}M_{\odot}$), medium ($7\times 10^{13}\leq
M_{200}\,h^{-1}M_{\odot}<5.5\times 10^{14}$) and massive ($M_{200}\geq
5.5\times 10^{14}\,h^{-1}M_{\odot}$) clusters.
In the left and right panels of Figure 3 we display the degree of relaxation
as a function of the mass and connectivity, respectively, for three
connectivity and mass bins. At a fixed mass, the left panel of Figure 3 shows
that there is no evident correlation between connectivity and the degree of
relaxation, as the three sub-samples overlap. On the other hand, at fixed
connectivity, the right panel of the figure shows a slight correlation between
the mass and $\chi$, indicating that less massive haloes are on average more
dynamically relaxed. This result is in disagreement with that found by
Gouin:2021 , where at fixed mass, weakly connected clusters are on average
more relaxed than highly connected clusters. These differences may depend on
the different data set and Cosmic Web extraction analysed in their work, to
which we refer the reader for more details.

Figure 3: Left panel: mean degree of relaxation as a function of the mass
$M_{200}$ for three connectivity $k_{200}$ bins ($k<4$, $k=4$ and $k>4$).
Right panel: mean degree of relaxation as a function of the connectivity
$k_{200}$ for three mass $M_{200}$ bins. In both panels, the errors bars
represent the errors on the mean values.

##  5 Conclusions

In this work we analysed the gas filamentary structures connected to The Three
Hundred hydrodynamical simulation clusters and their impact on galaxy cluster
properties. We extracted the gas skeletons at $z=0$ with the DisPerSE filament
finder in the 324 regions of the simulation and we estimated the connectivity
$k_{200}$ of haloes and clusters.
The main conclusions of this work can be summarized in the following:

  1. 1.

The connectivity is correlated with the mass of haloes and clusters, with more
massive clusters being on average more connected. This result is compatible
with previous results from literature, both from simulations and observations;

  2. 2.

We do not find a correlation between the connectivity and the dynamical state
of clusters, quantified in terms of the degree of relaxation $\chi$.

## References

  * (1) Bond J. R. et al., Nature 380, 603 (1996)
  * (2) Codis S. et al., Monthly Notices of the Royal Astronomical Society 479, 973 (2018)
  * (3) Cui W. et al., Monthly Notices of the Royal Astronomical Society 480, 2898 (2018)
  * (4) Klypin A. et al., Monthly Notices of the Royal Astronomical Society 457, 4340 (2016)
  * (5) Planck Collaboration et al., Astronomy and Astrophysics 594, A13 (2016)
  * (6) Sembolini F. et al., Monthly Notices of the Royal Astronomical Society 429, 323 (2013)
  * (7) Rasia E. et al., The Astrophysical Journal 813, L17 (2015)
  * (8) Cui W. et al., Monthly Notices of the Royal Astronomical Society 514, 977 (2022)
  * (9) Knollmann S. R. et al., Astrophysics Source Code Library (2011)
  * (10) Sousbie T., Monthly Notices of the Royal Astronomical Society 414, 350 (2011)
  * (11) Kuchner U. et al., Monthly Notices of the Royal Astronomical Society 494, 5473 (2020)
  * (12) Kuchner U. et al., Monthly Notices of the Royal Astronomical Society 503, 2065 (2021)
  * (13) AragÃn-Calvo M. A. et al., Monthly Notices of the Royal Astronomical Society 408, 2163 (2010)
  * (14) Sarron F. et al., Astronomy and Astrophysics 632, A49 (2019)
  * (15) Darragh Ford E. et al., Monthly Notices of the Royal Astronomical Society 489, 5695 (2019)
  * (16) Malavasi N. et al., Astronomy and Astrophysics 634, A30 (2020)
  * (17) Haggar R. et al., Monthly Notices of the Royal Astronomical Society 492, 6074 (2020)
  * (18) Cui W. et al., Monthly Notices of the Royal Astronomical Society 464, 2502 (2017)
  * (19) Gouin C. et al., Astronomy and Astrophysics 651, A56 (2021)


# QSpeech: Low-Qubit Quantum Speech Application Toolkit

1Corresponding author: Jianzong Wang<EMAIL_ADDRESS>Zhenhou Hong, Jianzong
Wang1, Xiaoyang Qu, Chendong Zhao, Wei Tao and Jing Xiao  Ping An Technology
(Shenzhen) Co., Ltd., Shenzhen, China
Emails:
<EMAIL_ADDRESS>
###### Abstract

Quantum devices with low qubits are common in the Noisy Intermediate-Scale
Quantum (NISQ) era. However, Quantum Neural Network (QNN) running on low-qubit
quantum devices would be difficult since it is based on Variational Quantum
Circuit (VQC), which requires many qubits. Therefore, it is critical to make
QNN with VQC run on low-qubit quantum devices. In this study, we propose a
novel VQC called the low-qubit VQC. VQC requires numerous qubits based on the
input dimension; however, the low-qubit VQC with linear transformation can
liberate this condition. Thus, it allows the QNN to run on low-qubit quantum
devices for speech applications. Furthermore, as compared to the VQC, our
proposed low-qubit VQC can stabilize the training process more. Based on the
low-qubit VQC, we implement QSpeech111Our implementation of QSpeech is
publicly available at https://github.com/zhenhouhong/QSpeech, a library for
quick prototyping of hybrid quantum-classical neural networks in the speech
field. It has numerous quantum neural layers and QNN models for speech
applications. Experiments on Speech Command Recognition and Text-to-Speech
show that our proposed low-qubit VQC outperforms VQC and is more stable.

###### Index Terms:

Quantum Neural Networks Library, Hybrid Quantum-classical Neural Networks,
Low-qubit, Variational Quantum Circuit, Speech Application

##  I Introduction

Recent progress in the development and commercialization of quantum
computation has significantly impacted the landscape of Quantum Neural Network
(QNN) [1, 2]. These quantum algorithms designed for the Noisy Intermediate-
Scale Quantum (NISQ) [3] computer are programmed with a Variational Quantum
Circuit (VQC) [4] and a classical optimizer. The approach is similar to the
parameter updating in many machine learning algorithms. The QNN has been
developed by some researchers in various applications [5, 6].

The number of qubits in the NISQ era [3] is $5$ to $50$ qubits. Several
quantum devices have less than 10 qubits. The application of quantum
algorithms in these quantum devices is a key challenge. Mainstream research
has focused on running the quantum algorithm in low-qubit ($2$ to $10$)
quantum devices [7, 8]. Using a scalable trapped atomic ion system [9], the
Grover Search algorithm [7] can be implemented on $3$ qubits quantum device.
The Quantum Principal Component Analysis (QPCA) [10] is developed to the
Resonant QPCA [11], which can be conducted in a low-qubit quantum device using
a resonant analysis algorithm.

However, it can be challenging when QNN or hybrid quantum-classical neural
networks are implemented in speech applications [12, 13]. First, these
algorithms with VQC in speech applications need many qubits, making them
difficult to implement. Second, with numerous qubits, QNN will have the Barren
Plateau problem [14, 15] in the training process. When restricting the number
of qubits, it is challenging to assure model performance.

(a) The framework of speech application with DNN.

(b) The framework of speech application with QNN.

Figure 1: The framework of speech application using DNN compared with QNN. The
blue boxes in (a) represents the DNN modules, the orange boxes in (b)
represents the QNN modules. The only difference is that DNN is replaced with
QNN, with the same preprocessing and postprocessing.

In this study, we propose a novel VQC called low-qubit VQC that allows the QNN
or hybrid quantum-classical neural networks to run in quantum devices with
fewer qubits. Our proposed low-qubit VQC can be used in the neural network
when high-dimensional units, such as 128 and 256, are required in the low-
qubit quantum devices. Furthermore, to stabilize the convergence, we add an
output clip operation in the first linear transformation in low-qubit VQC.
Adding output clip operation can avoid gradient vanishing problems. Compared
to VQC, the low-qubit VQC can enhance the convergence in training.
Additionally, we propose a framework called QSpeech, which uses the QNN module
for easy use in speech applications [16, 17]. QSpeech based on the PennyLane
[18] and PyTorch [19] as backend, can define a quantum circuit in the function
with a simple Python-based user interface. We implement the neural layer
functions in QSpeech such as Quantum Convolution (QConv), quantum long short
term memory [20], Quantum Gate Recurrent Unit (QGRU), and Quantum Attention
(QAttention). Thus, the QSpeech can easily implement the hybrid quantum-
classical neural network.

Our main contributions are as follows:

  * •

Low-qubit VQC has been proposed, allowing QNN to be used in a low-qubit
quantum device. It can stabilize the training process compared to VQC.

  * •

Based on low-qubit VQC, we design a library for the rapid prototyping of
hybrid quantum-classical neural networks in speech applications, named
QSpeech.

  * •

In the experiments, hybrid quantum-classical neural networks were successfully
applied in speech applications.

##  II Related Works

###  II-A Quantum Convolutional Neural Networks

Cong et al. [21] propose a quantum CNN based on a quantum circuit model that
incorporates the ideas of convolutional and pooling layers from classical
CNNs. The proposed architecture is similarly layered; however, it applies 1D
convolutions to the input quantum state (contrary to 2D/3D convolutions on
images). On the input state density $\rho-{in}$, the convolutional layer is
modeled as a quasi-local unitary operation. This unitary operator, denoted by
$U_{i}$, is applied to several successive sets of input qubits up to a
predefined depth. The pooling layer is implemented by measuring some qubits
and applying unitary rotations $V_{i}$ to the nearby qubits. The observations
on the qubits determine the rotation operation. This combines the
functionalities of dimensionality reduction (the output of $V_{i}$ is of lower
dimension) and nonlinearity (due to the partial measurement of qubits). After
the required number of convolutional blocks and pooling unitaries, the unitary
$F$ implements the fully connected layer. A final measurement on the output of
$F$ yields the network output.

Recently, Kerenidis et al. [22] identified the relationship between
convolutions and matrix multiplications and proposed the first quantum
algorithm to compute a CNN’s forward pass as a convolutional product. They
also proposed a quantum backpropagation algorithm to learn network parameters
through gradient descent. In an application of QCNNs, Zhang et al. [23] and
Melnikov et al. [24] propose special convolutional neural networks for
extracting features from graphs to identify those showing the quantum
advantage. Hong et al. [25] propose a novel structure of QCNN that outperforms
other DNN algorithms in predicting protein distance and contact.

###  II-B Quantum Recurrent Neural Networks

There have also been several interesting suggestions for developing quantum
variants of RNN. Hibat-Allah et al. [26] propose a quantum variant of RNNs
using variational wave functions to learn the approximate ground state of a
quantum Hamiltonian. Roth [27] propose an iterative retraining approach using
RNNs for simulating bulk quantum systems via mapping translations of lattice
vectors to the RNN time index. Hopfield Networks [28] were a popular early
type of a recurrent NN for which quantum variants have been developed in
several works [29].

###  II-C Quantum Learning for Speech Processing

Li et al. [30] propose a QNN that combines fuzzy theoretic principles with
continuous digits inference. Their results show that more than $15\%$ error
reduction is obtained on a speaker-independent continuous digits recognition
task. Li et al. [31] propose a QNN used in speech enhancement. The QNN can be
trained on the backpropagation algorithm. Its simulation results show that the
QNN outperforms the traditional spectral subtraction in speech enhancement. Fu
et al. [32] propose an improved particle swarm optimization (IPSO), called
IPSO-QNN, for training the QNN. Its experimental results in speech recognition
are better than other approaches. Yang et al. show that the QCNN can use in
privacy-preserving learning [6]. It uses the QCNN to encode the speech feature
to quantum state for data protection.

##  III Preliminaries

The entire quantum variational circuit [4] consists of three components: Data
Encoding Layer, Variational Layer, and Quantum Measurement Layer. Notably, the
number of qubits and measurements can be changed to fit the problem of
interest, and the variational layer can contain many dashed boxes to increase
the number of parameters, all subject to the capacity and capability of the
quantum machines used in the experiments.

Data Encoding Layer. We first encode a classical input vector into a quantum
state in this layer, which is required for further processing. A general
$N$-qubit quantum state can be represented as the following:

$\left|\phi\right\rangle=\sum_{(q_{1},q_{2},\ldots,q_{N})\in\left\\{0,1\right\\}^{N}}c_{q_{1},\ldots,q_{N}}\left|q_{1}\right\rangle\otimes\left|q_{2}\right\rangle\otimes\left|q_{3}\right\rangle\otimes\ldots\otimes\left|q_{N}\right\rangle$
(1)

where $c_{q_{1},\ldots,q_{N}}\in\mathbb{C}$ is the amplitude of each quantum
state and $q_{i}\in\left\\{0,1\right\\}$. The square of the amplitude
$c_{q_{1},\ldots,q_{N}}$ is the probability of measurement with the post-
measurement state in
$\left|q_{1}\right\rangle\otimes\left|q_{2}\right\rangle\otimes\left|q_{3}\right\rangle\otimes\ldots\otimes\left|q_{N}\right\rangle$,
and the total probability should sum to 1, i.e.,

$\sum_{(q_{1},q_{2},\ldots,q_{N})\in\left\\{0,1\right\\}^{N}}\left\|c_{q_{1},\ldots,q_{N}}\right\|^{2}=1$
(2)

For the quantum gates, the input will first be flattened and transformed into
rotation angles. In this study, $f_{e}(\cdot)$ represent the data encoder.

Variational Layer. After encoding the classical values into a quantum state,
it will be subject to a series of unitary transformations. The variational
layer consists of two parts, one for entanglement and the other for rotation,
which include numerous single-qubit unitary rotations $U$ parameterized by $1$
to $m$ parameters $\theta_{1},...,\theta_{m}$. The parameters labeled by
$\theta_{1},...,\theta_{m}$ are the ones that will be updated by the
optimization procedure. In this study,
$f_{u}(\theta_{1},...,\theta_{m};\cdot)$ represent the variational circuit.

Quantum Measurement(Decode) Layer. We performed a quantum measurement layer to
obtain the transformed data from VCs. The VCs’ ensemble samplings (expectation
values) were considered. Quantum simulation software (for example, PennyLane
or IBM Qiskit) can be used to calculate this value deterministically. While
implementing on a real quantum computer, it is necessary to prepare the same
system and repeatedly use the measurements to gain adequate statistics. In
this study, $f_{d}(\cdot)$ represent the measurement layer.

##  IV Proposed Methods

###  IV-A Challenges of VQC in Low-qubit Quantum Device

The VQC consists of data encoder, variational circuit, and measurement [33,
4]. It only processes the 1D array and necessitates qubits based on the size
of one dimension input. It is difficult to use the VQC in a high-dimensional
input $\mathbf{x}\in\mathbb{R}^{1\times n}$ ($n$ is 128 or 256 or even larger)
because it requires $n$ qubits to encode and transform the input in VQC, which
exceeds the number of qubits. Furthermore, with multiple qubits, it will be
easy to cause the Barren Plateaus problem [14, 15] during the training
process.

###  IV-B Problem Statement

Similar to deep learning [34, 35, 36] in speech application, using the QNN
$f_{Q}(\cdot)$ to model a function can generate $\mathbf{y}$, which represents
the classify of audio or generate speech from data input $\mathbf{x}$.

$\displaystyle\mathbf{y}=f_{q}(\mathbf{W},f_{v}(\theta),\mathbf{x})$ (3)
$\displaystyle\mathcal{L}(\mathbf{W},f_{v}(\theta),\mathbf{x},\mathbf{y}_{true})=\mathcal{L}(f_{q}(\mathbf{W},f_{v}(\theta),\mathbf{x}),\mathbf{y}_{true})$
(4)

In VQC $f_{v}(\theta)$, we use the RY gate $R_{y}(\cdot)$ in Data Encoding
Layer. $\mathbf{y}_{true}$ is the true label in data.

The optimization for QNN in speech application, is similar with DNN:

$\displaystyle\hat{\theta}=\mathop{\arg\min}_{\theta}\mathcal{L}(f_{q}(\mathbf{W},f_{V}(\theta),\mathbf{x}),\mathbf{y}_{true})$
(5)
$\displaystyle\hat{\mathbf{W}}=\mathop{\arg\min}_{\mathbf{W}}\mathcal{L}(f_{q}(W,f_{v}(\theta),\mathbf{x}),\mathbf{y}_{true})$

where the $\hat{\theta}$ and $\hat{\mathbf{W}}$ are optimized parameters in
$f_{q}(\mathbf{W},f_{v}(\theta),\mathbf{x})$. Based on PennyLane setting, VQC
$f_{v}(\theta)$ can perform the optimization using gradient descent
approaches. In this study, we follow the [37], using the Adam optimization
method [38] to minimize the loss
$\mathcal{L}(f_{q}(\mathbf{W},f_{v}(\theta),\mathbf{x}),and\mathbf{y}_{true})$.

###  IV-C Low-qubit Variational Quantum Circuit(VQC)

As shown in Fig. 2, the two purple blocks are the new parts that we added in
VQC. In low-qubit VQC, we add two linear transformations unlike in previous
VQC [21, 20],.

Figure 2: The proposed low-qubit VQC. The “LT” is the linear transformation,
“RY” is the rotation gates along the y axis, $U(\theta_{1},..,\theta_{m})$ is
the single-qubit unitary rotations $U$ parameterized by $1$ to $m$ parameters
$\theta_{1},...,\theta_{m}$.

Consider $\mathbf{x}$ as input, let
$\mathbf{x}=\\{x_{1},x_{2},\dots,x_{n}\\}$, $\mathbf{x}\in\mathbb{R}^{1\times
n}$. LT $f_{lt}(\mathbf{W};\cdot)$ is linear transformation, $\mathbf{W}$ is
the learnable parameter in $f_{lt}(\mathbf{W};\cdot)$.
$f_{lt}(\mathbf{W};\cdot):\mathbf{x}\rightarrow\mathbf{y}$. $\mathbf{q}_{1}$,
$\mathbf{q}_{2}$, and $\mathbf{q}_{3}$ are representative quantum states.
$\mathbf{y}_{1}$ and $\mathbf{y}_{2}$ are the outputs of first and second
$f_{lt}$. $\mathbf{y}=\\{y_{1},y_{2},\dots,y_{j}\\}$. Then the formula of low-
qubit VQC follows:

$\displaystyle\mathbf{y}_{1}=f^{1}_{lt}(\mathbf{W};\mathbf{x})$ (6)
$\displaystyle\mathbf{q}_{1}=f_{e}(\mathbf{y}_{1})$
$\displaystyle\mathbf{q}_{2}=f_{u}(\theta_{1},...,\theta_{m};c(\mathbf{q}_{1}))$
$\displaystyle\mathbf{q}_{3}=f_{d}(\mathbf{q}_{2})$
$\displaystyle\mathbf{y}_{2}=f^{2}_{lt}(\mathbf{W};\mathbf{q}_{3})$

Where $\theta_{1},...,\theta_{m}$ is the parameters in $U(\theta;\cdot)$, $m$
is a positive integer. $f_{e}$, $f_{u}$, and $f_{d}$ are the data encoder,
variational circuit, and measurement operation. $c(\cdot)$ is the CNOT gate.
$f^{1}_{lt}$ and $f^{2}_{lt}$ is the first and second linear transformation.

The advantage of adding the linear transformation is that it can squeeze high-
dimensional input in the data encoding layer or expand the output filters
after measurement. Several qubits are not required since they will be squeezed
with the linear transformation before the data encoding and variational
layers. It can make convergence stable in the training and need fewer
resources to simulate the quantum computation. Additionally, it makes low-
qubit VQC process multiple channels simultaneously while extracting the
features from input. It is beneficial for the QConv and QAttention layers. The
QAttention layer can allow high-dimensional units such as 128 or 256 with few
qubits. Thus, the Transformer with QAttention can be used in the low-qubit
quantum devices for speech applications [39] such as Text-to-Speech (TTS).

###  IV-D Analysis of the Low-qubit VQC

We found the $f_{e}(\mathbf{y}_{1})$ contains the $\arctan(\cdot)$, its
derivative may cause the gradient vanish.
$f_{e}(\mathbf{y}_{1})=R_{y}(\arctan(y_{11}))\left|0\right\rangle...R_{y}(\arctan(y_{1i}))\left|0\right\rangle$.
Its derivative results of $R_{y}(\arctan(y_{1i}))$ as follow:

$\displaystyle\frac{\partial R_{y}(\arctan(y_{1i}))}{\partial
w_{i}}=-R_{y}(\arctan(y_{1i}))^{\intercal}\frac{\partial\arctan(y_{1i})}{\partial
w_{i}}$ (7) $\displaystyle\frac{\partial\arctan(y_{1i})}{\partial
w_{i}}=\frac{1}{1+(y_{1i})^{2}}\frac{\partial y_{1i}}{\partial
w_{i}}=\frac{1}{1+(y_{1i})^{2}}x_{i}$

Where $y_{1i}$ is the element in $\mathbf{y}_{1}$.
$y_{1i}=\mathbf{w}\mathbf{x}$. $\mathbf{w}$ is a vector in $\mathbf{W}$.

We analysis the $x_{i}/(1+(y_{i})^{2})$. To observe the output, we fix the
$y_{i}$ and change the $x_{i}$. We can see that the $y_{i}$ cannot be too big
or too small, or the gradient will vanish. Therefore, we add the output clip
operation in the first linear transformation. It is possible to optimize the
weight of linear transformation. The output clip’s range is $[-3,3]$.

The LT with learnable parameters can be implemented using FC. We compare the
linear transformation with other subsampling approaches like bilinear
interpolation. The linear transformation has two advantages: first, it can
reduce the dimension and extract the features efficiently with learnable
parameters; second, it can extract the features adapt to the VQC with
learnable parameters. Bilinear interpolation simply reduces the dimension,
causing the gradient vanish based on the derivative result.

###  IV-E Framework of QSpeech

Fig. 3 shows the software architecture of QSpeech. In the QSpeech, the main
parts of QNN training and inference are written in python, which calls
PennyLane and PyTorch by switching the backend option. PennyLane allows
differentiable programming of quantum computers. We also provide complete
recipes for conducting speech application experiments. As shown in Fig. 3, the
QSpeech is based on the PyTorch and PennyLane. The functional or layers
backend are qcircuit, qconv, qlstm, qgru, and qattention. The functional
backend builds up the model libraries: Quantum Transformer-TTS (QTransformer-
TTS), Quantum M5 (QM5), and Quantum Tacotron (QTacotron). The Quantum Layers
(QLayers) are the functional backend, and the Quantum Models (QModels) are the
model libraries in the major libraries in QSpeech. Using the PennyLane as a
Quantum Machine Learning (QML) backend can easily embed the code into Python
and PyTorch.

Figure 3: The software architecture of QSpeech. All functional backends
consists of low-qubit VQC.

##  V Implemented Quantum Neural Network in Speech Application

###  V-A Quantum M5 For SCR

(a) M5 model.

(b) QM5 model.

Figure 4: The model architecture of QM5 compared with M5.

Figure 5: The process of QConv1d operation. The $a_{1},a_{2},...,a_{n}$ means
the value of input data. The $f_{1},f_{2},...,f_{n}$ means the feature maps.

The M5 [40] model is a deep convolutional neural network (CNN) that directly
uses time-domain waveforms as inputs. To extract the feature from the raw
speech data, the M5 model in speech command uses one dimension convolutional
layer (Conv1d). Thus, we first implement the quantum one dimension
convolutional (QConv1d) layer that can use in speech commands and other speech
applications [41, 42]. Then, we use the M5 model with the QConv1d layer in the
speech command task, called quantum M5(QM5).

Fig. 5 shows the QConv1d extract feature from raw speech data or feature map.
We propose a QConv1d that uses the low-qubit VQC and help extract from a
feature map, while avoiding source in the quantum simulator.

###  V-B Quantum Tacotron For TTS

Tacotron [43] is an end-to-end generative TTS model that synthesizes speech
directly from characters. It has an Attention Decoder part, which includes a
pre-net, an attention recurrent neural networks (RNN), and a decoder RNN. The
RNN component uses the GRU. As shown in Fig. 6, we change the original GRU to
QGRU in attention GRU and decoder GRU. The others are maintained.

(a) Attention decoder part in Tacotron.

(b) Attention decoder part in QTacotron.

Figure 6: The attention decoder with GRU compared with the attention decoder
with quantum GRU.

Figure 7: The proposed QGRU architecture. Each low-qubit VQC as detailed in
Figure 2. The sigmoid and tanh blocks represent the sigmoid and the hyperbolic
tangent activation function, respectively. $X_{t}$ is the input at time $t$,
$H_{t}$ is for the hidden state, $Y_{t}$ is the output. “DotMul” represents
element-wise multiplication.

To construct the basic unit of the proposed QGRU architecture, a QGRU cell, we
stack the low-qubit VQC blocks as mentioned earlier together. In Fig. 7, each
of the low-qubit VQC is described in the previous section (see also Fig. 2).
In a QGRU cell, there are five low-qubit VQC. For three low-qubit VQC in Fig.
7 below, the input is the concatenation $v_{t}$ of the hidden state $H_{t-1}$
from the previous time step and the current input vector $X_{t}$, and the
output is three vectors obtained from the measurements at the end of each low-
qubit VQC. The output value can then be further processed with two low-qubit
VQC in Fig. 7 top right, to get $Y_{t}$ and $H_{t}$.

###  V-C Quantum Transformer For TTS

Transformer-TTS [44] is an architecture that uses two parts (Encoder and
Decoder) to transform one sequence into another. Here we implement the quantum
self-attention with a single head and multihead. As shown in Fig. 8(a), adding
two VQC blocks in the input and output, others are the same with Scaled-Dot
Product Attention [44]. The quantum self-attention with multihead as shown in
Fig. 8(b), the $Q$, $V$, $K$ computed by the VQC, input to the quantum self-
attention with a single head. The final output is also computed by the VQC.
The others are the same as the Multihead Attention [44].

(a) Quantum self-attention with single head.

(b) Quantum self-attention with multihead.

Figure 8: The proposed quantum self-attention architecture with single head
and multihead. The low-qubit VQC as shown in Figure 2.

Figure 9: The proposed QTransformer-TTS architecture. The quantum self-
attention with multihead as shown in Figure 8(b).

With the quantum self-attention with multihead, the quantum Transformer-TTS
(QTransformer-TTS) will also be implemented. As shown in Fig. 9, change the
Multihead Attention to the quantum self-attention with multihead.

TABLE I: The Results of QM5 Compared with M5. In QM5, compared low-qubit VQC with VQC in different qubits. Model | Type of VQC | Number of qubits | Batch size | Test accuracy | Training time
---|---|---|---|---|---
M5 | - | - | 64 | 63.75% | -
- | - | 256 | 61.53% | -
QM5 | VQC | 2 | 64 | 30.57% | 741s
VQC | 2 | 256 | 29.13% | 2905s
Low-qubit VQC | 2 | 64 | 47.69% | 79s
Low-qubit VQC | 2 | 256 | 45.61% | 278s
VQC | 4 | 64 | 32.61% | 1659s
VQC | 4 | 256 | 30.18% | 6601s
Low-qubit VQC | 4 | 64 | 48.50% | 142s
Low-qubit VQC | 4 | 256 | 47.03% | 550s
VQC | 8 | 64 | 33.93% | 3521s
VQC | 8 | 256 | 31.74% | 13273s
Low-qubit VQC | 8 | 64 | 50.11% | 346s
Low-qubit VQC | 8 | 256 | 48.97% | 1265s

##  VI Experiments

###  VI-A Dataset

LJSpeech [45]. This dataset consists of 13,100 short audio clips of a single
speaker reading passages from seven nonfiction books. The texts are normalized
and inserted with the beginning character (space) and end character (a
period), e.g., “There are 16 apples” is converted to“ there are sixteen
apples.”. The LJSpeech dataset is randomly divided into two sets: 12,600
samples for training and 500 samples for testing.

Google Speech Command-V2 [46]. For spoken word recognition, we use the 35
class setting, including the following frequent speech commands: ‘left’, ‘go’,
‘yes’, ‘down’, ‘up’, ‘on’, ‘right’, ‘no’, ‘off’, ‘stop’, and so on with a
total of 84,843 training examples and 11,005 testing examples.

###  VI-B Settings

Quantum Environment. The simulation of quantum computation runs on the CPU. We
limit the number of qubits to small numbers (2 to 10). The $m$ of
$\theta_{1},...,\theta_{m}$ is equal to 3 in VQC and low-qubit VQC.

Model Architecture. The QTransformer-TTS has two layers. It replaces self-
attention with QAttention, the other settings follow [44]. The QGRU is a
substitute for GRU in QTacotron. QTacotron has same the same parameter as as
[43]. The QM5 replaces the 1D convolution (Conv1d) with 1D QConv (QConv1d).
The QTransformer-TTS, QM5, and QTacotron use the low-qubit VQC. As we
discussed in Sec. IV-C, VQC requires similar qubits as the input dimensions,
and it is considered large for quantum devices. Therefore, the low-qubit VQC
in QTransformer-TTS and QTacotron, cannot be compared with VQC. The QConv1d
with VQC or low-qubit VQC of QM5, in which the number of qubits is the same as
the kernel size, follows the [21, 20]. The linear transformation in low-qubit
VQC forms a fully connected layer.

Training. Adam is the optimizer in all experiments, with $\beta_{1}=0.9$,
$\beta_{2}=0.99$. The batch size in the training process of QTransformer-TTS
and QTacotron is $2$. QM5 can be trained in a large batch size (e.g., $64$,
$256$) using QConv1d with our proposed low-qubit VQC.

Compared Schemes. In the experiments, the QTransformer-TTS, QM5, and QTacotron
models are compared with the original Transformer-TTS, M5, and Tacotron. The
setups of training for the QTransformer-TTS, QM5, and QTacotron are the same
as the original Transformer-TTS, M5, and Tacotron. It will also compare the
results of different qubits in the experiments of QTransformer-TTS, QM5, and
QTacotron.

Evaluation. In Speech Command Recognition (SCR) experiments, the metric of
accuracy is the same as [40]. We use the mean opinion score (MOS) to measure
generated voices’ qualities in TTS experiments. Each sentence is judged by ten
males and ten females native English speakers.

###  VI-C Results and Analysis

QM5. In Fig. 10(a), the yellow curve is the training loss of QM5 with VQC.
Notably, it cannot stably converge as the iterative process increases. The
black curve, which is the QM5 with low-qubit VQC, can train in convergence and
stability, indicating that using low-qubit VQC can handle the Barren Plateaus
challenge. In Fig. 10(b), the training of QM5 with low-qubit VQC can improve
convergence and stability compared to the M5. With more qubits in QM5, the QM5
model will converge faster.

(a) Loss of VQC and low-qubit VQC in QM5.

(b) Loss of M5 and QM5(with low-qubit VQC).

Figure 10: The loss curves of M5 and QM5. (a) The training losses of QM5 with
VQC and low-qubit VQC. (b) The training loss curve of QM5 with different
numbers of qubits, compared with M5.

In Table I, with additional qubits, the QM5 obtains better performance.
Compared with VQC, the low-qubit VQC is $16\%$ higher. Table I shows that the
low-qubit VQC is faster than the VQC. For the same batch size, the low-qubit
VQC can reduce the training by $10$ times compared to the VQC. Table I also
shows the accuracy between the M5 and the QM5. The QM5, when using the low-
qubit VQC, can gain $16\%$ improvement in accuracy compared to the VQC,
indicating that the QM5 can work in the SCR task.

As shown in Fig. 11, the QConv1d can extract the feature from an input like
Conv1d. The feature pattern extracted by QConv1d is similar to that by Conv1d.

(a) Feature map of M5.

(b) Feature map of QM5.

Figure 11: The feature map output from M5 and QM5.(a) The feature map output
from Conv1d in M5.(b) The feature map output from QConv1d in QM5.

QTacotron. As shown in Fig. 12(a), the QTacotron with the low-qubit VQC
converges faster in training than Tacotron (the blue curve), indicating that
the low-qubit VQC can operate well in QGRU, like the low-qubit VQC in QConv1d.
The $0\%$ to $20\%$ of the training process indicates that the QTacotron
convergence is faster than that of the QTacotron. Following the training
steps, the QTacotron has less oscillation than the Tacotron. In Fig. 12(a),
QTacotron with different qubits has a similar loss convergence. However, it
shows that with more qubits, the QTacotron can converge better and with less
loss.

TABLE II: Results of Tacotron, QTacotron, Transformer-TTS, and QTransformer-TTS. QTacotron and QTransformer-TTS run in different numbers of qubits with low-qubit VQC. Model | Type of VQC | #Qubits | MOS
---|---|---|---
Tacotron | - | - | 3.09
QTacotron | Low-qubit VQC | 2 | 2.22
Low-qubit VQC | 4 | 2.56
Low-qubit VQC | 8 | 2.67
Transformer-TTS | - | - | 2.91
QTransformer-TTS | Low-qubit VQC | 2 | 2.19
Low-qubit VQC | 4 | 2.43
Low-qubit VQC | 8 | 2.52

Table II shows the MOS with Tacotron and QTacotron. QTacotron’s performance
will improve with more qubits of the low-qubit VQC. Therefore, the QTacotron
can achieve better results when using the low-qubit VQC. With eight qubits of
the low-qubit VQC, the MOS of QTacotron is about $2.67$, indicating that the
QTacotron with low-qubit VQC can run in the TTS task.

(a) Loss of Tacotron and QTacotron.

(b) Loss of Transformer-TTS and QTransformer-TTS.

Figure 12: The loss curves of Tacotron, QTacotron, Transformer-TTS, and
QTransformer-TTS. For QTacotron and QTransformer-TTS, both are compared with
different number of qubits. (a) The loss curves of Tacotron and QTacotron. (b)
The loss curves of Transformer-TTS and QTransformer-TTS.

In addition to the training process and performance, we provide the mel-
spectrogram obtained from the Tacotron and QTacotron. As shown in Fig.
13(a)(b), similar mel-spectrograms are obtained from the Tacotron and
QTacotron. We can observe that the QTacotron with low-qubit can achieve the
approximate performance as Tacotron.

QTransformer-TTS. In the training process, QTransformer-TTS outperforms
Transformer-TTS, as shown in Fig. 12(b). Compared to the training loss curve
in Fig. 12(b), the QTransformer-TTS converges faster than the Transformer-TTS
between $0\%$ and $40\%$ of the training process. Furthermore, QTransformer-
TTS, like the QTacotron, gives similar loss convergence in different qubits.
With more qubits of the low-qubit VQC, the QTransformer-TTS obtains lower loss
in the training process, indicating that the low-qubit VQC can work well in
QAttention, like the low-qubit VQC in QConv1d and QGRU.

Table II also shows the MOS with Transformer-TTS and QTransformer-TTS. We can
see that QTransformer-TTS performs better when using the low-qubit VQC. With
more qubits of the low-qubit VQC, QTransformer-TTS will get better
performance. With eight qubits of the low-qubit VQC, the MOS of QTransformer-
TTS is about $2.52$. It means that the QTransformer-TTS with low-qubit VQC can
run in the TTS task with a low-qubit quantum device.

Furthermore, in Fig. 13(c)(d), we show the mel-spectrogram generated by the
Transformer-TTS and QTransformer-TTS. The mel-spectrogram of the Transformer-
TTS and QTransformer-TTS is similar. We conclude that QTransformer-TTS with
low-qubit VQC can achieve the approximate performance as Transformer-TTS.

(a) Mel-spectrogram of Tacotron.

(b) Mel-spectrogram of QTacotron.

(c) Mel-spectrogram of Transformer-TTS.

(d) Mel-spectrogram of QTransformer-TTS.

Figure 13: The mel-spectrogram between the Tacotron, QTacotron, Transformer-
TTS, and QTransformer-TTS. (a) Mel-spectrogram generate by Tacotron. (b) Mel-
spectrogram generate by QTacotron. (c) Mel-spectrogram generate by
Transformer-TTS. (d) Mel-spectrogram generate by QTransformer-TTS.

Ablation Study. The low-qubit VQC is compared in two scenarios: first, the LT
is implemented by FC layer or bilinear interpolation; second, the LT layer
with or without output clip operation. Table III shows the results.

TABLE III: Compared results of QM5, QTacotron, and QTransformer-TTS. QTacotron and QTransformer-TTS run with low-qubit VQC in different settings. Model | Methods of LT | Results
---|---|---
QM5 | use bilinear interpolation | 37.51%
use FC | 43.41%
QTacotron | use bilinear interpolation | 2.11
use FC | 2.47
QTransformer-TTS | use bilinear interpolation | 2.24
use FC | 2.32
QM5 | without output clip | 43.21%
with output clip | 50.11%
QTacotron | without output clip | 2.41
with output clip | 2.67
QTransformer-TTS | without output clip | 2.39
with output clip | 2.52

Implemented methods of LT. In QM5, QTacotron, and QTransformer-TTS, the LT is
implemented by FC layer or the LT implemented by bilinear interpolation. As
shown in Table III, the results of using the FC layer are better than bilinear
interpolation. In low-qubit VQC, using the FC layer outperforms the bilinear
interpolation when implementing the LT.

With or without output clip. In QM5, QTacotron, and QTransformer-TTS, the
results of LT with output clip operation are compared with LT without output
clip operation. As shown in Table III shows, the results of using output clip
are better than without output clip. In low-qubit VQC, the results of LT with
output clip operation outperform LT without output clip operation.

##  VII Conclusion

This study introduced a novel VQC, called low-qubit VQC. It can address the
Barren Plateaus challenge and lower the requirement of qubits compared to VQC.
It is crucial because VQCs are used in low-qubit quantum devices in the NISQ
era. Low-qubit VQC can be used to implement QNN models in low-qubit quantum
devices. Based on low-qubit VQC, we propose the QSpeech toolkit. QSpeech fully
uses PyTorch and PennyLane as a QML engine and significantly simplifies the
implementation and training of the entire speech application. Our QSpeech
provided various quantum neural layers and hybrid quantum-classical neural
networks. The experimental results showed that the hybrid quantum-classical
neural networks with low-qubit VQC can be successfully used in speech
applications.

## Acknowledgment

This paper is supported by the Key Research and Development Program of
Guangdong Province under grant No.2021B0101400003. Corresponding author is
Jianzong Wang from Ping An Technology (Shenzhen) Co., Ltd (jzwang@188.com).

## References

  * [1] D. Wang, A. Sundaram, R. Kothari, A. Kapoor, and M. Roetteler, “Quantum algorithms for reinforcement learning with a generative model,” in _International Conference on Machine Learning_. PMLR, 2021, pp. 10 916–10 926.
  * [2] T. Birdal, V. Golyanik, C. Theobalt, and L. J. Guibas, “Quantum permutation synchronization,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2021, pp. 13 122–13 133.
  * [3] J. Preskill, “Quantum computing in the nisq era and beyond,” _Quantum_ , vol. 2, p. 79, 2018.
  * [4] M. Benedetti, E. Lloyd, S. Sack, and M. Fiorentini, “Parameterized quantum circuits as machine learning models,” _Quantum Science and Technology_ , vol. 4, no. 4, p. 043001, 2019.
  * [5] A. Bisarya, W. El Maouaki, S. Mukhopadhyay, N. Mishra, S. Kumar, B. K. Behera, P. K. Panigrahi, and D. De, “Breast cancer detection using quantum convolutional neural networks: A demonstration on a quantum computer,” _medRxiv_ , 2020.
  * [6] C.-H. H. Yang, J. Qi, S. Y.-C. Chen, P.-Y. Chen, S. M. Siniscalchi, X. Ma, and C.-H. Lee, “Decentralizing feature extraction with quantum convolutional neural network for automatic speech recognition,” in _ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2021, pp. 6523–6527.
  * [7] M. A. Nielsen and I. Chuang, “Quantum computation and quantum information,” 2002\.
  * [8] I. Tavernelli, A. Robert, P. Barkoutsos, and S. Woerner, “Resource-efficient quantum algorithm for protein folding,” _Bulletin of the American Physical Society_ , vol. 65, 2020.
  * [9] C. Figgatt, D. Maslov, K. A. Landsman, N. M. Linke, S. Debnath, and C. Monroe, “Complete 3-qubit grover search on a programmable quantum computer,” _Nature communications_ , vol. 8, no. 1, pp. 1–9, 2017.
  * [10] S. Lloyd, M. Mohseni, and P. Rebentrost, “Quantum principal component analysis,” _Nature Physics_ , vol. 10, no. 9, pp. 631–633, 2014.
  * [11] Z. Li, Z. Chai, Y. Guo, W. Ji, M. Wang, F. Shi, Y. Wang, S. Lloyd, and J. Du, “Resonant quantum principal component analysis.” _Science Advances_ , vol. 7, no. 34, 2021.
  * [12] X. Qu, J. Wang, and J. Xiao, “Evolutionary algorithm enhanced neural architecture search for text-independent speaker verification,” _Proc. Interspeech 2020_ , pp. 961–965, 2020.
  * [13] N. Zhang, J. Wang, W. Wei, X. Qu, N. Cheng, and J. Xiao, “Cacnet: Cube attentional cnn for automatic speech recognition,” in _2021 International Joint Conference on Neural Networks (IJCNN)_. IEEE, 2021, pp. 1–7.
  * [14] J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, “Barren plateaus in quantum neural network training landscapes,” _Nature communications_ , vol. 9, no. 1, pp. 1–6, 2018.
  * [15] C. Zhao and X.-S. Gao, “Analyzing the barren plateau phenomenon in training quantum neural networks with the zx-calculus,” _Quantum_ , vol. 5, p. 466, 2021.
  * [16] Z. Hong, J. Wang, X. Qu, J. Liu, C. Zhao, and J. Xiao, “Federated Learning with Dynamic Transformer for Text to Speech,” in _Proc. Interspeech 2021_ , 2021, pp. 3590–3594.
  * [17] C. Zhao, J. Wang, X. Qu, H. Wang, and J. Xiao, “r-g2p: Evaluating and enhancing robustness of grapheme to phoneme conversion by controlled noise introducing and contextual information incorporation,” in _ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2022, pp. 6197–6201.
  * [18] V. Bergholm, J. Izaac, M. Schuld, C. Gogolin, M. S. Alam, S. Ahmed, J. M. Arrazola, C. Blank, A. Delgado, S. Jahangiri _et al._ , “Pennylane: Automatic differentiation of hybrid quantum-classical computations,” 2018.
  * [19] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga _et al._ , “Pytorch: An imperative style, high-performance deep learning library,” _Advances in Neural Information Processing Systems_ , vol. 32, pp. 8026–8037, 2019.
  * [20] S. Y.-C. Chen, S. Yoo, and Y.-L. L. Fang, “Quantum long short-term memory,” _Bulletin of the American Physical Society_ , vol. 66, 2021.
  * [21] I. Cong, S. Choi, and M. D. Lukin, “Quantum convolutional neural networks,” _Nature Physics_ , vol. 15, no. 12, pp. 1273–1278, 2019.
  * [22] I. Kerenidis, J. Landman, and A. Prakash, “Quantum algorithms for deep convolutional neural networks,” in _International Conference on Learning Representations_ , 2019.
  * [23] Z. Zhang, D. Chen, J. Wang, L. Bai, and E. R. Hancock, “Quantum-based subgraph convolutional neural networks,” _Pattern Recognition_ , vol. 88, pp. 38–49, 2019.
  * [24] A. A. Melnikov, L. E. Fedichkin, and A. Alodjants, “Predicting quantum advantage by quantum walk with convolutional neural networks,” _New Journal of Physics_ , vol. 21, no. 12, p. 125002, 2019.
  * [25] Z. Hong, J. Wang, X. Qu, X. Zhu, J. Liu, and J. Xiao, “Quantum convolutional neural network on protein distance prediction,” in _2021 International Joint Conference on Neural Networks (IJCNN)_. IEEE, 2021, pp. 1–8.
  * [26] M. Hibat-Allah, M. Ganahl, L. E. Hayward, R. G. Melko, and J. Carrasquilla, “Recurrent neural network wave functions,” _Physical Review Research_ , vol. 2, no. 2, p. 023358, 2020.
  * [27] C. Roth, “Iterative retraining of quantum spin models using recurrent neural networks,” 2020.
  * [28] J. J. Hopfield, “Neural networks and physical systems with emergent collective computational abilities,” _Proceedings of the national academy of sciences_ , vol. 79, no. 8, pp. 2554–2558, 1982.
  * [29] H. Tang, Z. Feng, Y.-H. Wang, P.-C. Lai, C.-Y. Wang, Z.-Y. Ye, C.-K. Wang, Z.-Y. Shi, T.-Y. Wang, Y. Chen _et al._ , “Experimental quantum stochastic walks simulating associative memory of hopfield neural networks,” _Physical Review Applied_ , vol. 11, no. 2, p. 024020, 2019.
  * [30] F. Li, S. Zhao, and B. Zheng, “Quantum neural network in speech recognition,” in _6th International Conference on Signal Processing, 2002._ , vol. 2. IEEE, 2002, pp. 1267–1270.
  * [31] F. Li and G. Xu, “Quantum bp neural network for speech enhancement,” in _2009 Asia-Pacific Conference on Computational Intelligence and Industrial Applications (PACIIA)_ , vol. 2. IEEE, 2009, pp. 389–392.
  * [32] L. Fu and J. Dai, “A speech recognition based on quantum neural networks trained by ipso,” in _2009 International Conference on Artificial Intelligence and Computational Intelligence_ , vol. 2. IEEE, 2009, pp. 477–481.
  * [33] K. Mitarai, M. Negoro, M. Kitagawa, and K. Fujii, “Quantum circuit learning,” _Physical Review A_ , vol. 98, no. 3, p. 032309, 2018.
  * [34] Y. LeCun, Y. Bengio, and G. Hinton, “Deep learning,” _nature_ , vol. 521, no. 7553, pp. 436–444, 2015.
  * [35] X. Qu, J. Wang, and J. Xiao, “Enhancing data-free adversarial distillation with activation regularization and virtual interpolation,” in _ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2021, pp. 3340–3344.
  * [36] J. Liu, J. Wang, X. Qu, B. Chen, Z. Wei, and J. Xiao, “Automatic joint optimization of algorithm-level compression and compiler-based acceleration with reinforcement learning for dnn in edge devices,” in _2021 International Joint Conference on Neural Networks (IJCNN)_. IEEE, 2021, pp. 1–8.
  * [37] N. Killoran, T. R. Bromley, J. M. Arrazola, M. Schuld, N. Quesada, and S. Lloyd, “Continuous-variable quantum neural networks,” _Physical Review Research_ , vol. 1, no. 3, p. 033063, 2019.
  * [38] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in _ICLR (Poster)_ , 2015.
  * [39] L. Li, J. Wang, X. Qu, and J. Xiao, “Communication-memory-efficient decentralized learning for audio representation,” in _2021 International Joint Conference on Neural Networks (IJCNN)_. IEEE, 2021, pp. 1–8.
  * [40] W. Dai, C. Dai, S. Qu, J. Li, and S. Das, “Very deep convolutional neural networks for raw waveforms,” in _2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2017, pp. 421–425.
  * [41] Y. Tang, J. Wang, X. Qu, and J. Xiao, “Contrastive learning for improving end-to-end speaker verification,” in _2021 International Joint Conference on Neural Networks (IJCNN)_. IEEE, 2021, pp. 1–7.
  * [42] Z. Hong, J. Wang, W. Wei, J. Liu, X. Qu, B. Chen, Z. Wei, and J. Xiao, “When hearing the voice, who will come to your mind,” in _2021 International Joint Conference on Neural Networks (IJCNN)_. IEEE, 2021, pp. 1–6.
  * [43] Y. Wang, R. Skerry-Ryan, D. Stanton, Y. Wu, R. J. Weiss, N. Jaitly, Z. Yang, Y. Xiao, Z. Chen, S. Bengio _et al._ , “Tacotron: Towards end-to-end speech synthesis,” _Proc. Interspeech 2017_ , pp. 4006–4010, 2017.
  * [44] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, and I. Polosukhin, “Attention is all you need,” in _NIPS_ , 2017\.
  * [45] K. Ito _et al._ , “The lj speech dataset,” 2017.
  * [46] P. Warden, “Speech commands: A dataset for limited-vocabulary speech recognition,” 2018.