REDUCED DENSITY MATRIX ESTIMATION FOR PARTICLE-NUMBER-CONSERVING FERMION SYSTEM

Abstract

A computing device including one or more processing devices configured to receive classical shadow measurement results associated with a fermion wavefunction. The fermion wavefunction may describe a particle-number-conserving fermion system. The one or more processing devices may compute an estimated k-reduced density matrix (k-RDM), where k is a number of ladder operator pairs. Computing the estimated k-RDM may include computing, as terms of the estimated k-RDM, traces of respective products of a number operator of the particle-number-conserving fermion system, a single-particle basis rotation matrix dependent upon the classical shadow measurement result, elementary symmetric polynomials of Majorana bivectors, and a conjugate transpose of the single-particle basis rotation matrix. Computing each trace may include iteratively computing derivatives of a Pfaffian. The one or more processing devices may output the estimated k-RDM.

Description

BACKGROUND

Simulating the wavefunction of a system of fermions, such as in electronic structure problems in chemistry and material sciences, is one of the most promising applications of quantum computers. Once the fermion wavefunction is prepared on a quantum device, this wavefunction may be characterized in order to evaluate its properties, such as energy, density, and other many-body correlations.

In the case where the wavefunction ρ has particle number symmetry, meaning that the number of fermions η occupying n modes is fixed, the properties of the wavefunction are observables O that may be expressed as linear combinations of number-conserving k-Reduced Density Matrices (k-RDMs) D_{{right arrow over (q)}}^{{right arrow over (p)}}.

$\langle D_{\overrightarrow{q}}^{\overrightarrow{p}}\rangle =\mathrm{Tr}\left[D_{\overrightarrow{q}}^{\overrightarrow{p}}\rho \right]$

In the above equation, each k-RDM may be expressed as a product of k fermion creation operators and k fermion annihilation operators:

$D_{\overrightarrow{q}}^{\overrightarrow{p}}=a_{p_1}^{†}\dots a_{p_k}^{†}{a}_{q_k}\dots a_{q_1}$

The fermion creation and annihilation operators satisfy the standard anti-commutation relations:

$\left\{a_j,a_k^{†}\right\}={\delta }_{\mathrm{jk}}$ $\left\{a_j,a_k\right\}=\left\{a_j^{†},a_k^{†}\right\}=0$

In the above definition of D_{{right arrow over (q)}}^{{right arrow over (p)}}, the vectors {right arrow over (p)}, {right arrow over (q)} represent which modes are occupied by fermions. The integer elements of these vectors may be listed in ascending order, such that {right arrow over (p)}∈_n,η≐{(z₁, . . . , z_η): ∀j∈[η], 1≤z_j<z_j+1≤n}. In addition, the number operator for the fermion system is given by {circumflex over (n)}_j=a_j^†a_j.

Any observable then has the following form:

$O=\sum _{\overrightarrow{p},\overrightarrow{q}\in {𝒮}_{n,k}}{o}_{\overrightarrow{p},\overrightarrow{q}}{D}_{\overrightarrow{q}}^{\overrightarrow{p}}$

for some k and complex coefficients o_{{right arrow over (q)}}^{{right arrow over (p)}} For example, dipole moment and electron density are characterized by k=1, the electric structure Hamiltonian by k=2, the effective nuclear Hamiltonian by k=3, and the overlap with a single Slater determinant by k=η. Note that some observables may be linear combination of different k, but the equation for O is without loss of generality as the case of fixed particle number allows up-folding of any k-RDM to a sum of (k+1)-RDMs through the following identity:

$D_{\overrightarrow{q}}^{\overrightarrow{p}}=D_{\overrightarrow{q}}^{\overrightarrow{p}}\frac{{\sum }_{j\in \left[n\right]}{a}_j^{†}{a}_j}{\eta }$

Note that it suffices to consider at most k≤η, and RDMs with k>η are automatically zero.

SUMMARY

According to one aspect of the present disclosure, a computing device is provided, including one or more processing devices configured to receive a plurality of classical shadow measurement results associated with a fermion wavefunction. The fermion wavefunction may describe a plurality of fermions in a particle-number-conserving fermion system that occupy a plurality of modes. The one or more processing devices may be further configured to compute an estimated k-reduced density matrix (k-RDM), where k is a number of ladder operator pairs. Computing the estimated k-RDM may include computing, as terms of the estimated k-RDM, a plurality of traces of respective products of a number operator of the particle-number-conserving fermion system, a single-particle basis rotation matrix dependent upon the classical shadow measurement result, a plurality of elementary symmetric polynomials of a plurality of Majorana bivectors respectively associated with the plurality of fermions, and a conjugate transpose of the single-particle basis rotation matrix. Computing each of the traces may include iteratively computing a plurality of derivatives of a Pfaffian. The one or more processing devices may be further configured to output the estimated k-RDM.

This Summary is provided to introduce a selection of concepts in a simplified form that are further described below in the Detailed Description. This Summary is not intended to identify key features or essential features of the claimed subject matter, nor is it intended to be used to limit the scope of the claimed subject matter. Furthermore, the claimed subject matter is not limited to implementations that solve any or all disadvantages noted in any part of this disclosure.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 schematically shows a computing system at which one or more processing devices are configured to compute an estimated k-RDM according to a conventional k-RDM estimation protocol.

FIG. 2 schematically shows a computing system at which one or more processing devices are configured to compute an estimated k-RDM using an approach that has a lower time complexity than the conventional k-RDM estimation protocol.

FIG. 3 schematically shows the one or more processing devices when a plurality of terms of the estimated k-RDM are computed at least in part by computing a plurality of traces, according to the example of FIG. 2.

FIG. 4A shows a flowchart of a method for use with a computing device to compute an estimated k-RDM, according to the example of FIG. 2.

FIG. 4B shows additional steps of the method of FIG. 4A that may be performed for each of the classical shadow measurement results.

FIG. 4C shows additional steps of the method of FIG. 4A that may be performed when a plurality of traces and a plurality of derivates of Pfaffians are computed.

FIG. 5 shows a schematic view of an example computing system in which the computing device of FIG. 2 may be instantiated.

DETAILED DESCRIPTION

There exist prior protocols to evaluate all k-RDMs of a fermion system to some target error, or standard deviation, ϵ. These protocols include forming ρ, measuring it in some basis, and computing a noisy estimate of the k-RDMs from the measurement outcome. Such protocols may be made more efficient by minimizing the number of copies of ρ needed to achieve the target error ϵ. In the case where there is no particle symmetry, roughly

$\left(\begin{array}{c}n\\ k\end{array}\right)\sqrt{k}/{ϵ}^2$

measurements of the quantum state suffice. This expression for the sample efficiency has an explicit dependence on the number of modes but not the number of particles.

One prior protocol exploits particle number symmetry to estimate all number-conserving k-RDMs with a sample complexity of roughly

$\left(\begin{array}{c}\eta \\ k\end{array}\right)/{ϵ}^2$

that depends only on the number of fermions. This increase in sample complexity is substantial, as there are many natural applications where η«n, such as in simulating dynamical correlations or in a plane-wave representation of electronic structure. Although the estimator in this protocol is computationally efficient on a classical computer, it may still be impractical, requiring (n⁵) time to compute a single estimate of a single k-RDM. This (n⁵) time complexity may be a bottleneck for large-scale estimation considering that the electronic structure Hamiltonian already has (n⁴) RDMs in general. In addition, other applications such as unbiasing fermionic Monte Carlo may utilize millions of new estimates in every run.

FIG. 1 schematically shows a computing system 10 at which an estimated k-RDM {circumflex over (D)}_{{right arrow over (q)}}^{{right arrow over (p)}} is computed according to the prior protocol. The computing system 10 depicted in FIG. 1 includes a quantum computing device 12 that is configured to communicate with a classical computing device 20. The classical computing device 20 includes one or more processing devices 22 and memory 24.

As shown in the example of FIG. 1, the quantum computing device 12 receives a specification of a fermion wavefunction ρ over η fermions 14 in a particle-number-conserving fermion system 30 that occupy n modes 16. The fermions 14 may, for example, be electrons and atomic nuclei included in a molecule for which the fermion wavefunction ρ is simulated. The specification of the fermion wavefunction ρ may, for example, be received as an eigenstate of a Hamiltonian.

The one or more processing devices 22 are configured to compute the estimated k-RDMs {circumflex over (D)}_{{right arrow over (q)}}^{{right arrow over (p)}} as expectation values of a k-RDM operator {circumflex over (D)}_{{right arrow over (q)}}^{{right arrow over (p)}}. In the previous protocol, the expectation values {circumflex over (D)}_{{right arrow over (q)}}^{{right arrow over (p)}} of the k-RDM operators {circumflex over (D)}_{{right arrow over (q)}}^{{right arrow over (p)}} are estimated using classical shadow measurement results |{right arrow over (z)}. These classical shadow measurement results |{right arrow over (z)} may be measured as observables of a rotated fermion wavefunction obtained by applying a sampling rotation U(u_H) to the fermion wavefunction ρ. The sampling rotation U(u_H) is a rotation in a single-particle basis. Given a unitary matrix u_H sampled from the n×n Haar random distribution on the unitary group U_n, the sampling rotation U(u_H) may be defined as rotating each creation operator a_k^† to a linear combination of other creation operators as follows:

$U\left(u_H\right)a_k^{†}{U}^{†}\left(u_H\right)=\sum _j{u}_{kj}{a}_j^{†}$

The rotated fermion wavefunction may be given by UρU^†. The sampling rotations U(u_H) also have, within the η-particle subspace, matrix elements given by the following:

$\forall \overrightarrow{x},\overrightarrow{y}\in {𝒮}_{n,\eta },\langle \overrightarrow{x}❘U\left(u_H\right)❘\overrightarrow{y}\rangle ={\left({\wedge }^{\eta }{u}_H\right)}_{\overrightarrow{x},\overrightarrow{y}}=\det \left[u_{\overrightarrow{x},\overrightarrow{y}}\right]$

where ∧ is the exterior product and u_{{right arrow over (x)},{right arrow over (y)}} denotes the submatrix formed by taking rows x₁, x₂, . . . and columns y₁, y₂, . . . of u_H.

The quantum state of the rotated fermion wavefunction UρU^† may be measured in the computational basis to obtain a classical shadow measurement result |{right arrow over (z)}, as discussed above. The classical shadow measurement result |{right arrow over (z)} may be obtained with the following probability:

$\mathrm{Pr}[❘\overrightarrow{z}\rangle ]=\langle \overrightarrow{z}|U\left(u_H\right)\rho U^{†}\left(u_H\right)❘\overrightarrow{z}\rangle$

A permutation matrix that maps elements of {right arrow over (z)} to [η] may be defined as v_{{right arrow over (z)}}∈^n×n. Using this permutation matrix, a classical shadow may be computed as v=v_{{right arrow over (z)}}U_H.

Subsequently to computing the classical shadow v, an estimated k-RDM {circumflex over (D)}_{{right arrow over (q)}}^{{right arrow over (p)}} may be computed from the classical shadow v using the following formula:

$\langle {\overline{D}}_{\overrightarrow{q}}^{\overrightarrow{p}}\rangle =\langle \overrightarrow{q}❘U_k^{†}\left(v\right)E_{\eta ,k}{U}_k\left(v\right)❘\overrightarrow{p}\rangle$ $\mathrm{where}$ $E_{\eta ,k}\doteq \sum _{\overrightarrow{\tau }\in {𝒮}_{n,k}}❘\overrightarrow{r}\rangle \langle \overrightarrow{r}❘\frac{\left(\begin{array}{c}\eta -s^{\prime }\\ k-s^{\prime }\end{array}\right)\left(\begin{array}{c}\eta -\eta +s^{\prime }\\ s^{\prime }\end{array}\right)}{{\left(-1\right)}^{k+s^{\prime }}\left(\begin{array}{c}k\\ s^{\prime }\end{array}\right)}$

The above equation for k-RDM {circumflex over (D)}_{{right arrow over (q)}}^{{right arrow over (p)}} a involves matrix multiplication of dimension

$\left(\begin{array}{c}n\\ k\end{array}\right)\times \left(\begin{array}{c}n\\ k\end{array}\right)$

matrices and is therefore inefficient for large values of k.

The previous protocol also includes a method of computing {circumflex over (D)}_{{right arrow over (q)}}^{{right arrow over (p)}} or any k. Hermitian combinations of k-RDMs may be defined as follows:

$D_{\overrightarrow{q}}^{\overrightarrow{p}}=\frac{1}{2}\left(D_{\overrightarrow{q};0}^{\overrightarrow{p}}-i{D}_{\overrightarrow{q};\pi /2}^{\overrightarrow{p}}\right)$ $D_{\overrightarrow{q};\varphi }^{\overrightarrow{p}}\doteq e^{i\varphi }{D}_{\overrightarrow{q}}^{\overrightarrow{p}}+e^{-i\varphi }{D}_{\overrightarrow{p}}^{\overrightarrow{q}}$

Any Hermitian k-RDM D_{{right arrow over (q;ϕ)}}^{{right arrow over (p)}} is a linear combination of at most two diagonal k-RDMs. For example, when {right arrow over (p)}=[k] and {right arrow over (q)}=[2k]\[k]=(k+1 . . . 2k), then:

$D_{\left[2k\right]\setminus \left[k\right];\varphi }^{\left[k\right]}=U\left(\omega \left(\varphi \right)\right)\left({\hat{n}}_{\left[k\right]}-{\hat{n}}_{\left[2k\right]\setminus \left[k\right]}\right)U^{†}\left(\omega \left(\varphi \right)\right)$

where ω(ϕ) is a phase matrix defined as:

$\omega \left(\varphi \right)=\left(\frac{1}{\sqrt{2}}\left[\begin{array}{cc}{e}^{i\varphi /2}{I}_{\left[k\right]}& e^{i\varphi /2}{I}_{\left[k\right]}\\ e^{-i\varphi /2}{I}_{\left[k\right]}& e^{-i\varphi /2}{I}_{\left[k\right]}\end{array}\right]\right)\oplus I_{\left[n\right]\setminus \left[2k\right]}$

and where I_{{right arrow over (x)}}=Σ_{j∈{right arrow over (x)}}|j)j|. Accordingly, D_{{right arrow over (q;ϕ)}}^{{right arrow over (p)}} is a rotation of D_{{right arrow over (q)}}^{{right arrow over (p)}} in the complex plane.

Any rotated k-RDM D_{{right arrow over (q;ϕ)}}^{{right arrow over (p)}} may then be obtained by multiplying ω with a permutation that maps [k] to {right arrow over (p)} and [2k]\[k] to {right arrow over (q)}. The case where {right arrow over (p)} and {right arrow over (q)} have k′ elements not in common is treated analogously by the replacement ({circumflex over (n)}_[k]−{circumflex over (n)}_[2k]\[k])→({circumflex over (n)}_[k]−{circumflex over (n)}_[2k]\[k]){circumflex over (n)}_{[2k+k′]\[k′]}, since U(ω)a_j^†U(ω)^†=a_j^† acts as an identity on fermion operators with subscripts j>2k. In the above expressions, {circumflex over (n)}_[k]={circumflex over (n)}₁ . . . {circumflex over (n)}_k. Using this replacement, any k-RDM may be expressed as a linear combination of at most four diagonal estimated k-RDMs D_[k]^[k] conjugated by single-particle basis rotations U(ω_j),

$D_{\overrightarrow{q}}^{\overrightarrow{p}}=\sum _{j=1}^4{\alpha }_{j}U\left({\omega }_j\right){\hat{D}}_{\left[k\right]}^{\left[k\right]}{U}^{†}\left({\omega }_j\right)$

where the ω_j are permutations of the sparse matrix ω(ϕ).

In the previous protocol, the estimator for any k-RDM may be computed by combining the above equation for D_{{right arrow over (q)}}^{{right arrow over (p)}} with an alternate expression for {circumflex over (D)}_{{right arrow over (q)}}^{{right arrow over (p)}} given by the following formula:

$\langle U_k\left({\omega }_j\right){\hat{D}}_{\left[k\right]}^{\left[k\right]}{U}_k^{†}\left({\omega }_j\right)\rangle ={\mathrm{Tr}}_k\left[{\hat{n}}_{\left[k\right]}{U}_k\left({\omega }_j^{†}v\right)E_{\eta ,k}{U}_k^{†}\left({\omega }_j^{†}v\right)\right]$ $\mathrm{where}$ $E_{\eta ,k}\doteq \sum _{s=0}^k\frac{\left(\begin{array}{c}\eta -k+s\\ s\end{array}\right)\left(\begin{array}{c}n-\eta +k-s\\ k-s\end{array}\right)}{{\left(-1\right)}^{2k-s}\left(\begin{array}{c}k\\ s\end{array}\right)}{\mathrm{Sim}}_{k,s}\left({\prod }_{\left[\eta \right]}\right)$ $\mathrm{Si}{m}_{k,s}\left({\prod }_{\left[\eta \right]}\right)\doteq \sum _{j=0}^k{\left(-1\right)}^{j+s}\left(\begin{array}{c}j\\ s\end{array}\right)e_j\left({\hat{n}}_1\dots {\hat{n}}_{\eta }\right)={\left(-1\right)}^s\sum _{\alpha =0}^{\eta }{f}_{k,s}\left(\alpha \right)i^{\alpha }{e}_j\left({\beta }_1,\dots ,{\beta }_{\eta }\right)$ $f_{k,s}\left(\alpha \right)\doteq \sum _{j=\alpha }^k{\left(-1\right)}^j\left(\begin{array}{c}j\\ s\end{array}\right)\frac{1}{2^j}\left(\begin{array}{c}\eta -\alpha \\ \eta -j\end{array}\right)$

The trace may be computed over the k-particle subspace, such that:

$T{r}_k\left[X\right]=\sum _{\overrightarrow{x}\in {𝒮}_{n,k}}\langle \overrightarrow{x}❘X❘\overrightarrow{x}\rangle =\mathrm{Tr}\left[{\prod }_{k}X\right]$

where Π_k projects onto the k-particle subspace. β_j in the equation for Sim_k,s(Π_[η]) is a Majorana bivector, as discussed below. In the above equations, e_j are elementary symmetric polynomials. The elementary symmetric polynomials of degree k in d commuting variables are defined as follows:

$e_k\left(X_1,\dots ,X_d\right)\doteq \sum _{\overrightarrow{z}\in {𝒮}_{d,k}}{X}_{z_1}\dots X_{z_k}$

In addition, any element of the matrix product ω_j^†v may be computed using (1) arithmetic operations without matrix multiplication, since ω(ϕ) is sparse and ω_j^† is ω(ϕ) times some permutation.

Systems and methods of computing the estimated k-RDM {circumflex over (D)}_{{right arrow over (q)}}^{{right arrow over (p)}} with higher computational efficiency are provided below. FIG. 2 schematically shows the computing system 10 when the one or more processing devices 22 are configured to compute an estimated k-RDM {circumflex over (D)}_{{right arrow over (q)}}^{{right arrow over (p)}} of a k-RDM operator {circumflex over (D)}_{{right arrow over (q)}}^{{right arrow over (p)}} according to this more efficient approach. From each measurement of the quantum state, the following systems and methods allow any k-RDM to be estimated in (k²η) time. This increase in efficiency is particularly significant in examples in which k-RDM estimates are computed a large number of times.

The one or more processing devices 22 of the classical computing device 20 are configured to receive a plurality of classical shadow measurement results |{right arrow over (z)} associated with a fermion wavefunction ρ, as in the previous approach. The fermion wavefunction ρ describes a plurality of fermions 14 in a particle-number-conserving fermion system 30 that occupy a plurality of modes 16. The number of fermions is given by η and the number of modes 16 is given by n. Using the classical shadow measurement results |{right arrow over (z)}, the one or more processing devices 22 are further configured to compute an estimated k-RDM {circumflex over (D)}_{{right arrow over (q)}}^{{right arrow over (p)}}, where k is a number of ladder operator pairs (pairs of fermion creation and annihilation operators) that are used to define the k-RDM. The estimated k-RDM {circumflex over (D)}_{{right arrow over (q)}}^{{right arrow over (p)}} may, as in the previous approach, be computed as an expectation value of a k-RDM operator {circumflex over (D)}_{{right arrow over (q)}}^{{right arrow over (p)}}.

In the following discussion of the higher-efficiency approach to computing estimated k-RDMs {circumflex over (D)}_{{right arrow over (q)}}^{{right arrow over (p)}} the Majorana fermion representation γ_p of regular fermions is used. The Majorana representation γ_p has the following properties:

$a_p=\frac{{\gamma }_{2p-1}+i{\gamma }_{2p}}{2}$ $a_p^{†}=\frac{{\gamma }_{2p-1}-i{\gamma }_{2p}}{2}$ ${\hat{n}}_j=\frac{1}{2}\left(1+i{\beta }_j\right)$ ${\beta }_j\doteq {\gamma }_{2j-1}{\gamma }_{2j}$

In the above equations, α_p is an annihilation operator, α_p^† a creation operator, {circumflex over (n)}_j is a number operator, and β_j is a Majorana bivector. Thus, the one or more processing devices 22 are configured to compute the plurality of Majorana bivectors β_j as products of pairs of Majorana representations γ_p of the fermions 14.

For each of the classical shadow measurement results |{right arrow over (z)}, the computation of the estimated k-RDM {circumflex over (D)}_{{right arrow over (q)}}^{{right arrow over (p)}} may utilize a single-particle basis rotation matrix U(ω_j^†v) that is dependent upon the classical shadow measurement result |{right arrow over (z)}, as discussed in further detail below. When computing terms of the estimated k-RDM, the one or more processing devices 22 may be configured to compute the single-particle basis rotation matrix U(ω_j^†v) implicitly, as discussed below. Accordingly, the one or more processing devices 22 may avoid having to explicitly compute the single-particle basis rotation matrix, which may have a size of

$\left(\begin{array}{c}n\\ k\end{array}\right)\times \left(\begin{array}{c}n\\ k\end{array}\right)\mathrm{or}\left(\begin{array}{c}n\\ \eta \end{array}\right)\times \left(\begin{array}{c}n\\ \eta \end{array}\right).$

For each of the plurality of classical shadow measurement results |{right arrow over (z)}, the one or more processing devices 22 may be configured to compute a permutation matrix v_{{right arrow over (z)}}∈^n×n that maps elements of the classical shadow measurement result |{right arrow over (z)} to the fermions 14 included in the particle-number-conserving fermion system 30. In addition, the one or more processing devices 22 may be further configured to sample a unitary matrix u_H from a Haar random distribution. As in the previous approach shown in FIG. 1, the Haar random distribution may be the n×n Haar random distribution on the unitary group U_n.

When generating each of the classical shadow measurement results |{right arrow over (z)}, the one or more processing devices 22 may be further configured to utilize a sampling rotation matrix U(u_H) that encodes a rotation by the unitary matrix u_H. The one or more processing devices 22 may be further configured to sample each of the classical shadow measurement results |{right arrow over (z)} based at least in part on the fermion wavefunction ρ and on the sampling rotation matrix U(u_H). Since the sampling rotation matrix U(u_H) may have a size of

$\left(\begin{array}{c}n\\ k\end{array}\right)\times \left(\begin{array}{c}n\\ k\end{array}\right)\mathrm{or}\left(\begin{array}{c}n\\ \eta \end{array}\right)\times \left(\begin{array}{c}n\\ \eta \end{array}\right),$

the one or more processing devices 22 may be configured to compute the sampling rotation matrix U(u_H) implicitly when computing the classical shadow measurement results |{right arrow over (z)} rather than explicitly computing the sampling rotation matrix U(u_H). When the classical shadow measurement results |{right arrow over (z)} are sampled, as shown in FIG. 2, the rotations encoded in the corresponding sampling rotation matrices U(u_H) may be applied to the fermion wavefunction ρ at the quantum computing device 12 to generate a rotated fermion wavefunction UρU^† from which the classical shadow measurement results |{right arrow over (z)} are sampled. Although the example of FIG. 2 samples the classical shadow measurement |{right arrow over (z)} results from a quantum computing device 12, the fermion wavefunction ρ may instead be simulated at a classical analog simulator in other examples.

The one or more processing devices 22 may be further configured to compute the classical shadow v as a product of the permutation matrix v_{{right arrow over (z)}} and the unitary matrix u_H, such that v=V_{{right arrow over (z)}}u_H. The classical shadows v may accordingly be computed as discussed above with reference to the previous approach.

The one or more processing devices 22 may be further configured to compute a rotated unitary matrix ω_j^†v as a product of a phase matrix ω_j^† and the classical shadow v. The phase matrix ω_j^† in this example is the conjugate transpose of the phase matrix ω(ϕ) discussed above. Accordingly, in this example, the one or more processing devices 22 are further configured to compute the phase matrix ω_j^† based at least in part on the number of ladder operator pairs k, the total number of modes n, and a randomly selected phase angle ϕ.

The single-particle basis rotation matrix U(ω_j^†v) encodes a rotation by the rotated unitary matrix ω_j^†V. In the following discussion, ω_j^†v is indicated as u, such that the single-particle basis rotation matrix U(ω_j^†v) is represented as U(u). In the Majorana representation, the single-particle basis rotation matrix U(u) rotates Majorana operators to other linear combination of Majorana operators:

$U\left(u\right){\gamma }_p{U}^{†}\left(u\right)=\sum _q{\tilde{u}}_{p,q}{\gamma }_q$ $\mathrm{where}$ $\tilde{u}=\mathrm{Re}\left[u\right]\otimes \left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)+\mathrm{Im}\left[u\right]\otimes \left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right)$

The matrices in the above equation for ũ may be indicated as follows:

$\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)=I_{\left[2\right]}$ $\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right)=Y$

Using the indexing γ_p,x=γ_2p+x−1, the consistency of the Majorana representation with the above definition of the single-particle basis rotation may be verified as follows:

$\begin{array}{c}U\left(u\right)a_q^{†}{U}^{†}\left(u\right)=U\left(u\right)\frac{{\gamma }_{p,0}-i{\gamma }_{p,1}}{2}{U}^{†}\left(u\right)\\ =\frac{1}{2}\sum _q{\mathrm{Re}\left[u\right]}_{p,q}{\gamma }_{q,0}+{\mathrm{Im}\left[u\right]}_{p,q}{\gamma }_{q,1}-i\left({\mathrm{Re}\left[u\right]}_{p,q}{\gamma }_{q,1}-{\mathrm{Im}\left[u\right]}_{p,q}{\gamma }_{q,0}\right)\\ =\frac{1}{2}\sum _q\left({\mathrm{Re}\left[u\right]}_{p,q}+{i\mathrm{Im}\left[u\right]}_{p,q}\right){\gamma }_{q,0}-i\left({\mathrm{Re}\left[u\right]}_{p,q}+{i\mathrm{Im}\left[u\right]}_{p,q}\right){\gamma }_{q,1}\\ =\frac{1}{2}\sum _q{u}_{p,q}\left({\gamma }_{q,0}-i{\gamma }_{q,1}\right)=\sum _q{u}_{p,q}{a}_q^{†}\end{array}$

Given a real orthogonal matrix R∈_2n, the fermion operator γ_p and the single-particle basis rotation matrix U(R) associated with the real orthogonal matrix R satisfy the following property:

$U\left(R\right){\gamma }_p{U}^{†}\left(R\right)=\sum _q{R}_{p,q}{\gamma }_q$

where R_p,q are coefficients of the terms of the rotated fermion operator Uγ_pU^†.

Computing the estimated k-RDM includes computing a plurality of traces as terms of the estimated k-RDM {circumflex over (D)}_{{right arrow over (q)}}^{{right arrow over (p)}}. The terms U_k(u){circumflex over (D)}_[k]^[k]U_k^†(u) of the estimated k-RDM {circumflex over (D)}_{{right arrow over (q)}}^{{right arrow over (p)}}, (where the notation is simplified by absorbing ω_j into u) are linear combinations of the traces of various degrees of elementary symmetric polynomials of Majorana bivectors:

$\langle U_k\left(u\right){\hat{D}}_{\left[k\right]}^{\left[k\right]}{U}_k^{†}\left(u\right)\rangle =\sum _{j=0}^k{\alpha }_{n,\eta ,k,j}{\mathrm{Tr}}_k\left[{\hat{n}}_{\left[k\right]}{U}_k\left(u\right)e_j\left({\beta }_1,\dots ,{\beta }_{\eta }\right)U_k^{†}\left(u\right)\right]$

with some coefficients α_n,η,k,j. Accordingly, the traces computed as terms of the estimated k-RDM are computed as traces of respective products of a number operator {circumflex over (n)}_[k] of the particle-number-conserving fermion system 30, the single-particle basis rotation matrix U_k(u), the plurality of elementary symmetric polynomials e_j of a plurality of Majorana bivectors β_j respectively associated with the plurality of fermions 14, and the conjugate transpose U_k^†(u) of the single-particle basis rotation matrix U_k (u). The plurality of elementary symmetric polynomials e_j of the plurality of Majorana bivectors β_j may include the elementary symmetric polynomials with integer degrees j=0, . . . , η.

As discussed above, the estimated k-RDM {circumflex over (D)}_{{right arrow over (q)}}^{{right arrow over (p)}} may be computed as a sum of four or fewer terms U_k(u){circumflex over (D)}_[k]^[k]U_k^†(u) with corresponding coefficients. The speedup in the computation of the terms U_k(u){circumflex over (D)}_[k]^[k]U_k^†(u) arises from a faster method for computing traces of the following form:

$T{r}_r\left[{\hat{n}}_{\left[k\right]}U\left(u\right)e_j\left({\beta }_1,\dots ,{\beta }_{\eta }\right)U^{†}\left(u\right)\right]$

concurrently for each of the elementary symmetric polynomials j=0, . . . , η of the Majorana bivectors β_j. Once the traces are computed, each term U_k(u){circumflex over (D)}_[k]^[k]U_k^†(u) may be obtained by computing an appropriate linear combination of the traces.

FIG. 3 schematically shows the classical computing device 20 when the one or more processing devices 22 are configured to compute a plurality of rotated diagonal estimated k-RDM terms U_k(u){circumflex over (D)}_[k]^[k]U_k^†(u). The one or more processing devices 22, as shown in the example of FIG. 3, are configured to execute a trace computation module 40 at which the traces Tr_k[{circumflex over (n)}_[k]U(u)e_j(β₁, . . . , β_η)U^†(u)] are computed. A plurality of additional traces may also be computed at the trace computation module 40 when computing the terms U_k(u){circumflex over (D)}_[k]^[k]U_k^†(u), as discussed in further detail below.

The inner product Tr_k[{circumflex over (n)}_[k]U(u)e_j(β₁, . . . , β_η)U^†(u)] may be inferred from (n) calculations of fermionic linear optics. Fermionic linear optics is related to the strong simulability of matchgate circuits and non-interacting fermion distributions and is known to be possible in polynomial time (poly(n)). However, the approach discussed below is significantly faster.

Given any real anti-symmetric matrix M, which has block-diagonal form

$M=R\left[{\oplus }_{j\in \left[n\right]}{\lambda }_{j}Y\right]R^T,R\in {\mathrm{𝒮𝒪}}_{2n}$

the fermionic gaussian may be defined as follows:

$\rho \left(M\right)=\exp \left[\frac{i}{2}\sum _{j,k\in \left[2n\right]}{\theta }_j{M}_{j,k}{\theta }_k\right]=U\left(R\right)\prod _{j\in \left[n\right]}\left(1+i{\lambda }_j{\beta }_j\right)U^{†}\left(R\right)$

The trace of two products of fermionic gaussians evaluates to a Pfaffian as follows:

$Tr\left[\mathrm{XY}\right]={\left(-2\right)}^n\int D\theta D\mu e^{{\theta }^{T}u}\omega \left(X,\theta \right)\omega \left(Y,\mu \right)$ $\begin{array}{c}\mathrm{Tr}\left[\rho \left(M\right)\rho \left(K\right)\right]={\left(-2\right)}^n\int D\theta D\mu e^{{\theta }^{T}u}\omega \left(M,\theta \right)\omega \left(K,\mu \right)\\ ={\left(-2\right)}^n\int D\theta D\mu e^{{\theta }^{T}u}{e}^{\frac{i}{2}{\theta }^{T}M\theta }{e}^{\frac{i}{2}{\mu }^{T}K\mu }\\ ={\left(-2\right)}^n\int D\theta D\mu e^{{\theta }^{T}u+\frac{i}{2}{\mu }^{T}K\mu }{e}^{\frac{i}{2}{\theta }^{T}M\theta }\\ ={\left(-2\right)}^n{i}^n\mathrm{Pf}\left(K\right)\int D\theta e^{\frac{i}{2}{\theta }^T\left(-K^{-1}+M\right)\theta }=2^n\mathrm{Pf}\left(K\right)\mathrm{Pf}\left(M-K^{-1}\right)\\ =2^n\mathrm{Pf}\left(M\right)2^n\mathrm{Pf}\left(K-M^{-1}\right)\end{array}$

In a first step performed when computing the trace Tr_k[{circumflex over (n)}_[k]U(u)e_j(β₁, . . . , β_η)U^†(u)], the inner product {circumflex over (n)}_[k]U(u)e_j(β₁, . . . , β_η)U^†(u) in the above expression may be converted into an evaluation of at most 72 fermionic linear optics expressions Tr[ρ(M)ρ(K)] for some M and K. In a second step, the η fermionic linear optics expressions may be evaluated in a combined manner to obtain the value of the trace.

In the first step, any operator X, the trace Tr_k [{circumflex over (n)}_[k]X]=Tr_k [ρ_kX] over k particles is equal to the trace over all particles times a fermionic gaussian ρ_k. An invertible matrix ∧ may be defined as follows:

$\Lambda =I_{\left[k\right]}-I_{\left[n\right]\setminus \left[k\right]}=I_{\left[n\right]}-2I_{\left[n\right]\setminus \left[k\right]}=2I_{\left[k\right]}-I_{\left[n\right]}$

where I_[k] is an identity matrix with k rows and columns. From the definition of the fermionic gaussian,

${\rho }_k\doteq \rho \left(\Lambda \otimes y\right)=2^n\prod _{j=1}^k{\hat{n}}_j\prod _{j=k+1}^n\left(1-{\hat{n}}_j\right)$

which implies

$T{r}_k\left[{\hat{n}}_{\left[k\right]}X\right]=Tr\left[{\hat{n}}_{\left[k\right]}\prod _{j\in \left[n\right]\setminus \left[k\right]}\left(1-{\hat{n}}_j\right)X\right]=Tr\left[{\rho }_{k}X\right]$

The use of a fermionic gaussian in the above identity is a significant simplification over the previous approach, which instead inserted a generic projector Π_k onto the k-particle subspace as Tr_k[{circumflex over (n)}_[k]X]=Tr[Π_k{circumflex over (n)}_[k]X]. In contrast, in the current approach, Tr_k[{circumflex over (n)}_[k]X]=Tr[Π_kX].

The fermionic gaussian itself is a generating function for the elementary symmetric polynomials:

$\rho \left(M\right)=U\left(R\right)\prod _{j=1}^n\left(1+i\kappa {\lambda }_j{\beta }_j\right)U^{†}\left(R\right)=U\left(R\right)\sum _{j=0}^n{\left(i\kappa \right)}^j{e}_j\left({\lambda }_1{\beta }_1,\dots ,{\lambda }_n{\beta }_n\right)U^{†}\left(R\right)$

Thus, the following fermonic gaussian is a linear combination of the elementary symmetric polynomials e_j (β₁, . . . , β_η):

${\rho }_{\eta }\left(\kappa \right)\doteq \rho \left(\kappa I_{\left[n\right]}\otimes y\right)=\sum _{j=0}^{\eta }{\left(i\kappa \right)}^j{e}_j\left({\beta }_1,\dots {\beta }_{\eta }\right)$

A single evaluation of

$\mathrm{Tr}\left[{\rho }_k{U}_k\left(u\right){\rho }_{\eta }\left(\kappa \right)U_k^{†}\left(u\right)\right]=\sum _{j=0}^{\eta }{\left(i\kappa \right)}^j{\mathrm{Tr}}_k\left[U\left(u\right){\hat{n}}_{\left[k\right]}{U}^{†}\left(u\right)e_j\left({\beta }_1,{\mathrm{\dots \beta }}_{\eta }\right)\right]$

is therefore equal to a linear combination of the traces Tr_k[{circumflex over (n)}_[k]U(u)e_j(β₁, . . . , β_η)U^†(u)]. Evaluating Tr[ρ_kU_k(u)ρ_η(κ)U_k^†(u)] on (n) different values of κ then provides enough information to compute all the traces Tr_k[Π_[k]U(v)e_j({circumflex over (n)}₁, . . . , {circumflex over (n)}_η)U^†(v)], e.g., by polynomial interpolation.

In the second step, rather than computing the traces Tr_k [{circumflex over (n)}_[k]U(u)e_j (β₁, . . . , β_η)U^†(u)] via polynomial interpolation, a faster approach may instead be used. The traces may be expressed as Pfaffians of appropriately defined matrices:

$\begin{array}{c}\mathrm{Tr}\left[{\rho }_k{U}_k\left(u\right){\rho }_{\eta }\left(\kappa \right)U_k^{†}\left(u\right)\right]=\mathrm{Tr}\left[\rho \left(\Lambda \otimes Y\right)\rho \left({\kappa \mathrm{ũI}}_{\left[\eta \right]}\otimes {Yũ}^T\right)\right]\\ =\mathrm{Tr}\left[\rho \left({ũ}^T\Lambda \otimes Yũ\right)\rho \left(\kappa I_{\left[\eta \right]}\otimes Y\right)\right]\\ =2^n\mathrm{Pf}\left[{ũ}^T\Lambda \otimes Yũ\right]\mathrm{Pf}\left[\kappa I_{\left[\eta \right]}\otimes Y-{ũ}^T\Lambda \otimes Yũ\right]\\ =2^n{\left(-1\right)}^{n-k}\mathrm{Pf}\left[A\left(\kappa \right)\right]\end{array}$ $\mathrm{where}$ $A\left(\kappa \right)\stackrel{.}{=}\kappa I_{\left[\eta \right]}\otimes Y-{ũ}^T\Lambda \otimes Yũ$

Thus, the Pfaffian used to compute the traces is a Pfaffian of an intermediate matrix A(κ), which the one or more processing devices 22 are configured to compute based at least in part on the rotated unitary matrix ω_j^†v, the total number η of the plurality of fermions 14, the total number n of the plurality of modes 16, and the number of ladder operator pairs k.

Computing each of the traces includes iteratively computing a plurality of derivatives of the Pfaffian Pf[A]. By taking high-order derivatives with respect to K, the traces corresponding to specific elementary symmetric polynomials may be isolated. For example, the following derivative may be computed:

${\partial }_k^x\mathrm{Tr}\left[{\rho }_k{U}_k\left(u\right){\rho }_{\eta }\left(\kappa \right)U_k^{†}\left(u\right)\right]{|}_{\kappa =0}=x!i^x{\mathrm{Tr}}_k\left[{\hat{n}}_{\left[k\right]}U\left(u\right)e_x\left({\beta }_1,{\mathrm{\cdots \beta }}_{\eta }\right)U^{†}\left(u\right)\right]$

Accordingly, the trace may be computed for each of the elementary symmetric polynomials e_j for j=0, . . . , η. The one or more processing devices 22 may be configured to compute the derivatives of the Pfaffian Pf[A] with integer orders x=1, . . . , η.

An efficient method of computing the derivatives the Pfaffian Pf[A(K)] is provided below. The derivatives of a Pfaffian in general are given by:

$\partial \mathrm{Pf}\left(A\right)=\frac{1}{2}\mathrm{Pf}\left(A\right)\mathrm{Tr}\left[A^{-1}\partial A\right]$

For the intermediate matrix A discussed above, the derivatives are:

${\partial }_{\kappa }\mathrm{Pf}\left(A\right)=\frac{1}{2}\mathrm{Pf}\left(A\right)\mathrm{Tr}\left[A^{-1}\left(I_{\left[n\right]}\otimes Y\right)\right]$ ${\partial }_{\kappa }^2\mathrm{Pf}\left(A\right)=\frac{1}{2}{\partial }_{\kappa }\mathrm{Pf}\left(A\right)\mathrm{Tr}\left[A^{-1}\left(I_{\left[\eta \right]}\otimes Y\right)\right]-\frac{1}{2}\mathrm{Pf}\left(A\right)\mathrm{Tr}\left[A^{-1}\left(I_{\left[\eta \right]}\otimes Y\right)A^{-1}\left(I_{\left[\eta \right]}\otimes Y\right)\right]$ $\dots$ ${\partial }_{\kappa }^x\mathrm{Pf}\left(A\right)=\frac{1}{2}\sum _{j=0}^{x-1}\left(\begin{array}{c}x-1\\ j\end{array}\right){\partial }_{\kappa }^{x-j-1}{\left(-1\right)}^j\mathrm{Pf}\left(A\right)\mathrm{Tr}\left[{\left(A^{-1}\left(I_{\left[\eta \right]}\otimes Y\right)\right)}^{j+1}\right]$

This recursive computation allows the one or more processing devices 22 to compute higher-order derivatives from lower-order derivatives. Assuming that all the traces Tr[(A⁻¹(I,_[η]⊗Y))^j]|_κ=0 have been computed, the recursive computation of each of the derivatives ∀x∈[η], ∂_κ^xPf(A)|_κ=0 may be performed in (η²) time.

The one or more processing devices 22 are further configured to compute the traces included in the above equation for the derivatives ∂_κ^xPf(A) of the Pfaffian. As a base case for the intermediate matrix,

$A\left(0\right)=A^{-1}\left(0\right)=-{ũ}^T\Lambda \otimes Yũ$

Hence, the trace is given by:

$\begin{array}{c}\mathrm{Tr}\left[{\left(A^{-1}\left(I_{\left[\eta \right]}\otimes Y\right)\right)}^j\right]{❘}_{\kappa =0}=\mathrm{Tr}\left[{\left(\left(-{ũ}^T\Lambda \otimes Yũ\right)\left(I_{\left[\eta \right]}\otimes Y\right)\right)}^j\right]\\ ={\left(-1\right)}^j\mathrm{Tr}\left[{\left(\left(\Lambda \otimes Y\right)\left(ũI_{\left[\eta \right]}\otimes Y{ũ}^T\right)\right)}^j\right]\end{array}$

The trace Tr[(A⁻¹(I_[η]⊗Y))^j]|_κ=0 may be computed at least in part by computing each of the eigenvalues of (Λ⊗Y)(ũI_[n]⊗Yũ^T).

Computing eigenvalues of a 2n×2n matrix has a time complexity of (n^ω) and enables straightforward computation of the trace. Such an approach may be used in the non-particle-conserving case. However, when particle number is conserved, the above computation of the trace may be reduced to finding the eigenvalues of an even smaller 2k×2k matrix. Using the identity Λ=2I[_k]−I_[n], the trace may be computed as follows:

$\begin{array}{c}\mathrm{Tr}\left[{\left(A^{-1}\left(I_{\left[n\right]}\otimes Y\right)\right)}^j\right]{❘}_{\kappa =0}={\left(-1\right)}^j\left[{\left(\left({ũ}^T\left(2I_{\left[k\right]}-I_{\left[n\right]}\right)\otimes Yũ\right)\left(I_{\left[\eta \right]}\otimes Y\right)\right)}^j\right]\\ ={\left(-1\right)}^j[(2\left({ũ}^T{I}_{\left[k\right]}\otimes Yũ\right)\left(I_{\left[\eta \right]}\otimes Y\right)-\\ {\left({ũ}^T{I}_{\left[n\right]}\otimes Yũ\right)\left(I_{\left[\eta \right]}\otimes Y\right))}^j]\end{array}$

The following identities may be used to simplify the computation of the above trace:

$\begin{array}{c}\left({ũ}^T{I}_{\left[n\right]}\otimes Yũ\right)\left(I_{\left[\eta \right]}\otimes Y\right)=\left({\mathrm{Re}\left[u\right]}^T\otimes I_{\left[2\right]}-{\mathrm{Im}\left[u\right]}^T\otimes Y\right)\left(I_{\left[n\right]}\otimes Y\right)\\ \left(\mathrm{Re}\left[u\right]\otimes I_{\left[2\right]}+\mathrm{Im}\left[u\right]\otimes Y\right)\left(I_{\left[\eta \right]}\otimes Y\right)\\ =\left({\mathrm{Re}\left[u\right]}^T\otimes Y+{\mathrm{Im}\left[u\right]}^T\otimes I_{\left[2\right]}\right)\\ \left(\mathrm{Re}\left[u\right]\otimes I_{\left[2\right]}+\mathrm{Im}\left[u\right]\otimes Y\right)\left(I_{\left[\eta \right]}\otimes Y\right)\\ =({\mathrm{Re}\left[u\right]}^T\mathrm{Re}\left[u\right]\otimes Y-{\mathrm{Re}\left[u\right]}^T\mathrm{Im}\left[u\right]\otimes \\ I_{\left[2\right]}+{\mathrm{Im}\left[u\right]}^T\mathrm{Re}\left[u\right]\otimes I_{\left[2\right]}+{\mathrm{Im}\left[u\right]}^T\mathrm{Im}\left[u\right]\otimes Y)\\ \left(I_{\left[\eta \right]}\otimes Y\right)=\left(I_{\left[n\right]}\otimes Y\right)\left(I_{\left[\eta \right]}\otimes Y\right)\\ =-I_{\left[\eta \right]}\otimes I_{\left[2\right]}\end{array}$

In addition,

$\begin{array}{c}\left({ũ}^T{I}_{\left[k\right]}\otimes Yũ\right)\left(I_{\left[\eta \right]}\otimes Y\right)=\left({\mathrm{Re}\left[u\right]}^T\otimes I_{\left[2\right]}-{\mathrm{Im}\left[u\right]}^T\otimes Y\right)\left(I_{\left[k\right]}\otimes Y\right)\\ \left(\mathrm{Re}\left[u\right]\otimes I_{\left[2\right]}+\mathrm{Im}\left[u\right]\otimes Y\right)\left(I_{\left[\eta \right]}\otimes Y\right)\\ =\text{}\left({\mathrm{Re}\left[u\right]}^T\otimes Y+{\mathrm{Im}\left[u\right]}^T\otimes I_{\left[2\right]}\right)\\ \left(I_{\left[k\right]}\otimes I_{\left[2\right]}\right)\left(\mathrm{Re}\left[u\right]\otimes Y-\mathrm{Im}\left[u\right]\otimes I_{\left[2\right]}\right)\\ \left(I_{\left[\eta \right]}\otimes I_{\left[2\right]}\right)=-\left({\widetilde{\iota u}}^T{I}_{\left[k\right]}\otimes I_{\left[2\right]}\widetilde{\iota u}\right)\left(I_{\left[\eta \right]}\otimes I_{\left[2\right]}\right)\end{array}$ $\mathrm{where}$ $\widetilde{\iota u}=-\mathrm{Im}\left[u\right]\otimes I_{\left[2\right]}+\mathrm{Re}\left[u\right]\otimes Y$

Using the above identities, the binomial expansion, and the cyclic property of the trace, Tr[(A⁻¹(I_[η]⊗Y))^j]|_κ=0 may be computed as follows:

$\begin{array}{c}\mathrm{Tr}\left[{\left(A^{-1}\left(I_{\left[\eta \right]}\otimes Y\right)\right)}^j\right]{❘}_{\kappa =0}={\left(-1\right)}^j\mathrm{Tr}[(-2\left({\left(\widetilde{\iota u}\right)}^T{I}_{\left[k\right]}\otimes I_{\left[2\right]}\widetilde{\iota u}\right)\\ {\left(I_{\left[\eta \right]}\otimes I_{\left[2\right]}\right)+I_{\left[\eta \right]}\otimes I_{\left[2\right]})}^j]\\ ={\left(-1\right)}^j\sum _{y=0}^j{\left(-2\right)}^y\left(\begin{array}{c}j\\ y\end{array}\right)\mathrm{Tr}[I_{\left[\eta \right]}\otimes \\ {I_{\left[2\right]}[{\left(\widetilde{\iota u}\right)}^T{I}_{\left[k\right]}\otimes I_{\left[2\right]}\widetilde{\iota u})\left(I_{\left[\eta \right]}\otimes I_{\left[2\right]}\right)]}^y]\\ ={\left(-1\right)}^j[\mathrm{Tr}\left[I_{\left[\eta \right]}\otimes I_{\left[2\right]}\right]+\sum _{y=1}^j{\left(-2\right)}^y\left(\begin{array}{c}j\\ y\end{array}\right)\mathrm{Tr}\\ {[[{\left(\widetilde{\iota u}\right)}^T{I}_{\left[k\right]}\otimes I_{\left[2\right]}\widetilde{\iota u})\left(I_{\left[\eta \right]}\otimes I_{\left[2\right]}\right)]}^y]\\ ={\left(-1\right)}^j[2\eta +\sum _{y=1}^j{\left(-2\right)}^y\left(\begin{array}{c}j\\ y\end{array}\right)\\ \mathrm{Tr}\left[{\left[\left(I_{\left[k\right]}\otimes I_{\left[2\right]}\right)\left(\widetilde{\iota u}{I}_{\left[\eta \right]}\otimes I_{\left[2\right]}{\widetilde{\iota u}}^T\right)\right]}^y\right]]\end{array}$

A matrix M within the trace in the above equation may be defined as:

$M=\left(I_{\left[k\right]}\otimes I_{\left[2\right]}\right)\left(\widetilde{\iota u}{I}_{\left[\eta \right]}\otimes I_{\left[2\right]}{\widetilde{\iota u}}^T\right)$

The matrix M is non-zero on only a 2k×2k block. M may be rewritten as M=m·m^T, where

$m=\left(I_{\left[k\right]}\otimes I_{\left[2\right]}\right)\left(\widetilde{\iota u}{I}_{\left[\eta \right]}\otimes I_{\left[2\right]}\right)$

is a 2k×2η submatrix of . Thus, M may be obtained by matrix multiplication of a 2k×2η matrix with a 2η×2k matrix, which has a time complexity of (k²η^α) for some matrix multiplication exponent α≤1.

Subsequently to computing the matrix M, the one or more processing devices 22 may be further configured to compute each of the nonzero eigenvalues of M with a time complexity of (k^ω). The one or more processing devices 22 may be further configured to compute each of the powers of the eigenvalues of M (such as by repeated squaring) in (ηk) multiplications, thereby allowing each of the traces Tr[M^y] for y∈[η] to be computed with an overall time complexity of (k²η^α+k^ω+ηk). Using a cubic-time algorithm for matrix multiplication, the traces Tr[M^y] may be computed with a time complexity of (k²η), which is dominated by the time spent forming M.

Once the traces Tr[M^y] have been computed, the one or more processing devices 22 may be further configured to recursively compute the traces Tr [(A⁻¹(I_[η]⊗Y))^j]|_k=0 for j∈[η] with a time complexity of (η²). Accordingly, computing the derivatives ∂_k^xPf(A) of the Pfaffian may have an overall time complexity of (η²+k²η). The computation of the derivatives of the Pfaffian may be the highest-time-complexity step of the k-RDM estimation. Accordingly, the one or more processing devices 22 may be configured to compute the estimated k-RDM with a time complexity of (η²+k²η).

Subsequently to computing the estimated k-RDM {circumflex over (D)}_{{right arrow over (q)}}^{{right arrow over (p)}}, the one or more processing devices 22 are further configured to output the estimated k-RDM {circumflex over (D)}_{{right arrow over (q)}}^{{right arrow over (p)}} to an additional computing process. At the additional computing process, the one or more processing devices 22 may be further configured to estimate an observable O of the fermion wavefunction based at least in part on the estimated k-RDM {circumflex over (D)}_{{right arrow over (q)}}^{{right arrow over (p)}}. The one or more processing devices 22 may, for example, be configured to compute the estimated observable Ô as a linear combination of a plurality of estimated k-RDMs {circumflex over (D)}_{{right arrow over (q)}}^{{right arrow over (p)}}.

FIG. 4A shows a flowchart of a method 100 for use with a computing device to compute an estimated k-RDM. The computing device may include a classical computing device that is communicatively coupled to a quantum computing device.

The steps of the method 100 shown in FIG. 4A may be performed at the classical computing device.

At step 102, the method 100 includes receiving a plurality of classical shadow measurement results associated with a fermion wavefunction. The fermion wavefunction describes a plurality of fermions in a particle-number-conserving fermion system that occupy a plurality of modes.

At step 104, the method 100 further includes computing an estimated k-reduced density matrix (k-RDM), where k is a number of ladder operator pairs. Each of the ladder operator pairs includes a fermion creation operator and a fermion annihilation operator.

Step 104 includes step 106, which is performed for each of the classical shadow measurement results. At step 106, computing the estimated k-RDM further includes computing, as terms of the estimated k-RDM, a plurality of traces of respective products of: a number operator of the particle-number-conserving fermion system, a single-particle basis rotation matrix, a plurality of elementary symmetric polynomials of a plurality of Majorana bivectors respectively associated with the plurality of fermions, and a conjugate transpose of the single-particle basis rotation matrix. The single-particle basis rotation matrix is dependent upon the classical shadow measurement result, although the single-particle basis rotation matrix is not explicitly computed in the example of FIG. 4A. The estimated k-RDM may be computed as a sum of four or fewer of the terms computed in step 106. Computing each of the traces at step 106 includes, at step 108, iteratively computing a plurality of derivatives of a Pfaffian.

At step 110, the method 100 further includes outputting the estimated k-RDM. The estimated k-RDM may be output to an additional computing process at which one or more estimated values of an observable are computed using the estimated k-RDM. For example, the additional computing process may be a quantum chemistry simulation program. The observable may, for example, be an energy, a density, a dipole moment, or some other observable of the fermion system. Observables of the particle-number-conserving fermion system may be expressed as linear combinations of k-RDMs.

FIG. 4B shows additional steps of the method 100 that may be performed for each of the classical shadow measurement results. At step 112, the method 100 may further include generating the classical shadow measurement results. At step 114, step 112 may include sampling a unitary matrix from a Haar random distribution. The unitary matrix may have a dimension equal to the number of modes.

At step 116, step 112 may further include sampling the classical shadow measurement result from a rotated fermion wavefunction. The rotated fermion wavefunction may be computed based at least in part on the fermion wavefunction and on a sampling rotation matrix that encodes a rotation by the unitary matrix. The plurality of classical shadow measurement results may be received from the quantum computing device, at which the rotated fermion wavefunction may be prepared. The rotated fermion wavefunction may be a product of the sampling rotation matrix, the fermion wavefunction, and a conjugate transpose of the sampling rotation matrix. The classical shadow measurement results may be sampled without explicitly computing the sampling rotation matrix, which may be inefficient to compute.

At step 118, for each of the classical shadow measurement results, the method 100 may further include computing a permutation matrix that maps elements of the classical shadow measurement result to the fermions included in the particle-number-conserving fermion system. The method 100 may further include, at step 120, computing a classical shadow as a product of the permutation matrix and the unitary matrix.

At step 122, the method 100 may further include computing a rotated unitary matrix as a product of a phase matrix and the classical shadow. The single-particle basis rotation matrix may encode a rotation by the rotated unitary matrix. The phase matrix may be computed based at least in part on the number of ladder operator pairs, the total number of modes, and a randomly selected phase angle. Accordingly, the single-particle basis rotation matrix used when computing the k-RDM depends upon the classical shadow.

FIG. 4C shows additional steps of the method 100 that may be performed when performing steps 106 and 108. At step 124, the method 100 may further include computing an intermediate matrix based at least in part on the rotated unitary matrix, a total number of the plurality of fermions, a total number of the plurality of modes, and the number of ladder operator pairs. At step 126, the method 100 may further include computing the Pfaffian as a Pfaffian of the intermediate matrix. In addition, at step 128, the method 100 may further include computing the derivatives of the Pfaffian with integer orders from one to the number of fermions. For example, computing the derivatives of the Pfaffian may include iteratively computing traces of powers of an additional matrix generated from the rotated unitary matrix.

In examples in which the steps of FIG. 4C are performed, the estimated k-RDM may be computed with a time complexity of (η²+k²η), where η is a total number of the plurality of fermions. This time complexity represents a significant speedup over the (n⁵) time complexity of prior approaches to k-RDM estimation.

In some embodiments, the methods and processes described herein may be tied to a computing system of one or more computing devices. In particular, such methods and processes may be implemented as a computer-application program or service, an application-programming interface (API), a library, and/or other computer-program product.

FIG. 5 schematically shows a non-limiting embodiment of a computing system 200 that can enact one or more of the methods and processes described above. Computing system 200 is shown in simplified form. Computing system 200 may embody the computing system 10 described above and illustrated in FIG. 2. Computing system 200 may take the form of one or more personal computers, server computers, tablet computers, home-entertainment computers, network computing devices, video game devices, mobile computing devices, mobile communication devices (e.g., smart phone), and/or other computing devices, and wearable computing devices such as smart wristwatches and head mounted augmented reality devices.

Computing system 200 includes a logic processor 202 volatile memory 204, and a non-volatile storage device 206. Computing system 200 may optionally include a display subsystem 208, input subsystem 210, communication subsystem 212, and/or other components not shown in FIG. 5.

Logic processor 202 includes one or more physical devices configured to execute instructions. For example, the logic processor may be configured to execute instructions that are part of one or more applications, programs, routines, libraries, objects, components, data structures, or other logical constructs. Such instructions may be implemented to perform a task, implement a data type, transform the state of one or more components, achieve a technical effect, or otherwise arrive at a desired result.

The logic processor may include one or more physical processors configured to execute software instructions. Additionally or alternatively, the logic processor may include one or more hardware logic circuits or firmware devices configured to execute hardware-implemented logic or firmware instructions. Processors of the logic processor 202 may be single-core or multi-core, and the instructions executed thereon may be configured for sequential, parallel, and/or distributed processing. Individual components of the logic processor optionally may be distributed among two or more separate devices, which may be remotely located and/or configured for coordinated processing. Aspects of the logic processor may be virtualized and executed by remotely accessible, networked computing devices configured in a cloud-computing configuration. In such a case, these virtualized aspects are run on different physical logic processors of various different machines, it will be understood.

Non-volatile storage device 206 includes one or more physical devices configured to hold instructions executable by the logic processors to implement the methods and processes described herein. When such methods and processes are implemented, the state of non-volatile storage device 206 may be transformed—e.g., to hold different data.

Non-volatile storage device 206 may include physical devices that are removable and/or built-in. Non-volatile storage device 206 may include optical memory, semiconductor memory, and/or magnetic memory, or other mass storage device technology. Non-volatile storage device 206 may include nonvolatile, dynamic, static, read/write, read-only, sequential-access, location-addressable, file-addressable, and/or content-addressable devices. It will be appreciated that non-volatile storage device 206 is configured to hold instructions even when power is cut to the non-volatile storage device 206.

Volatile memory 204 may include physical devices that include random access memory. Volatile memory 204 is typically utilized by logic processor 202 to temporarily store information during processing of software instructions. It will be appreciated that volatile memory 204 typically does not continue to store instructions when power is cut to the volatile memory 204.

Aspects of logic processor 202, volatile memory 204, and non-volatile storage device 206 may be integrated together into one or more hardware-logic components. Such hardware-logic components may include field-programmable gate arrays (FPGAs), program- and application-specific integrated circuits (PASIC/ASICs), program- and application-specific standard products (PSSP/ASSPs), system-on-a-chip (SOC), and complex programmable logic devices (CPLDs), for example.

The terms “module,” “program,” and “engine” may be used to describe an aspect of computing system 200 typically implemented in software by a processor to perform a particular function using portions of volatile memory, which function involves transformative processing that specially configures the processor to perform the function. Thus, a module, program, or engine may be instantiated via logic processor 202 executing instructions held by non-volatile storage device 206, using portions of volatile memory 204. It will be understood that different modules, programs, and/or engines may be instantiated from the same application, service, code block, object, library, routine, API, function, etc. Likewise, the same module, program, and/or engine may be instantiated by different applications, services, code blocks, objects, routines, APIs, functions, etc. The terms “module,” “program,” and “engine” may encompass individual or groups of executable files, data files, libraries, drivers, scripts, database records, etc.

When included, display subsystem 208 may be used to present a visual representation of data held by non-volatile storage device 206. The visual representation may take the form of a graphical user interface (GUI). As the herein described methods and processes change the data held by the non-volatile storage device, and thus transform the state of the non-volatile storage device, the state of display subsystem 208 may likewise be transformed to visually represent changes in the underlying data. Display subsystem 208 may include one or more display devices utilizing virtually any type of technology. Such display devices may be combined with logic processor 202, volatile memory 204, and/or non-volatile storage device 206 in a shared enclosure, or such display devices may be peripheral display devices.

When included, input subsystem 210 may comprise or interface with one or more user-input devices such as a keyboard, mouse, touch screen, or game controller. In some embodiments, the input subsystem may comprise or interface with selected natural user input (NUI) componentry. Such componentry may be integrated or peripheral, and the transduction and/or processing of input actions may be handled on- or off-board. Example NUI componentry may include a microphone for speech and/or voice recognition; an infrared, color, stereoscopic, and/or depth camera for machine vision and/or gesture recognition; a head tracker, eye tracker, accelerometer, and/or gyroscope for motion detection and/or intent recognition; as well as electric-field sensing componentry for assessing brain activity; and/or any other suitable sensor.

When included, communication subsystem 212 may be configured to communicatively couple various computing devices described herein with each other, and with other devices. Communication subsystem 212 may include wired and/or wireless communication devices compatible with one or more different communication protocols. As non-limiting examples, the communication subsystem may be configured for communication via a wireless telephone network, or a wired or wireless local- or wide-area network. In some embodiments, the communication subsystem may allow computing system 200 to send and/or receive messages to and/or from other devices via a network such as the Internet.

The following paragraphs provide additional support for the claims of the subject application. According to one aspect, a computing device is provided, comprising one or more processing devices configured to receive a plurality of classical shadow measurement results associated with a fermion wavefunction, wherein the fermion wavefunction describes a plurality of fermions in a particle-number-conserving fermion system that occupy a plurality of modes. The one or more processing devices are further configured to compute an estimated k-reduced density matrix (k-RDM), where k is a number of ladder operator pairs, at least in part by, for each of the classical shadow measurement results, computing, as terms of the estimated k-RDM, a plurality of traces of respective products of: (1) a number operator of the particle-number-conserving fermion system; a single-particle basis rotation matrix dependent upon the classical shadow measurement result; (2) a plurality of elementary symmetric polynomials of a plurality of Majorana bivectors respectively associated with the plurality of fermions; and (3) a conjugate transpose of the single-particle basis rotation matrix, wherein computing each of the traces includes iteratively computing a plurality of derivatives of a Pfaffian. The one or more processors are further configured to output the estimated k-RDM.

In this aspect, the one or more processing devices can be configured to receive the plurality of classical shadow measurement results from a quantum computing device.

In this aspect, the one or more processing devices can be further configured to generate the plurality of classical shadow measurement results at least in part by, for each of the classical shadow measurement results: sampling a unitary matrix from a Haar random distribution; and sampling the classical shadow measurement result from a rotated fermion wavefunction, wherein the rotated fermion wavefunction is computed based at least in part on the fermion wavefunction and on a sampling rotation matrix that encodes a rotation by the unitary matrix.

In this aspect, for each of the plurality of classical shadow measurement results, the one or more processing devices can be further configured to: compute a permutation matrix that maps elements of the classical shadow measurement result to the fermions included in the particle-number-conserving fermion system; and compute a classical shadow as a product of the permutation matrix and the unitary matrix. Further in this aspect, the one or more processing devices can be further configured to compute a rotated unitary matrix as a product of a phase matrix and the classical shadow, wherein the single-particle basis rotation matrix encodes a rotation by the rotated unitary matrix. Further in this aspect, the Pfaffian can be a Pfaffian of an intermediate matrix; and the one or more processing devices can be configured to compute the intermediate matrix based at least in part on the rotated unitary matrix, a total number of the plurality of fermions, a total number of the plurality of modes, and the number of ladder operator pairs. Still further in this aspect, the one or more processing devices can be further configured to compute the phase matrix based at least in part on the number of ladder operator pairs, the total number of modes, and a randomly selected phase angle.

In this aspect, the one or more processing devices can be further configured to compute the plurality of Majorana bivectors as products of pairs of Majorana representations of the fermions.

In this aspect, the plurality of elementary symmetric polynomials of the plurality of Majorana bivectors can include the elementary symmetric polynomials with integer degrees from zero to the number of fermions.

In this aspect, the one or more processing devices can be configured to compute the derivatives of the Pfaffian with integer orders from one to the number of fermions.

In this aspect, the one or more processing devices can be configured to compute the estimated k-RDM with a time complexity of (η²+k²η), where η is a total number of the plurality of fermions.

In this aspect, the one or more processing devices can be further configured to estimate an observable of the fermion wavefunction based at least in part on the estimated k-RDM.

According to another aspect, a method for use with a computing device is provided, the method comprising: receiving a plurality of classical shadow measurement results associated with a fermion wavefunction, wherein the fermion wavefunction describes a plurality of fermions in a particle-number-conserving fermion system that occupy a plurality of modes. The method further comprises computing an estimated k-reduced density matrix (k-RDM), where k is a number of ladder operator pairs, at least in part by, for each of the classical shadow measurement results, computing, as terms of the estimated k-RDM, a plurality of traces of respective products of: (1) a number operator of the particle-number-conserving fermion system; (2) a single-particle basis rotation matrix dependent upon the classical shadow measurement result; (3) a plurality of elementary symmetric polynomials of a plurality of Majorana bivectors respectively associated with the plurality of fermions; and (4) a conjugate transpose of the single-particle basis rotation matrix. Each of the traces includes iteratively computing a plurality of derivatives of a Pfaffian. The method further comprises outputting the estimated k-RDM.

In this aspect, the plurality of classical shadow measurement results can be received from a quantum computing device.

In this aspect, the method can further comprise generating the plurality of classical shadow measurement results at least in part by, for each of the classical shadow measurement results: sampling a unitary matrix from a Haar random distribution; and sampling the classical shadow measurement result from a rotated fermion wavefunction, wherein the rotated fermion wavefunction is computed based at least in part on the fermion wavefunction and on a sampling rotation matrix that encodes a rotation by the unitary matrix.

In this aspect, the method can further comprise, for each of the classical shadow measurement results: computing a permutation matrix that maps elements of the classical shadow measurement result to the fermions included in the particle-number-conserving fermion system; and computing a classical shadow as a product of the permutation matrix and the unitary matrix.

In this aspect, the method can further comprise computing a rotated unitary matrix as a product of a phase matrix and the classical shadow, wherein the single-particle basis rotation matrix encodes a rotation by the rotated unitary matrix. Further in this aspect, the Pfaffian can be a Pfaffian of an intermediate matrix; and the method can further include computing the intermediate matrix based at least in part on the rotated unitary matrix, a total number of the plurality of fermions, a total number of the plurality of modes, and the number of ladder operator pairs.

In this aspect, the estimated k-RDM can be computed with a time complexity of (η²+k²η), where η is a total number of the plurality of fermions.

According to another aspect, a computing system is provided, comprising a quantum computing device and a classical computing device. The quantum computing device is configured to: prepare a fermion wavefunction, wherein the fermion wavefunction describes a plurality of fermions in a particle-number-conserving fermion system that occupy a plurality of modes; and measure a plurality of classical shadow measurement results associated with the fermion wavefunction. The classical computing device includes one or more processing devices configured to: receive the plurality of classical shadow measurement results; and compute an estimated k-reduced density matrix (k-RDM), where k is a number of ladder operator pairs. The estimated k-RDM can be computed at least in part by, for each of the classical shadow measurement results, computing, as terms of the estimated k-RDM, a plurality of traces of respective products of: (1) a number operator of the particle-number-conserving fermion system; (2) a single-particle basis rotation matrix dependent upon the classical shadow measurement result; (3) a plurality of elementary symmetric polynomials of a plurality of Majorana bivectors respectively associated with the plurality of fermions; and a conjugate transpose of the single-particle basis rotation matrix, wherein computing each of the traces includes iteratively computing a plurality of derivatives of a Pfaffian. The one or more processing devices are further configured to output the estimated k-RDM.

It will be understood that the configurations and/or approaches described herein are exemplary in nature, and that these specific embodiments or examples are not to be considered in a limiting sense, because numerous variations are possible. The specific routines or methods described herein may represent one or more of any number of processing strategies. As such, various acts illustrated and/or described may be performed in the sequence illustrated and/or described, in other sequences, in parallel, or omitted. Likewise, the order of the above-described processes may be changed.

The subject matter of the present disclosure includes all novel and non-obvious combinations and sub-combinations of the various processes, systems and configurations, and other features, functions, acts, and/or properties disclosed herein, as well as any and all equivalents thereof.

Claims

1. A computing device comprising: one or more processing devices configured to: receive a plurality of classical shadow measurement results associated with a fermion wavefunction, wherein the fermion wavefunction describes a plurality of fermions in a particle-number-conserving fermion system that occupy a plurality of modes;compute an estimated k-reduced density matrix (k-RDM), where k is a number of ladder operator pairs, at least in part by, for each of the classical shadow measurement results: computing, as terms of the estimated k-RDM, a plurality of traces of respective products of: a number operator of the particle-number-conserving fermion system;a single-particle basis rotation matrix dependent upon the classical shadow measurement result;a plurality of elementary symmetric polynomials of a plurality of Majorana bivectors respectively associated with the plurality of fermions; anda conjugate transpose of the single-particle basis rotation matrix,wherein computing each of the traces includes iteratively computing a plurality of derivatives of a Pfaffian; andoutput the estimated k-RDM.

2. The computing device of claim 1, wherein the one or more processing devices are configured to receive the plurality of classical shadow measurement results from a quantum computing device.

3. The computing device of claim 2, wherein the one or more processing devices are further configured to generate the plurality of classical shadow measurement results at least in part by, for each of the classical shadow measurement results: sampling a unitary matrix from a Haar random distribution; andsampling the classical shadow measurement result from a rotated fermion wavefunction, wherein the rotated fermion wavefunction is computed based at least in part on the fermion wavefunction and on a sampling rotation matrix that encodes a rotation by the unitary matrix.

4. The computing device of claim 3, wherein, for each of the plurality of classical shadow measurement results, the one or more processing devices are further configured to: compute a permutation matrix that maps elements of the classical shadow measurement result to the fermions included in the particle-number-conserving fermion system; andcompute a classical shadow as a product of the permutation matrix and the unitary matrix.

5. The computing device of claim 4, wherein the one or more processing devices are further configured to compute a rotated unitary matrix as a product of a phase matrix and the classical shadow, wherein the single-particle basis rotation matrix encodes a rotation by the rotated unitary matrix.

6. The computing device of claim 5, wherein: the Pfaffian is a Pfaffian of an intermediate matrix; andthe one or more processing devices are configured to compute the intermediate matrix based at least in part on the rotated unitary matrix, a total number of the plurality of fermions, a total number of the plurality of modes, and the number of ladder operator pairs.

7. The computing device of claim 6, wherein the one or more processing devices are further configured to compute the phase matrix based at least in part on the number of ladder operator pairs, the total number of modes, and a randomly selected phase angle.

8. The computing device of claim 1, wherein the one or more processing devices are further configured to compute the plurality of Majorana bivectors as products of pairs of Majorana representations of the fermions.

9. The computing device of claim 1, wherein the plurality of elementary symmetric polynomials of the plurality of Majorana bivectors includes the elementary symmetric polynomials with integer degrees from zero to the number of fermions.

10. The computing device of claim 1, wherein the one or more processing devices are configured to compute the derivatives of the Pfaffian with integer orders from one to the number of fermions.

11. The computing device of claim 1, wherein the one or more processing devices are configured to compute the estimated k-RDM with a time complexity of (η2+k2η), where η is a total number of the plurality of fermions.

12. The computing device of claim 1, wherein the one or more processing devices are further configured to estimate an observable of the fermion wavefunction based at least in part on the estimated k-RDM.

13. A method for use with a computing device, the method comprising: receiving a plurality of classical shadow measurement results associated with a fermion wavefunction, wherein the fermion wavefunction describes a plurality of fermions in a particle-number-conserving fermion system that occupy a plurality of modes;computing an estimated k-reduced density matrix (k-RDM), where k is a number of ladder operator pairs, at least in part by, for each of the classical shadow measurement results: computing, as terms of the estimated k-RDM, a plurality of traces of respective products of: a number operator of the particle-number-conserving fermion system;a single-particle basis rotation matrix dependent upon the classical shadow measurement result;a plurality of elementary symmetric polynomials of a plurality of Majorana bivectors respectively associated with the plurality of fermions; anda conjugate transpose of the single-particle basis rotation matrix,wherein computing each of the traces includes iteratively computing a plurality of derivatives of a Pfaffian; andoutputting the estimated k-RDM.

14. The method of claim 13, wherein the plurality of classical shadow measurement results are received from a quantum computing device.

15. The method of claim 13, further comprising generating the plurality of classical shadow measurement results at least in part by, for each of the classical shadow measurement results: sampling a unitary matrix from a Haar random distribution; andsampling the classical shadow measurement result from a rotated fermion wavefunction, wherein the rotated fermion wavefunction is computed based at least in part on the fermion wavefunction and on a sampling rotation matrix that encodes a rotation by the unitary matrix.

16. The method of claim 15, further comprising, for each of the classical shadow measurement results: computing a permutation matrix that maps elements of the classical shadow measurement result to the fermions included in the particle-number-conserving fermion system; andcomputing a classical shadow as a product of the permutation matrix and the unitary matrix.

17. The method of claim 16, further comprising computing a rotated unitary matrix as a product of a phase matrix and the classical shadow, wherein the single-particle basis rotation matrix encodes a rotation by the rotated unitary matrix.

18. The method of claim 17, wherein: the Pfaffian is a Pfaffian of an intermediate matrix; andthe method further includes computing the intermediate matrix based at least in part on the rotated unitary matrix, a total number of the plurality of fermions, a total number of the plurality of modes, and the number of ladder operator pairs.

19. The method of claim 13, wherein the estimated k-RDM is computed with a time complexity of (η2+k2η), where η is a total number of the plurality of fermions.

20. A computing system comprising: a quantum computing device configured to: prepare a fermion wavefunction, wherein the fermion wavefunction describes a plurality of fermions in a particle-number-conserving fermion system that occupy a plurality of modes; andmeasure a plurality of classical shadow measurement results associated with the fermion wavefunction; anda classical computing device including one or more processing devices configured to: receive the plurality of classical shadow measurement results;compute an estimated k-reduced density matrix (k-RDM), where k is a number of ladder operator pairs, at least in part by, for each of the classical shadow measurement results: computing, as terms of the estimated k-RDM, a plurality of traces of respective products of: a number operator of the particle-number-conserving fermion system;a single-particle basis rotation matrix dependent upon the classical shadow measurement result;a plurality of elementary symmetric polynomials of a plurality of Majorana bivectors respectively associated with the plurality of fermions; anda conjugate transpose of the single-particle basis rotation matrix,wherein computing each of the traces includes iteratively computing a plurality of derivatives of a Pfaffian; andoutput the estimated k-RDM.