# Efficient Computation of Map-scale Continuous Mutual Information on Chip in
Real Time

Keshav Gupta, Peter Zhi Xuan Li, Sertac Karaman, Vivienne Sze
Massachusetts Institute of Technology, Cambridge, Massachusetts 02139
Emails<EMAIL_ADDRESS>
###### Abstract

Exploration tasks are essential to many emerging robotics applications,
ranging from search and rescue to space exploration. The planning problem for
exploration requires determining the best locations for future measurements
that will enhance the fidelity of the map, for example, by reducing its total
entropy. A widely-studied technique involves computing the Mutual Information
(MI) between the current map and future measurements, and utilizing this MI
metric to decide the locations for future measurements. However, computing MI
for reasonably-sized maps is slow and power hungry, which has been a
bottleneck towards fast and efficient robotic exploration. In this paper, we
introduce a new hardware accelerator architecture for MI computation that
features a low-latency, energy-efficient MI compute core and an optimized
memory subsystem that provides sufficient bandwidth to keep the cores fully
utilized. The core employs interleaving to counter the recursive algorithm,
and workload balancing and numerical approximations to reduce latency and
energy consumption. We demonstrate this optimized architecture with a Field-
Programmable Gate Array (FPGA) implementation, which can compute MI for all
cells in an entire 201-by-201 occupancy grid (e.g., representing a
20.1m-by-20.1m map at 0.1m resolution) in 1.55 ms while consuming 1.7 mJ of
energy, thus finally rendering MI computation for the whole map real time and
at a fraction of the energy cost of traditional compute platforms. For
comparison, this particular FPGA implementation running on the Xilinx
Zynq-7000 platform is two orders of magnitude faster and consumes three orders
of magnitude less energy per MI map compute, when compared to a baseline GPU
implementation running on an NVIDIA GeForce GTX 980 platform. The improvements
are more pronounced when compared to CPU implementations of equivalent
algorithms.

##  I Introduction

Figure 1: (top) Top-level hardware architecture comprised of Address
Translation Units (ATU), crossbar interconnects, occupancy grid map, MI map,
and $16$ FCMI computation cores. (bottom) An example input occupancy grid,
corresponding MI map contours generated from a reference C++ implementation
and our FPGA implementation, and the absolute error between the two.

Exploration problems arise in an increasing set of applications, particularly
as robotic vehicles are being utilized in deep sea and space exploration
missions as well as search and rescue scenarios. In most applications, the
goal is to create a map of the environment, e.g., from range measurements, as
quickly as possible. A key problem, then, is to find the best locations for
future measurements in order to achieve this goal.

Consider an instance of this problem involving a robot equipped with a laser
scanner. A typical computation cycle in this process involves making a scan of
the robot’s surroundings, updating the partially-unknown map with the
information from the scan, determining the best location for the next scan and
moving to that location. An information-theoretic approach to determining the
best scan locations involves finding the location with the highest Mutual
Information (MI) between the potential future measurement and the existing
map. Unfortunately, computing MI at all map locations is slow and power
hungry, even for small-size maps [1], which makes it especially infeasible for
usage in medium to small scale robots with the energy budget of a few
kiloJoules. Most heuristics used to reduce the scale of this computation, such
as calculating the MI at sparse points instead of the entire map, fail to
provide guarantees, e.g., on the quality of the resulting motion [2, 3, 4].

To accelerate the computation of MI, several novel algorithms were proposed.
For instance, the efficient computation of Cauchy-Schwarz Quadratic Mutual
Information (CSQMI) is introduced in [5, 6]. In addition, the high-throughput
computation of Shannon Mutual information is proposed in the Fast Shannon
Mutual Information (FSMI) algorithm [7]. More recently, the computation of
Shannon MI is reformulated recursively and continuously along the range
measurements in the Fast Continuous Mutual Information (FCMI) algorithm [8],
which leads to a significant reduction in amount of compute needed for MI
across an entire map. Although these algorithms reduce the time complexity for
MI computation, their optimized implementations on CPUs and GPUs still achieve
a throughput that is orders-of-magnitude below what is required for autonomous
exploration, especially in large environments. This motivates the design of
specialized hardware accelerators for computing MI.

Hardware accelerators have been proposed for robotics applications such as
motion-planning [9, 10, 11, 12] and visual-inertial odometry [13, 14]. To our
best knowledge, the first and only hardware accelerator for computing MI was
proposed in [15], which contains 16 high-throughput computation cores that
compute FSMI. For moderately-sized environments (_i.e._ , a map size of
$200\times 200$ to $512\times 512$ cells), this accelerator computes MI for
the entire map in _near_ real time (_i.e._ , less than 10Hz), which is far
below the rate in which the robot makes range measurements (around
$30$-$60$Hz).

The computation of mutual information may seem like a task that can be easily
parallelized. One might be tempted to think, for example, that computational
cores dedicated to various parts of the map can work in parallel to compute
the mutual information for all cells in the map all at the same time. Indeed,
under certain widely-utilized assumptions, the computation decouples and cells
can be evaluated separately. Unfortunately, however, the computation does not
scale, not because it can not be decoupled, but because it becomes impossible
to feed all cores with the data that they need for their computations. More
specifically, all cores need to access different portions of the map, which
resides in memory, all at the same time. The memory bandwidth does not scale
to keep all cores busy. In fact, we show in this paper that naively scaling up
the number cores does not lead to a fast computation of the mutual information
metric, as most cores end up underutilized due to insufficient memory
bandwidth.

The main contribution of this paper is a scalable, multi-core hardware
architecture for computing mutual information for a 2D occupancy grid map.
This architecture has two key components: (i) a memory management system that
provides sufficient bandwidth to all cores, and (ii) a low-latency, energy-
efficient core (_i.e._ processing element) design that fully utilizes all
available bandwidth. Specifically, the memory management system exploits the
memory access pattern of the cores and splits the storage of the occupancy
grid map into smaller memory banks so that all cores access different banks
every cycle. The core contains a feedback architecture designed to handle the
dependencies that are introduced by the recursive nature of the FCMI algorithm
[8], so that all compute units within the core are fully utilized. Each core
is time-interleaved to compute the FCMI along multiple rays in order to avoid
stalls associated with the feedback architecture.

The proposed hardware architecture is implemented with 16 cores on a Xilinx
Zynq-7000 XC7Z045 FPGA. For a robot equipped with a 60-ray range sensor
utilizing an occupancy grid map of size $201\times 201$, our FPGA
implementation computes MI for the entire map in 1.55 ms (647 Hz, or much
faster than real time), which is $381\times$ faster than an optimized C++
implementation running on a ARM Cortex A57 (embedded cellphone processor) and
$71\times$ faster than an optimized CUDA implementation running on a NVIDIA
GTX 980 GPU. In addition, our FPGA implementation consumes $1.7$ mJ of energy
per MI map compute, which is $1314\times$ lower than ARM Cortex A57 and
$2650\times$ lower than NVIDIA GTX 980 GPU. A CPU on a robot can therefore
offload Mutual Information computation to such an accelerator, while
interfacing with the sensors and actuators and performing motion planning. The
extremely low latency and energy consumption of our accelerator allows robots
to explore efficiently without relying on heuristic approximations even with
the energy budget of a few Joules.

##  II Preliminaries

###  II-A The FCMI Algorithm

Occupancy grid map [16] is a common probabilistic representation of the 2D
environment and used in several well known autonomous exploration systems,
such as the frontier-exploration algorithm [17] and Shannon MI-based mapping
algorithms [18, 19, 20, 21, 22]. In an information theoretic framework, MI is
used to quantify the expected information gain from potential range
measurements in an unknown environment. The unknown environment is represented
by an occupancy grid map, which is updated continuously by the robot
throughout an exploration trial. The FCMI algorithm [8] computes the mutual
information between the map and potential range measurement rays emanating
from the sensors (_e.g._ , LiDAR) on such robot.

Let $Z=\\{Z_{\theta_{1}},Z_{\theta_{2}}\ldots\\}$ represent range measurements
made by sensor rays with angles $\Theta=\\{\theta_{1},\theta_{2},\ldots\\}$
emanating from robot’s location $x$. Let $\Delta_{\theta}$ denote the angular
resolution of the sensor rays, which is constant across all rays. The Shannon
MI $I_{x}(M;Z)$ between the measurements $Z$ made at location $x$ and map $M$
is defined as

$I_{x}(M;Z)=\sum_{\theta\in\Theta}I_{x}(M;Z_{\theta})\Delta_{\theta}.$ (1)

Note that the MI between a single range measurement $Z_{\theta}$ and map $M$
can be written as

$I_{x}(M;Z_{\theta})=h_{x}(Z_{\theta})-h_{x}(Z_{\theta}|M),$ (2)

where $h_{x}(Z_{\theta})$ is the entropy of the range measurement, and
$h_{x}(Z_{\theta}|M)$ is the conditional entropy of the range measurement
given the map.

Let $S$ denote the coordinates of cells in the occupancy grid map that are
intersected by a sensor ray of infinite length with angle $\theta$. Suppose
that the same sensor ray passes through location $x$ such that $x\in S$. For
the cell at location $x$, let $o$ denote its occupancy probability, which is
derived from previous range measurements. We define four expectations
$\alpha_{1},\beta_{1},\alpha_{0}$ and $\beta_{0}$ (as define in [8])
associated with the same cell as follows:

$\displaystyle\alpha_{1}$
$\displaystyle=e^{-\lambda_{m}w}\left((\alpha^{\prime}_{1}+\lambda
w\beta^{\prime}_{1})+w(\alpha^{\prime}_{0}+\lambda_{m}w\beta^{\prime}_{0})\right),$
$\displaystyle+\lambda_{m}^{-1}\left(\gamma_{3}(\lambda_{m}w)-\gamma_{2}(\lambda_{m}w)\log\lambda_{m}\right),$
(3a) $\displaystyle\beta_{1}$
$\displaystyle=e^{-\lambda_{m}w}\left(\beta^{\prime}_{1}+w\beta^{\prime}_{0}\right)+\lambda_{m}^{-1}\gamma_{2}(\lambda_{m}w),$
(3b) $\displaystyle\alpha_{0}$
$\displaystyle=e^{-\lambda_{m}w}\left(\alpha^{\prime}_{0}+\lambda_{m}w\beta^{\prime}_{0}\right)+\gamma_{2}(\lambda_{m}w)$
$\displaystyle-\gamma_{1}(\lambda_{m}w)\log\lambda_{m},$ (3c)
$\displaystyle\beta_{0}$
$\displaystyle=e^{-\lambda_{m}w}\beta^{\prime}_{0}+\gamma_{1}(\lambda_{m}w),$
(3d)

where
$\alpha^{\prime}_{1},\alpha^{\prime}_{0},\beta^{\prime}_{1},\beta^{\prime}_{0}$
denote the four expectations from the previous cell along the same sensor ray,
$w$ is the width of the current cell, $\Lambda$ is a large quantity used to
barely distort the scan measurements ($\Lambda=10^{7}$ for this work),
$\lambda=-\log(1-o)$, $\lambda_{m}=\min(\lambda,\Lambda)$ and
$\gamma_{s\in\\{1,2,3\\}}(\cdot)$ is the lower incomplete gamma function.

The entropy and conditional entropy in Eqn. (2) are approximated from the
aforementioned expectations as follows:

$\displaystyle h_{x}(Z_{\theta})$
$\displaystyle\approx\alpha_{1}\Delta_{\theta},$ (4a) $\displaystyle
h_{x}(Z_{\theta}|M)$
$\displaystyle\approx(1-\log\Lambda)\beta_{1}\Delta_{\theta}.$ (4b)

Due to the recursive computation of the four expectations along a sensor ray,
the entropy defined in Eqn. (4) for every cell within a map of size $H\times
H$ can be computed with a time complexity of $O(H^{2})$. Let $|\Theta|$ denote
the number of sensor rays emanating from the robot. Then, the mutual
information defined in Eqn. (1) for the entire map can be computed with a time
complexity of $O(|\Theta|H^{2})$.

###  II-B A Lower Bound for Latency

Let $n$ denote the number of cores in the hardware accelerator, which is
clocked from a single source with frequency $F$. Suppose that the robot makes
range measurements using $|\Theta|$ rays within an occupancy grid map of size
$H\times H$. Recall that the FCMI algorithm has a time complexity of
$O(|\Theta|H^{2})$ for computing the corresponding MI map of the same size. If
the hardware needs at least one cycle to compute the entropy defined in Eqn.
(4) for each cell along every sensor ray, the time between receiving an
occupancy grid and producing the resultant MI (the lower bound latency of the
system) can be determined from:

$t_{\text{LB}}=\frac{H^{2}|\Theta|}{nF},$ (5)

which has a unit of seconds, assuming unlimited memory bandwidth.

In this work, we present a hardware architecture with $n=16$ cores clocked at
$F=100$ MHz. In Section IV, we discuss how the achieved latency of the
proposed architecture scales with the number of cores $n$, and compare it with
the aforementioned lower bound.

##  III Hardware Architecture

###  III-A Architecture Overview

The proposed hardware architecture for FCMI computation is shown in Fig. 1
(top) for a 16-core implementation. The architecture contains a memory
subsystem that splits the occupancy grid map and MI map into 16 memory banks.
The Bresenham ray caster determines the cells that are intersected by each
ray. The Address Translation Unit (ATU) determines the memory storage location
for each cell in the memory subsystem. There are two crossbar interconnects
that function as switches for controlling the flow of data. The first crossbar
connects the ATUs to the memory subsystem for memory request generation. The
second crossbar connects the memory subsystem to the cores for data delivery
and MI map update. To achieve low-latency computation, we need to ensure that
(i) the memory subsystem can generate enough bandwidth to keep all cores
active, and (ii) the cores can effectively utilize all available bandwidth
without stalling.

###  III-B Memory Subsystem

(a)

(b)

(c)

Figure 2: (a) The micro-architecture of the FCMI core comprised of three
sections in the pipeline. (b) An example of the diagonal banking pattern with
4 banks, each in a different color. Any 4 consecutive cells here (examples
highlighted) are guaranteed to lie on distinct banks. (c) Time is represented
on the x axis and pipeline stages are represented on the y axis. Note that the
feedback section occupies Stages 11-18 of the pipeline. (left) Interleaving 8
rays (represented with distinct colors) into a core. (right) Computing only
one ray at a time, leading to a 7 cycle stall when waiting for results at the
output stage.

On the FPGA, memory units are implemented as dual-port SRAMs, which allows for
two memory accesses in each cycle. Since the MI map is updated every cycle
using one write port, only one read port from each SRAM is available for
delivering data to the cores. If multiple cores attempt to concurrently access
different cells stored within the same bank, a memory access conflict occurs
and some cores must stall.

Recall that the coordinates for each cell are computed using the Bresenham ray
casting algorithm [23], which we chose for its efficient hardware
implementations [24]. Since all 16 cores concurrently compute FCMI along rays
of the same angle as shown in Fig. 3, the cores are guaranteed to access
consecutive cells along a row or column of the map every cycle. This memory
access pattern is inherent to the algorithm and is highly predictable. In
order to exploit this property, the maps are partitioned and stored in 16
memory banks using a diagonal banking pattern, which ensures that any 16
consecutive cells within every row or column of the map are stored in distinct
banks [15]. An example of the diagonal banking pattern with 4 bank partitions
is shown in Fig. 2(b), which can be extended to 16 banks. The occupancy grid
map memory and the MI memory are both banked under this scheme.

Figure 3: Wrapping rays around map boundaries to balance core workload and
minimize latency.

###  III-C FCMI Core

The micro-architecture of the FCMI core is illustrated in Fig. 2(a). The cell
accepts a new occupancy value and a partial MI for one cell along the sensor
ray every cycle. The updated MI at the same cell is computed using the FCMI
formulation and written back to the MI map in the memory subsystem. The low-
latency, energy-efficient FCMI computation of each core is enabled by the
following:

### 1) Arithmetic Optimizations

We store and compute MI results using 32-bit fixed point numbers (20 integer,
12 fractional bits), since adders and multipliers implemented using fixed
point numbers tend to be faster and more energy efficient using the
specialized DSPs on the FPGA. The choice of 32-bit width comes from taking
into account the dynamic range of numbers needed to avoid under/over-flow in
this application. In addition, the core needs to evaluate the following
logarithms and reciprocals: $v=1-o$, $\lambda=-\log v$, $-\log\lambda$, and
$1/\lambda$ (from Eqn. (3)). They are all direct functions of the occupancy
probabilities $o$ which is uniformly quantized into $101$ values such that
$o\in\\{0,0.1,\ldots,1\\}$. Note that $o=0$ and $o=1$ are degenerate cases,
which are handled using pre-calculated results. A small look-up table is used
which store all values of logarithms and reciprocals. This look-up table
occupies only two Block RAM tiles per core on the FPGA.

Finally, the core also needs to evaluate an exponential to compute
$e^{-\lambda w}$, where $\lambda$ has one-to-one correspondence to the
occupancy probability $o$ (101 values) and $w$ is the cell width. To the best
of our knowledge, no such module exists that uses fixed point arithmetic that
would lead to a throughput of one result per clock cycle. To achieve such a
throughput, a look-up table or an approximation should be used. Note that $w$
is dependent on number of sensor rays $|\Theta|$ such that
$w=\left\lceil|\Theta|/8\right\rceil$ (accounting for symmetries w.r.t. the
axes). If $|\Theta|=60$, a look-up table mapping all possible combinations of
$(\lambda,w)$ to the corresponding $e^{-\lambda w}$ value would need to
contain $\lceil 101(60)/8\rceil=808$ values, which translates to $16$ Block
RAM tiles per core on the target FPGA. The FPGA does not have sufficient Block
RAM to support the look-up table for $16$ cores. In addition, since the size
of the look-up table is dependent on the angular resolution of the sensor
rays, the look-up table implementation would only be able to support a limited
set of sensor configurations. From formulation of FCMI, all $-\lambda w$
values lie inside $[-8,0]$. Thus, we use a piece-wise linear approximation to
the exponential function that splits the range into 16 uniform pieces and uses
a mean-squared error minimizing linear approximation to each. Such
approximation achieves a maximum relative error of $\pm 3\%$, while not
imposing a limit on sensor configurations, and consuming only one Block RAM
per core.

### 2) Hardware Pipeline

As evident from Eqn. (3), each cell update is a long chain of mathematical
operations. In order to reduce the critical path of the system to allow a
faster clock, we pipeline these operations into 19 stages and perform the
operations over consecutive clock cycles. The pipeline allows the system to be
clocked at $100\text{ MHz}$ instead of the combined combinational delay of the
system, which would be close to only $5\text{ MHz}$.

Since the quantities in Eqn. (3) are defined recursively, the pipeline
contains a feedback loop which increases the latency of the computation. A
naïve implementation of the pipeline executes the operations in Eqn. (3) for
each cell in 19 clock cycles such that each pipeline stage stalls for 18
cycles. As a solution, the 18 stages of the pipeline were split into three
sections: Preprocess, Feedback and Postprocess (as shown in Fig. 2(a)). This
reduces the feedback loop to only 8 stages in the Feedback section. In
addition, we introduce ray-interleaving to eliminate the remaining stalls
associated with the smaller feedback loop in Section III.

#### Preprocess

This stage uses the cell occupancy value to compute as many intermediate
values as possible without requiring information about the last cell. The
look-up tables are accessed to obtain $\lambda$, $-\log\lambda$ and
$1/\lambda$ defined in Section II-A. The exponential $e^{-\lambda w}$ and the
lower incomplete gamma function values $\gamma_{0}(\lambda w)$,
$\gamma_{1}(\lambda w)$, $\gamma_{2}(\lambda w)$ are also precalculated for
use in the feedback stage.

#### Feedback

This stage uses the intermediate values from the Preprocess stage and the
values of
$\\{\alpha^{\prime}_{1},\beta^{\prime}_{1},\alpha^{\prime}_{0},\beta^{\prime}_{0}\\}$
for the previous cell along the sensor ray to compute the updated
$\\{\alpha_{1},\beta_{1},\alpha_{0},\beta_{0}\\}$ values for the current cell.

#### Postprocess

This stage uses the updated $\\{\alpha_{1},\beta_{1},\alpha_{0},\beta_{0}\\}$
values for the current cell in addition to the current partial MI from the
memory subsystem to compute the updated MI of the same cell.

### 3) Ray Interleaving & Wrapping

The process of interleaving is represented visually in Fig. 2(c). From the
same figure, the recursion for the next cell along the sensor ray begins after
the previous cell along the same ray arrives at the final stage of the
Feedback section (stage 18). In the absence of ray interleaving, the core
needs to stall for 8 cycles for every cell along the sensor ray, which leads
to an increase in latency.

Due to ray interleaving, each core interleaves the computation of FCMI for 8
rays simultaneously to achieve a throughput of 1 cell per cycle. Therefore,
the amount of time a core spends processing each ray is proportional to its
length. In addition to ray interleaving, the workload is balanced across all
the cores that process rays with different lengths via wrapping. An example is
shown in Fig. 3 where the four cores that concurrently process the green rays
complete at different times. Note that orange rays wrap around the green rays
across the boundary of the map. By assigning the computation of the orange
rays to the same four cores after the green rays, the workload across all
cores becomes balanced.

##  IV Results

The proposed architecture is implemented on a Xilinx Zynq-7000 XC7Z045 FPGA
which is clocked at 100MHz and can store a maximum map size of $512\times
512$. We also evaluate the latency, energy per MI map compute and scalability
of the proposed architecture with respect to the number of cores. Finally, we
present the impacts that the proposed architecture can have in an exploration
experiment.

###  IV-A Accuracy

To evaluate the accuracy of our implementation, we generated contours lines
from MI maps for a variety of occupancy grid maps from real world exploration
trials.

In Fig. 1 (bottom), we compare the locations of high MI generated from our
FPGA implementation with those from a reference C++ implementation. The MI map
normalized to $[0,1]$ computed by the FPGA differs only marginally from the
one generated from the reference C++ implementation (by less than $0.05$
across all points), the major source of error being the piece-wise linear
approximation of the exponential function (Section III). We see that our
generated MI map is sufficiently accurate for real world exploration trials,
as seen from Section IV-E.

###  IV-B Latency

We evaluate the latency (_i.e._ the time elapsed between receiving an
occupancy grid and producing the corresponding MI map) of our FPGA
implementation that computes MI for an entire map of size $201\times 201$
cells with $60$ sensor rays. This map captures a robot in a $20.1\text{m
}\times 20.1\text{m}$ environment mapped at $10$ cm resolution. The FPGA
computes MI for the entire map in 1.55 ms using 16 cores.

Fig. 4 illustrates that the latency of our FPGA implementation reaches within
$2\%$ of the lower bound (Eqn. (5)) and scales linearly with the number of
cores. Without workload balancing or ray interleaving (Section III), the
latency would increase by $22.9\%$ or $8$ times respectively. Finally, if the
occupancy grid map and MI map are stored in one bank (_i.e._ no banking), the
latency does not decrease beyond what is achievable with a single core because
the memory cannot deliver enough bandwidth.

We compared the latency of our FPGA implementation that accelerates the FCMI
algorithm with the optimized C++ implementation on ARM Cortex A57 and Intel
Core i7-1065G7 (embedded mobile and desktop processors, respectively), and
optimized CUDA implementation on NVIDIA GTX 980 GPU. In addition, we compared
our FPGA implementation that accelerates the FCMI algorithm with a prior FPGA
implementation [15] that accelerates the FSMI algorithm. The results are shown
in Fig. 4. Our FPGA implementation is $71\times$ faster than the GTX 980,
$104\times$ faster than the FPGA implementation for FSMI, $190\times$ faster
than the Core i7 CPU, and $381\times$ faster than the Cortex A57 CPU.

The system latency scales in proportion to the number of pixels in the map
$(H^{2})$, and the number of sensor rays $(|\Theta|)$. For example, the system
computes the MI for a $256\times 256$ map in 2.51 ms, and for a $512\times
512$ map in 10.1 ms, and the relative speed-up over the baseline is
maintained.

###  IV-C Energy Consumption

We evaluate the energy per MI map compute (_i.e._ , the energy consumed by the
FPGA for every computation of an MI map from an occupancy grid) on the same
map as Section IV-B. Fig. 4 compares this metric of our FPGA implementation
with other hardware platforms, and illustrates how it scales with the number
of cores. The energy consumption is obtained post-place-and-route with
annotated switching activity.

With $16$ cores clocked at $100$ MHz, our FPGA implementation consumes only
$1.7$ mJ of energy per MI map compute, which is about $2650\times$ lower than
the GTX 980, $5360\times$ lower than the Core i7 CPU, $1314\times$ lower than
the Cortex A57 CPU and $187\times$ lower than the FPGA implementation for
FSMI.

Figure 4: (top, left) Impact of the number of cores on system latency. (top,
right) A comparison of latency achieved across platforms. (bottom, left)
Impact of the number of cores on energy consumed per MI compute. [Errata:
Careful re-simulation shows updated (better) numbers for Energy per MI Update.
The numbers from the original publication are shown by the black dashed line.]
(bottom, right) A comparison of the energy consumption per MI map compute
across platforms. FSMI (FPGA) is the previous work [15] and FCMI (FPGA) is
this work. Both architectures contain 16 cores and are implemented on the same
FPGA.

###  IV-D Resource Utilization

TABLE I: Resource utilization of the proposed architecture with 16 cores on a Xilinx Zynq-7000 XC7Z045 FPGA that stores a maximum map size of $512\times 512$. Module | LUT | LUT | FF | BRAM | DSP
---|---|---|---|---|---
(as Logic) | (as RAM)
Ray Caster | 1258 | 0 | 56 | 0 | 0
MI Memory | 10979 | 7 | 1123 | 232 | 0
Occ. Memory | 2407 | 0 | 282 | 56 | 0
Cores (16x) | 2970 | 170 | 2339 | 2 | 45
Total | 61729 | 2699 | 38991 | 320 | 720
Total% | 29% | 4% | 9% | 59% | 80%

Table I shows the resource utilization breakdown of our FPGA implementation
with 16 cores. The utilization of the entire system scales linearly with
respect to the number of cores, which offers great flexibility when
implemented on systems with different computational resources and energy
consumption. Although the size of the crossbar interconnect scales
quadratically with the number of cores, it occupies a very small fraction of
the design and this growth does not pose a problem for up to 16 cores.

###  IV-E Impact on Exploration

We simulate a model robot using the González-Baños-Latombe exploration
strategy [25] to explore a sample map sized at $201\times 201$ cells. The
robot makes a LiDAR scan and updates its occupancy grid at 30 Hz _i.e._ every
33 ms (a typical LiDAR sensor rate). Our system finishes the MI compute in
1.55 ms and remains idle while waiting for an update to the occupancy grid.
Since the other platforms require more than 33 ms to finish the MI compute,
they compute on the latest available occupancy grid. Fig. 5 (top) shows map
entropy of an example map as a function of trajectory length. It is evident
that being able to compute the MI map at a faster rate leads to a quicker drop
in map entropy.

More importantly, we measure the MI compute energy consumption during
exploration on several platforms while performing MI computes at the same rate
(_i.e._ , 1.7 Hz which matches the Cortex A57 CPU) in order to highlight the
energy efficiency of our architecture. From Fig. 5 (bottom), our system
completes the exploration task while spending around 5 J of energy on MI
compute, roughly three orders of magnitude lower than the compared platforms.
For context on the scale of these numbers, a typical cellphone battery holds
about 40 kJ of energy [26], which would drain after about 5-6 map exploration
trials running on its embedded processor. Our system, in comparison, would be
able to perform around 8,000 such trials on that energy budget.

Figure 5: (top) Map entropy of an example map as a function of trajectory
length using GBL exploration [25] with MI computed at the maximum rates
achievable on several platforms. (bottom) Map entropy on the same map as a
function of total energy spent on MI compute with MI being computed at the
rate matching the Cortex A57 CPU on all platforms in order to highlight the
energy efficiency of our system.

##  V Conclusions

In this work, we presented a novel scalable, multi-core architecture that
exploits the recursive nature of the FCMI algorithm. The key contributions in
this architecture are a memory subsystem that provides sufficient bandwidth to
all cores, and a low-latency, energy-efficient core design that interleaves
the computation across multiple sensor rays in order to utilize all available
memory bandwidth. We implemented the proposed architecture with 16 cores on an
FPGA, which computes the _entire_ mutual information map for a given occupancy
grid map in 1.55 ms (_i.e._ , 647 Hz) while consuming only 1.7 mJ of energy
per MI map compute. Not only does such a result lead to a significant
reduction in time, trajectory and energy consumption during autonomous
exploration, it also enables Mutual Information based exploration for medium
and small scale robots for the first time. In addition, the latency and FPGA
resource usage of the proposed architecture is scalable with the number of
cores, which makes the architecture re-configurable for both high performance
and lower energy applications, ranging from high-speed rovers to low-power
micro-drones.

Acknowledgements. This work was partially funded by the NSF RTML 1937501 and
the NSF CPS 1837212.

## References

  * [1] C. Cadena, L. Carlone, H. Carrillo, Y. Latif, D. Scaramuzza, J. Neira, I. Reid, and J. Leonard, “Past, present, and future of simultaneous localization and mapping: Towards the robust-perception age,” IEEE Transactions on Robotics, vol. 32, no. 6, p. 1309–1332, 2016.
  * [2] W. Burgard, M. Moors, C. Stachniss, and F. E. Schneider, “Coordinated multi-robot exploration,” IEEE Transactions on robotics, vol. 21, no. 3, pp. 376–386, 2005.
  * [3] H. H. González-Banos and J.-C. Latombe, “Navigation strategies for exploring indoor environments,” The International Journal of Robotics Research, vol. 21, no. 10-11, pp. 829–848, 2002.
  * [4] D. Holz, N. Basilico, F. Amigoni, S. Behnke, et al., “A Comparative Evaluation of Exploration Strategies and Heuristics to Improve Them,” in ECMR, pp. 25–30, 2011.
  * [5] B. Charrow, S. Liu, V. Kumar, and N. Michael, “Information-theoretic mapping using Cauchy-Schwarz Quadratic Mutual Information,” in IEEE International Conference on Robotics and Automation, 2015.
  * [6] B. Charrow, G. Kahn, S. Patil, S. Liu, K. Goldberg, P. Abbeel, N. Michael, and V. Kumar, “Information-Theoretic Planning with Trajectory Optimization for Dense 3D Mapping,” in Robotics: Science and Systems, vol. 6, 2015.
  * [7] Z. Zhang, T. Henderson, V. Sze, and S. Karaman, “FSMI: Fast computation of Shannon Mutual Information for Information Theoretic Mapping,” in IEEE International Conference on Robotics and Automation, 2019.
  * [8] T. Henderson, V. Sze, and S. Karaman, “An efficient and continuous approach to information-theoretic exploration,” in 2020 IEEE International Conference on Robotics and Automation (ICRA), pp. 8566–8572, IEEE, 2020.
  * [9] N. Atay and B. Bayazit, “A motion planning processor on reconfigurable hardware,” in Robotics and Automation, 2006. ICRA 2006. Proceedings 2006 IEEE International Conference on, pp. 125–132, 2006.
  * [10] S. Murray, W. Floyd-Jones, Y. Qi, D. J. Sorin, and G. Konidaris, “Robot Motion Planning on a Chip,” in Robotics: Science and Systems, 2016.
  * [11] Y. Kim, D. Shin, J. Lee, and H.-J. Yoo, “BRAIN: A Low-Power Deep Search Engine for Autonomous Robots,” IEEE Micro, vol. 37, no. 5, pp. 11–19, 2017.
  * [12] S. Xiao, N. Bergmann, and A. Postula, “Parallel RRT* architecture design for motion planning,” in Field Programmable Logic and Applications (FPL), 2017 27th International Conference on, pp. 1–4, 2017.
  * [13] Z. Zhang, A. A. Suleiman, L. Carlone, V. Sze, and S. Karaman, “Visual-Inertial Odometry on Chip: An Algorithm-and-Hardware Co-design Approach,” in Robotics: Science and Systems, 2017.
  * [14] A. Suleiman, Z. Zhang, L. Carlone, S. Karaman, and V. Sze, “Navion: A Fully Integrated Energy-Efficient Visual-Inertial Odometry Accelerator for Autonomous Navigation of Nano Drones,” in 2018 IEEE Symposium on VLSI Circuits, pp. 133–134, 2018.
  * [15] P. Z. X. Li, Z. Zhang, S. Karaman, and V. Sze, “High-throughput Computation of Shannon Mutual Information on Chip,” in Robotics: Science and Systems (RSS), 2019.
  * [16] A. Elfes, “Using occupancy grids for mobile robot perception and navigation,” Computer, no. 6, pp. 46–57, 1989.
  * [17] B. Yamauchi, “A frontier-based approach for autonomous exploration,” in IEEE International Symposium on Computational Intelligence in Robotics and Automation, pp. 146–151, 1997.
  * [18] B. J. Julian, S. Karaman, and D. Rus, “On mutual information-based control of range sensing robots for mapping applications,” The International Journal of Robotics Research, vol. 33, no. 10, pp. 1375–1392, 2014.
  * [19] F. Bourgault, A. A. Makarenko, S. B. Williams, B. Grocholsky, and H. F. Durrant-Whyte, “Information based adaptive robotic exploration,” in Intelligent Robots and Systems, 2002. IEEE/RSJ International Conference on, vol. 1, pp. 540–545, IEEE, 2002.
  * [20] G. M. Hoffmann and C. J. Tomlin, “Mobile sensor network control using mutual information methods and particle filters,” IEEE Transactions on Automatic Control, vol. 55, no. 1, pp. 32–47, 2010.
  * [21] T. Kollar and N. Roy, “Efficient Optimization of Information-Theoretic Exploration in SLAM,” in AAAI, vol. 8, pp. 1369–1375, 2008.
  * [22] C. Stachniss, G. Grisetti, and W. Burgard, “Information gain-based exploration using rao-blackwellized particle filters.,” in Robotics: Science and Systems, vol. 2, pp. 65–72, 2005.
  * [23] J. Bresenham, “Algorithm for computer control of a digital plotter,” IBM Syst. J., vol. 4, pp. 25–30, 1965.
  * [24] S. A. Edwards, “Bresenham’s Line Algorithm in Hardware,” 2012.  URL: http://www.cs.columbia.edu/~sedwards/classes/2012/4840/lines.pdf.
  * [25] H. H. González-Baños and J.-C. Latombe, “Navigation strategies for exploring indoor environments,” The International Journal of Robotics Research, vol. 21, no. 10-11, pp. 829–848, 2002.
  * [26] “Samsung s9 technical specifications: Gsmarena.” https://www.gsmarena.com/samsung_galaxy_s9-8966.php.


# Variational quantum circuits to prepare low energy symmetry states

Raja Selvarajan  Purdue University, Department of Physics and Astronomy, West
Lafayette, IN 47907, USA  Manas Sajjan  Purdue University, Department of
Chemistry, West Lafayette, IN 47907, USA  Sabre Kais  Purdue University,
Department of Physics and Astronomy, West Lafayette, IN 47907, USA  Purdue
University, Department of Chemistry, West Lafayette, IN 47907, USA  Purdue
Quantum Science and Engineering Institute, West Lafayette, IN 47907, USA
Corresponding author: Sabre Kais<EMAIL_ADDRESS>
###### Abstract

We explore how to build quantum circuits that compute the lowest energy state
corresponding to a given Hamiltonian within a Symmetry subspace by explicitly
encoding it into the circuit. We create an explicit unitary and a
variationally trained unitary that maps any vector output by ansatz
$A(\vec{\alpha})$ from a defined subspace to a vector in the symmetry space.
The parameters are trained varitionally to minimize the energy thus keeping
the output within the labelled symmetry value. The method was tested for a
spin XXZ hamiltonian using rotation and reflection symmetry and $H_{2}$
hamiltonian within $S_{z}=0$ subspace using $S^{2}$ symmetry. We have found
the variationally trained unitary surprisingly giving very good results with
very low depth circuits and can thus be used to prepare symmetry states within
near term quantum computers.

##  1 Introduction

One of the most important problems quantum computers have been envisioned to
solve is the simulation of Hamiltonian dynamics and computation of ground
state energies. Phase estimation algorithm [1, 2] despite solving this problem
requires circuits that run deep, something that cannot be afforded for within
the NISQ [3] era. Alternative algorithms based on hybrid models that make use
of variational methods, for instance variational quantum eigensolvers (VQE)
[4] and quantum imaginary time evolution [5] have been found to be more
resilient to noisy quantum devices [6]. The strength of a variational method
depends a lot on the variational ansatz that has been used for the simulation
of the state. Unitary coupled clusters [7] decomposed using trotterization
[8], RBM ansatz [9] based models and using hardware efficient gates to create
generalized layered models [10] are some of the circuit designs that have been
studied under variational methods.

Solving for the low energy states that lie within a Symmetry subspace ,i.e,
labelled by a specified symmetry value, forms a sub class within the
generalized constrained optimization problems, where the idea is to minimize a
cost function subject to a given set of constraints. In case where the
operators of the cost function commute with that of the constraints, one could
straightforwardly penalize the cost function with an additional term that
captures the error in symmetry value [11]. Alternatively one could design
circuits that variationally only explore the symmetry subspace by defining the
circuit using well defined structured features. These methods explore a
smaller Hilbert space for optimization and are likely to converge faster.
Barkoutsos et al [12] used a particle conserving gate alongside a particle
hole conserving representation to produce ground states of simple molecular
systems. Gard et al introduced a systematic way of preserving symmetry
subspaces for particle number, total spin, spin projection and time reversal
[13].

Unlike previous studies which aimed to provide algorithms for very specific
symmetries, we would like to demonstrate the efficacy of two new techniques
that can tackle any symmetry generically. We further compute the ground state
energy of a given Hamiltonian constrained to lie within the specified symmetry
subspace, i.e, the state is a eigenvector of the chosen symmetry operator with
a user-defined eigenvalue.

The organization of the paper is as follows. In Section 2 we illustrate in
detail the underlying theoretical framework of both the methods. In Section 3
we discuss the results using Heisenberg XXZ-Spin Hamiltonian and choose two
symmetry operators pertinent to the system. We also explore a real molecular
system ($H_{2}$) to filter singlet states of total spin angular momentum
squared ($S^{2}$) as the corresponding symmetry operator. We conclude
thereafter in Section 5 with a brief discussion of possible future extensions.

##  2 Method

Hybrid variational quantum algorithms have been extensively studied in the
context of solving unconstrained [14, 15, 16, 17] and even constrained
optimization problems [18, 19, 11]. The primary workhorse of such methods
revolve around iteratively minimizing a cost function through a usual
gradient-based optimization scheme to tweak the parameters of the quantum
circuit subsequently. The gradients of the cost-function are typically
computed directly from the quantum circuit [20], for instance using a
parameter shift method [21], while the succeeding parameter updates proceeds
classically. Updates are expressed as expectation values of output states
against some hermitian operator that can be re-expressed as a pauli string sum
and computed independently using measurements in the basis that diagonalizes
the operator or using the hadamard test [22]. The output state is typically
representative of an ansatz that is expressive enough to explore a sizeable
portion of the whole Hilbert space of dimension that scales as $2^{n}$ where
$n$ is the number of qubits used in the variational method. The ansatz so
chosen is generic and oblivious to the symmetry eigensector being sought.

In this work, we explore an alternative route. The cornerstone of our
technique lies in its ability to confine the state-ansatz apriori to a
particular eigenspace of the chosen symmetry operator even before the
optimization for minimal energy is attempted. This is attained by building
into the ansatz a variational state that selectively explores the symmetry-
subspace of interest only.

We now define the problem formally which shall be attempted to be solved in
this work. Let us consider a system characterized by a hamiltonian $H$ and let
the complete set of symmetry operators for this system be a set
$P=\\{O_{k}\\}_{k=1}^{N}$ wherein $[H,O_{m}]=0\>\>\forall\>\>O_{m}\in P$. Thus
the set $H\bigcup P$ completely characterizes the state-space of the physical
system. We select an operator $O_{k}\>\>\in P$ based on a user-defined choice
and find the target state $|\psi^{*}\rangle$ through the following
optimization scheme.

$\displaystyle|\psi^{*}\rangle$ $\displaystyle=$
$\displaystyle\underset{\psi(\vec{\theta})\>\in\>\Omega}{argmin_{\psi(\vec{\theta})}}\langle\psi(\vec{\theta})|H|\psi(\vec{\theta})\rangle$
$\displaystyle such\>that\>\>\>\>\Omega$ $\displaystyle=$
$\displaystyle\\{|x\rangle|\>\>O_{k}|x\rangle=S|x\rangle,\>\>\forall\>\>|x\rangle\in\mathbb{C}^{2^{n}}\\}$
(1)

where $S$ is the user-specified eigenvalue of $O_{k}$ which labels the desired
eigenspace and $\vec{\theta}$ are the variational parameters. Without loss of
generality one can envision that $|\Omega|=k$ i.e. the subspace is
$k$-dimensional.

We start by first constructing a generic $k$-dimensional variational ansatz
$A(\vec{\alpha})|0\rangle^{\bigotimes n}$ which allows us to subsequently
confine the acceptable state within the symmetry subspace $\Omega$ defined in
Eq.2, where $A(\vec{\alpha})$ is the unitary creating the variational ansatz
with $\vec{\alpha}$ being the corresponding parameters. If $k$ is expressible
as an exponentiation of 2, one could condition the variational circuit over
the first $n-m$ qubits, where $m=log_{2}(k)$. Otherwise we choose $m=\lfloor
log_{2}(k)\rfloor$ and then use a permutation gate to swap in the left over
basis states before applying another layer of variational ansatz.
Mathematically the action of the unitary $A(\vec{\alpha})$ can be envisioned
as

$A(\vec{\alpha})|0\rangle^{\bigotimes
n}=\sum_{i=1}^{2^{n}}a_{i}(\vec{\alpha})|i\rangle=\sum_{i=2^{n}-k}^{2^{n}}a_{i}(\vec{\alpha})|i\rangle$
(2)

wherein $|i\rangle$ denotes the computational basis states and
$a_{i}(\vec{\alpha})$ are the respective coefficients parameterized over
$\vec{\alpha}$. The last equality in Eq. 2 follows from the fact that
$a_{i}(\vec{\alpha})=0\>\>\forall\>\>i\>\>\leq 2^{n}-k$. This allows us to
create a variational ansatz over a $k$-dimensional subspace of
$\mathbb{C}^{2^{n}}$ as is required. Fig 1 provides a schematic view of the
circuit involved in construction of the ansatz $A(\alpha)$. We now present two
different methods that differs in the post-processing of the ansatz
$A(\vec{\alpha})|0\rangle^{\bigotimes n}$ once created. As mentioned before,
the ultimate goal of both the method will be to map the computational basis
states $\\{|i\rangle\\}_{i=2^{n}-k}^{2^{n}}$ onto the symmetry sub-space
$\Omega$ using a unitary transformation and then obtain the minimal energy
eigenstate in that subspace.

Figure 1: Ansatz ($A(\alpha)$) that has been used to build the state within a
k-dimensional space in Eq. 2

###  2.1 Method 1- Exact Unitary

This method uses an exact unitary that is composed of the eigenvectors of the
symmetry operator. Let $U$ be the unitary operator that diagonalizes $O_{k}$
i.e. $U=\sum_{i=1}^{2^{n}}|\phi_{i}\rangle\langle i|$ wherein
$|\phi_{i}\rangle$ are the eigenvectors of the symmetry operator $O_{k}$.
Using the unitary $U$ on the k-dimensional ansatz defined in Eq.2, we thus
get,

$UA(\vec{\alpha})\ket{0}^{\bigotimes
n}=\sum_{j=1}^{2^{n}}|\phi_{j}\rangle\langle
j|(\sum_{i=2^{n}-k}^{2^{n}}a_{i}(\vec{\alpha})\ket{i})=\sum_{i=2^{n}-k}^{2^{n}}a_{i}(\vec{\alpha})\ket{\phi_{i}}$
(3)

where in the last equation in Eq.3 only the eigenvectors
$|\phi_{i}\rangle$${}_{i=2^{n}-k}^{2^{n}}$ $\in\Omega$ survives. Operationally
the matrix $U$ is constructed by stacking the $k$ eigenvectors corresponding
to eigenvalue $S$ in the last $k$-columns. Note that the coefficients $a_{i}$
in Eq.3 are explicitly dependant on tunable parameters $\vec{\alpha}$ imparted
from the ansatz $A(\vec{\alpha})$. We thereafter train these parameters
$\vec{\alpha}$ by minimizing the energy of the output state using the
hamiltonian $H$ of the system as follows,

$\vec{\alpha}^{*}=argmin_{\vec{\alpha}}\>\>^{n\bigotimes}\bra{0}A^{\dagger}(\vec{\alpha})U^{\dagger}HUA(\vec{\alpha})\ket{0}^{\bigotimes
n}$ (4)

The minimization scheme in Eq.4 leads us to the minimal energy state within
the symmetry subspace as is required. As this method uses the exact unitary
$U$, it shall work irrespective of any initial ansatz with an output that is
restricted to the subspace. We describe schematically method 1 in Fig. 2(a).
This method however requires decomposing the matrix $U$ into a sequence of
elementary gates which depending on the symmetry may result in a quantum
circuit with very large depth. To circumvent this drawback, we propose an
alternative approach (Method 2) described next.

Figure 2: The circuit decomposition and essential steps in (a) Method 1 and
(b) Method 2

###  2.2 Method 2- Approximate Unitary construction

In this method, instead of using the exact unitary $U$ obtained from the
diagonalization of the symmetry operator $O_{k}$ as introduced in Method 1, we
approximate it using a parameterized ansatz $\tilde{U}$ which allows for a low
depth quantum circuit construction. Operationally we use the same ansatz as in
$A(\vec{\alpha})$ for constructing $\tilde{U}(\vec{\theta})$ and then learn
the parameters variationally by minimizing the cost function
$\langle(O_{k}-S)^{2}\rangle$ where $S$ is the desired eigenvalue. The state
on which the unitary $\tilde{U}(\vec{\theta})$ acts is
$A(\vec{\alpha})|0\rangle^{\bigotimes n}$ as introduced in the previous
section. The variational minimization over $\vec{\theta}$ is done over many
realization of the parameter set $\vec{\alpha}$ which is akin to an averaging
procedure over the underlying sampling distribution $p(\vec{\alpha})$.
Mathematically the parameter set $\vec{\theta}$ is learned as follows:

$\vec{\theta}^{*}=argmin_{\vec{\theta}}\>\>\langle^{n\bigotimes}\bra{0}A^{\dagger}(\vec{\alpha})\tilde{U}^{\dagger}(\vec{\theta})(O_{k}-S)^{2}\tilde{U}(\vec{\theta})A(\vec{\alpha})\ket{0}^{\bigotimes
n}\rangle_{p(\vec{\alpha})}$ (5)

where $\langle\rangle_{p(\vec{\alpha})}$ represents the averaging over the
distribution $p(\vec{\alpha})$. The parameters $\vec{\theta}^{*}$ are trained
so as to achieve a very low margin of error of allowing $\tilde{U}$ to
faithfully mimic the exact unitary $U$ and confine any subsequent operation to
the symmetry subspace $\Omega$ irrespective of the state prepared by the
ansatz $A$. With the parameters $\vec{\theta}$ known, one can proceed towards
doing a variational optimization on the parameters $\vec{\alpha}$ just as
before so as to minimize the energy according to the

$\vec{\alpha}^{*}=argmin_{\vec{\alpha}}\>\>^{n\bigotimes}\bra{0}A^{\dagger}(\vec{\alpha})\tilde{U}^{\dagger}(\vec{\theta^{*}})H\tilde{U}(\vec{\theta^{*}})A(\vec{\alpha})\ket{0}^{\bigotimes
n}$ (6)

We describe schematically method 2 in Fig. 2(b). In Fig.3 we describe entire
algorithm that is followed for both method 1 and method 2 and highlight the
essential steps. Since we have discussed in detail the first two key steps of
the algorithm which includes preparation of the ansatz using
$A(\vec{\alpha})\ket{0}^{\bigotimes n}$ (see Fig.1) and then construction and
operation by the unitary which confines the state within a symmetry subspace
(see Fig.2), we emphasize in somewhat detail the subsequent steps which
articulates the procedure for obtaining the ground state energy in Fig.3
(designated as (III)). The procedure is similar to any variational hybrid
algorithms routinely used nowadays. One computes the average energy using the
Hamiltonian of the system and the parameterized state constrained in the
symmetry sub-space (from step II in Fig.3). The gradient of this energy with
respect to the parameters of the state is constructed and the convergence of
the norm of the gradient is checked. If the desired threshold is not attained,
parameters of the state is further changed (by changing $\vec{\alpha}$ in step
I in Fig. 3) for the next iteration. The unitaries necessary to accomplish
this operation comes from the standard procedure of the conversion of the
system Hamiltonian into Pauli-strings. The accompanying circuit description of
such unitaries would then be highly system specific

Figure 3: Flowchart indicating the steps involved in the algorithm . (I)
Indicates preparation of the ansatz as illustrated in Fig. 1. (II) Preparation
of the unitary that dictates symmetry restriction as discussed in Fig.2. (III)
The key steps in the variational optimization of the energy.

##  3 Results

Both methods discussed above have been tested against $XXZ$ spin Hamiltonian
and $H_{2}$ molecule hamiltonian. We plot against chosen symmetries, the
energy error and state fidelity against the low energy state that respects the
symmetry. For method 2 that makes use of an approximate Unitary construction,
we show in addition, the deviation in the expected symmetry of output state
from the chosen value. The results have been obtained using Qiskit [23] state-
vector simulator. The parameter updates in every iteration has been globally
bound by 0.5 radians and reduced iteratively for convergence. The use of only
$R_{y}$ and $R_{x}$ gates in our ansatz allows for gradients to be calculated
with a shift of $\pi$ radians on the respective gate parameters.

###  3.1 $XXZ$ Spin Hamiltonian

The $XXZ$ spin hamiltonian is given by,

$H=\sum_{i}J(\sigma^{x}_{i}\sigma^{x}_{i+1}+\sigma^{y}_{i}\sigma^{y}_{i+1})+K\sigma^{z}_{i}\sigma^{z}_{i+1}$
(7)

where the summation is over the nearest neighbour spins with open boundaries.
We study the low energy states respected by the Reflection and Rotation
Symmetry [24].

Both symmetry operators have eigenvalues of $\pm 1$. We would like to
variationally probe the low energy states in each of these subspaces. We set
the number of spins to be $4$ and work with $J=1$ and $K=3$. This corresponds
to $16$ dimensional Hilbert space. For method 2 we have made use of D=5 layers
as shown in Fig 2 to train the unitary $\tilde{U}(\vec{\theta})$ up to a mean
error of 0.001 on the Symmetry value for 100 random samples generated by the
ansatz $A(\vec{\alpha})$

####  3.1.1 Reflection Symmetry

The reflection symmetry operator is given by,
$S=\Pi_{i}(\sigma^{x}_{i}\sigma^{x}_{i+1}+\sigma^{y}_{i}\sigma^{y}_{i+1}+\sigma^{z}_{i}\sigma^{z}_{i+1}+I)/2$.
The reflection operator has a 6 dimensional subspace that measures -1 and 10
dimensional subspace that measures +1. Fig 4 shows energy error and state
fidelity plots indicating the convergence of the variational methods to the
exact solution for each of the Symmetry values. Note that in method 2, for
$S=-1$, the method converges at the correct solution, despite the overall
ground state energy being much lower. We notice that the symmetry value
fluctuates within a very small interval around the exact value during the
training.

(a) Method 1: Energy error vs iterations

(b) Method 1: State fidelity vs iterations

(c) Method 2: Energy error vs iterations

(d) Method 2: State fidelity vs iterations

(e) Method 2: Symmetry value vs iterations

(f) Method 2: Symmetry value vs iterations

Figure 4: Computing the lowest energy state of XXZ Hamiltonian within subspace
labelled by the Reflection Symmetry Operator using Unitary constructed from
method 1 and method 2

####  3.1.2 Rotation Symmetry

The reflection symmetry operator is given by, $S=\Pi_{i}(\sigma^{x}_{i})$. The
symmetry subspace of the rotation operator with eigenvalue $\pm 1$ is 8
dimensional. Fig 5 shows energy error and state fidelity plots indicating the
convergence of the variational methods to the exact solution for each of the
Symmetry values. Note that in method 2, for $S=-1$, the method converges at
the right solution, despite the overall ground state energy being much lower.
We notice that the symmetry value fluctuates within a very small interval
around the exact value during the training.

(a) Method 1: Energy error vs iterations

(b) Method 1: State fidelity vs iterations

(c) Method 2: Energy error vs iterations

(d) Method 2: State fidelity vs iterations

(e) Method 2: Symmetry value vs iterations

(f) Method 2: Symmetry value vs iterations

Figure 5: Computing the lowest energy state of XXZ Hamiltonian within subspace
labelled by the Rotation Symmetry Operator using Unitary constructed from
method 1 and method 2

###  3.2 $H_{2}$ Hamiltonian

As an example of a real-molecular system we choose the prototypical $H_{2}$
molecule. The electronic Hamiltonian for the system is constructed in the
$STO-3G$ basis in the subspace spanned by vectors of the form
—$n_{2}(\downarrow),n_{1}(\downarrow),n_{2}(\uparrow),n_{1}(\uparrow)\rangle$
with $\sum_{i}S_{z}^{i}=0$ and $N=\sum_{i}a^{\dagger}_{i}a_{i}=2$ (the four
such states are $|0110\rangle,|0101\rangle,|1010\rangle,|1001\rangle$). The
corresponding matrix (in units of eV) at the equilibrium bond length of 0.725
Å is [25]:

$H=\begin{pmatrix}-1.06&0&0&0.18\\\ 0&-1.84&0.18&0\\\ 0&0.18&-0.23&0\\\
0.18&0&0&-1.06\end{pmatrix}$ (8)

Since the entire space is spanned by states with the zero eigenvalue of
$\sum_{i}S_{z}^{i}$, we have used the total angular momentum squared as the
symmetry operator of our choice ($O_{k}=\hat{S}^{2}$). The latter operator in
the aforesaid basis

$S^{2}=\frac{1}{2}\begin{pmatrix}1.0&0&0&-1.0\\\ 0&0&0&0\\\ 0&0&0&0\\\
-1.0&0&0&1.0\end{pmatrix}$ (9)

Only one of the eigenvalues of the $S^{2}=s(s+1)$ matrix is 1 and the
remaining ones are all 0. The techniques developed in this report requires the
symmetry subspace to be of dimension greater than 1. Thus we shall be
restricting to the subspace $S^{2}=s(s+1)=0$. For method 2 we have made use of
D=4 layers as shown in Fig 2 to train the unitary $U(\vec{\theta})$ up to a
mean error of 0.001 on the Symmetry value for 100 random samples generated by
the ansatz $A(\vec{\alpha})$. Fig 6 shows energy error and state fidelity
plots indicating the convergence of the variational methods to the exact
solution for $S^{2}=s(s+1)=0$. Here again, we notice that the symmetry value
fluctuates within a very small interval around the exact value during the
training.

(a) Method 1: Energy error vs iterations

(b) Method 1: State fidelity vs iterations

(c) Method 2: Energy error vs iterations

(d) Method 2: State fidelity vs iterations

(e) Method 2: Symmetry value vs iterations

Figure 6: Computing the lowest energy state of the hydrogen Hamiltonian within
$\sum_{i}S^{i}_{z}=0$ subspace labelled by $S^{2}=s(s+1)=0$ using Unitary
constructed from method 1 and method 2

##  4 Discussion

We discuss two general methods of restricting the exploration of a symmetry
labelled subspace for a variational ansatz in identifying the ground state
energy of a spin $XXZ$ Hamiltonian and $H_{2}$ Hamiltonian. Method 1 makes use
of an exact Unitary constructed from the Symmetry operator giving raise to a
large circuit decomposition. This is overcome through the use of variationally
trained Unitaries developed in method 2. The methods developed in here work
more generally where the constraint function need not commute with
Hamiltonian,i.e, need not be a symmetry, as reflected in construction of the
unitary $\tilde{U}(\vec{\theta})$. We would also like to point out that method
2 provides for a useful tool in developing Unitaries variationally that
satisfy a property of interest. The depth of the circuit and thus the number
of parameters to be trained is likely to increase with the number of qubits.
As most symmetry operators of interest are usually defined with the similar
operators acting over the entire qubit space, as a future work one might want
to investigate about minimal 2 qubit unitaries that leave the Symmetry value
unchanged and use them to create generalized ansatz over several qubits.

## References

  * [1] Daniel S. Abrams and Seth Lloyd.  Quantum algorithm providing exponential speed increase for finding eigenvalues and eigenvectors.  Physical Review Letters, 83(24):5162–5165, Dec 1999.
  * [2] M.A. Nielsen, I.L. Chuang, and I.L. Chuang.  Quantum Computation and Quantum Information.  Cambridge Series on Information and the Natural Sciences. Cambridge University Press, 2000.
  * [3] John Preskill.  Quantum computing in the nisq era and beyond.  Quantum, 2:79, Aug 2018.
  * [4] Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, Man-Hong Yung, Xiao-Qi Zhou, Peter J. Love, Alán Aspuru-Guzik, and Jeremy L. O’Brien.  A variational eigenvalue solver on a photonic quantum processor.  Nature Communications, 5(1), Jul 2014.
  * [5] Raja Selvarajan, Vivek Dixit, Xingshan Cui, Travis Humble, and Sabre Kais.  Prime factorization using quantum variational imaginary time evolution.  Scientific Reports, 11:20835, 10 2021.
  * [6] P. J. J. O’Malley, R. Babbush, I. D. Kivlichan, J. Romero, J. R. McClean, R. Barends, J. Kelly, P. Roushan, A. Tranter, N. Ding, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. G. Fowler, E. Jeffrey, E. Lucero, A. Megrant, J. Y. Mutus, M. Neeley, C. Neill, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. C. White, P. V. Coveney, P. J. Love, H. Neven, A. Aspuru-Guzik, and J. M. Martinis.  Scalable quantum simulation of molecular energies.  Phys. Rev. X, 6:031007, Jul 2016.
  * [7] Joonho Lee, William J. Huggins, Martin Head-Gordon, and K. Birgitta Whaley.  Generalized unitary coupled cluster wave functions for quantum computation.  Journal of Chemical Theory and Computation, 15(1):311–324, Nov 2018.
  * [8] Yuan Su, Hsin-Yuan Huang, and Earl T. Campbell.  Nearly tight trotterization of interacting electrons.  Quantum, 5:495, Jul 2021.
  * [9] Rongxin Xia and Sabre Kais.  Quantum machine learning for electronic structure calculations.  Nature Communications, 9, 10 2018.
  * [10] Abhinav Kandala, Antonio Mezzacapo, Kristan Temme, Maika Takita, Markus Brink, Jerry M. Chow, and Jay M. Gambetta.  Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets.  Nature, 549(7671):242–246, Sep 2017.
  * [11] Ilya G. Ryabinkin, Scott N. Genin, and Artur F. Izmaylov.  Constrained variational quantum eigensolver: Quantum computer search engine in the fock space, 2018.
  * [12] Panagiotis Barkoutsos, Jérôme Gonthier, Igor Sokolov, Nikolaj Moll, Gian Salis, Andreas Fuhrer, Marc Ganzhorn, Daniel Egger, Matthias Troyer, Antonio Mezzacapo, Stefan Filipp, and Ivano Tavernelli.  Quantum algorithms for electronic structure calculations: particle/hole hamiltonian and optimized wavefunction expansions.  05 2018.
  * [13] Bryan T. Gard, Linghua Zhu, George S. Barron, Nicholas J. Mayhall, Sophia E. Economou, and Edwin Barnes.  Efficient symmetry-preserving state preparation circuits for the variational quantum eigensolver algorithm.  npj Quantum Information, 6(1), Jan 2020.
  * [14] Yunseong Nam, Jwo-Sy Chen, Neal C. Pisenti, Kenneth Wright, Conor Delaney, Dmitri Maslov, Kenneth R. Brown, Stewart Allen, Jason M. Amini, Joel Apisdorf, Kristin M. Beck, Aleksey Blinov, Vandiver Chaplin, Mika Chmielewski, Coleman Collins, Shantanu Debnath, Andrew M. Ducore, Kai M. Hudek, Matthew Keesan, Sarah M. Kreikemeier, Jonathan Mizrahi, Phil Solomon, Mike Williams, Jaime David Wong-Campos, Christopher Monroe, and Jungsang Kim.  Ground-state energy estimation of the water molecule on a trapped ion quantum computer, 2019.
  * [15] J. I. Colless, V. V. Ramasesh, D. Dahlen, M. S. Blok, M. E. Kimchi-Schwartz, J. R. McClean, J. Carter, W. A. de Jong, and I. Siddiqi.  Computation of molecular spectra on a quantum processor with an error-resilient algorithm.  Phys. Rev. X, 8:011021, Feb 2018.
  * [16] Alexander J. McCaskey, Zachary P. Parks, Jacek Jakowski, Shirley V. Moore, T. Morris, Travis S. Humble, and Raphael C. Pooser.  Quantum chemistry as a benchmark for near-term quantum computers, 2019\.
  * [17] Yangchao Shen, Xiang Zhang, Shuaining Zhang, Jing-Ning Zhang, Man-Hong Yung, and Kihwan Kim.  Quantum implementation of the unitary coupled cluster for simulating molecular electronic structure.  Physical Review A, 95(2), Feb 2017.
  * [18] Kohdai Kuroiwa and Yuya O. Nakagawa.  Penalty methods for a variational quantum eigensolver.  Phys. Rev. Research, 3:013197, Feb 2021.
  * [19] Manas Sajjan, Shree Hari Sureshbabu, and Sabre Kais.  Quantum machine-learning for eigenstate filtration in two-dimensional materials.  Journal of the American Chemical Society, 143(44), 10 2021.
  * [20] Maria Schuld, Ville Bergholm, Christian Gogolin, Josh Izaac, and Nathan Killoran.  Evaluating analytic gradients on quantum hardware.  Physical Review A, 99(3), Mar 2019.
  * [21] K. Mitarai, M. Negoro, M. Kitagawa, and K. Fujii.  Quantum circuit learning.  Physical Review A, 98(3), Sep 2018.
  * [22] Kosuke Mitarai and Keisuke Fujii.  Methodology for replacing indirect measurements with direct measurements.  Phys. Rev. Research, 1:013006, Aug 2019.
  * [23] Héctor Abraham and AduOffei.  Qiskit: An open-source framework for quantum computing, 2019.
  * [24] Kira Joel, Davida Kollmar, and Lea F. Santos.  An introduction to the spectrum, symmetries, and dynamics of spin-1/2 heisenberg chains.  American Journal of Physics, 81(6):450–457, Jun 2013.
  * [25] Feng Zhang, Niladri Gomes, Noah F. Berthusen, Peter P. Orth, Cai-Zhuang Wang, Kai-Ming Ho, and Yong-Xin Yao.  Shallow-circuit variational quantum eigensolver based on symmetry-inspired hilbert space partitioning for quantum chemical calculations.  Physical Review Research, 3(1), Jan 2021.


# SelfRemaster: Self-Supervised Speech Restoration with
Analysis-by-Synthesis Approach Using Channel Modeling

###### Abstract

We present a self-supervised speech restoration method without paired speech
corpora. Because the previous general speech restoration method uses
artificial paired data created by applying various distortions to high-quality
speech corpora, it cannot sufficiently represent acoustic distortions of real
data, limiting the applicability. Our model consists of analysis, synthesis,
and channel modules that simulate the recording process of degraded speech and
is trained with real degraded speech data in a self-supervised manner. The
analysis module extracts distortionless speech features and distortion
features from degraded speech, while the synthesis module synthesizes the
restored speech waveform, and the channel module adds distortions to the
speech waveform. Our model also enables audio effect transfer, in which only
acoustic distortions are extracted from degraded speech and added to arbitrary
high-quality audio. Experimental evaluations with both simulated and real data
show that our method achieves significantly higher-quality speech restoration
than the previous supervised method, suggesting its applicability to real
degraded speech materials.

Index Terms: Speech restoration, self-supervised learning, neural vocoder,
historical speech

##  1 Introduction

Speech data used in modern speech processing research is usually recorded with
sophisticated digital recording equipment and in a well-designed environment.
However, most existing audio data contains acoustic distortions and have
different data properties from high-quality audio data. To expand the range of
applications of speech processing, it is necessary to make use of speech data
in varying qualities, e.g., speech material recorded with poor or outdated
audio devices [1, 2] and endangered languages’ voices [3, 4, 5]. However,
restoring or analyzing such degraded audio is challenging because paired high-
quality and degraded speech data is not generally available. The previous
general speech restoration method [6] relies on artificial paired corpora for
supervised learning. It cannot sufficiently represent acoustic distortions of
real data, which limits the applicability.

We propose a self-supervised speech restoration method without paired speech
data. Our proposed model consists of analysis, synthesis, and channel modules
based on the recording process of degraded speech. The analysis module
extracts time-variant distortionless speech features and time-invariant
acoustic distortion features from degraded speech (hereafter, “channel
feature”). The other modules carry out the generation process of distorted
speech; the synthesis module synthesizes distortionless speech (i.e., restored
speech), and the channel module adds acoustic distortion. The model is trained
in a self-supervised manner by minimizing the reconstruction loss between
degraded input and reconstructed speech. We also propose a dual-learning
method for stable self-supervised learning, where a backward process defines
distortionless speech feature loss with arbitrary high-quality speech corpora
while the forward process calculates the reconstruction loss. Experimental
evaluations show that our method achieves significantly higher-quality speech
restoration than the previous supervised method. Our
implementation111https://github.com/Takaaki-Saeki/ssl_speech_restoration and
audio samples222https://takaaki-saeki.github.io/ssl_remaster_demo/ are
publicly available. The contributions of this paper are as follows:

  * •

We propose self-supervised speech restoration that can learn various acoustic
distortions without a paired speech corpora and is applicable to real degraded
speech data.

  * •

Our method achieved significantly higher-quality speech restoration than the
previous supervised method on both simulated data and real historical
materials.

  * •

Our model enables audio effect transfer, in which only channel features are
extracted from degraded speech and applied to other high-quality speech.

##  2 Related work

Speech restoration. Many previous studies have focused on specific speech
restoration tasks such as bandwidth extension [7, 8], dereverberation [9, 10],
denoising [11, 12], and declipping [13, 14]. Similar to our work, some
previous studies [15, 16] jointly perform multiple speech restoration tasks,
and some methods leverage waveform synthesis [17, 6] to perform speech
restoration. Unlike these studies, our work does not focus on removing
additive noises and enables self-training. A previous study on general speech
restoration with a neural vocoder [6] is the closest work to ours. The most
significant difference is that our method does not require paired data, so it
can be trained with real audio materials in a self-supervised manner instead
of artificially creating degraded speech data. Furthermore, as our proposed
method also extracts channel features at the same time, it can transfer audio
effects, which adds the desired channel distortion to arbitrary audio.

Self-supervised learning for speech. Self-supervised speech representation
learning [18, 19, 20] and methods using those models for various downstream
tasks [21, 22] have widely been studied. There are several studies related to
our work which involve self-supervised learning on speech enhancement [23, 24,
25], though they only focus on denoising while our method restores speech from
various kinds of distorted speech. Similar to our method, the DDSP autoencoder
[26] learns disentangled features from acoustic signal in a self-supervised
manner and manipulates each feature. It uses a simple sinusoidal vocoder and
assumes only reverberation as acoustic distortion, while our method leverages
a more expressive waveform synthesis model to achieve high-quality speech
restoration and channel modules to capture various acoustic distortions.

##  3 Proposed method

###  3.1 Basic framework

Figure 1: Proposed self-supervised speech restoration method, consisting of
analysis, synthesis, and channel modules that simulate recording process of
degraded speech. Basic training process (top) minimizes reconstruction loss
between degraded and reconstructed speech. We also propose dual learning (top
& bottom) with arbitrary high-quality speech corpora.

We begin with an overview of the recording process of degraded speech
underlying our method. High-quality speech (i.e., distortionless speech) is
emitted from the mouth through a human speech production process, which can be
parameterized by time-variant speech features such as traditional source-
filter vocoder features or mel spectrograms. Then acoustic distortions (e.g.,
non-linear response of recording equipment and lossy audio coding) are added
to the high-quality speech, resulting in a final recorded audio. We assume the
distortions (i.e., channel features) to be time-invariant, where the recording
equipment or audio coding does not change in each audio sample.

The proposed speech restoration model consists of analysis, synthesis, and
channel modules based on the abovementioned process, as shown in Fig. 1 (top).
All of these modules are composed of neural networks, so they can be trained
in an end-to-end manner. We denote the degraded-speech waveform and the
degraded-speech features as $\vec{x}_{\mathrm{low}}$ and
$\vec{y}_{\mathrm{low}}$, respectively. The analysis module estimates the
speech features of the restored speech, $\hat{\vec{z}}_{\mathrm{res}}$, from
$\vec{y}_{\mathrm{low}}$. At the same time, the intermediate features of the
analysis module are input to other layers in order to finally estimate the
time-invariant channel features $\vec{c}$. Note that we use mel spectrogram
for $\vec{y}_{\mathrm{low}}$ and the two features described in Section 4.1 for
$\hat{\vec{z}}_{\mathrm{res}}$. The synthesis module simulates human speech
production; it takes $\hat{\vec{z}}_{\mathrm{res}}$ as input and estimates the
speech waveform of the restored speech $\hat{\vec{x}}_{\mathrm{res}}$. We use
HiFi-GAN [27] for the synthesis module, where it is trained on an arbitrary
high-quality speech corpus and the model parameters are frozen. The channel
module simulates the acoustic distortion; it is conditioned by the channel
feature $\vec{c}$ and distorts the restored speech
$\hat{\vec{x}}_{\mathrm{res}}$ to estimate the degraded speech
$\hat{\vec{x}}_{\mathrm{low}}$. The model is trained to minimize the
reconstruction loss, which can be defined as a multi-scale spectral loss [26]
as follows:

$\mathcal{L}_{\mathrm{recons}}=\sum_{i}\\{||\vec{s}_{i}-\hat{\vec{s}}_{i}||_{1}+\alpha||\log\vec{s}_{i}-\log\hat{\vec{s}}_{i}||_{1}\\},$
(1)

where $\vec{s}_{i}$ and $\hat{\vec{s}}_{i}$ are the amplitude spectrograms of
$\vec{x}_{\mathrm{low}}$ and $\hat{\vec{x}}_{\mathrm{low}}$, respectively. The
subscript $i$ is the window length of the Fourier transform and $\alpha$ is a
weight of the log term, where we used $i=(2048,1024,512,256,128,64)$ and
$\alpha=1.0$ in our experiment. Speech restoration is performed by driving the
trained analysis and synthesis modules. The analysis module receives the
acoustic feature of degraded speech, $\vec{y}_{\mathrm{low}}$, from the
analysis module, and the synthesis module synthesizes the restored-speech
waveform $\hat{\vec{x}}_{\mathrm{res}}$.

Furthermore, the trained model can also perform an “audio effect transfer” in
addition to speech restoration. The audio effect is transferred without the
synthesis module. We first extract the channel feature $\vec{c}$ from the
degraded speech $\vec{x}_{\mathrm{low}}$ using the trained analysis module.
The trained channel module conditioned by $\vec{c}$ receives arbitrary high-
quality audio and distorts it so that the resulting audio sounds distorted
like the referred degraded speech.

Fig. 2 shows the spectrograms obtained by the proposed method. The input low-
quality speech is simulated on the basis of the quantization and resampling
described in Section 4, and is missing high-frequency bands that are present
in the ground-truth high-quality speech and has distorted low-frequency bands.
In the restored speech, the missing and distorted bands are restored.
Furthermore, the reconstructed speech output by the channel module faithfully
reproduces the input speech.

Unfortunately, the above basic self-supervised training process is unstable
because each module has a high expressive power, and the modules are only
trained with the reconstruction loss between the degraded and reconstructed
waveform. Specifically, the analysis module can represent the effect of the
channel module, so there is no guarantee that the analysis module outputs the
speech features of high-quality speech. Therefore, we propose the dual-
learning method described in Section 3.2. In Section 3.3, we also introduce
supervised pre-training to improve the performance when only a very limited
amount of data is available for the model size, such as when we use rare
historical audio.

Figure 2: Spectrograms obtained with proposed method

###  3.2 Dual learning for stable self-supervised learning

Fig. 1 shows the proposed dual-learning method. In addition to the basic
training framework described in Section 3.1, we introduce a training task that
propagates information in the reverse direction. We denote the two training
tasks in the forward direction and backward direction as
$\mathcal{T}_{\mathrm{forward}}$ and $\mathcal{T}_{\mathrm{backward}}$,
respectively. $\mathcal{T}_{\mathrm{forward}}$ is the basic learning framework
described in Section 3.1; degraded speech is sent to the analysis module to
reconstruct the degraded waveform. Meanwhile, in
$\mathcal{T}_{\mathrm{backward}}$, high-quality speech is input to the channel
module to estimate the features of high-quality speech. Our method is
analogous to the forward and backward learning between two domains in machine
translation [28].

Here, we use arbitrary high-quality speech, which does not align with the
degraded speech and only consists of other speaker’s utterances. The high-
quality speech waveform $\vec{X}_{\mathrm{high}}$ is first passed through the
channel module to generate the degraded speech waveform
$\hat{\vec{X}}_{\mathrm{low}}$ with channel features. Next, the speech
features of the degraded speech, $\hat{\vec{Y}}_{\mathrm{low}}$, are input to
the analysis module to estimate the speech feature of the restored speech,
$\hat{\vec{Z}}_{\mathrm{res}}$. Then a feature loss is defined as the mean
squared error between $\hat{\vec{Z}}_{\mathrm{res}}$ and the speech feature
$\vec{Z}_{\mathrm{high}}$ is obtained from the high-quality speech as:

$\mathcal{L}_{\mathrm{feature}}=||\vec{Z}_{\mathrm{high}}-\hat{\vec{Z}}_{\mathrm{res}}||_{2}.$
(2)

The final loss function is given as the weighted sum of the loss function
between the reconstruction loss and the feature loss as:

$\mathcal{L}=(1-\beta)*\mathcal{L}_{\mathrm{recons}}+\beta*\mathcal{L}_{\mathrm{feature}},$
(3)

where $\beta$ is a weight of the feature loss. Here, the gradient due to
feature loss is not propagated to the channel module, and the channel module
is only trained on the reconstruction loss. Intuitively,
$\mathcal{T}_{\mathrm{forward}}$ and $\mathcal{T}_{\mathrm{reverse}}$ can be
interpreted as:

  * •

$\mathcal{T}_{\mathrm{forward}}$ obtains the analysis module that adapts to
degraded speech and the channel module that expresses the acoustic distortion.

  * •

$\mathcal{T}_{\mathrm{backward}}$ obtains the analysis module to estimate
high-quality speech features.

This dual learning method enables self-supervised learning of the analysis
module and channel module without paired data.

###  3.3 Supervised pretraining for low-resource settings

Data scarcity is an obstacle in restoring historical audio materials; only a
small amount of audio is available due to the limited amount of recording
storage size available in the past. Thus, we introduce a supervised pre-
training method to obtain the initial weights for such low data resource
conditions as:

  1. 1.

We supervisedly train the analysis and channel modules using paired data
artificially created from arbitrary high-quality speech data.

  2. 2.

With the initialization from supervised pre-training, we then carry out the
self-supervised learning described in Section 3.1 or Section 3.2 with real
degraded speech data.

In the supervised pretraining, we generate pseudo-degraded speech waveform
$\vec{X}_{\mathrm{pseudo-low}}$ by randomly applying distortion to the
original high-quality speech. When we denote the restored speech feature
estimated from the pseudo-degraded speech feature $\vec{Y}_{\mathrm{pseudo-
low}}$ as $\hat{\vec{Z}}_{\mathrm{res}}$, we define the feature loss as
$\mathcal{L}_{\mathrm{feature}}=||\vec{Z}_{\mathrm{high}}-\hat{\vec{Z}}_{\mathrm{res}}||_{2}$.
Note that we input the high-quality speech waveform to the channel module
instead of using the output of the synthesis module. Then we train the model
to minimize the weighted sum of the reconstruction loss and feature loss in
the same manner as Eq. (3).

##  4 Experimental evaluations

###  4.1 Experimental condition

Database. We used both simulated and real degraded speech datasets with the
sampling rate set to 22.05 kHz. We created simulated datasets by applying
various distortions to the JSUT corpus [29], which consists of around six
hours of high-quality speech utterances of a Japanese female speaker. We
randomly selected 25 and 25 sentences for the validation and test sets,
respectively. For the real data, we used a Japanese historical speech dataset
[30]. This dataset was recorded in the 1960s–1970s and consists of nine
speakers telling folktales, and the whole duration was about 22-minutes. The
analog audio clips recorded on the cassette tape were digitized using a radio
cassette player. Because this data contains a lot of additive noise, iZotope
RX9 was used for noise reduction as preprocessing. We randomly selected two-
minute audio clips for the validation and test sets. We used the JVS corpus
[29] as a high-quality speech corpus for training the synthesis module,
supervised pretraining, and $\mathcal{T}_{\mathrm{backward}}$ of dual
learning.

Speech feature analysis. We used two types of speech features
$\hat{\vec{z}}_{\mathrm{res}}$: mel spectrogram and source-filter vocoder
features. The 80-dimensional mel spectrogram was extracted with a frame size
of 1024 and frame shift of 256. The vocoder features consists of F0 and
41-dimensional mel cepstrum coefficients with a 5-ms frame shift, which were
extracted by the WORLD vocoder [31]. Hereafter, we denote these types as
“Melspec” and “SourceFilter”, respectively.

Table 1: Evaluation results with simulated data

| (a) Band-limited | (b) Clipped | (c) Quantized & Resampled | (d) Overdrive
---|---|---|---|---
| MCD | MOS | MCD | MOS | MCD | MOS | MCD | MOS
Ground-truth | - | $4.51\pm 0.092$ | - | $4.58\pm 0.083$ | - | $4.67\pm 0.098$ | - | $4.65\pm 0.098$
Input | $17.85$ | $2.38\pm 0.115$ | $11.21$ | $2.45\pm 0.124$ | $22.25$ | $1.73\pm 0.115$ | $16.73$ | $1.54\pm 0.105$
Liu, et al. [6] | $6.98$ | $3.74\pm 0.115$ | $9.33$ | $3.01\pm 0.119$ | $9.34$ | $2.80\pm 0.109$ | $10.23$ | $2.00\pm 0.112$
SSL-dual (MelSpec) | 5.61 | $\textbf{4.20}\pm 0.090$ | 7.20 | $\textbf{3.49}\pm 0.114$ | 8.07 | $\textbf{3.27}\pm 0.111$ | 8.93 | $\textbf{2.68}\pm 0.108$
SSL-dual (SourceFilter) | $5.88$ | $3.46\pm 0.115$ | $7.59$ | $2.49\pm 0.190$ | $9.60$ | $2.66\pm 0.109$ | $9.06$ | $2.58\pm 0.125$

Model details. We used the U-Net architecture [32] for the analysis and
channel modules. In each down sampling of the U-Net, the temporal resolution
was reduced by half with four layers of residual convolution blocks and
average pooling. Each residual convolution block consists of convolutional
layers and batch normalization layers [33] with a skip connection [34]. During
up-sampling, the temporal resolution is doubled by applying deconvolution,
resulting in time-variant high-quality speech features with the same temporal
resolution as that of the input features. For more details, refer to the
published implementation. For the synthesis module, we used the pretrained
multi-speaker model333https://github.com/jik876/hifi-gan under the “MelSpec”
condition and we trained the model with the JVS corpus under the
“SourceFilter” condition.

Training settings. The batch size was four and the maximum number of epochs
was set to 50, and the model used for evaluation was selected based on the
validation results. Adam [35] was used as the optimizer with the initial
learning rate set to $0.001$. We applied a learning rate scheduling that
multiplies the learning rate by 0.5 if the validation loss did not decrease
over three epochs. In dual learning and supervised pretraining, we set $\beta$
of Eq. (3) to $0.1$ and $0.001$, respectively. For the evaluation with real
data, we used the supervised pretraining described in Section 3.3, in which
quantization and resampling were randomly applied to each speech waveform.
After $\mu$-law quantization [36] with random quantization bits ranging from 6
to 10 bits, the data was resampled with random sampling frequencies of {8,
11.25, 12, 16} kHz.

###  4.2 Evaluation on speech restoration with simulated data

We evaluated the proposed method by adding various artificial distortions to
the high-quality speech corpus described in Section 4.1. The purpose of this
evaluation is to investigate the performance in an experimental setting
similar to that of the previous supervised speech restoration model [6]. We
compared our proposed method described in Section 3.2 (“SSL-dual”) with the
previous supervised model (“Liu, et al.”) for generalized speech restoration
tasks [6], where we used a pretrained model provided by the
author444https://github.com/haoheliu/voicefixer. This supervised model is
trained with artificial paired data generated by randomly applying noise,
reverberation, clipping, and low-pass filtering to multi-speaker and multi-
lingual clean speech corpora.

We created datasets simulating two in-domain and two out-of-domain
distortions. “In-domain” and “out-of-domain” indicate whether the training
data of the supervised model includes the distortion or not. The in-domain
distortions are (a) Band-limited: Applying biquad-lowpass filtering to the
high-quality speech with a cut-off frequency of 4 kHz; and (b) Clipped:
Clipping the high-quality speech waveform with an absolute-value threshold of
0.25. The out-domain ones are more significant distortions than in-domain
ones, which are (c) Quantized & resampled: Quantizing the high-quality speech
to 8-bit $\mu$-law and resampling it to $8$ kHz; and (d) Overdrive:
Overdriving the high-quality speech using the SoX implementation [37].

We conducted an objective evaluation using mel cepstral distortion (MCD) [38]
with the dimension of the cepstrum set to 24. We also evaluated with
subjective mean opinion score (MOS) test, where we recruited 40 Japanese
speakers for each evaluation and calculated MOS with 95 % confidence
intervals. Table 1 lists the results, with the best results shown in bold.

As shown in the table, “SSL-dual (Melspec)” significantly improves the quality
of input degraded speech in all cases, indicating that it can greatly improve
the speech quality. We also verified that the proposed method significantly
outperforms the previous supervised method in both out-of-domain and in-domain
cases. Particularly, in the simple band-limited case, the proposed method
achieves a high MOS of $4.20$. In the proposed methods, the quality of
“Melspec” is higher than that of “SourceFilter”. The difference in MOS may be
due to the lower naturalness of the neural vocoder in the SourceFilter, while
MCD shows a significant improvement over the previous supervised model.

###  4.3 Evaluation on speech restoration with real data

We also evaluated the proposed method with real data. Note that there is no
ground-truth speech in this data. Because this is a small dataset, we applied
supervised pre-training as shown in Section 4.1. Table 2(a) lists the results
of the MOS test. The suffix “-pre” indicates the use of the supervised pre-
training.

The proposed method “SSL–dual–pre (Melspec),” achieved the highest score and
significantly outperformed the previous supervised model. Here, the
performance of the previous method and some of the proposed methods were lower
than the input, which is most likely due to the quality degradation in the
neural vocoder. We also determined that the MOS of both MelSpec and
SourceFilter can be improved by supervised pre-training. “SSL–dual–pre
(Melspec)” showed a higher MOS than the input but there was not significant
difference. Therefore, we additionally conducted an AB test with 40 Japanese
evaluators to clarify the difference in quality. Table 2(b) indicate that the
proposed method achieved significantly higher preference score than the input.
The above evaluation results show that our proposed method is feasible when
only a small amount of real data is available.

Table 2: Evaluation results with real data

| MOS
---|---
Input | $2.98\pm 0.136$
Liu, et al. [6] | $2.80\pm 0.138$
SSL–pre (MelSpec) | $2.70\pm 0.127$
SSL–dual (MelSpec) | $2.96\pm 0.135$
SSL–dual–pre (MelSpec) | $\textbf{3.06}\pm 0.140$
SSL–pre (SourceFilter) | $2.14\pm 0.123$
SSL–dual (SourceFilter) | $2.41\pm 0.124$
SSL–dual–pre (SourceFilter) | $2.55\pm 0.135$
(a) MOS tests of all methods

SSL–dual–pre (MelSpec) vs. Input |  $p$-value
---|---
0.565 \- 0.435 | 0.047
(b) Preference AB test

Table 3: Evaluation results of audio effect transfer | Simulated | Real
---|---|---
Reference | $3.98\pm 0.152$ | $2.99\pm 0.146$
High-quality | $1.16\pm 0.069$ | $1.30\pm 0.093$
Mean spec. diff. | $1.68\pm 0.113$ | -
Proposed | $\textbf{3.44}\pm 0.159$ | $\textbf{2.12}\pm 0.155$

###  4.4 Evaluation on audio effect transfer

To investigate the effectiveness of the audio effect transfer, we retrieved
the acoustic distortion from the training data and added it to different high-
quality speech data. There were two kinds of training data: “Simulated”
generated with quantization and resampling as in Section 4.2 and “Real” used
in Section 4.3. The high-quality speech to be modified was sampled from the
JVS corpus. We conducted a five-level similarity MOS (SMOS) test with 40
listeners to evaluate the degree to which the output speech sounded distorted
like the training data. “Reference” is a sample that has the same channel
features as the training data. For comparison, we prepared a method (“Mean
spec. diff”) that multiplies the amplitude spectrum of the high-quality speech
by the time-averaged amplitude spectrum difference of training data and high-
quality data.

The results show that the proposed method is closer to the target channel
features than the original high-quality speech on both simulated and real
data. Regarding the results on the simulated data, the SMOS is significantly
higher than that of the case corrected by the mean amplitude spectral
difference. Although the score of “Proposed” is 0.5–0.8 smaller than the
reference, it performed sufficiently.

##  5 Conclusion

We proposed a self-supervised speech restoration method consisting of
analysis, synthesis, and channel modules. Our method achieved significantly
higher-quality speech restoration than the previous supervised method in
several simulated data and was effective even with low-resource real data. Our
future work includes applying our method to a wider variety of real data, such
as where multiple speakers are speaking at the same time.

Acknowledgements: Part of this work was supported by JSPS KAKENHI Grant Number
21H04900.

## References

  * [1] J. Godfrey, E. Holliman, and J. McDaniel, “SWITCHBOARD: Telephone speech corpus for research and development,” in _Proc. ICASSP_ , 1992, pp. 517–520.
  * [2] J. H. L. Hansen, A. Sangwan, A. Joglekar, A. E. Bulut, L. Kaushik, and C. Yu, “Fearless Steps: Apollo-11 corpus advancements for speech technologies from earth to the moon,” in _Proc. INTERSPEECH_ , 2018, pp. 2758–2762.
  * [3] B. Abraham, D. Goel, D. Siddarth, K. Bali, M. Chopra, M. Choudhury, P. Joshi, P. Jyoti, S. Sitaram, and V. Seshadri, “Crowdsourcing speech data for low-resource languages from low-income workers,” in _Proc. LREC_ , 2020, pp. 2819–2826.
  * [4] S. Takamichi and H. Saruwatari, “CPJD corpus: Crowdsourced parallel speech corpus of Japanese dialects,” in _Proc. LREC_ , 2018, pp. 434–437.
  * [5] K. Matsuura, S. Ueno, M. Mimura, S. Sakai, and T. Kawahara, “Speech corpus of ainu folklore and end-to-end speech recognition for ainu language,” in _Proc. LREC_ , 2020, pp. 2622–2628.
  * [6] H. Liu, Q. Kong, Q. Tian, Y. Zhao, D. Wang, C. Huang, and Y. Wang, “Voicefixer: Toward general speech restoration with neural vocoder,” _arXiv preprint arXiv:2109.13731_ , 2021.
  * [7] J. Kontio, L. Laaksonen, and P. Alku, “Neural network-based artificial bandwidth expansion of speech,” _IEEE Transactions on Audio, Speech, and Language Processing_ , vol. 15, no. 3, pp. 873–881, 2007.
  * [8] V. Kuleshov, S. Z. Enam, and S. Ermon, “Audio super resolution using neural networks,” _arXiv preprint arXiv:1708.00853_ , 2017.
  * [9] B. Schwartz, S. Gannot, and E. A. P. Habets, “Online speech dereverberation using kalman filter and em algorithm,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 23, no. 2, pp. 394–406, 2015.
  * [10] O. Ernst, S. E. Chazan, S. Gannot, and J. Goldberger, “Speech dereverberation using fully convolutional networks,” in _Proc. EUSIPCO_ , 2018, pp. 390–394.
  * [11] Y. Luo and N. Mesgarani, “Conv-TasNet: Surpassing ideal time–frequency magnitude masking for speech separation,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 27, pp. 1256–1266, 2019.
  * [12] Y. Koizumi, K. Yatabe, M. Delcroix, Y. Masuyama, and D. Takeuchi, “Speech enhancement using self-adaptation and multi-head self-attention,” _Proc. ICASSP_ , pp. 181–185, 2020.
  * [13] F. Bie, D. Wang, J. Wang, and T. F. Zheng, “Detection and reconstruction of clipped speech for speaker recognition,” _Speech Communication_ , vol. 72, pp. 218–231, 2015.
  * [14] P. Záviška, P. Rajmic, A. Ozerov, and L. Rencker, “A survey and an extensive evaluation of popular audio declipping methods,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 15, pp. 5–24, 2021.
  * [15] K. Han, Y. Wang, D. Wang, W. S. Woods, I. Merks, and T. Zhang, “Learning spectral mapping for speech dereverberation and denoising,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 23, no. 6, pp. 982–992, 2015.
  * [16] X. Shu, Y. Zhu, Y. Chen, L. Chen, H. Liu, C. Huang, and Y. Wang, “Joint echo cancellation and noise suppression based on cascaded magnitude and complex mask estimation,” _arXiv preprint arXiv:2107.09298_ , 2021.
  * [17] A. Polyak, L. Wolf, Y. Adi, O. Kabeli, and Y. Taigman, “High fidelity speech regeneration with application to speech enhancement,” _Proc. ICASSP_ , pp. 7143–7147, 2021.
  * [18] A. Baevski, H. Zhou, A. Mohamed, and M. Auli, “wav2vec 2.0: A framework for self-supervised learning of speech representations,” _arXiv preprint arXiv:2006.11477_ , 2020.
  * [19] W.-N. Hsu, B. Bolte, Y.-H. H. Tsai, K. Lakhotia, R. Salakhutdinov, and A. Mohamed, “Hubert: Self-supervised speech representation learning by masked prediction of hidden units,” _arXiv preprint arXiv:2106.07447_ , 2021\.
  * [20] S. Chen, C. Wang, Z. Chen, Y. Wu, S. Liu, Z. Chen, J. Li, N. Kanda, T. Yoshioka, X. Xiao, J. Wu, L. Zhou, S. Ren, Y. Qian, Y. Qian, J. Wu, M. Zeng, X. Yu, and F. Wei, “Wavlm: Large-scale self-supervised pre-training for full stack speech processing,” _arXiv preprint arXiv:2110.13900_ , 2021\.
  * [21] C. Yi, J. Wang, N. Cheng, S. Zhou, and B. Xu, “Applying wav2vec2.0 to speech recognition in various low-resource languages,” _arXiv preprint arXiv:2012.12121_ , 2020.
  * [22] H. Nguyen, F. Bougares, N. A. Tomashenko, Y. Estève, and L. Besacier, “Investigating self-supervised pre-training for end-to-end speech translation,” in _Proc. INTERSPEECH_ , 2020, pp. 1466–1470.
  * [23] Y. Wang, S. Venkataramani, and P. Smaragdis, “Self-supervised learning for speech enhancement,” _arXiv preprint arXiv:2006.10388_ , 2020.
  * [24] R. E. Zezario, T. Hussain, X. Lu, H.-M. Wang, and Y. Tsao, “Self-supervised denoising autoencoder with linear regression decoder for speech enhancement,” in _Proc. ICASSP_ , 2020, pp. 6669–6673.
  * [25] Z. Du, M. Lei, J. Han, and S. Zhang, “Self-supervised adversarial multi-task learning for vocoder-based monaural speech enhancement,” in _Proc. INTERSPEECH_ , 2020, pp. 3271–3275.
  * [26] J. Engel, L. H. Hantrakul, C. Gu, and A. Roberts, “DDSP: Differentiable digital signal processing,” in _Proc. ICLR_ , 2020.
  * [27] J. Kong, J. Kim, and J. Bae, “HiFi-GAN: Generative adversarial networks for efficient and high fidelity speech synthesis,” _arXiv preprint arXiv:2010.05646_ , 2020.
  * [28] D. He, Y. Xia, T. Qin, L. Wang, N. Yu, T.-Y. Liu, and W.-Y. Ma, “Dual learning for machine translation,” in _Proc. NIPS_ , 2016, pp. 820––828.
  * [29] S. Takamichi, R. Sonobe, K. Mitsui, Y. Saito, T. Koriyama, N. Tanji, and H. Saruwatari, “JSUT and JVS: Free japanese voice corpora for accelerating speech synthesis research,” _Acoustical Science and Technology_ , vol. 41, pp. 761–768, 2020.
  * [30] T. Sasaki, _Tono no mukashibanashi / Sasaki Tokuo hen_. Ohfusha, 1985.
  * [31] M. Morise, F. Yokomori, and K. Ozawa, “WORLD: A vocoder-based high-quality speech synthesis system for real-time applications,” _IEICE transactions on information and systems_ , vol. E99-D, no. 7, pp. 1877–1884, 2016\.
  * [32] O. Ronneberger, P. Fischer, and T. Brox, “U-Net: Convolutional networks for biomedical image segmentation,” in _Proc. MICCAI_ , 2015, pp. 234–241.
  * [33] S. Ioffe and C. Szegedy, “Batch normalization: Accelerating deep network training by reducing internal covariate shift,” in _Proc ICML_ , 2015, p. 448–456.
  * [34] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” _Proc. CVPR_ , pp. 770–778, 2016.
  * [35] D. Kingma and B. Jimmy, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * [36] C. Recommendation, “Pulse code modulation (pcm) of voice frequencies,” in _ITU_ , 1988.
  * [37] C. Bagwell, et al., “SoX - Sound eXchange.” [Online]. Available: {http://sox.sourceforge.net/}
  * [38] T. Fukada, K. Tokuda, T. Kobayashi, and S. Imai, “An adaptive algorithm for mel-cepstral analysis of speech,” in _Proc. ICASSP_ , 1992, pp. 137–140.