Classical and Quantum Gravity PAPER • OPEN ACCESS Parameter estimation of gravitational waves with a quantum metropolis algorithm To cite this article: Gabriel Escrig et al 2023 Class. Quantum Grav. 40 045001   View the article online for updates and enhancements. You may also like From single layer to multilayer networks in mild cognitive impairment and Alzheimer’s disease Ignacio Echegoyen, David López-Sanz, Fernando Maestú et al. - PDRs4All: A JWST Early Release Science Program on Radiative Feedback from Massive Stars PI Team:, Olivier Berné, Émilie Habart et al. - Origin of the magnetic transition at 100 K in -Fe2O3 nanoparticles studied by x-ray absorption fine structure spectroscopy J López-Sánchez, A Muñoz-Noval, C Castellano et al. - This content was downloaded from IP address 147.96.28.131 on 16/02/2023 at 18:32 https://doi.org/10.1088/1361-6382/acafcf /article/10.1088/2632-072X/ac3ddd /article/10.1088/2632-072X/ac3ddd /article/10.1088/2632-072X/ac3ddd /article/10.1088/1538-3873/ac604c /article/10.1088/1538-3873/ac604c /article/10.1088/1538-3873/ac604c /article/10.1088/1361-648X/aa904b /article/10.1088/1361-648X/aa904b /article/10.1088/1361-648X/aa904b /article/10.1088/1361-648X/aa904b /article/10.1088/1361-648X/aa904b /article/10.1088/1361-648X/aa904b /article/10.1088/1361-648X/aa904b /article/10.1088/1361-648X/aa904b Classical and Quantum Gravity Class. Quantum Grav. 40 (2023) 045001 (15pp) https://doi.org/10.1088/1361-6382/acafcf Parameter estimation of gravitational waves with a quantum metropolis algorithm Gabriel Escrig1,∗, Roberto Campos1,2, Pablo A M Casares1 and M A Martin-Delgado1,3 1 Departamento de Física Teórica, Universidad Complutense de Madrid, Madrid, Spain 2 Quasar Science Resources, SL, Madrid, Spain 3 CCS-Center for Computational Simulation, Universidad Politécnica de Madrid, Madrid, Spain E-mail: gescrig@ucm.es Received 22 September 2022; revised 21 November 2022 Accepted for publication 3 January 2023 Published 20 January 2023 Abstract After the first detection of a gravitational wave in 2015, the number of suc- cesses achieved by this innovative way of looking through the Universe has not stopped growing. However, the current techniques for analyzing this type of events present a serious bottleneck due to the high computational power they require. In this article we explore how recent techniques based on quantum algorithms could surpass this obstacle. For this purpose, we propose a quantiz- ation of the classical algorithms used in the literature for the inference of grav- itational wave parameters based on the well-known quantum walks technique applied to a Metropolis–Hastings algorithm. Finally, we develop a quantum environment on classical hardware, implementing ametric to compare quantum versus classical algorithms in a fair way. We further test all these developments in the real inference of several sets of parameters of all the events of the first detection period GWTC-1 and we find a polynomial advantage in the quantum algorithms, thus setting a first starting point for future algorithms. Keywords: gravitational waves, quantum walks, Metropolis–Hastings algorithms (Some figures may appear in colour only in the online journal) ∗ Author to whom any correspondence should be addressed. Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. 1361-6382/23/045001+15$33.00 © 2023 The Author(s). Published by IOP Publishing Ltd Printed in the UK 1 https://doi.org/10.1088/1361-6382/acafcf https://orcid.org/0000-0003-2881-085X https://orcid.org/0000-0002-2527-4177 https://orcid.org/0000-0001-5500-9115 https://orcid.org/0000-0003-2746-5062 mailto:gescrig@ucm.es http://crossmark.crossref.org/dialog/?doi=10.1088/1361-6382/acafcf&domain=pdf&date_stamp=2023-1-20 https://creativecommons.org/licenses/by/4.0/ https://creativecommons.org/licenses/by/4.0/ Class. Quantum Grav. 40 (2023) 045001 G Escrig et al 1. Introduction After 100 years of its prediction and an international effort spanning several decades, on 14 September 2015, the two advanced Laser Interferometer Gravitational-Wave Observatory (LIGO) detectors simultaneously detected the first signal of a gravitational wave (GW) [1]. This made it possible to open a new window through which to look at the Universe, com- plementary to the observation of the electromagnetic spectrum. This new method promises the observation of astrophysical events never seen before, but its processing is costly. The computational bottleneck arises in the Bayesian inference of the parameters of the sources that generated the gravitational radiation and will grow when a larger number of parameters, higher simulation accuracy and more events to analyze are obtained due to the improvement of the detectors between observation periods. Historically, the computational advances have been very large. Over the last 20 years, dif- ferent types of algorithms have been proposed (see [2, 3]), but those based on the formalism of Bayesian inference have stood out from the rest, making Bayesian parameter estimation the language of GW astronomy. Among these algorithms the most important are the Metropolis– Hastings algorithms based on Monte Carlo Markov Chains [4]. In addition, a large number of software libraries have been developed (see [5–7]), making it clear that the field of GW is growing. Currently, the computational analysis cost is very large and will get larger for several reas- ons. First, in the near future, and with the continuous improvement of the detectors, many more black hole and neutron star mergers will be measured. For instance, in the next observing period O4, the LIGO-Virgo-KAGRA collaboration is expected to observe about 79+89 −44 binary black hole (BBH) events and about 10+52 −10 binary neutron star (BNS) events over one year [8]. Since neutron stars are more complex systems than black holes, they require more parameters to model and larger computational times. Specifically, the estimated time to analyze a BBH merger event is on the order of days, while the estimated time to analyze a BNS merger event is on the order of weeks. As a result, the overall analysis time for the next run could be longer than one year. Second, with the addition of new observatories such as KAGRA and in the longer term the space interferometer LISA, there is hope to be able to observe different events such as supernovae explosions, the first candidate considered for detecting GWs, or massive black holes at galactic centers. These highlights the importance of developing new algorithms that accelerate the ana- lysis of these astrophysical events. In this manuscript, we review the most commonly used algorithms for the analysis of GWs and explore how quantum walk-based algorithms can be used to search for the set of parameters that best fits the data obtained in the detectors. In par- ticular, we compare the performance of a quantum Metropolis–Hasting algorithm against its classical counterpart using the data from the events detected so far, which can be found in the PyCBC catalog at [9]. Quantum walks are a cornerstone family of quantum algorithms that exhibit a quadratically larger eigenvalue gap than the associated randomwalk. Since the inverse of the eigenvalue gap is often amultiplicative factor in the complexity of the algorithm, quantumwalks are capable to outperform their classical counterparts, as we will explain. Specifically, quantum Metropolis– Hastings algorithms exhibiting this quadratic advantage have been developed [10–14], and have also been proposed as a method to tackle hard optimization problems [15–17]. Although the use of quantum walks is quite extensive throughout the literature [18], to the best of our knowledge there are no previous analysis of the use of this type of methods in the study of gravitational wave events, only in their detection [19]. Thus, we hope this work will be of interest to the gravitational astronomical community. 2 Class. Quantum Grav. 40 (2023) 045001 G Escrig et al The rest of the manuscript is organized as follows. In the next section 2 we introduce the basics of classical and quantum Metropolis–Hastings algorithms, and how they compare with each other. Then, in section 3 we describe our simulation setup, whose results are explained in 4. We end the article with the conclusions and future work in section 5. 2. Bayesian inference: classical and quantum Metropolis–Hastings algorithms The standard formalism for tackling the task of estimating the unknown parameters of a grav- itational wave source classically is the formalism of Bayesian inference based on Markov Chain Monte Carlo methods [20], as previously stated, despite the fact that several algorithms have been proposed in the literature. Among all the algorithms of this type, the Metropolis– Hastings algorithm stands out, and will consequently be the topic of study of this work, both in its classical and quantum versions. 2.1. Classical Metropolis–Hastings The formalism of Bayesian inference describes the available knowledge about a parameter which we want to estimate from observations, θ, as a probability density that we do not know a priori. In this framework, the Bayes’ Theoremwill be ourmain tool, where a prior probability distribution p(θ) is updated as new data d is received from the experiment to give the posterior distribution p(θ|d) [21], p(θ|d) = p(θ)p(d|θ) p(d) . (1) The models we will explore have many parameters, in particular the models used for the BBH and BNS events contain 15 parameters [9], among which the distance to the source, the masses, the spin, the coalescence time or the position in the sky can be highlighted. We will indicate collectively this set of parameters as θ = {θ1,θ2, . . . ,θN}. Thus, in order to describe the collective knowledge about all parameters and their relationships, the joint probability distri- bution on the multidimensional space p(θ|d)will be used, being able to recover the probability for a specific parameter by marginalizing over the other unwanted parameters. Since the evidence, p(d), only depends on the data obtained and is common for the whole parameter space, it can be taken as a normalization factor [22], so that the posterior probability reduces to p(θ|d)∝ p(θ)L(d|θ), (2) where we have introduced the likelihood function as the conditional probability L(d|θ) := p(d|θ). The likelihood function is something that we choose, it is a description of the meas- urement. Thus, it is appropriate to choose a right likelihood function with the noise in our data. For this purpose, it is common in the literature of gravitational wave astronomy to assume a Gaussian noise in the detectors [23], so that the Gaussian-noise likelihood function has the form: L(d|θ)∝ exp ( −1 2 ˆ ∞ 0 |d̃( f)− µ̃(θ, f)|2 σ2( f) df ) . (3) All these functions are computed in the Fourier domain, with f denoting frequencies, where d̃ represents the data set obtained in the usual way at the detector, µ̃ the waveform for a cer- tain relativistic numerical model calculated from the parameter set θ, and σ the detector noise 3 Class. Quantum Grav. 40 (2023) 045001 G Escrig et al (spectral density). It is interesting to note that despite having an integral at all possible frequen- cies, in practice we will limit this interval to those frequencies at which ground-based detectors work, also removing some frequencies such as the harmonics of domestic electricity, ensuring that at all these frequencies σ is non-zero, and therefore there cannot be any divergence in (3). Since the likelihood function is equal to a negative exponential function, we can define an effective energy function from the argument of that exponential, E(θ) := 1 2 ˆ ∞ 0 |d̃( f)− µ̃(θ, f)|2 σ2( f) df. (4) This energy is a measure of how different is the available data, d̃, from our relativistic numerical model of gravitational wave, µ̃, calculated from parameters θ. The use of numerical models of gravitational waveforms is due to the complexity of finding analytical solutions for some GW sources to the Einstein equations. For our purposes, we will use a single numerical model depending on the event (BBH or BNS). Then, the closer the parameters are to the source value, the closer the wave predicted by the relativistic numerical model resembles the data obtained, the lower the energy and therefore the higher the value of the likelihood. The goal in the parameter estimation of GWs is to find the set of parameters that best fits the data obtained from the detectors. For this purpose, the Metropolis–Hastings algorithm is the standard and best-performing technique for parameter inference in gravitational wave astronomy, which performs a random walk over the gravitational wave model parameter space Ω. This algorithm is specially designed to reach the equilibrium state as quickly as possible, that is, to reach the probability distribution ϕ such that Wϕ= ϕ, where W is the transition matrix of our stochastic process. To reach this state we will start from an initial, usually uni- form, distribution ϕ(0) and a random transition will be proposed to obtain the next state. In this way: (a) If at the current iteration i, the chain has a location in the parameter space of θi, we propose a random transition∆θi to the new parameter state θi+1 = θi+∆θi. (b) We now calculate the function likelihood found in (3) for the new state, and accept this transition with probability min ( 1, p(θi+1) p(θi) ( L(d|θi+1) L(d|θi) ) 1 T ) . (5) The Metropolis–Hastings algorithm works by iterating over these steps. It is important to note that the variable T introduced in the algorithm is what is known in the literature as tem- pering [4], playing a role analogous to temperature, so that by varying this parameter we can obtain different speeds of convergence to the optimal parameters. We can go a step further and introduce simulated annealing, so we can start the optimization with a high temperature and lowering it with some annealing scheme, β, so that at first the algorithm is more likely to accept larger steps and gradually accept smaller steps until finally converging to the optimal set parameters. Note that if the prior probability is uniform, we only need the function likelihood to decide whether to accept the next step in the random walk. We can finally formulate in matrix terms the random walk, and thus leave an algorithm ready for quantization. To do this, starting from a state i, the Metropolis–Hastings algorithm proposes a random change to one of the configurations, j, connected to i. We will call T ij the 4 Class. Quantum Grav. 40 (2023) 045001 G Escrig et al probability of such a proposal. If we have a uniform prior, it is customary to choose a uniform probability distribution, Tij = { 1 N there is a transition that connects i with j, 0 if there isn’t, (6) forN the number of possible transitions from state i. Alternatively, we could choose Tij ∝ p(θj) p(θi) . Then, this change is accepted with a probability Aij =min ( 1, ( L(d|θj) L(d|θi) )β ) =min ( 1,eβ[E(θj)−E(θi)] ) , (7) resulting in an overall transition probability for a given step of i→ j: Wij = { TijAij if i 6= j, 1− ∑ k ̸=jTkjAkj if i= j. (8) As such, we have explained how the Metropolis–Hastings algorithm can be used to solve our Bayesian inference problem. 2.2. Quantum Metropolis–Hastings Similarly, the objective for the Metropolis–Hastings algorithm is sampling the Boltzmann dis- tribution ρβ = |ϕβ〉〈ϕβ |= Z−1(β) ∑ θ∈Ω e−βE(θ)|θ〉〈θ|. E(θ) is the energy function defined above for gravitational wave parameters θ, β is the annealing schedule, playing the role of the inverse of the temperature and Z(β) = ∑ θ∈Ω e−βE(θ) is the partition function, a normaliza- tion factor. Here, β plays an important role in finding the lowest energy state because if it is large enough, only the configuration with the lowest possible energy will appear with a high probability when sampling the state. 2.2.1. Construction of the quantum Metropolis–Hastings algorithm. There are several ways of constructing quantum walks, but the best known for its results for this purpose is the one proposed by Szegedy [24, 25]. Szegedy’s algorithm is based on two key features. First, it duplicates the registers encoding the state, so the operation space is twice the original Hil- bert space. Given the acceptance probabilities Wij = TijAij, defined previously in (8), for the transition from state i to state j, one defines the unitary operator U| j〉|0〉 := | j〉 ∑ i∈Ω √ Wji|i〉= | j〉|pj〉. (9) The second key feature is that Szegedy quantum walk operator consists of two rotations similar to those performed in the well-known Grover’s algorithm [26, 27]. Taking R0 := 1− 2Π0 = 1− 2(1⊗ |0〉〈0|) (10) the reflection over the state |0〉 in the second Hilbert subspace, and S the swap gate that exchanges both subspaces, we define the quantum walk step as W := U†SUR0U †SUR0. (11) The first R0 represents the first rotation, while U†SUR0U†SU is the second. The Metropolis–Hastings algorithm we are going to use is a modification of Szegedy’s quantum walk [14], replacing the bipartite space with the use of a coin. As a result, we will substitute the use of W with a coin version of the evolution operator, UW . This operator will 5 Class. Quantum Grav. 40 (2023) 045001 G Escrig et al act on 3 quantum registers: a register for the system, which will indicate the current state of the system, a register encoding the transition, and a register for the Boltzmann coin encoding probabilities Aij. We will denote them respectively by |x〉S|m〉M|b〉C. The operator UW is equi- valent to half of the Szegedy operatorW in (8) under conjugation by the isomorphism operator Y, which maps the states in the second register of the original quantum walk to the transitions needed to reach them [14]: Y†(R0U †SU)Y= UW (12) with Y : |x〉S|y〉S ′ 7→ |x〉S|m〉M such that m(x) = y. (13) The coined version operator UW is then constructed as UW = RV†B†FBV. (14) In this equation, V prepares a superposition in the register |·〉M over all possible transitions that the state may undergo in the next step. B rotates the coin qubit |·〉C to have an amplitude of |1〉C corresponding to the probability of acceptance Aij indicated in (7). F changes the register of the system |·〉S to the new configuration, conditioned by the value of |·〉M and |·〉C. Finally, R is nothing more than a reflection on the |0〉MC state. If additional auxiliary registers |·〉A are present, then the rotation R will reflect on those registers too. It is straightforward to construct these operators from the previous definitions. As V pre- pares the possible transitions in the register |·〉M, these will be given by T ij. The operator V consequently prepares a uniform superposition V : |0〉M → N− 1 2 ∑ m |m〉M. (15) Then B prepares the coin qubit so that the probability of being in state |1〉C corresponds to the probability of acceptance, Aij B : |x〉S|m〉M|0〉C → |x〉S|m〉M (√ 1−Ax,m(x)|0〉C+ √ Ax,m(x)|1〉C ) . (16) F will update the system registry to the new configuration, F : |x〉S|m〉M|b〉C → |b ·m(x)+ (1− b) · x〉S|m〉M|b〉C. (17) And finally R is the reflection operator on the state |0〉MC R : |0〉M|0〉C →−|0〉M|0〉C | j〉M|b〉C → | j〉M|b〉C if ( j,b) 6= (0,0). (18) There have been several proposals to quantize the Metropolis–Hastings algorithm using these quantumwalks. These are often based on a slow evolution of |θt〉 → |θt+1〉 using quantum phase estimation to project the state into the stationary state of the corresponding Markov Chain. Since the eigenvalue gap of the quantum walk operator is quadratically larger than the classical one, this technique often allows achieving the result with a quadratic quantum advantage [10–13]. However, classically the Metropolis–Hastings algorithm is usually used with a rapidly varying annealing schedule, where no complexity-theoretical guarantees can be 6 Class. Quantum Grav. 40 (2023) 045001 G Escrig et al provided. As such, instead of relying on the complex quantum phase estimation procedure, we opt for a simpler and heuristic unitary procedure introduced first in [14]. It consists of executing L times the evolution operator on an initial state U |ψ(L)〉 := UWL . . .UW2UW1 |ϕ(0)〉 (19) where t= 1,2, . . . ,L will also define the annealing schedule, for the chosen values of β(t) at each step. This is, in some ways, the simplest approach one could think of to quantize the Metropolis–Hastings algorithm, as it is very similar to the classical way of performing many steps of the random walk. Thus, the procedure of our quantum algorithm is clear. Given an initial state |ϕ(0)〉 which represents a uniform superposition, the algorithm consists of applying the evolution operator defined in (14) iteratively, obtaining the state (19). Hopefully, the probability of measuring the set of parameters that minimizes the energy, |θ∗〉, will increase until it is close to 1. So, having introduced the quantumMetropolis–Hastings algorithm, we will next explore how to compare it with its classical counterpart. 2.3. Figure of merit: the total time to solution Given the two algorithms, the classical algorithm and its quantum version, we are now inter- ested in finding a tool capable of comparing both algorithms in a fair way. When we talk about a good algorithm for parameter inference in gravitational wave astronomy we are interested in two aspects: on the one hand, we want a model that finds the maximum likelihood parameters with a high probability. On the other, we want this algorithm to be fast. For example, comput- ing all possible parameter configurations would be quite accurate, but extremely expensive. When a random walk is employed to minimize some function E(θ), the minimum θ∗ is reached only with probability p< 1. Starting from a certain distribution q(θ) and apply- ing the W walk sequentially t iterations, the probability of success is p(t) = (W tq)(θ∗). To increase this probability to a constant value of 1− δ, it is sufficient to repeat the procedure L= log(1−δ) log(1−p(t)) times. Thus, we can find a natural metric given these conditions as the number of quantum walk steps t times the number of attempts L, introducing the figure of merit called TTS [14], defined in other terms as the expected average time it would take for the algorithm to find the solution if we could repeat the procedure in case of failure: TTS(t) = t log(1− δ) log(1− p(t)) (20) where t ∈ N is the number of quantum/classical walk steps performed in one execution of the quantum/classical Metropolis–Hastings algorithm, p(t) is the probability of hitting the correct parameter set after those iterations in each run of the algorithm, and δ is an arbitrary target probability of success, which we set to 0.9. This metric represents the tradeoff between longer walks and the increased probability of success. Longer walks may achieve a higher probability of success and thus be repeated fewer times, but increasing the duration t of the walk beyond a certain point could have a negligible impact on its probability of success p(t) and may not be worth it. Thus we define the minimum TTS over the total of all iterations as min(TTS) = mintTTS(t), thus finding a way to compare quantum and classical walks by comparing the lowest achieved values of the TTS, i.e. comparing the min(TTS) for different simulations. It remains to explain how to calculate p(t) in each of the algorithms. For the quantum case, the heuristic method of using the quantum walk explained in the previous section is proposed. 7 Class. Quantum Grav. 40 (2023) 045001 G Escrig et al Startingwith a |ϕ(0)〉 state, it consists of applying the quantumwalk operatorUW j sequentially, obtaining the state (19). Thus, the probability of obtaining the searched state is: p(θ = θ∗) = |〈θ∗|UWL . . .UW 2UW1 |ϕ(0)〉|2. (21) So the construction of the TTS is straightforward, obtaining for the Lth iteration of the algorithm the expression TTS(L) = L log(1− δ) log(1− |〈θ∗|UWL . . .UW2UW1 |ϕ(0)〉|2) . (22) Similarly, for the classical algorithm, we estimate the probability of obtaining θ∗ after t steps of the Metropolis–Hastings algorithm. 3. Simulation setup In contrast to standard Bayesian inference analysis, the objective of this paper is to benchmark the usefulness of both classical and quantum Metropolis–Hastings algorithms more than to obtain the most likely parameter values. To carry out this task, the fairest way to compare the two algorithms would be to run the classical Metropolis–Hastings algorithm on classical hardware and the quantum Metropolis–Hastings algorithm on quantum hardware. However, in the current era, only noisy intermediate-scale quantum computers are available, and these are only capable of running simple algorithms. Current quantum computers are not capable of executing the kind of algorithms proposed in this paper. This can be seen just by studying the cost of implementing the oracle (14). With the amount of qubits needed for each component of the Walk Operator [15] and the amount of parameters to be estimated with the desired precision, current quantum computers with a few hundred qubits would be unable to do it. Moreover, this has also been proven heuristically in [16]. Under the same oracle construction as in this work, they have only managed to execute the calculation of two angles at the same time with only 1 bit accuracy and with a high amount of noise, so at this time only a simulation is achievable. In this work, we have used the available data for GW mergers from the first GW catalog GWTC-1 [29]. The data and the posterior distributions inferred by [9] can be found in the data availability statement section. Since the quantum simulators are restricted in the number of qubits that we can simulate, we will fix most of the parameters to their inferred value, and only run the Metropolis–Hastings algorithm in a small subset of them. This limitation is due to the use of quantum simulators, that scale exponentially on classical computers, but it should not be an issue once we have early fault-tolerant quantum algorithms available. Likewise, the computation of the energy deltas going into the likelihood function is abstracted away to a readout table, which simplifies the quantum circuits required by our algorithm. While this requires exhaustively precomputing these values in advance, this is again a reflection of the limitations of current quantum simulators. When provided with a fault-tolerant quantum computer, any model to compute these energies can be implemented in a quantum computer at a constant multiplicative overhead. This is because quantum computers can perform classical calculations using exactly the same algorithms and therefore have the same scaling, differing only in a constant multiplicative overhead, which will be very high because an operation on a quantum computer takes ten orders of magnitude longer than on a classical computer [30]. The quantum advantage does not happen in this energy calculation. Rather, the polynomial speedup associated with quantum walks will be reflected in the number of steps of the algorithm, the factor we measure with the TTS introduced in the previous section. 8 Class. Quantum Grav. 40 (2023) 045001 G Escrig et al Figure 1. Flowchart of the algorithm implemented. This algorithm consists of three stages: a first stage where, once the parameters to be inferred and their number of pre- cision bits have been chosen, the energies are calculated with the help of the GW lib- rary PyCBC [6]. A second stage, where the two Metropolis–Hastings algorithms are executed, the quantum algorithm being programmed from the Qiskit library [28], and a third stage where the results are stored and the values for the chosen figure of merit are calculated. Finally, the value of min total time to solution (TTS) is obtained for each of the classical/quantum parts for later comparison. The pipeline of our library is explained in figure 1, which we will explain in the following sections. 3.1. The energy calculator module The algorithm starts once we have chosen the set of parameters we want to infer from the con- figuration file and the number of precision bits with which we want to discretize the parameter space. For example, for a parameter such as the inclination, which takes values in the interval [0,π], taking 1 bit of precision means that this parameter can only take the values {0,π}, while a precision of 2 bits indicates discretization to the values {0, π3 , 2π 3 ,π}. Then, the real LIGO data of all GW events from GWTC-1 are downloaded. For each of these events, all possible combinations of the chosen parameters and number of precision bits are calculated and the energy of each of these combinations is stored. We have used the PyCBC gravitational wave library to obtain the energy of each combination using a model that marginalizes the polarization, in the same way as [6], so we are left with a total of 14 parameters to vary for the BBH events listed in GWTC-1. We further need the range of possible values for these parameters, which we extract from the configuration files of the repository of the catalog [9]. Thus, this leaves a clear scheme for calculating the energy for any set of parameters for any BBH event for all observation periods. 3.2. The Metropolis–Hastings solver module The Metropolis–Hastings solver module makes use of the software tool QMS [17]. It receives a description of the problem as a list of tuples (state, energy) and it calculates the difference 9 Class. Quantum Grav. 40 (2023) 045001 G Escrig et al of energy between positions, called deltas. QMS tries, simultaneously, to find the minimum energy state and to avoid getting stuck in a local minimum energy state. The main functionality of QMS is the quantum search to find the minimum energy state. However, it also serves as a test-bed for classical–quantum algorithm comparison. For this reason, QMS performs a classical search, in order to compare both results, classical and quantum. This software tool can return the resulting amplitude probabilities, as real quantum hardware does, or just sample from this target distribution. QMS executes the quantum search in a quantum simulator running in a classical computer. The quantum simulator used in this work is QASM Simulator of Qiskit [28]. Delta energies are introduced in the algorithm as a precalculated oracle, or table, queried via a serial QRAM. 3.3. The output processor module Finally, the last module consists of generating all the results obtained by the Metropolis– Hastings algorithms for each of the GW events. As we have mentioned, our purpose is not to develop a quantum algorithm that estimates the parameters of the GW sources, but to demon- strate whether it is possible to achieve a quantum advantage in GW inference, thus creating a starting point for possible future developments. For this, we have previously introduced the figure of merit min(TTS). Such will be the output of this module of the algorithm. Thus, in this module the computation of the TTS will be carried out, storing the smallest value in TTS, the min(TTS), in order to later be able to study the behavior of this figure of merit as the complexity of the simulations increases. This complexity can be increased either by increasing the number of precision bits or by increasing the set of parameters to be inferred. We want to study the behavior of the classical min(TTS) versus the quantum min(TTS) as we increase the complexity of the simulation in order to find the exponent that will indicate the possible quantum advantage. In logarithmic scales, this reduces to finding the slope of the best linear least square fit. There will be a quantum advantage when the quantum min(TTS) grows as a power of the classical min(TTS) with an exponent less than 1. For example, a quadratic quantum advantage implies that this exponent is 0.5. 4. Simulation results In this section, we will test the classical and quantum Metropolis–Hastings algorithms. To achieve that goal we will study all the 11 events found in the first GW catalogue, the GWTC- 1. In this catalog only BBH events have been detected [29], so taking into account that we use the PyCBC model that marginalizes the polarization in order to reproduce the results of the catalogue of [9], we are left with 14 parameters available for the BBH events. A good starting point is to make inference in pairs of parameters that make physical sense. Among these two-parameter sets we may highlight the two masses, reparameterized to chirp mass and mass ratio. Also the dimensionless spin-magnitude of the two black holes. In addi- tion, and to complete the set of parameters, the luminosity distance, reparametrised to the comoving volume together with the inclination and the coalescense phase, are also two inter- esting pair of parameters. As what we want to study is the slope, and this is obtained by increasing the complexity of the simulation, we will represent the previous pairs of parameters by increasing the num- ber of precision bits. The results of these simulations will then be the plot of the quantum min(TTS) versus the classical one on a logarithmic scale to obtain in a simple and visual 10 Class. Quantum Grav. 40 (2023) 045001 G Escrig et al Figure 2. Comparison of the minimum value of the classical and quantumTTS achieved for the 2-parameter inference simulations. Top left the source-frame chirp mass and the mass ratio have been inferred. Top right the comoving volume and the inclination have been inferred. Below left the dimensionless spin-magnitude of the larger and smaller object have been inferred. Below right the comoving volume and the coalescence phase have been inferred. All results have been obtained for 3–6 bits of precision with a con- stant value of β. way the exponent. Starting with the inference of two parameters at a time, this is shown in figure 2. The classical and quantum min(TTS) are on the x-axis and y-axis respectively on a logarithmic scale, so a slope less than one would imply an advantage of the quantum algorithms. Conversely, a slope greater than one would imply an advantage of the classical algorithms. To distinguish this behavior, in all these graphs there is a dashed gray line sep- arating the space in which the quantum min(TTS) is less than the classical min(TTS) and vice versa. For the small sizes we have considered (low complexity) the quantum Metropolis provides modest results when compared to those obtained by the classical Metropolis. How- ever, the key point of these plots is to note that the fact that the exponent is less than 1 leads us to expect that the quantum advantage dominates as the complexity increases and makes the quantum Metropolis more useful than its classical counterpart. Furthermore, the next logical step to increase the complexity of the simulation is to infer more parameters at the same time. The results for this inference of four parameters is shown in figure 3. In summary, since we cannot run a quantum Metropolis–Hastings algorithm at high levels of complexity, we run the quantum and classical algorithms for simple cases and see how it 11 Class. Quantum Grav. 40 (2023) 045001 G Escrig et al Figure 3. Comparison of the minimum value of the classical and quantumTTS achieved for the 4-parameter inference simulations. On the left the dimensionless spin-magnitude of the larger and smaller object, the comoving volume and the inclination have been inferred. On the right the source-frame chirp mass, the mass ratio, the coalescence time and the comoving volume have been inferred. All results have been obtained for 2–3 bits of precision with a constant value of β. Table 1. Fitted exponents for different sets of variable parameters. 2 PARAMETERS INFERENCE Parameters inferred Fit Exponents Chirp mass and mass ratio 0.95 Dimensionless spin1 and spin2 0.89 Comoving volume and inclination 0.87 Comoving volume and coalesence phase 0.88 4 PARAMETERS INFERENCE Parameters inferred Fit Exponents Coalesence time, masses and comoving volume 0.68 Dimensionless spin1, spin2, comoving volume and inclination 0.67 behaves as the complexity increases. In this way, the fit for the calculation of the exponent is shown in table 1. Although only inferred on a small data set, this result may give us some hints as to whether our technique based on the use of a quantumMetropolis–Hastings algorithm, will be useful for estimating GW source parameters if a sufficiently large and error tolerant quantum computer is available, a point that will be discussed in the conclusions. 5. Conclusions and outlook Throughout this work we have studied how quantum computation could complement classical techniques in the inference of gravitational wave parameters. To this aim, we have explored the feasibility of a quantumMetropolis–Hastings algorithm based on quantumwalks for parameter inference in GW astronomy. 12 Class. Quantum Grav. 40 (2023) 045001 G Escrig et al This algorithm has been built to run both the classical and quantum Metropolis–Hastings algorithm, since it is our aim to compare these two through the TTS figure of merit. For this purpose we have studied the first 11 gravitational wave events of the GWTC-1 catalogue. It is currently not possible to run the Metropolis–Hastings quantum algorithms on real quantum hardware to compare them under the same conditions as their classical versions. How- ever, it is possible to simulate the behavior of the quantumMetropolis–Hastings algorithmwith classical hardware and compare their behaviors using an appropriate figure of merit, although we will effectively be limited in the size of the simulation, due to the number of complex num- bers needed to describe a quantum system growing exponentially with the size of the system. Thus, we have studied the algorithms for simple but magnitude relevant examples in a black hole merger. Specifically, by inferring 2 and 4 event parameters at the same time, increasing the complexity of the simulation by increasing the number of precision bits to test the behavior of the quantum algorithm against the classical one. In the case of inference of two parameters at the same time, it has been possible to run the algorithm up to 6 bits of precision, which implies a discretization comparable to the uncertainty of classical parameter estimation estimates [9, 29]. In the case of 4-parameter inference at the same time, the limit on the number of precision bits is reduced to half of the 2-parameter case above, due to the large computational cost of simulating a quantum processor on classical hardware. The maximum number of bits achieved in these simulations has been 3. It is important to note that this is not a limitation in the imple- mentation of quantum Metropolis–Hastings algorithms to gravitational wave astronomy, it is only due to the large computational cost of simulating quantum algorithms on classical hard- ware, which will not be a problem when run on sufficiently large and fault tolerant quantum hardware. In this way, in the 2-parameter inference a clear polynomial advantage has been obtained, achieving exponents of around 0.8∼ 0.9. This result is very interesting because, although a considerable number of points are already in the quantum advantage regime, this trend indic- ates that the advantage will become clearer the higher the complexity of the simulation. Although we have dealt with problems of reduced sizes (complexity) where a polynomial quantum advantage is not yet useful, however the importance of our results is that they allow us to extrapolate this advantage to situations where the complexity of the simulation is very high, as it is in the inference of a larger set of parameters. The trend we found with this polynomial quantum advantage would imply an improvement of the quantum algorithms with respect to the classical ones, thus being able to reduce quantitatively and sensibly the execution times. In the case of 4-parameter inference, the results indicate a higher polynomial advantage of quantum algorithms over classical algorithms than the 2-parameter case above. This could be due to the fact that we are increasing the complexity of the simulation by inferring more parameters at the same time. Still, it is interesting the fact that the vast majority of the points are already at a quantum advantage. Although it may seem that the number of events with 2 bits in the quantum advantage zone is higher than for 3 bits, which might lead one to believe that the quantum advantage is now smaller as the precision bits increase, the decisive factor that marks the possible quantum advantage while increasing complexity is the scaling exponent of the fitting. Then, we find that the scaling exponent is indeed smaller even than for the 2-parameter inference, finding an even more remarkable quantum advantage. There are dozens of combinations in which to infer 2 parameters at the same time and hundreds in which to infer 4 parameters at the same time. There is also a high computational cost required for both the extraction of the energy files and the execution of the Metropolis– Hastings algorithms. Therefore, wewanted to perform the simulations on parameter sets where 13 Class. Quantum Grav. 40 (2023) 045001 G Escrig et al there is a good physical motivation to infer at the same time. The results obtained are very encouraging about the use of this type of quantum techniques in the study of gravitational radiation. Thus, in this work we have proposed for the first time a quantum algorithm based on quantum walks in gravitational wave astronomy, achieving polynomial advantages over clas- sical algorithms for the inference of a reduced set of physical parameters. This represents a first starting point for possible future developments and providing convincing arguments to motivate the use of quantum computation in the field of astrophysics. Data availability statement The data that support the findings of this study are openly available at the follow- ing URL/DOI: PyCBC, Data used of 4-OGC 2021 https://github.com/gwastro/4-ogc (also available at: www.atlas.aei.uni-hannover.de/work/yifan.wang/4ogc/release_prod/convertsnr/ 1212_posterior/) (Accessed 10 August 2022). Acknowledgment We would like to thank Yifan Wang, Alex Nitz and Collin Capano on the usage of PyCBC. We acknowledge support from the CAM/FEDER Project No. S2018/TCS-4342 (QUITEMAD-CM), Spanish MINECO Grants MINECO/FEDER Projects, PGC2018- 099169-BI00 FIS2018,MCINwith funding fromEuropeanUnionNextGenerationEU (PRTR- C17.I1) and Ministry of Economic Affairs Quantum ENIA project. M A M-D has been par- tially supported by the U.S. Army Research Office through Grant No. W911NF-14-1-0103. P A M C thanks the support of a MECD Grant FPU17/03620, and R C the support of a CAM Grant IND2019/TIC17146. ORCID iDs Gabriel Escrig https://orcid.org/0000-0003-2881-085X Roberto Campos https://orcid.org/0000-0002-2527-4177 Pablo A M Casares https://orcid.org/0000-0001-5500-9115 M A Martin-Delgado https://orcid.org/0000-0003-2746-5062 References [1] Abbott B P et al (LIGO Scientific Collaboration and Virgo Collaboration) 2016 Observation of gravitational waves from a binary black hole merger Phys. Rev. Lett. 116 061102 [2] Christensen N and Meyer R 2022 Parameter estimation with gravitational waves Rev. Mod. Phys. 94 025001 [3] DÁlvares J et al 2021 Exploring gravitational-wave detection and parameter inference using deep learning methods Class. Quantum Grav. 38 155010 [4] Röver C 2007 Bayesian inference on astrophysical binary inspirals based on gravitational-wave measurements PhD Thesis Auckland U [5] Price-Whelan A M et al The Astropy Collaboration 2018 The astropy project: building an open- science project and status of the v2.0 core package Astron. J. 156 123 [6] Usman S A et al 2016 The PyCBC search for gravitational waves from compact binary coalescence Class. Quantum Grav. 33 215004 [7] Ashton G et al 2019 Bilby: a user-friendly Bayesian inference library for gravitational-wave astro- nomy Astrophys. J. Suppl. Ser. 241 27 14 https://github.com/gwastro/4-ogc http://www.atlas.aei.uni-hannover.de/work/yifan.wang/4ogc/release_prod/convertsnr/1212_posterior/ http://www.atlas.aei.uni-hannover.de/work/yifan.wang/4ogc/release_prod/convertsnr/1212_posterior/ https://orcid.org/0000-0003-2881-085X https://orcid.org/0000-0003-2881-085X https://orcid.org/0000-0002-2527-4177 https://orcid.org/0000-0002-2527-4177 https://orcid.org/0000-0001-5500-9115 https://orcid.org/0000-0001-5500-9115 https://orcid.org/0000-0003-2746-5062 https://orcid.org/0000-0003-2746-5062 https://doi.org/10.1103/PhysRevLett.116.061102 https://doi.org/10.1103/PhysRevLett.116.061102 https://doi.org/10.1103/RevModPhys.94.025001 https://doi.org/10.1103/RevModPhys.94.025001 https://doi.org/10.1088/1361-6382/ac0455 https://doi.org/10.1088/1361-6382/ac0455 https://doi.org/10.3847/1538-3881/aabc4f https://doi.org/10.3847/1538-3881/aabc4f https://doi.org/10.1088/0264-9381/33/21/215004 https://doi.org/10.1088/0264-9381/33/21/215004 https://doi.org/10.3847/1538-4365/ab06fc https://doi.org/10.3847/1538-4365/ab06fc Class. Quantum Grav. 40 (2023) 045001 G Escrig et al [8] Abbott B P et al 2020 Prospects for observing and localizing gravitational-wave transients with advanced LIGO, advanced Virgo and KAGRA Living Rev. Relativ. 23 3 [9] Nitz A H, Kumar S, Wang Y-F, Kastha S, Wu S, Schäfer M, Dhurkunde R and Capano C D 2021 4-OGC: catalog of gravitational waves from compact-binary mergers (arXiv:2112.06878) [10] Somma R D, Boixo S and Barnum H 2007 Quantum simulated annealing (arXiv:0712.1008) [11] Somma R D, Boixo S, Barnum H and Knill E 2008 Quantum simulations of classical annealing processes Phys. Rev. Lett. 101 130504 [12] Wocjan P and Abeyesinghe A 2008 Speedup via quantum sampling Phys. Rev. A 78 042336 [13] Yung M-H and Aspuru-Guzik A 2012 A quantum–quantum metropolis algorithm Proc. Natl Acad. Sci. 109 754 [14] Lemieux J, Heim B, Poulin D, Svore K and Troyer M 2020 Efficient quantum walk circuits for Metropolis-Hastings algorithm Quantum 4 287 [15] Lemieux J, Duclos-Cianci G, Sénéchal D and Poulin D 2021 Resource estimate for quantummany- body ground-state preparation on a quantum computer Phys. Rev. A 103 052408 [16] Casares P AM, Campos R andMartin-DelgadoMA 2022 Qfold: quantumwalks and deep learning to solve protein folding Quantum Sci. Technol. 7 025013 [17] Campos R, Casares P A M and Martin-Delgado M A 2022 Quantum metropolis solver: a quantum walks approach to optimization problems (arXiv:2207.06462) [18] Chakraborty S, Novo L, Ambainis A and Omar Y 2016 Spatial search by quantum walk is optimal for almost all graphs Phys. Rev. Lett. 116 100501 [19] Gao S, Hayes F, Croke S,Messenger C andVeitch J 2022Quantum algorithm for gravitational-wave matched filtering Phys. Rev. Res. 4 023006 [20] Veitch J et al 2015 Parameter estimation for compact binaries with ground-based gravitational-wave observations using the lalinference software library Phys. Rev. D 91 042003 [21] van der Sluys M, Raymond V, Mandel I, Röver C, Christensen N, Kalogera V, Meyer R and Vec- chio A 2008 Parameter estimation of spinning binary inspirals using Markov Chain Monte Carlo Class. Quantum Grav. 25 184011 [22] Röver C, Meyer R and Christensen N 2007 Coherent Bayesian inference on compact binary inspir- als using a network of interferometric gravitational wave detectors Phys. Rev. D 75 062004 [23] Thrane E and Talbot C 2019 An introduction to Bayesian inference in gravitational-wave astro- nomy: parameter estimation, model selection and hierarchical models Publ. Astron. Soc. Aust. 36 e010 [24] Szegedy M 2004 Quantum speed-up of Markov Chain based algorithms 45th Annual IEEE Symp. on Foundations of Computer Science pp 32–41 [25] Portugal R 2018 Quantum Walks and Search Algorithms (New York: Springer) (https://doi.org/ 10.1007/978-1-4614-6336-8) [26] Galindo A andMartín-DelgadoMA 2000 Family of Grover’s quantum-searching algorithms Phys. Rev. A 62 062303 [27] Galindo A and Martín-Delgado M A 2002 Information and computation: classical and quantum aspects Rev. Mod. Phys. 74 347 [28] ANIS M S et al 2021 Qiskit: an open-source framework for quantum computing (available at: https://github.com/Qiskit/qiskit/tree/0.36.2) [29] Abbott B P et al (LIGO Scientific Collaboration and Virgo Collaboration) 2019 GWTC-1: a gravitational-wave transient catalog of compact binary mergers observed by LIGO and Virgo during the first and second observing runs Phys. Rev. X 9 031040 [30] Babbush R, McClean J R, Newman M, Gidney C, Boixo S and Neven H 2021 Focus beyond quad- ratic speedups for error-corrected quantum advantage PRX Quantum 2 010103 15 https://doi.org/10.1007/s41114-020-00026-9 https://doi.org/10.1007/s41114-020-00026-9 https://arxiv.org/abs/2112.06878 https://arxiv.org/abs/0712.1008 https://doi.org/10.1103/PhysRevLett.101.130504 https://doi.org/10.1103/PhysRevLett.101.130504 https://doi.org/10.1103/PhysRevA.78.042336 https://doi.org/10.1103/PhysRevA.78.042336 https://doi.org/10.1073/pnas.1111758109 https://doi.org/10.1073/pnas.1111758109 https://doi.org/10.22331/q-2020-06-29-287 https://doi.org/10.22331/q-2020-06-29-287 https://doi.org/10.1103/PhysRevA.103.052408 https://doi.org/10.1103/PhysRevA.103.052408 https://doi.org/10.1088/2058-9565/ac4f2f https://doi.org/10.1088/2058-9565/ac4f2f https://arxiv.org/abs/2207.06462 https://doi.org/10.1103/PhysRevLett.116.100501 https://doi.org/10.1103/PhysRevLett.116.100501 https://doi.org/10.1103/PhysRevResearch.4.023006 https://doi.org/10.1103/PhysRevResearch.4.023006 https://doi.org/10.1103/PhysRevD.91.042003 https://doi.org/10.1103/PhysRevD.91.042003 https://doi.org/10.1088/0264-9381/25/18/184011 https://doi.org/10.1088/0264-9381/25/18/184011 https://doi.org/10.1103/PhysRevD.75.062004 https://doi.org/10.1103/PhysRevD.75.062004 https://doi.org/10.1017/pasa.2019.2 https://doi.org/10.1017/pasa.2019.2 https://doi.org/10.1007/978-1-4614-6336-8 https://doi.org/10.1007/978-1-4614-6336-8 https://doi.org/10.1103/PhysRevA.62.062303 https://doi.org/10.1103/PhysRevA.62.062303 https://doi.org/10.1103/RevModPhys.74.347 https://doi.org/10.1103/RevModPhys.74.347 https://github.com/Qiskit/qiskit/tree/0.36.2 https://doi.org/10.1103/PhysRevX.9.031040 https://doi.org/10.1103/PhysRevX.9.031040 https://doi.org/10.1103/PRXQuantum.2.010103 https://doi.org/10.1103/PRXQuantum.2.010103 Parameter estimation of gravitational waves with a quantum metropolis algorithm 1. Introduction 2. Bayesian inference: classical and quantum Metropolis–Hastings algorithms 2.1. Classical Metropolis–Hastings 2.2. Quantum Metropolis–Hastings 2.2.1. Construction of the quantum Metropolis–Hastings algorithm. 2.3. Figure of merit: the total time to solution 3. Simulation setup 3.1. The energy calculator module 3.2. The Metropolis–Hastings solver module 3.3. The output processor module 4. Simulation results 5. Conclusions and outlook References