Quantum estimation bound of Gaussian matrix permanent

Joonsuk Huh joonsukhuh@gmail.com Department of Chemistry, Sungkyunkwan University, Suwon 16419, Korea SKKU Advanced Institute of Nanotechnology (SAINT), Sungkyunkwan University, Suwon 16419, Korea Institute of Quantum Biophysics, Sungkyunkwan University, Suwon 16419, Korea

Abstract

Exact calculation and even multiplicative error estimation of matrix permanent are challenging for both classical and quantum computers. Regarding the permanents of random Gaussian matrices, the additive error estimation is closely linked to boson sampling, and achieving multiplicative error estimation requires exponentially many samplings. Our newly developed formula for matrix permanents and its corresponding quantum expression have enabled better estimation of the average additive error for random Gaussian matrices compared to Gurvits’ classical sampling algorithm. The well-known Ryser formula has been converted into a quantum permanent estimator. When dealing with real random Gaussian square matrices of size $N$ , the quantum estimator can approximate the matrix permanent with an additive error smaller than $\epsilon(\sqrt{\mathrm{e}N})^{N}$ , where $\epsilon$ is the estimation precision. In contrast, Gurvits’ classical sampling algorithm has an estimation error of $\epsilon(2\sqrt{N})^{N}$ , which is exponentially larger ( $1.2^{N}$ ) than the quantum method. As expected, the quantum additive error bound fails to reach the multiplicative error bound of $(2\pi N)^{1/4}\epsilon(\sqrt{N/\mathrm{e}})^{N}$ . Additionally, the quantum permanent estimator can be up to quadratically faster than the classical estimator when using quantum phase estimation-based amplitude estimation.

^†^†preprint: APS/123-QED

I Introduction

Almost three decades ago, Peter Shor [1] discovered a celebrated polynomial-time quantum algorithm for integer factorization. Although the best-known classical algorithm for integer factorization scales quasi-exponentially, the problem’s computational complexity still needs to be settled. Thus, claiming quantum supremacy with Shor’s quantum algorithm is unclear [2]. Moreover, implementing Shor’s quantum algorithm for factoring large integer numbers takes a lot of work.

In this context, Aaronson and Arkhipov [3] proposed a quantum sampling problem known as boson sampling, which has the potential to challenge the extended Church-Turing thesis and demonstrate quantum supremacy. The difficulty of boson sampling is based on the computational complexity of matrix permanents. Approximating permanents of arbitrary complex or real square matrices within a constant factor is a #P-hard problem. And boson sampling is a scattering process of non-interacting single photons in a linear optical network [4, 5]. The output photon distribution of the boson sampling device is expressed as a function of permanents of random Gaussian matrices (RGMs). Aaronson and Arkhipov demonstrated that if an efficient classical simulation of boson sampling, specifically an additive error matrix permanent estimator, were possible, it would lead to the collapse of the polynomial hierarchy (PH) to the third level, which is widely believed not to be the case. Therefore, the classical hardness of boson sampling relies on the complexity of calculating the permanents of complex RGMs. Still, the #P-hardness of permanents of Gaussian matrices has yet to be proven that the permanent of Gaussian conjecture (PGC) asserts that the multiplicative error estimation of permanents of RGMs is #P-hard. Additionally, assuming a non-uniform distribution of permanents of RGMs by the permanent anti-concentration conjecture (PAC), the multiplicative error estimation is essentially equivalent to the additive error estimation within an error proportional to their average magnitude.

Troyansky and Tishby [4] demonstrated a connection between a Hermitian matrix permanent and the expectation value of a linear optical quantum observable. They found that computing the permanent requires exponentially many samples due to significant observable variance. Therefore, the boson sampling device is impractical in estimating permanents, moreover, the magnitudes are generally exponentially small. Although there have been limited attempts to link quantum circuit models to matrix permanents, most studies have focused on proving the computational hardness of computing permanents using quantum circuit models [6, 7, 8, 9], rather than developing quantum algorithms for estimating matrix permanents.

Even though there is no proof, it is believed that quantum computers cannot solve the #P-hard problem in polynomial time. It is important to note that estimating the matrix permanent, in this paper, involves additive error approximation, and not exact calculation or relative error estimation, which are widely believed to be impossible even for quantum computers. This paper aims to explore whether quantum computers can estimate matrix permanents and how their performance compares to classical sampling algorithms.

We have developed a new formula to computing an $N$ -dimensional matrix permanent based on Ryser’s formula [10]. Our formula expresses the permanent as a sum of $N$ -th powers of a binary function, which allows us to reformulate the power expansion as an overlap integral using the quantum Ising Hamiltonian propagator. We have proposed additive error protocols using the Hadamard test [11, 12] and quantum phase estimation (QPE) [13]. We found a recently published paper [14], which connects matrix permanents to spin models as Fefferman et al. [8] and Peropadre et al. [15]. We compared our results of real Gaussian matrices using the method of Fefferman et al. [8] to the result of Gurvits’ classical sampling algorithm [16, 3], and found an improved upper bound of additive error estimation of Gaussian matrix permanents compared to the existing classical sampling algorithm. However, it cannot achieve the expected relative error limit, as this would lead to the collapse of PH to BQP (bounded-error quantum polynomial time). Additionally, we note that the quantum estimator can achieve a quadratic speedup when using QPE-based amplitude estimation.

The paper’s structure is as follows: First, we transform the Ryser formula into a sum of $N$ -th power of binary functions. Then, we demonstrate that a matrix permanent can be approximately represented as a quantum overlap integral, followed by the additive error analysis with real RGMs. Finally, we discuss the implications of our findings and conclude the paper.

II Results

II.1 Classical formulas

Recall the definition of permanent of a square matrix $A\in\mathbb{C}^{N\times N}$ :

\displaystyle\mathrm{Per}(A)=\sum_{\vec{{\sigma}}\in S_{\mathrm{p}}}\prod_{j=1}^{N}A_{j,\sigma_{j}},

(1)

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. (1) requires $N!$ terms in the summation, and each term is composed of products of $N$ elements of the matrix. In 1963, Ryser [10] discovered a formula for matrix permanent, which involves only $2^{N}$ terms in the summation,

\displaystyle\textrm{Per}(A)=(-1)^{N}\sum_{S\subseteq\{1,\ldots,N\}}(-1)^{|S|}\prod_{j=1}^{N}\sum_{k\in S}A_{jk},

(2)

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. We can rewrite the Ryser formula (2) by introducing an $N$ -dimensional binary vector $\vec{y}$ as follows:

\displaystyle\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}),

(3)

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 [17, 18, 19, 20]:

\displaystyle\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),

(4)

where we used the identity $\sum_{\vec{x}\in\{-1,1\}^{N}}\prod_{k}^{N}x_{k}=0$ . Eq. 4 shows the magnitude is bounded by $\|A\|^{N}$ . Gurvits [16, 3] suggested an additive error algorithm $\pm\epsilon\|A\|^{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}$ . $\|\cdot\|$ is the matrix 2-norm; otherwise, the order is indicated.

Using the identity of a monomial by Lemma 1 of Kan [21], we have

\displaystyle\prod_{j=1}^{N}\vec{A}_{j}\cdot\vec{x}=\frac{1}{N!2^{N}}\sum_{\vec{x}^{\prime}\in\{-1,1\}^{N}}\left(\prod_{k=1}^{N}x^{\prime}_{k}\right)\left(\vec{x^{\prime}}^{\top}A\vec{x}\right)^{N}.

(5)

Now, we have a new expression for the matrix permanent by inserting Eq. (5) into Eq. (4) as follows:

\displaystyle\mathrm{Per}(A)=\frac{1}{N!2^{2N}}\sum_{\vec{x},\vec{x}^{\prime}\in\{-1,1\}^{N}}\left(\prod_{k=1}^{N}x_{k}x_{k}^{\prime}\right)\left(\vec{x^{\prime}}^{\top}A\vec{x}\right)^{N},

(6)

whose magnitude is bounded by $\|A\|_{1,1}^{N}/N!$ and $\|A\|_{1,1}=\sum_{j,k}|A_{jk}|$ .

We name our new finding in Eq. (6) the Glynn-Kan-Huh (GKH) formula in this paper. Thus, Per( $A$ ) can be interpreted as an average of $2^{2N}$ $N$ -th power of binary functions by Eq. (6). The implication of the GKH formula is discussed in Appendix A.

At a glance, our new formula for permanents in Eq. (6) is exponentially more inefficient compared with the original Ryser formula, Eq. (2). However, our new formula is suitable for developing quantum algorithms for permanents because they are composed of an expectation value of $N$ -th powers of the binary function, i.e. $\vec{x^{\prime}}^{\top}A\vec{x}$ . In the following section, we develop a quantum algorithm for additive error approximation of matrix permanents.

II.2 Quantum estimator

From the GKH formula in Eq. (6), we immediately identify the following quantum expectation value as a matrix permanent:

	$\displaystyle\textrm{Per}(A)=$
	$\displaystyle\tfrac{1}{N!2^{2N}}\sum_{y=0}^{2^{2N}-1}\langle y\|\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}\|y\rangle,$		(7)

where $|y\rangle=|y_{1}\rangle\cdots|y_{2N}\rangle$ , $y=y_{1}\cdots y_{2N}$ is given as a $2N$ -bit string and $Z_{k}$ is a Pauli-Z operator of the $k$ -th qubit. Using the identity, $\frac{1}{2^{N}}\sum_{y=0}^{2^{2N}-1}|y\rangle=\bigotimes_{k=1}^{2N}H_{k}|\mathbf{0}\rangle=|\phi\rangle$ , where $H_{k}$ is a Hadamard operator of the $k$ -th qubit and $|\mathbf{0}\rangle=|0\rangle^{\otimes 2N}$ , we rewrite the matrix permanent as an expectation value of a quantum operator, i.e., $\mathrm{Per}(A)=\frac{1}{N!2^{2N}}\mathrm{Tr}[G_{Z}(A)]=\frac{1}{N!}\langle\phi|G_{Z}(A)|\phi\rangle=\frac{1}{N!}\langle\mathbf{0}|G_{X}(A)|\mathbf{0}\rangle$ , where $G_{X}(A)=(\bigotimes_{k=1}^{2N}H_{k})G_{Z}(A)(\bigotimes_{k=1}^{2N}H_{k})$ and $X_{k}$ = $H_{k}Z_{k}H_{k}$ . Let us call $G_{X}(A)$ the GKH operator, defined as

\displaystyle G_{X}(A)=\left(\bigotimes_{k=1}^{N}X_{k}X_{N+k}\right)\left(\sum_{j,k=1}^{N}A_{jk}X_{N+j}X_{k}\right)^{N},

(8)

which is factored with two components–parity operator and Ising Hamiltonian operator, respectively denoted as $P$ and $\mathcal{H}$ . That is, $G_{X}(A)=P\mathcal{H}^{N}(A)$ , where $\mathcal{H}(A)=\sum_{j,k=1}^{N}A_{jk}X_{N+j}X_{k}$ and $P=\bigotimes_{k=1}^{N}X_{k}X_{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 non-unitary and generally non-Hermitian. The remaining task is to approximately implement the power operator in a quantum circuit. We can perform the approximation through various approaches [22]. For example, we can consider quantum signal processing with quantum singular value transformation and the linear combination of unitaries (LCU) technique [23, 24, 25]. However, the repeated application of LCU for the power operator requires a prefactor for the normalization of the quantum states grows as the power of the Hamiltonian norm, $\|\mathcal{H}(A)\|^{N}$ ( $\|\mathcal{H}(A)\|\leq\|A\|_{1,1}$ ), resulting in large errors in the end. The quantum power method of Seki et al. [26] can also be used. The method is based on the finite difference to approximate the Hamiltonian power operator as the high-order time differentiation of the time propagator. However, it also suffers from a large error due to a convergence condition related to the Hamiltonian norm. Instead, we adopt the method by Fefferman et al. [8], where the hardness of the exact boson sampling was discussed with permanents of integer matrices because it produces less error than the quantum power method. The matrix permanent is approximated as a single overlap integral of the time-propagator of quantum Ising Hamiltonian. The overlap integral is evaluated within an additive error using the Hadamard test or the QPE-based amplitude estimation.

II.2.1 Matrix permanent and time-propagator

Let’s rewrite the matrix permanent with the GKH operator and leave only the Hamiltonian power operator as follows,

	$\displaystyle\mathrm{Per}(A)$	$\displaystyle=\frac{1}{N!}\langle\mathbf{0}\|P\mathcal{H}(A)^{N}\|\mathbf{0}\rangle$
		$\displaystyle=\langle\mathbf{1}\|\frac{\mathcal{H}(A)^{N}}{N!}\|\mathbf{0}\rangle,$		(9)

where $P|\mathbf{0}\rangle=|\mathbf{1}\rangle=|1\cdots 1\rangle$ . Eq. (9) can be interpreted as a transition amplitude from $|\mathbf{0}\rangle$ to $|\mathbf{1}\rangle$ . A single application of $\mathcal{H}(A)$ results in the pairwise bit flips with amplitudes $A_{jk}$ , and thus, the matrix permanent is a counting amplitudes problem. On the other hand, we can consider an overlap integral of the time-propagator, i.e.,

\displaystyle\langle\mathbf{1}|\exp(-\mathrm{i}\mathcal{H}(A)t)|\mathbf{0}\rangle=\langle\mathbf{1}|\sum_{l=0}^{\infty}\frac{(-\mathrm{i}t\mathcal{H}(A))^{N+2l}}{(N+2l)!}|\mathbf{0}\rangle,

(10)

where the terms of Taylor series in a lower order than $N$ vanish because the scattering operations of $\mathcal{H}(A)$ on $|\mathbf{0}\rangle$ at least $N$ times to form $|\mathbf{1}\rangle$ and the pairwise selection should be made an even number of times for the amplitude to survive that the additional odd order (from the $N$ -th order) contributions vanish [8]. Now, we can relate the overlap integral of the time-propagator with the matrix permanent as follows:

\displaystyle\langle\mathbf{1}|\exp(-\mathrm{i}\mathcal{H}(A)t)|\mathbf{0}\rangle

\displaystyle=(-\mathrm{i}t)^{N}\mathrm{Per}(A)+\sum_{l=1}^{\infty}(-\mathrm{i}t)^{N+2l}R_{l},

(11)

where the overlap integral is expanded as a time-series that the leading order’s coefficient is the matrix permanent, and the remaining error terms follow: $R_{l}=\langle\mathbf{1}|\frac{\mathcal{H}(A)^{N+2l}}{(N+2l)!}|\mathbf{0}\rangle$ . Here, we assume $A$ is a real matrix to have the time-propagator as a unitary operator, and $N$ is an even number for brevity. So only the real part matters. For odd $N$ , we need to take care of the imaginary part instead of the real part. After rearranging Eq. (11), we can express the matrix permanent as follows:

\displaystyle\mathrm{Per}(A)=(-t^{2})^{-N/2}\left[\langle\mathbf{1}|\exp(-\mathrm{i}\mathcal{H}(A)t)|\mathbf{0}\rangle+R\right],

(12)

where $R=-\sum_{l=1}^{\infty}(-t^{2})^{N/2+l}R_{l}$ and the equation implies if the remaining error term is small enough, we can approximate the matrix permanent via estimating the overlap integral of the time-propagator. Note that we only need the real part of the overlap integral because only the real part survives. The quantum estimation of the overlap integral can be made by the Hadamard test [11, 27, 28] or the QPE-based method [13], whose costs scale $\mathcal{O}(1/\epsilon^{{}^{\prime}2})$ and $\mathcal{O}(1/\epsilon^{\prime})$ , respectively. We note here that the computational cost of the Hadamard test of the Ising Hamiltonian scales identically to the classical samplings because the Hamiltonian is provided in a computational basis. Thus, the computational advantage can be achieved only in the QPE approach, which scales quadratically better.

Now the matrix permanent estimation problem becomes an evaluation problem of the real part of the quantum overlap integral, i.e., $\mathrm{Re}\braket{\phi}{U}{\phi}=\mathrm{Re}\braket{\mathbf{0}}{(\bigotimes_{k=1}^{2N}H_{k})U(\bigotimes_{k=1}^{2N}H_{k})}{\mathbf{0}}$ , for a unitary operator $U$ , which includes the parity operator and the time-propagator. We present a Hadamard test quantum circuit for the overlap integral evaluation in Appendix B. The Hadamard test circuit requires a polynomial number of quantum gates and depth. We refer readers to Knill et al. [13] for the QPE-based algorithm.

II.2.2 Average additive error upper bound

This paper focuses on the real RGMs to discuss the algorithm’s remarkable implications. The average upper bound of $|R_{l}|$ in Eq. (12) (i.e., $\mathcal{E}(|R_{l}|$ ) with respect to real RGMs $\{W\}$ , where $W\sim\mathcal{N}(0,1)^{N\times N}_{\mathbb{R}}$ of i.i.d. real Gaussians is

$\displaystyle\mathcal{E}(\|R_{l}\|)$	$\displaystyle\leq\frac{1}{2}\mathcal{E}\left(\langle\mathbf{0}\|\frac{\mathcal{H}(W)^{N+2l}}{(N+2l)!}\|\mathbf{0}\rangle+\langle\mathbf{1}\|\frac{\mathcal{H}(W)^{N+2l}}{(N+2l)!}\|\mathbf{1}\rangle\right)$
	$\displaystyle=\frac{1}{2(N+2l)!}\frac{\partial^{N+2l}}{\partial\beta^{N+2l}}\mathcal{E}(\langle\mathbf{0}\|\exp(\beta\mathcal{H}(W))\|\mathbf{0}\rangle$
	$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\langle\mathbf{1}\|\exp(\beta\mathcal{H}(W))\|\mathbf{1}\rangle)\Big{\|}_{\beta=0},$	(13)

where we used the fact that $\mathcal{H}^{N+2l}$ is a positive semidefinite operator (by assuming $N$ as an even number) for the inequality and we introduced a generating function with a parameter $\beta$ to perform the Gaussian integral for averaging $\mathcal{E}$ with respect to the Gaussian distribution. The resulting identity is a simple exponential function of $\beta^{2}$ :

$\displaystyle\mathcal{E}\left(\langle\mathbf{0}\|\exp(\beta\mathcal{H}(W))\|\mathbf{0}\rangle\right)$	$\displaystyle=\langle\mathbf{0}\|\mathcal{E}\left(\mathrm{e}^{\beta\sum_{j,k=1}^{N}W_{jk}X_{N+j}X_{k}}\right)\|\mathbf{0}\rangle$
	$\displaystyle=\exp\left(\frac{\beta^{2}N^{2}}{2}\right)$
	$\displaystyle=\mathcal{E}\left(\langle\mathbf{1}\|\exp(\beta\mathcal{H}(W))\|\mathbf{1}\rangle\right),$	(14)

where, in the second line, the one-dimensional Gaussian integrals were performed regarding the corresponding normal distributions of $\{W_{jk}\}$ , and we used $X_{N+j}^{2}X_{k}^{2}=1$ . An alternative proof of the above relation is given in Appendix C. Finally, we have an upper bound of $\mathcal{E}(|R_{l}|)$ :

	$\displaystyle\mathcal{E}(\|R_{l}\|)$	$\displaystyle\leq\frac{1}{(N+2l)!}\frac{\partial^{N+2l}}{\partial\beta^{N+2l}}\exp\left(\frac{\beta^{2}N^{2}}{2}\right)\Big{\|}_{\beta=0}$
		$\displaystyle=\frac{(N^{2}/2)^{N/2+l}}{(N/2+l)!}.$		(15)

Consequently, the remaining error term in Eq. (12) has an averaged upper bound,

$\displaystyle\mathcal{E}\left(\left\|R\right\|\right)$	$\displaystyle\leq(Nt^{2})^{N/2}\sum_{l=1}^{\infty}\frac{(N/2)^{N/2+l}(Nt^{2})^{l}}{(N/2+l)!}$
	$\displaystyle\leq(Nt^{2})^{N/2}\left(\sum_{l=1}^{\infty}\frac{(N/2)^{N/2+l}}{(N/2+l)!}\right)\left(\sum_{l=1}^{\infty}(Nt^{2})^{l}\right)$
	$\displaystyle<(Nt^{2})^{N/2}\exp(N/2)\left(\sum_{l=1}^{\infty}(Nt^{2})^{l}\right)$
	$\displaystyle=\frac{Nt^{2}}{(1-Nt^{2})}(\mathrm{e}Nt^{2})^{N/2},$	(16)

where we assumed $Nt^{2}<1$ for the convergence and used

\displaystyle\left(\sum_{l=1}^{\infty}\frac{(N/2)^{N/2+l}}{(N/2+l)!}\right)=\exp(N/2)-\sum_{l=1}^{N/2}\frac{(N/2)^{l}}{l!}.

(17)

To make the averaged remaining error bounded by an exponentially small number ( $\epsilon_{\mathrm{R}}=(\mathrm{e}-c)^{-1}c^{N/2+1}$ ), we need $c=\mathrm{e}Nt^{2}<1$ , which is tighter than $Nt^{2}<1$ .

The quantum amplitude is estimated using the Hadamard test or the QPE-based amplitude estimation algorithm with the additive error precision $\epsilon^{\prime}$ . Therefore, the total error $\epsilon_{\mathrm{T}}$ in Eq. (12) has an average upper bound as follows:

\displaystyle\epsilon_{\mathrm{T}}

\displaystyle\leq t^{-N}\left(\epsilon^{\prime}+\epsilon_{\mathrm{R}}\right),

(18)

where we set $\epsilon=\epsilon^{\prime}+\epsilon_{\mathrm{R}}$ and the upper bound becomes $\epsilon(\sqrt{\mathrm{e}N/c})^{N}$ when $t=\sqrt{c/(N\mathrm{e})}$ . Note that $\epsilon^{\prime}$ is polynomially small and $\epsilon_{\mathrm{R}}$ is exponentially small. Here, we also note that the quantum algorithm’s sampling cost can scale $\mathcal{O}(1/\epsilon^{\prime})$ , which is quadratically better than the classical sampler if we use the QPE-based amplitude estimation [13] instead of the Hadamard test. Now, we compare our error bound with the Gurvits classical sampling algorithm [16, 3] for real RGMs.

Gurvits algorithm [16, 3]– The Gurvits classical sampling algorithm generates an additive error $\epsilon\|W\|^{N}$ with a sampling cost $\mathcal{O}(1/\epsilon^{2})$ and the average error is $\epsilon(2\sqrt{N})^{N}$ , where we used $\mathcal{E}(\|W\|)=2\sqrt{N}$ [3, 29, 30]. This error is an exponential factor $(2/\sqrt{\mathrm{e}/c})^{N}<1.2^{N}$ (where $\mathrm{e}/4<c<1$ ) larger than the quantum version.

II.2.3 Tighter average error bound using GKH formula

The upper bound in Eq. (16) can be improved by taking the alternating sign changes in the summation into account as follows:

	$\displaystyle\|R\|$	$\displaystyle=\left\|-\sum_{k=0}^{\infty}\left[(-t^{2})^{N/2+2k+1}R_{2k+1}+(-t^{2})^{N/2+2k+2}R_{2k+2}\right]\right\|$
		$\displaystyle=\left\|-(-t^{2})^{N/2}\sum_{k=0}^{\infty}(-t^{2})^{2k+1}(R_{2k+1}-t^{2}R_{2k+2})\right\|.$		(19)

Note that $R_{l}\geq 0$ ( $R_{l}<0$ ) for all $l$ when $\mathrm{Per}(A)\geq 0$ ( $\mathrm{Per}(A)<0$ ); this is clearly seen when the structure of $R_{l}$ in eigenvalues is considered. Without losing the generality, we assume $\mathrm{Per}(A)\geq 0$ and thus $R_{l}\geq 0$ for all $l$ . Similar to the GKH formula (6), we can express $R_{l}$ as follows:

\displaystyle R_{l}=\frac{1}{(N+2l)!2^{2N}}\sum_{\vec{x},\vec{x}^{\prime}\in\{-1,1\}^{N}}\left(\prod_{k=1}^{N}x_{k}x_{k}^{\prime}\right)\left(\vec{x^{\prime}}^{\top}A\vec{x}\right)^{N+2l}.

(20)

Accordingly, $R_{l}$ can be rewritten as a difference between two positive semidefinite terms, i.e.,

	$\displaystyle R_{l}$	$\displaystyle=\sum_{\{\vec{x},\vec{x^{\prime}}\}\in\{-1,1\}^{N}}\frac{\mathcal{S}_{+}\left(\vec{x^{\prime}}^{\top}A\vec{x}\right)^{N+2l}}{(N+2l)!2^{2N}}-\frac{\mathcal{S}_{-}\left(\vec{x^{\prime}}^{\top}A\vec{x}\right)^{N+2l}}{(N+2l)!2^{2N}}$
		$\displaystyle=R_{l}^{(+)}-R_{l}^{(-)},$		(21)

where the first and second terms are denoted $R_{l}^{(\pm)}$ and $\mathcal{S}_{\pm}$ refer to $\left(\prod_{k=1}^{N}x_{k}\pm\prod_{k=1}^{N}x_{k}^{\prime}\right)^{2}/4$ , respectively. With the new notation, we can recast $R$ as follows:

	$\displaystyle\|R\|=\Big{\|}-(-t^{2})^{N/2}\sum_{k=0}^{\infty}(-t^{2})^{2k+1}$	$\displaystyle[R_{2k+1}^{(+)}-t^{2}R_{2k+2}^{(+)}$
	$\displaystyle-($	$\displaystyle R_{2k+1}^{(-)}-t^{2}R_{2k+2}^{(-)})]\Big{\|}.$		(22)

When $t\leq N$ , $R_{2k+1}^{(\pm)}\geq t^{2}R_{2k+2}^{(\pm)}$ because $0\leq R_{2k+1}^{(\pm)}\leq R_{2k+2}^{(\pm)}(N+4k+4)(N+4k+3)$ , which can be shown by the power mean inequality:

\displaystyle R_{l+1}^{\pm}(N+2l+2)!\geq[R_{l}^{\pm}(N+2l)!]^{\frac{N+2l+2}{N+2l}}\geq R_{l}^{\pm}(N+2l)!.

(23)

The last inequality is generally true for random Gaussian matrices because $\|A\|_{1,1}=\mathcal{O}(N^{2})$ that $R_{l}^{\pm}(N+2l)!\geq 1$ . Then, we conclude as follows:

\displaystyle R_{l+1}^{\pm}(N+2l+2)(N+2l+1)\geq R_{l}^{\pm}.

(24)

Thus, with a small $t$ ( $\leq N$ ), $R$ is bounded as follows:

$\displaystyle\|R\|$	$\displaystyle\leq(t^{2})^{N/2}\sum_{k=0}^{\infty}(t^{2})^{2k+1}[R_{2k+1}^{(+)}-t^{2}R_{2k+2}^{(+)}$
	$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+(R_{2k+1}^{(-)}-t^{2}R_{2k+2}^{(-)})]$
	$\displaystyle=-(t^{2})^{N/2}\sum_{l=1}^{\infty}(-t^{2})^{l}(R_{l}^{(+)}+R_{l}^{(-)}).$	(25)

Now we can set the average upper bound of $R$ , i.e.,

$\displaystyle\mathcal{E}(\|R\|)$	$\displaystyle\leq-(t^{2})^{N/2}\sum_{l=1}^{\infty}(-t^{2})^{l}\mathcal{E}(R_{l}^{(+)}+R_{l}^{(-)})$
	$\displaystyle=-(-1)^{N/2}\sum_{l=1}^{\infty}\frac{(-N^{2}t^{2}/2)^{N/2+l}}{(N/2+l)!}$
	$\displaystyle=\left\|\exp(-N^{2}t^{2}/2)-\sum_{l=0}^{N/2}\frac{(-N^{2}t^{2}/2)^{l}}{l!}\right\|$
	$\displaystyle=\frac{\exp(-\theta N^{2}t^{2}/2)}{(N/2+1)!}\left(\frac{N^{2}t^{2}}{2}\right)^{N/2+1},$	(26)

where $0\leq\theta<1$ in the Lagrange form of the remainder in Taylor’s theorem for the Gaussian function. We used the Stirling approximation,

\displaystyle\frac{1}{(N/2+1)!}<\frac{1}{\sqrt{2\pi(N/2+1)}}\left(\frac{\mathrm{e}}{N/2+1}\right)^{N/2+1},

(27)

for the following inequality,

	$\displaystyle\frac{\mathrm{e}^{-\theta N^{2}t^{2}/2}\left(\frac{N^{2}t^{2}}{2}\right)^{N/2+1}}{(N/2+1)!}$	$\displaystyle<\frac{\mathrm{e}^{-\theta N^{2}t^{2}/2}\left(\frac{\mathrm{e}N^{2}t^{2}}{N+2}\right)^{N/2+1}}{\sqrt{2\pi(N/2+1)}}$
		$\displaystyle=\frac{\left[\exp\left(1-\theta\tau\right)\tau\right]^{N/2+1}}{\sqrt{2\pi(N/2+1)}},$		(28)

where $\tau=N^{2}t^{2}/(N+2)$ . To make the average upper bound of $|R|$ exponentially small, we set $\exp\left(1-\theta\tau\right)\tau<1$ that $\frac{\mathrm{ln}(\mathrm{e}\tau)}{\tau}<\theta$ . Therefore, when $\tau=1/\mathrm{e}$ , the constant $\theta$ always exists in $(0,1)$ . Finally, we have an exponentially small upper bound,

\displaystyle\mathcal{E}(R)<\frac{\mathrm{e}^{-\tfrac{\theta}{\mathrm{e}}}}{\sqrt{2\pi(N/2+1)}}\exp\left(-\tfrac{\theta}{2\mathrm{e}}N\right)=\epsilon_{\mathrm{R}}^{\mathrm{tight}},

(29)

with

\displaystyle t=\left(\sqrt{\mathrm{e}\left(\frac{N}{N+2}\right)N}\right)^{-1},

(30)

and

\displaystyle t^{-N}=\left(\sqrt{\mathrm{e}\left(\frac{N}{N+2}\right)N}\right)^{N}.

(31)

Consequently, we have the total error bound as

$\displaystyle\epsilon_{\mathrm{T}}$	$\displaystyle\leq t^{-N}\left(\epsilon^{\prime}+\epsilon_{\mathrm{R}}^{\mathrm{tight}}\right)$
	$\displaystyle=\epsilon\left(\sqrt{\mathrm{e}\left(\frac{N}{N+2}\right)N}\right)^{N}$
	$\displaystyle<\epsilon\left(\sqrt{\mathrm{e}N}\right)^{N},$	(32)

where we set $\epsilon=\epsilon^{\prime}+\epsilon_{\mathrm{R}}^{\mathrm{tight}}$ .

III Conclusion and Outlook

We have introduced a new classical formula for matrix permanent called the GKH (Glynn-Kan-Huh) formula. This formula was derived by combining the monomial expansion of Kan [21] with the Glynn formula [17]. Recently, it was rederived using quantum optical language [31]. While the GKH formula is more inefficient than Ryser’s formula, it has a simple structure of the $N$ -th power of a binary function, allowing for the potential use of matrix symmetry and linear algebra techniques, such as singular value decomposition, to reduce computational costs. Unlike the Ryser and Glynn formulas, linear algebra can be introduced to the GKH formula. For instance, it is known that the multiplicative error estimations of a 0-1 entry matrix and a Hermitian positive semidefinite matrix can be obtained in polynomial time. The GKH formula could facilitate the development of the BPP^NP algorithm for the Hermitian positive semidefinite matrix. It is also noteworthy that the GKH formula can be extended to matrices with repeated rows and columns [18, 19, 20, 31] and can be used to develop corresponding quantum formulas.

The classical GKH formula can be converted to the quantum expectation value of the GKH operator, allowing for the estimation of matrix permanents using quantum amplitudes. Quantum algorithms based on the GKH operator can take different forms depending on how the Hamiltonian power operator is approximated. The Hamiltonian power operator can be incorporated into a large unitary operator, which may include ancilla qubits, such as LCU or quantum singular value transformation. In either case, the permanent can be expressed as a single overlap integral (denoted as $a$ ) of the large unitary operator. The resulting permanent expression contains a significant prefactor, i.e., $\mathrm{Per}(A)=\|A\|^{N}_{1,1}/N!\times a$ , resulting in generally superexponentially large (for $\|A\|_{1,1}\geq\mathcal{O}(N^{2})$ ) additive error in the additive error estimation protocol (e.g., Hadamard test) for the overlap integral or a multiplicative error with exponentially expensive cost due to the exponentially small amplitude [13]. However, it is unknown whether some matrices give polynomially small amplitudes for efficient multiplicative error quantum amplitude estimation [13]. When $\|A\|_{1,1}=\mathcal{O}(N)$ and the size of permanent is only exponentially large, the amplitude can be in a polynomial size. If such matrices exist, the relative error permanent estimation of such matrices may not be #P-hard, because we strongly believe BQP cannot access the #P-hardness [3]; the possibility is still open, but, so far, it was shown by Lloyd et al. [32] that it requires nonlinearity for a quantum device to access a #P-hard problem efficiently.

Instead, in this paper, we approximated the matrix permanent by representing it as an overlap integral of a time-propagator of the quantum Ising Hamiltonian. The first surviving term in the Taylor expansion with respect to time $t$ is proportional to the matrix permanent, while the higher-order terms are considered as an additive error. In this paper, we examine the real RGMs and discuss the implications of our findings. By imposing a condition on the size of the time parameter, $t=\mathcal{O}(1/\sqrt{N})$ , we showed that the averaged remaining error can be exponentially small. As a result, we can approximate the permanents of real RGMs with an average additive error smaller than $\epsilon(\sqrt{\mathrm{e}N})^{N}$ , which represents an exponential factor ( $1.2^{N}$ ) improvement compared to the Gurvits classical sampler. It is important to note that our upper bound might be more optimal, and there could be potential for improving the quantum estimator’s error. However, we caution that achieving the relative error limit ( $(2\pi N)^{1/4}\epsilon(\sqrt{N/\mathrm{e}})^{N}$ ) would potentially lead the Polynomial Hierarchy to collapse to BQP, due to conjectures (PGA and PAC) [3].

Interestingly, the complexity of the circuit depends on the matrix elements of the quantum Ising Hamiltonian, as different matrices result in various circuit complexities [27, 33, 34]. Thus, the quantum circuit model can be used to analyze the computational complexity of a matrix permanent (cf. [9]).

Although the estimation error of the quantum estimator did not meet the relative error limit of $(2\pi N)^{1/4}\epsilon(\sqrt{N/\mathrm{e}})^{N}$ , our quantum analysis represents a significant breakthrough in the quantum estimation of matrix functions such as (loop)-hafnians, Torontonians, determinants, and immanants with minor modifications. More details about this work will be available soon.

Acknowledgements.

The author acknowledges Raymond Kan, Hyukjoon Kwon, Dominik Hangleiter, Gwonhak Lee, Chae-Yeun Park, and Taewan Kim for the fruitful discussions. 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, NRF-2022M3H3A106307411, NRF-2023M3K5A1094805, NRF-2023M3K5A1094813). This work was also partly supported by Institute for Information & communications Technology Promotion (IITP) grant funded by the Korea government(MSIP) (No. 2019-0-00003, Research and Development of Core technologies for Programming, Running, Implementing and Validating of Fault-Tolerant Quantum Computing System).

Appendix

Appendix A Classical implication of of the GKH formula (Eq. (6))

Interestingly, Eq. (6) can be rewritten as a difference between two positive semidefinite functions for a real square matrix of $A$ with an even dimension,

	$\displaystyle\mathrm{Per}(A)=$	$\displaystyle\frac{1}{N!2^{2N}}\Huge[\sum_{\{\vec{x},\vec{x^{\prime}}\}\in\{-1,1\}^{N}}\mathcal{S}_{+}\left(\vec{x^{\prime}}^{\top}A\vec{x}\right)^{N}$
		$\displaystyle-\sum_{\{\vec{x},\vec{x^{\prime}}\}\in\{-1,1\}^{N}}\mathcal{S}_{-}\left(\vec{x^{\prime}}^{\top}A\vec{x}\right)^{N}\Huge],$		(33)

where $\mathcal{S}_{\pm}$ refers to $\left(\prod_{k=1}^{N}x_{k}\pm\prod_{k=1}^{N}x_{k}^{\prime}\right)^{2}/4$ , respectively. For the odd-dimensional case, we can use the following identity to have the even-dimensional matrix, $\mathrm{Per}(A)=\mathrm{Per}\left(A\bigoplus 1\right)$ . Eq. (33) directly shows #P-hard is GapP-hard by the definition of GapP [35]. On the other hand, the multiplicative error estimation of a permanent of real matrix is also #P-hard. The positive semidefinite functions in Eq. (33) can be evaluated with multiplicative error estimation using Stockmeyer’s BPP^NP counting algorithm [3, 36], where BPP and NP stand for bounded error probabilistic polynomial time and nondeterministic polynomial time, respectively.

Meanwhile, the multiplicative error estimation of the permanent of a Hermitian positive semidefinite matrix is classified as $\mathrm{BPP^{NP}}$ [37, 38, 39]. Although it is known that permanents of non-negative matrices can be evaluated within a constant factor by Jerrum-Sinclair-Vigoda (JSV) algorithm [40] in polynomial time, the JSV algorithm has no practical implementations so far [41, 42]. Therefore, further computational complexity analysis of estimating permanents of real matrices using Eq. (33) is required. We would be able to discover new efficient classical algorithms with the GKH formula.

On the other hand, we have a condition for the superiority of the GKH sampler (Monte Carlo sampling with the GKH formula) with respect to the Gurvits algorithm as

\displaystyle\epsilon\frac{\|A\|_{1,1}^{N}}{N!}\leq\epsilon\|A\|^{N},

(34)

which implies

\displaystyle\|A\|_{1,1}\leq\frac{N}{\mathrm{e}}\|A\|,

(35)

where we used $N!\simeq(N/\mathrm{e})^{N}$ .

Based on this condition (35), we may consider three cases:

1.

$\|A\|_{1,1}\leq\frac{N}{e}\|A\|\leq\frac{N}{e}$ , which implies $\|A\|\leq 1$ : The GKH sampler produces exponentially small errors, smaller than Gurvits’ exponentially small error.
2.

$\|A\|_{1,1}\leq\frac{N}{e}\leq\frac{N}{e}\|A\|$ , which implies $\|A\|\geq 1$ : The GKH sampler still produces exponentially small errors, but the Gurvits algorithm makes an exponentially large error.
3.

$\frac{N}{e}\leq\|A\|_{1,1}\leq\frac{N}{e}\|A\|$ , which implies $\|A\|\geq 1$ : Both classical samplers produce exponentially large errors, but the GKH sampler’s error is smaller than the Gurvits’ one.

We can set a tighter condition with the 1-norm and 2-norm relation as

\displaystyle\|A\|_{1,1}\leq N\|A\|_{1}\leq\frac{N}{e}\|A\|,

(36)

where we set the upper bound, as in Eq. (35), to find a sufficient condition of matrix norms. Recalling the well-known inequality, $1/\sqrt{N}\leq\|A\|/\|A\|_{1}\leq\sqrt{N}$ , we have

\displaystyle e\leq\frac{\|A\|}{\|A\|_{1}}\leq\sqrt{N},

(37)

which implies $e^{2}\leq N$ . We observe that the domain in Eq. (37) is non-empty for the sufficiently large $N\geq 8$ and even it is dominant (50 %) when $N\geq 28$ . However, we note, here, that the random Gaussian matrices do not belong to the domain in average because the average 1-norm and 2-norm are $\mathcal{O}(N)$ and $\mathcal{O}(\sqrt{N})$ , respectively. Therefore, the Gurvits sampler performs better than the GKH sampler for random Gaussian matrices.

Appendix B Circuit synthesis for Hadamard test

Refer to caption — Figure 1: Quantum circuits for evaluating a real part of the overlap integral using the Hadamard test. (a) The Hadamard test circuit for the additive error estimation of the real part of the quantum overlap integral $\mathrm{Re}\braket{\phi}{U}{\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. 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)$ , where $R_{Z}(\theta)=\exp(-\mathrm{i}\tfrac{\theta}{2}Z)$ .

Here, we present the quantum circuit synthesis for the additive error protocol based on the Hadamard test.

The time-propagator can be factorized into two-qubit gates as

\displaystyle U(A^{\prime};t/2)=(-1)^{N}\prod_{p,q=1}^{N}\exp\left(-\mathrm{i}A^{\prime}_{pq}\frac{t}{2}Z_{N+p}Z_{q}\right),

(38)

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. 1(a), and the required controlled- $U$ gate, $V$ , is decomposed as a simple product of controlled- $R_{ZZ}$ gates as

	$\displaystyle V(A^{\prime};\Delta t/2)=$
	$\displaystyle\prod_{p,q=1}^{N}(\ket{0^{(c)}}\bra{0^{(c)}}+\ket{1^{(c)}}\bra{1^{(c)}}\otimes R_{ZZ}^{(p,q)}(A^{\prime}_{pq}t)),$		(39)

where $(c)$ denotes the control ancilla qubit. Accordingly, Fig. 1 shows the detailed quantum circuit corresponding to the Hadamard test, where $A^{\prime}=A+\pi t^{-1}I$ and $I$ is an identity matrix. We used the following identity for the parity operator, $(\bigotimes_{k=1}^{2N}H_{k})P(\bigotimes_{k=1}^{2N}H_{k})=P_{Z}=(-1)^{N}\bigotimes_{k=1}^{N}\exp\left(-\frac{\pi}{2}\mathrm{i}Z_{k}Z_{N+k}\right)$ , to integrate it in the unitary propagator $U$ of the quantum Ising Hamiltonian. After all, it requires $2N+1$ qubits with $\mathcal{O}(N^{2})$ circuit depth and CNOT gates per single overlap integral evaluation and $\mathcal{O}(1/\epsilon^{{}^{\prime}2})$ samples.

Appendix C Alternative proof of Eq. (14)

	$\displaystyle\mathcal{E}\left(\langle\mathbf{0}\|\exp(\beta\mathcal{H}(W))\|\mathbf{0}\rangle\right)$
	$\displaystyle=\mathcal{E}\left(\langle\mathbf{0}\|\exp(\beta\sum_{j,k=1}^{N}W_{jk}X_{N+j}X_{k})\|\mathbf{0}\rangle\right)$
	$\displaystyle=\mathcal{E}\left(\langle\phi\|\exp(\beta\sum_{j,k=1}^{N}W_{jk}Z_{N+j}Z_{k})\|\phi\rangle\right)$
	$\displaystyle=\frac{1}{2^{2N}}\sum_{\vec{x},\vec{x}^{\prime}\in\{-1,1\}^{N}}\mathcal{E}\left(\exp(\beta\sum_{j,k=1}^{N}W_{jk}x^{\prime}_{j}x_{k})\right)$
	$\displaystyle=\frac{1}{2^{2N}}\sum_{\vec{x},\vec{x}^{\prime}\in\{-1,1\}^{N}}\prod_{j,k=1}^{N}\left[\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\mathrm{d}W_{jk}\mathrm{e}^{-\frac{1}{2}W_{jk}^{2}+\beta W_{jk}x^{\prime}_{j}x_{k}}\right]$
	$\displaystyle=\exp\left(\frac{\beta^{2}N^{2}}{2}\right),$		(40)

where we used $x_{j}^{{}^{\prime}2}x_{k}^{2}=1$ .

References

Shor [1999] P. W. Shor, Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer, SIAM Rev. 41, 303 (1999).
Harrow and Montanaro [2017] A. W. Harrow and A. Montanaro, Quantum computational supremacy, Nature 549, 203 (2017).
Aaronson and Arkhipov [2013] S. Aaronson and A. Arkhipov, The computational complexity of linear optics, Theory Comput. 9, 143 (2013).
Troyansky and Tishby [1996] L. Troyansky and N. Tishby, Permanent Uncertainty : On the quantum evaluation of the determinant and the permanent of a matrix, in Proc. 4th Workshop on Physics and Computation PhysComp96, 1 (1996).
Scheel [2004] S. Scheel, Permanents in linear optical networks, arXiv:quant-ph/0406127 (2004).
Rudolph [2009] T. Rudolph, Simple encoding of a quantum circuit amplitude as a matrix permanent, Phys. Rev. A 80, 4 (2009).
Aaronson [2011] S. Aaronson, A linear-optical proof that the permanent is #P-hard, Proc. R. Soc. A: Math. Phys. Eng. Sci. 467, 3393 (2011).
Fefferman et al. [2017] B. Fefferman, M. Foss-Feig, and A. V. Gorshkov, Exact sampling hardness of Ising spin models, Phys. Rev. A 96, 032324 (2017).
Park et al. [2023] C.-Y. Park, P. A. M. Casares, J. M. Arrazola, and J. Huh, The hardness of quantum spin dynamics, arXiv:quant-ph/2312.07658 (2023).
Ryser [2014] H. J. Ryser, Combinatorial Mathematics, the Carus mathematical monographs (The Mathematical Association of America, 2014).
Aharonov et al. [2007] D. Aharonov, I. Arad, E. Eban, and Z. Landau, Polynomial Quantum Algorithms for Additive approximations of the Potts model and other Points of the Tutte Plane, arXiv:quant-ph/0702008 (2007).
Aharonov et al. [2009] D. Aharonov, V. Jones, and Z. Landau, A polynomial quantum algorithm for approximating the jones polynomial, Algorithmica (New York) 55, 395 (2009).
Knill et al. [2007] E. Knill, G. Ortiz, and R. D. Somma, Optimal quantum measurements of expectation values of observables, Phys. Rev. A 75, 012328 (2007).
Alase et al. [2023] A. Alase, O. Doty, and D. L. Feder, Matrix permanent and determinant from a spin system, Phys. Rev. A 108, 012207 (2023).
Peropadre et al. [2017] B. Peropadre, A. Aspuru-Guzik, and J. J. García-Ripoll, Equivalence between spin Hamiltonians and boson sampling, Phys. Rev. A 95, 032327 (2017).
Gurvits [2005] L. Gurvits, On the Complexity of Mixed Discriminants and Related Problems, In 30th Internat. Symp. Mathematical Foundations of Computer Science MFCS’05, 447 (2005).
Glynn [2010] D. G. Glynn, The permanent of a square matrix, Eur. J. Comb. 31, 1887 (2010).
Aaronson and Hance [2012] S. Aaronson and T. Hance, Generalizing and Derandomizing Gurvits’s Approximation Algorithm for the Permanent, arXiv:quant-ph/1212.0025 (2012).
Chin and Huh [2018] S. Chin and J. Huh, Generalized concurrence in boson sampling, Sci. Rep. 8, 1 (2018).
Yung et al. [2019] M.-H. Yung, X. Gao, and J. Huh, Universal bound on sampling bosons in linear optics and its computational implications, Nat. Sci. Rev. 6, 719 (2019).
Kan [2008] R. Kan, From moments of sum to moments of product, J. Multivar. Anal. 99, 542 (2008).
Aulicino et al. [2022] J. C. Aulicino, T. Keen, and B. Peng, State preparation and evolution in quantum computing: A perspective from Hamiltonian moments, Int. J. Quantum Chem. 122 (2022).
Childs and Wiebe [2012] A. M. Childs and N. Wiebe, Hamiltonian Simulation Using Linear Combinations of Unitary Operations, Quantum Inf. Comput. 12, 901 (2012).
Berry et al. [2015] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, Simulating hamiltonian dynamics with a truncated taylor series, Phys. Rev. Lett. 114, 1 (2015).
Low and Chuang [2019] G. H. Low and I. L. Chuang, Hamiltonian Simulation by Qubitization, Quantum 3, 163 (2019).
Seki and Yunoki [2021] K. Seki and S. Yunoki, Quantum Power Method by a Superposition of Time-Evolved States, PRX Quantum 2, 010333 (2021).
Matsuo et al. [2014] A. Matsuo, K. Fujii, and N. Imoto, Quantum algorithm for an additive approximation of Ising partition functions, Phys. Rev. A 90, 1 (2014).
Huang et al. [2019] H.-Y. Huang, K. Bharti, and P. Rebentrost, Near-term quantum algorithms for linear systems of equations, arXiv:quant-ph/1909.07344 (2019).
Vershynin [2018] R. Vershynin, High-Dimensional Probability: An Introduction with Applications in Data Science, Cambridge Series in Statistical and Probabilistic Mathematics (Cambridge University Press, 2018).
Rebrova and Vershynin [2018] E. Rebrova and R. Vershynin, Norms of random matrices: Local and global problems, Advances in Mathematics 324, 40 (2018).
Chabaud et al. [2022] U. Chabaud, A. Deshpande, and S. Mehraban, Quantum-inspired permanent identities, Quantum 6, 877 (2022).
Abrams and Lloyd [1998] D. S. Abrams and S. Lloyd, Nonlinear Quantum Mechanics Implies Polynomial-Time Solution for NP-Complete and #P Problems, Phys. Rev. Lett. 81, 3992 (1998).
Fujii and Morimae [2017] K. Fujii and T. Morimae, Commuting quantum circuits and complexity of Ising partition functions, New J. Phys. 19, 033003 (2017).
Mann and Bremner [2019] R. L. Mann and M. J. Bremner, Approximation algorithms for complex-valued ising models on bounded degree graphs, Quantum 3, 162 (2019).
Hangleiter and Eisert [2023] D. Hangleiter and J. Eisert, Computational advantage of quantum random sampling, Rev. Mod. Phys. 95, 035001 (2023).
Stockmeyer [1985] L. J. Stockmeyer, On approximation algorithms for #P, SIAM J. Comput. 14, 849 (1985).
Lund et al. [2014] A. P. Lund, A. Laing, S. Rahimi-Keshari, T. Rudolph, J. L. O’Brien, and T. C. Ralph, Boson Sampling from a Gaussian State, Phys. Rev. Lett. 113, 100502 (2014).
Rahimi-Keshari et al. [2015] S. Rahimi-Keshari, A. P. Lund, and T. C. Ralph, What Can Quantum Optics Say about Computational Complexity Theory?, Phys. Rev. Lett. 114, 060501 (2015).
Chakhmakhchyan and Cerf [2017] L. Chakhmakhchyan and N. J. Cerf, Quantum-inspired algorithm for estimating the permanent of positive semidefinite matrices, Phys. Rev. A 022329, 1 (2017).
Jerrum et al. [2004] M. Jerrum, A. Sinclair, and E. Vigoda, A polynomial-time approximation algorithm for the permanent of a matrix with nonnegative entries, J. ACM 51, 671 (2004).
Bezáková et al. [2008] I. Bezáková, D. Štefankovič, V. V. Vazirani, and E. Vigoda, Accelerating Simulated Annealing for the Permanent and Combinatorial Counting Problems, SIAM J. Comput. 37, 1429 (2008).
Newman and Vardi [2020] J. E. Newman and M. Y. Vardi, FPRAS Approximation of the Matrix Permanent in Practice, arXiv:cs/2012.03367 (2020).

$\displaystyle\mathcal{E}(\|R_{l}\|)$	$\displaystyle\leq\frac{1}{2}\mathcal{E}\left(\langle\mathbf{0}\|\frac{\mathcal{H}(W)^{N+2l}}{(N+2l)!}\|\mathbf{0}\rangle+\langle\mathbf{1}\|\frac{\mathcal{H}(W)^{N+2l}}{(N+2l)!}\|\mathbf{1}\rangle\right)$
	$\displaystyle=\frac{1}{2(N+2l)!}\frac{\partial^{N+2l}}{\partial\beta^{N+2l}}\mathcal{E}(\langle\mathbf{0}\|\exp(\beta\mathcal{H}(W))\|\mathbf{0}\rangle$
	$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\langle\mathbf{1}\|\exp(\beta\mathcal{H}(W))\|\mathbf{1}\rangle)\Big{\|}_{\beta=0},$	(13)

$\displaystyle\mathcal{E}\left(\langle\mathbf{0}\|\exp(\beta\mathcal{H}(W))\|\mathbf{0}\rangle\right)$	$\displaystyle=\langle\mathbf{0}\|\mathcal{E}\left(\mathrm{e}^{\beta\sum_{j,k=1}^{N}W_{jk}X_{N+j}X_{k}}\right)\|\mathbf{0}\rangle$
	$\displaystyle=\exp\left(\frac{\beta^{2}N^{2}}{2}\right)$
	$\displaystyle=\mathcal{E}\left(\langle\mathbf{1}\|\exp(\beta\mathcal{H}(W))\|\mathbf{1}\rangle\right),$	(14)

	$\displaystyle\|R\|$	$\displaystyle=\left\|-\sum_{k=0}^{\infty}\left[(-t^{2})^{N/2+2k+1}R_{2k+1}+(-t^{2})^{N/2+2k+2}R_{2k+2}\right]\right\|$
		$\displaystyle=\left\|-(-t^{2})^{N/2}\sum_{k=0}^{\infty}(-t^{2})^{2k+1}(R_{2k+1}-t^{2}R_{2k+2})\right\|.$		(19)