# Multifractality in non-unitary random dynamics

Jason Iaconis  Department of Physics and Center for Theory of Quantum Matter,
University of Colorado, Boulder CO 80309, USA  Xiao Chen  Department of
Physics, Boston College, Chestnut Hill, MA 02467, USA

###### Abstract

We explore the multifractality of the steady state wave function in non-
unitary random quantum dynamics in one dimension. We focus on two classes of
random systems: the hybrid Clifford circuit model and the non-unitary free
fermion dynamics. In the hybrid Clifford model, we map the measurement driven
transition to an Anderson localization transition in an effective graph space
by using properties of the stabilizer state. We show that the volume law phase
with nonzero measurement rate is non-ergodic in the graph space and exhibits
weak multifractal behavior. We apply the same method to the hybrid Clifford
quantum automaton circuit and obtain similar multifractality in the volume law
phase. For the non-unitary random free fermion system with a critical steady
state, we compute the moments of the probability distribution of the single
particle wave function and demonstrate that it is also weakly multifractal and
has strong variations in real space.

##  1 Introduction

Non-unitary dynamics have attracted a lot of attention in the past few years.
It has been shown that for a generic many-body unitary quantum dynamics
subject to local projective measurement, there exists an entanglement phase
transition at the level of the quantum trajectories Li  _et al._ (2018a);
Skinner  _et al._ (2019); Gullans and Huse (2020a); Chan  _et al._ (2019);
Zabalo  _et al._ (2020); Gullans and Huse (2020b); Li  _et al._ (2019);
Iaconis  _et al._ (2020); Jian  _et al._ (2020); Bao  _et al._ (2020); Tang
and Zhu (2020). In the steady state wave function, by increasing the
measurement rate $p$, the entanglement entropy changes from a highly entangled
volume-law scaling to a short-range entangled area-law scaling. In particular,
when $p$ is nonzero, the volume-law phase has a non-trivial subleading
correction term which is absent in the conventional thermal phase Li  _et al._
(2019); Gullans and Huse (2020a); Fan  _et al._ (2020); Li and Fisher (2020).
The stability of this non-thermal volume law phase has interesting
interpretations in the language of quantum error correction Choi  _et al._
(2020); Gullans and Huse (2020a); Fan  _et al._ (2020); Li and Fisher (2020).

For this measurement induced phase transition, the dynamics is random in both
space and time. The randomness comes from various sources including the choice
of unitary gate, the position of the measurement gate and the measurement
outcome. The presence of the randomness leads to an emergent critical point
which is distinct from any conventional critical point in a clean
systemSkinner  _et al._ (2019); Li  _et al._ (2020). Furthermore, in the
highly entangled non-thermal volume law phase, a subleading correction to the
entanglement entropy is caused by the fluctuation of the random dynamics and
is different from the prediction given by simple mean-field theory estimation
Li and Fisher (2020); Fan  _et al._ (2020); Li  _et al._ (2021). This
observation implies the random fluctuation effect becomes dominant in the low
dimensional quantum dynamics and renders the physics significantly different
from that in the clean systems.

In this paper, we will go beyond this entanglement picture and investigate
multifractal behavior in the non-unitary random dynamics. Multifractality has
been observed in many random systems and historically has played an important
role in identifying the Anderson localization Wegner (1980); Evers and Mirlin
(2008) and spin glass phase transitions Derrida (1981). For instance, the
Anderson localization phase transition in disordered systems is a continuous
phase transition separating the extensive metallic state from a localized
state. At the critical point, the single particle wave function has strong
spatial fluctuations and is multifractal in nature. This can be characterized
by using the inverse participation ratio (IPR):

$\displaystyle I_{q}=\int d^{d}{\bf r}|\psi({\bf r})|^{2q}\sim L^{-\tau_{q}},$
(1)

where $\psi({\bf r})$ is the normalized single particle wave function in real
space and the exponent $\tau_{q}$ is an infinite set of critical exponents
describing the moments of $|\psi({\bf r})|^{2}$ Wegner (1980); Evers and
Mirlin (2008). We can further introduce the fractional dimension $D_{q}$ via
the relation $\tau_{q}=D_{q}(q-1)$, in order to quantify how extended the wave
function is. At the critical point, $D_{q}$ takes a fractal value and has a
non-trivial dependence on $q$. In contrast, in the metallic phase, the wave
function is uniform in the space with $D_{q}=d$ (the spatial dimension of the
system), while in the localized phase, the wave function is exponentially
localized with $D_{q}=0$.

In this paper, we first consider hybrid random Clifford circuits in which the
wave functions can be represented using the stabilizer formalism Li  _et al._
(2018b, 2019). Using this formalism, the steady state wave functions can be
transformed into a so called graph state by applying only local unitary
operations. We analyse the IPR of the eigenvectors of the adjacency matrix
associated with this graph. We find that throughout the non-thermal volume law
phase, the graph has a high connectivity and $D_{q}$ takes a fractional value
between 0 and 1. Furthermore, this fractional dimension, $D_{q}$, has a non-
trivial dependence on $q$, demonstrating the multifractal nature of these
graph states. When the measurement rate $p>p_{c}$, the random graph obtains a
local structure with low connectivity and has $D_{q}=0$, analogous to the
Anderson localized phase. We further apply the same method to the hybrid
Clifford quantum automaton circuit Iaconis  _et al._ (2020) and find similar
multifractal behavior in the volume law phase.

In addition, we consider non-unitary random free fermion dynamics in one
dimension. Previous studies indicate the existence of a critical phase in this
model, which enjoys emergent two dimensional conformal symmetry with a
spacelike time direction Chen  _et al._ (2020). We observe that in the
critical steady state, the single particle wave function has strong
fluctuations in space and is multifractal in nature. This provides strong
evidence that the wave function is qualitatively different from the critical
state of the clean free fermion system, in which the single particle wave
function’s amplitude is uniform in space and the criticality comes from
quantum coherence effects.

The rest of the paper is organized as follows. In Sec. 2, we study the
multifractal behavior in the hybrid random Clifford circuit. We first review
both the stabilizer and the graph state formalism in Sec. 2.1. We then compute
the IPR of the corresponding graph state in the steady state of the hybrid
random Clifford circuit model in Sec. 2.2. We apply the same method to the
hybrid Clifford quantum automaton circuit in Sec. 2.3. In Sec. 3, we analyze
the multifractal behavior of the non-unitary free fermion dynamics. In Sec. 4
, we summarize our results and discuss possible directions for future work.

##  2 hybrid random Clifford circuit

In this section, we study the dynamics of hybrid quantum circuits which are
composed solely of 2-site random unitary gates drawn uniformly from the
Clifford group interspersed with a layer of projective Pauli measurements (See
Fig. 1). This model exhibits an entanglement phase transition from a volume
law phase to an area law phase as we vary the measurement rate $p$ . Using the
stabilizer formalism, it is possible to efficiently simulate this Clifford
dynamics and analyze the scaling of the entanglement in both of these phases
as well as at the critical point $p_{c}$. In this work, we study the so called
graph states, which can always be obtained from the steady state stabilizer
wave function Raussendorf  _et al._ (2003); Van den Nest  _et al._ (2004). We
will show that this ensemble of random graphs, which can be specified by a
corresponding adjacency matrix, possess eigenstates which exhibit a
localization transition at $p_{c}$. Furthermore, we will show that in the
regime $0<p<p_{c}$, these eigenstates show the characteristic properties of
multifractality.

###  2.1 Stabilizer formalism

In Clifford circuits, the dynamics are conveniently described using the notion
of stabilizer operators. These stabilizers can be used to completely define a
class of quantum states and simulate their quantum dynamics. In this
subsection, we give a brief review of the stabilizer formalism Gottesman
(1998); Aaronson and Gottesman (2004).

An $L$-qubit stabilizer state, $|\psi\rangle$, is completely defined as the
simultaneous eigenstate of $L$ commuting and independent Pauli string
operators $M_{i}$ with eigenvalue +1. The $M_{i}$ operators are the generators
of the stabilizer group $\mathcal{S}$, which is a subgroup of the $L$-qubit
Pauli group. We can write each generator as
$M_{i}=X_{1}^{a^{i}_{1}}Z_{1}^{b^{i}_{1}}X_{2}^{a^{i}_{2}}Z_{2}^{b^{i}_{2}}\cdots
X_{L}^{a^{i}_{L}}Z_{L}^{b^{i}_{L}}$, with $a^{i}_{n},b^{i}_{n}$ taking values
0 or 1. The information contained in these $L$ vectors,
$(a^{i}_{1},a^{i}_{2},\cdots,a^{i}_{L},b^{i}_{1},b^{i}_{2},\cdots,b^{i}_{L})$,
can be arranged into the so called “stabilizer tableau”. This is a $L\times
2L$ binary matrix $T=[T_{X},T_{Z}]$, where the first square matrix, $T_{X}$,
stores the information of $\\{a\\}$ and the second matrix, $T_{Z}$, stores the
information contained in $\\{b\\}$. For instance, a trivial product state in
the z-direction $|\psi\rangle=|00\cdots 0\rangle$ has $M_{i}=Z_{i}$ for all
$i=1\dots L$. In the corresponding stabilizer tableau, $T_{X}$ is a zero
matrix and $T_{Z}$ is the identity matrix. Under the hybrid Clifford dynamics
described in Fig. 1), the stabilizer operators, $M_{i}$, are mapped to a new
set of stabilizer operators $M_{i}^{\prime}$. Therefore, the wave function
remains a stabilizer state, with the stabilizer tableau being updated
accordingly 111We disregard the phase information since it is not important in
our paper..

The Rényi entanglement entropy for a subsystem A is defined as
$S_{n}=\frac{1}{1-n}\log_{2}\mbox{Tr}\rho_{A}^{n}$ where $\rho_{A}$ is the
reduced density matrix for region A. For the stabilizer states, $S_{n}$ is
independent of the Rényi index $n$ and obeys the form Hamma  _et al._ (2005)

$\displaystyle S_{n}=L_{A}-\log_{2}|\mathcal{S}_{A}|,$ (2)

where $\mathcal{S}_{A}$ is a stabilizer group defined in subsystem A, which is
a subgroup of $\mathcal{S}$. $|\mathcal{S}_{A}|$ counts the number of
independent stabilizers supported only on A.

In the hybrid random Clifford dynamics, when the measurement rate $p$ is
small, the steady state entanglement entropy has volume law scaling. Roughly
speaking, this implies that the corresponding stabilizer generators of this
state span the entire system. A single projective measurement replaces one
stabilizer with a local operator. As we increase $p$, the long stabilizers are
gradually replaced by a series of short stabilizers. Eventually, at a high
enough measurement rate $p>p_{c}$, the steady state entanglement entropy has
area law scaling.

Figure 1: The cartoon for the one dimensional hybrid random Clifford circuit.
The random two-qubit unitary gates are arranged in a brick-wall fashion, while
the single-qubit Z measurements are randomly applied with probability $p$ at
each time. Figure 2: The emergent graph state which results from the random
Clifford dynamics with a measurement rate corresponding to the (left) volume
law, (middle) critical and (right) area law entangled phases. Each node of the
graph represents a single qubit, with the edges corresponding to the
associated CZ gates which are applied to a reference state to create the graph
state. The placement of the nodes are determined according to the Fruchterman-
Reingold force-directed drawing algorithm which treats each node as a particle
and introduces forces between them Fruchterman and Reingold (1991). There are
attractive forces between adjacent nodes and repulsive forces between all
pairs of nodes. The position of the nodes is fixed by minimizing the energy of
the entire system. In the non-thermal volume law phase, a finite fraction of
the nodes are connected with $O(L)$ other nodes and forms the “core” of the
graph. The outer layers of the graph is composed of the nodes which are
connected with only a few number of other nodes. As we increase the
measurement rate $p$, the core of the graph becomes smaller and more nodes are
pushed to the outer layers of the graph. At the critical point, the core
disappears and the graph has a self-similar structure. In the area law phase,
most of the nodes are living on the outermost circle and are only connecting
with nearby nodes.

####  2.1.1 Graph states

In this subsection, we introduce an important subclass of stabilizer states
known as graph states. We briefly review the definition of the graph states
and their connection with the more general stabilizer states.

An undirected graph, $G$, is defined by a set of vertices $V$ and edges $E$.
For any graph, $G$, we can define a corresponding graph state
$|\psi(G)\rangle$ as

$\displaystyle|\psi\rangle=\prod_{(ij)\in
E}\Lambda_{ij}(Z)|+\rangle^{\otimes|V|}$ (3)

where $\Lambda(Z)_{ij}={\rm I}-2|11\rangle\langle 11|$ is the two-qubit
control Z (CZ) unitary operator and $|+\rangle^{\otimes|V|}$ is the reference
state with $|+\rangle=(|0\rangle+|1\rangle)/\sqrt{2}$ polarized in x
direction. The stabilizer generators of such a graph state are given by the
set

$\displaystyle M_{i}=X_{i}\prod_{j|(ij)\in E}Z_{j}.$ (4)

In the language of the stabilizer tableau, $T_{X}$ is an identity matrix.
Since all $M_{i}$ operators commute with each other, $T_{Z}$ is required to be
a symmetric binary matrix. We further require that all the diagonal elements
of $T_{Z}$ are zeros. This Z stabilizer tableau is exactly the adjacency
matrix of the graph $G$. When $T^{ij}_{Z}=1$, we have $(i,j)\in E$ and
$T^{ij}_{Z}=0$ otherwise.

A key feature for our analysis is that all stabilizer states are equivalent to
a graph state up to the application of single qubit unitary rotations S (phase
gate) or H (Hadamard gate) Van den Nest  _et al._ (2004). In the language of
the stabilizer tableau, this mapping onto a graph state can be done in two
steps: (1) We swap/add the rows of $T$ (Gaussian elimination over finite
$\mathbb{Z}_{2}$ field) to transform $T_{X}$ into an upper triangular matrix.
In this process, we may also apply H gates to enforce that all the diagonal
elements of $T_{X}$ are equal to one. (2) We add rows in $T$ to transform
$T_{X}$ into an identity matrix. We further apply S gates to enforce the
condition that $T_{Z}$ contains only zeros along its diagonal. Note that row
operations which are applied to the stabilizer tableaux do not change the
stabilizer wave function. Since the only nontrivial operations we apply are
single qubit unitary gates, the entanglement entropy is invariant under this
transformation.

As a consequence, for a stabilizer state evolved under the hybrid Clifford
dynamics, at any time, it can be transformed into a graph state with the
complete quantum information contained in the corresponding adjacency matrix.
Its entanglement entropy is closely related to the connectivity properties of
the underlying graph. In particular, if we bipartition a graph into two
subsets A and $\overline{\rm A}$ with

$\displaystyle T_{Z}=\begin{pmatrix}T_{Z}^{AA}&T_{Z}^{A\bar{A}}\\\
T_{Z}^{\bar{A}A}&T_{Z}^{\bar{A}\bar{A}}\end{pmatrix},$ (5)

the connectivity between A and $\overline{\rm A}$ can be quantified by ${\rm
rank}_{2}(T_{Z}^{A\bar{A}})$, which is exactly the entanglement entropy of A
Hein  _et al._ (2006). Stabilizer wave functions in the volume law phase are
characterized by graphs with very high connectivity which lack locality
between the connected vertices of the graph. On the other hand, in the area
law phase, the adjacency matrix is more sparse and vertices are only connected
to other vertices which are spatially nearby in the original circuit
construction.

In the following section, we will show that there exists a structural change
of the adjacency matrix across $p_{c}$ which is similar to the Anderson
localization transition of random matrix models. One of the most prominent
examples of such a model is the power-law random banded matrix ensemble, in
which the off-diagonal elements of a random matrix have zero mean and variance
$1/r^{2\alpha}$, where $r$ is the distance from the diagonal element Evers and
Mirlin (2000, 2008). This matrix describes a one dimensional random free
fermion Hamiltonian with long-range hopping. As we increase $\alpha$, the
hopping becomes increasingly local and the single particle eigenstate
undergoes an Anderson localization transition from an extended state to a
localized state. Precisely at the transition point $\alpha=1$, the critical
eigenstate is neither localized or extended. The calculation of moments
indicates that the eigenstate has strong fluctuations and is multifractal –
characterized by a infinite set of fractal dimensions Evers and Mirlin (2008).
Motivated by these studies, we map the hybrid circuit measurement driven
transition to an Anderson localization transition in the effective graph
space. We diagonalize the adjacency matrix and analyze the possible
multifractal behavior in its eigenstates.

In Fig. 2, we show three examples of characteristic graphs which describe the
random Clifford circuit stabilizer states $i)$ deep within the volume phase
$0<p<p_{c}$, $ii)$ at the critical point $p_{c}$, and $iii)$ in the area law
phase $p>p_{c}$. In the volume law phase, when there is a nonzero measurement
rate, the graph consists of a core of very highly connected nodes plus an
outer shell of nodes with low connectivity. The highly connected core leads to
volume law entanglement scaling and the same structure has also been observed
in the pure unitary evolution with $p=0$. The presence of the outer shell
nodes is due to the measurement and as we will show later, gives rise to a
multifractal structure. At exactly the critical point $p=p_{c}$, the inner
core disappears. The graph at this point has a complex self similar structure
and the entanglement entropy scales logarithmically in the subsystem size.
Finally in the area law phase, the connectivity is dramatically reduced and
the graph has an emergent local structure. In the following subsection, in
order to quantitatively characterize the structure of these graphs, we will
analyze the properties of the eigenvectors and eigenvalues of the adjacency
matrix.

###  2.2 Adjacency matrix and multifractality

As we discussed before, each graph state is specified by an adjacency matrix
$T_{Z}$. Since it is a hermitian binary matrix, we can treat it as a
Hamiltonian for a free fermion system with

$\displaystyle H=\sum_{ij}T_{Z}^{ij}c_{i}^{\dagger}c_{j}.$ (6)

This model describes a free fermion hopping on an ensemble of random graphs,
$G$, generated by the corresponding graph state of the steady state of the
hybrid Clifford circuit. We expect that in this free fermion system, there
exists an Anderson localization transition with respect to the node degrees.
When $p>p_{c}$, in the area law phase, the random graph has a local structure
and fermions can only hop to the sites which are spatially nearby. Since we
are considering the one dimensional random system, the single particle wave
function is always spatially localized. In contrast, in the volume law phase
$p<p_{c}$, the Hamiltonian has long range hopping terms and the wave function
can become delocalized.

To quantitatively characterize how extensive the wave functions are, we
consider the IPR defined as

$\displaystyle I_{q}(L)=\sum_{i=1}^{L}|\Psi(x_{i})|^{2q}.$ (7)

where $\Psi(x_{i})$ is a random eigenstate of $T_{Z}$. This quantity computes
the $q$-th moment of the eigenstate coefficients and scales as $I_{q}\sim
L^{-\tau_{q}}$ with $\tau_{q}=D_{q}(q-1)$. As we mentioned in the
introduction, $D_{q}$ is the fractional dimension and has been used to
distinguish between the extensive and localized states in the Anderson
localization transition. In our model, when $p>p_{c}$, the system is localized
and we have $D_{q}=0$ for all $q>0$. We are going to study $D_{q}$ for the
volume law phase with $p<p_{c}$.

To numerically compute $\tau(q)$ and $D_{q}$ in the random system at finite
system size, we need to take an ensemble average over $\Psi(x_{i})$. The
correct way to do this is to consider the quenched average over $\log
I_{q}(L)$, i.e.,

$\displaystyle\tau_{q}=-\frac{\langle\log I_{q}(L)\rangle}{\log L}.$ (8)

The average $\langle\cdot\rangle$ is taken over both different eigenvectors,
$\Psi_{k}$, in one realization of $T_{Z}$ and different instances of the
random Clifford circuit. We also compute the annealed average defined as

$\displaystyle\tau_{q}^{*}=-\frac{\log\langle I_{q}(L)\rangle}{\log L}.$ (9)

In random systems, this quantity is much easier to obtain analytically.
Previous experience in random matrices and spin glass systems tell us that
when the system is ergodic, the quenched average and annealed average give the
same result Evers and Mirlin (2008); Derrida (1981); Chamon  _et al._ (1996).
On the other hand, in many non-ergodic systems, $\tau_{q}$ and $\tau^{*}_{q}$
can be quite different.

By studying $I_{q}(L)$, we find that at any finite measurement rate
$0<p<p_{c}$, the behavior of the graph states generated in the hybrid random
Clifford circuits is dramatically different from the $p=0$ limit. To see this
we numerically measure both $\tau_{q}(L)$ and $\tau^{*}_{q}(L)$ at general $q$
using the ensemble of graphs generated by the random Clifford circuits with up
to $L=4000$ sites. We furthermore use these scaling exponents to extract the
value of the fractal dimension $D_{q}$ in the thermodynamic limit
$L\rightarrow\infty$. In both cases, we will see that the usual volume law
Clifford wave function in the absence of measurements ($p=0$) behaves as a
fully extended wave function, while for any finite measurement rate there is
evidence of multifractal scaling behavior.

Figure 3: Comparison of the multifractal exponents $\tau_{q}$ (circles) and
$\tau^{*}_{q}$ (squares) of the IPR as a function of $q$. (a) In the special
case when $p=0$, both $\tau_{q}$ and $\tau^{*}_{q}$ are the same regardless of
the value of $q$. (b) For non-zero measurement rate $p=0.05\ll p_{c}$,
$\tau_{q}$ and $\tau^{*}_{q}$ show drastically different behavior when $q>2$.
Their difference increases slightly as we increase the system size. For other
values of $p<p_{c}$, similar behaviors are observed, providing evidence of the
non-ergodic nature of the full volume law phase.

We first consider the limit $p=0$. We measure both $\tau_{q}$ and
$\tau^{*}_{q}$, using Eq.’s 8 and 9. As shown in Fig. 3, we find that
$\tau_{q}$ and $\tau^{*}_{q}$ take the same value for arbitrary $q$ (on top of
each other in the plot). This indicates that when the measurement rate $p=0$,
the single particle wave function is ergodic and $I_{q}(L)$ is a self-
averaging quantity. We can further use $\tau_{q}$ to extract the fractal
dimension $D_{q}$. We find that the finite size scaling of the fractal
dimension very closely follows the form

$\displaystyle D_{q}(L)\sim 1-\frac{f(q)}{\log L}.$ (10)

The fractal dimension $D_{q}$ approaches 1 for all $q$ (See Fig. 5 for the
plot of $D_{2}(L)$ at finite $L$). This is again consistent with a single
particle wave function which is fully ergodic and extensive in the
thermodynamic limit. Note that the same finite size sub-leading correction of
$D_{q}$ at finite $L$ is also observed in free fermion models which use the
Gaussian orthogonal ensemble (GOE) random matrices as the Hamiltonians Bäcker
_et al._ (2019). Finally, we also examine the level spacing statistics of the
eigenvalues in $T_{Z}$ by computing the probability distribution $P(s_{i})$ of
$s_{i}=e_{i}-e_{i+1}$, the spacing between adjacent eigenvalues. It is known
that in GOE random matrices, $P(s)$ takes the following form

$\displaystyle P(s)=\frac{\pi}{2}se^{-\frac{\pi}{4}s^{2}}.$ (11)

In Fig. 7, we plot the level spacing distribution for the ensemble of
adjacency matrices with $p=0$. To numerically obtain $P(s)$, we perform the
unfolding procedure Brody  _et al._ (1981); Avishai  _et al._ (2002), which
compensates for the non-constant density of states in the eigenvalue
distribution. We find that $P(s)$ very closely follows the same GOE form.
Overall, we find that graphs states generated from random Clifford circuits
without measurements are very well behaved. The associated adjacency matrix
shares similar properties with a GOE random matrix and its eigenstates show
behavior consistent with fully extended ergodic wave functions.

Figure 4: The derivative of the fractal dimension $D_{2}$ as a function of
measurement rate $p$, for a fixed system size $L=2000$. We see that there is a
peak in the derivative $\partial D_{2}(p)/\partial p$, at exactly the
transition point $p_{c}\approx 0.16$.

Figure 5: (Top) Finite size scaling of the fractal dimension
$D_{q}=\tau_{q}/(q-1)$. The data appears to follow the form
$D_{q}(L)\sim\frac{c}{\log(L+L_{0})}+D_{q}(\infty)$. (Bottom) The fractal
dimension extrapolated to the thermodynamic limit as a function of $p$ for
different values of $q$. Note that the non-linearity as a function of $q$ is
present throughout the volume law phase, providing evidence of a full multi-
fractal phase.

We now consider the case of a finite measurement rate. The behavior of the
wave functions for any $p>0$ is significantly different from the $p=0$ case.
The quenched exponents $\tau_{q}$ behave very differently than the annealed
exponents $\tau^{*}_{q}$ and the fractional dimension $D_{q}$ shows
complicated multifractal behaviour. First consider the behavior of $\tau_{q}$
vs that of $\tau^{*}_{q}$, for some nonzero value of $p$ deep in the volume
law phase. As shown in Fig. 3 (b), for $p=0.05$, $\tau_{q}$ and $\tau^{*}_{q}$
are equal only when the moment $q\lesssim q_{c}=2$. For $q\gtrsim 2$,
$\tau_{q}$ grows linearly with $q$ and there is a large gap between $\tau_{q}$
and $\tau^{*}_{q}$. Furthermore the gap $|\tau_{q}-\tau^{*}_{q}|$ grows with
system size, indicating that this effect will persist in the thermodynamic
limit. This discrepancy between the quenched and annealed average exists for
all $0<p<p_{c}$ and indicates that the steady state for the non-thermal volume
law phase is non-ergodic in the graph space.

We now focus on the behavior of the fractal dimension $D_{q}$ as a function of
both $q$ and measurement rate $p$, within this non-thermal phase. For any
finite size system, we measure $\tau_{q}$ using Eq. 8, and determine the
fractional dimension using the relationship $D_{q}=\tau_{q}/(q-1)$. Note that
for $q=2$, $D_{q}=\tau_{q}$. We find that the fractal dimension at fixed $q$
decreases monotonically as we increase $p$. In particular, as we cross the
phase transition point, we observe a sharp peak exactly at $p_{c}\approx 0.16$
in the derivative of $D_{q}(p)$ (See Fig. 4) for $\partial D_{2}(p)/\partial
p$ with fixed system size.

In order to calculate the precise value of the fractal dimension $D_{q}$ in
the thermodynamic limit we must perform finite size scaling. The finite size
results for $q=2$ are shown in Fig. 5. We again find that there are
logarithmic corrections in the finite size limit, and so a careful
extrapolation to the thermodynamic limit must be performed. In fact, we find
that the finite size effects are more significant for nonzero measurement
rates. We empirically find a very good fit to the form

$\displaystyle D_{q}(L)=D_{q}(\infty)+\frac{f(q)}{\log(L+L_{0})},$ (12)

where we include the additional fitting parameter $L_{0}$. The data follows
this scaling form for all systems sizes we measured from $L=256$ to $L\sim
4000$ sites. Note that $L_{0}$ is much smaller than the largest system size
$L=4000$ in the fitting and as $L\rightarrow\infty$, we recover the same form
as for the $p=0$ case. After extrapolating to the $L=\infty$ limit, we notice
that in the range $0<p<p_{c}$, $D_{2}(p)$ has a non-integer fractal value.

Using this extrapolation method, we plot $D_{q}$ as a function of p in the
thermodynamic limit. We show the results in Fig. 5, for $q=1,2$ and $8$. Note
that in the $q=1$ case, $D_{1}(q)$ is defined by the limit of the equation
$D_{q}=\tau_{q}/(q-1)$, where $\tau_{1}(p)=0$ due to the normalization of the
wave fucntion. We find that for all $q$, $D_{q}(p)$ forms a continuous curve
which interpolates between $D_{q}=1$ at $p=0$ and $D_{q}=0$ at $p=p_{c}$.
Importantly, at all finite measurement rates in the volume law phase, the
fractal dimension exhibits a strong $q$ dependence. That is, for any nonzero
measurement rate in the volume law phase, the resultant graph states show
multifractal behavior which is not present in the usual volume law phase
without measurement. This multifractal behavior is reminiscent of the critical
behavior of wave functions near an Anderson localization transition.

Figure 6: (Top) The anomalous dimension $\Delta_{q}=q-1-\tau_{q}$ (solid line)
and $\Delta_{1-q}$ (dashed line) for small $q<1$, at different measurement
rates $p$. (Bottom) We plot $\Delta_{q}/q$ vs $q$. We find that for small $p$,
this quantity is linear in $q$, ($\Delta_{q}/q\sim 1-q$) consistent with
predictions for a weak multifractal system.

We also look at the non-linearity of the multifractal exponent $\tau_{q}$ for
$0<q<1.0$. Previously, the field theory calculation at the critical point of
the Anderson localization transition suggests that the anomalous dimension
$\Delta_{q}\equiv(q-1)-\tau_{q}$, which is defined as the deviation of
$\tau_{q}$ from the fully ergodic case, is symmetric around $q=1/2$ Mirlin
_et al._ (2006); Evers and Mirlin (2008); Gruzberg  _et al._ (2011). In Fig. 6
(a), we plot $\Delta_{q}$ and $\Delta_{1-q}$ as a function of $q$ in the
thermodynamic limit and we find that they are close to each other. The
difference between them is small for small $p$ and slightly increases as we
increase $p$. When $p$ is close to zero, we find that $\Delta_{q}/q$ is a
linear function of $q$, indicating that $\Delta_{q}\sim q(1-q)$. This
parabolic form of $\Delta_{q}$ has also been observed in the critical wave
function with weak disorder Chamon  _et al._ (1996); Evers and Mirlin (2008);
Mirlin  _et al._ (2006).

Finally, we once again consider the eigenvalue spacing statistics of the
adjacency matrix, $T_{Z}$, for nonzero measurement rate. For free-fermion
models which undergo an Anderson localization transition, there is a
qualitative change of the level spacing statistics as one moves across the
critical point. In our model, as we mentioned previously, when $p=0$, the
nearest neighbor level spacing distribution $P(s)$ is described by GOE. On the
other hand, as shown in Fig. 7, near the critical point at $p=p_{c}$, the
distribution $P(s)$ exhibits significant difference from GOE and is close to a
semi-Poisson distribution Geraedts  _et al._ (2016); Bogomolny  _et al._
(1999),

$\displaystyle P(s)=4se^{-2s}.$ (13)

Namely, there exists level repulsion between adjacent energy levels as
$P(s\rightarrow 0)=0$, and the tail of the distribution of $P(s)$ appears to
decay exponentially as $\sim e^{-s}$. There appear to be small deviations from
the exact semi-Poisson distribution in the intermediate regime. We also
compute $P(s)$ in the volume law phase with $0<p<p_{c}$. Deep inside the
volume law phase, $P(s)$ appears to be described by the GOE distribution. For
measurement rates closer to the critical point, the tail of the spacing
distribution decays with some form between the GOE and semi-Poisson
distribution. This deviation from GOE might be due to finite size effects and
we expect that in the volume law phase, $P(s)$ becomes GOE in the
thermodynamic limit.

Figure 7: The eigenvalue level-spacing probability distribution of the
adjacency matrices which result from the hybrid random Clifford dynamics. The
unfolding procedure is applied only to eigenvalues in the range
$\lambda_{i}\in[-20,-3]$, where the density of states is relatively smooth.
$P(s)$ appears close to the GOE distribution deep in the volume law phase, but
appears to approach the semi-Poisson distribution near the critical point. The
top and bottom plots are the same results but on a linear scale and log scale
respectively.

###  2.3 Hybrid random Clifford quantum automaton circuit

We now consider a hybrid circuit in which the unitary dynamics is composed
solely of gates which preserve the computational basis. These are known as
quantum automaton (QA) circuits, and have been studied in Ref’s Iaconis  _et
al._ , 2020; Iaconis, 2020. The general form of a QA gate is given by

$\displaystyle U_{QA}|m\rangle=e^{i\theta_{m}}|\pi(m)\rangle,$ (14)

where $\pi(m)$ is the permutation group on the $2^{N}$ basis states
$|m\rangle$. When acting on an initial product state which has all spins
perpendicular to the computational basis, $U_{QA}$ can create complex highly
entangled wave functions Iaconis  _et al._ (2020); Iaconis (2020). In
particular, we have

$\displaystyle|+\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2^{N}}\sum_{n}|n\rangle$ (15) $\displaystyle
U_{QA}|+\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2^{N}}\sum_{n}e^{i\theta_{n}}|n\rangle.$ (16)

When a finite rate of non-unitary composite measurements (explained below) are
added to $U_{QA}$, there is again a phase transition between a volume law and
and area law phase. In this case, the universality class of the critical point
is distinct from that of the hybrid random Clifford circuit discussed in Sec.
2.2 which possesses an emergent conformal symmetry and therefore has critical
exponent $z=1$ Li  _et al._ (2020). It was shown in Ref. Iaconis  _et al._ ,
2020 that the critical point of the hybrid QA circuit falls exactly in the
directed percolation universality class. This is a well known non-equilibrium
critical point which has a dynamical critical exponent $z=1.581$.

Figure 8: The cartoon for the Clifford QA circuit subject to Z measurement
followed by Hadamard gate. The dashed box denotes the time evolution in one
time step. The unitary QA circuit is composed of CNOT and CZ gates. There are
two types of CNOT gates (control qubit on the left or the right qubit) and we
apply them randomly with equal probability.

In what follows, we consider the case where the hybrid QA circuit contains
unitary gates which belong to the Clifford group. We include a finite rate of
composite measurements, whereby a spin is projectively measured in the Pauli Z
basis followed by a Hadamard rotation (see Fig. 8). A Clifford QA circuit of
this form was also studied in Ref. Iaconis  _et al._ , 2020, and found to
display the same characteristics as the generic QA hybrid circuit. In this
model, there exists an entanglement phase transition between a volume law and
area law phase at $p_{c}=0.075$. We now study the multifractal properties of
the adjacency matrices formed by the steady state graph states of the hybrid
Clifford QA circuits. Notice that in this model, the wave function is a graph
state up to single qubit unitary S rotations (H rotation is not required).

Figure 9: The IPR scaling exponent $\tau_{q}$ and the annealed scaling
exponent $\tau^{*}_{q}$ for the QA circuit at $p=0.06$ in the non-thermal
volume law phase. Once again, we find that for any finite measurement rate,
$\tau$ and $\tau^{*}$ show dramatically different behavior at large $q$.

When the measurement rate $p=0$, the resulting eigenstates of the adjacency
matrix with respect to the node degrees are fully extended and ergodic. In
particular the quenched and annealed multifractal exponents are the same,
i.e., $\tau_{q}=\tau^{*}_{q}$. Using the same scaling form as in Eq. 10, in
the thermodynamic limit, we find that $\tau_{q}\approx q-1$ with
$D_{q}=\tau_{q}/(q-1)\approx 1$ for all $q$. In Fig. 12, we plot the
eigenvalue level spacing distribution. We find that for QA circuits without
measurements, the distribution follows the GOE form given by Eq. 11 with the
wave function being fully ergodic and extended.

On the other hand, for any finite measurement rate, we again find that all
steady state wave functions in the volume law phase display multifractal
behavior. In addition, at large $q$, there is a large discrepancy between the
quenched exponent $\tau_{q}$, and the annealed exponent $\tau_{q}^{*}$. This
can be seen in Fig. 9, where at $p=0.06$, there is a large gap between
$\tau_{q}$ and $\tau^{*}_{q}$ at large $q$, and this gap increases as we move
towards the thermodynamic limit. This signals that there is a breakdown of
ergodicity also in this QA model.

We further extract the fractal dimension $D_{q}$ as a function of measurement
rate for different values of $q$. We use the scaling form in Eq. 12, to
extrapolate $D_{q}$ to the thermodynamic limit. We show these results in Fig.
10, where we see very similar behavior as with the random Clifford circuit in
Sec. 2.2. In particular, $D_{q}$, when scaled to the thermodynamic limit,
takes a non-integer fractal value for $0<p<p_{c}$. This fractal dimension goes
to zero at the critical point and remains zero in the area law phase for all
$q$. Importantly, in the volume law phase, the value of $D_{q}$ has a strong
$q$ dependence. This indicates that, again, the wave functions show
complicated multifractal behavior throughout the entire volume law phase. Note
that the finite size effects are somewhat larger in the QA model, and
therefore the extracted curves for $D_{q}(L=\infty)$ are not as smooth as in
the random Clifford model in Sec. 2.2.

In Fig. 11, we plot both the anomalous dimension $\Delta_{q}$, and
$\Delta_{q}/q$ for $0<q<1$. For small measurement rates,
$\Delta_{q}\approx\gamma q(1-q)$, for some small constant $\gamma$. This
nonlinear functional form also follows the predictions for weak multifractal
systems.

ummWe finally also study the level spacing statistics near the critical point
$p_{c}$ in the Clifford QA model. The results are show in Fig. 12. We find
that the ensemble of adjacency matrices show level spacing statistics which
are very close the the semi-Poisson distribution. The tail of the distribution
decays like $e^{-s}$, and due to level repulsion $P(s)$ goes to zero as
$s\rightarrow 0$. In summary, we also find multifractal behavior in the steady
state wave functions of the non-thermal volume law phase of the random
Clifford QA circuit. In future works, the more general QA model might provide
a tractable platform for extending these results to non-Clifford systems.

Figure 10: The fractal dimension $D_{q}$ for different measurement rates, $p$,
in the thermodynamic limit, for the QA circuit. Note that the finite size
effects appear to be larger in the QA circuit and performing the scaling to
the thermodynamic limit is less precise than in the random Clifford circuit.

Figure 11: (Top) The anomalous dimension $\Delta_{q}$ and $\Delta_{1-q}$ for
small $q<1$ in the QA circuit. For measurement rate $p$ close to zero,
$\Delta_{q}\approx\Delta_{1-q}$. (Bottom) We also show $\Delta_{q}/q$, and
find that for small $p$, this is linear is $q$, which is again consistent with
predictions for weak multifractal behavior.

Figure 12: The unfolded level spacing statistics of the adjacency matrix for
the QA Clifford wave functions at $p=0.0$ and at the critical point $p=0.075$.
We again find that the level spacing distribution of the adjacency matrix
appears to follow a GOE distribution for $p=0$ and is close to the semi-
Poisson distribution at the critical point.

##  3 non-unitary random free fermion dynamics

We consider the one dimensional discrete non-unitary free fermion circuit
defined in Ref. Chen  _et al._ , 2020. The (unnormalized) non-unitary time
evolution operator is defined as

$\displaystyle U=\prod_{t=1}^{T}U_{\beta}(t)U_{\tau}(t),$ (17)

where $U_{\tau}(t)=\exp(-2i\tau H_{1}(t))$ denotes the unitary evolution with
$H_{1}(t)=\sum_{x}\kappa_{x,t}c^{\dagger}_{x}c_{x+1}+H.C.$ and
$U_{\beta}(t)=\exp(-2\beta H_{2}(t))$ denotes the imaginary evolution governed
by a random onsite potential
$H_{2}(t)=\sum_{x}\lambda_{x,t}c^{\dagger}_{x}c_{x}$. Both $\kappa_{x,t}$ and
$\lambda_{x,t}$ are random in space and time with the following simple
distribution:

$\displaystyle P(\kappa_{x,t})$
$\displaystyle=\frac{1}{2}\delta(\kappa_{x,t}-1)+\frac{1}{2}\delta(\kappa_{x,t}+1)$
$\displaystyle P(\lambda_{x,t})$
$\displaystyle=\frac{1}{2}\delta(\lambda_{x,t}-1)+\frac{1}{2}\delta(\lambda_{x,t}).$
(18)

Under this non-unitary random evolution, the wave function evolves as

$\displaystyle|\psi(T)\rangle=\frac{U}{\sqrt{Z}}|\psi(0)\rangle,\quad\text{with}\
\ Z=\langle\psi(0)|U^{\dagger}U|\psi(0)\rangle.$ (19)

It is shown in Ref. Chen  _et al._ , 2020 that this non-unitary dynamics has
an emergent two dimensional conformal symmetry for any arbitrary $\beta>0$.
The steady state $|\psi(T\to\infty)\rangle$ is critical and has an
entanglement entropy which scales logarithmically in the subsystem size, the
same as for ground states of critical systems.

If we start with an initial pure Gaussian state, under the non-unitary
evolution, $|\psi(T)\rangle$ remains a pure Gaussian state and can be simply
written as a product state in some suitable basis, i.e.,

$\displaystyle|\psi(T)\rangle=\prod_{n=1}^{N}f_{n}^{\dagger}(T)|0\rangle,$
(20)

where $N$ is the total number of the fermions and is conserved under the time
evolution. $\\{f_{n}^{\dagger}(T)\\}$ with $n=1,\cdots N$ form a canonical
basis for fermion creation operators at time $T$ and satisfy
$\\{f_{m}^{\dagger},f_{n}^{\dagger}\\}=0$. They can be expanded in the
$\\{c_{x}^{\dagger}\\}$ basis defined in the spatial direction, i.e.,

$\displaystyle f_{n}^{\dagger}=\sum_{x=1}^{L}u_{x}c_{x}^{\dagger},$ (21)

where $L$ is the system size and $u_{x}$ satisfies the normalization
constraint $\sum_{x}|u_{x}|^{2}=1$.

We study the spatial distribution of $f_{n}^{\dagger}$ by computing the q-th
moment of the wave function defined as

$\displaystyle I_{q}(L)=\sum_{x=1}^{L}|u_{x}|^{2q}\sim L^{-\tau_{q}}.$ (22)

As we have done for the Clifford circuit, we compute both $\tau_{q}$ and
$\tau_{q}^{*}$ for the steady state and we find that they are different when
$q>q_{c}$, indicating the nonergodicity of the single particle wave function
(See Fig. 13 (a)). Numerically, we find that $q_{c}$ depends on $\beta$ and
decreases as we increase $\beta$. We also observe that $\tau_{q}$ has a non-
trivial dependence on $q$. When $q$ is large, $\tau_{q}=\alpha q$ with
$\alpha<1$, while for small $q$, $\tau_{q}$ is non-linear in $q$. We further
analyze the finite size effect in $\tau_{q}$ and present the numerical results
for $D_{2}(L)$ in Fig. 13 (b). Again, we observe that $D_{2}(L)$ has a
logarithmic correction at finite $L$ and slowly converges to
$D_{2}(L\to\infty)$, which is between 0 and 1 and depends on the value of
$\beta$. Similar results are obtained for other $q>1$ and are presented in
Fig. 13 (c). In the unitary evolution limit $\beta=0$, we expect that
$D_{q}\to 1$ and the single particle wave function is uniform in space. On the
other hand, in the limit $\beta=\infty$, $D_{q}=0$ for $q>0$ and the wave
function is localized in the space. At finite $\beta$, $D_{q}$ interpolates
between 0 and 1 and forms a set of continuous exponents depending on the
parameter $q$, implying a multifractal structure of the wave function.

Figure 13: (a) $\tau_{q}$ and $\tau_{q}^{*}$ in the non-unitary free fermion
dynamics at various $\beta$. The solid yellow curve is a guideline with
$\tau_{q}=q-1$. (b) $D_{2}(L)$ vs $1/\log L$ at various $\beta$. The crossing
of the curves with the vertical axis gives $D_{2}(L\to\infty)$. (c)
$D_{q}(L\to\infty)$ as a function $\beta$ for various $q$. (d) The anomalous
dimension $\Delta_{q}$ and $\Delta_{1-q}$ as a function of $q$ at various
$\beta$. The solid curve is for $\Delta_{q}$ and the dashed curve is for
$\Delta_{1-q}$. (e) $\Delta_{q}/q$ vs $q$ at various $\beta$. The curve with
$\beta=1$ is close to a straight line.

In addition, we also compute $\tau_{q}$ in the limit $L\to\infty$ in the
regime $0<q<1$ by using the same extrapolation method. We numerically extract
$\tau_{q}(L\to\infty)$ and present the result for $\Delta_{q}$ in Fig. 13 (d).
It is approximately symmetric around $1/2$. In particular, we observe that
$\Delta_{q}\approx\gamma q(1-q)$ when $\gamma$ is small (See the $\beta=1$
curve in Fig. 13 (e)), similar to what we have observed for the Clifford
circuit.

The non-trivial dependence of $\tau_{q}$ on $q$ implies that this critical
wave function has strong spatial variation and is distinct from the critical
state in the free fermion model without disorder, in which the single particle
wave function is extended and has a uniform distribution in space. The
criticality in these clean systems is caused by the delicate quantum coherence
which is fragile and can be easily destroyed when randomness is introduced.

##  4 Discussion and Conclusion

We investigated the multifractal behavior in two classes of non-unitary random
dynamics by numerically examining the inverse participation ratio. In the
hybrid random Clifford circuit, we transform the steady state wave function to
a graph state characterized by an adjacency matrix. We compute the eigenstates
of the adjacency matrix and observe multifractal behavior in the graph space
in the non-thermal volume law phase. We further obtain similar multifractal
behavior in the volume law phase of the hybrid Clifford QA circuit. We expect
that the multifractality can also be observed in the volume law phase of other
random non-unitary Clifford circuit with discrete symmetry Lavasani  _et al._
(2021); Sang and Hsieh (2021) or generated with only measurement gates
Ippoliti  _et al._ (2021).

On the other hand, in the non-unitary random free fermion dynamics, the steady
state is critical and can be written as $\prod_{n}f^{\dagger}_{1}\cdots
f^{\dagger}_{N}|0\rangle$ in some suitable basis. We numerically confirm that
these single particle wave functions have strong fluctuations and are
multifractal in real space. Previously, the multifractal exponents have been
analytically computed in various disordered free fermion systems, including
the two dimensional Dirac fermion in a random potentialChamon  _et al._
(1996); Castillo  _et al._ (1997); Gruzberg  _et al._ (2011) and the power-law
random banded matrix Evers and Mirlin (2008). It would be interesting to
generalize and apply these techniques to the non-unitary dynamics in order to
analytically compute multifractal exponents, perhaps in some large N non-
unitary modelsZhang  _et al._ (2021).

In the non-unitary free fermion dynamics, the random imaginary potential can
be replaced by a continuous weak measurementAlberton  _et al._ (2020); Cao
_et al._ (2019). In such hybrid dynamics, when the measurement strength is
small, the weak measurement is analogous to a random imaginary potential and
leads to a similar critical phase with multifractality. As we further increase
the measurement rate, a phase transition to an area law entangled phase
occurs. The fermions are now fully localized due to the measurement and there
is no longer multifractality.

In a broad sense, the phase transitions in both the hybrid Clifford circuit
and the hybrid free fermion circuit with weak measurement are “Anderson
localization”-like. In the former, there is a localization transition in the
associated graph space while in the latter, the localization occurs in real
space. In both models, before it enters into the Anderson localized phase,
there exists a phase in which the wave function has strong fluctuations and is
multifractal in nature. This is different from conventional Anderson
localization, where multifractality appears only at the critical point.
Recently, some disordered free fermion models exhibiting Anderson localization
transitions have been constructed, in which a non-ergodic metallic phase with
multifractality is identified, similar to what we have found in this paper De
Luca  _et al._ (2014); Kravtsov  _et al._ (2015). In particular, the
multifractal behavior observed in both the non-unitary Clifford and free
fermion models belong to weak mulitfractality class, characterized by a linear
growth of $\tau_{q}$ with large $q$ and an approximate parabolic form in small
$q$ Evers and Mirlin (2008); Derrida (1981). It would be interesting to
discover non-unitary random dynamics with strong multifractal behavior in
which $\tau_{q}$ becomes zero above $q>q_{c}$.

It is widely believed that there is a generic measurement driven transition in
an interacting system which has neither a stabilizer representation nor can be
described by a simple free fermion dynamics. We expect that this non-thermal
volume law phase has strong random fluctuations and is still multifractal in
nature. However, we are unaware of any good approach to characterize the
multifractality correctly and it is also unclear if we can map the measurement
induced transition to an Anderson localization transition in a proper basis.
We leave these interesting problems for future study.

Note added: During the completion of this work, we became aware of a work
investigating multifractality in the hyrbid Haar random circuit at the
critical point in a different context Zabalo  _et al._ (2021).

###### Acknowledgements.

This work is partially supported (J.I.) by the National Science Foundation
under Grant Number 1734006. J.I. is supported by a Simons Investigator Award
to Leo Radzihovsky from the Simons Foundation. This work is performed in part
(X.C.) at Aspen Center for Physics, which is supported by National Science
Foundation grant PHY-1607611.

## References

  * Li  _et al._ (2018a) Y. Li, X. Chen, and M. P. Fisher, Physical Review B 98, 205136 (2018a).
  * Skinner  _et al._ (2019) B. Skinner, J. Ruhman, and A. Nahum, Physical Review X 9, 031009 (2019).
  * Gullans and Huse (2020a) M. J. Gullans and D. A. Huse, Physical Review X 10, 041020 (2020a).
  * Chan  _et al._ (2019) A. Chan, R. M. Nandkishore, M. Pretko, and G. Smith, Physical Review B 99, 224307 (2019).
  * Zabalo  _et al._ (2020) A. Zabalo, M. J. Gullans, J. H. Wilson, S. Gopalakrishnan, D. A. Huse, and J. Pixley, Physical Review B 101, 060301 (2020).
  * Gullans and Huse (2020b) M. J. Gullans and D. A. Huse, Physical review letters 125, 070606 (2020b).
  * Li  _et al._ (2019) Y. Li, X. Chen, and M. P. A. Fisher, Phys. Rev. B 100, 134306 (2019).
  * Iaconis  _et al._ (2020) J. Iaconis, A. Lucas, and X. Chen, Physical Review B 102, 224311 (2020).
  * Jian  _et al._ (2020) C.-M. Jian, Y.-Z. You, R. Vasseur, and A. W. Ludwig, Physical Review B 101, 104302 (2020).
  * Bao  _et al._ (2020) Y. Bao, S. Choi, and E. Altman, Physical Review B 101, 104301 (2020).
  * Tang and Zhu (2020) Q. Tang and W. Zhu, Phys. Rev. Research 2, 013022 (2020).
  * Fan  _et al._ (2020) R. Fan, S. Vijay, A. Vishwanath, and Y.-Z. You, arXiv preprint arXiv:2002.12385 (2020).
  * Li and Fisher (2020) Y. Li and M. P. A. Fisher, “Statistical mechanics of quantum error-correcting codes,” (2020), arXiv:2007.03822 [quant-ph] .
  * Choi  _et al._ (2020) S. Choi, Y. Bao, X.-L. Qi, and E. Altman, Physical Review Letters 125, 030505 (2020).
  * Li  _et al._ (2020) Y. Li, X. Chen, A. W. W. Ludwig, and M. P. A. Fisher, “Conformal invariance and quantum non-locality in hybrid quantum circuits,” (2020), arXiv:2003.12721 [quant-ph] .
  * Li  _et al._ (2021) Y. Li, S. Vijay, and M. P. A. Fisher, “Entanglement domain walls in monitored quantum circuits and the directed polymer in a random environment,” (2021), arXiv:2105.13352 [cond-mat.stat-mech] .
  * Wegner (1980) F. Wegner, Zeitschrift für Physik B Condensed Matter 36, 209 (1980).
  * Evers and Mirlin (2008) F. Evers and A. D. Mirlin, Rev. Mod. Phys. 80, 1355 (2008).
  * Derrida (1981) B. Derrida, Phys. Rev. B 24, 2613 (1981).
  * Li  _et al._ (2018b) Y. Li, X. Chen, and M. P. A. Fisher, Phys. Rev. B 98, 205136 (2018b).
  * Chen  _et al._ (2020) X. Chen, Y. Li, M. P. A. Fisher, and A. Lucas, Physical Review Research 2, 033017 (2020).
  * Raussendorf  _et al._ (2003) R. Raussendorf, D. E. Browne, and H. J. Briegel, Physical Review A 68 (2003), 10.1103/physreva.68.022312.
  * Van den Nest  _et al._ (2004) M. Van den Nest, J. Dehaene, and B. De Moor, Physical Review A 69 (2004), 10.1103/physreva.69.022316.
  * Gottesman (1998) D. Gottesman, “The heisenberg representation of quantum computers,” (1998), arXiv:quant-ph/9807006 [quant-ph] .
  * Aaronson and Gottesman (2004) S. Aaronson and D. Gottesman, Physical Review A 70, 052328 (2004).
  * Note (1) We disregard the phase information since it is not important in our paper.
  * Hamma  _et al._ (2005) A. Hamma, R. Ionicioiu, and P. Zanardi, Phys. Rev. A 71, 022315 (2005).
  * Fruchterman and Reingold (1991) T. M. Fruchterman and E. M. Reingold, Software: Practice and experience 21, 1129 (1991).
  * Hein  _et al._ (2006) M. Hein, W. Dür, J. Eisert, R. Raussendorf, M. Nest, and H.-J. Briegel, arXiv preprint quant-ph/0602096 (2006).
  * Evers and Mirlin (2000) F. Evers and A. D. Mirlin, Phys. Rev. Lett. 84, 3690 (2000).
  * Chamon  _et al._ (1996) C. d. C. Chamon, C. Mudry, and X.-G. Wen, Phys. Rev. Lett. 77, 4194 (1996).
  * Bäcker  _et al._ (2019) A. Bäcker, M. Haque, and I. M. Khaymovich, Physical Review E 100 (2019), 10.1103/physreve.100.032117.
  * Brody  _et al._ (1981) T. A. Brody, J. Flores, J. B. French, P. A. Mello, A. Pandey, and S. S. M. Wong, Rev. Mod. Phys. 53, 385 (1981).
  * Avishai  _et al._ (2002) Y. Avishai, J. Richert, and R. Berkovits, Phys. Rev. B 66, 052416 (2002).
  * Mirlin  _et al._ (2006) A. D. Mirlin, Y. V. Fyodorov, A. Mildenberger, and F. Evers, Phys. Rev. Lett. 97, 046803 (2006).
  * Gruzberg  _et al._ (2011) I. A. Gruzberg, A. W. W. Ludwig, A. D. Mirlin, and M. R. Zirnbauer, Phys. Rev. Lett. 107, 086403 (2011).
  * Geraedts  _et al._ (2016) S. D. Geraedts, R. Nandkishore, and N. Regnault, Phys. Rev. B 93, 174202 (2016).
  * Bogomolny  _et al._ (1999) E. B. Bogomolny, U. Gerland, and C. Schmit, Phys. Rev. E 59, R1315 (1999).
  * Iaconis (2020) J. Iaconis, “Quantum state complexity in computationally tractable quantum circuits,” (2020), arXiv:2009.05512 [quant-ph] .
  * Lavasani  _et al._ (2021) A. Lavasani, Y. Alavirad, and M. Barkeshli, Nature Physics 17, 342–347 (2021).
  * Sang and Hsieh (2021) S. Sang and T. H. Hsieh, Physical Review Research 3 (2021), 10.1103/physrevresearch.3.023200.
  * Ippoliti  _et al._ (2021) M. Ippoliti, M. J. Gullans, S. Gopalakrishnan, D. A. Huse, and V. Khemani, Physical Review X 11 (2021), 10.1103/physrevx.11.011030.
  * Castillo  _et al._ (1997) H. E. Castillo, C. de C. Chamon, E. Fradkin, P. M. Goldbart, and C. Mudry, Phys. Rev. B 56, 10668 (1997).
  * Zhang  _et al._ (2021) P. Zhang, S.-K. Jian, C. Liu, and X. Chen, arXiv preprint arXiv:2104.04088 (2021).
  * Alberton  _et al._ (2020) O. Alberton, M. Buchhold, and S. Diehl, arXiv preprint arXiv:2005.09722 (2020).
  * Cao  _et al._ (2019) X. Cao, A. Tilloy, and A. D. Luca, SciPost Phys. 7, 24 (2019).
  * De Luca  _et al._ (2014) A. De Luca, B. Altshuler, V. Kravtsov, and A. Scardicchio, Physical Review Letters 113 (2014), 10.1103/physrevlett.113.046806.
  * Kravtsov  _et al._ (2015) V. E. Kravtsov, I. M. Khaymovich, E. Cuevas, and M. Amini, New Journal of Physics 17, 122002 (2015).
  * Zabalo  _et al._ (2021) A. Zabalo, M. J. Gullans, J. H. Wilson, R. Vasseur, A. W. W. Ludwig, S. Gopalakrishnan, D. A. Huse, and J. H. Pixley, “Operator scaling dimensions and multifractality at measurement-induced transitions,” (2021), arXiv:2107.03393 [cond-mat.dis-nn] .


# Integrated Speech Enhancement Method Based on Weighted Prediction Error and
DNN for Dereverberation and Denoising

###### Abstract

Both reverberation and additive noises degrade the speech quality and
intelligibility. the weighted prediction error (WPE) performs well on
dereverberation but with limitations. First, The WPE doesn’t consider the
influence of the additive noise which degrades the performance of
dereverberation. Second, it relies on a time-consuming iterative process, and
there is no guarantee or a widely accepted criterion on its convergence. In
this paper, we integrate deep neural network (DNN) into WPE for
dereverberation and denoising. DNN is used to suppress the background noise to
meet the noise-free assumption of WPE. Meanwhile, DNN is applied to directly
predict spectral variance of the target speech to make the WPE work without
iteration. The experimental results show that the proposed method has a
significant improvement in speech quality and runs fast.

Index Terms: Speech enhancement, DNN, WPE

##  1 Introduction

In real-world environments, speech signals are often contaminated by additive
noises. In an enclosed space, such as a living room, the signals are also
corrupted by their reflections from walls and other surfaces. Both the
reverberation and additive noises degrade the audible quality of speech signal
and the performance of automatic speech recognition (ASR)[1, 2, 3]. Therefore,
a speech enhancement processing is requested to deal with the reverberation
and noises simultaneously.

Over the past decades, several single- and multi-channel dereverberation
approaches have been proposed, which can be broadly categorized into acoustic
channel equalization[4, 5, 6], spectral enhancement[7] and probabilistic
model-based approaches[8, 9]. Acoustic channel equalization techniques remove
the reverberation by reconstructing the room impulse responses (RIRs) between
the acoustic source and the microphone. Although perfect dereverberation is
achievable in theory, the performance of such methods heavily depends on the
RIRs estimation, which is a tough problem in practice[10, 11]. Most spectral
enhancement methods, assuming that the early and later reflections are
mutually independent, are derived from speech denoising methods. These methods
always introduce disruptive speech distortion. Other speech dereverberation
approaches are based on statistical acoustics model. One of the representative
methods is WPE[12, 13] which models the reverberation with an autoregressive
(AR) process and uses the Maximum likelihood (ML) estimation for
dereverberation. To obtain the satisfactory results, the WPE relies on an
iterative procedure to optimize the AR weights and desired speech spectral
variance.

Although the WPE performs well, it has several limitations. First, since it is
hard to satisfy the noise-free assumption, the performance of WPE is always
influenced by background noise in real environments. Second, the original WPE
employs an iterative procedure which is time consuming. Third, there is no
guarantee on the convergence of the prediction weights and the performance may
reduce when more iterations are applied[13]. Therefore, several studies (e.g
in[14]) focus on replacing the iterative procedure.

Recently, deep neural networks (DNNs) show their power on the speech
enhancement[15, 16, 17, 18]. The DNN-based methods usually predict the
magnitude spectrogram of interest signal or a mask to remove the undesired
parts[19, 20, 21]. Although, the DNN can build the non-linear relationships
between the mixture and the target by training on large amounts of data, it is
totally data-driven without considering the speech signal processing theory.
In this paper, the DNN and WPE are combined. The DNN is applied to remove the
background noise and estimate the desired speech spectral variance directly,
which can reduce the influence of background noise to WPE and make it work
without iteration.

The rest of this paper is organized as follows. In Section 2, the problem
formulation and the WPE algorithm are introduced. Section 3 presents the
proposed model and algorithm in detail. Experiments and evaluations are given
in Section 4. Finally, Section 5 provides the conclusions.

##  2 WPE METHOD

In this section, a brief description of the WPE method from a statistical
viewpoint is presented. A scenario where a single speech source is captured by
$M$ microphones is considered. In the STFT domain, $s_{n,k}$ denotes the clean
speech signal with time frame index $n\in\\{1,...,N\\}$ and frequency bin
index $k\in\\{1,...,K\\}$. The speech signal at the $m$-th
$\left(m\in\\{1,...,M\\}\right)$ microphone, $x^{m}_{n,k}$, can be modeled as

$x^{m}_{n,k}=\sum_{l=0}^{L_{h}-1}\left(h^{m}_{l,k}\right)^{*}s_{n-l,k}+e^{m}_{n,k},\\\
$ (1)

where $h^{m}_{l,k}$ is an approximation of the acoustic transfer function
(ATF) between the speech source and the $m$-th microphone with length of
$L_{h}$, and $\left(\cdot\right)^{*}$ denotes the complex conjugate operator.
The additive term $e^{m}_{n,k}$ represents the sum of modeling errors and
background noise. In WPE method, the $e_{n,k}^{m}(\forall n,k,m)$ is always
assumed to 0. The convolutive model in (1) is often rewritten as

$x^{m}_{n,k}=\sum_{l=0}^{D-1}\left(h^{m}_{l,k}\right)^{*}s_{n-l,k}+\sum_{l=D}^{L_{h}-1}\left(h^{m}_{l,k}\right)^{*}s_{n-l,k},\\\
$ (2)

where the $\sum_{l=0}^{D-1}\left(h^{m}_{l,k}\right)^{*}s_{n-l,k}$ indicates
the sum of the direct signal and the early reflection at the $m$-th
microphone. And $\sum_{l=D}^{L_{h}-1}\left(h^{m}_{l,k}\right)^{*}s_{n-l,k}$
denotes the later reflection. $D$ corresponds to the duration of the early
reflections. To simplify the expression, we rewritten (2) as

$x^{m}_{n,k}=d^{m}_{n,k}+\sum_{l=D}^{L_{h}-1}\left(h^{m}_{l,k}\right)^{*}s_{n-l,k}.\\\
$ (3)

Because the early reflections signal can actually improve speech
intelligibility[22], most dereverberation methods aim to reconstruct
$d^{m}_{n,k}$ as the desired speech. By replacing the convolutive model in (3)
with an AR model, the signal observed at the first microphone ($m=1$) can be
rewritten in the well-known multi-channel linear prediction (MCLP) form

$x^{1}_{n,k}=d^{1}_{n,k}+\sum_{m=0}^{M}\left(\textbf{g}^{m}_{k}\right)^{H}\textbf{x}^{m}_{n-D,k}\\\
$ (4)

where $d^{1}_{n,k}$ is the desired signal, and $(\cdot)^{H}$ denotes the
conjugate transposition operator. The vector
$\textbf{g}^{m}_{k}\in\mathbb{C}^{L_{k}}$ is the regression vector of order
$L_{k}$ for the $m$-th channel. $\textbf{x}^{m}_{n,k}$ and
$\textbf{g}^{m}_{k}$ are defined as

$\textbf{x}^{m}_{n-D,k}=\left[x^{m}_{n-D,k},x^{m}_{n-D-1,k},...,x^{m}_{n-D-(L_{k}-1),k}\right]^{T},\\\
$
$\begin{split}\textbf{g}^{m}_{k}=\left[g^{m}_{0,k},g^{m}_{1,k},...,g^{m}_{L_{k}-1,k}\right]^{T},\end{split}$
(5)

where $(\cdot)^{T}$ denotes the transposition operator. The MCLP model (4) can
be written in a compact form using the multi-channel regression vector
$\textbf{g}_{k}\in\mathbb{C}^{ML_{k}}$

$x^{1}_{n,k}=d_{n,k}+\textbf{g}^{H}_{k}\textbf{x}_{n-D,k}\\\ $ (6)

where

$d_{n,k}\equiv d^{1}_{n,k},\\\ $
$\textbf{g}_{k}=\left[(\textbf{g}^{1}_{k})^{T},...,(\textbf{g}^{M}_{k})^{T}\right]^{T},\\\
$
$\begin{split}\textbf{x}_{n,k}=\left[(\textbf{x}^{1}_{n,k})^{T},...,(\textbf{x}^{M}_{n,k})^{T}\right]^{T}.\end{split}$
(7)

From (6), the desired signal can be rewritten as

$d_{n,k}=x^{1}_{n,k}-\textbf{g}^{H}_{k}\textbf{x}_{n-D,k}.\\\ $ (8)

In WPE method, the desired signal in each frequency bin can be modeled as a
circular complex Gaussian distribution with zero-mean and frequency-dependent
variance $\sigma^{2}_{d_{n,k}}$. Assuming independence across time frames, by
using the maximum likelihood (ML) estimation of the desired speech $d_{n,k}$
at each frequency, the joint distribution of the desired speech coefficients
at frequency bin, $k$, is given by

$p\left(\textbf{d}_{k}\right)=\prod^{N}_{n=1}\frac{1}{\pi\sigma^{2}_{d_{n,k}}}\mathbf{\exp}\left(-\frac{\left|d_{n,k}\right|^{2}}{\sigma^{2}_{d_{n,k}}}\right),\\\
$ (9)

where $\sigma^{2}_{d_{n,k}}$ is the time-varying spectral variance of the
desired speech. By inserting $d_{n,k}$ from (8) into (9) and taking the
negative of logarithm of $p(\textbf{d}_{k})$ in (9), the objective function
can be written as

$\begin{split}\mathcal{J}(\Theta_{k})&=-\log p(\textbf{d}_{k}|\Theta_{k})\\\
&=\sum_{n=1}^{N}(\log\sigma^{2}_{d_{n,k}}+\frac{|x^{1}_{n,k}-\textbf{g}^{H}_{k}\textbf{x}_{n-D,k}|^{2}}{\sigma^{2}_{d_{n,k}}}),\\\
\end{split}$ (10)

where constant terms are ignored.
$\Theta_{k}=\\{\textbf{g}_{k},\sigma^{2}_{d_{1,k}},\sigma^{2}_{d_{2,k}},...,\sigma^{2}_{d_{N,k}}\\}$
are the unknown parameters for the $k$-th frequency bin. These parameters can
be split into two groups: $\\{g_{k}\\}$, the AR weights, and
$\\{\sigma^{2}_{d_{2,k}},...,\sigma^{2}_{d_{N,k}}\\}$, the spectral variance
of the desired speech. A two-step algorithm to minimize the objective function
is adopted by optimize the AR weights and desired speech spectral variance,
alternatively. First, fix the $\textbf{g}_{k}$,
$\\{\sigma^{2}_{d_{1,k}},\sigma^{2}_{d_{2,k}},...,\sigma^{2}_{d_{N,k}}\\}$ are
adopted to minimize objective function. The estimated ${\sigma}^{2}_{d_{n,k}}$
can be obtained by

${\hat{\sigma}^{2}_{d_{n,k}}}=|{d_{n,k}}|^{2},\quad n=1,2,...,N.\\\ $ (11)

It is the square of the desired speech spectrum of the first channel, see (8).
Then, fix the
$\\{\sigma^{2}_{d_{1,k}},\sigma^{2}_{d_{2,k}},...,\sigma^{2}_{d_{N,k}}\\}$,
and $\textbf{g}_{k}$ is applied to minimize the objective function. The
estimated $\textbf{g}_{k}$ can be obtained by

$\hat{\textbf{g}}_{k}=\left(\sum_{n=1}^{N}\frac{\textbf{x}_{n-D,k}\textbf{x}^{H}_{n-D,k}}{\sigma^{2}_{d_{n,k}}}\right)^{-1}\sum_{n=1}^{N}\frac{\textbf{x}_{n-D,k}(x^{1}_{n,k})^{*}}{\sigma^{2}_{d_{n,k}}}.\\\
$ (12)

The estimated $\hat{\textbf{g}}_{k}$ is then used in (8) to obtain $d_{n,k}$.
The above two steps are repeated until some convergence criterions are
satisfied or a maximum number of iteration is exceeded.

##  3 PROPOSED METHOD

As described in section 1, the WPE method has three limitations: the noiseless
assumption, high time complexity of iteration and no guarantee on the
convergence. In this paper, DNN is introduced to overcome these three
limitations. First, we generate noiseless signals to meet the WPE assumption.
Second, the iteration is removed to reduce the time complexity by estimating
the desired speech spectral variance directly with DNN. In original WPE, the
iteration is required by the alternative optimization of the AR weights and
the desired speech spectral variance. We obtain the desired speech spectral
variance directly from the DNN, and only need to get the AR weights by (12).
Therefore, the iteration is not needed any more. Without iterations, there is
no need to consider the convergence problem.

We illustrate the proposed method in Figure 1, and describe it in three steps:
In the first step, the noise is removed and the desired speech spectral
variance is estimated. Then, the output of WPE is obtained with the estimation
of the AR weights. Finally, the residual noise is removed by applying the
estimated mask received in the first step.

Figure 1: Schematic diagram of proposed method.

###  3.1 Denoising and Variance Estimation

A single channel DNN is used to predict the noiseless reverberation spectrum
and the desired speech spectral variance estimation simultaneously. The
structure of the DNN is shown in Figure 2. Eq.(11) shows that the desired
speech spectral variance $\sigma^{2}_{d_{n,k}}$ is the square of the desired
speech spectrum of the first channel (refer to Eq.(8)). We train a DNN to
predict the noiseless reverberation spectrum $X_{R}^{m}$ and desired speech
spectrum $X_{S}^{m}$ using ideal ratio mask $IRM^{m}_{R}$ and
$IRM^{m}_{S}$[23] as the targets which are defined as

$\begin{split}IRM^{m}_{R}&=\min\\{\frac{X^{m}_{R}}{X^{m}+\epsilon},1\\}\\\
IRM^{m}_{S}&=\min\\{\frac{X^{m}_{S}}{X^{m}+\epsilon},1\\}.\end{split}$ (13)

After obtaining the estimated mask, the estimated noiseless reverberation
spectrum $\hat{X}^{m}_{R}$ and the desired speech spectrum $\hat{X}^{m}_{S}$
are obtained by

$\begin{split}\hat{X}^{m}_{R}&=X^{m}\otimes IRM^{m}_{R}\\\
\hat{X}^{m}_{S}&=X^{m}\otimes IRM^{m}_{S},\end{split}$ (14)

where $\otimes$ denotes an element-wise multiplication.

Figure 2: Structure of the DNN. The inputs are the log spectra of the current
frame and its neighboring frames. The outputs are $IRM_{R}$ and $IRM_{S}$
respectively.

###  3.2 Dereverberation

After getting $\hat{X}_{R}$ and $\hat{X}^{m}_{S}$, the phase of mixture
($\angle{X}$) is used to reconstruct the noise-free phase spectrum
$\tilde{X}_{R}$. The WPE is adopted to reconstruct the noise-free desired
signal $\hat{d}_{n,k}$. $\tilde{X}_{R}$ and $|\hat{X}^{1}_{S}|^{2}$ are used
as the input of WPE and $\sigma^{2}_{d_{n,k}}$ in Eq.(11), respectively. The
AR weights, $\textbf{g}_{k}$, can be obtained from (12). Then the output of
WPE, $\hat{d}_{n,k}$, can be calculated from (8). The procedure is outlined in
Algorithm (1).

Algorithm 1 Dereverberation algorithm.

0: $L_{k}$, $D$ initialization:
$\sigma^{2}_{d_{n,k}}\leftarrow|\hat{X}^{1}_{S}|^{2}$,
$X\leftarrow\tilde{X}_{R}$

$\textbf{A}_{k}\leftarrow\sum_{n=1}^{N}\frac{\textbf{x}_{n-D,k}\textbf{x}^{H}_{n-D,k}}{\sigma^{2}_{d_{n,k}}}$

$\textbf{b}_{k}\leftarrow\sum_{n=1}^{N}\frac{\textbf{x}_{n-D,k}(x^{1}_{n,k})^{*}}{\sigma^{2}_{d_{n,k}}}$

$\hat{\textbf{g}}_{k}\leftarrow\textbf{A}^{-1}_{k}\textbf{b}_{k}$

$\hat{d}_{n,k}\leftarrow
x^{1}_{n,k}-\hat{\textbf{g}}^{H}_{k}\textbf{x}_{n-D,k}$

###  3.3 Post processing

The performance can be improved further by masking the $IRM^{1}_{S}$ to the
output $\hat{d}_{n,k}$ to form the final enhanced spectrum. Finally the
enhanced signal is converted from frequency domain to time domain using the
inverse short time Fourier transform (ISTFT).

##  4 Experiments

###  4.1 Experimental data and Metrics

The official set of the 2014 REVERB Challenge2[24] is used in the experiment.
The data set consists of a training, a development and a (final) evaluation
test set. The training and development sets are constructed using separated
parts of the WSJCAM0[25] corpus convolved with 24 RIRs and corrupted by
various types of noises at 10 dB. The reverberation time (RT) of the 24 RIRs
ranges roughly from 0.2 to 0.8 sec. Each RIR includes 8 channels. The test set
includes simulated and real world data sets. The simulated data is generated
by convolving WSJCAM0 corpus with 6 RIRs (3 rooms, 2 types of distance between
a speaker and a microphone array), and mixing with various types of noises at
10 dB. It should be mentioned that all the RIRs are recorded in real rooms and
different for the training and simulated sets.

We randomly select 5000 sentences (only using the first channel) from the
official training set to learn the DNN weights. 400 sentences are selected
from the official simulated data set to evaluate the proposed method.

A fixed 50-ms frame size was used with $80\%$ overlap between frames. The
discrete Fourier transform (DFT) is applied on each frame. The length of the
DFT is 800.

The performance is evaluated with perceptual evaluation of speech quality
(PESQ)[26] and cepstral distance (CD)[27] which reflect the quality of the
objective speech. The cepstral distance between two signals is defined as

$CD=\frac{10}{\log
10}\sqrt{(c_{0}-\hat{c}_{0})^{2}+2\sum_{k=1}^{12}(c_{k}-\hat{c}_{k})^{2}}\\\ $
(15)

where $c_{k}$ and $\hat{c}_{k}$ are the cepstral coefficients of the anechoic
speech signal and the estimated desired signal, respectively.

###  4.2 Experimental Setup and Parameter Selection

To evaluate the performance, the single channel DNN in proposed method uses
magnitude spectrogram with 5-frame context window (2 before and 2 after) as
the input features. The DNN has three hidden layers with 1024 rectified linear
units (ReLUs)[28] for each. Sigmoidal units are used in the output layer since
the IRM is bounded between 0 and 1. The input is normalized to zero mean and
unity variance over all feature vectors in the training set. The DNN weights
are randomly initialized without pretrain. Root Mean Square Propagation
(RMSprop)[29] is utilized for optimization. The order of the regression vector
$L_{k}$ and the prediction delay $D$ in WPE of the proposed method are set to
15 and 3 respectively.

Since the proposed algorithm (denoted as ‘Proposed’) is independent on the
number of microphones, we evaluate and compare the performances in single and
multi-channel scenario, separately. To evaluate the single-channel enhancement
performance, DNN-Han[16] and original WPE are used to compare with the
proposed method. DNN-Han use a DNN to learn the spectral mapping from
reverberant and noisy signal to clean signal directly. In multi-channel speech
enhancement experiment, we compare the proposed method with original WPE. In
both single- and multi-channel experiments, the $L_{k}$ and the $D$ in
original WPE are the same as the ones in the proposed method.

The original WPE cannot reduce the background noise. In order to compare
fairly, we apply the $IRM^{1}_{R}$, obtained in the first step, to the output
of the WPE for noise reduction. This is denoted as ‘WPE-mask’.

###  4.3 Experimental results

Figure 3 shows a comparison of the proposed method and the WPE for multi-
channel speech enhancement. “Ref” indicates the noise and reverberation
speech. The results show that the proposed method achieves best results on
PESQ and CD. The WPE-mask outperforms the WPE because of noise reduction.
However, it is still worse than the proposed method. It means that our method
has better dereverberation ability due to the noise reduction and precise
variance estimation by DNN.

The performance of single channel speech enhancement is shown in Figure 4. We
can observe that the proposed method is still better than WPE-mask. In
addition, the proposed method outperforms the DNN-Han. The possible reason is
that the DNN-Han is a totally data-driven model which does not consider the
statistical characteristics of speech.

The computational costs of the WPE and the proposed method are also evaluated.
In WPE, the best performance is achieved with 5 iterations. The comparison
result presented in Figure 5 shows that the computational effort has been
considerably reduced by eliminating the iterative process.

Finally, an example of reconstructed spectrogram is shown in Figure 6.
Obviously, our method has a better performance on denoising. As shown in the
black ellipse, our method can also suppress more reverberations.

Figure 3: PESQ and CD measures versus the numbers of iterations with various
source algorithms for multi-channel speech enhancement. The reported values
are obtained as the average of all test utterances.

Figure 4: PESQ and CD measures versus the numbers of iterations with various
source algorithms for single-channel speech enhancement. The reported values
are obtained as the average of all test utterances. Figure 5: Processing time
required from a 10 secs. speech for WPE and proposed method. Intel(R) Xeon(R)
CPU E5-2690 v2 @ 3.00GHz with RAM: 32G and Matlab environment. Figure 6: Top
left: The spectrogram of the mixture speech; Bottom left: Clean speech
spectrogram; Top right: Speech separation using WPE; Bottom right: Speech
separation using proposed method.

##  5 Conclusion

In this paper, a novel algorithm for speech denoising and dereverbration is
proposed. The experimental results show that the proposed method outperforms
the DNN and the conventional WPE. The proposed method takes advantage of the
merits of the machine learning method: DNN and the conventional speech
processing method: WPE. Using DNN, the influence of the background noise is
reduced in the process of WPE and more accurate parameters are obtained
directly. In addition, the proposed method can be used in single- or multi-
channel speech enhancement, which is very flexible in practice. The proposed
method can be easily extended to Recurrent Neural Networks (RNNs) and Long
Short Term Memory Networks (LSTMs) for better performance. We are also
investigating this avenue.

## References

  * [1] Oriol Vinyals, Suman V Ravuri, and Daniel Povey.  Revisiting recurrent neural networks for robust asr.  In Acoustics, Speech and Signal Processing (ICASSP), 2012 IEEE International Conference on, pages 4085–4088. IEEE, 2012.
  * [2] Gibak Kim, Yang Lu, Yi Hu, and Philipos C Loizou.  An algorithm that improves speech intelligibility in noise for normal-hearing listeners.  The Journal of the Acoustical Society of America, 126(3):1486–1494, 2009.
  * [3] Cemil Demir, Murat Saraclar, and Ali Taylan Cemgil.  Single-channel speech-music separation for robust asr with mixture models.  IEEE Transactions on Audio, Speech, and Language Processing, 21(4):725–736, 2013.
  * [4] Benjamin Cauchi, Patrick A Naylor, Timo Gerkmann, Simon Doclo, and Stefan Goetze.  Late reverberant spectral variance estimation using acoustic channel equalization.  In Signal Processing Conference (EUSIPCO), 2015 23rd European, pages 2481–2485. IEEE, 2015.
  * [5] Ina Kodrasi, Ante Jukić, and Simon Doclo.  Robust sparsity-promoting acoustic multi-channel equalization for speech dereverberation.  In Acoustics, Speech and Signal Processing (ICASSP), 2016 IEEE International Conference on, pages 166–170. IEEE, 2016.
  * [6] Takuya Yoshioka, Tomohiro Nakatani, and Masato Miyoshi.  Integrated speech enhancement method using noise suppression and dereverberation.  IEEE Transactions on Audio, Speech, and Language Processing, 17(2):231–246, 2009.
  * [7] Emanuël AP Habets, Sharon Gannot, and Israel Cohen.  Late reverberant spectral variance estimation based on a statistical model.  IEEE Signal Processing Letters, 16(9):770–773, 2009.
  * [8] Emanuel AP Habets.  Multi-channel speech dereverberation based on a statistical model of late reverberation.  In Acoustics, Speech, and Signal Processing, 2005. Proceedings.(ICASSP’05). IEEE International Conference on, volume 4, pages iv–173. IEEE, 2005.
  * [9] Tomohiro Nakatani, Biing-Hwang Juang, Takuya Yoshioka, Keisuke Kinoshita, Marc Delcroix, and Masato Miyoshi.  Speech dereverberation based on maximum-likelihood estimation with time-varying gaussian source model.  IEEE transactions on audio, speech, and language processing, 16(8):1512–1527, 2008.
  * [10] Alfred Mertins, Tiemin Mei, and Markus Kallinger.  Room impulse response shortening/reshaping with infinity-and $p$-norm optimization.  IEEE Transactions on Audio, Speech, and Language Processing, 18(2):249–259, 2010.
  * [11] Ina Kodrasi, Stefan Goetze, and Simon Doclo.  Regularization for partial multichannel equalization for speech dereverberation.  IEEE Transactions on Audio, Speech, and Language Processing, 21(9):1879–1890, 2013.
  * [12] Tomohiro Nakatani, Takuya Yoshioka, Keisuke Kinoshita, Masato Miyoshi, and Biing-Hwang Juang.  Blind speech dereverberation with multi-channel linear prediction based on short time fourier transform representation.  In Acoustics, Speech and Signal Processing, 2008. ICASSP 2008. IEEE International Conference on, pages 85–88. IEEE, 2008.
  * [13] Tomohiro Nakatani, Takuya Yoshioka, Keisuke Kinoshita, Masato Miyoshi, and Biing-Hwang Juang.  Speech dereverberation based on variance-normalized delayed linear prediction.  IEEE Transactions on Audio, Speech, and Language Processing, 18(7):1717–1731, 2010.
  * [14] Mahdi Parchami, Wei-Ping Zhu, and Benoit Champagne.  Speech dereverberation using linear prediction with estimation of early speech spectral variance.  pages 504–508, 2016.
  * [15] Xiao-Lei Zhang and DeLiang Wang.  A deep ensemble learning method for monaural speech separation.  IEEE/ACM Transactions on Audio, Speech and Language Processing (TASLP), 24(5):967–977, 2016.
  * [16] Kun Han, Yuxuan Wang, DeLiang Wang, William S Woods, Ivo Merks, and Tao Zhang.  Learning spectral mapping for speech dereverberation and denoising.  IEEE Transactions on Audio, Speech, and Language Processing, 23(6):982–992, 2015.
  * [17] Hao Li, Shuai Nie, Xueliang Zhang, and Hui Zhang.  Jointly optimizing activation coefficients of convolutive nmf using dnn for speech separation.  In Nelson Morgan, editor, INTERSPEECH, pages 550–554. ISCA, 2016\.
  * [18] Xueliang Zhang, Hui Zhang, Shuai Nie, Guanglai Gao, and Wenju Liu.  A pairwise algorithm using the deep stacking network for speech separation and pitch estimation.  IEEE/ACM Transactions on Audio, Speech, and Language Processing, 24(6):1066–1078, 2016.
  * [19] Yuxuan Wang, Arun Narayanan, and DeLiang Wang.  On training targets for supervised speech separation.  IEEE/ACM Transactions on Audio, Speech and Language Processing (TASLP), 22(12):1849–1858, 2014.
  * [20] Yong Xu, Jun Du, Li-Rong Dai, and Chin-Hui Lee.  An experimental study on speech enhancement based on deep neural networks.  IEEE Signal processing letters, 21(1):65–68, 2014.
  * [21] Arun Narayanan and DeLiang Wang.  Ideal ratio mask estimation using deep neural networks for robust speech recognition.  In Acoustics, Speech and Signal Processing (ICASSP), 2013 IEEE International Conference on, pages 7092–7096. IEEE, 2013.
  * [22] Israel Cohen, Jacob Benesty, Sharon Gannot, and Kommunikationstechnik.  Speech processing in modern communication.  Journal of the Acoustical Society of America, 129(1):535–535, 2010\.
  * [23] Soundararajan Srinivasan, Nicoleta Roman, and DeLiang Wang.  Binary and ratio time-frequency masks for robust speech recognition.  Speech Communication, 48(11):1486–1501, 2006.
  * [24] Keisuke Kinoshita, Marc Delcroix, Sharon Gannot, Emanuel A. P. Habets, Reinhold Haeb-Umbach, Walter Kellermann, Volker Leutnant, Roland Maas, Tomohiro Nakatani, and Bhiksha Raj.  A summary of the reverb challenge: state-of-the-art and remaining challenges in reverberant speech processing research.  Journal on Advances in Signal Processing, 2016(1):1–19, 2016.
  * [25] T. Robinson, J. Framsen, D. Pye, J. Foote, and S. Renals.  WSJCAM0: A British English speech corpus for large vocabulary continuous speech recognition.  In Proceedings IEEE International Conference on Acoustics, Speech and Signal Processing, 1995.
  * [26] Antony W Rix, John G Beerends, Michael P Hollier, and Andries P Hekstra.  Perceptual evaluation of speech quality (pesq)-a new method for speech quality assessment of telephone networks and codecs.  In Acoustics, Speech, and Signal Processing, 2001. Proceedings.(ICASSP’01). 2001 IEEE International Conference on, volume 2, pages 749–752. IEEE, 2001.
  * [27] Sadaoki Furui.  Digital Speech Processing, Synthesis, and Recognition (Revised and Expanded).  Marcel Dekker, Inc., 2000.
  * [28] Xavier Glorot, Antoine Bordes, and Yoshua Bengio.  Deep sparse rectifier neural networks.  15(106):275, 2011.
  * [29] Tijmen Tieleman and Geoffrey Hinton.  Lecture 6.5-rmsprop: Divide the gradient by a running average of its recent magnitude.  COURSERA: Neural networks for machine learning, 4(2), 2012.