# Partial and complete qubit estimation using a single observable:
optimization and quantum simulation

Cristian A. Galvis Florez, J. Martínez-Cifuentes, K. M. Fonseca-Romero∗
∗Departamento de Física, Universidad Nacional de Colombia - Sede Bogotá,
Facultad de Ciencias, Grupo de Óptica e Información Cuántica, Carrera 30 Calle
45-03, C.P. 111321, Bogotá, Colombia<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Quantum state estimation is an important task of many quantum information
protocols. We consider two families of unitary evolution operators, one with a
one-parameter and the other with a two-parameter, which enable the estimation
of a single spin component and all spin components, respectively, of a two-
level quantum system. To evaluate the tomographic performance, we use the
quantum tomographic transfer function (qTTF), which is calculated as the
average over all pure states of the trace of the inverse of the Fisher
information matrix. Our goal is to optimize the qTTF for both estimation
models. We find that the minimum qTTF for the one-parameter model is achieved
when the entangling power of the corresponding unitary operator is at its
maximum. The models were implemented on an IBM quantum processing unit, and
while the estimation of a single-spin component was successful, the whole spin
estimation displayed relatively large errors due to the depth of the
associated circuit. To address this issue, we propose a new scalable circuit
design that improves qubit state tomography when run on an IBM quantum
processing unit.

  *

Keywords: Quantum state estimation, Quantum tomographic transfer function,
Quantum Computing, IBM Quantum Experience.

##  1 Introduction

The concept of quantum mechanical computers was first proposed in the works of
Feynman [1, 2]. Soon, several quantum algorithms were developed that are more
efficient than their classical counterparts [3, 4]. These algorithms have been
implemented in current quantum computers [5, 6] using a few qubits. One of the
key tasks in assessing the performance of these quantum algorithms is the
estimation of the computer’s quantum state [7, 8]. This is particularly
important in quantum cryptography, where the qubit state needs to be
accurately specified both at the source and after transmission [9]. Due to
errors in state preparation and processing, it is essential to have accurate
state estimates at various stages of quantum protocols. Various methods for
state estimation [10] have been implemented, such as maximal entropy [11, 12],
maximum likelihood [13], Bayesian methods [14], and linear regression
estimation [15]. These methods may vary in their resource requirements, such
as the number of measurements needed or the knowledge of the system’s
evolution. For example, usual tomography assume access to $n^{2}$ different
measurements; adaptive methods, the possibility to measure arbitrary
observables; stroboscopic methods [16, 17, 18] require the knowledge of the
evolution of the system. Like the single-observable state estimation [19],
stroboscopic tomography requires only a single observable. However, single-
observable tomography also requires an auxiliary system of dimension at least
equal to that of the system whose state is estimated and the joint evolution
at a single instant of time.

Peres’ model [20] , a particular instance of the von Neumann model of
measurement, involves an auxiliary two-level system to measure a qubit
observable (sec. 2.1). Its extension [21], which employs two auxiliary two-
level systems (sec. 2.2), is an example of a variant of the single-observable
tomographic method, in which the system is not directly measured. This would
be the case of cavity QED, for example, where two-level atoms are used as
probes to measure the state of the confined photons.

Peres’ model and its extension are parameterized by one or two parameters,
respectively. Here, we consider and solve the problem of finding the
parameters which optimize the tomographic performance of these models, as
measured by the quantum tomographic transfer function (qTTF) [22]. The qTTF is
the average, over all pure states of the system, of the trace of the inverse
of the Fisher information matrix. The models considered in this work, are run
on a five-qubit IBM quantum processing unit (QPU). The experimental qTTF is
found by averaging the error over a set of six points, equally spaced on the
Bloch sphere. This experimental quantifier is satisfactory for Peres’ model
and for an alternative model of qubit state estimation, but it is large for
the circuit implementation of the extended Peres’ model (sec. 3) . The
alternative model, which employs only two CNOT gates, is proposed, optimized
and run on an IBM QPU (subsection 3.3). An alternative proposal for multiqubit
states, implemented on IBM quantum processors, use entangled bases [23].
Scalability of our model is discussed in the last section of this paper
(section 4).

##  2 Partial and complete spin estimation methods

###  2.1 Estimation of a single spin component

The spin component $s_{z}$ of a qubit $S$ can be estimated by using an
auxiliary qubit (meter) $A$. The interaction between the system and meter is
described by the Hamiltonian
$H_{A}=g\ket{1^{S}}\bra{1^{S}}\otimes\ket{1^{A}}\bra{1^{A}}$ [20], where $g$
is a coupling constant. The associated evolution operator is

$U_{A}(\theta)=e^{-iH_{A}T/\hbar}=\begin{pmatrix}1&0&0&0\\\ 0&1&0&0\\\
0&0&1&0\\\ 0&0&0&e^{-i\theta}\\\ \end{pmatrix}.$ (1)

$U(\theta),$ which has been written in two-qubit basis in the lexicographic
order $\ket{00},\ket{01},\ket{10},\ket{11},$ has the form of a controlled
phase gate ($C_{P}(\theta)$). The evolution operator depends on the single
parameter $\theta=gT/\hbar,$ where $T$ represents the interaction time. Due to
the periodicity of $U_{A}(\theta)$, we consider $\theta\in[0,2\pi).$ The
initial state of the qubit-meter system is assumed to be
$\ket{\psi^{S}_{0}}\otimes\ket{+^{A}}$, where the initial states of system and
meter are, respectively, $\ket{\psi^{S}_{0}}=c_{0}\ket{0}+c_{1}\ket{1},$ and
$\ket{+^{A}}=\frac{1}{\sqrt{2}}\left(\ket{0}+\ket{1}\right).$ Coefficients
$c_{0}$ and $c_{1}$ satisfy the normalization condition
$|c_{0}|^{2}+|c_{1}|^{2}=1.$ At time $T$, the meter is measured in the
$\sigma_{x}$ basis. It is detected in the state $\ket{+}$
($\sigma_{x}\ket{+}=\ket{+}$) with probability $P_{+}=P_{0}$ or in the state
$\ket{-}$ ($\sigma_{x}\ket{-}=-\ket{-}$) with probability $P_{-}=P_{1}$. Here,

$\displaystyle P_{0}$
$\displaystyle=|c_{0}|^{2}+|c_{1}|^{2}\cos^{2}\frac{\theta}{2}=\frac{1}{2}\left(1+\cos^{2}\frac{\theta}{2}\right)+\frac{1}{2}\sin^{2}\frac{\theta}{2}s_{z},$
(2) $\displaystyle P_{1}$
$\displaystyle=|c_{1}|^{2}\sin^{2}\frac{\theta}{2}=\frac{1}{2}\sin^{2}\frac{\theta}{2}\left(1-s_{z}\right),$
(3)

and
$s_{z}=\braket{\psi_{0}^{S}}{\sigma_{z}^{S}}{\psi_{0}^{S}}=|c_{0}|^{2}-|c_{1}|^{2}.$
In terms of the probabilities of detection, the parameter $s_{z}$ can be
estimated using the relation

$s_{z}=\csc^{2}\left(\frac{\theta}{2}\right)\left(P_{0}-P_{1}-\cos^{2}\left(\frac{\theta}{2}\right)\right).$
(4)

$\large\Qcircuit@C=1.0em@R=0.3em@!R{\lstick{{A}:}&\gate{H}\barrier[0em]{1}\qw\ctrl{1}\barrier[0em]{1}\qw\gate{H}\barrier[0em]{1}\qw\meter\qw\qw\\\
\lstick{{S}:}\gate{U_{\xi}}\qw\gate{P_{\theta}}\qw\qw\qw\qw\qw\qw\\\
\lstick{c:}\cw\cw\cw\cw\cw\cw\dstick{}\cw\cwx[-2]\cw\cw\\\ }$ Figure 1:
Quantum circuit for the estimation of a single spin component of a qubit $S$
using the auxiliary qubit $A$.

The circuit that represents this process, shown in Figure 1, is divided into
four sections: preparation, interaction, change of basis, and measurement.
Taking into account that the qubits’ default state in quantum circuits is
$\ket{0},$ the initial state of the qubit-meter system is
$U_{\xi}\ket{0^{S}}\otimes U_{H}\ket{0^{A}},$ where $U_{\xi}$ denotes a
general single-qubit gate rotation defined by a set of angles $\xi$ such that
$U_{\xi}\ket{0^{S}}=\ket{\psi_{0}^{S}},$ and $U_{H}$ a Hadamard gate.

The evolution operator (1), is implemented by a $C_{P}(\theta)$ gate. Since
only measurements on the $\sigma_{z}$ basis are available, it is necessary to
perform a change of basis (the Hadamard gate of the third section) to measure
the auxiliary qubit in the $\sigma_{x}$ basis.

It is essential to find an estimation for $s_{z},$ and to characterize the
associated error (the square root of the variance). Any unbiased estimator
satisfies the Cramér-Rao bound
Cov$(\mathbf{s},\mathbf{s})\geq\mathbb{F}^{-1}(\mathbf{s})$, where
$\mathbf{s}$ is the vector of parameters to be estimated,
Cov$(\mathbf{s},\mathbf{s})$ the corresponding covariance matrix, and
$\mathbb{F}(\mathbf{s})$ the associated Fisher information matrix. The
elements of the Fisher matrix are given by [24]

$F_{\mu\nu}(\mathbf{s})=-\mathbb{E}\left(\frac{\partial^{2}\ln
P(\mathbf{s})}{\partial s_{\mu}\partial s_{\nu}}\right),$

were $\mathbb{E}$ denotes the statistical expectation value. We consider that
the multinomial distribution
$P(\mathbf{n}|\mathbf{s})=(\sum_{k}n_{k})!\prod_{k}p_{k}^{n_{k}}/(n_{k}!),$
gives the probability that $n_{k}$ events corresponding to the $k$-th outcome
were detected. Assuming the probabilities $p_{k}$ to linearly depend on the
parameters $s_{\mu}$ we get

$F_{\mu\nu}(\mathbf{s})=\mathbb{E}\left(\sum_{k}\frac{n_{k}}{p_{k}^{2}}\frac{\partial
p_{k}}{\partial s_{\mu}}\frac{\partial p_{k}}{\partial s_{\nu}}\right).$

Using $n_{k}=Np_{k},$ where $N=\sum_{k}n_{k},$ and taking into account that
the resulting expression is not longer a function of $n_{k},$ we obtain

$F_{\mu\nu}(\mathbf{s})=N\sum_{j=0}^{1}\frac{1}{p_{j}}\frac{\partial
p_{j}}{\partial s_{\mu}}\frac{\partial p_{j}}{\partial s_{\nu}}.$ (5)

The $N$ factor, which describes a variance proportional to $N^{-1}$, will be
ignored, unless we explicitly consider the effect of the sample size. Here,
$s_{z}$ is the only parameter to be estimated. Hence, the Fisher matrix
becomes just a number, whose inverse is

$F^{-1}(c_{0},c_{1},\theta)=\frac{4P_{0}(c_{0},c_{1})P_{1}(c_{0},c_{1})}{\sin^{4}\frac{\theta}{2}}.$
(6)

Since $F^{-1}(c_{0},c_{1},\theta)$ is the lower limit to the variance
associated to $s_{z}$, it is defined as the estimation error. In A.1, we show
that, in the limit $N\rightarrow\infty$, the variance of the average of
$s_{z}$ over $N$ shots coincides with the Fisher error

$\operatorname{Var}(\hat{s}_{z})\rightarrow\frac{1}{N-1}F^{-1}\rightarrow\frac{1}{N}F^{-1}.$
(7)

To study the dependence of the error on the initial state of qubit $S$, we
parameterize the state of the system using $\xi=\\{\alpha_{1},\alpha_{2}\\}$
as $c_{0}=e^{i\alpha_{2}}\cos\alpha_{1}$ and
$c_{1}=e^{-i\alpha_{2}}\sin\alpha_{1}$, where $\alpha_{1}\in[0,\pi/2],$ and
$\alpha_{2}\in[0,\pi]$. The resulting expression for the error,

$\tilde{F}^{-1}(\alpha_{1},\alpha_{2},\theta)=F^{-1}(c_{0}(\alpha_{1},\alpha_{2}),c_{1}(\alpha_{1},\alpha_{2}),\theta)=4\sin^{2}(\alpha_{1})(\csc^{2}(\theta/2)-\sin^{2}(\alpha_{1})),$

is independent of $\alpha_{2}.$ This error is plotted in Figure 2 for
different values of $\theta$ in the interval $[\pi/2,\pi]$. Values of $\theta$
between 0 and $\frac{\pi}{2}$, not depicted in Figure 2, have larger
estimation errors. It can be seen that, for a given fixed initial state,
$\theta=\pi$ produces the least error. Hence, for the interaction considered
in this work, it is convenient to use $\theta=\pi$ to obtain the best
estimation of $s_{z}.$

The maximum error, for fixed $\theta$, satisfy

$\frac{\partial
F^{-1}}{\partial\alpha_{1}}=4\sin(2\alpha_{1})\left(\csc^{2}(\theta/2)-2\sin^{2}\alpha_{1}\right)=0.$
(8)

If $0\leq\theta\leq\pi/2$ ($\pi/2\leq\theta\leq\pi$) the maximum error occurs
at $\alpha_{1}=\pi/2$ ($\alpha_{1}=\arcsin(\frac{1}{\sqrt{2}}\csc(\theta/2))$)
where $F^{-1}=4(\csc^{2}(\theta/2)-1)$ ($F^{-1}=\csc^{4}(\theta/2)$). Hence,
the maximum error for constant $\theta$ is
$F^{-1}_{max}=\max(4(\csc^{2}(\theta/2)-1),\csc^{4}(\theta/2)).$

Since the system-meter state is pure, its entanglement can be quantified by
the entropy of either reduced density matrix,
$\rho_{S}=\operatorname{tr}_{A}\left(\ket{\psi_{AS}}\bra{\psi_{AS}}\right),$
or
$\rho_{A}=\operatorname{tr}_{S}\left(\ket{\psi_{AS}}\bra{\psi_{AS}}\right).$
Here
$\ket{\psi_{AS}}=(1/\sqrt{2})\left(c_{0}\ket{00}+c_{0}\ket{01}+c_{1}\ket{10}+c_{1}e^{-i\theta}\ket{11}\right).$
The von Neumann entropy [25] of the system, after the interaction with the
meter, is

$E(\rho_{S})=-\operatorname{Tr}\left(\rho_{S}\log\rho_{S}\right)=-\lambda_{+}\log\lambda_{+}-\lambda_{-}\log\lambda_{-},$
(9)

where $\lambda_{\pm}=(1/2)(1\pm\sqrt{x})$ and
$x={1-\sin^{2}(2\alpha_{1})\sin^{2}(\theta/2)}$. System’s entropy
$E(\rho_{S})$ is not only qualitatively similar to the error $F^{-1}$, it also
vanishes for the same combinations of $\alpha_{1}$ and $\theta:$
$(\alpha_{1}=0,\textrm{any }\theta)$, and $(\alpha_{1}=\pi/2,\theta=\pi).$ The
initial state corresponding to $\alpha_{1}=0$, $\ket{0^{S}}\ket{+^{A}}$, is
left unchanged by the evolution operator (1). On the other hand, the initial
state when $\alpha_{1}=\frac{\pi}{2}$ and $\theta=\pi,$
$\ket{1^{S}}\ket{+^{A}},$ changes to $\ket{1^{S}}\ket{-^{A}}.$ The error
vanishes in both cases because only one outcome is available.

Figure 2: a) Error $F^{-1}$ as a function of $\alpha_{1}$ for
$\theta=\frac{\pi}{2},\frac{2\pi}{3},\frac{3\pi}{4}$ and $\pi$. Error
decreases as $\theta\rightarrow\pi$. b) Entropy of the system state, after the
interaction, as a function of $\alpha_{1}$ for
$\theta=\frac{\pi}{2},\frac{2\pi}{3},\frac{3\pi}{4}$ and $\pi$. The entropy
and the error vanish for the same values of $\alpha_{1}$ and $\theta$.

To study the dependence of the error of the estimation only on the coupling
parameter $\theta$, we define the average of
$F^{-1}(\alpha_{1},\alpha_{2},\theta)$ (we have dropped the tilde to lighten
the notation) over the set of initial pure states of the system,

$\overline{F^{-1}}(\theta)=\frac{1}{\pi}\int_{0}^{\pi/2}d\alpha_{1}\int_{0}^{\pi}d\alpha_{2}\,F^{-1}(\alpha_{1},\alpha_{2},\theta)\sin(2\alpha_{1}).$
(10)

This average is equivalent to the natural generalization of the quantum
tomographic transfer function (qTTF), first defined in [22], for the case
where only one parameter is estimated. Optimization of the qTTF over $\theta$
allows us to determine which $U(\theta)$ performs better for the estimation of
$s_{z}$. We obtain
$\overline{F^{-1}}(\theta)=\frac{2}{3}(2+\cos\theta)\csc^{2}(\theta/2),$ which
is a monotonic decreasing function in $[0,\pi].$ $\overline{F^{-1}}(\theta)$
diverges for $\theta\to 0$ and attains its minimum at $\theta=\pi,$
$\overline{F^{-1}}(\pi)=2/3.$

Figure 3: qTTF $\overline{F^{-1}}$ (solid line) and maximum error
$F_{max}^{-1}$ (dashed line) as functions of the coupling parameter $\theta$.

Figure 3 shows $\overline{F^{-1}}(\theta)$ and the maximum error for fixed
$\theta$, $F^{-1}_{max}$, as a function of $\theta.$ As can be seen, the value
$\theta=\pi$ corresponds to the coupling value with the minimum average error,
which agrees with the results shown in Figure 2, where the configuration with
$\theta=\pi$ has the lowest values of the error
$F^{-1}(\alpha_{1},\alpha_{2},\theta)$. This value coincides with the maximum
entangling power of the unitary evolution $U$ [26].

###  2.2 Complete qubit state estimation

The complete initial state of qubit $S$ can be estimated using two auxiliary
meters $A$ and $B$. We consider the interaction Hamiltonian

$H_{AB}=\left(\lambda_{A}\Pi^{S}_{1}\otimes\Pi^{A}_{1}\otimes
I^{B}+\lambda_{B}\Pi^{S}_{+}\otimes I^{A}\otimes\Pi^{B}_{1}\right),$ (11)

where $\Pi_{k}^{M}=|k^{M}\rangle\langle k^{M}|$, with $k=1,+$ and $M=A,B$.
Here, $I_{A}$ and $I_{B}$ are the identity operators associated to the Hilbert
spaces of meters $A$ and $B$, respectively, and $\lambda_{A}$, $\lambda_{B}$
are constants. The eigenstates of $\sigma_{x}$ are denoted by $\ket{+}$ and
$\ket{-}$. The time evolution operator at time $T$,
$U_{AB}(\theta_{A},\theta_{B}),$

$U_{AB}(\theta_{A},\theta_{B})=e^{-iH_{AB}T/\hbar}=\exp\left(-i\theta_{A}\Pi^{S}_{1}\otimes\Pi^{A}_{1}\otimes
I^{B}-i\theta_{B}\Pi^{S}_{+}\otimes I^{A}\otimes\Pi^{B}_{1}\right)$ (12)

depends on the parameters $\theta_{A}=(\lambda_{A}T/\hbar)$ and
$\theta_{B}=(\lambda_{B}T/\hbar)$.

As in the previous section, the $A$ and $B$ meters are prepared in the state
$\ket{+}=\left(\ket{0}+\ket{1}\right)/\sqrt{2}$ and qubit $S$ is prepared in
an arbitrary state $\ket{\psi_{0}^{S}}=c_{0}\ket{0}+c_{1}\ket{1}$. The joint
initial state
$|\Psi(0)\rangle=\ket{\psi_{0}^{S}}\otimes|+^{A}\rangle\otimes|+^{B}\rangle$
is transformed by the interaction into the final state

$\ket{\Psi_{f}}=U_{AB}(\theta_{A},\theta_{B})\ket{\Psi(0)}=\frac{1}{2}\sum_{ij=0,1}U_{ij}|\psi_{0}^{S}\rangle\otimes|i^{A},j^{B}\rangle,$
(13)

where $|i^{A},j^{B}\rangle=|i^{A}\rangle\otimes|j^{B}\rangle$ and

$\displaystyle U_{00}$ $\displaystyle=I,$ $\displaystyle U_{01}$
$\displaystyle=e^{-i\theta_{B}\Pi_{+}^{S}},$ (14) $\displaystyle U_{10}$
$\displaystyle=e^{-i\theta_{A}\Pi_{1}^{S}},$ $\displaystyle U_{11}$
$\displaystyle=e^{-\theta_{A}\Pi_{1}^{S}-i\theta_{B}\Pi_{+}^{S}}.$ (15)

The probability $p_{kl}$ that meters A and B are detected in states
$\ket{k^{A}}$ and $\ket{l^{A}}$, respectively, $k,l=\pm,$ is

$p_{kl}=\frac{1}{4}s_{0}+\sum_{\mu=1}^{3}\left(a_{\mu}k+b_{\mu}l+c_{\mu}kl\right)s_{\mu}.$
(16)

Here, $s_{0}=1$ and ${s_{\mu}}$, $\mu=1,2,3$, are the components of the Bloch
vector of the initial state of the system
$\rho_{0}^{S}=\ket{\psi_{0}}\bra{\psi_{0}}$,
$s_{\mu}=\operatorname{tr}(\rho_{0}^{S}\sigma_{\mu})$. The matrix $\sigma_{0}$
is the identity matrix, and $\sigma_{i},$ with $i=1,2,3$ are the Pauli
matrices. The coefficients $a_{\mu}$, $b_{\mu},$ and $c_{\mu}$ are given by
$a_{\mu}=\frac{1}{2}\operatorname{Tr}\left[A\sigma_{\mu}\right]$,
$b_{\mu}=\frac{1}{2}\operatorname{Tr}\left[B\sigma_{\mu}\right]$ and
$c_{\mu}=\frac{1}{2}\operatorname{Tr}\left[C\sigma_{\mu}\right]$, where

$\displaystyle A$
$\displaystyle=\frac{1}{16}\sum_{ij}U_{ij}^{\dagger}U_{1-i,j},\qquad
B=\frac{1}{16}\sum_{ij}U_{ij}^{\dagger}U_{i,1-j},\qquad
C=\frac{1}{16}\sum_{ij}U_{ij}^{\dagger}U_{1-i,1-j}.$ (17)

The explicit expressions for coefficients $a_{\mu}$, $b_{\mu},$ and $c_{\mu}$
can be found in the B. Although we have assumed that qubit $S$ is prepared in
a pure initial state, equation (16) holds for pure and mixed initial states
$\rho_{0}^{S}$.

By defining the probability vector
$\textbf{P}=(p_{00},p_{01},p_{10},p_{11})^{T}$, and the state vector
$\textbf{S}=(s_{0},s_{1},s_{2},s_{3})^{T}$, we can write equation (16) as
$\textbf{P}=\mathbb{T}\textbf{S}$, where $\mathbb{T}$ is an implicitly defined
matrix. If this matrix is invertible, we can estimate vector S, associated
with the initial state of the system, as

$\textbf{S}=\mathbb{T}^{-1}\textbf{P}.$ (18)

This equation lets us estimate the vector S using the probabilities
(frequencies) measured in meters $A$ and $B$.

The linear inversion estimation, which we have just described, sometimes
predicts non-positive semidefinite states. The discrete maximum-likelihood
estimator [27, 28], which can also be expressed in terms of the matrix
$\mathbb{T}$ and the probabilities P, always predicts physical states. This
estimator is computed iteratively using the $R\rho R$ Algorithm [29]. The
$\rho^{(n+1)}$ estimation is computed from the the previous estimation,
$\rho^{(n)},$ estimation

${\rho}^{(n+1)}=\mathcal{N}\left[{R}\left({\rho}^{(n)}\right){\rho}^{(n)}{R}\left({\rho}^{(n)}\right)\right],$
(19)

where $\mathcal{N}[\cdot]$ denotes normalization to trace unity of the
corresponding operator. Here,

${R}=\sum_{\mu=0}^{3}{r}_{\mu}\sigma_{\mu}\qquad\text{and}\qquad
r_{\mu}^{(n)}=\sum_{q=1}^{4}\frac{P_{q}}{\check{P}_{q}^{(n)}}T_{q\mu}.$

The components of the vector of experimental frequencies P are denoted by
$P_{q},$ and $\check{P}_{q}^{(n)}=\sum_{\mu=0}^{3}T_{q\mu}s_{\mu}^{(n)}$,
where $s_{\mu}^{(n)}$ are the components of $\textbf{S}^{(n)},$ the Bloch
components of ${\rho}^{(n)}.$ The initial state, $\rho^{(0)}=\sigma_{0}/2$ is
the maximally mixed state.

To simulate the quantum state estimation we propose the three-qubit circuit
shown in Figure 4. One of the qubits corresponds to the system $S$, whose
initial state is going to be estimated; the other two correspond to the $A$
and $B$ meters.

${\large\Qcircuit@C=1.0em@R=0.2em@!R{\lstick{{A}:}&\gate{H}\barrier[0em]{2}\qw\multigate{2}{U_{AB}\
\textrm{operator}}\barrier[0em]{2}\qw\gate{H}\barrier[0em]{2}\qw\meter\qw\qw\qw\\\
\lstick{{S}:}\gate{U_{\xi}}\qw\ghost{U_{AB}\
\textrm{operator}}\qw\qw\qw\qw\qw\qw\qw\\\
\lstick{{B}:}\gate{H}\qw\ghost{U_{AB}\
\textrm{operator}}\qw\gate{H}\qw\qw\meter\qw\qw\\\
\lstick{c:}\cw\cw\cw\cw\cw\cw\cw\cwx[-3]\cw\cwx[-1]\cw\cw\\\ \vspace{2em}}}$
Figure 4:  General structure of the quantum circuit used for the quantum state
estimation.

The first part of the circuit corresponds to the initialization of the qubits,
using a Hadamard gate for each of the meters, and a general $U_{\xi}$ single-
qubit gate rotation for the system. The second part is the core of the
algorithm. The actual implementation of $U_{AB}(\theta_{A},\theta_{B})$ was
found employing the transpilation algorithm of IBM-Q.

The Hadamard gates, at the end of the circuit, perform the transformation to
the $\sigma_{x}$ basis of the meters.

The Fisher information matrix for this model,

$\mathbb{F}=\mathbb{D}^{T}\left(\mathbb{V}\mathbb{P}^{-1}\mathbb{V}^{T}\right)\mathbb{D},$
(20)

depends on the measurement matrix $\mathbb{D}$, and the auxiliary matrices
$\mathbb{V}$ and $\mathbb{P}$

$\mathbb{D}=\begin{pmatrix}a_{1}&b_{1}&c_{1}\\\ a_{2}&b_{2}&c_{2}\\\
a_{3}&b_{3}&c_{3}\end{pmatrix},\quad\mathbb{V}=\begin{pmatrix}1&1&-1&-1\\\
1&-1&1&-1\\\
1&-1&-1&1\end{pmatrix},\quad\mathbb{P}=\begin{pmatrix}p_{00}&0&0&0\\\
0&p_{01}&0&0\\\ 0&0&p_{10}&0\\\ 0&0&0&p_{11}\end{pmatrix}.$

To obtain this result, it is convenient to write the equation
$\mathbf{P}=\mathbb{T}\mathbf{S}$ as
$\mathbf{P}=(\frac{1}{4}\mathbf{1}+\mathbb{V}^{T}\mathbf{D}_{0})s_{0}+\mathbb{V}^{T}\mathbb{D}\mathbf{s},$
where $\mathbf{1}=(1,1,1,1)^{T},\mathbf{D}_{0}=(a_{0},b_{0},c_{0})^{T}$, and
$\mathbf{s}=(s_{1},s_{2},s_{3})^{T}.$ Direct use of (5) gives
$F_{\mu\nu}=\sum_{k=1}^{4}\frac{1}{p_{k}}(\mathbb{V}^{T}\mathbb{D})_{k\mu}(\mathbb{V}^{T}\mathbb{D})_{k\nu},$
which can be written in matrix form as (20).

We define the error of the estimator as the trace of the inverse of the Fisher
information matrix

$\Delta(\alpha_{1},\alpha_{2},\theta_{A},\theta_{B})=\operatorname{Tr}(\mathbb{F}^{-1}).$
(21)

In the limit of large $N$, the variance of the estimator ($\hat{s}^{i}$) of
the $i$-th component of the Bloch vector is related to an element of the
Fisher error matrix

$\operatorname{Var}\left(\hat{{s}}^{i}\right)=\frac{1}{N}(\mathbb{F}^{-1})_{ii},$

as shown in A.2.

The mean error over the whole set of initial pure states of the system, the
qTTF, is used to characterize the error of the setup. It depends only on the
parameters $\theta_{A}$ and $\theta_{B}.$ As mentioned in Section 2.1, an
optimization of the average error over $\theta_{A}$ and $\theta_{B}$ allows
the determination of which $U_{AB}(\theta_{A},\theta_{B})$ performs better for
the estimation of the initial state of qubit $S$. Using the scipy.optimize
package [30], we found that $\theta_{A}^{\textrm{min}}\approx 3.45$ and
$\theta_{B}^{\textrm{min}}\approx-8.42$ minimize the mean error, which takes
the value $\bar{\Delta}(\theta_{A},\theta_{B})\approx 17.0$.

The error
$\Delta(\alpha_{1},\alpha_{2},\theta_{A}^{\textrm{min}},\theta_{B}^{\textrm{min}})$,
shown in Figure (5), displays moderate fluctuations as a function of the
initial state.

Figure 5: Fisher error,
$\Delta(\alpha_{1},\alpha_{2},\theta_{A}^{\textrm{min}}\approx
3.45,\theta_{B}^{\textrm{min}}\approx-8.42)$ as a function of the initial
state (parameterized by $\alpha_{1}$ and $\alpha_{2}$).

##  3 Results and Discussion

###  3.1  Estimation of a single spin component

We first simulate the estimation of $s_{z}$ with the circuit in Figure 1. The
results of the simulation for different values of $\theta$ with
$c_{0}=\frac{1}{2}$ and $c_{1}=\frac{\sqrt{3}}{2}$ are shown in Figure 6.

Figure 6: a) Ideal simulation of the circuit in figure 1 in the
qasm_simulator, performed in a local (not quantum) simulator. b) Simulation in
the ibmq_lima QPU exhibits some bias with respect to the theoretical result,
indicating the presence of hardware errors that affect the measurement
probabilities. c) Simulation using the qiskit measurement error mitigator does
not significantly improve error.

The computer simulation, made with $1024$ shots, displays some fluctuations
around the theoretical results of $P_{0}$ and $P_{1}$, given by equations (2)
and (3). The results of a simulation in the ibmq_lima QPU are slightly
separated from the theoretical curve, due to imperfect evolution (imperfect
quantum gates, decoherence, dissipation) and measurement errors. Qiskit
measurement error mitigator, which takes into account the errors due to the
state preparation and measurement (SPAM), does not significantly improve the
estimation. State estimation is biased because our model does not faithfully
describe the QPUs. A better model would incorporate corrections to the
expressions for the probabilities $P_{0}$ and $P_{1}.$ In this work, however,
we do not construct an improved model because our main objective is to assess
the performance of the estimator of $s_{z}$ as a function of the parameter
$\theta$, and not of its implementation. Hence, we proceed as if our estimator
were unbiased.

The simulations, performed in the ibmq_lima QPU for three different values of
$\theta$ (see Table 1), qualitatively agree with the expected results,
including the fact that the estimations with $\theta=\pi$ have lower
variances.

Table 1: Estimation ($\hat{s}_{z}$) of the $z$ spin component for the six eigenstates of the Pauli matrices, in qasm_simulator classical simulator and the IBM ibmq_lima QPU, for the angles $\theta=\pi,2\pi/3,\pi/2$. Each state was simulated five times with $5\times 1024$ shots each. Simulator | State | $s_{z}$ |  $\hat{s}_{z}$ ($\theta=\pi$) |  $\hat{s}_{z}$ ($\theta=2\pi/3$) |  $\hat{s}_{z}$ ($\theta=\pi/2$)
---|---|---|---|---|---
QASM | $\ket{\psi_{z0}}$ | $1.0$ | $1.00\pm 0.00$ | $1.00\pm 0.00$ | $1.00\pm 0.00$
| $\ket{\psi_{z1}}$ | $-1.0$ | $-1.00\pm 0.00$ | $-1.03\pm 0.02$ | $-0.99\pm 0.06$
| $\ket{\psi_{x0}}$ | $0.0$ | $0.00\pm 0.04$ | $-0.01\pm 0.04$ | $-0.01\pm 0.05$
| $\ket{\psi_{x1}}$ | $0.0$ | $0.01\pm 0.02$ | $0.00\pm 0.01$ | $0.01\pm 0.06$
| $\ket{\psi_{y0}}$ | $0.0$ | $0.00\pm 0.02$ | $-0.02\pm 0.04$ | $0.03\pm 0.03$
| $\ket{\psi_{y1}}$ | $0.0$ | $0.01\pm 0.04$ | $0.01\pm 0.03$ | $0.01\pm 0.04$
IBMQ | $\ket{\psi_{z0}}$ | $1.0$ | $0.95\pm 0.06$ | $0.96\pm 0.04$ | $0.94\pm 0.06$
| $\ket{\psi_{z1}}$ | $-1.0$ | $-0.92\pm 0.08$ | $-0.8\pm 0.2$ | $-0.8\pm 0.2$
| $\ket{\psi_{x0}}$ | $0.0$ | $-0.04\pm 0.05$ | $0.00\pm 0.01$ | $0.03\pm 0.09$
| $\ket{\psi_{x1}}$ | $0.0$ | $0.1\pm 0.1$ | $0.1\pm 0.1$ | $0.2\pm 0.2$
| $\ket{\psi_{y0}}$ | $0.0$ | $0.00\pm 0.02$ | $-0.03\pm 0.05$ | $-0.01\pm 0.02$
| $\ket{\psi_{y1}}$ | $0.0$ | $0.06\pm 0.07$ | $0.1\pm 0.1$ | $0.1\pm 0.1$

The $\operatorname{qTTF}$ can also be experimentally measured, using a state
2-design, which is set of states X satisfying

$\displaystyle\int
d\psi\braket{\psi}{O_{1}}{\psi}\braket{\psi}{O_{2}}{\psi}=\frac{1}{|X|}\sum_{\psi\in
X}\braket{\psi}{O_{1}}{\psi}\braket{\psi}{O_{2}}{\psi},$ (22)

for all operators $O_{1},\,O_{2}$ [31]. Here, $|X|$ is the number of elements
of $X$. It is easy to see that $F^{-1}$ can be written as
$\braket{\psi}{O_{1}}{\psi}\braket{\psi}{O_{2}}{\psi},$ where
$O_{1}=\frac{1}{2}\left(1+\cos^{2}\frac{\theta}{2}\right)+\frac{1}{2}\sin^{2}\frac{\theta}{2}\sigma_{z}$
and $O_{2}=\frac{1}{2}\csc^{2}\frac{\theta}{2}\left(1-\sigma_{z}\right).$ We
use
$X=\left\\{\ket{\psi_{z_{0}}},\ket{\psi_{z_{1}}},\ket{\psi_{x_{0}}},\ket{\psi_{x_{1}}},\ket{\psi_{y_{0}}},\ket{\psi_{y_{1}}}\right\\}=\left\\{|0\rangle,|1\rangle,\frac{|0\rangle+|1\rangle}{\sqrt{2}},\frac{|0\rangle-|1\rangle}{\sqrt{2}},\frac{|0\rangle+i|1\rangle}{\sqrt{2}},\frac{|0\rangle-i|1\rangle}{\sqrt{2}}\right\\},$
that is, the set of eigenvectors of the three Pauli matrices.

From the simulation results, we can build $N$ samples to make an estimation of
$s_{z}$ as it is described in A.1. In figure 7 we compare the average of
variances of the estimations and the qTTF parameter over $N$ for the six
states of the three Pauli matrices. We can see how the Cramér-Rao bound is
satisfied since
$\overline{\operatorname{Var}(\hat{s}_{z})}\geq\overline{\operatorname{qTTF}}/N.$
As $N\rightarrow\infty$ we see that we are reaching the equality of the
Cramér-Rao bound.

Figure 7:  Graph of the qTTF parameter $\overline{F^{-1}}$ and variance of the
estimations of $N$ samples, for the angles $\theta=\pi,2\pi/3,\pi/2$.

###  3.2 Complete qubit state estimation

For the complete qubit state estimation model, the interaction (11) with the
optimal parameters was simulated using the qasm_simulator and the ibmq_lima
QPU. We consider the six eigenstates of the three Pauli spin matrices,
$\left\\{\ket{\psi_{z_{0}}},\ket{\psi_{z_{1}}},\ket{\psi_{x_{0}}},\ket{\psi_{x_{1}}},\ket{\psi_{y_{0}}},\ket{\psi_{y_{1}}}\right\\}.$
Each state was estimated by averaging the results of five estimations; each
estimation used 1024 shots. The fidelity between the initial state of the
system $\rho_{0}^{S}$ and the estimated state $\hat{\rho}_{0}^{S}$,
$\mathcal{F}(\rho_{0}^{S},\hat{\rho}_{0}^{S})=[\operatorname{Tr}\sqrt{\sqrt{\rho_{0}^{S}}\hat{\rho}_{0}^{S}\sqrt{\rho_{0}^{S}}}]^{2}$,
is low (reaching just 33% for $\ket{\psi_{0}}=\ket{\psi_{x0}}$), as shown in
Table 2.

This error in the estimation of the total spin Bloch vector is a consequence
of the number of gates needed to simulate this system (around $44$ CNOT gates
and $82$ single-qubit gates). Another quantity that helps to understand the
low fidelity of the estimated states is the circuit depth, the largest number
of gates, from the input to the output, always moving ahead along qubit wires.
The transpiled circuit that represents the evolution operator, has a circuit
depth of $83.$

Table 2:  Estimation $(\hat{s}_{x},\hat{s}_{y},\hat{s}_{z})$ of the Bloch vector $\vec{s}=(s_{x},s_{y},s_{z})$ corresponding to the eigenstates of the three Pauli operators, using the optimal values of the circuit parameters $\theta_{A}$ and $\theta_{B}$ approximately $3.45$ and $-8.42$ respectively. Estimations were executed in qasm_simulator classical simulator and the IBM ibmq_lima QPU, using $5\times 1024$ shots.  Simulator | State | $s_{x}$ | $\hat{s}_{x}$ | $s_{y}$ | $\hat{s}_{y}$ | $s_{z}$ | $\hat{s}_{z}$ | Fidelity
---|---|---|---|---|---|---|---|---
QASM | $\ket{\psi_{z0}}$ | 0.0 | $-0.01\pm 0.06$ | 0.0 | $-0.03\pm 0.08$ | 1.0 | $0.94\pm 0.04$ | $99.7\%$
| $\ket{\psi_{z1}}$ | 0.0 | $-0.02\pm 0.07$ | 0.0 | $-0.02\pm 0.03$ | -1.0 | $-0.97\pm 0.03$ | $99.8\%$
| $\ket{\psi_{x0}}$ | 1.0 | $0.99\pm 0.02$ | 0.0 | $-0.01\pm 0.04$ | 0.0 | $-0.00\pm 0.06$ | $99.8\%$
| $\ket{\psi_{x1}}$ | -1.0 | $-0.992\pm 0.007$ | 0.0 | $0.04\pm 0.01$ | 0.0 | $0.02\pm 0.08$ | $99.8\%$
| $\ket{\psi_{y0}}$ | 0.0 | $-0.01\pm 0.04$ | 1.0 | $0.94\pm 0.06$ | 0.0 | $0.05\pm 0.04$ | $99.8\%$
| $\ket{\psi_{y1}}$ | 0.0 | $0.01\pm 0.08$ | -1.0 | $-0.994\pm 0.004$ | 0.0 | $0.05\pm 0.06$ | $99.7\%$
IBMQ | $\ket{\psi_{z0}}$ | 0.0 | $-0.60\pm 0.08$ | 0.0 | $0.2\pm 0.1$ | 1.0 | $0.75\pm 0.07$ | $85\%$
| $\ket{\psi_{z1}}$ | 0.0 | $-0.7\pm 0.2$ | 0.0 | $0.34\pm 0.07$ | -1.0 | $0.3\pm 0.3$ | $37\%$
| $\ket{\psi_{x0}}$ | 1.0 | $-0.3\pm 0.3$ | 0.0 | $0.2\pm 0.1$ | 0.0 | $0.7\pm 0.2$ | $32\%$
| $\ket{\psi_{x1}}$ | -1.0 | $-0.75\pm 0.07$ | 0.0 | $0.1\pm 0.1$ | 0.0 | $0.5\pm 0.2$ | $92\%$
| $\ket{\psi_{y0}}$ | 0.0 | $-0.75\pm 0.04$ | 1.0 | $0.17\pm 0.08$ | 0.0 | $0.63\pm 0.06$ | $55\%$
| $\ket{\psi_{y1}}$ | 0.0 | $-0.5\pm 0.2$ | -1.0 | $0.2\pm 0.1$ | 0.0 | $0.74\pm 0.08$ | $41\%$

Even if this is already a low probability for having a good simulation, we
have not taken into account other error sources, such as state preparation and
measurements errors (SPAM), as well as those related to coherence time [32,
33]; which justify even more the low fidelities obtained in Table 2.

In figure 8 we take $N$ samples of spin estimations for each of the six
eigenstates of the Pauli matrices and compare the mean of the qTTF parameter
over $N$ and the average of the covariance matrix trace. The Cramér-Rao bound
$\overline{\operatorname{Cov}(s,s^{\prime})}\geq\overline{\operatorname{qTTF}}/N$
is not only always satisfied, but it also provides a good estimation for the
variance.

Figure 8: Mean of the qTTF error parameter over $N$ and trace of the
covariance matrix means over the six eigenstates of the Pauli matrix against
the $N$ samples.

###  3.3 Optimal complete qubit state estimation

Low fidelity results found before are not a characteristic of the tomographic
method employed, but from the implementation in noisy quantum computers. We
modify the $U_{AB}$ operator in equation (12), minimizing the number of CNOTs
gates used. The resulting circuit is shown in fig. 9.

$\Qcircuit@C=1.0em@R=0.2em@!R{\\\
\nghost{{A}:}&\lstick{{A}:}\gate{{H}}\barrier[0em]{2}\qw\gate{{U}_{\vec{\theta}_{A_{1}}}}\targ\gate{{U}_{\vec{\theta}_{A_{2}}}}\qw\qw\barrier[0em]{2}\qw\gate{{H}}\barrier[0em]{2}\qw\meter\qw\qw\qw\\\
\nghost{{S}:}\lstick{{S}:}\gate{{U_{\xi}}}\qw\qw\ctrl{-1}\gate{{H}}\ctrl{1}\gate{{H}}\qw\qw\qw\qw\qw\qw\qw\\\
\nghost{{B}:}\lstick{{B}:}\gate{{H}}\qw\qw\qw\gate{{U}_{\vec{\theta}_{B_{1}}}}\targ\gate{{U}_{\vec{\theta}_{B_{2}}}}\qw\gate{{H}}\qw\qw\meter\qw\qw\\\
\nghost{c:}\lstick{c:}\lstick{}\cw\cw\cw\cw\cw\cw\cw\cw\cw\cw\dstick{}\cw\cwx[-3]\dstick{}\cw\cwx[-1]\cw\cw\\\
\\\ }$ Figure 9: Quantum circuit proposed for optimal complete qubit state
estimation in NISQ quantum computers.

The circuit in figure 9 depends on twelve parameters, three for each $U$ gate,
$U(\vec{\theta})=U(\theta,\phi,\lambda)=\begin{pmatrix}\cos\left(\theta\right)&-e^{i\lambda}\sin\left(\theta\right)\\\
e^{i\phi}\sin\left(\theta\right)&e^{i(\phi+\lambda)}\cos\left(\theta\right)\end{pmatrix}$
. Optimization of qTTF over these parameters is made using random initial
parameters. The best parameters

$\displaystyle\vec{\theta}_{A_{1}}=(0.59,1.58,2.52),\quad\vec{\theta}_{A_{2}}=(2.55,1.94,0.31),$
$\displaystyle\vec{\theta}_{B_{1}}=(0.70,4.31,3.46),\quad\vec{\theta}_{B_{2}}=(0.67,6.47,4.47),$

correspond to an average error $\bar{\Delta}\approx 8.$

The estimation of the quantum state is good enough by linear inversion, which
relies on linear regression estimation algorithms [34]. Thus, it has a lower
computational cost with respect to the model in subsection 3.2 which performs
tomography using a maximum likelihood estimation algorithm. The results of
state estimation for the six eigenstates of the Pauli operators to perform the
estimation protocol can be found in Table 3. The fidelities of the estimated
states are above $98\%$ when the model of figure 9 is implemented in an IBM
QPU.

Table 3:  Estimation $(\hat{s}_{x},\hat{s}_{y},\hat{s}_{z})$ of the Bloch vector $\vec{s}=(s_{x},s_{y},s_{z})$ corresponding to the eigenstates of the three Pauli operators, using the optimal circuit described in figure 9. Estimations were executed in qasm_simulator classical simulator and the IBM ibmq_lima QPU, using $5\times 1024$ shots.  Simulator | State | $s_{x}$ | $\hat{s}_{x}$ | $s_{y}$ | $\hat{s}_{y}$ | $s_{z}$ | $\hat{s}_{z}$ | Fidelity
---|---|---|---|---|---|---|---|---
QASM | $\ket{\psi_{z0}}$ | 0.0 | $0.00\pm 0.07$ | 0.0 | $-0.02\pm 0.02$ | 1.0 | $1.00\pm 0.06$ | $99.8\%$
| $\ket{\psi_{z1}}$ | 0.0 | $0.01\pm 0.05$ | 0.0 | $0.03\pm 0.05$ | -1.0 | $-0.99\pm 0.05$ | $99.8\%$
| $\ket{\psi_{x0}}$ | 1.0 | $0.95\pm 0.04$ | 0.0 | $0.02\pm 0.01$ | 0.0 | $0.03\pm 0.01$ | $99.9\%$
| $\ket{\psi_{x1}}$ | -1.0 | $-1.00\pm 0.02$ | 0.0 | $-0.02\pm 0.05$ | 0.0 | $0.00\pm 0.05$ | $99.9\%$
| $\ket{\psi_{y0}}$ | 0.0 | $-0.04\pm 0.03$ | 1.0 | $1.02\pm 0.03$ | 0.0 | $-0.05\pm 0.07$ | $99.7\%$
| $\ket{\psi_{y1}}$ | 0.0 | $0.00\pm 0.05$ | -1.0 | $-1.01\pm 0.07$ | 0.0 | $0.02\pm 0.05$ | $99.9\%$
IBMQ | $\ket{\psi_{z0}}$ | 0.0 | $0.12\pm 0.04$ | 0.0 | $0.05\pm 0.03$ | 1.0 | $0.96\pm 0.03$ | $99.5\%$
| $\ket{\psi_{z1}}$ | 0.0 | $-0.14\pm 0.03$ | 0.0 | $0.01\pm 0.03$ | -1.0 | $-0.81\pm 0.04$ | $99.2\%$
| $\ket{\psi_{x0}}$ | 1.0 | $0.93\pm 0.04$ | 0.0 | $0.14\pm 0.03$ | 0.0 | $0.04\pm 0.04$ | $99.2\%$
| $\ket{\psi_{x1}}$ | -1.0 | $-0.88\pm 0.04$ | 0.0 | $-0.09\pm 0.04$ | 0.0 | $0.15\pm 0.04$ | $98.9\%$
| $\ket{\psi_{y0}}$ | 0.0 | $-0.00\pm 0.06$ | 1.0 | $0.87\pm 0.03$ | 0.0 | $0.14\pm 0.04$ | $99.1\%$
| $\ket{\psi_{y1}}$ | 0.0 | $0.10\pm 0.05$ | -1.0 | $-0.75\pm 0.03$ | 0.0 | $0.02\pm 0.09$ | $99.1\%$

##  4  Conclusions

We reviewed qubit tomographic models that make use of one or two parameters,
enabling the estimation of a single component or the entire Bloch vector
associated with the initial state of the qubit, respectively. We use
tomography techniques that rely on either linear regression estimators or the
discrete maximum-likelihood estimator. These models are based on the
interaction between the qubit under study and one or two auxiliary qubits that
act as meters. We used the quantum tomographic transfer function (qTTF), the
trace of the inverse of the Fisher information matrix averaged over all pure
states, as a measure of the global error of the models. The qTTF, which
displays enormous variations as a function of the parameters of the models, is
significantly higher for the models which aim to completely estimate the
state. The best complete estimation model has a global error about 25 times
larger than its single-component counterpart.

For the case of one component estimation, the smallest global error is
obtained for the value of the parameter which maximizes the entangling power
of the associated unitary transformation. This suggests that entanglement
between the system and the meter might play a role in the quality of the
estimation. The estimation error for complete state estimation model, before
averaging over the possible initial states of the system, presents moderate
fluctuations as a function of the initial state.

We implemented both estimation models on an IBM QPU. Experimental results for
several values of the parameter of the single-component model agrees well with
the theoretical results, with small deviations due to imperfect evolution and
measurement errors. The experimental errors corresponding to the best complete
state estimation model, due to the number of CNOT and single-qubit gates of
the implementation, are so large that the fidelity between the estimated and
the prepared initial state of the qubit can be as low as $30\%$.

The fidelities obtained with a modified model, taking into account the natural
interactions of the QPU, were above $98\%$. Additionally, the corresponding
qTTF is half the error for the original complete state estimation optimal
model.

The methodology employed in this work can be applied to other estimation
models (states in larger Hilbert spaces, different system-meter interactions,
unitary or dissipative evolution), although the optimization becomes
increasingly difficult as the number of parameters (state and evolution)
increases. These modifications of the estimation model can be performed taking
into account the characteristics of the QPU, for example, its basis gates or
topology. This makes the model practical for being adapted to any quantum
device or any other research purpose, such as simulations of experimental
setups.

Our estimation method can be extended to $n$ qubits employing a strategy
similar to the one reported in [35], in which each qubit interacts with only
one auxiliary qubit. In our case, two auxiliary qubits and two CNOTs are
necessary for each qubit, exactly as shown in Figure 9. In other words, state
estimation of $n$ qubits, $2n$ auxiliary qubits and $2n$ CNOT gates are
required. Notice that the total number of outcomes corresponding to the
measurement of $\sigma_{z}$ for each auxiliary qubit is $2^{2n},$ exactly the
number of parameters of a general unnormalized density matrix for $n$ qubits.
Notice that single-observable tomography for $n$ qubits requires $n$ auxiliary
qubits and $3n$ CNOT gates. In effect, for each qubit an auxiliary qubit is
required [35] as mentioned before. An optimal circuit for general two qubit
interaction requires at most 3 CNOT gates and 15 elementary one-qubit gates
[36]. In practice, the topology of the QPUs pose additional constraints, which
might considerably increase the number of required CNOT gates.

As a final comment, we would like to highlight the need for a general theory
of resource management for state estimation, addressing how different
resources can be optimally substituted for one another. For example, our
results suggest that $n^{2}$ local measurements are equivalent to a single
local observable, $n$ auxiliary qubits, $15n$ elementary qubit gates and $3n$
CNOT gates.

##  5 Acknowledgements

This research received no specific grant from any funding agency.

##  6 Conflict of interest

The authors have no conflicts of interest to declare.

##  7 Data Availability

The data as well as the codes that support the findings of this study are
available from the corresponding authors on reasonable request.

##  Appendix A Equivalence between the variance of the estimator and the
Fisher Error parameter.

###  A.1 Estimation of $s_{z}$

From the experiment, we have a set of observations $S=\\{s_{k}\\}_{k=1}^{N}.$
Here $s_{k}$ takes the values $s_{0}=s_{z}(1,0,\theta)$ with probability
$P_{0}=p$ and $s_{1}=s_{z}(0,1,\theta)$ with probability $P_{1}=1-p$, where
$s_{z}$ is given by eq. (4),
$s_{z}(P_{0},P_{1},\theta)=\csc^{2}\left(\frac{\theta}{2}\right)\left(P_{0}-P_{1}-\cos^{2}\left(\frac{\theta}{2}\right)\right).$

Each random variable $s_{k}$ has a Bernoulli probability distribution, while
the set $S$ has a binomial probability distribution.

The mean $\bar{s}$ and variance $\sigma^{2}$ of the observations are given by

$\displaystyle\bar{s}$
$\displaystyle=\frac{1}{N}\sum_{k=1}^{N}s_{k}=\frac{n_{0}}{N}s_{0}+\frac{n_{1}}{N}s_{1},\quad\textrm{and}$
$\displaystyle\sigma^{2}$
$\displaystyle=\frac{1}{N-1}\sum_{k=1}^{N}(s_{k}-\bar{s})^{2}=\frac{n_{0}(s_{0}-\bar{s})^{2}}{N-1}+\frac{n_{1}(s_{1}-\bar{s})^{2}}{N-1},$

where $n_{0}$ is the number of $s_{0}$ observations and $n_{1}$ is the number
of $s_{1}$ observations.

The variance can be written as

$\sigma^{2}=\frac{N}{N-1}f_{0}f_{1}(s_{0}-s_{1})^{2}$

where $f_{0}=\frac{n_{0}}{N}$ and $f_{1}=\frac{n_{1}}{N}$. When
$N\rightarrow\infty$, $f_{0}\rightarrow P_{0}=p$ and $f_{1}\rightarrow
P_{1}=1-p$. Therefore, $\bar{s}\rightarrow s_{0}p+s_{1}(1-p)$ and
$\sigma^{2}\rightarrow p(1-p)(s_{0}-s_{1})^{2},$ which are the theoretical
values of the mean and the variance of a Bernoulli distribution respectively.

When $\bar{s}$ is taken as a stochastic variable, its variance is

$\sigma^{2}(\bar{s})=\frac{1}{N^{2}}\sum_{k=1}^{N}\sigma=\frac{1}{N-1}f_{0}f_{1}(s_{0}-s_{1})^{2}=\frac{1}{N-1}F^{-1},$

where we have used $\operatorname{Var}(\alpha
X)=\alpha^{2}\operatorname{Var}(X)$ with $\alpha$ constant, and
$\operatorname{Var}(\sum_{k}s_{k})=\sum_{k=1}^{N}\operatorname{Var}(s_{k})=N\sigma.$
Here, $F^{-1}$ is the inverse of the (scaled) Fisher parameter, given by eq.
(6).

###  A.2 Total spin estimation

Similarly, as we have done for the estimation of $s_{z}$, here we have a set
of observations $S^{i}=\\{s_{k}^{i}\\}_{k=1}^{N}$ for the $i$-th component of
the spin. Those observations come from the measurement of the meters $A$ and
$B.$ Then, we can define the mean value for each direction of the spin as

$\bar{{s}^{i}}=\frac{1}{N}\sum_{k=1}^{N}s_{k}^{i}=\frac{n_{00}}{N}s_{00}^{i}+\frac{n_{01}}{N}s_{01}^{i}+\frac{n_{10}}{N}s_{10}^{i}+\frac{n_{11}}{N}s_{11}^{i},$

where $i=x,y,z.$ Here we have to take into account that
$s^{i}=s^{i}(n_{00},n_{01},n_{10},n_{11}).$ Then, $s_{00}^{i}=s^{i}(N,0,0,0),$
$s_{01}^{i}=s^{i}(0,N,0,0),$ $s_{10}^{i}=s^{i}(0,0,N,0),$ and
$s_{11}^{i}=s^{i}(0,0,0,N).$ For convenience, we change the notation
$n_{00}\rightarrow n_{1}$, $n_{01}\rightarrow n_{2}$, $n_{10}\rightarrow
n_{3}$ y $n_{11}\rightarrow n_{4}.$ The variances of these random variables,
in each direction, are

$\operatorname{Var}(s^{i}_{k})={\sigma^{i}}^{2}=\frac{1}{N-1}\sum_{k=1}^{N}(s^{i}_{k}-\bar{{s}^{i}})^{2}.$

This can be written, taking into account the explicit form of the mean, as

${\sigma^{i}}^{2}=\frac{N}{N-1}\left(\sum_{l=1}^{4}f_{l}(s_{l}^{i}-\sum_{m=1}^{4}s_{m}^{i}f_{m})^{2}\right).$

Indices $n,l$ correspond to the $4$ possible results of the meter’s
measurements. Similarly to the previous case, we have $f_{l}=\frac{n_{l}}{N}.$
When $N\rightarrow\infty,$ we have

${\sigma^{i}}^{2}\rightarrow\sum_{l=1}^{4}p_{l}(s_{l}^{i}-\sum_{m=1}^{4}p_{m}s_{m}^{i})^{2}.$

In this limit, the mean value of the estimation in the $i$-th direction can be
written as

$\bar{{s}^{i}}\rightarrow
p_{1}s_{1}^{i}+p_{2}s_{2}^{i}+p_{3}s_{3}^{i}+p_{4}s_{4}^{i}.$

This relation can be written in the matrix form
$\bar{\vec{{s}}}=\mathbb{S}\vec{p},$ where

$\bar{\vec{{s}}}=\begin{pmatrix}1\\\ \bar{{s}^{x}}\\\ \bar{{s}^{y}}\\\
\bar{{s}^{z}}\end{pmatrix},\hskip
10.00002pt\mathbb{S}=\begin{pmatrix}1&1&1&1\\\
{{s}_{1}^{x}}&{{s}_{2}^{x}}&{{s}_{3}^{x}}&{{s}_{4}^{x}}&\\\
{{s}_{1}^{y}}&{{s}_{2}^{y}}&{{s}_{3}^{y}}&{{s}_{4}^{y}}&\\\
{{s}_{1}^{z}}&{{s}_{2}^{z}}&{{s}_{3}^{z}}&{{s}_{4}^{z}}&\end{pmatrix},\hskip
10.00002pt\vec{p}=\begin{pmatrix}p_{1}\\\ p_{2}\\\ p_{3}\\\
p_{4}\end{pmatrix}.$

Here we have to take into account that $\sum_{k=1}^{4}p_{k}=1.$

Equation (18), $\vec{s}=\mathbb{T}^{-1}\vec{p},$ allows us to build an
estimator $\hat{\vec{s}},$

$\hat{\vec{s}}=\mathbb{T}^{-1}\mathbb{S}^{-1}\bar{\vec{{s}}}=\mathbb{G}\bar{\vec{{s}}}.$

The $i$-th component of its variance
$\operatorname{Var}(\hat{\vec{s}})=\operatorname{Var}\left(\mathbb{G}\bar{\vec{{s}}}\right)$
is

$\displaystyle\operatorname{Var}\left(\hat{{s}}^{i}\right)$
$\displaystyle=\operatorname{Var}\left(\sum_{j=0}^{4}G_{ij}\bar{s}_{j}\right)=\sum_{j=0}^{4}G_{ij}^{2}\operatorname{Var}\left(\frac{1}{N}\sum_{k=1}^{N}s^{j}_{k}\right)$
$\displaystyle=\sum_{j=0}^{4}G_{ij}^{2}\frac{1}{N^{2}}\sum_{k=1}^{N}{\sigma^{j}}^{2}.$

In terms of the experimental results we have

$\operatorname{Var}\left(\hat{{s}}^{i}\right)=\frac{1}{N}\sum_{j=0}^{4}G_{ij}^{2}{\sigma^{j}}^{2}=\frac{1}{N}\left(\sum_{j=0}^{4}G_{ij}^{2}\sum_{l=1}^{4}p_{l}(s_{l}^{j}-\sum_{m=1}^{4}p_{m}s_{m}^{j})^{2}.\right),$

where we have used $\operatorname{Var}(\alpha
X)=\alpha^{2}\operatorname{Var}(X)$ with $\alpha$ constant, and
$\operatorname{Var}(\sum_{k}s_{k}^{i})=\sum_{k=1}^{N}\operatorname{Var}(s_{k}^{i})=N\sigma^{i}.$
Here $i,j=0,x,y,z.$ i.e. they label the components of the spin vector. The
component $0$ satisfy the normalization $s^{0}=1,$ then, it is clear that
$\operatorname{Var}\left(\hat{{s}}^{0}\right)=0.$ On the other hand, indices
$l,m$ correspond to the four possible results of the meter’s measurements,
$00,01,10$ and $11.$
It can be numerically shown that (see the supplementary material)

$\sum_{j=0}^{4}G_{ij}^{2}\sum_{l=1}^{4}p_{l}(s_{l}^{j}-\sum_{m=1}^{4}p_{m}s_{m}^{j})^{2}=\left(\mathbb{F}^{-1}\right)_{ii},$

which corresponds to the Fisher error (21).

##  Appendix B Explicit expressions for coefficients $a_{\mu}$, $b_{\mu},$ and
$c_{\mu}$

Setting $\theta_{C}=\sqrt{\theta_{A}^{2}+\theta_{B}^{2}},$ we can write the
coefficients $a_{\mu}$, $b_{\mu},$ and $c_{\mu}$ as follows.

$\displaystyle
a_{0}=\frac{1}{8}\cos\left(\frac{\theta_{A}}{2}\right)\left[\cos\left(\frac{\theta_{A}}{2}\right)+\frac{\theta_{B}\sin\left(\frac{\theta_{B}}{2}\right)\sin\left(\frac{\theta_{C}}{2}\right)}{\theta_{C}}+\cos\left(\frac{\theta_{B}}{2}\right)\cos\left(\frac{\theta_{C}}{2}\right)\right],$
$\displaystyle
a_{1}=\frac{1}{8}\sin\left(\frac{\theta_{A}}{2}\right)\left[\sin\left(\frac{\theta_{B}}{2}\right)\cos\left(\frac{\theta_{C}}{2}\right)-\frac{\theta_{B}\cos\left(\frac{\theta_{B}}{2}\right)\sin\left(\frac{\theta_{C}}{2}\right)}{\theta_{C}}\right],$
$\displaystyle
a_{2}=-\frac{\theta_{A}\sin\left(\frac{\theta_{A}}{2}\right)\sin\left(\frac{\theta_{B}}{2}\right)\sin\left(\frac{\theta_{C}}{2}\right)}{8\theta_{C}},$
$\displaystyle
a_{3}=-\frac{1}{8}\sin\left(\frac{\theta_{A}}{2}\right)\left[\frac{\theta_{A}\cos\left(\frac{\theta_{B}}{2}\right)\sin\left(\frac{\theta_{C}}{2}\right)}{\theta_{C}}+\sin\left(\frac{\theta_{A}}{2}\right)\right],$
$\displaystyle
c_{0}=\frac{1}{32}\left[4\cos\left(\frac{\theta_{C}}{2}\right)\cos\left(\frac{\theta_{A}}{2}+\frac{\theta_{B}}{2}\right)+\cos\left(\theta_{A}-\theta_{B}\right)+\cos\left(\theta_{A}\right)+\cos\left(\theta_{B}\right)+1\right],$
$\displaystyle
c_{1}=\frac{1}{8}\left[\cos\left(\frac{\theta_{A}}{2}\right)\sin\left(\frac{\theta_{B}}{2}\right)\sin\left(\frac{\theta_{A}-\theta_{B}}{2}\right)-\frac{\theta_{B}\sin\left(\frac{\theta_{C}}{2}\right)\sin\left(\frac{\theta_{A}+\theta_{B}}{2}\right)}{\theta_{C}}\right],$
$\displaystyle
c_{2}=-\frac{1}{8}\sin\left(\frac{\theta_{A}}{2}\right)\sin\left(\frac{\theta_{B}}{2}\right)\sin\left(\frac{\theta_{A}-\theta_{B}}{2}\right).$

The remaining coefficients can be written in terms of the previous ones as
$b_{0}(\theta_{A},\theta_{B})=a_{0}\left(\theta_{B},\theta_{A}\right)$,
$b_{1}\left(\theta_{A},\theta_{B}\right)=a_{3}\left(\theta_{B},\theta_{A}\right),$
$b_{2}\left(\theta_{A},\theta_{B}\right)=-a_{2}\left(\theta_{B},\theta_{A}\right),$
$b_{3}\left(\theta_{A},\theta_{B}\right)=a_{1}\left(\theta_{B},\theta_{A}\right),$
$c_{3}\left(\theta_{A},\theta_{B}\right)=c_{1}\left(\theta_{B},\theta_{A}\right).$

## References

## References

  * [1] Feynman RP.  Simulating physics with computers.  International Journal of Theoretical Physics. 1982;21.
  * [2] Feynman RP.  Quantum mechanical computers.  Foundations of Physics. 1986;16.
  * [3] Shor PW.  Algorithms for quantum computation: discrete logarithms and factoring.  In: Proceedings 35th Annual Symposium on Foundations of Computer Science; 1994. p. 124-34.
  * [4] Grover LK.  A Fast Quantum Mechanical Algorithm for Database Search.  In: Proceedings of the Twenty-Eighth Annual ACM Symposium on Theory of Computing. STOC ’96. New York, NY, USA: Association for Computing Machinery; 1996. p. 212-9.  Available from: https://doi.org/10.1145/237814.237866.
  * [5] Kiktenko EO, Kublikova DN, Fedorov AK.  Estimating the precision for quantum process tomography.  Optical Engineering. 2020;59(6):1 9.  Available from: https://doi.org/10.1117/1.OE.59.6.061614.
  * [6] Gaikwad A, Shende K, Arvind, Dorai K.  Implementing efficient selective quantum process tomography of superconducting quantum gates on IBM quantum experience.  Scientific Reports. 2022 Mar;12(1):3688.  Available from: https://doi.org/10.1038/s41598-022-07721-3.
  * [7] Banaszek K, Cramer M, Gross D.  Focus on quantum tomography.  New Journal of Physics. 2013 dec;15(12):125020.  Available from: https://doi.org/10.1088/1367-2630/15/12/125020.
  * [8] James DFV, Kwiat PG, Munro WJ, White AG.  Measurement of qubits.  Phys Rev A. 2001 Oct;64:052312.  Available from: https://link.aps.org/doi/10.1103/PhysRevA.64.052312.
  * [9] Bechmann-Pasquinucci H, Tittel W.  Quantum cryptography using larger alphabets.  Phys Rev A. 2000 May;61:062308.  Available from: https://link.aps.org/doi/10.1103/PhysRevA.61.062308.
  * [10] Dong D, Petersen IR.  Quantum estimation, control and learning: opportunities and challenges.  Annual Reviews in Control. 2022.
  * [11] Gupta R, Xia R, Levine RD, Kais S.  Maximal Entropy Approach for Quantum State Tomography.  PRX Quantum. 2021 Feb;2:010318.  Available from: https://link.aps.org/doi/10.1103/PRXQuantum.2.010318.
  * [12] Gupta R, Levine RD, Kais S.  Convergence of a Reconstructed Density Matrix to a Pure State Using the Maximal Entropy Approach.  The Journal of Physical Chemistry A. 2021;125(34):7588-95.  PMID: 34410718.  Available from: https://doi.org/10.1021/acs.jpca.1c05884.
  * [13] Lohani S, Searles TA, Kirby BT, Glasser RT.  On the Experimental Feasibility of Quantum State Reconstruction via Machine Learning.  IEEE Transactions on Quantum Engineering. 2021;2:1-10.
  * [14] Lukens JM, Law KJH, Jasra A, Lougovski P.  A practical and efficient approach for Bayesian quantum state estimation.  New Journal of Physics. 2020 jun;22(6):063038.  Available from: https://dx.doi.org/10.1088/1367-2630/ab8efa.
  * [15] Qi B, Hou Z, Wang Y, Dong D, Zhong HS, Li L, et al.  Adaptive quantum state tomography via linear regression estimation: Theory and two-qubit experiment.  npj Quantum Information. 2017 Apr;3(1):19.  Available from: https://doi.org/10.1038/s41534-017-0016-4.
  * [16] Jamiołkowski A.  Minimal number of operators for observability ofN-level quantum systems.  International Journal of Theoretical Physics. 1983;22(4):369-76.
  * [17] Czerwiński A.  Optimal evolution models for quantum tomography.  Journal of Physics A: Mathematical and Theoretical. 2016;49(7):075301.
  * [18] Czerwinski A.  Selected Concepts of Quantum State Tomography.  Optics. 2022;3(3):268-86.
  * [19] D’Ariano GM.  Universal quantum observables.  Physics Letters A. 2002;300(1):1-6.
  * [20] Peres A.  When is a quantum measurement?  American Journal of Physics. 1986;54(8):688-92.  Available from: https://doi.org/10.1119/1.14505.
  * [21] Saavedra D, Fonseca-Romero KM.  Complete and incomplete state estimation via the simultaneous unsharp measurement of two incompatible qubit operators.  Phys Rev A. 2019 Apr;99:042130.  Available from: https://link.aps.org/doi/10.1103/PhysRevA.99.042130.
  * [22] Řeháček J, Teo YS, Hradil Z.  Determining which quantum measurement performs better for state estimation.  Phys Rev A. 2015 Jul;92:012108.  Available from: https://link.aps.org/doi/10.1103/PhysRevA.92.012108.
  * [23] Pereira L, Zambrano L, Delgado A.  Scalable estimation of pure multi-qubit states.  npj Quantum Information. 2022 May;8(1):57.  Available from: https://doi.org/10.1038/s41534-022-00565-9.
  * [24] Kay SM.  Fundamentals of Statistical Signal Processing: Estimation Theory.  Prentice Hall PTR; 1993.
  * [25] Benenti G, Casati G, Strini G.  Principles of Quantum Computation and Information. vol. 2.  World Scientific; 2004.  Available from: https://www.worldscientific.com/doi/abs/10.1142/5528.
  * [26] Zanardi P, Zalka C, Faoro L.  Entangling power of quantum evolutions.  Phys Rev A. 2000 Aug;62:030301.  Available from: https://link.aps.org/doi/10.1103/PhysRevA.62.030301.
  * [27] Paris M, Řeháček J.  Quantum State Estimation.  Springer, Berlin, Heidelberg; 2004.  Available from: https://link.springer.com/book/10.1007/b98673.
  * [28] Teo YS.  Introduction to Quantum-State Estimation.  World Scientific; 2015.  Available from: https://www.worldscientific.com/doi/abs/10.1142/9617.
  * [29] Řeháček J, Hradil Z, Knill E, Lvovsky AI.  Diluted maximum-likelihood algorithm for quantum tomography.  Phys Rev A. 2007 Apr;75:042108.  Available from: https://link.aps.org/doi/10.1103/PhysRevA.75.042108.
  * [30] Virtanen P, Gommers R, Oliphant TE, Haberland M, Reddy T, Cournapeau D, et al.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  Nature Methods. 2020;17:261-72.
  * [31] Bendersky A, Pastawski F, Paz JP.  Selective and efficient quantum process tomography.  Phys Rev A. 2009 Sep;80:032116.  Available from: https://link.aps.org/doi/10.1103/PhysRevA.80.032116.
  * [32] Jattana MS, Jin F, De Raedt H, Michielsen K.  General error mitigation for quantum circuits.  Quantum Information Processing. 2020 Nov;19(11):414.  Available from: https://doi.org/10.1007/s11128-020-02913-0.
  * [33] Unruh WG.  Maintaining coherence in quantum computers.  Phys Rev A. 1995 Feb;51:992-7.  Available from: https://link.aps.org/doi/10.1103/PhysRevA.51.992.
  * [34] Qi B, Hou Z, Li L, Dong D, Xiang G, Guo G.  Quantum State Tomography via Linear Regression Estimation.  Scientific Reports. 2013 Dec;3(1):3496.  Available from: https://doi.org/10.1038/srep03496.
  * [35] Wang H, Zheng W, Yu Y, Jiang M, Peng X, Du J.  Determining an n-qubit state by a single apparatus through a pairwise interaction.  Physical Review A. 2014;89(3):032103.
  * [36] Vatan F, Williams C.  Optimal quantum circuits for general two-qubit gates.  Physical Review A. 2004;69(3):032315.


Xianzhao Chen, Yist Y. Lin, Kang Wang, Yi He, Zejun Ma

# Improving Frame-level Classifier for Word Timings with Non-peaky CTC in End-
to-End Automatic Speech Recognition

###### Abstract

End-to-end (E2E) systems have shown comparable performance to hybrid systems
for automatic speech recognition (ASR). Word timings, as a by-product of ASR,
are essential in many applications, especially for subtitling and computer-
aided pronunciation training. In this paper, we improve the frame-level
classifier for word timings in E2E system by introducing label priors in
connectionist temporal classification (CTC) loss, which is adopted from prior
works, and combining low-level Mel-scale filter banks with high-level ASR
encoder output as input feature. On the internal Chinese corpus, the proposed
method achieves 95.68%/94.18% compared to the hybrid system 93.0%/90.22% on
the word timing accuracy metrics. It also surpass a previous E2E approach with
an absolute increase of 4.80%/8.02% on the metrics on 7 languages. In
addition, we further improve word timing accuracy by delaying CTC peaks with
frame-wise knowledge distillation, though only experimenting on LibriSpeech.

Index Terms: automatic speech recognition, connectionist temporal
classification, word timing, end-to-end

##  1 Introduction

Word timing is the start and end time of a word. It is widely used in many ASR
applications, such as computer-aided pronunciation training, subtitling and
filler word removal. Word timings can also be converted into frame-level
classifications for pre-training acoustic model in ASR systems [1].
Correspondingly, the word timing metric can be used to measure the performance
of the frame-level classifier [2].

Recently, E2E systems adopting Deep Neural Networks (DNN) have demonstrated
excellent ASR performance [3, 4]. Unlike traditional hybrid systems, E2E
systems directly align inputs and outputs without modeling phonetic alignment.
In E2E systems, Listen, Attend and Spell (LAS) [5] uses attention mechanisms,
CTC [6] and RNN Transducer (RNN-T) [7] introduce blank token to align
variable-length inputs and outputs. Though E2E systems trained with soft
alignments achieve superior ASR performance compared to traditional hybrid
systems, the estimated word timing in such systems is usually not comparable
with that in traditional ones.

There are several methods for estimating word timings in hybrid systems and
E2E systems. In a hybrid system, such as DNN-Hidden Markov Models (HMM)
system, a DNN model is trained with frame-level classification using
alignments from a pretrained GMM-HMM system, and the HMMs are trained using
the posterior probabilities of acoustic features provided by DNN [8]. Accurate
word timings are obtained from forced alignment of the recognition results.
However, a hybrid system requires many steps of training, making it less
desirable.

In E2E systems, word timings can be estimated by the forced alignment results
of character-level CTC models, where the CTC peak of the first character
indicate the word start time and the CTC peak of the last character indicate
the word end time [9]. The CTC model cannot estimate word timings well when
the duration of the modeling unit is relatively long, e.g., Chinese
characters. Because the blank probability of CTC model is dominant in almost
all frames, and the non-blank probability is only relatively high in few
frames. This is called the peaky behavior [10]. CTC-based alignments for word
timings can be improved by alleviating the peaky behavior [11, 12], but these
methods have complicated regularization terms which require modifications of
the model architecture or CTC loss topology respectively. A two-pass model
trained with frame-level classification, where alignments are from traditional
hybrid systems, predicts word timings by constrained attention weights [2].
The method achieves state-of-the-art performance on word timing metrics, but
its dependence on hybrid systems makes it less favored. A purely E2E system
was proposed where alignments from a hybrid system are not needed [13]. But it
requires training a frame-level classifier with CTC loss, manually expanding
CTC peaks to word timings, and training another classifier using the word
timings with cross-entropy (CE) loss. Continuous integrate-and-fire (CIF) is a
soft and monotonic alignment well suited for estimating word timings in E2E
non-autoregressive ASR systems. The CIF-based system obtains good acoustic
boundaries by accumulating weights of frames, until exceeding a threshold
value, and firing tokens. With several steps of post-processing, its
performance is comparable to HMM-based system and only need to be verified in
languages other than Chinese [14].

To the best of our knowledge, we are the first to adopt label priors to
eliminate the CTC peaky behavior [15] for pure E2E word timing estimation and
the proposed system surpasses the performance of HMM-based system on Chinese
language. Our contributions in this paper are:

  1. 1.

We adopt label priors to eliminate the peaky behavior of CTC and train frame-
level classifiers to predict word timings, which can be easily extended to
other languages.

  2. 2.

We introduce Mel features as low-level local information and combine with
Transformer encoder outputs as high-level utterance-wise information, which
can be seen as feature fusion and is found to be beneficial for word timing
estimation.

  3. 3.

We further improve word timing accuracy by a frame-level knowledge
distillation method to shift word boundaries adaptively.

E2E speech recognition and word timing estimation systems are introduced in
Section 2. In Section 3, methods for optimizing the frame-level classifier
used in this paper are presented. Section 4 is the experimental setup. Results
and analyses are presented in Section 5, and a conclusion is given in Section
6.

##  2 System Overview

The E2E system used in this paper is illustrated in Figure 1. The ASR system
employs LAS [5] which is an attention-based encoder-decoder architecture. The
word timing system employs a frame-level classifier following the encoder and
the Mel feature is introduced optionally. The Mel filter bank sequence is
given as $\mathbf{x}=(x_{1},...,x_{T})$ and the wordpiece sequence is given as
$\mathbf{y}=(y_{1},...,y_{U})$, where $y_{u}$ is a wordpiece from a vocabulary
of the size $V$. The output of the decoder is denoted as
$P(y_{u}|\mathbf{x},y_{<u})$ at the step $u$. The output of the encoder is
denoted as $\mathbf{h}=(h_{1},...,h_{T})$ and the encoder output with Mel
feature denoted as $\mathbf{h}^{\prime}=(h_{1};x_{1},...,h_{T};x_{T})$
(represented by +Mel). The output of the frame-level classifier is denoted as
$P(v,t|\mathbf{x})$ of the wordpiece $v$ at the time $t$. The training stage
consists of two steps. First, the encoder-decoder is trained to maximize
$P(y_{u}|\mathbf{x},y_{<u})$. Second, the parameters in the LAS are fixed and
the frame-level classifier is trained using different methods which are
described in the next section. In the decoding stage, the encoder-decoder
produces ASR results $\mathbf{y}^{\prime}$ and the word timings come from the
forced alignment between $\mathbf{y}^{\prime}$ and $P(v,t|\mathbf{x})$.

Figure 1: system workflow

##  3 Methodology

The frame-level classifier trained with CTC loss has peaky behavior where the
peak roughly falls within the ground truth range of a token which can be
character, wordpiece, or word. The exact start and end times of a token can be
obtained by estimating the left and right boundaries based on CTC peak
position. In this section, word timing estimation including peak position
estimation and boundary estimation will be presented separately. In addition,
we also propose two methods to further polish word timing accuracy.

###  3.1 Peak position estimation

This subsection presents two methods combined to estimate the CTC peak
position.

CTC loss is proposed for training frame-level classifiers in E2E systems [6].
In order to align arbitrary frame-level and token-level sequences, the CTC
model introduces a blank token (represented by $\phi$). The CTC loss is
computed by Equation 1, where $P(\pi_{t}|\mathbf{x})$ is the posterior
probability of the token $\pi_{t}$ from the frame-level classifier at the time
$t$, and the function $\mathcal{B}$ removes consecutive equivalent tokens as
well as blank tokens so that $\mathcal{B}(\boldsymbol{\pi})=\mathbf{y}$. There
are many valid paths $\boldsymbol{\pi}$ that satisfy
$\mathcal{B}(\boldsymbol{\pi})=\mathbf{y}$, so the forward-backward algorithm
is used to optimize the calculation of the CTC valid paths. Since the E2E
system does not rely on the frame-level classifier to output ASR results, it
must use forced alignment to obtain word timings. Forced alignment, computed
by Equation 2, aims to find the most likely valid path $\boldsymbol{\pi}$. The
start and end time of each token $y_{u}$ is obtained through the path
$\boldsymbol{\pi}$. Many experiments have shown that frame-level classifiers
trained with CTC loss always predict blank tokens, known as the peaky behavior
[12, 16, 17]. As shown in Figure 2, forced alignment of non-blank tokens last
for few frames and peaks appear near the ground truths. To summarize, the
peaks can be obtained as anchor points for the position estimation by forced
alignment from the CTC model.

$\mathcal{L}_{\text{CTC}}=-\log
P(\textbf{y}|\textbf{x})=-\log\sum_{\boldsymbol{\pi}\in\mathcal{B}^{-1}(y)}\prod_{t}P(\pi_{t}|\mathbf{x})$
(1)
$\boldsymbol{\pi}=\operatorname*{argmax}_{\mathcal{B}(\boldsymbol{\pi})=\mathbf{y}}P(\boldsymbol{\pi}|\mathbf{x})=\operatorname*{argmax}_{\mathcal{B}(\boldsymbol{\pi})=\mathbf{y}}\prod_{t}P(\pi_{t}|\mathbf{x})$
(2)

Figure 2: CTC peaks

###  3.2 Boundary estimation

This subsection introduces two methods of peak expansion to estimate
boundaries.

A previous method to estimate boundaries is based on neighboring CTC peaks,
and teaches the frame-level classifier with a new posterior probability using
cross-entropy [13]. For convenience, this method is abbreviated to CETC. The
position of CTC peak of the token $y_{u}$ is denoted as $t_{\text{peak}}^{u}$.
The start and end times are computed by Equation 3. The probability of the
token $y_{u}$ at the time $t$ is $P(y_{u},t|\mathbf{x})$ from the classifier
trained with $\mathcal{L}_{\text{CETC}}=\text{Cross-
Entropy}(P,P_{\text{guided}})$, where $P_{\text{guided}}$ is computed by
Equation 4. The hyperparameters $\alpha_{\text{left}}$ and
$\alpha_{\text{right}}$ control the left and right boundaries, and $\beta$
controls that the further away from the CTC peak, the lower the probability.
According to the CETC paper, the
$(\alpha_{\text{left}},\alpha_{\text{right}},\beta)$ are set to
$(0.2,0.7,0.5)$. The CETC method estimating boundaries by neighboring peaks is
used to estimate word timings.

$\begin{split}t_{\text{start}}^{u}&=t_{\text{peak}}^{u}-\alpha_{\text{left}}\cdot(t_{\text{peak}}^{u}-t_{\text{peak}}^{u-1})\\\
t_{\text{end}}^{u}&=t_{\text{peak}}^{u}+\alpha_{\text{right}}\cdot(t_{\text{peak}}^{u+1}-t_{\text{peak}}^{u})\end{split}$
(3)
$\begin{split}P_{\text{guided}}(y_{u},t|\mathbf{x})&=\left\\{\begin{array}[]{ccl}(\frac{t-{t}_{\text{peak}}^{u}}{{t}_{\text{peak}}^{u}-{t}_{\text{start}}^{u}})^{\beta},&{t_{\text{start}}^{u}\leq
t<t_{\text{spike}}^{u}}\\\
(\frac{{t}_{\text{end}}^{u}-t}{{t}_{\text{end}}^{u}-{t}_{\text{peak}}^{u}})^{\beta},&{t_{\text{spike}}^{u}\leq
t\leq t_{\text{end}}^{u}}\end{array}\right.\\\
P_{\text{guided}}(\phi,t|\mathbf{x})&=1-\sum
P(\neq\phi,t|\mathbf{x})\end{split}$ (4)

A CTC model without peaky behavior is well suited for boundary estimation. A
recent paper explains and elegantly solves the CTC peaky behavior by
introducing the label prior, which is abbreviated to the non-peaky CTC (NPC)
method here [18]. The label prior must be calculated using logits before
softmax function, computed by Equation 5. The modified logits are computed by
Equation 6 and the NPC paper suggests using different $\gamma$ in training and
inference, which will be reported in the following section. The CTC loss is
calculated on the modified logits. The NPC method extending a peak to a span
by reducing the probability of blank token is employed to estimate the token
boundaries.

$\begin{split}O_{\text{prior}}(y_{u}|\mathbf{x})=\frac{1}{T}\sum_{t=1}^{T}O(y_{u},t|\mathbf{x})\end{split}$
(5)
$\begin{split}O^{\prime}(y_{u},t|\mathbf{x})=O(y_{u},t|\mathbf{x})-\gamma\cdot
O_{\text{prior}}(y_{u}|\mathbf{x})\end{split}$ (6)

###  3.3 Polishing word timing accuracy

In our experiments, we found that the NPC method tends to estimate word
boundaries earlier in time, which leaves room for improvement in word timing
accuracy. Having the manually labelled development set, the best shifting for
predicted boundaries can be grid searched and applied in post-processing. On
the other hand, to shift word timings adaptively, the CTC peaks are shifted
using a frame-by-frame knowledge distillation method called Peak First
Regularization (PFR) [10]. The PFR loss is computed by Equation 7, where
$\mu=-1$ delays CTC peaks, $\mu=1$ advances CTC peaks and $P_{\text{PFR}}$ is
the processed classifier probability. A temperature coefficient $\tau$ in
Equation 8 smooths the output of the frame-level classifier
$O(v,t|\mathbf{x})$ of the token $v$ at the time $t$. According to the
independence assumption of CTC, the PFR loss must be employed before freezing
the LAS and the loss is computed by Equation 9. It is necessary to prevent the
gradient from passing through the teacher frame
$P_{\text{PFR}}(v,t+\mu|\mathbf{x})$. The default
$(\mu,\tau,\lambda_{\text{CE}})$ is $(-1,10.0,0.95)$ when using the PFR loss.

$\mathcal{L}_{\text{PFR}}=\sum_{t}\sum_{v}P_{\text{PFR}}(v,t+\mu|\mathbf{x})\log\frac{P_{\text{PFR}}(v,t+\mu|\mathbf{x})}{P_{\text{PFR}}(v,t|\mathbf{x})}$
(7) $P_{\text{PFR}}(v,t|\mathbf{x})=\text{softmax}(O(v,t|\mathbf{x})/\tau)$
(8)
$\mathcal{L}=\lambda_{\text{CE}}\cdot\mathcal{L}_{\text{CE}}+(1-\lambda_{\text{CE}})\cdot(\mathcal{L}_{\text{CTC}}+\lambda_{\text{PFR}}\cdot\mathcal{L}_{\text{PFR}})$
(9)

##  4 Experiments

###  4.1 Datasets

We use the LibriSpeech 960h for base experiment and apply the word timing
estimation method to the internal data in 7 languages including Chinese (zh),
English (en), German (de), Russian (ru), Spanish (es), Indonesian (id) and
Vietnamese (vi) where the total durations are 100k, 40k, 9k, 10k, 20k, 6k and
10k hours respectively. The word timing test set for LibriSpeech is the forced
alignment result from HMM-based hybrid system with
script111https://github.com/kaldi-
asr/kaldi/blob/master/egs/librispeech/s5/run.sh and the internal test sets are
annotated by the experts. The E2E system employs 80-dimensional Mel-filter
bank features as input and wordpiece sequences obtained with SentencePiece
tokenizer [19] as output.

###  4.2 Modeling

In this paper, an encoder-decoder model outputs ASR results and a frame-level
classifier of 2 dense layers with 2048 hidden units following the encoder
outputs word timings. The 80-dimensional Mel-filter bank is computed on a
$25\text{\,}\mathrm{ms}$ window with a $10\text{\,}\mathrm{ms}$ shift. In the
encoder, two CNN kernels of width, height, input channel size, and output
channel size $(3,3,1,64)$ and $(3,3,64,128)$ with $(2,2)$ stride are used to
extract high-level feature. Then the CNN feature projected to 512-dimension is
fed into 18 transformer blocks. The transformer blocks have a multi-head self-
attention layer with 8 heads and a feed-forward layer adopting GLU activation
function with 2048 hidden units. The output of the last transformer block is
projected from 512 to 2048. The decoder has a multi-head attention layer with
2 heads, 4 LSTM layers with 1024 hidden units and an embedding layer with 512
hidden units. The encoder output optionally combines 4 consecutive frames of
the Mel features into 2368-dimensional features, which are fed to the frame-
level classifier as acoustic features. Overall the model has 168M parameters.

The E2E system is trained in Tensorflow [20] using the Lingvo [21] toolkit and
RETURNN [22] toolkit (for calculating the modified CTC loss) on 16 Tesla A100
GPUs for around 5 days. The hybrid system is trained in Kaldi [8].

###  4.3 Configuration

Using Adam optimizer to train the LAS model, the highest learning rate is
$7.5\mathrm{e}{-4}$ at the end of the warm-up at 4k steps. The learning rate
starts decaying at 60k steps decaying by half every 50k steps. A model is
trained 200k steps, and models for high-resource languages are trained up to
300k steps. Audios of similar length are grouped into same batches, and each
step of batch size multiplied by the longest audio in the batch is
approximately 1224 seconds. The training technique is from the LAS paper [5],
using Specaugment [23] where $(F,m_{F},T,p,m_{T})=(27,2,100,0.1,1)$. The
variational weight noise [24] is turned on optionally where standard deviation
is 0.075 and the uniform label smoothing [25] is introduced with uncertainty
0.1. In addition MWER [26] and CTC fine tuning with fixed learning rates
$1\mathrm{e}{-6}$ are trained for 20k and 10k steps respectively.

###  4.4 Measuring word timings

The word timing metrics are adopted from [2], which are the average offset
(ms) of the predicted word start to the ground truth (Ave. ST$\Delta$) and the
percentage of those word start offsets smaller than $200\text{\,}\mathrm{ms}$
(%WS<200ms). Similarly, the Ave. ED$\Delta$ and %WE<200ms for word end. In our
experiments, word timing metrics are only calculated on the matched token
pairs according to the alignment of minimum edit distance. For higher
granularity we employ %WS<80ms and %WE<80ms in our experiments.

##  5 Results

###  5.1 Analyzing the label prior

Table 1 shows the word timing metrics with different label priors used in
training and inference on LibriSpeech test-clean. The frame-level classifier
which is trained with CTC loss predicts mostly blanks and rarely the other
tokens, known as peaky behavior. The model with peaky behavior has bad word
timing metrics even though inference with the label prior when
$\gamma_{\text{train}}=0,\gamma_{\text{inf}}\neq 0$. It shows that the word
timing metrics have a significant improvement because of training with the
label prior as a normalization term when $\gamma_{\text{train}}\neq 0$. In the
following experiments, the default
$(\gamma_{\text{train}},\gamma_{\text{inf}})$ are (0.25, 1.0) because the
summation of the %WS<80ms and the %WE<80ms is the largest.

Table 1: Effect of the label prior ($\gamma$) in training and inference on word timing metrics on LibriSpeech $\gamma_{\text{train}}$ | $\gamma_{\text{inf}}$ | Ave. ST$\Delta$ | Ave. ED$\Delta$ | %WS<80ms | %WE<80ms
---|---|---|---|---|---
0.0 | 0.0 | 117.51 | 158.06 | 60.57 | 44.59
0.0 | 1.0 | 103.21 | 148.25 | 69.05 | 50.31
0.25 | 0.0 | 41.36 | 52.92 | 96.71 | 90.19
0.25 | 1.0 | 40.28 | 50.54 | 96.75 | 91.18
0.5 | 0.0 | 42.4 | 50.35 | 95.48 | 91.35
0.5 | 1.0 | 42.54 | 49.18 | 95.14 | 91.87
0.75 | 0.0 | 47.36 | 51.4 | 93.01 | 90.51
0.75 | 1.0 | 48.11 | 50.77 | 92.4 | 90.63
1.0 | 0.0 | 49.89 | 52.61 | 91.66 | 89.57
1.0 | 1.0 | 50.49 | 51.55 | 91.0 | 89.64

###  5.2 Improving frame-level classifiers in 7 languages

Table 2 shows the performance of 4 frame-level classifiers in Chinese and
compares the proposed method NPC+Mel to the previous work CETC in 6 other
languages. HMM denotes the frame-level classifier from HMM-based hybrid
system, and Offset denotes the offset in $\text{\,}\mathrm{ms}$ applied to all
timestamps. It shows that NPC outperforms HMM in Chinese when
$40\text{\,}\mathrm{ms}$ offset is applied. Further, the introduction of Mel
feature is to provide accurate local information to reduce the amount of
offset needed. Nevertheless it doesn't meet the purpose, it improves
performance. Hence, we suspect the cause of the word timing shift is not from
feature extraction but the use of CTC loss. NPC+Mel obtains an average of
86.79%/78.99% for %WS<80ms and %WE<80ms in all 7 languages, greatly
outperforms CETC's 82.14%/70.94% and slightly surpasses NPC's
86.53%/78.10%222Due to space limitation, the details of NPC results for 6
languages are not listed in Table 2.

Table 2: Word timing metrics in 7 languages in internal datasets

Method | Lang | Ave. ST$\Delta$ | Ave. ED$\Delta$ | %WS<80ms | %WE<80ms | Offset
---|---|---|---|---|---|---
HMM | zh | 41.47 | 47.05 | 93.0 | 90.22 | 0
CETC | zh | 39.62 | 44.53 | 91.86 | 89.72 | 0
NPC | zh | 32.69 | 35.39 | 95.5 | 93.07 | 40
NPC+Mel | zh | 31.93 | 33.78 | 95.68 | 94.18 | 40
CETC | en | 55.66 | 105.43 | 78.72 | 58.33 | 0
NPC+Mel | en | 44.05 | 70.35 | 88.64 | 77.91 | 60
CETC | ru | 53.6 | 92.63 | 82.18 | 61.62 | 0
NPC+Mel | ru | 48.49 | 80.49 | 85.2 | 69.68 | 60
CETC | es | 180.58 | 189.14 | 59.72 | 55.5 | 0
NPC+Mel | es | 173.7 | 181.4 | 62.53 | 60.62 | 70
CETC | de | 44.78 | 63.0 | 86.01 | 75.94 | 0
NPC+Mel | de | 33.62 | 52.64 | 93.61 | 83.7 | 60
CETC | id | 49.3 | 80.85 | 85.53 | 69.98 | 0
NPC+Mel | id | 46.23 | 80.62 | 88.11 | 76.28 | 60
CETC | vi | 38.48 | 47.92 | 90.99 | 85.49 | 0
NPC+Mel | vi | 30.09 | 39.62 | 94.85 | 90.36 | 40

###  5.3 Analyzing CTC peak positions

The introduction of CTC loss brings the peaky behavior, so we first analyze
the CTC peak positions. As shown in Figure 3, the number of CTC peaks is
measured as y-axis and the relative position of the peak time to the ground
truth word timings is measured as x-axis. $x=0$ means that CTC peaks are
emitted on the word start time and $x=-1$ means that CTC peaks appear a whole
word duration early before the start time. The red dotted dashed line in the
figure shows the average relative position of the peak to the ground truth is
0.249. The grey dashed lines indicate the average word start and end times of
the NPC+Mel prediction, and the CTC peaks are located on average at 0.406 of
the relative position of the prediction. The reason that applying positive
offset to the predicted timestamps increases accuracy is because the predicted
word start and end times are usually earlier than they should be, as shown in
Figure 3, the grey dashed lines do not lie exactly on $x=0$ and $x=1$. Table 3
shows that the average word duration of the prediction is very close to the
ground truth. In conclusion, the NPC+Mel method accurately predicts the word
duration, and by delaying the emission time of CTC peaks we might reduce the
amount of offset needed for increasing accuracy.

Figure 3: Distribution of the relative position of CTC peaks Table 3: Average
word duration of the ground truth and prediction ($\text{\,}\mathrm{ms}$) in
internal datasets

| zh | en | ru | es | de | id | vi
---|---|---|---|---|---|---|---
Ground truth | 163.3 | 287.1 | 376.6 | 287.7 | 289.5 | 368.0 | 201.0
Prediction | 162.8 | 257.2 | 331.6 | 288.7 | 279.3 | 337.5 | 195.8

###  5.4 Delaying CTC peaks and timestamps

We conduct the experiment in LibriSpeech. The CTC peak is delayed by
introducing PFR loss during training. Table 4 shows the performance with
different $\lambda_{\text{PFR}}$. The E2E system is trained from scratch with
CE, CTC and PFR loss. Because of the independence assumption in CTC loss, PFR
loss cannot change the position of CTC peaks in the case of fixed encoder
parameters. Without fixing the encoder parameters, the ASR result is affected
and we observe WER increases greatly as $\lambda_{\text{PFR}}$ surpasses
$2.5$. The result shows that setting $\lambda_{\text{PFR}}=1.5$ not only
reduces offset to $0\text{\,}\mathrm{ms}$ but also achieves best word timing
metrics.

Table 4: Effect of $\lambda_{\text{PFR}}$ on the offset needed to achieve the
highest accuracy on LibriSpeech

$\lambda_{\text{PFR}}$ | Ave. ST$\Delta$ | Ave. ED$\Delta$ | %WS<80ms | %WE<80ms | Offset | WER
---|---|---|---|---|---|---
0.0 | 40.28 | 50.54 | 96.75 | 91.18 | 40 | 5.3
0.5 | 35.06 | 42.04 | 97.35 | 93.84 | 20 | 4.1
1.0 | 34.73 | 40.85 | 98.0 | 95.18 | 10 | 4.7
1.5 | 38.69 | 42.9 | 97.99 | 96.04 | 0 | 5.0
2.0 | 32.57 | 35.96 | 97.89 | 94.82 | -10 | 4.2
2.5 | 40.36 | 41.81 | 97.84 | 95.59 | -20 | 4.7
3.0 | 51.62 | 55.69 | 97.67 | 93.96 | -30 | 6.3

##  6 Conclusions

In this paper, we improve frame-level classifiers with non-peaky CTC loss in a
purely E2E system. The proposed method outperforms the HMM-based hybrid system
in Chinese and the CETC method in 7 languages. We observe the timestamp
accuracy can be further improved, which requires manually labelled development
set for gridsearching the offset of timestamp in post-processing. Also we
employ PFR method to adaptively shift CTC peaks in training, improving word
timing performance. In the future, the proposed approach can be applied to a
multilingual E2E ASR system for word timings.

## References

  * [1] L. Huang, J. Sun, Y. Tang, J. Hou, J. Chen, J. Zhang, and Z. Ma, ``HMM-Free Encoder Pre-Training for Streaming RNN Transducer,'' in _Proc. Interspeech 2021_ , 2021, pp. 1797–1801.
  * [2] T. N. Sainath, R. Pang, D. Rybach, B. Garcıa, and T. Strohman, ``Emitting word timings with end-to-end models,'' _Proc. Interspeech 2020_ , pp. 3615–3619, 2020.
  * [3] S. Karita, N. Chen, T. Hayashi, T. Hori, H. Inaguma, Z. Jiang, M. Someki, N. E. Y. Soplin, R. Yamamoto, X. Wang _et al._ , ``A comparative study on transformer vs rnn in speech applications,'' in _2019 IEEE Automatic Speech Recognition and Understanding Workshop (ASRU)_. IEEE, 2019, pp. 449–456.
  * [4] A. Gulati, J. Qin, C.-C. Chiu, N. Parmar, Y. Zhang, J. Yu, W. Han, S. Wang, Z. Zhang, Y. Wu _et al._ , ``Conformer: Convolution-augmented transformer for speech recognition,'' _arXiv preprint arXiv:2005.08100_ , 2020.
  * [5] W. Chan, N. Jaitly, Q. Le, and O. Vinyals, ``Listen, attend and spell: A neural network for large vocabulary conversational speech recognition,'' in _2016 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2016, pp. 4960–4964.
  * [6] A. Graves, S. Fernández, F. Gomez, and J. Schmidhuber, ``Connectionist temporal classification: labelling unsegmented sequence data with recurrent neural networks,'' in _Proceedings of the 23rd international conference on Machine learning_ , 2006, pp. 369–376.
  * [7] A. Graves, ``Sequence transduction with recurrent neural networks,'' _arXiv preprint arXiv:1211.3711_ , 2012.
  * [8] D. Povey, A. Ghoshal, G. Boulianne, L. Burget, O. Glembek, N. Goel, M. Hannemann, P. Motlicek, Y. Qian, P. Schwarz, J. Silovsky, G. Stemmer, and K. Vesely, ``The kaldi speech recognition toolkit,'' in _IEEE 2011 Workshop on Automatic Speech Recognition and Understanding_. IEEE Signal Processing Society, Dec. 2011, iEEE Catalog No.: CFP11SRW-USB.
  * [9] A. Baevski, Y. Zhou, A. Mohamed, and M. Auli, ``wav2vec 2.0: A framework for self-supervised learning of speech representations,'' _Advances in neural information processing systems_ , vol. 33, pp. 12 449–12 460, 2020.
  * [10] Z. Tian, H. Xiang, M. Li, F. Lin, K. Ding, and G. Wan, ``Peak-first ctc: Reducing the peak latency of ctc models by applying peak-first regularization,'' _arXiv preprint arXiv:2211.03284_ , 2022.
  * [11] Y. Teytaut, B. Bouvier, and A. Roebel, ``A study on constraining connectionist temporal classification for temporal audio alignment,'' _Proc. Interspeech 2022_ , pp. 5015–5019, 2022.
  * [12] H. Liu, S. Jin, and C. Zhang, ``Connectionist temporal classification with maximum entropy regularization,'' _Advances in Neural Information Processing Systems_ , vol. 31, 2018.
  * [13] X. Chen, H. Ni, Y. He, K. Wang, Z. Ma, and Z. Xie, ``Emitting Word Timings with HMM-Free End-to-End System in Automatic Speech Recognition,'' in _Proc. Interspeech 2021_ , 2021, pp. 2571–2575.
  * [14] X. Shi, Y. Chen, S. Zhang, and Z. Yan, ``Achieving timestamp prediction while recognizing with non-autoregressive end-to-end asr model,'' _arXiv preprint arXiv:2301.12343_ , 2023.
  * [15] A. Zeyer, R. Schlüter, and H. Ney, ``Why does ctc result in peaky behavior?'' _arXiv preprint arXiv:2105.14849_ , 2021.
  * [16] Y. Miao, M. Gowayyed, and F. Metze, ``Eesen: End-to-end speech recognition using deep rnn models and wfst-based decoding,'' in _2015 IEEE Workshop on Automatic Speech Recognition and Understanding (ASRU)_. IEEE, 2015, pp. 167–174.
  * [17] H. Sak, A. Senior, K. Rao, O. Irsoy, A. Graves, F. Beaufays, and J. Schalkwyk, ``Learning acoustic frame labeling for speech recognition with recurrent neural networks,'' in _2015 IEEE international conference on acoustics, speech and signal processing (ICASSP)_. IEEE, 2015, pp. 4280–4284.
  * [18] A. Zeyer, R. Schlüter, and H. Ney, ``Why does CTC result in peaky behavior?'' _CoRR_ , vol. abs/2105.14849, 2021. [Online]. Available: https://arxiv.org/abs/2105.14849
  * [19] T. Kudo and J. Richardson, ``SentencePiece: A simple and language independent subword tokenizer and detokenizer for neural text processing,'' in _Proceedings of the 2018 Conference on Empirical Methods in Natural Language Processing: System Demonstrations_. Brussels, Belgium: Association for Computational Linguistics, Nov. 2018, pp. 66–71. [Online]. Available: https://www.aclweb.org/anthology/D18-2012
  * [20] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng, ``TensorFlow: Large-scale machine learning on heterogeneous systems,'' 2015, software available from tensorflow.org. [Online]. Available: https://www.tensorflow.org/
  * [21] J. Shen, P. Nguyen, Y. Wu, Z. Chen _et al._ , ``Lingvo: a modular and scalable framework for sequence-to-sequence modeling,'' 2019.
  * [22] A. Zeyer, T. Alkhouli, and H. Ney, ``RETURNN as a generic flexible neural toolkit with application to translation and speech recognition,'' in _Proceedings of ACL 2018, System Demonstrations_. Association for Computational Linguistics, 2018. [Online]. Available: https://doi.org/10.18653%2Fv1%2Fp18-4022
  * [23] D. S. Park, W. Chan, Y. Zhang, C.-C. Chiu, B. Zoph, E. D. Cubuk, and Q. V. Le, ``SpecAugment: A Simple Data Augmentation Method for Automatic Speech Recognition,'' in _Proc. Interspeech 2019_ , 2019, pp. 2613–2617.
  * [24] A. Graves, ``Practical variational inference for neural networks,'' in _Advances in neural information processing systems_. Citeseer, 2011, pp. 2348–2356.
  * [25] C. Szegedy, V. Vanhoucke, S. Ioffe, J. Shlens, and Z. Wojna, ``Rethinking the inception architecture for computer vision,'' in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 2818–2826.
  * [26] R. Prabhavalkar, T. N. Sainath, Y. Wu, P. Nguyen, Z. Chen, C.-C. Chiu, and A. Kannan, ``Minimum word error rate training for attention-based sequence-to-sequence models,'' in _2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2018, pp. 4839–4843.