Parameter Estimation of Gravitational Waves with a Quantum Metropolis Algorithm

After the first detection of a gravitational wave in 2015, the number of successes 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 quantization of the classical algorithms used in the literature for the inference of gravitational 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 a metric 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.


I. INTRODUCTION
After 100 years of its prediction and an international effort spanning several decades, on September 14th, 2015, the two Advanced 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, complementary 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, different types of algorithms have been proposed (see Refs [2,3]), but those based on the formalism of Bayesian inference have stood out from the rest, making Bayesian parameter estimation the language of gravitational-wave 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 Refs [5][6][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 reasons. 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 * gescrig@ucm.es † robecamp@ucm.es ‡ pabloamo@ucm.es § mardel@ucm.es 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 analysis 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 particular, 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,10].
Quantum walks are a cornerstone family of quantum algorithms that exhibit a quadratically larger eigenvalue gap than the associated random walk. Since the inverse of the eigenvalue gap is often a multiplicative factor in the complexity of the algorithm, quantum walks are capable to outperform their classical counterparts, as we will explain. Specifically, quantum Metropolis-Hastings algorithms exhibiting this quadratic advantage have been developed [11][12][13][14][15], and have also been proposed as a method to tackle hard optimization problems [16][17][18]. Although the use of quantum walks is quite extensive throughout the literature [19], to the best of our knowledge there are no previous analysis of the use of this type of methods in arXiv:2208.05506v2 [gr-qc] 10 May 2023 the study of gravitational wave events, only in their detection [20]. Thus, we hope this work will be of interest to the gravitational astronomical community.
The rest of the manuscript is organized as follows. In the next section II we introduce the basics of classical and quantum Metropolis-Hastings algorithms, and how they compare with each other. Then, in section III we describe our simulation setup, whose results are explained in IV. We end the article with the conclusions and future work in section V.

II. BAYESIAN INFERENCE: CLASSICAL AND QUANTUM METROPOLIS-HASTINGS ALGORITHMS
The standard formalism for tackling the task of estimating the unknown parameters of a gravitational wave source classically is the formalism of Bayesian inference based on Markov Chain Monte Carlo methods (MCMC) [21], 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.

A. 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' Theorem will be our main 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) [22], The models we will explore have many parameters, in particular the models used for the BBH and BNS events contain 15 parameters [10], 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 distribution on the multidimensional space p(θ|d) will be used, being able to recover the probability for a specific parameter by marginalising 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 normalisation factor [23], so that the posterior probability reduces to 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 measurement. 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 [24], so that the Gaussian-noise likelihood function has the form: All these functions are computed in the Fourier domain, with f denoting frequencies, whered represents the data set obtained in the usual way at the detector,μ the waveform for a certain relativistic numerical model calculated from the parameter set θ, and σ the detector noise (spectral density). It is interesting to note that despite having an integral at all possible frequencies, in practice we will limit this interval to those frequencies at which groundbased 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, 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 uniform, distribution φ (0) and a random transition will be proposed to obtain the next state. In this way: 1. 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 .
2. We now calculate the function likelihood found in (3) for the new state, and accept this transition with probability 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 tempering [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 optimisation 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 probability of such a proposal. If we have a uniform prior, it is customary to choose a uniform probability distribution, for N the number of possible transitions from state i. Alternatively, we could choose T ij ∝ p(θj ) p(θi) . Then, this change is accepted with a probability resulting in an overall transition probability for a given step of i → j: As such, we have explained how the Metropolis-Hastings algorithm can be used to solve our Bayesian inference problem.

B. Quantum Metropolis-Hastings
Similarly, the objective for the Metropolis-Hastings algorithm is sampling the Boltzmann distribution ρ β = |φ β φ β | = 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 normalisation 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.

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 [25,26]. 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 Hilbert space. Given the acceptance probabilities W ij = T ij A ij , defined previously in (8), for the transition from state i to state j, one defines the unitary operator 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 [27,28]. Taking 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 The first R 0 represents the first rotation, while The Metropolis-Hastings algorithm we are going to use is a modification of Szegedy's quantum walk [15], 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, U W . This operator will 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 A ij . We will denote them respectively by |x S |m M |b C . The operator U W is equivalent to half of the Szegedy operator W 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 [15]: with Y : |x S |y S → |x S |m M such that m(x) = y. (13) The coined version operator U W is then constructed as 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 A ij 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 M C 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 prepares the possible transitions in the register |· M , these will be given by T ij . The operator V consequently prepares a uniform superposition Then B prepares the coin qubit so that the probability of being in state |1 C corresponds to the probability of acceptance, A ij F will update the system registry to the new configuration, And finally R is the reflection operator on the state |0 M C There have been several proposals to quantise the Metropolis-Hastings algorithm using these quantum walks. 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 [11][12][13][14]. However, classically the Metropolis-Hastings algorithm is usually used with a rapidly varying annealing schedule, where no complexity-theoretical guarantees can be 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 Ref. [15]. It consists of executing L times the evolution operator on an initial state U 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 quantum Metropolis-Hastings algorithm, we will next explore how to compare it with its classical counterpart.

C. Figure of Merit: the Total Time to Solution
Given the two algorithms, the classical algorithm and its quantum version, we are now interested 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, computing 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 applying the W walk sequentially t iterations, the probability of success is p(t) = (W t q)(θ * ). 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 Total Time to Solution (TTS) [15], 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: where t ∈ N is the number of quantum/classical walk steps performed in one execution of the quan- tum/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 Total Time to Solution over the total of all iterations as min(TTS) = min t TTS(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. Starting with a |φ (0) state, it consists of applying the quantum walk operator U W j sequentially, obtaining the state (19). Thus, the probability of obtaining the searched state is: So the construction of the TTS is straightforward, obtaining for the L-th iteration of the algorithm the ex- Similarly, for the classical algorithm, we estimate the probability of obtaining θ * after t steps of the Metropolis-Hastings algorithm.

III. 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 [16] 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 [17]. 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 [30]. The data and the posterior distributions inferred by [10] can be found in [9]. 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 10 orders of magnitude longer than on a classical computer [31]. The quantum advantage does not happen in this 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 β. 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 Total Time to Solution introduced in the previous section.
The pipeline of our library is explained in figure 1, which we will explain in the following sections.

A. The Energy Calculator Module
The algorithm starts once we have chosen the set of parameters we want to infer from the configuration 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 Py-CBC 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 [10]. Thus, this leaves a clear scheme for calculating the energy for any set of parameters for any BBH event for all observation periods.

B. The Metropolis-Hastings Solver Module
The Metropolis-Hastings Solver Module makes use of the software tool QMS [18]. It receives a description of the problem as a list of tuples (state, energy) and it calculates the difference 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 [29]. Delta energies are introduced in the algorithm as a precalculated oracle, or table, queried via a serial QRAM.

C. 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 demonstrate 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.

IV. 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 [30], so taking into account that we use the PyCBC model that marginalises the polarisation in order to reproduce the results of the catalogue of reference [10], 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 twoparameter sets we may highlight the two masses, reparameterised to chirp mass and mass ratio. Also the dimensionless spin-magnitude of the two black holes. In addition, 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 interesting 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 number 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 way the exponent. Starting with the inference of two parameters at a time, this is shown in Fig. 2.
The classical and quantum min(TTS) are on the xaxis 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 separating 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. However, 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 4 parameters is shown in Fig. 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 behaves as the complexity increases. In this way, the fit for the calculation of the exponent is shown in Table I.
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 quantum Metropolis-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.

V. 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 Quantum Metropolis-Hastings Algorithm based on quantum walks for parameter inference in GW astronomy.
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. However, it is possible to simulate the behavior of the quantum Metropolis-Hastings algorithm with 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 numbers 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 [10,30].
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 implementation of quantum Metropolis-Hastings algorithms to gravitational wave astronomy, it is only due to the large computational cost of simulating quantum algorithms on classical hardware, 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 indicates 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, we wanted to perform the simulations on parameter sets where 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 classical 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.