# On sampling determinantal and Pfaffian point processes on a quantum computer

Rémi Bardenet, Michaël Fanuel, Alexandre Feller

<figure><img src="image_1.png" /></figure>

► To cite this version:

Rémi Bardenet, Michaël Fanuel, Alexandre Feller. On sampling determinantal and Pfaffian point processes on a quantum computer. Journal of Physics A: Mathematical and Theoretical, 2024, 57 (5), pp.055202. 10.1088/1751-8121/ad1b75. hal-04107251

HAL Id: hal-04107251

https://hal.science/hal-04107251

Submitted on 26 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On sampling determinantal and Pfaffian point processes on a quantum computer

Rémi Bardenet, Michaël Fanuel, and Alexandre Feller

Université de Lille, CNRS, Centrale Lille
UMR 9189 - CRISTAL, F-59000 Lille, France

## Abstract

DPPs were introduced by Macchi as a model in quantum optics the 1970s. Since then, they have been widely used as models and subsam- pling tools in statistics and computer science. Most applications require sampling from a DPP, and given their quantum origin, it is natural to wonder whether sampling a DPP on a quantum computer is easier than on a classical one. We focus here on DPPs over a finite state space, which are distributions over the subsets of {1, ..., N} parametrized by an $N \times N$ Hermitian kernel matrix. Vanilla sampling consists in two steps, of respec- tive costs $O(N^3)$ and $O(Nr^2)$ operations on a classical computer, where r is the rank of the kernel matrix. A large first part of the current pa- per consists in explaining why the state-of-the-art in quantum simulation of fermionic systems already yields quantum DPP sampling algorithms. We then modify existing quantum circuits, and discuss their insertion in a full DPP sampling pipeline that starts from practical kernel specifica- tions. The bottom line is that, with P (classical) parallel processors, we can divide the preprocessing cost by P and build a quantum circuit with $O(Nr)$ gates that sample a given DPP, with depth varying from $O(N)$ to $O(r \log N)$ depending on qubit-communication constraints on the target machine. We also connect existing work on the simulation of supercon- ductors to Pfaffian point processes, which generalize DPPs and would be a natural addition to the machine learner's toolbox. Finally, the circuits are empirically validated on a classical simulator and on 5-qubit machines.

**Keywords**— Determinantal and Pfaffian point processes, fermionic systems, quantum circuits, Givens rotations.

## Contents

<table><tr><td>1 Introduction</td><td>2</td></tr><tr><td>2 Determinantal and Pfaffian point processes</td><td>5</td></tr><tr><td>2.1 Determinantal point processes</td><td>5</td></tr><tr><td>2.2 Pfaffian point processes</td><td>7</td></tr></table>

<footnote>*Corresponding author: remi.bardenet@univ-lille.fr</footnote>
<table><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>3 From qubits to fermions</td><td>8</td></tr><tr><td>3.1 Models in quantum physics</td><td>8</td></tr><tr><td>3.1.1 An observable and a state define a random variable</td><td>8</td></tr><tr><td>3.1.2 Commuting observables and a state define a random vector</td><td>9</td></tr><tr><td>3.2 The canonical anti-commutation relations</td><td>9</td></tr><tr><td>3.3 Fermionic operators acting on qubits</td><td>10</td></tr><tr><td>4 From fermions to Pfaffian point processes</td><td>11</td></tr><tr><td>4.1 Building a DPP from a quantum measurement</td><td>11</td></tr><tr><td>4.2 The case of L-ensembles</td><td>12</td></tr><tr><td>4.3 The case of projection DPPs</td><td>14</td></tr><tr><td>4.4 Building a PfPP from the BdG Hamiltonian</td><td>15</td></tr><tr><td>5 Quantum circuits to sample DPPs and PfPPs</td><td>17</td></tr><tr><td>5.1 Quantum circuits</td><td>17</td></tr><tr><td>5.2 The case of a projection DPP</td><td>18</td></tr><tr><td>5.2.1 QR by Givens rotations</td><td>18</td></tr><tr><td>5.2.2 Parallel QR algorithms and qubit operation constraints</td><td>20</td></tr><tr><td>5.2.3 A hybrid parallel/quantum algorithm for projection DPPs</td><td>21</td></tr><tr><td>5.3 Reducing general DPPs to mixture DPPs</td><td>23</td></tr><tr><td>5.4 Dilating a general DPP to a projection DPP</td><td>24</td></tr><tr><td>5.5 Sampling the PfPP corresponding to the BdG Hamiltonian</td><td>24</td></tr><tr><td>5.5.1 QR decomposition with double Givens rotations</td><td>25</td></tr><tr><td>5.5.2 Quantum gates</td><td>28</td></tr><tr><td>6 Numerical experiments</td><td>28</td></tr><tr><td>6.1 Projection DPPs</td><td>28</td></tr><tr><td>6.2 Pfaffian PPs</td><td>32</td></tr><tr><td>7 Discussion</td><td>34</td></tr><tr><td>A Mathematical details</td><td>36</td></tr><tr><td>A.1 More about Givens operators</td><td>36</td></tr><tr><td>A.2 Proof of Lemma 4.7</td><td>36</td></tr><tr><td>A.3 Proof of Proposition 4.8</td><td>36</td></tr><tr><td>B Gate details</td><td>37</td></tr><tr><td>B.1 Givens rotation using a controlled unitary gate</td><td>38</td></tr><tr><td>B.2 Qiskit 0.42.1 implementation without controlled unitary gate</td><td>39</td></tr><tr><td>C Sources of error in quantum computers</td><td>39</td></tr></tbody></table>

# 1 Introduction

Determinantal point processes (DPPs) were introduced in the thesis of Macchi, 1972, recently translated and reprinted as (Macchi, 2017). Macchi's motivation was the design of probabilistic models for free fermions in quantum optics; see the preface of (Macchi, 2017) for a history of DPPs, and (Bardenet, Feller, et al., 2022) for an extended discussion of the links between free fermions and DPPs. DPPs have known another surge of interest since the 90s for their application to random matrix theory
(Johansson, 2005). More recently, they have been adopted as models or sampling tools in fields like spatial statistics (Lavancier, Møller, and Rubak, 2015), Monte Carlo methods (Bardenet and Hardy, 2020), machine learning (Kulesza and Taskar, 2012), or numerical linear algebra (Derezinski and Mahoney, 2021). In the latter two fields, the considered DPPs are often finite, in the sense that a DPP is a probability measure over subsets of a ground set of finite cardinality $N \gg 1$. Such a finite DPP is specified by an $N \times N$ matrix called its kernel matrix, which we assume here to be Hermitian.

In machine learning as in numerical linear algebra, it is crucial to be able to sample from the considered finite DPPs. For instance, a famous DPP is the subset of edges of a uniform spanning tree in a connected graph (Pemantle, 1991). Sampling these uniform spanning trees is a necessary step for building the randomized preconditioners of Laplacian systems in (Kyng and Song, 2018). As another example, DPPs have been used as randomized summaries of large collections of items, such as a corpus of texts. Sampling the corresponding DPP then allows to extract a small representative subset of sentences (Kulesza and Taskar, 2011, Section 4 and references therein). Other machine learning applications include constructing coresets (Tremblay, Barthelmé, and Amblard, 2019), kernel matrix approximation (Derezinski, Khanna, and Mahoney, 2020a; Fanuel, Schreurs, and Suykens, 2021) and feature extraction for linear regression (Belhadji, Bardenet, and Chainais, 2020).

Much research has thus been devoted to sampling finite DPPs on a classical computer, either exactly or approximately. The default generic exact sampler is the ‘HKPV’ sampler (Hough, Krishnapur, Peres, and Virág, 2006). To fix ideas, when applied to a *projection* DPP, i.e., a DPP that puts all its mass on subsets of some fixed cardinality $r \le N$, and assuming the kernel matrix is given in diagonalized form, HKPV has complexity $O(Nr^2)$. For DPPs on graphs such as uniform spanning trees, there also exist dedicated algorithms, such as the loop-erased random walks of Wilson (1996), with an expected number of steps equal to the mean commute time to a chosen root node of the graph.

Given that DPPs are originally a model in quantum electronic optics, and are still the default mathematical object used to describe a quantum physical system known as *free fermions* (Dean, Le Doussal, Majumdar, and Schehr, 2019), it is natural to ask whether finite DPPs can be sampled more efficiently on a quantum computer than on a classical computer. Somewhat implicitly, the question has actually already been tackled in a string of physics papers whose goal is the more ambitious quantum simulation of fermionic systems (Ortiz, Gubernatis, Knill, and Laflamme, 2001; Wecker, Hastings, Wiebe, Clark, Nayak, and Troyer, 2015; Kivlichan, McClean, Wiebe, Gidney, Aspuru-Guzik, Chan, and Babbush, 2018; Jiang, Sung, Kechedzhi, Smelyanskiy, and Boixo, 2018). In a reverse cross-disciplinary direction, and still rather implicitly, the quantum algorithms therein are reminiscent of parallel QR decompositions, a key topic in numerical algebra (Sameh and Kuck, 1978; Demmel, Grigori, Hoemmen, and Langou, 2012). While finishing this work, we also realized that in the computer science literature, and independently of the aforementioned physics works, Kerenidis and Prakash (2022) recently proposed similar quantum algorithms to sample projection DPPs, as a building block for quantum data analysis pipelines. For our purpose, their main contribution is a quantum circuit with depth logarithmic in $N$, when Jiang et al. (2018) only discuss depths linear in $N$.

On our side, motivated by applications of finite DPPs in data science, we initiated in (Bardenet, Feller, et al., 2022) a programme of reconnection of DPPs to their physical fermionic roots, to foster cross-disciplinary research between mathematics, computer science, and physics on the topic, even if our languages and lore are quite
different. In particular, physicists have developed tools for the analysis and the con-
struction of fermionic systems that we would like to apply to DPPs in data science,
without reinventing the wheel. The current paper is part of this programme, and sums
up our understanding of what the state of the art in physics tells us on sampling finite
DPPs, after we follow in the footsteps of (Macchi, 1972) and map a given DPP to a
fermionic density operator. This cross-disciplinary, self-contained survey is our first
contribution.

As an example of what our disciplines can bring to each other, our second con-
tribution is to relate the Pfaffian point processes as defined by Koshida (2021) – a
generalization of DPPs that is natural in theoretical physics but has not yet been
used in data science – to a quantum algorithm by Jiang et al. (2018) for solving
the Bogoliubov-de Gennes Hamiltonian. As another example of the fertility of cross-
disciplinary work, after we make the link between the quantum circuits of (Wecker
et al., 2015; Jiang et al., 2018) and parallel QR decompositions (Demmel et al., 2012),
many new variants of the quantum circuits in (Jiang et al., 2018) become immediately
available, adapting to a range of qubit-communication and hardware constraints. In
particular, we exhibit a variant of the quantum circuits in (Jiang et al., 2018) with
the same dimensions as the best circuit in (Kerenidis et al., 2022).

Overall, our conclusions on quantum DPP sampling are that if a projection kernel is given in diagonalized form **K** = **Q**\***Q**, with **Q** ∈ **C**<sup>r×N</sup> a matrix with orthonormal rows, one can build quantum circuits that sample DPP(**K**) with O(rN) one- and two-qubit gates, and depth depending on what hardware constraints we put on which qubits can be jointly operated. Acting only on neighbouring qubits, depth is O(N) (Wecker et al., 2015; Jiang et al., 2018), while acting on arbitrary pairs of qubits can take the depth down to O(r log N); see our variant of (Jiang et al., 2018) in Section 5.2.2 and the logarithmic depth Clifford loaders of Kerenidis et al., 2022. Such depths (i.e., the largest number of gates applied to any single qubit) favourably compare to the time complexity O(Nr<sup>2</sup>) of the classical HKPV algorithm, or the expected complexity in O(Nr + r<sup>3</sup> log r) of the randomized version of HKPV in (Dereziński, Clarkson, Mahoney, and Warmuth, 2019; Barthelmé, Tremblay, and Amblard, 2023).

That being said, diagonalizing **K** on a classical computer as a preprocessing step
remains a O(N<sup>3</sup>) bottleneck, or at least O(Nd<sup>2</sup>) in the common case where the diag-
onalization of **K** can be reduced to the SVD of an N × d matrix. This bottleneck thus
seems to cancel the advantage of using a quantum circuit if one insists on starting
with **K** stored on a classical computer. Yet, while the projection kernel **K** may not be
available in diagonalized form, it is common in data science applications (Kyng et al.,
2018; Belhadji et al., 2020) to specify it implicitly, as a set of vectors spanning its
range. As noted by Barthelmé et al., 2023, using a (classical) parallel QR algorithm
and *P* processors, we can reduce the classical preprocessing step to O(Nd<sup>2</sup>/P) flops.
Importantly for our quantum pipeline, we discuss here how to further reuse the in-
termediate steps of this preprocessing in the design of the quantum circuit to apply
next. This yields a hybrid parallel/quantum algorithm to sample projection DPPs.
Compared to the classical HKPV sampler, our pipeline thus provides a linear speedup.
Compared to the expected complexity of the randomized classical algorithm discussed
in (Dereziński et al., 2019; Barthelmé et al., 2023), we show a gain in the sampling
step, but we arguably share the same bottleneck of classical parallel QR preprocessing.
Finally, the necessity for classical preprocessing may disappear in the future, once **K**
can be assumed to be initially available as a quantum state, stored on a quantum
computer.

The rest of the paper is organized as follows. In Section 2, we define DPPs and
one of their generalizations, Pfaffian PPs (PfPPs). In Section 3, we introduce the vocabulary of quantum field theory, at the basis of the connection between PfPPs and free fermions. By sticking to the case of a finite-dimensional state space, we can avoid technical difficulties and provide a rigorous, stand-alone introduction, mostly following (Nielsen, 2005). Section 4 is devoted to building a Hamiltonian starting from a DPP or a PfPP, so that a simple quantum measurement yields a sample from the corresponding point process. In Section 5, we show how Wecker et al., 2015; Jiang et al., 2018 build circuits to simulate the fermionic systems corresponding to our point processes. Our presentation insists on the implicit links with parallel QR algorithms, which allow us to introduce variants of the circuits with smaller complexity under assumptions on the qubit communication constraints of the target machine. Finally, we investigate in Section 6 the implementation of the circuits with the library Qiskit (Qiskit contributors, 2023), and the noise when running the circuits on 5-qubit IBMQ machines (IBM Quantum, 2021).

Appendix A contains a few detailed proofs that we extracted from the main body of the paper. Appendix B contains a discussion on gate details to implement the basic operations introduced in the circuits of Section 5. Appendix C is an overview of the sources of error in current quantum computers and their orders of magnitude.

**Notations.** The complex conjugate of a complex number $z$ is denoted by $\bar{z}$. Similarly, $\bar{M}$ denotes the entrywise complex conjugate of a matrix $M$. The transpose of $M$ reads $M^\top$ and its Moore-Penrose pseudo-inverse is $M^+$. The adjoint of an operator $A$ is written $A^*$. Also, we denote the canonical basis elements of $\mathbb{C}^N$ by $e_k$, $1 \le k \le N$. The sigmoid function is $\sigma(x) = 1/(1 + \exp(-x))$. For $N \in \mathbb{N}^*$, let $[N] = \{1, \dots, N\}$. A point process on $[N]$ is a probability measure over subsets of $[N]$.

# 2 Determinantal and Pfaffian point processes

In this section, we introduce discrete determinantal point processes (DPPs), and refer to (Kulesza et al., 2012) for their elementary properties. We also introduce Pfaffian point processes (PfPPs, Rains, 2000; Soshnikov, 2003), a generalization of DPPs that has not yet been used in machine learning, to the best of our knowledge. As we shall see in Section 4.4, both DPPs and PfPPs naturally appear when modeling physical particles known as fermions.

## 2.1 Determinantal point processes

**Definition 2.1 (DPP).** Let $K \in \mathbb{C}^{N \times N}$. A random subset $Y \subseteq [N]$ is drawn from the DPP of marginal kernel $K$, denoted by $Y \sim \text{DPP}(K)$ if and only if

$$ \forall S \subseteq [N], \quad \mathbb{P}(S \subseteq Y) = \det(K_S), \qquad (1) $$

where $K_S = [K_{i,j}]_{i,j \in S}$. We take as convention $\det(K_\emptyset) = 1$.

A priori, it is not obvious that a given complex matrix $K$ defines a DPP.

**Theorem 2.2 (Macchi, 1972, Soshnikov, 2000).** When $K$ is Hermitian, existence of DPP($K$) is equivalent to the spectrum of $K$ being included in $[0, 1]$.

In particular, when the spectrum of $K$ is included in $\{0, 1\}$, we call $K$ a projection kernel, and the corresponding DPP a projection DPP. Letting $r$ be the number of unit
eigenvalues of its kernel, samples from a projection DPP have fixed cardinality $r$ with probability 1 (Hough et al., 2006, Lemma 17). In applications, projection kernels of rank $r$ are often available in one of two forms: either

$$ \mathbf{K} = \mathbf{A}(\mathbf{A}^*\mathbf{A})^+\mathbf{A}^*, \qquad (2) $$

where $\mathbf{A} \in \mathbb{C}^{N \times M}$ is any matrix with rank $r \le \min(N, M)$, or in diagonalized form

$$ \mathbf{K} = \mathbf{U}\mathbf{U}^*, \qquad (3) $$

where $\mathbf{U} \in \mathbb{C}^{N \times r}$ has orthonormal columns. We give an example application for each form.

**Example 2.3** (Uniform spanning trees). Consider a finite connected graph $G$ with $M$ vertices and $N$ edges, encoded by its oriented edge-vertex incidence matrix $\mathbf{A} \in \{-1, 0, 1\}^{N \times M}$. There are a finite number of spanning trees of $G$, and we draw one uniformly at random. The edges in that random tree correspond to a subset $Y$ of the indices $[N]$ of the rows of $\mathbf{A}$. It turns out (Pemantle, 1991) that $Y$ is a projection DPP with kernel (2).

Uniform spanning trees are useful in many contexts, e.g. to build preconditioners for certain linear systems (Kyng et al., 2018, Section 5). Another example of application of DPPs is column-subset selection.

**Example 2.4** (Column subset selection). Belhadji et al., 2020 propose to select $k$ columns of a “fat” matrix<fnref n="1" /> $\mathbf{X} \in \mathbb{R}^{n \times N}$, $N \gg n$, using the projection DPP with rank-k kernel

$$ \mathbf{K} = \mathbf{V}_{:,1:k} \mathbf{V}_{:,1:k}^{\top}, \qquad (4) $$

where $\mathbf{X} = \mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}$ is the singular value decomposition of $\mathbf{X}$. This is an example of DPP with a kernel specified by (3). Belhadji et al., 2020 prove that the projection of $\mathbf{X}$ onto the subspace spanned by the selected columns is essentially an optimal low-rank approximation of $\mathbf{X}$. This ensures statistical guarantees in sketched linear regression.

Because we assume that the kernel is Hermitian, a DPP can be seen as a repulsive distribution, in the sense that for all distinct $i, j \in [N]$,

$$ \begin{aligned} \mathbb{P}(\{i, j\} \subseteq Y) &= \mathbf{K}_{i,i} \mathbf{K}_{j,j} - \mathbf{K}_{i,j} \overline{\mathbf{K}}_{i,j} \\ &= \mathbb{P}(\{i\} \subseteq Y)\mathbb{P}(\{j\} \subseteq Y) - |\mathbf{K}_{i,j}|^2 \\ &\le \mathbb{P}(\{i\} \subseteq Y)\mathbb{P}(\{j\} \subseteq Y). \end{aligned} $$

This repulsiveness enforces diversity in samples, and is particularly adequate in applications where a DPP is used to extract a small diverse subset of a large collection of $N$ items. Beyond column subset selection, this diversity is natural in machine learning tasks such as summary extraction (Kulesza and Taskar, 2012) or experimental design (Derezinski, Liang, and Mahoney, 2020b; Poinas and Bardenet, 2023).

<footnote n="1">¹Note our different notation compared to (Belhadji et al., 2020), who use $N$ for the number of rows.</footnote>
## 2.2 Pfaffian point processes

Similarly to the determinant, the Pfaffian of a $2k \times 2k$ skew-symmetric matrix is a polynomial of the matrix entries

$$ \mathrm{Pf}(\mathbf{A}) = \sum_{\sigma \text{ contraction}} \mathrm{sgn}(\sigma) \mathbf{A}_{\sigma(1)\sigma(2)} \cdots \mathbf{A}_{\sigma(2k-1)\sigma(2k)}. $$

Recall that a contraction of order $m$ ($m$ even) is a permutation such that $\sigma(1) < \sigma(3) < \dots < \sigma(m-1)$, and $\sigma(2i-1) < \sigma(2i)$ for $i \le m/2$. To relate to determinants, note that the Pfaffian of a skew-symmetric matrix $\mathbf{A}$ of even size satisfies $\det \mathbf{A} = (\mathrm{Pf}\mathbf{A})^2$.

**Definition 2.5 (PfPP).** Let $\mathbb{K}: [N] \times [N] \to \mathbb{C}^{2 \times 2}$ satisfy $\mathbb{K}(i, j)^{\top} = -\mathbb{K}(j, i)$ for all $1 \le i, j \le N$. A random subset $Y \subseteq [N]$ is drawn from the PfPP of marginal kernel $\mathbb{K}$, denoted by $Y \sim \mathrm{PfPP}(\mathbb{K})$ if and only if

$$ \forall S \subseteq [N], \quad \mathbb{P}(S \subseteq Y) = \mathrm{Pf}(\mathbf{K}_S), \qquad (5) $$

where $\mathbf{K}_S = [\mathbb{K}(i, j)]_{i,j \in S}$ is a complex matrix made of $|S|$ blocks of size $2 \times 2$.

Sufficient conditions on $\mathbb{K}$ for the existence of PfPP($\mathbb{K}$) were given by Kargin (2014, Theorem 1.3) when $(\begin{smallmatrix} 0 & -1 \\ 1 & 0 \end{smallmatrix})\mathbb{K}(i, j)$ can be mapped to a self-adjoint quaternionic kernel taking values in the set of $2 \times 2$ complex matrices. Later Kassel (2015, Theorem 2.3) gave an equivalent of the Macchi-Soshnikov Theorem 2.2 for this type of processes; see (Kassel and Lévy, 2022, Theorem 7.6 and Proposition 7.11). This class of Pfaffian PPs was also studied by Bufetov, Deelan Cunden, and Qiu (2021) in the continuous setting.

More recently, Koshida (2021) gave another sufficient condition for the existence of a Pfaffian point process on a discrete ground set, which is well-suited to the processes considered in our paper. The PfPPs of Koshida (2021) correspond to the case where $(\begin{smallmatrix} 0 & 1 \\ 1 & 0 \end{smallmatrix})\mathbb{K}(i, j)$ is a self-adjoint complex kernel. The intersection of the classes of PfPPs studied by Kargin (2014) and Koshida (2021) is simply the set of PfPPs for which the $2 \times 2$ matrix $\mathbb{K}(i, j)$ has a vanishing diagonal, i.e., DPPs with Hermitian kernels; see Example 2.7 below.

Before going further, we introduce a few useful notations. Consider a $2N \times 2N$ complex matrix $\mathbf{S}$, viewed as made of four $N \times N$ blocks. Define the following transformation of $\mathbf{S}$, called here *particle-hole transformation*, which consists in taking the complex conjugation and exchanging blocks along diagonals, i.e.

$$ \mathrm{ph}(\mathbf{S}) = \mathbf{C}\bar{\mathbf{S}}\mathbf{C}, \text{ with } \mathbf{C} = \begin{pmatrix} \mathbf{0} & \mathbf{I} \\ \mathbf{I} & \mathbf{0} \end{pmatrix}. $$

**Proposition 2.6 (Koshida, 2021).** Let $\mathbf{S} = (\begin{smallmatrix} \mathbf{S}_{11} & \mathbf{S}_{12} \\ \mathbf{S}_{21} & \mathbf{S}_{22} \end{smallmatrix})$ be an Hermitian $2N \times 2N$ matrix such that $0 \preceq \mathbf{S} \preceq \mathbf{I}$ and $\mathrm{ph}(\mathbf{S}) = \mathbf{I} - \mathbf{S}$. There exists a Pfaffian point process with the marginal kernel

$$ \mathbb{K}(i, j) = \begin{pmatrix} \mathbf{S}_{21}(i, j) & \mathbf{S}_{22}(i, j) \\ \mathbf{S}_{11}(i, j) - \delta_{ij} & \mathbf{S}_{12}(i, j) \end{pmatrix}, \quad 1 \le i, j \le N. $$

A few remarks are in order. First, the properties of $\mathbf{S}$ allow to simplify the expression of the kernel in

$$ \mathbb{K}(i, j) = \begin{pmatrix} \mathbf{S}_{21}(i, j) & \mathbf{S}_{22}(i, j) \\ -\mathbf{S}_{22}(j, i) & \mathbf{S}_{21}(j, i) \end{pmatrix}, \qquad (6) $$

where $\mathbf{S}_{21}$ is skew-symmetric and $\mathbf{S}_{22}$ is Hermitian. Second, DPPs with Hermitian kernels appear as particular instances of the PfPPs of Proposition 2.6.
**Example 2.7** (vanishing diagonal). Let $\mathbf{S}$ satisfy the conditions of Proposition 2.6, and let $\mathbb{K}$ be the corresponding Pfaffian kernel. If $\mathbf{S}_{21} = 0$, $Y \sim \text{PfPP}(\mathbb{K})$ is distributed according to DPP($\mathbf{S}_{22}$).

Third, for $Y \sim \text{PfPP}(\mathbb{K})$ and $i \neq j$, the 2-point correlation function is

$$ \begin{aligned} \mathbb{P}(\{i, j\} \subseteq Y) &= \mathbf{S}_{22}(i, i)\mathbf{S}_{22}(j, j) - |\mathbf{S}_{22}(i, j)|^2 + |\mathbf{S}_{21}(i, j)|^2 \\ &= \mathbb{P}(\{i\} \subseteq Y)\mathbb{P}(\{j\} \subseteq Y) - |\mathbf{S}_{22}(i, j)|^2 + |\mathbf{S}_{21}(i, j)|^2. \end{aligned} \quad (7) $$

Compared with DPPs with Hermitian kernels, Equation (7) suggests that a Pfaffian point process as in Proposition 2.6 is less repulsive than the related determinantal process DPP($\mathbf{S}_{22}$) – an intuition for this fact is given in Section 4.4 below. Relatedly, note that the 1-point correlation functions of DPP($\mathbf{S}_{22}$) and PfPP($\mathbb{K}$) for $\mathbb{K}$ given in (6) are the same. In particular, the expected cardinality of $Y \sim \text{PfPP}(\mathbb{K})$ is simply $\mathbb{E}|Y| = \text{Tr}(\mathbf{S}_{22})$.

# 3 From qubits to fermions

The content of this section is standard; see e.g., the reference textbook (Nielsen and Chuang, 2010) for quantum computing basics and (Nielsen, 2005) for the Jordan-Wigner transform. We also refer to (Bardenet, Feller, et al., 2022), which presents all the basic elements required in this section in the context of optical measurements and the resulting point processes.

## 3.1 Models in quantum physics

A quantum model is given by (i) a Hilbert space $(\mathbb{H}, \langle \cdot | \cdot \rangle)$ called the state space, and (ii) a collection of self-adjoint operators $\mathbb{H} \to \mathbb{H}$ called observables, of which one particular observable $H: \mathbb{H} \to \mathbb{H}$ is singled out and called the Hamiltonian. Let $\psi \neq 0$ be an element of $\mathbb{H}$. All elements of the form $z\psi$ for a complex $z \neq 0$ represent the same quantum state, called a pure quantum state, as opposed to more general states to be defined later. To simplify expressions, it is conventional to only consider elements $\psi$ of unit norm, and to denote a unit-norm pure state by the “ket” $|\psi\rangle$, keeping in mind that, as long as $|z| = 1$, all vectors $z|\psi\rangle \in \mathbb{H}$ represent the same state. The corresponding “bra” $\langle\psi|$ is the linear form $|x\rangle \mapsto \langle\psi|x\rangle$.

### 3.1.1 An observable and a state define a random variable

Henceforth, we assume that $\mathbb{H}$ has finite dimension $d$. Take an observable $A$. By the spectral theorem, $A$ can be diagonalized in an orthonormal basis, say with eigenpairs $(\lambda_i, u_i)$ with $1 \le i \le d$. For simplicity, we momentarily assume that all the eigenvalues of $A$ have multiplicity 1. Together with a state $|\psi\rangle$, the observable $A = \sum_i \lambda_i u_i u_i^*$ describes a random variable $X_{A,\psi}$ on spec($A$) = $\{\lambda_1, \dots, \lambda_d\}$, through

$$ \mathbb{P}(X_{A,\psi} = \lambda_i) = |\langle\psi|u_i\rangle|^2. \quad (8) $$

When modeling statistical uncertainty on a state, like when describing the noisy output of an experimental device, physical states are not modelled as unit-norm vectors of $\mathbb{H}$, but rather as positive trace-one operators. To see how, we first map $|\psi\rangle$ to the
rank-one projector $\rho = |\psi\rangle\langle\psi|$. Then the distribution (8) can be equivalently defined as

$$ \mathbb{P}(X_{A,\rho} = \lambda_i) = \mathrm{tr} [\rho 1_{\{\lambda_i\}}(A)], \quad (9) $$

where for any $f: \mathbb{R} \to \mathbb{R}$, we have $f(A) = \sum_i f(\lambda_i) u_i u_i^*$. In particular, the expectation of $X_{A,|\psi\rangle\langle\psi|}$ is $\langle\psi|A|\psi\rangle$. Note that Definition (9) generalizes to operators A with eigenvalues with arbitrary multiplicity, and to states $\rho$ beyond projectors. In particular, for $\mu$ a probability measure on $\mathbb{H}$,

$$ \rho = \mathbb{E}_{\psi \sim \mu} |\psi\rangle\langle\psi| \quad (10) $$

still defines a probability measure on the spectrum of A through (9). The expectation of that distribution, also called the *expectation value* of operator A, is $\langle A \rangle_{\rho} \triangleq \mathrm{tr} \, \rho A$. In physics, the association (9) of a state-observable pair $(\rho, A)$ with the random variable $X_{A,\rho}$ is known as *Born's rule*.

Any $\rho$ that is not a rank-one projector is called a *mixed state*, by opposition to rank-one projectors, which are *pure states*. Mixed states like (10) are commonly used to describe any uncertainty in the actual state of an experimental entity.<fnref n="2" />

### 3.1.2 Commuting observables and a state define a random vector

When all pairs of a set of observables $A_1, \dots, A_p$ commute, then these observables can be diagonalized in the same orthonormal basis $(u_i)$, and (9) can be naturally generalized to describe a random vector $(X_{A_j,\rho})$ of dimension $p$, with values in the Cartesian product

$$ \mathrm{spec}(A_1) \times \dots \times \mathrm{spec}(A_p). $$

More precisely, the law of $(X_{A_j,\rho})$ is given by

$$ \mathbb{P}(X_{A_1, \rho} = x_1, \dots, X_{A_p, \rho} = x_p) = \mathrm{tr} [\rho 1_{\{x_1\}}(A_1) \dots 1_{\{x_p\}}(A_p)], \quad (11) $$

see e.g. (Bouten, Van Handel, and James, 2007) for an introduction. In this paper, we will associate a Pfaffian point process to a particular mixed state and a set of commuting observables, which can respectively be efficiently prepared and measured on a quantum computer. To define these objects, we first need to explain how physicists build Hamiltonians of fermionic systems.

## 3.2 The canonical anti-commutation relations

The Hamiltonian and its structure are often the key part in specifying a model, much like the factorization of a joint distribution in a probabilistic model. In the case of fermions, a family of physical particles that includes electrons, Hamiltonians are typically built as polynomials of fermionic creation-annihilation operators, i.e., operators that satisfy the so-called *canonical anti-commutation relations* (CAR).

**Definition 3.1 (CAR).** Let $\mathbb{H}$ be a Hilbert space. The operators $c_j : \mathbb{H} \to \mathbb{H}$, $j = 1, \dots, p$ and their adjoints are said to satisfy the canonical anti-commutation relations if

$$ \{c_i, c_j\} = \{c_i^*, c_j^*\} = 0 \quad \text{and} \quad \{c_i, c_j^*\} = \delta_{ij} \mathbb{I}, \quad (\text{CAR}) $$

where $\{u, v\} := uv + vu$ is the anti-commutator of operators u and v.

<footnote n="2">²Statisticians would say *epistemic* uncertainty, i.e., imprecise knowledge of the state.</footnote>
Assuming existence<fnref n="3" /> for a moment, and limiting ourselves to a finite dimensional Hilbert space $\mathbb{H}$ of dimension $d = 2^N$, one can say many things on $\mathbb{H}$ from the fact that there are $N$ operators $c_1, \dots, c_N$ satisfying (CAR). On that topic, we recommend reading (Nielsen, 2005), from which we borrow the following lemma.

**Lemma 3.2** (Fock basis; see e.g. Nielsen, 2005). Let $\dim\mathbb{H} = 2^N$, $N \in \mathbb{N}^*$, and assume that $c_1, \dots, c_N$ are distinct operators on $\mathbb{H}$ that satisfy (CAR). First, there is a vector $|\emptyset\rangle \in \mathbb{H}$, called the vacuum, which is a simultaneous eigenvector of all $c_i^* c_i$, $i = 1, \dots, N$, always with eigenvalue 0. Second, for $\mathbf{n} = (n_1, \dots, n_N) \in \{0, 1\}^N$, consider

$$|\mathbf{n}\rangle := \prod_{i=1}^{N} (c_i^*)^{n_i} |\emptyset\rangle.$$

Then $\mathcal{B}_{\text{Fock}} = (|\mathbf{n}\rangle)$ is an orthonormal basis of $\mathbb{H}$. Third, for all $1 \le i \le N$,

$$c_i^* c_i |\mathbf{n}\rangle = n_i |\mathbf{n}\rangle, \qquad (12)$$

and, for $i \ne j$,

$$c_i^* c_j |\mathbf{n}\rangle = \pm n_j (1 - n_i) |\tilde{\mathbf{n}}\rangle, \qquad (13)$$

where $\tilde{n}_i = 1$, $\tilde{n}_j = 0$, and $\tilde{n}_k = n_k$ for $k \ne i, j$.

The basis built in the lemma is called the Fock basis. Its construction depends on the choice of the operators $c_1, \dots, c_N$. When there is a risk of confusion, we shall thus further denote $|\mathbf{n}\rangle$ and $|\emptyset\rangle$ as $|\mathbf{n}_c\rangle$ and $|\emptyset_c\rangle$, respectively. Second, because applying $c_i$ to a vector of the Fock basis has the effect of zeroing the $i$th component if it was 1, and mapping to zero otherwise, we call the $c_i$'s annihilation operators. Similarly, we call their adjoints creation operators.

## 3.3 Fermionic operators acting on qubits

Henceforth, we let $\mathbb{H} := (\mathbb{C}^2)^{\otimes N}$. This is the state space describing $N$ qubits, of dimension $2^N$. A qubit corresponds to any physical system, the state of which is described by one out of two levels. We associate these two levels with a distinguished orthonormal basis ($|0\rangle, |1\rangle$) of $\mathbb{C}^2$, commonly called the computational basis.

Consider the so-called Pauli operators on $\mathbb{C}^2$, given in the computational basis as

$$\mathbf{I} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}, \quad \sigma_x = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}, \quad \sigma_y = \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix}, \quad \sigma_z = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}. \qquad (14)$$

Note that $\sigma_x^2 = \sigma_y^2 = \sigma_z^2 = \mathbf{I}$, and that $\sigma_x \sigma_y = i\sigma_z$.

**Definition 3.3** (The Jordan-Wigner (JW) transformation). Define the JW annihilation operators $a_j$ on $\mathbb{H}$, for $j \in \{1, \dots, N\}$, by their matrices in the computational basis

$$\underbrace{\sigma_z \otimes \cdots \otimes \sigma_z}_{j-1 \text{ times}} \otimes \left( \frac{\sigma_x + i\sigma_y}{2} \right) \otimes \mathbf{I} \otimes \cdots \otimes \mathbf{I}. \qquad (15)$$

The creation operator $a_j^*$, adjoint of $a_j$, has matrix

$$\underbrace{\sigma_z \otimes \cdots \otimes \sigma_z}_{j-1 \text{ times}} \otimes \left( \frac{\sigma_x - i\sigma_y}{2} \right) \otimes \mathbf{I} \otimes \cdots \otimes \mathbf{I}.$$

<footnote n="3">³Definition 3.3 gives a realization of (CAR); see Proposition 3.4.</footnote>
It is easy to check the following properties of the Jordan-Wigner operators.

**Proposition 3.4.** Let $a_1, \dots, a_N$ be the JW operators from Definition 3.3. They satisfy (CAR). Moreover, the so-called number operators $N_j = a_j^* a_j$ have the simple matrix form

$$ \underbrace{\mathbf{I} \otimes \dots \otimes \mathbf{I}}_{j-1 \text{ times}} \otimes \frac{\mathbf{I} - \sigma_z}{2} \otimes \mathbf{I} \otimes \dots \otimes \mathbf{I}. $$

Finally, the Fock basis $|n\rangle$ obtained by setting $c_j = a_j$ in Lemma 3.2 is actually the computational basis, i.e.,

$$ |n\rangle = |n_1\rangle \otimes \dots \otimes |n_N\rangle. $$

There are alternative ways to define fermionic operators on $\mathbb{H}$ using Pauli matrices, like the parity or Bravyi-Kitaev transformations (Bravyi and Kitaev, 2002); or the Ball-Verstraete-Cirac transformation (Ball, 2005; Verstraete and Cirac, 2005). In particular, and at the expense of simplicity, the Bravyi-Kitaev transformation yields operators with smaller weight than Jordan-Wigner: only $O(\log N)$ terms are allowed to differ from the identity in the equivalent of (15).

# 4 From fermions to Pfaffian point processes

We first show in Section 4.1 how to build any discrete DPP with Hermitian kernel from a mixed state corresponding to a particle number preserving Hamiltonian and a set of commuting observables. This connection between DPPs and quasi-free states was observed by Olshanski (2020). Then we do the same for a class of discrete PfPPs in Section 4.4 associated to Hamiltonians without particle number conservation.

## 4.1 Building a DPP from a quantum measurement

Let $\mathbf{H} \in \mathbb{C}^{N \times N}$ be Hermitian, and consider the operator

$$ H = \sum_{i,j=1}^{N} c_i^* \mathbf{H}_{ij} c_j, \qquad (16) $$

which we think of as a Hamiltonian acting on $\mathbb{H}$. This Hamiltonian preserves the particle number, i.e., $\mathbf{H}$ commutes with the number operator $\sum_i c_i^* c_i$. Let $(\lambda_k, \mathbf{u}_k)_{1 \le k \le N}$ be the eigenpairs of $\mathbf{H}$ where the eigenvalues satisfy $\lambda_1 \le \dots \le \lambda_N$. The Hamiltonian (16) can thus be rewritten in “diagonalized” form

$$ H = \sum_{i,j=1}^{N} c_i^* \left( \sum_{k=1}^{N} \lambda_k \mathbf{u}_k \mathbf{u}_k^* \right)_{ij} c_j = \sum_{i,j,k=1}^{N} \lambda_k c_i^* \mathbf{u}_{ki} \mathbf{u}_{kj}^* c_j = \sum_{k=1}^{N} \lambda_k b_k^* b_k, $$

where we defined the operators

$$ b_k = \sum_{j=1}^{N} \mathbf{u}_{kj}^* c_j, \qquad (17) $$

which can be checked to satisfy (CAR).

In what follows, we consider in Section 4.2 DPPs for which the correlation kernel has a spectrum in $[0, 1)$, referred to as L-ensembles in the literature. Next, in Section 4.3, we discuss projection DPPs.
## 4.2 The case of L-ensembles

Below, Proposition 4.1 states that there is a natural DPP associated with a Hamiltonian like (16), whose marginal kernel $K$ is obtained by applying a sigmoid to the spectrum of $H$; formally $K = \sigma(-\beta\bar{H})$ with $\beta > 0$ being a parameter called inverse temperature.

**Proposition 4.1** (DPP kernel by taking the sigmoid). Let $\beta > 0$ and $\mu \ge 0$, $H$ and $(b_k)$ be respectively defined by (16) and (17). Consider the mixed state

$$ \rho = \frac{1}{Z} e^{-\beta(H - \mu \sum_{j=1}^{N} b_j^* b_j)}, \qquad (18) $$

where the normalization constant $Z$ ensures that $\text{tr} \rho = 1$. For $i \in [N]$, the observable $N_i = c_i^* c_i$ is a projector, so that the random variable $X_{N_i, \rho}$ associated to $N_i$ and the state (18) by (9) takes values in $\{0, 1\}$. Moreover, all $N_i$'s commute pairwise, and thus define a joint Boolean vector $(X_{N_i, \rho})_{i \in [N]}$ through (11). Consider the point process

$$ S = \{i \in [N] : X_{N_i, \rho} = 1\}, $$

corresponding to jointly observing all operators $N_i$ in the state (18), and reporting the indices corresponding to a 1. Then $S$ is determinantal with correlation kernel

$$ K = \bar{U} \text{Diag} \left( \frac{e^{-\beta(\lambda_k - \mu)}}{1 + e^{-\beta(\lambda_k - \mu)}} \right) \bar{U}^*, $$

where $U$ and $(\lambda_k)_k$ are obtained by the diagonalization $H = U \text{Diag} (\lambda_k) U^*.$

*Proof.* All number operators $N_i = c_i^* c_i$ commute pairwise and have spectrum in $\{0, 1\}$; see Lemma 3.2. Consequently, their joint measurement indeed describes a random binary vector $X = (X_{N_i, \rho})$. By Born's rule (11), the correlation function of the point process corresponding to the indices of the 1s in $X$ are given by

$$ \mathbb{P}(\{i_1, \dots, i_k\} \subseteq X) = \text{tr}(\rho N_{i_1} \dots N_{i_k}) \triangleq \langle N_{i_1} \dots N_{i_k} \rangle_{\rho}. $$

Because of the anti-commutation relations (CAR), an explicit computation known as Wick's theorem implies that for any $i_1, \dots, i_k \in [N]$,

$$ \langle N_{i_1} \dots N_{i_k} \rangle_{\rho} = \det \left( \langle c_{i_m}^* c_{i_n} \rangle_{\rho} \right)_{m,n}. \qquad (19) $$

Wick's theorem is a standard result in quantum physics; see e.g. (Bardenet et al., 2022, Section 3.5) for a rewriting with our notations of one of the canonical references (Cohen-Tannoudji, Diu, and Laloë, 1977).

Now, Equation (19) implies that the point process $X$, consisting of simultaneously measuring the occupation of all qubits using $(N_i)$, is determinantal with correlation kernel

$$ K_{ij} = \langle c_i^* c_j \rangle_{\rho}. \qquad (20) $$

In order to provide an explicit computation of $K_{ij}$, we use Lemma 3.2 with the operators $b_1, \dots, b_N$, to obtain a basis ($|n\rangle$) = ($|n\rangle_b$) of eigenvectors of all operators of the
form $b_i^*b_i$, $i = 1, \dots, N$. We now proceed to computing the trace in (20) by summing over that basis. We write

$$
\begin{align*}
\mathbf{K}_{ij} &= \operatorname{tr}[\rho c_i^* c_j] \\
&= \operatorname{tr}\left[\rho c_i^* c_j \sum_{\mathbf{n}} |\mathbf{n}\rangle\langle\mathbf{n}|\right] \\
&= Z^{-1} \sum_{\mathbf{n}} \langle\mathbf{n}|e^{-\beta \sum_{p}(\lambda_p - \mu)b_p^* b_p} c_i^* c_j |\mathbf{n}\rangle.
\end{align*}
$$

Now note that $c_j = \sum_{\ell=1}^{N} \mathbf{u}_{j\ell} b_{\ell}$, so that $c_i^* c_j = \sum_{k,\ell=1}^{N} \mathbf{u}_{ik}^* \mathbf{u}_{j\ell} b_k^* b_{\ell}$, and

$$
\mathbf{K}_{ij} = Z^{-1} \sum_{k,l=1}^{N} \mathbf{u}_{ik}^{*} \mathbf{u}_{j\ell} \sum_{\mathbf{n}} \langle \mathbf{n} | e^{-\beta \sum_{p} (\lambda_{p} - \mu) b_{p}^{*} b_{p}} b_{k}^{*} b_{\ell} | \mathbf{n} \rangle .
$$

Because of (13) and the orthonormality of the basis ($|\mathbf{n}\rangle$), for $k \neq \ell$,

$$
\langle \mathbf{n} | e^{-\beta \sum_p (\lambda_p - \mu) b_p^* b_p} b_k^* b_\ell | \mathbf{n} \rangle = 0.
$$

Therefore,

$$
\mathbf{K}_{ij} = Z^{-1} \sum_{k=1}^{N} \mathbf{u}_{ik}^{*} \mathbf{u}_{jk} \sum_{\mathbf{n}} \langle \mathbf{n} | e^{-\beta \sum_{p} (\lambda_p - \mu) b_p^* b_p} b_k^* b_k | \mathbf{n} \rangle \quad (21)
$$

$$
= Z^{-1} \sum_{k=1}^{N} \mathbf{u}_{ik}^{*} \mathbf{u}_{jk} \sum_{\mathbf{n}} n_k e^{-\beta \sum_{p} (\lambda_p - \mu) n_p}. \quad (22)
$$

Now, we pause and compute the normalization constant $Z$ of $\rho$. Starting from $\operatorname{tr} \rho = 1$, we write

$$
Z = \operatorname{tr} e^{-\beta \sum_{p} (\lambda_p - \mu) b_p^* b_p} = \sum_{\mathbf{n}} \langle \mathbf{n} | e^{-\beta \sum_{p} (\lambda_p - \mu) b_p^* b_p} | \mathbf{n} \rangle = \sum_{\mathbf{n}} e^{-\beta \sum_{p} (\lambda_p - \mu) n_p}.
$$

Rewriting the sum as a product, we obtain

$$
Z = \prod_{p=1}^{N} \sum_{n_p \in \{0,1\}} e^{-\beta(\lambda_p - \mu)n_p} = \prod_{p=1}^{N} (1 + e^{-\beta(\lambda_p - \mu)}). \quad (23)
$$

Now, explicitly writing the dependence of $Z = Z(\lambda)$ to $\lambda = (\lambda_1, \dots, \lambda_N)$, Equations (22) and (23) together yield

$$
\mathbf{K}_{ij} = \sum_{k=1}^{N} \mathbf{u}_{ik}^{*} \mathbf{u}_{jk} \frac{\sum_{\mathbf{n}} n_k e^{-\beta \sum_{p} (\lambda_p - \mu) n_p}}{\sum_{\mathbf{n}} e^{-\beta \sum_{p} (\lambda_p - \mu) n_p}} \quad (24)
$$

$$
= \sum_{k=1}^{N} \mathbf{u}_{ik}^{*} \mathbf{u}_{jk} \left[ -\frac{1}{\beta} \frac{d}{d\lambda_k} \log Z(\lambda) \right] \quad (25)
$$

$$
= \sum_{k=1}^{N} \mathbf{u}_{ik}^{*} \mathbf{u}_{jk} \frac{e^{-\beta(\lambda_k - \mu)}}{1 + e^{-\beta(\lambda_k - \mu)}}, \quad (26)
$$

for any $1 \le i, j \le N$.

□
**Corollary 4.2** (Hamiltonian by taking the logit). Let $\mathbf{K} = \mathbf{VDV}^*$ be an Hermitian DPP kernel with $0 \prec \mathbf{D} \prec \mathbf{I}$. Then Proposition 4.1 yields a DPP with kernel $\mathbf{K}$ provided we choose $\mathbf{V} = \bar{\mathbf{U}}$ and the eigenvalues $\lambda$ so that

$$ \beta(\lambda_i - \mu) = \log \frac{1 - d_i}{d_i}. $$

Furthermore, assuming $d_1 \ge \dots \ge d_N \ge 0$, and $d_r < \mu < d_{r+1}$, $\mathbf{K}$ converges to the projection kernel onto the first $r$ columns of $\mathbf{V}$.

*Proof.* The first part is a consequence of (21). The second statement is maybe easiest to see in terms of Frobenius norm, to wit $\|\mathbf{A}\|^2 = \sum_{i=1}^N \sigma_i^2(\mathbf{A})$, where $\sigma_i(\mathbf{A})$ are the singular values of $\mathbf{A}$. If $\mathbf{P}$ denotes the projector into the first $r$ columns of $\mathbf{V}$, then the absolute values of the eigenvalues of the Hermitian matrix $\mathbf{K} - \mathbf{P}$ are its singular values, and all of them converge to 0. $\square$

## 4.3 The case of projection DPPs

With the notation of Section 4.2, we note that $(|\mathbf{n}\rangle) = (|\mathbf{n}\rangle_b)$ is a basis of eigenvectors of $\rho$, and that

$$ \rho |\mathbf{n}\rangle = \frac{e^{-\beta \sum_{p} (\lambda_p - \mu) n_p}}{\sum_{\mathbf{n}'} e^{-\beta \sum_{p} (\lambda_p - \mu) n_p'}} |\mathbf{n}\rangle. $$

In particular, the only eigenvalue that does not vanish when $\beta \to \infty$ is that of $|1, \dots, 1, 0, \dots, 0\rangle$, where the 1s occupy the first $r$ components, and $r \in \{1, \dots, N\}$ is such that with $\lambda_r < \mu < \lambda_{r+1}$. Indeed, the ratio of the eigenvalue of $|1, \dots, 1, 0, \dots, 0\rangle$ with that of any other eigenvector diverges to $+\infty$, while all eigenvalues remain in $[0, 1]$ and the trace is fixed to 1. Thus, denoting $|\psi\rangle = |1, \dots, 1, 0, \dots, 0\rangle$,

$$ \rho \to |\psi\rangle\langle\psi| $$

in Frobenius norm as $\beta \to \infty$. In particular,

$$ \mathrm{tr} N_{i_1} \dots N_{i_k} \rho \to \mathrm{tr} N_{i_1} \dots N_{i_k} |\psi\rangle\langle\psi| $$

as $\beta \to \infty$. Combined with the second statement of Corollary 4.2, we know that the correlation functions of the point process corresponding to preparing $|\psi\rangle\langle\psi|$ and measuring the occupation of all qubits and recording where the 1s occur is the projection DPP of kernel $\mathbf{P} = \mathbf{V}_{:,1:r} \mathbf{V}_{:,1:r}^*$.

**Remark 4.3** (Beyond Gaussian states). We have shown that the kernels of projection DPPs are obtained as limits of kernels associated with Gaussian states (18). Actually, a DPP kernel can be associated with a density matrix of a quasi-free state generalizing the Gaussian state, for which Wick's theorem also holds; see e.g. (Koshida, 2021, Section 1.3). In particular, every quasi-free state is a convex combination of pure states. We give now a few details by following (Dierckx, Fannes, and Pogorzelska, 2008). Let $\mathcal{C} \subseteq [N]$ and let $|\psi_C\rangle = \prod_{i \in \mathcal{C}} c_i^* |0\rangle$. Let $\mathbf{K}$ be an Hermitian matrix with eigenvalues $(\nu_p)$ in $[0, 1]$. The density matrix associated with $\mathbf{K}$ is the quasi-free state

$$ \rho = \sum_{\mathcal{C} \subseteq [N]} \alpha_{\mathcal{C}} |\psi_{\mathcal{C}}\rangle\langle\psi_{\mathcal{C}}| \quad \text{with} \quad \alpha_{\mathcal{C}} = \prod_{p \in \mathcal{C}} \nu_p \prod_{q \in [N] \setminus \mathcal{C}} (1 - \nu_q). $$
## 4.4 Building a PfPP from the BdG Hamiltonian

We shall see that Pfaffian point processes appear using slightly more sophisticated Hamiltonians. Unlike Hamiltonian (16), the Hermitian operator

$$H = \frac{1}{2} \sum_{i,j=1}^{N} \mathbf{M}_{ij} (c_i^* c_j - c_j^* c_i) + \frac{1}{2} \sum_{i,j=1}^{N} (\Delta_{ij} c_i^* c_j^* + \Delta_{ij}^* c_i c_j), \quad (27)$$

with complex matrices $\mathbf{M} = \mathbf{M}^*$ and $\Delta$, does not commute with the total number operator $N = \sum_i c_i^* c_i$. Physicists say that $H$ does not “preserve” the total number of particles.

Note that, due to the anti-commutation relation $c_i^* c_j^* = -c_j^* c_i^*$, we can simply redefine $\Delta$ so that $\Delta = -\Delta^\dagger$. The Hamiltonian (27) becomes the so-called Bogoliubov-de Gennes (BdG) Hamiltonian

$$H_{\text{BdG}} = \frac{1}{2} \sum_{i,j=1}^{N} \mathbf{M}_{ij} (c_i^* c_j - c_j^* c_i) + \frac{1}{2} \sum_{i,j=1}^{N} (\Delta_{ij} c_i^* c_j^* - \overline{\Delta}_{ij} c_i c_j), \quad (28)$$

which is a model for superconductors; see e.g. (Schnyder, Ryu, Furusaki, and Ludwig, 2008) for a modern overview.

We now investigate the point process associated to the occupation numbers of the Gaussian state $Z^{-1} \exp(-\beta H_{\text{BdG}})$, for $\beta > 0$, as we did for (18). It is convenient to stack the operators $c_1, \dots, c_N$ and their adjoints in column vectors, and write $\mathbf{c} = [c_1 \dots c_N]^\top$ and $\mathbf{c}^* = [c_1^* \dots c_N^*]^\top$. We can then rewrite the Hamiltonian more compactly as

$$H_{\text{BdG}} = \frac{1}{2} (c_c^* c_c)^* \mathbf{H}_{\text{BdG}} (c_c^* c_c) \quad \text{with} \quad \mathbf{H}_{\text{BdG}} = \begin{pmatrix} -\bar{\mathbf{M}} & -\bar{\Delta} \\ \bar{\Delta} & \mathbf{M} \end{pmatrix}. \quad (29)$$

By construction, the matrix $\mathbf{H}_{\text{BdG}}$ obeys the particle-hole symmetry, namely $\overline{\mathbf{C}\mathbf{H}_{\text{BdG}}\mathbf{C}} = -\mathbf{H}_{\text{BdG}}$ where we recall the definition of the involution $\mathbf{C} = \begin{pmatrix} 0 & \mathbf{I} \\ \mathbf{I} & 0 \end{pmatrix}$. The name particle-hole symmetry comes from the fact that

$$\mathbf{C}(c_c^*) = (c_c^*)$$

flips the role of creation and annihilation operators.

Before going further, we introduce a group of transformations preserving (CAR) that are used to diagonalize the Bogoliubov-de Gennes Hamiltonian; see (Moore, 2014, Section 18.4.3).

**Definition 4.4** (orthogonal complex matrix). A complex invertible $2N \times 2N$ matrix $\mathbf{W}$ satisfying $\mathbf{W}^\top \mathbf{C} \mathbf{W} = \mathbf{C}$ is called an orthogonal complex matrix. The group of orthogonal complex matrices is denoted by $O(\mathbf{C}, \mathbb{C})$.

For any orthogonal complex matrix $\mathbf{W}$ and any $N$ operators $c_1, \dots, c_N$ satisfying (CAR), another set of creation-annihilation operators satisfying (CAR) is given by the so-called *Bogoliubov transformation*

$$(b_b^*) = W(c_c^*), \quad (30)$$

and by further requiring that $\mathbf{W}$ is unitary, $b_k^*$ is the adjoint of $b_k$ for all $k \in [N]$. Hence, in what follows, we only consider transformations $\mathbf{W} \in O(\mathbf{C}, \mathbb{C}) \cap U(2N)$.
**Lemma 4.5.** A unitary orthogonal complex matrix is of the form $(\frac{U}{V}\frac{V}{U})$, with $UU^* + VV^* = I$ and $UV^\top + VU^\top = 0$.

Next, following Jiang et al. (2018), we diagonalize (28) by a convenient change of variables, which consists in finding a suitable Bogoliubov transformation. The upshot is that we have the decomposition

$$H_{\text{BdG}} = W^* \begin{pmatrix} \text{Diag}(-\epsilon_k) & 0 \\ 0 & \text{Diag}(\epsilon_k) \end{pmatrix} W,$$

as shown in Lemma 4.6 below.

**Lemma 4.6 (Diagonalization of BdG Hamiltonian).** Let $\Omega = \frac{1}{\sqrt{2}}(I_i I_i)$ and define $A = -i\Omega(\frac{\Delta}{-M} - \frac{M}{-\Delta})\Omega^*$. The following statements hold.

(i) $H_{\text{BdG}} = \frac{1}{2}f^\top A f$ where $f = \Omega(\begin{pmatrix} c^* \\ c \end{pmatrix})$ and $A$ is real skew-symmetric.

(ii) There exists a real orthogonal matrix $R$ and a real vector $\epsilon = [\epsilon_1, \dots, \epsilon_N]^\top$ such that $RAR^\top = (\begin{pmatrix} 0 & \text{Diag}(\epsilon_k) \\ -\text{Diag}(\epsilon_k) & 0 \end{pmatrix})$.

(iii) Another set of creation-annihilation operators satisfying (CAR) is given by $(\begin{pmatrix} b^* \\ b \end{pmatrix}) = W(\begin{pmatrix} c^* \\ c \end{pmatrix})$, where $W = \Omega^* R \Omega \in O(C, C) \cap U(2N)$.

(iv) We have the diagonalization $H_{\text{BdG}} = \frac{1}{2}(\begin{pmatrix} b^* \\ b \end{pmatrix})^\top (\begin{pmatrix} 0 & \text{Diag}(\epsilon_k) \\ -\text{Diag}(\epsilon_k) & 0 \end{pmatrix}) (\begin{pmatrix} b^* \\ b \end{pmatrix})$.

We refer to (Jiang et al., 2018) for a proof sketch.

Now, we leverage these results to compute the expectation value of bilinears under $\rho_{\text{BdG}}$. Lemma 4.7 states a result analogous to Proposition 4.1, and its proof, given in Appendix A.2, relies on similar techniques such as Wick's theorem.

**Lemma 4.7.** Let $\epsilon = [\epsilon_1, \dots, \epsilon_N]^\top$ be the eigenvalues of $H_{\text{BdG}}$ as given by Lemma 4.6 and let $\rho = \frac{1}{2}\exp(-\beta H_{\text{BdG}})$. We have

$$S \triangleq \begin{pmatrix} (\langle c_i c_j^* \rangle_\rho)_{i,j} & (\langle c_i c_j \rangle_\rho)_{i,j} \\ (\langle c_i^* c_j^* \rangle_\rho)_{i,j} & (\langle c_i^* c_j \rangle_\rho)_{i,j} \end{pmatrix} = \overline{W}^* \begin{pmatrix} \text{Diag}(\sigma(\beta \epsilon_k)) & \mathbf{0} \\ \mathbf{0} & \text{Diag}(\sigma(-\beta \epsilon_k)) \end{pmatrix} \overline{W},$$

where $\sigma(x)$ is the sigmoid function; see Section 1. Furthermore, we have that $S = \sigma(-\beta \overline{H}_{\text{BdG}})$ and $S$ satisfies $0 \preceq S \preceq I$, as well as $CSC = I - S$.

For convenience, we introduce the following notations, for $1 \le i, j \le N$,

$$\begin{pmatrix} \langle c_i c_j^* \rangle_\rho & \langle c_i c_j \rangle_\rho \\ \langle c_i^* c_j^* \rangle_\rho & \langle c_i^* c_j \rangle_\rho \end{pmatrix} \triangleq \begin{pmatrix} \delta_{ij} - \mathbf{K}_{ij}^T & -\overline{\mathbf{P}}_{ij} \\ \overline{\mathbf{P}}_{ij} & \mathbf{K}_{ij} \end{pmatrix}.$$

where $0 \preceq \mathbf{K} \preceq I$ is Hermitian and $\mathbf{P}$ is skew-symmetric. By definition, we have

$$S = \begin{pmatrix} I - \overline{\mathbf{K}} & \mathbf{P}^* \\ \mathbf{P} & \mathbf{K} \end{pmatrix}.$$

The particle-hole transformation amounts to replacing in $S$ each $c_i$ by $c_i^*$ and vice-versa. As a consequence of (CAR), $S$ satisfies $CSC = I - S$; see Section 2.2.

Proposition 4.8 can be seen as a generalization to mixed states of the analysis of Terhal and DiVincenzo (2002), and restates the results of Koshida (2021).
**Proposition 4.8.** Let $N_i = c_i^* c_i$ for $1 \le i \le N$ and let $\rho = Z^{-1} \exp(-\beta H_{BdG})$. For any $i_1, \dots, i_k \in [N]$, we have

$$ \langle N_{i_1} \dots N_{i_k} \rangle_{\rho} = \mathrm{Pf}(\mathbb{K}(i_m, i_n))_{1 \le m, n \le k}, \qquad (31) $$

where each block of the above $2k \times 2k$ matrix is given in terms of the following kernel valued in $2 \times 2$ matrices:

$$ \mathbb{K}(i, j) = \begin{pmatrix} \mathbf{P}_{ij} & \mathbf{K}_{ij} \\ -\mathbf{K}_{ji} & -\mathbf{P}_{ij} \end{pmatrix}. \qquad (32) $$

The latter satisfies $\mathbb{K}(i, j)^{\top} = -\mathbb{K}(j, i)$ for $1 \le i, j \le N$.

The proof of this result, given in Appendix A.3, is also based on Wick's theorem.

# 5 Quantum circuits to sample DPPs and PfPPs

Armed with the connections between PfPPs and fermions in Section 4, it remains to connect fermions and quantum circuits. Quantum circuits are briefly introduced in Section 5.1 for self-containedness. In Section 5.2, we describe the quantum circuit of Wecker et al., 2015, later modified by Jiang et al. (2018), that corresponds to a projection DPP given the diagonalized form of its kernel. In Section 5.2.2, we depart from (Jiang et al., 2018) and highlight the connections of their construction with a classical parallel QR algorithm in numerical algebra (Sameh et al., 1978). Consequently, we propose to take inspiration from more recent parallel QR algorithms, such as (Demmel et al., 2012), to construct circuits with different constraints on the communication between qubits. In particular, we recover circuits with depth the shortest depth reported in (Kerenidis et al., 2022), but with QR-style arguments rather than sophisticated data loaders. Moreover, and maybe more importantly for DPP sampling, while sampling from a DPP with *non-diagonalized* kernel remains limited by the initial (classical) diagonalization step, we argue that if one chooses the right avatar in the available distributed QR algorithms, we can even give a hybrid pipeline of a classical parallel and a quantum algorithm to sample some projection DPPs with *non-diagonalized* kernel, with a linear speedup compared to the vanilla classical DPP sampler of Hough et al. (2006). The covered DPPs include practically relevant cases such as the uniform spanning trees and column subsets of Examples 2.3 and 2.4. In Sections 5.3 and 5.4, we give two standard arguments, respectively due to (Hough et al., 2006) and (Lyons, 2003), to reduce the treatment of (non-projection) Hermitian DPPs to projection DPPs. This concludes our treatment of DPPs. In Section 5.5, we go back to connecting point processes to the work of Jiang et al. (2018), showing how one can use their circuit corresponding to the BdG Hamiltonian to sample a Pfaffian PP.

## 5.1 Quantum circuits

We refer to (Nielsen et al., 2010, Part II) for a description of quantum circuits and the basic building blocks. In short, in the quantum circuit model of quantum computation, one describes a computation by the initialization of a set of $N$ qubits, a sequence of unitary operators among a small set of physically-implementable operators called *gates*, and a physically-implementable observable called *measurement*. Not all gates can be implemented on every quantum computer, but software development kits like Qiskit
(Qiskit contributors, 2023) allow the user to define a quantum circuit using a large enough set of gates. The latter usually include any tensor product of identity matrices and Pauli matrices (14) (the so-called X, Y, and Z gates, respectively corresponding to $\sigma_x$, $\sigma_y$, and $\sigma_z$ in (14)) and a few two-qubit gates such as the ubiquitous CNOT. Qiskit then “transpiles” the resulting circuit into a sequence of actually-implementable gates for a given machine. For instance, with the current technology, not all qubits can be jointly operated on a given machine, e.g., one may not be able to apply a gate that acts on two qubits that are physically too far from each other in the actual quantum machine, or only be able to do so with a significant chance of error. Assuming the transpiling process preserves the dimensions of the circuit, one measures the complexity of a quantum circuit by quantities such as its total number of gates and its depth, i.e., the largest number of elementary gates applied to any single qubit.

Finally, in order to judge the possibility of sampling DPPs today, and be able to estimate what we can do in the future as hardware and software improve, it is important to keep in mind the main sources of error and their order of magnitude in current quantum hardware; see Appendix C.

## 5.2 The case of a projection DPP

Consider again $(a_j)$ to be the Jordan-Wigner operators, which satisfy (CAR), and $(|n\rangle)$ to be the corresponding Fock basis. Let $r \in \mathbb{N}^*$, $\mathbf{Q} \in \mathbb{C}^{r \times N}$ have orthonormal rows, and $b_j^* = \sum_{\ell=1}^N q_{j\ell} a_\ell^*$ for $1 \le j \le r$. Let further

$$|\psi\rangle = b_1^* \dots b_r^* |0\rangle. \qquad (33)$$

From Section 4.3, we know that simultaneously measuring $N_i = a_i^* a_i$, $i = 1, \dots, N$ on $\rho = |\psi\rangle\langle\psi|$ samples the projection DPP with kernel $\mathbf{K} = \mathbf{Q}^* \mathbf{Q}$. Measuring $N_i$ is easy, because the Fock basis of the JW operators is the computational basis, see Proposition 3.4, and measurement in the computational basis is a basic operation on any quantum computer. For the same reason, implementing $a_1^* \dots a_r^* |0\rangle$ simply amounts to initializing the first $r$ qubits to $|1\rangle$, and the rest to $|0\rangle$. It remains to be able to prepare the state (33), where the $b_i$'s intervene, not the $a_i$'s. This is done by a procedure akin to the QR algorithm by successive Givens rotations; a standard reference on QR algorithms and matrix computations in general is (Golub and Van Loan, 2012).

### 5.2.1 QR by Givens rotations

**Proposition 5.1** (Wecker et al., 2015; Kivlichan et al., 2018; Jiang et al., 2018). There is a unitary operator $\mathcal{U}(\mathbf{Q})$ on $\mathbb{H}$ such that

$$b_1^* \dots b_r^* |0\rangle = \mathcal{U}(\mathbf{Q}) a_1^* \dots a_r^* |0\rangle,$$

and $\mathcal{U}(\mathbf{Q})$ is a product of unitaries corresponding to elementary quantum gates.

*Proof.* The Givens rotation with parameters $\theta, \phi$ and indices $\ell^1, \ell^2 \in [N]$ is the unitary matrix

$$\mathbf{G} = \mathbf{G}_{\ell^1, \ell^2}(\theta, \phi) = \mathbf{P} \begin{pmatrix} \gamma & \mathbf{0} \\ \mathbf{0} & \mathbb{I} \end{pmatrix} \mathbf{P}^{-1} \quad \text{with} \quad \gamma = \begin{pmatrix} \cos \theta & e^{-i\phi} \sin \theta \\ -\sin \theta e^{i\phi} & \cos \theta \end{pmatrix}, \qquad (34)$$
where **P** is the matrix of the permutation (1 $\ell^1$)(2 $\ell^2$), allowing to select the vectors on which the rotation is applied.<fnref n="4" /> Choosing $\theta, \phi$ relevantly, one can make sure that $QG^*$ has a zero in position ($\ell^1, \ell^2$), while all columns other than ($\ell_1, \ell_2$) are left unchanged. Iteratively choosing ($\ell^1, \ell^2$), we can sequentially introduce zeros in **Q** by multiplying it by $n_G$ Givens rotations $G_1, \dots, G_{n_G}$, never changing an entry that we previously set to zero, until

$$QG_1^* \dots G_{n_G}^* = (\Lambda|0), \qquad (35)$$

where $\Lambda$ is an $r \times r$ diagonal unitary matrix, and the right block is the $r \times (N-r)$ zero matrix.

Now, to each Givens rotation $G = G_{\ell_1, \ell_2}(\theta, \phi)$, we associate a unitary operator $\mathcal{G} = G_{\ell_1, \ell_2}(\theta, \phi)$ on $\mathbb{H}$. We call $\mathcal{G}$ a Givens operator, as it realizes the Givens rotation $G$ by a conjugation of creation operators, i.e.

$$\mathcal{G}a_{\ell_1}^* \mathcal{G}^* = \cos \theta a_{\ell_1}^* + e^{-i\phi} \sin \theta a_{\ell_2}^*,$$

$$\mathcal{G}a_{\ell_2}^* \mathcal{G}^* = -e^{i\phi} \sin \theta a_{\ell_1}^* + \cos \theta a_{\ell_2}^*,$$

$$\mathcal{G}a_i^* \mathcal{G}^* = a_i^*, \quad \forall i \notin \{\ell^1, \ell^2\}.$$

We refer to Appendix A.1 for an explicit construction. For future convenience, we note that this action can be compactly summarized if one stacks up the creation operators in a vector $\mathbf{a}^* = (a_1^* \cdots a_N^*)^\top$, and writes

$$\mathcal{G}\mathbf{a}^*\mathcal{G}^* = \mathbf{G} \cdot \mathbf{a}^*, \qquad (36)$$

with $\cdot$ defined by matrix-vector multiplication.

Finally, consider the unitary

$$\mathcal{U}(\mathbf{Q}) \triangleq \mathcal{G}_{n_G} \dots \mathcal{G}_1, \qquad (37)$$

where $\mathcal{G}_i$ is the Givens operator corresponding to $\mathbf{G}_i$ in (35). Then, by construction and up to a phase factor, for all $1 \le k \le r$,

$$\mathcal{U}(\mathbf{Q}) a_k^* \mathcal{U}(\mathbf{Q})^* = b_k^*, \quad \text{and} \quad \mathcal{U}(\mathbf{Q}) |\emptyset\rangle = |\emptyset\rangle.$$

In particular,

$$|\psi\rangle = b_1^* \dots b_r^* |\emptyset\rangle = \mathcal{U}(\mathbf{Q}) a_1^* \mathcal{U}(\mathbf{Q})^* \dots \mathcal{U}(\mathbf{Q}) a_r^* \mathcal{U}(\mathbf{Q})^* |\emptyset\rangle \\ = \mathcal{U}(\mathbf{Q}) a_1^* \dots a_r^* |\emptyset\rangle,$$

again up to a phase factor, which is irrelevant in the resulting state $|\psi\rangle\langle\psi|$.

<figure><img src="image_3.png" /></figure>

The operator $\mathcal{U}(\mathbf{Q})$ in Proposition 5.1 is indeed a product of elementary two-qubit gates, because any Givens operator $\mathcal{G}_{\ell_1, \ell_2}(\theta, \phi)$ can be implemented as such, as first put forward by (Wecker et al., 2015). We discuss gate details in Appendix B.

To see how many gates are required for a given DPP kernel, we need to discuss a final degree of freedom: the Givens chain of rotations (34) is not unique. This is where constraints on which qubits can be jointly acted upon enter the picture.

<footnote n="4">⁴The parametrization of Givens rotations slightly differs from (Jiang et al., 2018).</footnote>
### 5.2.2 Parallel QR algorithms and qubit operation constraints

The product of Givens rotations in (34) can give rise to a quantum circuit of short depth if there are subproducts of rotations that can be performed on disjoint pairs of qubits. Independently of quantum computation, this is exactly the same kind of constraint that numerical algebraists have been studying in a long string of works on parallel QR factorization;<fnref n="5" /> see (Demmel et al., 2012) and references therein.

As a first example, after a preprocessing phase, the preparation of (33) in (Jiang et al., 2018, Section III) implicitly implements a parallel QR algorithm known as Sameh-Kuck (Sameh et al., 1978). More precisely, in a preprocessing phase, Jiang et al., 2018 zero out $r(r-1)/2$ entries in the upper-right corner of $Q$ by pre-multiplying $Q$ by a product of (complex) Givens rotations. Let us write this product of unitary matrices by $V$. This preprocessing requires $r(r-1)/2$ Givens rotations which are implemented thanks to a classical computer. To fix ideas, when $r = 3$ and $N = 6$, this results in a matrix of the following form

$$VQ = \begin{pmatrix} * & * & * & * & 0 & 0 \\ * & * & * & * & * & 0 \\ * & * & * & * & * & * \end{pmatrix}. \qquad (38)$$

Note that replacing $Q$ by $VQ$ does not change the kernel of the resulting DPP since $V^*V = I$. Now, Jiang et al., 2018 apply QR to $VQ$ as in the proof of Proposition 5.1, applying Givens rotations on disjoint pairs of neighbouring columns, as they become available. Rather than giving a formal description of the algorithm, available in (Sameh et al., 1978), we follow Jiang et al., 2018 and depict its application on Example (38). We use the convenient notation of Jiang et al., 2018, i.e. bold characters to show the most recently actively updated entries of the matrix, and underlined characters to show entries that automatically result from the rows of $Q$ being orthonormal. The successive parallel rounds of the algorithm are then

$$ \begin{aligned} VQ &= \begin{pmatrix} * & * & * & * & 0 & 0 \\ * & * & * & * & 0 & 0 \\ * & * & * & * & * & * \end{pmatrix} \rightarrow \begin{pmatrix} * & * & * & \mathbf{0} & 0 & 0 \\ * & * & * & * & * & 0 \\ * & * & * & * & * & * \end{pmatrix} \\ & \rightarrow \begin{pmatrix} * & * & \mathbf{0} & 0 & 0 & 0 \\ * & * & * & * & \mathbf{0} & 0 \\ * & * & * & * & * & * \end{pmatrix} \rightarrow \begin{pmatrix} \lambda_1 & \mathbf{0} & 0 & 0 & 0 & 0 \\ \mathbf{0} & * & * & * & \mathbf{0} & 0 \\ * & * & * & * & * & * \end{pmatrix} \\ & \rightarrow \begin{pmatrix} \lambda_1 & 0 & 0 & 0 & 0 & 0 \\ 0 & \lambda_2 & \mathbf{0} & 0 & 0 & 0 \\ 0 & \mathbf{0} & * & * & \mathbf{0} & 0 \end{pmatrix} \rightarrow \begin{pmatrix} \lambda_1 & 0 & 0 & 0 & 0 & 0 \\ 0 & \lambda_2 & 0 & 0 & 0 & 0 \\ 0 & 0 & \lambda_3 & \mathbf{0} & 0 & 0 \end{pmatrix}. \end{aligned} \qquad (39) $$

The factors $\lambda_i$ are of unit modulus, as in the construction behind Proposition 5.1. The upshot is that only a lower triangular matrix remains whose only non-zero entries are on the diagonal, while the other entries of the lower triangle automatically vanish due to row orthonormality.

The resulting circuit, an example of which is given in Figure 1, requires $O(Nr)$ gates, and has depth $O(N)$. The $O(r^2)$ preprocessing, which is optional and could also be performed in the quantum circuit, allows getting rid of a handful of (necessarily

<footnote n="5">⁵Note that, unlike the common QR decomposition, the matrix $R$ obtained here has more structure: it contains only one non-zero element per row. To some extent, the decomposition we discuss is close to a singular value decomposition; see the end of Section 5.2.3 for a short discussion.</footnote>
faulty) quantum gates at a small classical cost. Finally, constraining Givens operators to act on neighbouring qubits, or equivalently, Givens rotations to act on neighbouring columns, is practically relevant if the actual quantum computer at disposal favours two-qubit gates acting on neighbouring qubits.

As a second extreme example, assume that we have no constraint on which pairs of qubits can be jointly operated. We can then perform parallel QR on $Q$ using only $O(r \log N)$ parallel rounds, simply by acting on all available disjoint pairs in a single row until all but one entry are zeroed, keeping in mind that the leftmost $r$ entries of each row will be automatically updated because of orthonormality constraints. On an example with $N = 8$ and $r = 3$, this would yield

$$ \begin{align*} Q &\rightarrow \begin{pmatrix} * & \mathbf{0} & * & \mathbf{0} & * & \mathbf{0} & * & \mathbf{0} \\ * & * & * & * & * & * & * & * \\ * & * & * & * & * & * & * & * \end{pmatrix} \rightarrow \begin{pmatrix} * & 0 & \mathbf{0} & 0 & * & 0 & \mathbf{0} & 0 \\ * & * & * & * & * & * & * & * \\ * & * & * & * & * & * & * & * \end{pmatrix} \\ &\rightarrow \begin{pmatrix} \frac{\lambda_1}{0} & 0 & 0 & 0 & \mathbf{0} & 0 & 0 & 0 \\ 0 & * & * & * & * & * & * & * \\ 0 & * & * & * & * & * & * & * \end{pmatrix} \rightarrow \begin{pmatrix} \frac{\lambda_1}{0} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & * & * & \mathbf{0} & * & \mathbf{0} & * & \mathbf{0} \\ 0 & * & * & * & * & * & * & * \end{pmatrix} \\ &\rightarrow \begin{pmatrix} \lambda_1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & * & \mathbf{0} & 0 & * & 0 & \mathbf{0} & 0 \\ 0 & * & * & * & * & * & 0 & * \end{pmatrix} \rightarrow \begin{pmatrix} \lambda_1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{\lambda_2}{0} & 0 & 0 & \mathbf{0} & 0 & 0 & 0 \\ 0 & * & * & * & * & * & * & * \end{pmatrix} \\ &\rightarrow \begin{pmatrix} \lambda_1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \lambda_2 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & * & \mathbf{0} & * & \mathbf{0} & * & \mathbf{0} \end{pmatrix} \rightarrow \begin{pmatrix} \lambda_1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \lambda_2 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & * & 0 & * & 0 & * & 0 \end{pmatrix} \\ &\rightarrow \begin{pmatrix} \lambda_1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \lambda_2 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{\lambda_3}{0} & 0 & \mathbf{0} & 0 & 0 & 0 \end{pmatrix}. \end{align*} $$

The resulting quantum circuit has $O(Nr)$ gates as the one of Jiang et al., 2018, but a depth of only $O(r \log_2 N)$. This depth is similar to the shortest depth obtained by Kerenidis et al., 2022 with a similar tree-like structure called “parallel Clifford loaders”.

Just like in parallel QR, we argue that the choice of the chain (34) of rotations should depend on the particular hardware constraints, like which qubits we allow to be jointly operated. Parallel QR algorithms allow quite arbitrary dependency structures, see e.g. (Demmel et al., 2012, Section 4). While initially motivated by communication- or storage constraints in parallel classical computing, these dependency structures are actually tailored to designing quantum circuits to prepare states like (33) depending on a given quantum architecture. Taking the analogy with parallel QR even one step further, we now show that the initial bottleneck of diagonalizing the kernel of a DPP can be combined with the design of the quantum circuit in a single run of a parallel QR algorithm.

### 5.2.3 A hybrid parallel/quantum algorithm for projection DPPs

The link of Givens-based quantum circuits with parallel QR algorithms suggests pipelines to sample from projection DPPs, even when the kernel is not yet in diagonalized form.

**Proposition 5.2.** Assume that we have access only to $\mathbf{A} \in \mathbb{R}^{d \times N}$, and that we want to sample from DPP($\mathbf{K}$), where $\mathbf{K} = \mathbf{V}_{:, [r]} \mathbf{V}_{:, [r]}^*$ and $\mathbf{V}$ is defined by the singular value decomposition $\mathbf{A} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^*$. Then, given $P$ classical processors, for a run time
<figure><img src="image_4.png" /><figcaption>Figure 1: A circuit sampling a DPP with projection kernel of rank $r = 3$, with $N = 5$ items. On the left-hand side, Pauli X gates are used to create fermionic modes on the three first qubits. Note also the parallel QR Givens rotations on neibouring qubits indicated by parametrized $XX + YY$ gates. See Appendix B for gate details. On the right-hand side, measurements of occupation numbers are denoted by black squares.</figcaption></figure>

in $O(Nd^2/P)$ up to logarithmic factors in $P$, we can design a quantum circuit with depth $O(N \log P/P)$ and $O(Nd)$ gates, such that simultaneous measurements in the computational basis at the output of the circuit yield a sample of DPP($\mathbf{K}$).

Three comments are in order. First, specifying the kernel $\mathbf{K}$ as in Proposition 5.2, where $\mathbf{K}$ is the projection kernel onto the first $r$ principal components of a full-rank rectangular matrix $\mathbf{A}$, is common in practice. In particular, it covers the column subset selection of Example 2.4 and the uniform spanning forests of Example 2.3. Second, for simplicity we have omitted the additional $O(r^2)$ preprocessing introduced by Jiang et al., 2018. Third, at least if we neglect physical sources of noise in the quantum circuit (see Appendix C), the bottleneck in the hybrid sampler remains the initial classical cost $O(Nd^2/P)$. Compared with the vanilla classical approach, we have gained a linear speedup thanks to parallelization. Note that it is not clear how such a linear speedup can be gained in the vanilla classical algorithm, but that it can be gained in non-quantum randomized variants (Dereziński et al., 2019; Barthelmé et al., 2023). Dequantizing Proposition 5.2, to obtain a non-quantum algorithm with a similar cost, is actually an interesting question, which we leave to future work.

*Proof.* To prove Proposition 5.2, we first use a Givens-based parallel QR algorithm to compute $\mathbf{A}^* = \mathbf{QR}$; see (Demmel et al., 2012) and references therein for a recent entry. One can incorporate here communication constraints that will turn into qubit-communication constraints in the final algorithm. For simplicity, we assume $N \gg d$ and use the parallel tall-skinny QR of (Demmel et al., 2012, Section 2.1), with $P \le N/d$ processors. Without going into details, for $P$ a power of two, the idea is to partition $\mathbf{A}^*$ into $P$ equal $N/P \times d$ blocks, perform QR for each block using a user-chosen QR algorithm, and for $\log P$ stages, group the resulting $R$ matrices in pairs and perform QR. The run time is $O(Nd^2/P)$, up to logarithmic factors, which is a linear speedup compared to classical, non-parallel QR. Now, perform the singular value decomposition
of the $d \times d$ matrix $\mathbf{R} = \tilde{\mathbf{U}}\tilde{\mathbf{\Sigma}}\tilde{\mathbf{V}}^*$; this costs $\mathcal{O}(d^3)$. Then

$$ \mathbf{A}^*\mathbf{A} = \mathbf{QRR}^*\mathbf{Q}^* = (\mathbf{Q}\tilde{\mathbf{U}})\tilde{\mathbf{\Sigma}}^2(\mathbf{Q}\tilde{\mathbf{U}})^*, $$

so that the first $r$ principal components of $\mathbf{A}$ are $\tilde{\mathbf{U}}_{:, [r]}$. Essentially, we have used the QR algorithm to compute the principal directions: this detour is valuable because we simultaneously (i) benefit from parallelization and (ii) obtain $\mathbf{Q}$ as a product of Givens rotations by construction. If we now decompose $\tilde{\mathbf{U}}_{:, [r]}$ as a product of Givens rotations, say using the Sameh-Kuck algorithm like in (Jiang et al., 2018), we paid a (classical) preprocessing cost of $\mathcal{O}(Nd^2/P + d^3)$, up to logarithmic factors in $P$, to obtain all the information needed to create a quantum circuit with Givens gates that samples DPP($\mathbf{K}$) starting from the knowledge of $\mathbf{A}$. The depth depends on the blockwise QR algorithm used in the initial parallel tall-skinny QR step. If we use Givens rotations à la Sameh-Kuck to enforce only nearest-neighbour qubit interactions, we obtain a depth of $\mathcal{O}(N \log P/P)$, and a number of gates of $\mathcal{O}(Nd)$. The limiting factor of the overall pipeline remains the QR factorization of $\mathbf{A}$, in $\mathcal{O}(Nd^2/P)$ flops, but we have gained a linear speedup. $\square$

## 5.3 Reducing general DPPs to mixture DPPs

In Section 5.2, we gave quantum circuits to sample projection DPPs. If the kernel $\mathbf{K}$ is not that of a projection, a standard argument by Hough et al. (2006) allows writing the corresponding DPP as a statistical mixture of projection DPPs, thus yielding a quantum sampler for DPP($\mathbf{K}$) with a further classical preprocessing step, the latter implementing the choice of the mixture component.

More precisely, assume $\mathbf{K}$ is available in diagonalized form

$$ \mathbf{K} = \mathbf{U} \text{Diag}(\nu_k) \mathbf{U}^*, \quad (40) $$

with $\nu_1, \dots, \nu_N \in [0, 1]$. Then sampling DPP($\mathbf{K}$) can be done in two steps. First, sample $N$ independent Bernoulli random variables with respective parameters $\nu_1, \dots, \nu_N$ on a classical computer. Let $\mathcal{C} \subseteq [N]$ be the set of indices of successful Bernoulli trials. Second, measure all the occupation numbers $N_i$ for $1 \le i \le N$ in the quantum circuit of Section 5.2; in other words, sample from the projection DPP of correlation kernel $\mathbf{Q}^*\mathbf{Q}$ with $\mathbf{Q} = (\mathbf{U}_:\mathcal{C})^*$, where $\mathbf{U}_:\mathcal{C}$ is the set of columns indexed by $\mathcal{C}$.

To conclude this section, we note how one can arrive at this mixture construction by inspecting the quantum state “above” DPP($\mathbf{K}$). Consider the factorization of the state (18) obtained thanks to the diagonalization of $\mathbf{H}$, i.e.,

$$ \rho = \prod_{k=1}^{N} \frac{\exp(-\beta(\lambda_k - \mu)b_k^*b_k)}{1 + \exp(-\beta(\lambda_k - \mu))}, \quad (41) $$

where we used the formula (23) for the normalization constant. Note that the eigenvalues of the correlation kernel (40) are obtained by taking the sigmoid of the Hamiltonian eigenvalues, i.e., $\nu_k = \sigma(-\beta(\lambda_k - \mu))$, in the light of Proposition 4.1. By a simple inspection, the factor $k \in [N]$ in the product (41), is given by

$$ \frac{\exp(-\beta(\lambda_k - \mu)b_k^*b_k)}{1 + \exp(-\beta(\lambda_k - \mu))} = (1 - \nu_k) \times \pi_{\text{Ker } b_k} + \nu_k \times \pi_{\text{Ker } b_k^*}, $$

where $\pi_{\text{Ker } b_k}$ (resp. $\pi_{\text{Ker } b_k^*}$) is the orthogonal projector onto the null space of $b_k$ (resp. $b_k^*$). Inspecting Born's rule reveals that $\text{tr} \rho N_{i_1} \dots N_{i_\ell}$ is equal to the correlation
function of a statistical mixture of projection DPPs, each coming with a weight of the form

$$ \prod_{i \in I} \nu_i \prod_{j \notin I} (1 - \nu_j) $$

for some $I \subset [N]$. These weights correspond to the independent Bernoulli trials introduced by Hough et al., 2006 and discussed immediately after (40).

## 5.4 Dilating a general DPP to a projection DPP

In principle, there exists another strategy to sample from a general DPP by only using a quantum circuit such as described in Section 5.2 – thus bypassing the first step of the sampling algorithm consisting of the classical sampling of Bernoulli random variables. Any DPP on the ground set $\{1, ..., N\}$ can be dilated to a projection DPP on $\{1, ..., N, N + 1, ..., 2N\}$ by the transformation

$$ K \mapsto K_{\text{dil}} = \begin{pmatrix} K \\ (K(I - K))^{1/2} \end{pmatrix} \begin{pmatrix} (K(I - K))^{1/2} \\ I - K \end{pmatrix}, $$

see (Lyons, 2003, Section 8). The points of a sample of DPP($K_{\text{dil}}$) which belong to $\{1, ..., N\}$ yield a sample of DPP($K$).

## 5.5 Sampling the PfPP corresponding to the BdG Hamiltonian

We now turn to PfPPs. As above, we can write the mixed state of the Bogoliubov-de Gennes Hamiltonian at inverse temperature $\beta > 0$ as

$$ \frac{\exp(-\beta H_{\text{BdG}})}{Z} = \prod_{k=1}^{N} \frac{\exp(-\beta \epsilon_k b_k^* b_k)}{1 + \exp(-\beta \epsilon_k)}, \quad (42) $$

where $\epsilon_k$ is the energy of the eigenmode $k$. Hence, in complete analogy with the case of DPPs in Section 5.3, the product in (42) entails a statistical mixture formula. Therefore, sampling the corresponding Pfaffian point process amounts to proceed as follows:

1. sample $N$ independent Bernoullis with success probability $\sigma(-\beta\epsilon_k)$ on a classical computer, to obtain the set of successful indices $C \subseteq \{1, ..., N\}$,

2. jointly measure the occupation numbers $N_i = a_i^* a_i$ for all $1 \le i \le N$ in the pure state $\prod_{i \in C} b_i^* |\emptyset_b\rangle$ with a quantum circuit,

where $|\emptyset_b\rangle$ is the joint Fock vacuum of the $b_k$'s. Note that the preparation of $\prod_{i \in C} b_i^* |\emptyset_b\rangle$ is discussed in Section 5.5.1.

The output of this sampling algorithm is the set of indices for which the measure has given an occupation number equal to 1 rather than equal to 0. Again, this can be done easily thanks to the representation of the occupation numbers $N_i$'s in Proposition 3.4.

**Remark 5.3 (Projective S).** In the light of Proposition 4.8, the Pfaffian point process associated with the pure state $\prod_{i \in C} b_i^* |\emptyset\rangle$ is determined by the orthogonal projection matrix

$$ S = \overline{\mathbf{W}}^* \begin{pmatrix} \mathbf{I}_{\bar{C}} & \mathbf{0} \\ \mathbf{0} & \mathbf{I}_C \end{pmatrix} \overline{\mathbf{W}}, $$
where $\mathbf{I}_C$ is the diagonal matrix with entry $(i, i)$ equal to 1 if $i \in C$ and zero otherwise, whereas $\mathbf{I}_{\bar{c}} = \mathbf{I} - \mathbf{I}_C$. Notice that $\mathbf{C}\bar{\mathbf{S}}\mathbf{C} = \mathbf{I} - \mathbf{S}$ holds.

### 5.5.1 QR decomposition with double Givens rotations

We now describe how to adapt the discussion of Section 5.2 to the case where the particle number is not conserved. To simplify notation, we suppose that $\mathcal{C} = \{1, \dots, r\}$, and we prepare the state

$$b_1^* \dots b_r^* |\emptyset_b\rangle, \quad (43)$$

where the creation-annihilation operators $b_k$'s are obtained from the $a_k$'s thanks to a Bogoliubov transformation$^6$ (30) with a unitary orthogonal complex matrix $\mathbf{W}$ given in Lemma 4.6. Due to particle number non-conservation, the Fock vacuum of the $b_i$'s, denoted by $|\emptyset_b\rangle$ and given below, is not annihilated by the $a_i$'s. Also, recall that the Bogoliubov transformation is given by a matrix of the form

$$\mathbf{W} = \begin{pmatrix} \bar{\mathbf{W}}_1 & \bar{\mathbf{W}}_2 \\ \mathbf{W}_2 & \mathbf{W}_1 \end{pmatrix}, \quad (44)$$

where the blocks are such that this transformation is unitary; see Definition 4.4. Explicitly, as explained in (Jiang et al., 2018, Section IV), the full transformation is determined by the lower blocks of $\mathbf{W}$, which encode the expression of $\mathbf{b}$ in terms of $\mathbf{a}^*$ and $\mathbf{a}$. Now, we give the following result that we formalize and adapt from (Jiang et al., 2018, Equation (31)).

**Lemma 5.4.** There exists a $2N \times 2N$ orthogonal complex matrix $\mathbf{O}$ and an $N \times N$ unitary matrix $\mathbf{V}$ such that

$$\mathbf{V}(\mathbf{W}_2|\mathbf{W}_1)\mathbf{O}^* = (\mathbf{0}|\mathbf{D}),$$

where $\mathbf{D} = \text{Diag}(z_k)$ is a diagonal matrix with $|z_k| = 1$ for $1 \le k \le N$. Furthermore, the matrix $\mathbf{O}$ is here the following product of Givens rotations and particle-hole transformations

$$\mathbf{O} = \mathbf{B}\mathbf{M}_{N-1}\mathbf{B} \dots \mathbf{B}\mathbf{M}_1\mathbf{B}, \quad (45)$$

the terms of which we now explain. The matrix $\mathbf{M}_i$ is a product of $N-i$ double Givens rotation matrices. The double Givens rotation $\mathbf{\Gamma}$ associated with $(\theta, \phi)$ and vector indices $1 \le \ell^1 < \ell^2 \le N$ is defined as

$$\mathbf{\Gamma} = \begin{pmatrix} \mathbf{G}_{\ell^1, \ell^2}(\theta, \phi) & \mathbf{0} \\ \mathbf{0} & \bar{\mathbf{G}}_{\ell^1, \ell^2}(\theta, \phi) \end{pmatrix} \quad \text{with } \mathbf{G} \text{ defined in (34).}$$

The matrix $\mathbf{B}$ in (45) is a permutation matrix exchanging the last vector of the first block with the last vector of the second block. Formally, it is given by the so-called particle-hole matrix

$$\mathbf{B} = \begin{pmatrix} \mathbf{I} - \mathbf{e}_N \mathbf{e}_N^\top & \mathbf{e}_N \mathbf{e}_N^\top \\ \mathbf{e}_N \mathbf{e}_N^\top & \mathbf{I} - \mathbf{e}_N \mathbf{e}_N^\top \end{pmatrix}.$$

Finally, note that because we use the Jordan-Wigner transform, $\mathbf{B}$ can be implemented by a Pauli X gate.

<footnote>$^6$Recall that the $a_k$'s are the Jordan-Wigner operators representing the $c_k$'s in (30).</footnote>
Let us give a few details that sketch the proof of Lemma 5.4. On the one hand, the matrix $V$ is a product of single Givens rotations applied to the rows of $(\mathbf{W}_2 \mathbf{W}_1)$, to yield an upper triangle of zeros in the left block. For instance, when $N = 3$, we obtain

$$V (\mathbf{W}_2 | \mathbf{W}_1) = \begin{pmatrix} 0 & 0 & * & * & * & * \\ 0 & * & * & * & * & * \\ * & * & * & * & * & * \end{pmatrix}. \qquad (46)$$

On the other hand, the matrix $O$ is a product of double Givens rotations and particle-hole transforms, and further turns the left block into a block of zeros. More precisely, the next step consists in filling with zeros the remainder of the left block of (46), one anti-diagonal at a time, starting with the one in blue. For each anti-diagonal, a right-multiplication with the particle-hole transformation $B$ is used to zero the top right element, while right-multiplication with double Givens matrices fill with remainder of the anti-diagonal with zeros.

Before illustrating the procedure on our example, we pause to explain how to leverage the structure of the unitary orthogonal complex matrix $W$ to guess the structure of the entries in the right block. Denote by $\mathbf{r}_{1,i}$ and $\mathbf{r}_{2,i}$, the $i$th row vector of the first and second block, respectively. Thanks to Lemma 4.5, we have the following inner product identities between rows

$$\mathbf{r}_{1,i} \mathbf{r}_{2,j}^{\top} = - \mathbf{r}_{2,i} \mathbf{r}_{1,j}^{\top}, \qquad (47)$$

$$\mathbf{r}_{1,i} \mathbf{r}_{1,j}^{*} = \delta_{ij} - \mathbf{r}_{2,i} \mathbf{r}_{2,j}^{*}. \qquad (48)$$

Property (47) explains why we necessarily already have a zero in the right-hand block of (46).

Now, we show the right multiplication by $O$ on our example (46), starting with a particle-hole transformation on (46), namely

$$\begin{align*}
& \mathbf{V} (\mathbf{W}_2 | \mathbf{W}_1) \\
& \xrightarrow{\times \mathbf{B}} \begin{pmatrix} 0 & 0 & \mathbf{0} & * & * & * \\ 0 & * & * & * & * & * \\ * & * & * & * & * & * \end{pmatrix} \xrightarrow{\times \Gamma^*} \begin{pmatrix} 0 & 0 & 0 & * & * & * \\ 0 & \mathbf{0} & * & * & * & * \\ * & * & * & * & * & * \end{pmatrix} \quad \text{by (47)} \\
& \xrightarrow{\times \Gamma^*} \begin{pmatrix} 0 & 0 & 0 & \frac{\alpha}{\beta_{(2)}} & \frac{\alpha}{\beta_{(1)}} & 0 \\ 0 & 0 & * & * & * & * \\ \mathbf{0} & * & * & * & * & * \end{pmatrix} \quad \text{by (47) at (1), (48) at (2), (47) at (3)} \\
& \xrightarrow{\times \mathbf{B}} \begin{pmatrix} 0 & 0 & 0 & \alpha & 0 & 0 \\ 0 & 0 & \mathbf{0} & \frac{\alpha}{\beta_{(2)}} & * & * \\ 0 & * & * & * & * & * \end{pmatrix} \quad \text{by (48)} \\
& \xrightarrow{\times \Gamma^*} \begin{pmatrix} 0 & 0 & 0 & \alpha & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{\beta}{\beta_{(3)}} & 0 \\ 0 & \mathbf{0} & \frac{\gamma}{\beta_{(4)}} & 0 & \frac{\alpha}{\beta_{(2)}} & \frac{\alpha}{\beta_{(1)}} \end{pmatrix} \quad \text{by (47) at (1), by (48) at (2), (3), (4)} \\
& \xrightarrow{\times \mathbf{B}} \begin{pmatrix} 0 & 0 & 0 & \alpha & 0 & 0 \\ 0 & 0 & 0 & 0 & \beta & 0 \\ 0 & 0 & 0 & 0 & 0 & \gamma \end{pmatrix} = \mathbf{V} (\mathbf{W}_2 | \mathbf{W}_1) \mathbf{O}^*,
\end{align*}$$

with unit modulus complex numbers $\alpha, \beta, \gamma$. To help the reader understand which of properties (47) and (48) applied, an index was added to matrix elements which are automatically set to zero.
**Proposition 5.5** (Formalization of a result of Jiang et al., 2018). The number of double Givens gates as defined in Lemma 5.4 to achieve the factorization (45) is $N(N-1)/2$, whereas the number of particle-hole transformations is $N$.

Lemma 5.4 is now leveraged to yield the factorization of $\mathbf{W}$ in (44).

**Lemma 5.6.** The matrix $\mathbf{O}$ in (45) is a unitary orthogonal complex matrix, and we have the factorization

$$ \begin{pmatrix} \mathbf{D}\bar{\mathbf{V}} & \mathbf{0} \\ \mathbf{0} & \bar{\mathbf{D}}\mathbf{V} \end{pmatrix} \mathbf{W}\mathbf{O}^* = \begin{pmatrix} \mathbf{I} & \mathbf{0} \\ \mathbf{0} & \mathbf{I} \end{pmatrix}. $$

*Proof.* The discussion above yields the expression of the lower block. Since $\mathbf{O}$ and $\mathbf{W}$ are orthogonal complex matrices and unitary as well, the upper block is determined by the lower block in the light of Definition 4.4. $\square$

Furthermore, the $\mathbf{W}$ matrix of Lemma 5.6 yields a Bogoliubov transformation of the following factorized form

$$ \begin{pmatrix} \mathbf{b}^* \\ \mathbf{b} \end{pmatrix} = \begin{pmatrix} \mathbf{V}^\top \bar{\mathbf{D}} & \mathbf{0} \\ \mathbf{0} & \mathbf{V}^* \mathbf{D} \end{pmatrix} \begin{pmatrix} \tilde{\mathbf{b}}^* \\ \tilde{\mathbf{b}} \end{pmatrix} \quad \text{with} \quad \begin{pmatrix} \tilde{\mathbf{b}}^* \\ \tilde{\mathbf{b}} \end{pmatrix} \triangleq \mathbf{O} \begin{pmatrix} \mathbf{a}^* \\ \mathbf{a} \end{pmatrix}, $$

which gives rise to two remarks.

**Remark 5.7.** By observing that the creation operators $\mathbf{b}^*$ only depend on $\tilde{\mathbf{b}}^*$ and not $\tilde{\mathbf{b}}$, we can directly use the factorization given in Section 5.2. To do so, for $k \in [r]$, we write

$$ b_k^* = \sum_{j=1}^{N} \mathbf{Q}_{kj} \tilde{b}_j^*, $$

where $\mathbf{Q} = (\mathbf{V}^\top \bar{\mathbf{D}})_{[r]}$ is the matrix made of the r first rows of $\mathbf{V}^\top \bar{\mathbf{D}}$. The product of Givens transformations factorizing this matrix is $\mathbf{U}(\mathbf{Q})$ as given in Section 5.2.1, which is realized by the operator $\mathcal{U}(\mathbf{Q}) : \mathbb{H} \to \mathbb{H}$; see (37). Finally, the desired state (43) reads

$$ b_1^* \dots b_r^* |\emptyset_b\rangle = \mathcal{U}(\mathbf{Q}) \tilde{b}_1^* \dots \tilde{b}_r^* |\emptyset_{\tilde{b}}\rangle, \quad (49) $$

where the Fock vacuum is $|\emptyset_b\rangle = |\emptyset_{\tilde{b}}\rangle$ and with $\mathbf{Q} = (\mathbf{V}^\top \bar{\mathbf{D}})_{[r]}$.

**Remark 5.8.** We can define a unitary $\mathcal{W} = \mathcal{W}(\mathbf{O}) : \mathbb{H} \to \mathbb{H}$ representing the multiplication by $\mathbf{O}$ as follows

$$ \begin{pmatrix} \mathcal{W}\mathbf{a}^*\mathcal{W}^* \\ \mathcal{W}\mathbf{a}\mathcal{W}^* \end{pmatrix} = \mathbf{O} \begin{pmatrix} \mathbf{a}^* \\ \mathbf{a} \end{pmatrix} = \begin{pmatrix} \tilde{\mathbf{b}}^* \\ \tilde{\mathbf{b}} \end{pmatrix}, $$

where the action of $\mathcal{W}$ is entrywise on the vectors $\mathbf{a}^*$ and $\mathbf{a}$. By construction, we have

$$ \mathcal{W}(\mathbf{O}) = \mathcal{B}\mathcal{M}_1\mathcal{B}\mathcal{M}_2\mathcal{B}\dots\mathcal{B}\mathcal{M}_{N-1}\mathcal{B}, \quad (50) $$

where $\mathcal{B}$ is such that $\mathcal{B}\mathbf{a}_N^*\mathcal{B}^* = \mathbf{a}_N$ while leaving all other creation-annihilation operators unchanged, while $\mathcal{M}_i$ is a composition of $N - i$ Givens operators. Each Givens operator $\Gamma : \mathbb{H} \to \mathbb{H}$ simply represents multiplication by the Givens matrix $\Gamma$ as follows

$$ \begin{pmatrix} \Gamma \mathbf{a}^* \Gamma^* \\ \Gamma \mathbf{a} \Gamma^* \end{pmatrix} = \Gamma \cdot \begin{pmatrix} \mathbf{a}^* \\ \mathbf{a} \end{pmatrix}. $$
Again, the Givens operators in (50) appear in reverse order compared with Givens matrices in (45). Also, contrary to the particle number preserving case, the Fock vacuum of the $a_i$'s does not coincide with the Fock vacuum of the operators $\tilde{b}_1, \dots, \tilde{b}_N$, but we rather have $|\emptyset_{\tilde{b}}\rangle = W|\emptyset_a\rangle$.

### 5.5.2 Quantum gates

We now combine the two remarks above. Let $V$ and $O$ be given by Lemma 5.4. For short, denote $Q = (V^\top D)_{[r]:}$. The factorization of the fermionic Gaussian state (43) as a composition of double Givens gates and $X$ gates reads

$$b_1^* \dots b_r^* |\emptyset_b\rangle = U(Q)W(O) \cdot a_1^* \dots a_r^* |\emptyset_a\rangle.$$

This result is obtained from (49) by using $\tilde{b}_k^* = Wa_k^*W^*$ with $1 \le k \le N$.

# 6 Numerical experiments

In this section, we demonstrate the circuits discussed in Section 5 on a Qiskit simulator (Qiskit contributors, 2023), and on a few 5-qubit IBMQ machines (IBM Quantum, 2021). Python code to reproduce the experiments in available on GitHub.<fnref n="7" />

## 6.1 Projection DPPs

We start with a projection DPP on $\{1, 2, 3, 4, 5\}$ with rank 3. The kernel is $K = Q_{:, [k]} Q_{:, [k]}^\top$, where $Q$ is obtained by the (Householder-)QR decomposition of an $N \times N$ matrix full of i.i.d. real Gaussians. We show the kernel in Figure 2a. Following the circuit construction of Jiang et al., 2018 with 2-qubit gates acting only on neighbouring qubits on a line, we obtain the circuit shown in Figure 1, where each gate labeled as “XX+YY” is a Givens gate, i.e., a sequence of CNOT and $Z$ gates as discussed in Appendix B.

Now, before it can be run on a particular machine, the circuit needs to be transpiled, i.e. written as an equivalent sequence of gates that correspond to what can be physically implemented on the machine. We show in Figure 3 the calibration data for three 5-qubit IBMQ machines: `lima`, `quito`, and `manila`. This calibration takes the form of a graph, where nodes represent qubits, and edges represent the possibility of applying a CNOT gate, the only two-qubit gate in our original circuit in Figure 1. While `manila` can implement CNOT gates between neighbouring qubits on a line, as implicitly assumed in the construction of Jiang et al., 2018, the two other machines have a T-structured graph that will force the transpiled circuit to look quite different from Figure 1. Indeed, we show the transpiled circuits in the first three panels of Figure 4. Note how `manila` uses CNOT gates between neighbouring qubits, as in the original circuit, resulting in overall fewer CNOT gates than on `quito` and `lima`. The latter two machines cannot afford CNOT gates jointly acting on qubits 2 and 3, for instance, and end up compensating for this by using more CNOT gates in total. Since CNOT gates are among the most error-prone manipulations, minimizing the number of CNOT gates is an important feature. Intuitively, had we known in advance that we would run the circuit on a machine with a particular graph, we should have designed

<footnote n="7">⁷https://github.com/For-a-few-DPPs-more/quantum-sampling-DPPs</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 3: Calibration data for three IBMQ machines: lima, quito, and manila. Dark blue is low, light violet is high, but colors are not comparable across subfigures, although the orders of magnitude are similar. Node colors indicate readout errors (all of the order of $10^{-2}$), edge colors indicate CNOT error, of the same $10^{-2}$ order of magnitude for lima and quito, but going down to $5 \times 10^{-3}$ for manila. All machines have 5 qubits, but their volumes are respectively 8, 16, and 32.</figcaption></figure>

the circuit in Figure 1 differently, by rather running parallel QR with Givens rotations only along edges of the machine’s graph.

Before we observe the results, we note that the three machines come with different characteristics. For instance, manila has overall the lowest readout errors, and lima the largest. It is common to summarize the characteristics of a machine in a single number $2^{d_{max}}$, called *quantum volume*, where $d_{max}$ is –loosely speaking– the depth of a square circuit that we can expect to run reliably. The volume reported by IBM is obtained numerically, using a procedure known as randomized benchmarking; see Appendix C. The machines manila, quito, and lima respectively have reported volumes 32, 16, and 8. Thus, even on manila, the transpiled circuit is much larger than the “guaranteed” $5 \times 5$ square circuit, and we should expect noise in our measurements, as we shall now see.

Figure 4d shows the empirical distributions corresponding to independently preparing the input and measuring the output of the transpiled circuits in Figure 4, 20,000 times each. In blue and orange, we show for reference the probability under DPP(K) of the corresponding subset, as well as the empirical frequencies coming from sampling the output of the circuit using a simulator, which amounts to independently drawing 20,000 samples of the DPP on a classical computer.
<figure><img src="image_6.png" /><figcaption>Figure 4: Transpiled circuits corresponding to the input in Figure 1 for three IBMQ machines: lima, quito, and manila. Figure 4d shows the corresponding empirical distribution.</figcaption></figure>
The empirical measure of the classical samples in orange is close to the underlying DPP, as testified by the total variation (TV) distance<fnref n="8" /> between the two distributions, which is below $10^{-2}$. Actually, if we repeatedly and independently draw sets of 20,000 samples, we obtain the empirical distribution of the TV distance shown in Figure 2b. In contrast, the TV distance between the empirical measure obtained from 20,000 runs on manila of the corresponding transpiled circuit is 0.27, while it is 0.39 and 0.40 for the T-structured machines lima and quito, respectively. Looking at Figure 2b again, a test based on the estimated TV distance would easily reject even the hypothesis that at least one of the quantum circuits samples from the correct distribution. Moreover, we confirm (not shown) that the difference in TV between manila and its competitors is significant, which confirms the intuition coming from its larger volume and smaller number of CNOT gates, due to a QR decomposition adapted to its qubit communication graph. Finally, although a test would reject that the quantum circuits sample from the correct DPP, the resulting distributions are still close to the target DPP, especially for manila, as confirmed by their respective TV. Interestingly, the quantum circuits actually yield point processes that are supported on (almost) all subsets of {1, ..., 5}, while the target DPP, being of projection, only charges subsets of cardinality 3. The noise does respect the structure of the DPP, somehow: two-item subsets that (wrongly) appear in the support of the empirical measures correspond to items that are marginally favoured by the DPP, as can be seen on the diagonal of the kernel in Figure 2a. Intuitively, the appearance of a subset of cardinality 2 is partly due to readout error on one of the qubits supposed to output 1. This intuition is confirmed by the calibration data: for lima, for instance, Figure 3a shows that the qubit labeled '3' has a large readout error; simultaneously, there is a deficit of appearance of the triplet with labels {0, 3, 4} in its empirical distribution in Figure 4d, while {0, 4} wrongly appears.

**Optimal QR for a T-structured communication graph.** The QR-inspired fermionic circuits implemented in Qiskit follow (Jiang et al., 2018), and thus use Givens rotations between neighbouring columns. While this suits the calibration data of manila, which has a linear qubit-communication graph, lima and quito rather have a communication graph shaped as a $T$. As a result, transpilation is less straightforward and yields a bigger circuit than for manila. In particular, the transpiled quito circuit in Figure 4b has 15 CNOT gates, while both the original circuit in Figure 1 and the transpiled manila circuit have only 12 CNOT gates.

As we discussed in Section 5, a shorter (and less error-prone) transpiled circuit would intuitively result from a QR decomposition that respects the qubit-communication graph. For concreteness, we give a QR decomposition that is better suited to the communication graph of quito, with the mapping to the columns of $Q$ given in Figure 5.

<footnote n="8">⁸Recall that the TV is the maximum absolute difference between the probabilities assigned to a subset, where the maximum is taken across subsets.</footnote>
Assuming no preprocessing, we find

$$
\begin{equation}
\begin{split}
\mathbf{Q} &\rightarrow \begin{pmatrix} * & * & \mathbf{0} & * & \mathbf{0} \\ * & * & * & * & * \\ * & * & * & * & * \end{pmatrix} \rightarrow \begin{pmatrix} * & * & 0 & \mathbf{0} & 0 \\ * & * & * & * & * \\ * & * & * & * & * \end{pmatrix} \\
&\rightarrow \begin{pmatrix} \frac{\lambda_1}{0} & \mathbf{0} & 0 & 0 & 0 \\ 0 & * & * & * & \mathbf{0} \\ 0 & * & * & * & 0 \end{pmatrix} \rightarrow \begin{pmatrix} \lambda_1 & 0 & 0 & 0 & 0 \\ 0 & * & * & \mathbf{0} & 0 \\ 0 & * & * & * & * \end{pmatrix} \\
&\rightarrow \begin{pmatrix} \lambda_1 & 0 & 0 & 0 & 0 \\ 0 & \mathbf{0} & \frac{\lambda_2}{0} & 0 & 0 \\ 0 & * & \frac{\lambda_3}{0} & * & \mathbf{0} \end{pmatrix} \rightarrow \begin{pmatrix} \lambda_1 & 0 & 0 & 0 & 0 \\ 0 & 0 & \lambda_2 & 0 & 0 \\ 0 & \frac{\lambda_3}{0} & 0 & \mathbf{0} & 0 \end{pmatrix}.
\end{split}
\tag{51}
\end{equation}
$$

While not necessary optimal in any sense, our guiding principle for the decomposition (51) is to fill the matrix with zeros such that, for each row, the final complex phase appears at a node which, if removed from the graph, leaves the resulting graph connected. Note that to fall back onto a “diagonal” matrix, although this last step is unnecessary for DPP sampling, a final permutation between the second and third columns can be realized by an extra Givens gate with $\theta = 0$ and $\phi = \pi/2$ (i.e., a so-called ISwap gate) between qubit 1 and qubit 2. Overall, the sequence of Givens rotations corresponding to (51) can be transpiled on *quito* or *lima* with only the expected 2-per-rotation CNOT gates, since all rotations are applied to neighbours in the graph. We leave the characterization and benchmarking of the optimal QR decomposition for a given qubit communication graph for future work.

<figure><img src="image_7.png" /><figcaption>Figure 5: Allowed CNOT gates between qubits in the T-structured communication graph of *lima* or *quito*, for a map sending column 1 $\mapsto$ qubit 0, column 2 $\mapsto$ qubit 1, etc.</figcaption></figure>

## 6.2 Pfaffian PPs

In this section, we illustrate the second step of the algorithm of Section 5.5 to sample PfPPs of the type described by Koshida (2021) and associated with a pure state of the form $\prod_{i \in C} b_i^* |0_b\rangle$, namely, for which the matrix $S$ is projective as explained in Remark 5.3.

To construct $S$, we consider the quadratic form $H_{BdG}$ given in (28) with $N = 5$
with the Hermitian and skew-symmetric part given respectively by

$$ M = \begin{pmatrix} 1 & 0.5 & 0.2 & 0.2 & 0.2 \\ 0.5 & 1 & 0.5 & 0.2 & 0.2 \\ 0.2 & 0.5 & 1 & 0.5 & 0.2 \\ 0.2 & 0.2 & 0.5 & 1 & 0.5 \\ 0.2 & 0.2 & 0.2 & 0.5 & 1 \end{pmatrix} \quad \text{and} \quad \Delta = \begin{pmatrix} 0 & 1 & 0 & 0 & 0 \\ -1 & 0 & 1 & 0 & 0 \\ 0 & -1 & 0 & 1 & 0 \\ 0 & 0 & -1 & 0 & 1 \\ 0 & 0 & 0 & -1 & 0 \end{pmatrix}. $$

Next, we select the subset $\mathcal{C}$ of the 3 smallest positive eigenvalues of $H_{BdG}$. Let $K$ be the kernel associated with the projection

$$ S = \overline{\mathbf{W}}^* \begin{pmatrix} \mathbf{1}_c & \mathbf{0} \\ \mathbf{0} & \mathbf{1}_c \end{pmatrix} \overline{\mathbf{W}}, $$

by virtue of Proposition 2.6. The Hamiltonian quadratic form and the associated Pfaffian kernel are displayed in Figure 6. In this case, PfPP($K$) is a Pfaffian L-ensemble (Borodin and Rains, 2005, Proposition 1.6), for which the probability mass function has the following simple expression:

$$ \mathbb{P}(Y = S) = \frac{\mathrm{Pf}(\mathbf{L}_S)}{\mathrm{Pf}(\mathbf{J} + \mathbf{L})}, \qquad (52) $$

with the Gram matrix $\mathbf{L}_S = [\mathbb{L}(i, j)]_{i, j \in S}$ and where $\mathbb{L} = (K - \mathcal{J})^{-1} - \mathcal{J}$ is an $N \times N$ skew-symmetric matrix made of 2 × 2 blocks. For this numerical experiments, we only

<figure><img src="image_8.png" /><figcaption>Figure 6: Visualization of the Hamiltonian $H_{BdG}$ and the components defining the Pfaffian kernel $K$.</figcaption></figure>

compare (classical) simulations of the quantum circuits to the true point process. As in Section 6.1, 20,000 samples are drawn independently from $Y \sim \mathrm{PfPP}(K)$ thanks to the Qiskit simulator. Next, the empirical frequencies of subsets of $[N]$ are computed and compared with the probabilities (52).

In Figure 7, we observe that the empirical probabilities match the expected case, with a TV distance of 0.009. Furthermore, it is manifest from Figure 6b that there is a weak repulsion between items 2 and item 4 – see the entry (1, 3) in the grayscale matrix – and that each of these two elements has a large marginal probability. Hence, they can be expected to be sampled together. This is confirmed in Figure 7 where the subset {2, 4} (corresponding to the label (1, 3)) corresponds to a large mass under the empirical measure. Also, note that all subsets naturally have the same parity. For completeness, the histogram of the total variation distance between the ground truth distribution and the estimated distribution, computed over 1000 runs, is displayed in Figure 8.
<figure><img src="image_9.png" /><figcaption>Figure 7: Histogram of the empirical frequencies of the subsets sampled by PfPP($K$) vs subset probabilities as given by (52).</figcaption></figure>

# 7 Discussion

Inspired by the pioneering work of Macchi (1972) on point processes in quantum physics, we have studied quantum algorithms for DPP sampling. We did so by reducing DPP sampling to a fermion sampling problem, and then leveraging and modifying recent algorithms for fermion sampling (Wecker et al., 2015; Jiang et al., 2018). While many of the steps are either common lore in one or the other field, or recently published material, we believe that there is value in a self-contained survey of how to reduce a finite DPP to a fermion sampling problem, all the way from the mathematical object to the implementation on the quantum machine. We hope that this paper can help spark further cross-disciplinary work. Moreover, writing down all the steps from a DPP as specified in machine learning to its quantum implementation has allowed us to make contributions on top of the survey, like the extension of the argument to a class of Pfaffian PPs, and the first steps in adapting the QR-decomposition behind fermion sampling algorithms to qubit communication constraints. This opens several research avenues.

First, as mentioned in the introduction, projective DPPs can also be sampled thanks to the Clifford loaders defined in (Kerenidis et al., 2022), introduced independently of the physics literature that we cover in this paper. Yet the structure of the arguments is related, and it would be enlightening to explicitly compare, on the one hand, the parallel implementation we give in Section 5.2.2 with a circuit depth logarithmic in $N$ and, on the other hand, data loaders of Kerenidis et al. (2022) with a similar depth.
<figure><img src="image_10.png" /><figcaption>Figure 8: Empirical distribution of the total variation distance between the true probability mass function of the PfPP given by Figure 6b and its simulation.</figcaption></figure>

Second and in the same spirit, an interesting extension of this work would be to develop an algorithm optimally matching any qubit communication graph to a QR decomposition scheme, where by *optimal* we mean minimizing e.g. the total variation distance between the output of the circuit and the original point process. This would generalize the case of the *T* graph described in Section 6.1, and lead to transpilers with fewer noisy gates. A potential strategy would be to follow the approach of Frerix and Bruna, 2019, who optimize a sparsity-inducing objective to approximate a unitary matrix as a product of Givens rotations.

Third, a natural improvement of the proposed classical preprocessing followed by circuit construction would be the inclusion of the kernel diagonalization *in* the quantum circuit, using for instance the recent developments about quantum SVD (Rebentrost, Steffens, Marvian, and Lloyd, 2018). The combination of this quantum preprocessing and a QR-based circuit would constitute a turn-key sampling pipeline.

A fourth and maybe more speculative research perspective would be to leverage our knowledge of all the correlation functions of point processes such as DPPs, PfPPs, and permanental PPs (Jahangiri, Arrazola, Quesada, and Killoran, 2020) to develop a statistical test of quantum decoherence in a given machine. In particular, our aim would be to design statistics of PPs which are sensitive to the different kinds of errors affecting a quantum computer.

Finally, from a mathematical perspective, we think it is worth exploring in more depth the structure of PfPPs. While potentially offering more modeling power in machine learning applications, they have received little attention, likely due to their high sampling cost on a classical computer. Since the sampling overhead on a quantum computer is minor, they are likely to become a valuable addition to the machine learner's toolbox. For starters, we are unaware of a formula for the probability mass function generalizing (52) to the case where $K-J$ is not invertible. Such a construction would generalize the *extended* L-ensemble construction in (Tremblay, Barthelmé, Usevich, and Amblard, 2023).

## Acknowledgements

We acknowledge support from ERC grant Blackjack (ERC-2019-STG-851866) and ANR AI chair Baccarat (ANR-20-CHIA-0002). Furthermore, we acknowledge the use of IBM Quantum services for this work. The views expressed are those of the authors, and do not reflect the official policy or position of IBM or the IBM Quantum team.
# A Mathematical details

## A.1 More about Givens operators

We briefly explain here how the unitary operator $\mathcal{G}$ representing a Givens rotation in (36) can be constructed. A unitary operator $\mathcal{V}: \mathbb{H} \to \mathbb{H}$ such that

$$ \begin{pmatrix} \mathcal{V}a_1^* \mathcal{V}^* \\ \mathcal{V}a_2^* \mathcal{V}^* \end{pmatrix} = \begin{pmatrix} \cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{pmatrix} \begin{pmatrix} a_1^* \\ a_2^* \end{pmatrix} $$

is given by $\mathcal{V} = \exp(\theta(a_2^*a_1 - a_1^*a_2))$, where we used the BCH formula (Hall, 2010). Thanks to this observation, we readily have that

$$ \mathcal{G} = \mathcal{T}\mathcal{V}\mathcal{T}^*, $$

with $\mathcal{T} = \exp(i\phi(a_1^*a_1 - a_2^*a_2)/2)$, realizes the Givens transformation

$$ \begin{pmatrix} \mathcal{G}a_1^*\mathcal{G}^* \\ \mathcal{G}a_2^*\mathcal{G}^* \end{pmatrix} = \begin{pmatrix} \cos \theta & e^{-i\phi} \sin \theta \\ -e^{i\phi} \sin \theta & \cos \theta \end{pmatrix} \begin{pmatrix} a_1^* \\ a_2^* \end{pmatrix}. $$

## A.2 Proof of Lemma 4.7

We begin by a remark about the diagonalization of $H_{\text{BdG}}$. Let $\mathbf{C} = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}$. As a consequence of Lemma 4.6, we can write

$$ H_{\text{BdG}} = \frac{1}{2} \begin{pmatrix} \mathbf{c}^* \\ \mathbf{c} \end{pmatrix}^* \begin{pmatrix} -\overline{\mathbf{M}} & -\overline{\Delta} \\ \Delta & \overline{\mathbf{M}} \end{pmatrix} \begin{pmatrix} \mathbf{c}^* \\ \mathbf{c} \end{pmatrix}. $$

The diagonalization of $H_{\text{BdG}}$ reads

$$ H_{\text{BdG}} = \frac{1}{2} \begin{pmatrix} \mathbf{b}^* \\ \mathbf{b} \end{pmatrix}^* \begin{pmatrix} -\text{Diag}(\epsilon_k) & \mathbf{0} \\ \mathbf{0} & \text{Diag}(\epsilon_k) \end{pmatrix} \begin{pmatrix} \mathbf{b}^* \\ \mathbf{b} \end{pmatrix}. $$

Note that the expectation value of the bilinears are

$$ \begin{align*} \begin{pmatrix} \langle c_i c_j^* \rangle_\rho & \langle c_i c_j \rangle_\rho \\ \langle c_i^* c_j^* \rangle_\rho & \langle c_i^* c_j \rangle_\rho \end{pmatrix}_{i,j} &= \mathbf{W}^\top \begin{pmatrix} \langle b_i b_j^* \rangle_{i,j} & \langle b_i b_j \rangle_{i,j} \\ \langle b_i^* b_j^* \rangle_{i,j} & \langle b_i^* b_j \rangle_{i,j} \end{pmatrix}_{i,j} \overline{\mathbf{W}} \\ &= \mathbf{W}^\top \begin{pmatrix} (1-d_i)\delta_{ij})_{i,j} & (0)_{i,j} \\ (0)_{i,j} & (d_i\delta_{ij})_{i,j} \end{pmatrix} \overline{\mathbf{W}} \\ &= \mathbf{W}^\top \begin{pmatrix} \mathbf{I} - \text{Diag}(d_k) & \mathbf{0} \\ \mathbf{0} & \text{Diag}(d_k) \end{pmatrix} \overline{\mathbf{W}} \end{align*} $$

with $d_k = \exp(-\beta\epsilon_k)/(1 + \exp(-\beta\epsilon_k))$ for all $1 \le k \le N$. The proof is completed if we recall that the sigmoid $\sigma$ satisfies $\sigma(x) = 1 - \sigma(-x)$.

## A.3 Proof of Proposition 4.8

Assume that $i_1, \dots, i_k$ are pairwisely distinct. The object of interest is

$$ \langle N_{i_1} \dots N_{i_k} \rangle_\rho = \langle c_{i_1}^* c_{i_1} \dots c_{i_k}^* c_{i_k} \rangle_\rho. $$
We now use a direct consequence of Wick's theorem for expectations in a Gaussian state, see (Bardenet et al., 2022, Theorem 3): let $m$ even and let $\beta_1, \dots, \beta_m$ be linear combinations of the $c_i$'s and $c_i^*$'s for $1 \le i \le N$, then

$$ \langle \beta_1 \dots \beta_m \rangle_{\rho} = \sum_{\sigma \text{ contraction}} \operatorname{sgn}(\sigma) \langle \beta_{\sigma(1)} \beta_{\sigma(2)} \rangle_{\rho} \dots \langle \beta_{\sigma(m-1)} \beta_{\sigma(m)} \rangle_{\rho}. \quad (53) $$

Take $m = 2k$ and use the following definition: $\beta_{2\ell-1} = c_{i\ell}^*$ and $\beta_{2\ell} = c_{i\ell}$ for $1 \le \ell \le k$. Thus, we now show that (53) is the Pfaffian of a skewsymmetric matrix made of the $2 \times 2$ blocks.

Let us construct this skewsymmetric matrix. In details, for the pair $(\ell, \ell')$ with $1 \le \ell < \ell' \le k$, we denote the $2 \times 2$ block by $\mathbb{K}_{\ell\ell'}$. In order to make any block matrix $(\mathbb{K}_{i_m i_n})_{1 \le m, n \le k}$ skewsymmetric, we need to have

$$ \mathbb{K}_{\ell\ell'}^{\top} = -\mathbb{K}_{\ell'\ell}. \quad (54) $$

The form of the $(\ell, \ell')$ block with $1 \le \ell < \ell' \le k$ is

$$ \mathbb{K}_{\ell\ell'} = \begin{pmatrix} \langle \beta_{2\ell-1} \beta_{2\ell'-1} \rangle_{\rho} & \langle \beta_{2\ell-1} \beta_{2\ell'} \rangle_{\rho} \\ \star & \langle \beta_{2\ell} \beta_{2\ell'} \rangle_{\rho} \end{pmatrix} = \begin{pmatrix} \langle c_{i\ell}^* c_{i\ell'}^* \rangle_{\rho} & \langle c_{i\ell}^* c_{i\ell'} \rangle_{\rho} \\ \star & \langle c_{i\ell} c_{i\ell'} \rangle_{\rho} \end{pmatrix}. $$

Note that $\langle \beta_{2\ell-1} \beta_{2\ell'} \rangle_{\rho}$ does not appear in (53) since $\ell < \ell'$. Now, the $\star$ is completed by $-\langle c_{i\ell}^* c_{i\ell'} \rangle_{\rho}$ to ensure the property (54) which guarantees skewsymmetry of the block kernel matrix. Recalling (CAR), the Pfaffian kernel is defined as the skewsymmetric matrix whose $(\ell, \ell')$-block is

$$ \mathbb{K}_{\ell\ell'} = \begin{pmatrix} \langle c_{i\ell}^* c_{i\ell'}^* \rangle_{\rho} & \langle c_{i\ell}^* c_{i\ell'} \rangle_{\rho} \\ \langle c_{i\ell} c_{i\ell'} \rangle_{\rho} - \delta_{i\ell i\ell'} & \langle c_{i\ell} c_{i\ell'} \rangle_{\rho} \end{pmatrix}. $$

Now, we define the kernel $\mathbb{K}(i_\ell, i_{\ell'}) = \mathbb{K}_{\ell\ell'}$. Then, the expression (53) matches the definition of the Pfaffian, i.e.,

$$ \langle \beta_1 \dots \beta_m \rangle_{\rho} = \mathrm{Pf} (\mathbb{K}(i, j))_{1 \le i, j \le m}. $$

By using Lemma 4.7, we obtain the expression (32) and (31) follows.

## B Gate details

We give here details about the implementation of a Givens rotation using elementary gates. Note first that zeroing out a matrix entry $y \in \mathbb{C} \setminus \{0\}$ can be done with a Givens rotation matrix as follows

$$ \begin{pmatrix} \cos \theta & e^{-i\phi} \sin \theta \\ -e^{i\phi} \sin \theta & \cos \theta \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} = \begin{pmatrix} r \\ 0 \end{pmatrix}, $$

with $\exp i\phi = x^*y/|xy|$, $\cos \theta = |x|/(|x|^2 + |y|^2)^{1/2}$ and $\sin \theta = |y|/(|x|^2 + |y|^2)^{1/2}$.

Let us list a few elementary gates to implement the corresponding operation on a pair of qubits. For ease of comparison with the online documentation of Qiskit, we denote in this section the Pauli matrices (14) by $X$, $Y$, and $Z$.
## B.2 Qiskit 0.42.1 implementation without controlled unitary gate.

The Givens rotation is implemented in Qiskit, as in Figure 1, thanks to a gate

$$R_{XX+YY}(2\theta, \phi - \pi/2) = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & c\theta & e^{-i\phi}s\theta & 0 \\ 0 & -e^{i\phi}s\theta & c\theta & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix}$$

where the 2-qubit parameterized $XX + YY$ interaction is defined as

$$\begin{aligned} R_{XX+YY}(\theta, \beta) &= (\exp(-i\beta Z/2) \otimes I) \exp\left(-i\frac{\theta}{2} \frac{X \otimes Y + Y \otimes X}{2}\right) (\exp(i\beta Z/2) \otimes I) \\ &= \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & c\theta/2 & -ie^{-i\beta}s\theta/2 & 0 \\ 0 & -ie^{i\beta}s\theta/2 & c\theta/2 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix}. \end{aligned}$$

This gate is implemented in Qiskit 0.42.1 as a composition of elementary gates containing 12 single-qubit gates and 2 CNOT gates. Also, note that some details in their gate definitions may vary from the one used in this section.

## C Sources of error in quantum computers

Different quantities are usually specified by the constructor to assess the efficiency of a quantum machine. First, the error rate (or readout error) is an estimated frequency of getting an undesired measurement (a “bit flip”) when measuring a state drawn independently from a fixed, user-defined prior distribution. Second, the de-excitation (or relaxation) of a qubit prepared in its excited state, and thus supposed to represent $|1\rangle$, is a natural physical process. It is usually understood to be caused by the coupling of the qubit to electromagnetic radiation or, more abstractly, its environment. The typical timescale of the relaxation processes is called the relaxation time, and usually denoted by $T_1$. Moreover, the presence of an environment is not only the source of relaxation (bit-flip error), but also an additional source of phase errors. This type of error destroys the quantum coherence properties in and between the qubits, i.e., it modifies the correlation structure in Boolean vectors built using Born’s rule (9), possibly to the point of making the qubits behave as simple classical bits. This decoherence process is usually characterized by a typical timescale denoted $T_2$. The output of any circuit with depth longer than $T_1$ or $T_2$ is likely to be strongly contaminated with the corresponding noises.

To give concrete numbers, the IBM 127 qubits Washington platform<fnref n="9" /> has a median $T_1$ (the median is over all qubits) of $T_1 = 95.22$ µs, and the median $T_2$ is $T_2 = 92.17$ µs. This means that, in practice, after about 100 µs, more than half of the qubits have either decohered (become classical bits) or have relaxed to their ground state. Moreover, the median readout error is $p_{\text{err, ro}} = 1.350 \times 10^{-2}$, and the median CNOT gate error is $p_{\text{err, CNOT}} = 1.287 \times 10^{-2}$. Ideally, characterizing the noise of all the gates

<footnote n="9">⁹See https://quantum-computing.ibm.com/services/resources?services=systems</footnote>
would require a quantum process tomography as well as quantum state tomography, to characterize the state of all qubits. In practice, this is out of reach for large systems. To circumvent this issue, Elben, Flammia, Huang, Kueng, Preskill, Vermersch, and Zoller, 2022 propose to partially characterize the noise with a scalable and robust algorithm called *randomized benchmarking*. This is how the cited numbers for IBM machines have been benchmarked.<fnref n="10" />

The long-term goal of the quantum computer race is to build a fault-tolerant quantum computer, based on quantum error corrections codes. These are techniques that build up on redundancies of the logical qubits to be robust to noise-induced errors according; see the so-called *threshold theorem* (Kitaev, 2003; Knill, Laflamme, and Zurek, 1998). As in the classical case, the required redundancy depends of the actual values of the error rates, which are to this day still too high to have a realistic implementation. In the meantime, other techniques are being developed to alleviate the influence or errors, either by directly eliminating errors on the hardware itself, like with so-called *spin echos* or *dynamic decoupling* techniques (Preskill, 1998), or by statistical postprocessing, with so-called *error mitigation techniques* (Cai, Babbush, Benjamin, Endo, Huggins, Li, McClean, and O'Brien, 2022).

## References

Ball, R. C. (2005). "Fermions without Fermion Fields". In: *Phys. Rev. Lett.* 95 (17), p. 176407. URL: https://link.aps.org/doi/10.1103/PhysRevLett. 95.176407.

Bardenet, R. and Hardy, A. (2020). "Monte Carlo with Determinantal Point Processes". In: *Annals of Applied Probability*. URL: https://doi.org/10.1214/19-AAP1504.

Bardenet, R. et al. (2022). "From point processes to quantum optics and back". In: *arXiv preprint arXiv:2210.05522*. URL: https://arxiv.org/abs/2210.05522.

Barthelmé, S., Tremblay, N., and Amblard, P.-O. (2023). "A Faster Sampler for Discrete Determinantal Point Processes". In: *Proceedings of The 26th International Conference on Artificial Intelligence and Statistics*. Vol. 206. Proceedings of Machine Learning Research. PMLR, pp. 5582-5592. URL: https://proceedings.mlr.press/v206/barthelme23a.html.

Belhadji, A., Bardenet, R., and Chainais, P. (2020). "A determinantal point process for column subset selection". In: *Journal of Machine Learning Research (JMLR)*. URL: http://jmlr.org/papers/v21/19-080.html.

Borodin, A. and Rains, E. M. (2005). "Eynard-Mehta theorem, Schur process, and their Pfaffian analogs". In: *Journal of statistical physics* 121.3, pp. 291-317.

Bouten, L., Van Handel, R., and James, M. R. (2007). "An introduction to quantum filtering". In: *SIAM Journal on Control and Optimization* 46.6, pp. 2199-2241. URL: https://arxiv.org/abs/math/0601741.

<footnote n="10">¹⁰https://qiskit.org/textbook/ch-quantum-hardware/randomized-benchmarking.html</footnote>
Bravyi, S.B. and Kitaev, A.Y. (2002). "Fermionic Quantum Computation". In: *Annals of Physics* 298.1, pp. 210–226. URL: https://www.sciencedirect.com/science/article/pii/S0003491602962548.

Bufetov, A. I., Deelan Cunden, F., and Qiu, Y. (2021). "Conditional measures for Pfaffian point processes: Conditioning on a bounded domain". In: *Annales de l'Institut Henri Poincaré, Probabilités et Statistiques* 57.2, pp. 856–871. DOI: 10.1214/20-AIHP1099. URL: https://doi.org/10.1214/20-AIHP1099.

Cai, Z., Babbush, R., Benjamin, S. C., Endo, S., Huggins, W. J., Li, Y., McClean, J. R., and O'Brien, T. E. (2022). "Quantum error mitigation". In: *arXiv preprint arXiv:2210.00921*. URL: https://arxiv.org/abs/2210.00921.

Cohen-Tannoudji, C., Diu, B., and Laloë, F. (1977). *Quantum mechanics; 1st ed*. New York, NY: Wiley.

Dean, D. S., Le Doussal, P., Majumdar, S. N., and Schehr, G. (2019). "Non-interacting fermions in a trap and random matrix theory". In: *Journal of Physics A: Mathematical and Theoretical* 52.14, p. 144006.

Demmel, J., Grigori, L., Hoemmen, M., and Langou, J. (2012). "Communication-optimal parallel and sequential QR and LU factorizations". In: *SIAM Journal on Scientific Computing* 34.1, A206–A239. URL: https://arxiv.org/abs/0806.2159.

Derezinski, M., Khanna, R., and Mahoney, M. W. (2020a). "Improved guarantees and a multiple-descent curve for Column Subset Selection and the Nystrom method". In: *Advances in Neural Information Processing Systems*. Vol. 33. Curran Associates, Inc., pp. 4953–4964. URL: https://tinyurl.com/nhj5wexa.

Derezinski, M., Liang, F., and Mahoney, M. (2020b). "Bayesian experimental design using regularized determinantal point processes". In: *Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics*. Vol. 108. Proceedings of Machine Learning Research. PMLR, pp. 3197–3207. URL: https://proceedings.mlr.press/v108/derezinski20a.html.

Derezinski, M. and Mahoney, M. W. (2021). "Determinantal point processes in randomized numerical linear algebra". In: *Notices of the American Mathematical Society* 68.1, pp. 34–45. URL: https://arxiv.org/abs/2005.03185.

Dereziński, M., Clarkson, K. L., Mahoney, M. W., and Warmuth, M. K. (2019). "Minimax experimental design: Bridging the gap between statistical and worst-case approaches to least squares regression". In: *Conference on Learning Theory*. PMLR, pp. 1050–1069.

Dierckx, B., Fannes, M., and Pogorzelska, M. (2008). "Fermionic quasifree states and maps in information theory". In: *Journal of Mathematical Physics* 49.3. URL: https://doi.org/10.1063/1.2841326.

Elben, A., Flammia, S. T., Huang, H.-Y., Kueng, R., Preskill, J., Vermersch, B., and Zoller, P. (2022). "The randomized measurement toolbox". In: *Nature Reviews Physics*, pp. 1–16. URL: https://arxiv.org/abs/2203.11374.

Fanuel, M., Schreurs, J., and Suykens, J.A.K. (2021). "Diversity Sampling is an Implicit Regularization for Kernel Methods". In: *SIAM Journal on Math-
matics of Data Science 3.1, pp. 280–297. URL: https://doi.org/10.1137/20M1320031.

Frerix, T. and Bruna, J. (Sept. 2019). “Approximating Orthogonal Matrices with Effective Givens Factorization”. In: *Proceedings of the 36th International Conference on Machine Learning*. Vol. 97. Proceedings of Machine Learning Research. PMLR, pp. 1993–2001. URL: https://proceedings.mlr.press/v97/frerix19a.html.

Golub, G. H. and Van Loan, C. F. (2012). *Matrix computations*. JHU Press.

Hall, B. C. (2010). *Lie Groups, Lie Algebras, and Representations: An Elementary Introduction*. Vol. 222. Graduate Texts in Mathematics. New York: Springer.

Hough, J. B., Krishnapur, M., Peres, Y., and Virág, B. (2006). “Determinantal processes and independence”. In: *Probability surveys*. URL: https://arxiv.org/abs/math/0503110.

IBM Quantum (2021). URL: https://quantum-computing.ibm.com/.

Jahangiri, S., Arrazola, J.M., Quesada, N., and Killoran, N. (2020). “Point processes with Gaussian boson sampling”. In: *Phys. Rev. E* 101 (2), p. 022134. URL: https://link.aps.org/doi/10.1103/PhysRevE.101.022134.

Jiang, Z., Sung, K. J., Kechedzhi, K., Smelyanskiy, V. N., and Boixo, S. (2018). “Quantum algorithms to simulate many-body physics of correlated fermions”. In: *Physical Review Applied* 9.4, p. 044036. URL: https://arxiv.org/abs/1711.05395.

Johansson, K. (2005). “Random matrices and determinantal processes”. In: *ArXiv Mathematical Physics e-prints*. eprint: math-ph/0510038.

Kargin, V. (2014). “On Pfaffian random point fields”. In: *Journal of Statistical Physics* 154.3, pp. 681–704. URL: https://doi.org/10.1007/s10955-013-0900-z.

Kassel, A. (2015). “Learning About Critical Phenomena from Scribbles and Sandpiles”. In: *ESAIM: Proc.* 51, pp. 60–73. URL: https://doi.org/10.1051/proc/201551004.

Kassel, A. and Lévy, T. (2022). “Determinantal probability measures on Grassmannians”. In: *Annales de l’Institut Henri Poincaré D*. URL: https://doi.org/10.4171/aihpd/152.

Kerenidis, I. and Prakash, A. (2022). “Quantum machine learning with subspace states”. In: *arXiv preprint arXiv:2202.00054*. URL: https://arxiv.org/abs/2202.00054.

Kitaev, A. Yu. (2003). “Fault-tolerant quantum computation by anyons”. In: *Annals of physics* 303.1, pp. 2–30. URL: https://doi.org/10.1016/S0003-4916(02)00018-0.

Kivlichan, I. D., McClean, J., Wiebe, N., Gidney, C., Aspuru-Guzik, A., Chan, G. K.-L., and Babbush, R. (2018). “Quantum Simulation of Electronic Structure with Linear Depth and Connectivity”. In: *Phys. Rev. Lett.* 120 (11), p. 110501. URL: https://link.aps.org/doi/10.1103/PhysRevLett.120.110501.

Knill, E., Laflamme, R., and Zurek, W. H. (1998). “Resilient quantum computation: error models and thresholds”. In: *Proceedings of the Royal Soci-
ety of London. Series A: Mathematical, Physical and Engineering Sciences
454.1969, pp. 365-384. URL: https://www.jstor.org/stable/i203157.

Koshida, S. (2021). "Pfaffian Point Processes from Free Fermion Algebras: Per-
fectness and Conditional Measures". In: *SIGMA. Symmetry, Integrability
and Geometry: Methods and Applications* 17, p. 008. URL: https://doi.
org/10.3842/SIGMA.2021.008.

Kulesza, A. and Taskar, B. (2011). "Learning Determinantal Point Processes". In: *Proceedings of the Twenty-Seventh Conference on Uncertainty in Artificial Intelligence*. UAI'11. Barcelona, Spain: AUAI Press, pp. 419-427. ISBN: 9780974903972. URL: https://arxiv.org/abs/1202.3738.

(2012). "Determinantal point processes for machine learning". In: *Foundations and Trends in Machine Learning*. URL: http://dx.doi.org/10.1561/2200000044.

Kyng, R. and Song, Z. (2018). "A Matrix Chernoff Bound for Strongly Rayleigh Distributions and Spectral Sparsifiers from a few Random Spanning Trees". In: *2018 IEEE 59th Annual Symposium on Foundations of Computer Science (FOCS)*, pp. 373-384. URL: https://arxiv.org/abs/1810.08345.

Lavancier, F., Møller, J., and Rubak, E. (2015). "Determinantal point process models and statistical inference". In: *Journal of the Royal Statistical Society, Series B*. B. URL: https://www.jstor.org/stable/24775312.

Lyons, R. (2003). "Determinantal probability measures". In: *Publications Mathematiques de l'IHES* 98, pp. 167-212. URL: https://doi.org/10.1007/s10240-003-0016-0.

Macchi, O. (1972). "Processus ponctuels et coincidences – Contributions à l'étude théorique des processus ponctuels, avec applications à l'optique statistique et aux communications optiques". PhD thesis. Université Paris-Sud.

— (2017). *Point processes and coincidences – Contributions to the theory, with applications to statistical optics and optical communication, augmented with a scholion by Suren Poghosyan and Hans Zessin*. Walter Warmuth Verlag.

Moore, G. W. (2014). "Quantum symmetries and compatible hamiltonians". URL: http://www.physics.rutgers.edu/gmoore/QuantumSymmetryBook.pdf.

Nielsen, M. A. (2005). *The fermionic canonical commutation relations and the Jordan-Wigner transform*. Tech. rep. The University of Queensland. URL: https://futureofmatter.com/assets/fermions_and_jordan_wigner.pdf.

Nielsen, M. A. and Chuang, I. L. (2010). *Quantum computation and quantum information*. 10th anniversary edition. Cambridge University Press. URL: https://doi.org/10.1017/CBO9780511976667.

Olshanski, G. (2020). "Determinantal point processes and fermion quasifree states". In: *Communications in Mathematical Physics* 378.1, pp. 507-555. URL: https://doi.org/10.1007/s00220-020-03716-1.

Ortiz, G., Gubernatis, J. E., Knill, E., and Laflamme, R. (2001). "Quantum algorithms for fermionic simulations". In: *Physical Review A* 64.2, p. 022319. URL: https://doi.org/10.1103/PhysRevA.64.022319.
Pemantle, R. (1991). "Choosing a Spanning Tree for the Integer Lattice Uniformly". In: *The Annals of Probability* 19.4, pp. 1559-1574. URL: http://www.jstor.org/stable/2244527.

Poinas, A. and Bardenet, R. (2023). "On proportional volume sampling for experimental design in general spaces". In: *Statistics and Computing* 33.1, p. 29. URL: https://doi.org/10.1007/s11222-022-10115-0.

Preskill, J. (1998). "Lecture notes for physics 229: Quantum information and computation". In: *California Institute of Technology* 16.1, pp. 1-8. URL: http://theory.caltech.edu/~preskill/ph229/.

Qiskit contributors (2023). *Qiskit: An Open-source Framework for Quantum Computing*. DOI: 10.5281/zenodo.2573505.

Rains, E. M. (2000). "Correlation functions for symmetrized increasing sequences". In: *arXiv preprint math/0006097*. URL: https://arxiv.org/abs/math/0006097.

Rebentrost, P., Steffens, A., Marvian, I., and Lloyd, S. (2018). "Quantum singular-value decomposition of nonsparse low-rank matrices". In: *Phys. Rev. A* 97 (1), p. 012327. URL: https://link.aps.org/doi/10.1103/PhysRevA.97.012327.

Sameh, A. H. and Kuck, D. J. (1978). "On stable parallel linear system solvers". In: *Journal of the ACM (JACM)* 25.1, pp. 81-91. URL: https://doi.org/10.1145/322047.322054.

Schnyder, A.P., Ryu, S., Furusaki, A., and Ludwig, A.W.W. (2008). "Classification of topological insulators and superconductors in three spatial dimensions". In: *Phys. Rev. B* 78 (19), p. 195125. URL: https://link.aps.org/doi/10.1103/PhysRevB.78.195125.

Soshnikov, A. (2000). "Determinantal random point fields". In: *Russian Mathematical Surveys* 55, pp. 923-975. URL: https://dx.doi.org/10.1070/RM2000v055n05ABEH000321.

— (2003). "Janossy Densities. II. Pfaffian Ensembles". In: *Journal of statistical physics* 113, pp. 611-622. URL: https://doi.org/10.1023/A:1026077020147.

Terhal, B. M. and DiVincenzo, D. P. (2002). "Classical simulation of noninteracting-fermion quantum circuits". In: *Phys. Rev. A* 65 (3), p. 032325. URL: https://link.aps.org/doi/10.1103/PhysRevA.65.032325.

Tremblay, N., Barthelmé, S., and Amblard, P.-O. (2019). "Determinantal Point Processes for Coresets". In: *Journal of Machine Learning Research* 20.168, pp. 1-70. URL: http://jmlr.org/papers/v20/18-167.html.

Tremblay, N., Barthelmé, S., Usevich, K., and Amblard, P.-O. (2023). "Extended L-ensembles: a new representation for Determinantal Point Processes". In: *The Annals of Applied Probability* 33.1, pp. 613-640.

Verstraete, F. and Cirac, J.I. (2005). "Mapping local Hamiltonians of fermions to local Hamiltonians of spins". In: *Journal of Statistical Mechanics: Theory and Experiment* 2005.09, P09012. URL: https://dx.doi.org/10.1088/1742-5468/2005/09/P09012.

Wecker, D., Hastings, M. B., Wiebe, N., Clark, B. K., Nayak, C., and Troyer, M. (2015). "Solving strongly correlated electron models on a quantum computer".
puter”. In: *Physical Review A* 92.6, p. 062318. URL: https://link.aps.org/doi/10.1103/PhysRevA.92.062318.

Wilson, D. B. (1996). “Generating Random Spanning Trees More Quickly than the Cover Time”. In: *Proceedings of the Twenty-Eighth Annual ACM Symposium on Theory of Computing*. STOC ’96. Philadelphia, Pennsylvania, USA: Association for Computing Machinery, pp. 296–303. URL: https://doi.org/10.1145/237814.237880.