\section{Introduction}
Almost three decades ago, Peter Shor~\cite{Shor1999} discovered a celebrated polynomial time quantum algorithm for integer factorisation. Although the best-known classical algorithm for integer factorisation scales quasi-exponentially, the problem's computational complexity has not been settled yet. Therefore, we cannot claim quantum supremacy~\cite{Harrow2017} clearly with Shor's quantum algorithm. For factoring an $L$-bit number, it requires $\mathcal{O}(L)$ qubits and $\mathcal{O}(L^{3})$ quantum gate operations mainly for the quantum Fourier transform (QFT). The current classical factoring limit is several hundred bits. Therefore, implementing Shor's quantum algorithm for factoring large integer numbers requires a fault-tolerant error-correcting universal quantum computer, which cannot be built in the near term due to technical challenges. Nevertheless, we cannot show quantum supremacy with Grover's search algorithm~\cite{Grover1996} because it can only achieve a quadratic speedup compared to the best-known classical algorithm. Consequently, other problems, which can provide clear computational complexity proofs and are implementable without expensive quantum subroutines like the QFT, are expected to demonstrate quantum computers' computational superiority against their classical counterparts. 

In this context, Aaronson and Arkhipov (AA)~\cite{Aaronson2011} suggested a quantum sampling problem~\cite{Lund2017,Harrow2017}, called boson sampling, that could potentially challenge the extended Church Turing thesis and demonstrate the quantum supremacy. The hardness of boson sampling is established on the computational complexity of matrix permanents. Approximating permanents of arbitrary complex or real square matrices with a constant factor is a \#P-hard problem. Boson sampling can be interpreted directly as a scattering process of non-interacting single photons in a linear optical network~\cite{Troyansky1996,Scheel2004}. The output photon distribution of the boson sampling device is given as a function of permanents of complex Gaussian submatrices. Using two conjectures, namely, permanent of Gaussian conjecture (PGC) and permanent anti-concentration conjecture (PAC), AA showed an efficient classical simulation of boson sampling, i.e. an additive error matrix permanent estimator, would collapse the polynomial hierarchy (PH) to the third level; this is strongly believed not to happen. Specifically, the classical hardness of boson sampling relies on the hardness of the permanents of random complex Gaussian matrices. The Gaussian matrices are hidden in a big Haar random unitary matrix of a linear optical network. Nevertheless, the \#P-hardness of permanents of Gaussian matrices has not yet been proved. The multiplicative error estimation of permanents of Gaussian matrices is conjectured as \#P-hard in PGC. Assuming a non-uniform distribution of permanents of random Gaussian matrices by PAC, it was shown that the multiplicative error estimation is equivalent to the additive error estimation for the permanents of Gaussian matrices. Recently, Pan et al. ~\cite{Zhongeabe8770} demonstrated the quantum advantage with the Gaussian boson sampling~\cite{Hamilton2017} setup that uses Gaussian photon sources instead of single-photon sources. Their quantum advantage arguments exploited the classical hardness of a matrix function called a Torontonian~\cite{Quesada2020}, a generating function of hafnians (or, equivalently, multivariate normal moments). Gaussian boson sampling is also believed to be a classically hard problem. It has some practical applications, contrary to the original boson sampling, in areas like molecular spectroscopy~\cite{huh2015}, graph~\cite{Bradler2018}, and molecular docking ~\cite{Banchieaax1950}. 

Similar complexity logics of boson sampling were applied to a quantum circuit model, instantaneous quantum polynomial time (IQP) problem~\cite{Shepherd2009}. The IQP circuit is composed of commuting gates. The resulting output qubit distribution cannot be efficiently determined using classical computers: the corresponding quantum overlap integrals can be interpreted as real-time partition functions of the quantum Ising Hamiltonian (Heisenberg model)~\cite{Matsuo2014,Fujii2017,Mann2019}. Google research~\cite{Arute2019} demonstrated the quantum supremacy with the random circuit sampling problem using their 53-superconducting-qubit device by introducing non-commuting quantum gates, unlike the IQP setup. However, the verification and performance remain  arguable~\cite{Pednault2019,Huang2020,NVIDIA}.

Regardless of these quantum sampling efforts, quantum algorithms for well-classified decision problems are needed for practical applications and quantum computational supremacy or advantage issues. The development of such quantum algorithms for decision problems can naturally lead to practical applications, unlike quantum sampling problems. Hopefully, we do not need to rely on the conjectures as to the sampling problems and the expensive QFT subroutine. 

Troyansky and Tishby~\cite{Troyansky1996} showed how a Hermitian matrix permanent is expressed with an expectation value of a linear optical quantum observable. They concluded that exponentially many samples are needed to compute a permanent because the observable variance is significant. Note that the boson sampling device cannot be used to estimate permanents as the magnitudes are generally exponentially small. There are only a few efforts, as yet, to relate quantum circuit models to matrix permanents. Most of the studies focused on the computational hardness proof of computing permanents via the quantum circuit models~\cite{Rudolph2009,Aaronson2011hard}, which did not result in practical quantum algorithms for estimating matrix permanents.  

Computing permanents of matrices with a constant factor is classically intractable (\#P-complete for complex matrices and \#P-hard for real matrices)~\cite{Aaronson2011}. It is believed, although without proof, that even quantum algorithms are not expected to access the \#P-hard problem. Nevertheless, developing a quantum algorithm for approximating the matrix permanent is essential to locating bounded error quantum polynomial time (BQP) relative to the PH, which is a significant complexity question because of its apparent computational complexity. To the best of our knowledge, there is no practical quantum algorithm reported for computing permanents so far.

We transformed the well-known classical algorithm, Ryser's formula~\cite{Ryser1963}, computing an $N$-dimensional  matrix permanent with an exponential cost $\mathcal{O}(N^{2}2^{N})$, to sums of $N$-th powers of a binary function. Then, we reformulated the power expansion as an $N$-th derivative of the quantum Ising Hamiltonian propagator. Consequently, we clearly showed how the boson sampling and the IQP problems are connected~\cite{Matsuo2014,Fujii2017,Mann2019} via expressing the matrix permanent as a linear combination of the real-time partition functions of the Ising Hamiltonian. Using the quantum power method of Seki et al.~\cite{Seki2021}, which exploits the finite difference scheme for the time derivatives, the matrix permanent estimation problem becomes an overlap integral evaluation problem. \emph{This paper presents the first efficient quantum algorithm for approximating the permanents of arbitrary real and complex matrices with multiplicative error and additive error protocols.} The multiplicative error protocol is suggested by combining the quantum amplitude estimation methods of Knill, Ortiz, and Somma (KOS)~\cite{Knill2007} and Aaronson and Rall (AR)~\cite{Aaronson2020}. The multiplicative error estimation of permanent implies that quantum computers can access the \#P-hard problem~\cite{Aaronson2011} (see Theorem 1 and Corollary 1). We also propose an additive error protocol using the Hadamard test~\cite{Aharonov2007,Aharonov2009} for the quantum overlap integral estimation. We compared our quantum algorithm to the best classical random sampling algorithm, i.e. the Gurvits algorithm~\cite{gurvits2005,Aaronson2011}, and found criteria for the quantum advantage. Our additive error quantum protocol requires a polynomial number of quantum gates and depth,  $\mathcal{O}(N^{2})$ for $2N+1$ qubits, that it is practically implementable in noisy intermediate-scale quantum (NISQ) devices. We also note that our quantum algorithm is potentially applicable to other matrix functions such as (loop)-hafnians~\cite{kan:2007,bjorklund2019}, Torontonians~\cite{Quesada2018}, determinants, and immanants, which are related to symmetries of quantum particles~\cite{Somma2003}.  

The rest of the paper is structured as follows: First, we transform the Ryser formula to a sum of $N$-th power of binary functions. Then, we show that a matrix permanent can be approximately cast into the sum of quantum overlap integrals. The multiplicative error and the additive error protocols are followed. Finally, after presenting quantum circuits for estimating a matrix permanent with an additive error using the Hadamard test, we discuss the implications of our findings and conclude the paper.    

\section{Classical Formulas}
Recall the definition of permanent of an $N$-dimensional square matrix $A$ (Per($A$)):
\begin{align}
\mathrm{Per}(A)=\sum_{\vec{{\sigma} }\in S_{\mathrm{p}}}\prod_{i=1}^{N}A_{j,\sigma_{j}} ,
\label{eq:PERA}
\end{align}
where $A_{jk}$ are the entries of $A$ and the set $S_{\mathrm{p}}$ includes all the  permutations of $(1,2,\ldots,N)$, $\{\vec{\sigma} \}$ ($\vec{\sigma} =(\sigma_1,\ldots \sigma_N)$ is an $N$-dimensional vector).  
Brute force computation of a matrix permanent in Eq. \eqref{eq:PERA} requires $N!$ terms in the summation, and each term is composed of products of  $N$ elements of the matrix.
In 1963, Ryser~\cite{Ryser1963} discovered a formula for matrix permanent, which involves only $2^{N}$ terms in the summation,  
\begin{align}
\textrm{Per}(A) = (-1)^{N}\sum_{S\subseteq \{1,\ldots,N\}}(-1)^{\vert S \vert}
\prod_{j=1}^N \sum_{k\in S}A_{k,j} ,
\label{eq:Ryser}
\end{align}
where $|S|$ is the cardinality of subset $S$. Although, this is a dramatic improvement, the algorithm is still not classically efficient as it scales exponentially $\mathcal{O}(N^22^{N})$. 
We can rewrite the Ryser formula~\eqref{eq:Ryser} by introducing an $N$-dimensional binary vector $\vec{y}$ as follows:  
\begin{align}
\textrm{Per}(A) = \sum_{\vec{y}\in \{0,1\}^{N}}(-1)^{\sum_{k=1}^{N}(y_{k}+1)}
\prod_{j=1}^N (\vec{A}_{j}\cdot\vec{y}) ,
\label{eq:Ryser2}
\end{align}
where $\vec{A}_{j}$ is the $j$-th row of $A$. 
By introducing a variable $x_{k}=2y_{k}-1$, which maps from $\{0,1\}^{N}$ to $\{-1,1\}^{N}$, we can find the Glynn  formula~\cite{Glynn2010,aaronson2012,Chin2018,Yung2019}: 
\begin{align}
\textrm{Per}(A) = \frac{1}{2^{N}}\sum_{\vec{x}\in \{-1,1\}^N}\left(\prod_{k=1}^{N}x_{k}\right)\left(\prod_{j=1}^N\vec{A}_{j}\cdot\vec{x}\right) ,
\label{eq:Glynn}
\end{align}
where we used the identity $\sum_{\vec{x}\in\{-1,1\}^{N}}\prod_{k}^{N}x_{k}=0$.  
Gurvits~\cite{gurvits2005,Aaronson2011} suggested an additive error algorithm $\pm\epsilon\Vert A \Vert^{N}$ based on Glynn's formula by interpreting it as an expectation value $E(\prod_{k=1}^{N}x_{k}\prod_{j=1}^N\vec{A}_{j}\cdot\vec{x})$ over the random binary vector $\vec{x}$. $\Vert \cdot \Vert$ is the matrix 2-norm; otherwise, the order is indicated.   

Using the identity of a monomial by Kan~\cite{kan:2007}, we have     
\begin{align}
    \prod_{j=1}^{N}\vec{A}_{j}\cdot\vec{x}=\frac{1}{N!2^{N}}\sum_{\vec{x}'\in\{-1,1\}^{N}}\left(\prod_{k=1}^{N}x'_{k}\right)\left(\vec{x'}^{\top}A\vec{x}\right)^{N} .
    \label{eq:monomial}
\end{align}
Now, we have a new expression for the matrix permanent by inserting Eq.~\eqref{eq:monomial} into Eq.~\eqref{eq:Glynn} as follows: 
\begin{align}
    \mathrm{Per}(A)=\frac{1}{N!2^{2N}}\sum_{\vec{x},\vec{x'}\in\{-1,1\}^{N}}\left(\prod_{k=1}^{N}x_{k}x_{k}'\right)\left(\vec{x'}^{\top}A
    \vec{x}\right)^{N} .
    \label{eq:Kanexpression}
\end{align}
We name our new finding in  Eq.~\eqref{eq:Kanexpression} the Glynn-Kan formula in this paper. 
Thus, Per($A$) can be interpreted as an average of $2^{2N}$ $N$-th power of binary functions by  Eq.~\eqref{eq:Kanexpression}. Interestingly, Eq.~\eqref{eq:Kanexpression} can be rewritten as a difference between two positive semidefinite functions for a real square matrix of $B$ with an even dimension, 
\begin{align}
    \mathrm{Per}(B)=\frac{1}{N!2^{N}}\left[\frac{1}{2^{N}}\sum_{\{\vec{x},\vec{x'}\}\in \mathcal{S}_{+}}\left(\vec{x'}^{\top}B
    \vec{x}\right)^{N}
    -\frac{1}{2^{N}}\sum_{\{\vec{x},\vec{x'}\}\in \mathcal{S}_{-}}\left(\vec{x'}^{\top}B
    \vec{x}\right)^{N}\right],
\end{align}
where $\mathcal{S}_{\pm}$ refers to a subset of $\left(\prod_{k=1}^{N}x_{k}x_{k}'\right)=\pm 1$, respectively, and each subset is composed of $2^{N}$ elements. 
Because  $\left(\prod_{k=1}^{N}x_{k}x_{k}'\right)$ can be evaluated efficiently, and we know $\mathcal{S}_{\pm}$ are evenly distributed, we could evaluate permanents of real matrices of even dimension with multiplicative error according to Stockmeyer's counting algorithm~\cite{Aaronson2011,Stockmeyer1985}, which would imply that this problem is classified as  BPP$^{\mathrm{NP}}$, where BPP and NP stand for bounded error probabilistic polynomial time and nondeterministic polynomial time, respectively. For the odd-dimensional case, we can use the following identity to have the even-dimensional matrix, $\mathrm{Per}(B)= \mathrm{Per}\left(B\bigoplus 1\right)$. Here, we express the matrix permanent as a difference between the two expectation values, 
\begin{align}
    \mathrm{Per}(B)=\frac{1}{N!2^{N}}\left[\underset{\{\vec{x},\vec{x'}\}\in \mathcal{S}_{+}}{E}\left[\left(\vec{x'}^{\top}B
    \vec{x}\right)^{N}\right]  
    -\underset{\{\vec{x},\vec{x'}\}\in \mathcal{S}_{-}}{E}\left[\left(\vec{x'}^{\top}B
    \vec{x}\right)^{N}\right]\right] .
    \label{eq:diffexpectation}
\end{align}
According to Stockmeyer's theorem, the even and odd phase expectation values can be evaluated with the multiplicative errors $h_{\pm}=1+1/\mathrm{poly}_{\pm}(N)$, respectively, where ‘poly’ denotes a polynomial function. If we can efficiently find an identical multiplicative factor for the even and odd phase expectation values, i.e. $h_{+}=h_{-}=h$, permanents of real matrices can be efficiently estimated with a multiplicative factor $h$ using a $\mathrm{BPP^{NP}}$ algorithm. 
However, its implication would result in the collapse of the PH to the third level, i.e. $\mathrm{P^{\#P}}=\mathrm{BPP^{NP}}$. 
Therefore, we generally expect that $h_{+}\neq h_{-}$. 
For exceptional cases, however, such as real matrices having non-negative entries and Hermitian positive semidefinite matrices, we could find a common multiplicative factor $h$. Meanwhile, the multiplicative error estimation of the permanent of  a Hermitian positive semidefinite matrix is classified as   $\mathrm{BPP^{NP}}$~\cite{Lund2014,Rahimi-Keshari2015,Chakhmakhchyan2017}. Although it is known that permanents of non-negative matrices can be evaluated within a constant factor by Jerrum-Sinclair-Vigoda (JSV) algorithm~\cite{Jerrum2004} in polynomial time, the JSV algorithm has no practical implementations so far~\cite{Bezakova2008,Newman2020}. 
Therefore, further computational complexity analysis of estimating permanents of real matrices using Eq.~\eqref{eq:diffexpectation} is required because it is known that computing a permanent of an arbitrary real matrix is a \#P-hard problem~\cite{Aaronson2011} and we would be able to discover new efficient classical algorithms with it.     

When $A$ is a complex matrix, i.e.   $A=B+\mathrm{i}C$, where  $B$ and $C$ are real square matrices, we can expand Eq.~\eqref{eq:Kanexpression} as follows: 
\begin{align}
\mathrm{Per}(A)
    =\frac{1}{N!2^{2N}}\sum_{l=0}^{N}\mathrm{i}^{l}
    \begin{pmatrix}
    N\\l
    \end{pmatrix}
    \sum_{\vec{x},\vec{x'}\in\{-1,1\}^{N}}\left(\prod_{k=1}^{N}x_{k}x_{k}'\right)\left(\vec{x'}^{\top}
    B
    \vec{x}\right)^{N-l}
    \left(\vec{x'}^{\top}
    C
    \vec{x}\right)^{l} .
   
   
   
   
   
   
   
   
   
   
   
   
    \label{eq:Kanexpressioncomplex}    
\end{align}
The exact evaluation of Eq.~\eqref{eq:Kanexpressioncomplex} for the permanent of a complex matrix requires the summation of $(N+1)2^{2N}$ real binary functions with complex weights. At a glance, our new formulas for permanents,  Eqs.~\eqref{eq:Kanexpression} and \eqref{eq:Kanexpressioncomplex}, are exponentially more inefficient compared with the original Ryser's formula, Eq.~\eqref{eq:Ryser}. However, our new formulas are suitable for developing a BQP quantum algorithm for permanents because they are composed of an expectation value of $N$-th powers of the binary function, i.e. $\vec{x'}^{\top}A
    \vec{x}$. 
We use these formulas to develop quantum algorithms for estimating the permanents of arbitrary square matrices. 
In the following section, we develop quantum algorithms for multiplicative error and additive error approximation of permanents of real and complex matrices. Eqs.~\eqref{eq:Kanexpression} and \eqref{eq:Kanexpressioncomplex} are approximately evaluated using a finite difference method and quantum overlap integrals of the quantum Ising Hamiltonian propagators. Thus, our quantum algorithm for permanents exhibits a direct connection between the boson sampling and IQP. We propose two overlap integral evaluation protocols of multiplicative and additive errors.    

\section{Quantum algorithm}
From the Glynn-Kan formula in Eq.~\eqref{eq:Kanexpression}, we  immediately identify the following quantum expectation value as a matrix permanent: 
\begin{align}
\textrm{Per}(A) = 
\frac{1}{N!2^{2N}}\sum_{y=0}^{2^{2N}-1}
\langle y\vert\left(\bigotimes_{k=1}^{N}Z_{k}Z_{N+k}\right)
\left(\sum_{j,k=1}^{N}A_{jk}Z_{N+j}Z_{k}\right)^{N}\vert y \rangle ,
\label{eq:Glynnexp}
\end{align}
where $\vert y\rangle=\vert y_{1}\rangle\cdots\vert y_{2N}\rangle$, $y=y_{1}\cdots y_{2N}$ is given as a $2N$-bit string and $Z_{k}$ is a Pauli-Z operator on the $k$-th qubit. Using the identity, $\frac{1}{2^{2N}}\sum_{y=0}^{2^{2N}-1}
\vert y\rangle=\bigotimes_{k=1}^{2N}H_{k}\vert 0\cdots 0\rangle=\vert\phi\rangle$, where $H_{k}$ is a Hadamard operator of the $k$-th qubit, we rewrite the matrix permanent as an expectation value of a quantum operator, i.e. $\mathrm{Per}(A)=\frac{1}{N!}\langle\phi\vert G(A)\vert\phi\rangle$. 
Let us call $G(A)$ the Glynn-Kan operator, defined as   
\begin{align}
G(A)=\left(\bigotimes_{k=1}^{N}Z_{k}Z_{N+k}\right)
\left(\sum_{j,k=1}^{N}A_{jk}Z_{N+j}Z_{k}\right)^{N},  
\end{align}
which is factored with two components--parity operator and Ising operator, respectively denoted as $P$ and $\mathcal{H}$. That is, 
$G(A)=P\mathcal{H}^{N}(A)$, where 
$\mathcal{H}(A)=\sum_{j,k=1}^{N}A_{jk}Z_{N+j}Z_{k}$ and 
$P=\bigotimes_{k=1}^{N}Z_{k}Z_{N+k}$.
Thus far, we have expressed the matrix permanent as an expectation value of qubit operators and qubit states. The resulting expression involves a power operator of $\mathcal{H}$, which is generally non-Hermitian. The remaining task is to approximately implement the power operator in a quantum circuit. We can perform the approximation through different approaches. For example, we can consider quantum signal processing with quantum singular value transformation and the linear combination of unitaries (LCU) technique~\cite{Childs2012,Berry2015,Low2019}. However, the repeated application of LCU for the power operator requires many ancilla qubits, and the success rate scales badly. 
In this paper, we adopt the approach of Seki et al.~\cite{Seki2021}, the quantum power method, a quantum version of the finite difference method for time-derivatives of a quantum propagator, because it results in simpler quantum circuits. 

\subsection{Unitary Approximation}

\noindent\textbf{Real Matrix}:
When $A$ is a real matrix, i.e. $A=B$, $\mathcal{H}(B)$ is a Hermitian operator that we can express the power operator as the time derivative of the corresponding time propagator $U(B;t)=\exp(-\mathrm{i}\mathcal{H}(B)t)$ directly as  
\begin{align}
  \mathcal{H}(B)^{N}=\mathrm{i}^{N}\frac{\partial^{N}U(B;t)}{\partial t^{N}}\Big\vert_{t=0}.   
\end{align}
We can approximately implement the time derivatives using the finite difference method by Seki et al.~\cite{Seki2021}, as follows, with the time step $\Delta t$, 
\begin{align}
    \frac{\mathrm{i}^{N}}{\Delta t^{N}}\left[U\left(B;\Delta t/2\right)-U\left(B;-\Delta t/2\right)\right]^{N}
    &=\frac{\mathrm{i}^{N}}{\Delta t^{N}}
    \sum_{l=0}^{N}(-1)^{l}\binom{N}{l}U\left((N-2l)B;\Delta t/2\right)  \\
    &=
    \mathcal{H}(B)^{N}-\frac{N\Delta t^{2}}{24}\mathcal{H}(B)^{N+2}+\cdots
    , 
    \label{eq:FDforReal}
\end{align}
where $\Vert \mathcal{H}(B) \Vert \Delta t/2\le 1$ is the convergence condition for the expansion. Thus, the Glynn-Kan operator is approximately given as a linear combination of unitary gates as follows:
\begin{align}
G(B)\simeq\frac{(-1)^{N}}{\Delta t^{N}}\sum_{l=0}^{N}(-1)^{l}\binom{N}{l}U\left(B^{(l)};\Delta t/2 \right),
\end{align}
where $B^{(l)}=(N-2l)B+\pi\Delta t^{-1} I$ and $I$ is an identity matrix. We used the following identity for the parity operator,  $P=\mathrm{i}^{N}\bigotimes_{k=1}^{N}\exp\left(-\frac{\pi}{2}\mathrm{i} Z_{k}Z_{N+k}\right)$, 
 to hide it in the unitary propagator $U$. By hiding the parity operator, we can simplify the required quantum circuits.  
Finally, the permanent of the real matrix $B$ is approximated as the sum of $N+1$ quantum overlap integrals,
\begin{align}
    \mathrm{Per}(B)&=\frac{(-1)^{N}}{N!\Delta t^{N}}\sum_{l=0}^{N}(-1)^{l}\binom{N}{l}\braket{\phi\vert U\left(B^{(l)};\Delta t/2 \right)\vert \phi} +
     \mathcal{O}(\frac{N\Delta t^{2}}{24N!}\Vert\mathcal{H}(B)\Vert^{N+2})
    .
    \label{eq:permanentofreal}
\end{align}
For the quantum algorithm's practical application, we can systematically improve  
the error introduced by the finite difference with the Richardson extrapolation method~\cite{Seki2021} for the given time step $\Delta t$.   
Eq.~\eqref{eq:permanentofreal} can be further simplified using the time reversal symmetry, $[\mathrm{i}^{N}\hat{U}(-\Delta t/2)]=[\mathrm{i}^{N}\hat{U}(\Delta t/2)]^{^{\dagger}}$. We only need to estimate a half ($\floor{(N+1)/2}$) of the real parts of the quantum overlap integrals, respectively, as follows: 
\begin{align}
    \mathrm{Per}(B)
    &\simeq
    \frac{2(-1)^{N}}{N!\Delta t^{N}}\sum_{l=0}^{\floor{\tfrac{N-1}{2}}} (-1)^{l}\binom{N}{l}
    \mathrm{Re}\braket{\phi\vert U\left(B^{(l)};\Delta t/2 \right)\vert \phi}
    .
    \label{eq:permanentofrealpart}
\end{align}



\noindent\textbf{Complex Matrix}: When $A$ is a complex matrix, i.e.  $A=B+\mathrm{i}C$, we exploit the classical expression in Eq.~\eqref{eq:Kanexpressioncomplex} for the quantum version as follows:
\begin{align}
\mathrm{Per}(A)
    &=\frac{1}{N!}\sum_{l=0}^{N}\mathrm{i}^{l}
    \begin{pmatrix}
    N\\l
    \end{pmatrix}
    \braket{\phi\vert
    P
    \mathcal{H}^{N-l}(B)
    \mathcal{H}^{l}(C)\vert\phi} 
    \\
    &=
    \frac{(-1)^{N}}{N!\Delta t^{N}}\sum_{l=0}^{N}
    \sum_{j=0}^{N-l}
    \sum_{k=0}^{l}
    \mathrm{i}^{l}
    (-1)^{j+k}
    \binom{N}{l}
    \binom{N-l}{j}
    \binom{l}{k}
    \braket{\phi\vert
    U(A^{(l,j,k)};\Delta t/2)\vert \phi} \nonumber \\
    &+\mathcal{O}\left(\frac{N\Delta t^{2}}{24N!}\Vert\mathcal{H}(A)\Vert^{N-1}(\Vert \mathcal{H}(B)\Vert^{3}+\Vert \mathcal{H}(C)\Vert^{3} )\right),
    \label{eq:Kanexpressioncomplexquantum}    
\end{align}
where $A^{(l,j,k)}=(N-l-2j)B+(l-2k)C+\pi\Delta t^{-1}I$; we used the following expansion:
\begin{align}
    &\frac{\mathrm{i}^{N}}{\Delta t^{N}}\left[U\left(B;\frac{\Delta t}{2}\right)-U\left(B;-\frac{\Delta t}{2}\right)+\mathrm{i}\left(U\left(C;\frac{\Delta t}{2}\right)-U\left(C;-\frac{\Delta t}{2}\right)\right)\right]^{N}\\
    &=
    \mathcal{H}(A)^{N}-\frac{N\Delta t^{2}}{24}\mathcal{H}(A)^{N-1}\left(\mathcal{H}(B)^{3}+\mathrm{i}\mathcal{H}(C)^{3}\right)+\cdots    , 
     \label{eq:FDforComplex}
\end{align}
which has the convergence condition 
\begin{align}
    \Vert \mathcal{H}(A) \Vert \Delta t/2\le
    \left(\Vert \mathcal{H}(B) \Vert+\Vert \mathcal{H}(C) \Vert\right) \Delta t/2
    \le 2 ;
\end{align}
thus, $\Vert \mathcal{H}(A) \Vert \Delta t/4\le 1$. 
Approximating the permanent of a complex matrix requires $(N^{3}+6N^{2}+11N+6)/6$ quantum overlap integrals in the summation, which is still a polynomial order but two orders more expensive than the real matrix case, whereas the circuit complexity is identical to the real case.  
    Here, we also only need to estimate half of the real part of the quantum overlap integrals:  
     $(N^{3}+6N^{2}+11N+6)/12$ and  $(N^{3}+6N^{2}+8N)/12$ overlap  integrals are to be evaluated for the odd and even dimensions, respectively. It is straightforward to reduce the number of terms for the complex matrix case as the real matrix one, but here, we do not present a simple expression for the complex case as in  Eq.~\eqref{eq:permanentofrealpart}.

\subsection{Overlap Integral Estimation}
Now the matrix permanent estimation problem becomes the evaluation problem of the real part of the quantum overlap integral, i.e. $\mathrm{Re}\braket{\phi\vert U \vert \phi}$, for both real and complex matrices. In the following two subsections, we present two overlap integral estimation protocols. One is a multiplicative error protocol based on KOS~\cite{Knill2007}, and AR~\cite{Aaronson2020}. The modified protocol in this paper can estimate the overlap integrals with a constant factor, which is essential for small overlap integrals. However, it requires more extensive quantum resources than the additive error estimation method. Another protocol uses the Hadamard test~\cite{Aharonov2007,Aharonov2009}, which generates additive errors. 

\subsubsection{multiplicative error protocol}
Brassard et al.~\cite{Brassard2002} proposed a quantum amplitude estimation algorithm. Here, the quantum amplitude refers to the absolute value of an overlap integral, $f$. The amplitude can be obtained within an additive error precision ($\epsilon=1/\mathrm{poly(N)}$) by estimating the eigenvalues of the Grover diffusion operator of $U$ using the quantum phase estimation technique. It requires $\mathcal{O}(1/(\epsilon f))$ queries to the Grover diffusion operator. Other Grover search-like algorithms were suggested by multiple groups~\cite{Wie2019,Grinko2021,Suzuki2020, Aaronson2020} to avoid the expensive QFT subroutines. AR developed such an algorithm, which can perform multiplicative error estimation of quantum amplitude with a similar number of quires as in Brassard et al.~\cite{Brassard2002}. On the other hand, KOS proposed a quantum protocol to estimate quantum overlap integrals as a function of multiple quantum amplitudes using Brassard et al.'s quantum algorithm. Here, we propose a modified quantum protocol to obtain the multiplicative error estimation of the quantum overlap integral by combining the KOS and AR methods. In generalizing the KOS method, we define a quantum amplitude ($f(\theta)$) involving an ancilla qubit $c$ with its rotation angle $\theta$ as follows:
\begin{align}
    f(\theta)=\vert\braket{\psi\vert R_{Z}^{(c)}(\theta)\otimes V\vert \psi}\vert,
\end{align}
where $\ket{\psi}=H^{(c)}\ket{0^{(c)}}\otimes \ket{\phi}$, $R_{Z}^{(c)}(\theta)=\exp(-\mathrm{i} \tfrac{\theta}{2} Z^{(c)})$, and $V=\ket{0^{(c)}}\bra{0^{(c)}}+\ket{1^{(c)}}\bra{1^{(c)}}\otimes U$. Then, we use the following function to estimate the real part of the quantum overlap integrals: 
\begin{align}
    \Delta f(\theta_{1},\theta_{2})=f(\theta_{1})^{2}-f(\theta_{2})^{2}
    =\frac{1}{2}\left[\mathrm{Re}\braket{\phi\vert U \vert \phi}(\mathrm{cos}(\frac{\theta_{1}}{2})-\mathrm{cos}(\frac{\theta_{2}}{2}))-\mathrm{Im}\braket{\phi\vert U \vert \phi}(\mathrm{sin}(\frac{\theta_{1}}{2})-\mathrm{sin}(\frac{\theta_{2}}{2}))\right].
\end{align}
Now, we can obtain the real part of the quantum overlap integral as $\Delta f(0,2\pi)=\mathrm{Re}\braket{\phi\vert U \vert \phi}$, where $f(\theta_{1})$ and $f(\theta_{2})$ are estimated using the AR algorithm with a multiplicative error constant $g/2$, yielding a multiplicative error factor $(1\pm g/2)^{2}\simeq 1\pm g$ for small $g$. The query complexity of the AR algorithm is $\mathcal{O}(g^{-1} f(\theta)^{-1}\log (1/\delta))$ with failure probability $\delta$. The AR algorithm does not work when the quantum amplitude $f(\theta)$ is too small as the required number of queries diverges. The failure occurs for evaluating $\Delta f(0,2\pi)$ when $\vert\mathrm{Re}\braket{\phi\vert U \vert \phi}\vert~\simeq 1$. Fortunately, we can avoid this unwanted situation efficiently. We can set a threshold for the oracle iteration number to test whether $f(\theta)$ is too small. When $f(\theta)<<1$, we can evaluate $\sqrt{2}\Delta f(\pi/2,3\pi/2)$ for $\mathrm{Re}\braket{\phi\vert U \vert \phi}$ instead of $\Delta f(0,2\pi)$. Therefore, in any case, we can obtain the real part of the overlap integrals with a constant factor $1\pm g$ via the AR amplitude estimation algorithm efficiently in a polynomial time. Note that the corresponding Grover diffusion operator oracle for $U$ can be implemented efficiently in polynomial depth quantum circuits of $2N+3$ qubits~\cite{Aaronson2020}. This implies that the BQP quantum algorithm can estimate the overlap integral with a constant factor.
Without loss of generality and for simplicity, we mainly focus on estimating the permanent of a real matrix with a multiplicative error since the protocol is directly applicable to complex matrices. 
To transform the error introduced by the finite difference approximation as a multiplicative form, we reexpress Eq.~\eqref{eq:FDforReal} as follows:
\begin{align}
    \frac{\mathrm{i}^{N}}{\Delta t^{N}}\left[U\left(B;\Delta t/2\right)-U\left(B;-\Delta t/2\right)\right]^{N}
    =
    \mathcal{H}(B)^{N}\left[1-\frac{N}{6}\left(\frac{\Delta t\mathcal{H}(B)}{2}\right)^{2}+\cdots\right]
    .
\end{align}
Consequently, we have the following expression for the matrix permanent scaled with a constant error factor caused by the finite difference: 
\begin{align}
    \mathrm{Per}(B)
    \simeq
    \frac{2(-1)^{N}}{N!\Delta t^{N}}\sum_{l=0}^{\floor{\tfrac{N-1}{2}}} (-1)^{l}\binom{N}{l}
    \mathrm{Re}\braket{\phi\vert U\left(B^{(l)};\Delta t/2 \right)\vert \phi}
    (1+g_{\mathrm{FD}})
    ,
\end{align}
where $g_{\mathrm{FD}}=N\vert\mathrm{sinc} x -1\vert$ and $x=\Delta t\Vert\mathcal{H}(B)\Vert/2\le 1$. $g_{\mathrm{FD}}$ is bounded as  $g_{\mathrm{FD}}< Nx^{2}/6$ within the domain $x\le 1$, which means $g_{\mathrm{FD}}< 1/\mathrm{poly}(N)$ when $x\le \sqrt{6}/\sqrt{N \mathrm{poly}(N)}$. Therefore, with a polynomially small value of $x$, we have the polynomially small error factor $g_{\mathrm{FD}}$ from the finite difference scheme. Finally, we conclude that the matrix permanent can be estimated as $\tilde{P}$ with a multiplicative error, which is known to be a \#P-hard problem~\cite{Aaronson2011}, via the KOS and AR algorithms, i.e. 
\begin{align}
    \vert\mathrm{Per}(B)\vert(1+g_{\mathrm{FD}})(1-\tfrac{g}{2})^{2}\le 
    \vert&\tilde{P}\vert
    \le
    \vert\mathrm{Per}(B)\vert(1+g_{\mathrm{FD}})(1+\tfrac{g}{2})^{2} \nonumber \\
    \Longrightarrow\vert\mathrm{Per}(B)\vert(1-g_{\times})\le 
    \vert&\tilde{P}\vert
    \le
    \vert\mathrm{Per}(B)\vert(1+g_{\times}),  
\end{align}
where $g_{\times}=(1+g_{\mathrm{FD}})(1+g/2)^{2}-1$. 
There is no restriction on the size of the matrix norm for efficient estimation of matrix permanents because we can control the error with $\Delta t \le 2\sqrt{6}/(\Vert \mathcal{H}(B)\Vert\sqrt{N \mathrm{poly}(N)})$.
A similar conclusion can be made for complex matrices. The existence of such a multiplicative error quantum algorithm for matrix permanent would lead to the surprising conclusion that BQP would equal $\mathrm{P^{\#P}}$ because we know  $\mathrm{BQP}\subseteq \mathrm{P^{\#P}}$~\cite{Bernstein1997,Aaronson2010}. Here, we do not definitively claim that $\mathrm{BQP}=\mathrm{P^{\#P}}$; further computational complexity analysis is required with our quantum algorithm for the conclusion. The bottom line is that the constant factor estimation of the real parts of the overlap integrals will lead us to the constant factor estimation of matrix permanents.  
We note here that the constant factor algorithm assumes the ideal quantum device, and only the statistical error was considered in the analysis. 
Armed with our multiplicative error protocol, we have the following theorem and corollary: 
\begin{theorem}
$\mathrm{(PH\subseteq BQP)}$
Given matrix $A\in \mathbb{R}^{N\times N}$ or $A\in \mathbb{C}^{N\times N}$, and an ideal BQP machine, if the estimation of $\mathrm{Per}(A)$ with a multiplicative factor $(1+ g_{\times})$ for $0<g_{\times}<1$ is \#P-hard, then  $\mathrm{PH}\subseteq\mathrm{BQP}$.
\end{theorem}
\emph{Proof.} We can estimate the permanent of an arbitrary matrix with a constant factor with the quantum algorithm suggested above in polynomial time. If the problem is \#P-hard, $\mathrm{BQP}=\mathrm{P^{\#P}}$. According to Toda's theorem ($\mathrm{PH}\subseteq\mathrm{P^{\#P}}$)~\cite{Toda}, the BQP contains the PH. 
\begin{corollary}
$\mathrm{(BPP\neq BQP)}$
If~$\mathrm{BPP}=\mathrm{BQP}$, then PH collapses to the second level according to Theorem 1 and Toda's theorem.  
\end{corollary}
Since we strongly believe that the PH does not collapse~\cite{Aaronson2011}, Corollary 1 tells us that   $\mathrm{BPP}\neq\mathrm{BQP}$. This implies that quantum computing is more powerful than classical computing if the PH does not collapse.  

\subsubsection{additive error protocol} 
We proposed a multiplicative error protocol for estimating a matrix permanent in the previous subsection. However, the proposed quantum protocol based on the AR algorithm is impractical due to the large polynomial-size overhead; therefore, not suitable for NISQ devices. We expect that other Grover-like quantum amplitude estimation algorithms~\cite{Wie2019,Grinko2021,Suzuki2020}, which require fewer quantum resources, would be helpful in practice, but they all perform the additive error estimation. In this subsection, we propose an additive error protocol, which can be implemented in simpler quantum circuits and is suitable for NISQ devices.    
We can obtain the real parts of quantum overlap integrals via the ancilla qubit measurement of the Hadamard test ~\cite{Aharonov2007,Aharonov2009} in Fig.~\ref{fig:quantumcircuits}(a). 
\begin{figure}[t]
\centering
    \includegraphics[width=1.0\textwidth]{quantumcircuits.jpg}
\caption{Quantum circuits for evaluating a real part of the overlap integral using the Hadamard test. (a) The Hadamard test circuit for the additive error  ($1/\mathrm{poly(N)}$) estimation of the real part of the quantum overlap integral $\mathrm{Re}\braket{\phi\vert U\vert\phi}$. The real part of the quantum overlap integral can be obtained via the difference between the probabilities of getting $0$ and $1$ in the ancilla qubit (denoted as $c$) measurement, i.e. $\mathrm{Pr}(0)-\mathrm{Pr}(1)$. The detailed quantum circuit  implementation is presented: (b) the circuit synthesis of the controlled-$R_{ZZ}(\theta)=\exp(-\mathrm{i}\tfrac{\theta}{2} Z_{1}Z_{2})$ gate and   
    (c) the circuit synthesis of controlled-$R_{Z}(\theta)=\exp(-\mathrm{i}\tfrac{\theta}{2} Z)$. 
    }
    \label{fig:quantumcircuits}
\end{figure}
The number of samplings for each real part is   $\mathcal{O}\left(\epsilon_{\mathrm{HT}}^{-2}\log(1/\delta)\right)$, where  $\epsilon_{\mathrm{HT}}$ and  $\delta$ are the required precision and a failure probability of the Hadamard test~\cite{Aharonov2007,Matsuo2014,Huang2019}, respectively. The errors in our quantum algorithm are caused by the finite difference method for the differentiation and the Hadamard test for the quantum overlap integral additive error estimation. In the following, we analyse the additive errors for (i) arbitrary real and complex matrices compared with the Gurvits classical sampling algorithm, and (ii) real and complex Gaussian matrices in connection to the boson sampling problem. The total additive error of a real matrix permanent $\epsilon_{\mathrm{Tot}}$ has an upper bound as follows: 
 \begin{align}
\epsilon_{\mathrm{Tot}}&\le\frac{1}{N!}\left(\epsilon_{\mathrm{HT}}\left(\frac{2}{\Delta t}\right)^{N} +\epsilon_{\mathrm{FD}}\frac{N\Delta t^{2}}{24}\Vert\mathcal{H}(B)\Vert^{N+2}\right),
\end{align}
where 
$0<\epsilon_{\mathrm{HT}}$,    $\epsilon_{\mathrm{FD}}=1/\mathrm{poly}(N)<1$. Here,  
$\epsilon_{\mathrm{HT}}$ and $\epsilon_{\mathrm{FD}}$ are  constant error factors caused by the Hadamard test and the finite difference method, respectively, and we used the following identity, 
$\tfrac{1}{\Delta t^{N}}\sum_{l=0}^{N}\binom{N}{l}=\left(\frac{2}{\Delta t}\right)^{N}$. Note that the upper bound is set for the worst-case scenario, although we expect that the error cancellation can occur among the additive errors having different signs. Using the upper bound of the inverse factorial, $\frac{1}{N!} < \frac{1}{\sqrt{2\pi N}}\left(\frac{e}{N}\right)^{N}$, we found an upper bound of the additive error expressed with powers,    
\begin{align}
\epsilon_{\mathrm{Tot}}\le\frac{1}{\sqrt{2\pi N}}\left(\frac{2e}{N\Delta t}\right)^{N}
\left(\epsilon_{\mathrm{HT}}+\frac{\epsilon_{\mathrm{FD}}N}{6}\left(\frac{\Vert\mathcal{H}(B)\Vert\Delta t}{2}\right)^{N+2}\right).
\end{align}
Because of the convergence condition, we have 
\begin{align}
   \epsilon_{\mathrm{Tot}}\le\epsilon_{\mathrm{R}}\left(\frac{2e}{N\Delta t}\right)^{N},
\end{align}
where $\epsilon_{\mathrm{R}}=1/\mathrm{poly}(N)$. 
The total error of the complex matrix case is bounded as   
\begin{align}
\epsilon'_{\mathrm{Tot}}&\le\frac{1}{N!}\left(\epsilon_{\mathrm{HT}}\left(\frac{4}{\Delta t}\right)^{N} +\epsilon_{\mathrm{FD}}\frac{N\Delta t^{2}}{24}\Vert\mathcal{H}(A)\Vert^{N-1}(\Vert \mathcal{H}(B)\Vert^{3}+\Vert \mathcal{H}(C)\Vert^{3})\right).
\end{align}
We obtain a similar upper bound as the real matrix case using the identity,   $\tfrac{1}{\Delta t^{N}}\sum_{l=0}^{N}\sum_{j=0}^{N-l}
    \sum_{k=0}^{l}\binom{N}{l}\binom{N-l}{j}\binom{l}{k}=\left(\frac{4}{\Delta t}\right)^{N}$, as  
\begin{align}
\epsilon'_{\mathrm{Tot}}\le\frac{1}{\sqrt{2\pi N}}\left(\frac{4e}{N\Delta t}\right)^{N}
\left(\epsilon_{\mathrm{HT}}+\frac{4\epsilon_{\mathrm{FD}}N}{3}\left(\frac{\Vert\mathcal{H}(A)\Vert\Delta t}{4}\right)^{N+2}\right),
\end{align}
where we used $\Vert\mathcal{H}(B)\Vert^{3}+\Vert\mathcal{H}(C)\Vert^{3}\le 2\Vert\mathcal{H}(A)\Vert^{3}$.
Again, due to the convergence condition, we have 
\begin{align}
   \epsilon'_{\mathrm{Tot}}\le\epsilon_{\mathrm{C}}\left(\frac{4e}{N\Delta t}\right)^{N},
   \label{eq:upperboundcomplex}
\end{align}
where $\epsilon_{\mathrm{C}}=1/\mathrm{poly}(N)$. 
    
\textbf{(i) Exponentially small additive error estimation:}
Using the upper bounds of the total errors, we can achieve an exponentially small additive error with the following conditions: 
\begin{align}
&\frac{2e}{N}<\Delta t\le\frac{2}{\Vert \mathcal{H}(B) \Vert}, \\
&\frac{4e}{N}<\Delta t\le \frac{4}{\Vert \mathcal{H}(A) \Vert}, 
\label{eq:expbound}
\end{align}
for real and complex matrices, respectively. Note that $\Vert \mathcal{H}(A) \Vert=\sum_{jk}\vert A_{jk}\vert$. 
We may consider a bound condition for the 1-norm of the matrix,  
\begin{align}
&\Vert \mathcal{H}(B) \Vert \le N\Vert B\Vert_{1}\le \frac{N}{e}, \\
&\Vert \mathcal{H}(A) \Vert \le N\Vert A\Vert_{1}\le \frac{N}{e}, 
\end{align}
therefore, $\Vert B\Vert_{1}\le 1/e$ and $\Vert A\Vert_{1}\le 1/e$ for real and complex matrices, respectively.  Because $1/e\ge\Vert A\Vert_{1}\ge\Vert A\Vert/\sqrt{N}$, the sufficient condition of the matrix 2-norm for exponentially small error is $\Vert A\Vert \le \sqrt{N}/e$.














\textbf{(ii) Additive error bound condition for beating the Gurvits classical sampling algorithm:}
Gurvits' classical sampling algorithm can estimate a matrix permanent with an additive error $\epsilon \Vert A \Vert^{N}$, where $0<\epsilon=1/\mathrm{poly}(N)<1$. It can produce exponentially small, polynomially small, or exponentially large errors depending on the size of the matrix spectral norm. Here, we test our quantum algorithm against Gurvits' classical sampling algorithm, considering what condition can yield a better additive error. Without loss of generality, we consider only complex matrices here because the real matrices carry identical bound conditions as the complex matrices.    

We devise an inequality using Eq.~\eqref{eq:upperboundcomplex} and Gurvits' additive error to find a condition for the superiority of our quantum algorithm, i.e. 
\begin{align}
   \epsilon\left(\frac{4e}{N\Delta t}\right)^{N}\le\epsilon \Vert A \Vert^{N},
\end{align}
where we set $0<\epsilon_{\mathrm{C}}=\epsilon<1$. The resulting bound for $\Delta t$ is given as follows: 
\begin{align}
    \frac{4e}{N\Vert A\Vert}\le \Delta t \le \frac{4}{\Vert\mathcal{H}(A)\Vert}.
    \label{eq:deltatbound}
\end{align}
Thus, we have a sufficient condition for the superiority of our quantum algorithm with respect to the Gurvits algorithm as  \begin{align}
    \Vert\mathcal{H}(A)\Vert\le\frac{N}{e}\Vert A \Vert.
    \label{eq:betterbound}
\end{align}
Based on this condition~\eqref{eq:betterbound}, we may consider three cases:  
\begin{enumerate}
    \item $\Vert\mathcal{H}(A)\Vert\le\frac{N}{e}\Vert A \Vert\le\frac{N}{e}$, which implies $\Vert A \Vert \le 1$: We can obtain exponentially small errors, smaller than Gurvits' exponentially small error. 
    \item $\Vert\mathcal{H}(A)\Vert\le\frac{N}{e}\le\frac{N}{e}\Vert A \Vert$, which implies $\Vert A \Vert \ge 1$: We can still obtain exponentially small errors, but the Gurvits algorithm produces an   exponentially large  error.
     \item $\frac{N}{e}\le\Vert\mathcal{H}(A)\Vert\le\frac{N}{e}\Vert A \Vert$, which implies $\Vert A \Vert \ge 1$: Both of our quantum algorithm and Gurvits' classical algorithm produce an  exponentially large error, but our quantum algorithm's error is smaller than the classical one.  
\end{enumerate}
The quantum algorithm can perform exponentially better than its classical counterpart if the condition in Eq.~\eqref{eq:betterbound} is satisfied. Here, we emphasise that case 2 above shows a domain where the quantum algorithm produces an exponentially small error,  whereas the classical algorithm carries an exponentially large error. 
Now, let us see what Eq.~\eqref{eq:betterbound} means.  
We know the Hamiltonian norm is upper-bounded by the 1-norm as
\begin{align}
    \Vert\mathcal{H}(A)\Vert\le N\Vert A \Vert_{1} \le\frac{N}{e}\Vert A \Vert,
\end{align}
where we set the upper bound, as in Eq.~\eqref{eq:betterbound}, to find a sufficient condition of matrix norms. 
Recalling the known inequality, $1/\sqrt{N}\le\Vert A\Vert/\Vert A\Vert_{1}\le\sqrt{N}$, we have
\begin{align}
    e\le\frac{\Vert A\Vert}{\Vert A\Vert_{1}}\le \sqrt{N}, 
    \label{eq:ratiobound}
\end{align}
which implies $e^{2}\le N$. We observe that the domain in Eq.~\eqref{eq:ratiobound} is non-empty for the sufficiently large $N\ge 8$. Let us call it the quantum advantage domain. As obviously shown in the bound condition in Eq.~\eqref{eq:ratiobound}, the quantum advantage domain becomes dominant as the matrix size increases and becomes larger than 50~\% for $N\ge 20$; see  Fig.~\ref{fig:quantumadvantage}.  
\begin{figure*}[t]
    \centering
    \includegraphics[width=0.7\textwidth]{quantumadvantageratio.jpg}
    \caption{Quantum advantage domain ratio ($Q(N)$) in Eq.~\eqref{eq:ratiobound}.}
    \label{fig:quantumadvantage}
\end{figure*}

\textbf{(iii) Average hardness for Gaussian matrices:}
Now, we discuss the average hardness of additive error estimation of Gaussian matrices. According to AA's two conjectures, the additive error estimation of a Gaussian matrix within $\epsilon \sqrt{N!}$ is a \#P-hard problem. If we can achieve this error with our quantum algorithm, then $\mathrm{BQP=P^{\#P}}$, which would be surprising again. 
We test the following inequality for this purpose: 
\begin{align}
   \epsilon\left(\frac{4e}{N\Delta t}\right)^{N}\le\epsilon \left(N!\right)^{k}<\epsilon \left(\frac{N}{2}\right)^{Nk},
\end{align}
where $k$ is an order parameter that $k\le 1/2$ would lead us to the surprising conclusion. 
Using the upper bound condition of $\Delta t$ in Eq.~\eqref{eq:deltatbound}, we have 
\begin{align}
    \frac{e}{2}\Vert\mathcal{H}(A)\Vert<\left(\frac{N}{2}\right)^{k+1}.  
\end{align}
Here, we introduce the average value of the Hamiltonian norm over random Gaussian matrices, $E(\Vert \mathcal{H}(A)\Vert)=\sqrt{2\pi^{-1}}N^{2}$, to the inequality, that is, 
\begin{align}
2e\sqrt{\frac{2}{\pi}}
<\left(\frac{N}{2}\right)^{k-1}.  
\end{align}
Therefore, $k$ has to be greater than 1. Our additive error quantum protocol based on the Hadamard test cannot reach the requirement for the \#P-hardness, unlike our multiplicative error protocol in the previous subsection.    

























    























\textbf{(iv) Quantum circuit synthesis:} 
\begin{table*}[tb]
\centering
\begin{tabular}{|m{5cm}||m{3cm}|m{5cm}|} 
 \hline
 \backslashbox[5cm]{Criteria}{Matrices} 
 & Arbitrary real ($B$)& Arbitrary complex ($A$) \\ 
 \hline  \hline
 Number of $\mathrm{Re}\braket{\phi\vert U\vert \phi}$ to be evaluated &
  $\floor{\frac{N+1}{2}}$& 
 Odd $N$: 
 $(N^{3}+6N^{2}+11N+6)/12$ 
 Even $N$: 
 $(N^{3}+6N^{2}+8N)/12$\\ 
 \hline
 Required number of qubits per $\mathrm{Re}\braket{\phi\vert U\vert \phi}$ estimation& $2N+1$& $2N+1$\\ 
 \hline
 Number of CNOT gates per $\mathrm{Re}\braket{\phi\vert U\vert \phi}$ estimation & $4N^{2}+2N$&$4N^{2}+2N$ \\ 
 \hline
 Circuit depth per $\mathrm{Re}\braket{\phi\vert U\vert \phi}$ estimation& $9N^2+1$&$9N^2+1$  \\ 
 \hline
 Total number of samples for estimating a permanent with the Hadamard test& $\mathcal{O}\left(\frac{N}{\epsilon_{\mathrm{HT}}^{2}}\log(\frac{1}{\delta})\right)$ & $\mathcal{O}\left(\frac{N^{3}}{\epsilon_{\mathrm{HT}}^{2}}\log(\frac{1}{\delta})\right)$\\ 
 \hline
 Total error upper bound & $\epsilon\left(\frac{2e}{N\Delta t}\right)^{N}$
 & $\epsilon\left(\frac{4e}{N\Delta t}\right)^{N}$\\ 
 \hline
 $\Delta t$ bound for efficient exponentially small additive error estimation& $\frac{2e}{N}<\Delta t < \frac{2}{\Vert \mathcal{H}(B) \Vert}$& $\frac{4e}{N}<\Delta t < \frac{4}{\Vert \mathcal{H}(A) \Vert}$\\ 
 \hline
 Sufficient condition for beating the Gurvits algorithm & $e\le \frac{\Vert B\Vert}{\Vert B\Vert_{1}}\le \sqrt{N}$ & $e\le \frac{\Vert A\Vert}{\Vert A\Vert_{1}}\le \sqrt{N}$\\ [\smallskipamount]
 \hline
\end{tabular}
\caption{Summary table for the additive error estimation protocol.}
\label{table:summary}
\end{table*}
As mentioned previously, the multiplicative error quantum protocol is efficient from a quantum complexity perspective, but its implementation is impractical. Therefore, here, we present the quantum circuit synthesis for the additive error protocol based on the Hadamard test for the practical application of our quantum algorithm in NISQ devices. Without loss of generality, we present the quantum circuit for complex matrices; for the real matrix case, we can simply set $C=0$. 
The time-propagator 
can be factorised into two-qubit gates as 
\begin{align}
    U(A^{(l,j,k)};\Delta t/2)=
    \prod_{p,q=1}^{N} \exp\left(-\mathrm{i}A_{pq}^{(l,j,k)}\frac{\Delta t}{2}Z_{N+p}Z_{q}\right), 
\end{align}
where the two-qubit gate is a simple Ising $ZZ$-gate $R_{ZZ}^{(p,q)}(\theta_{pq})=\exp(-\mathrm{i}\tfrac{\theta_{pq}}{2} Z_{N+p}Z_{q})$. 
We estimate the overlap integrals using the Hadamard test method in Fig.~\ref{fig:quantumcircuits}(a), and the required controlled-$U$ gate, $V$, is decomposed as a simple product of controlled-$R_{ZZ}$ gates as   
\begin{align}
&V(A^{(l,j,k)};\Delta t/2)=\prod_{p,q=1}^{N}
(\ket{0^{(c)}}\bra{0^{(c)}}+\ket{1^{(c)}}\bra{1^{(c)}}\otimes R_{ZZ}^{(p,q)}(A_{pq}^{(l,j,k)}\Delta t)), 
\end{align}
where $(c)$ denotes the control ancilla qubit. Accordingly, Fig.~\ref{fig:quantumcircuits}(a) shows the detailed quantum circuit corresponding to the Hadamard test. Thus, we need to evaluate  $\mathcal{O}(N)$ and $\mathcal{O}(N^{3})$ real part of the overlap integrals to estimate the permanents of real and complex matrices, respectively. In any case, it requires $2N+1$ qubits with $\mathcal{O}(N^{2})$ circuit depth and CNOT gates per single overlap integral evaluation. Tab.~\ref{table:summary} summarises the additive error protocol.  




















\section{Discussion and Conclusion}
This paper reported a new classical formula for matrix permanent named the Glynn-Kan formula, which was discovered by introducing the monomial expansion of Kan~\cite{kan:2007} to the Glynn formula~\cite{Glynn2010}. Although the Glynn-Kan formula is exponentially more inefficient than Rayser's formula, it has a simple structure of the $N$-th power of binary functions. Interestingly, for real matrices, it can be cast as a difference between two positive semidefinite functions that Stockmeyer's counting algorithm would lead it to the  $\mathrm{BPP^{NP}}$ algorithm. The implication would be surprising, i.e. $\mathrm{BPP^{NP}}=\mathrm{P^{\#P}}$, that the PH would collapse to the third level according to Toda's theorem~\cite{Aaronson2011}. However, this paper cannot conclude because the two multiplicative error factors of the positive semidefinite functions should be identical and efficiently obtainable to become a $\mathrm{BPP^{NP}}$ algorithm. It would be worth exploring to the study of the computational complexity of the permanents of various matrices using the Glynn-Kan formula. For example, the multiplicative error estimations of a matrix having 0-1 entries and the Hermitian positive semidefinite matrix are known to be obtainable in polynomial times.  The Glynn-Kan formula would help develop the $\mathrm{BPP^{NP}}$ algorithm for the Hermitian positive semidefinite matrix. It is worth mentioning that the Glynn-Kan formula can be generalized for matrices with repeated rows and columns~\cite{aaronson2012,Chin2018,Yung2019} and can be used to develop the corresponding quantum algorithms.   

By converting the classical Glynn-Kan formula to the quantum expectation value of the Glynn-Kan operator, we developed quantum algorithms for computing matrix permanents for the first time. Using the finite difference method~\cite{Seki2021}, we then approximated the quantum expectation value as the sum of the quantum overlap integrals. Consequently, the matrix permanent evaluation problem becomes a quantum overlap integral evaluation problem. We suggested two quantum protocols for the quantum overlap integral evaluation. A multiplicative error protocol was proposed by combining  KOS'~\cite{Knill2007} and AR's~\cite{Aaronson2020} amplitude estimation algorithms. It requires a polynomial number of queries to the Grover diffusion operator, invoking a polynomial number of quantum gates that a BQP quantum circuit can be used for the estimation. Accordingly, quantum computers can efficiently compute the matrix permanents with multiplicative errors. 

Because the multiplicative error estimation of matrix permanent is known to be a \#P-hard problem, the implication of our finding is again remarkable, $\mathrm{BQP}=\mathrm{P^{\#P}}$. Although that is strongly believed not to be the case~\cite{Aaronson2011}, the possibility is still open. So far, it has been shown by Lloyd et al. ~\cite{Abrams1998} that it requires nonlinearity for a quantum device to access a \#P-hard problem efficiently. Our multiplicative error quantum algorithm provides additional evidence to Refs.~\cite{Aaronson2010,Fefferman2010,Raz2019} that the BQP has capabilities beyond the PH. This paper proposes to show how we can adequately locate the BQP concerning the PH using quantum algorithms for computing matrix permanents. 
Another implication of our multiplicative error quantum algorithm is that if $\mathrm{BPP}=\mathrm{BQP}$, then according to our quantum algorithm, $\mathrm{BPP}=\mathrm{P^{\#P}}$ and the PH collapses. This would imply that quantum computing is more powerful than classical computing ($\mathrm{BPP}\neq\mathrm{BQP}$; see Theorem 1 and Corollary 1). However, we emphasise that we do not definitively claim that  $\mathrm{BQP}=\mathrm{P^{\#P}}$ in this paper; instead, we provide a new tool to investigate computational complexity questions about the relationships among the BPP, BQP, and \#P. More careful computational complexity analyses are required, including quantum device errors, not only with the statistical sampling errors. It would be interesting to test our multiplicative error protocol with the QFT-based algorithm by Brassard et al.~\cite{Brassard2002} if it can perform the multiplicative error amplitude estimation. Moreover, a more efficient quantum algorithm can be developed with quantum signal processing and LCU techniques~\cite{Childs2012,Berry2015,Low2019}. We believe that our new findings in classical and quantum algorithms for computing matrix permanents would play an essential role in identifying the unclear relationship between the BQP and PH. 
It is known so far only that $\mathrm{BQP \subseteq P^{\#P}=P^{PP}}$~\cite{Bernstein1997}. 
Hence, thorough computational complexity analyses are anticipated in subsequent studies. 

We also proposed an additive error quantum protocol for computing matrix permanents. Compared to the multiplicative error protocol, this algorithm can be practically implementable in NISQ devices. It requires $2N+1$ qubits and $4N^2+2N$ CNOT gates per single overlap integral evaluation of its real part. Interestingly, the required circuit complexities of the real and complex matrices are not different; only the required numbers of overlap integral evaluations are different. The additive error quantum protocol is comparable with the Gurvits classical sampling algorithm. We found a quantum advantage domain of the matrix norm where the quantum algorithm can provide a better estimation than the classical one. When the system size is greater than 20, we expect the quantum advantage domain to be more than 50 \%.              

Lastly, we note that our quantum algorithm directly connects the boson sampling and the IQP to approximate matrix permanents as a linear combination of the real-time partition function of quantum Ising Hamiltonian~\cite{Matsuo2014,Fujii2017,Mann2019}. The circuit complexity differs depending on the quantum Ising Hamiltonian's  matrix elements since different matrices would give different circuit complexities. Therefore, the quantum circuit model can analyse the computational complexity of a matrix permanent. Further investigations of the connection between boson sampling and qubit sampling (IQP) should be explored. 

Our quantum algorithm for computing matrix permanent attained a milestone in the quantum computing and the computational complexity theory and opened new computational complexity questions. It could be used to answer those questions by itself.  
Our quantum algorithm can be modified to compute other matrix functions like (loop)-hafnians, Torontonians, determinants, and immanants with minor changes; the related works will appear elsewhere soon. The subsequent paper will explore various numerical procedures, including the Gaussian process and interpolation techniques, applying our additive error quantum algorithm using cloud quantum computers like IBM Quantum and IonQ. 
With the help of our findings, we believe that the BQP would serve as an essential bridge in answering the Millennium question $\mathrm{P}=\mathrm{NP}$ by adequately positioning BPP, BQP, and \#P complexity classes.   










       


    
 






\begin{acknowledgments}
This work was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF), funded by the Ministry of Education, Science and Technology (NRF-2021M3E4A1038308, NRF-2021M3H3A1038085, NRF-2019M3E4A1079666). 
\end{acknowledgments}


\section*{References}
