Toward Prediction of Financial Crashes with a D-Wave Quantum Annealer

The prediction of financial crashes in a complex financial network is known to be an NP-hard problem, which means that no known algorithm can efficiently find optimal solutions. We experimentally explore a novel approach to this problem by using a D-Wave quantum annealer, benchmarking its performance for attaining a financial equilibrium. To be specific, the equilibrium condition of a nonlinear financial model is embedded into a higher-order unconstrained binary optimization (HUBO) problem, which is then transformed into a spin-1/2 Hamiltonian with at most, two-qubit interactions. The problem is thus equivalent to finding the ground state of an interacting spin Hamiltonian, which can be approximated with a quantum annealer. The size of the simulation is mainly constrained by the necessity of a large number of physical qubits representing a logical qubit with the correct connectivity. Our experiment paves the way for the codification of this quantitative macroeconomics problem in quantum annealers.

Prediction of financial crashes in a complex financial network is known to be an NP-hard problem, which means that no known algorithm can guarantee to find optimal solutions efficiently. We experimentally explore a novel approach to this problem by using a D-Wave quantum annealer, benchmarking its performance for attaining a financial equilibrium. To be specific, the equilibrium condition of a nonlinear financial model is embedded into a higher-order unconstrained binary optimization (HUBO) problem, which is then transformed to a spin-1/2 Hamiltonian with at most two-qubit interactions. The problem is thus equivalent to finding the ground state of an interacting spin Hamiltonian, which can be approximated with a quantum annealer. The size of the simulation is mainly constrained by the necessity of a large quantity of physical qubits representing a logical qubit with the correct connectivity. Our experiment paves the way to codify this quantitative macroeconomics problem in quantum annealers.

I. INTRODUCTION
Economics is a complex science in which the agents' psychology plays an essential role that is often hardly grasped by mathematical models. However, economists do not relinquish to try to predict market behavior employing sophisticated models, leading to the field of quantitative finance. Following this idea, quantitative finance and economics emerged. They were applied to understand the evolution of financial markets and economies, as well as to forecast their possible future. A realistic question in risk management is: would there be a drastic drop in the market values if the prices of assets suffer some small perturbation? The cross-holdings and nonlinear character of financial network dynamics will cause chain reactions, implying that sudden drops of a market value might affect other nodes in the network resulting in a financial crisis. Presently, the prediction of crashes is mainly performed by studying previous cases in history and comparing with the current configuration [1][2][3][4][5][6]. While this empirical approach has been succesful [7], the economic environment is constantly evolving. Hence, we cannot limit ourselves to predicting economic disasters which are qualitatively similar to past events. Therefore, ab initio simulations of financial networks will become essential for avoiding financial crises. This problem was recently shown to be NP-Hard [8]. Therefore, given the current standpoint on complexity theory, this problem is not expected to be efficiently solvable in a classical computer. Indeed, given the global knowledge of a financial network, the time to compute the consequences of a perturbation would by far exceed the age of the universe.
An alternative approach to this problem is presented in Refs. [9,10], where they depict how to tackle this type of problems by using quantum annealers. In particular, a mathematically identical problem is simulated and the corresponding result measured [11][12][13][14]. Specifically, it was shown that obtaining the equilibrium configuration of a financial network is equivalent to solving a higher-order unconstrained binary optimization (HUBO) problem, which should be feasible for a quantum annealer allowing for multi-qubit interactions. Unfortunately, this hardware has not been realized yet, as stateof-the-art quantum annealers are restricted to two-qubit interactions [15]. A possible work-around, which comes at the price of introducing ancillary qubits, is to find an effective Hamiltonian with the same low-energy sub-arXiv:1904.05808v4 [quant-ph] 16 Feb 2023 space and two-qubit interactions at most. This leaves us with the problem of solving a quadratic unconstrained binary optimization (QUBO) problem whose optimum encodes the equilibrium configuration of a financial network. This problem can be addressed employing a quantum annealer. The D-Wave 2000Q quantum annealer, equipped with a Chimera architecture, requires a large quantity of physical qubits to obtain the desired connectivity and limits the number of institutions and assets considered. An analysis of the changes experienced by the financial network to reach its equilibrium configuration will tell whether a crash has occurred.
In this paper, we experimentally implement the study presented in Refs. [9,10]. Specifically, we compute the equilibrium configuration of a financial network before and after perturbation with a D-Wave 2000Q quantum annealer, and compare the result to alternative methods. Although the D-Wave machine has been successfully used to solve problems in engineering [16], cryptography [17], biology [18], and quantitative finance [19,20] among others, it is the first time that quantum annealing is applied to solve a problem of macroeconomics. This should attract more attention from the finance and economic disciplines towards quantum computing [21][22][23][24][25][26][27].
The contents are organized as follows: in Sec. II, we introduce the model of financial network that will be considered. Sec. III reviews the quantum annealing algorithm to find financial equilibrium. Sec. IV experimentally proves the validity of the scheme by finding the financial equilibrium of a random network of the largest implementable size with a D-Wave 2000Q quantum annealer; for this network, we also show experimentally how the scheme allows to compute the financial equilibrium. Sec. V analyzes the achieved results and discusses further possible improvements. The conclusions drawn from the work are shown in Sec. VI.

II. FORMULATION OF THE MODEL
A nonlinear network model for financial markets is proposed in Ref. [9]. It is made up of n institutions and FIG. 1. Example of a financial network: the yellow nodes and green nodes denote institutions and assets, respectively. Links denote ownerships and cross-holdings, described by the ownership matrices D and C, respectively. Diagonal matrix C represents the self-ownership of institutions, which would be plotted as self-loops in the graph representation. The equity value Vi of institution i is defined by summing up its ownership of all assets and cross-holdings. m assets, and aims at representing the market values of institutions by mapping it onto a graph, as shown in Fig. 1. We codify the prices of the m assets by an m−dimensional vector p ∈ R m , where the element p k represents the price of asset k. Moreover, an n × m ownership matrix D can be defined such that the element D ik ≥ 0 corresponds to the percentage of asset k owned by institution i. There is also an n × n ownership matrix C that describes the cross-holdings and self-ownerships between institutions. The coefficients C ij denote the percentage of institution j owned by institution i. By considering all self-ownerships (i.e., the diagonal elements) from C one forms a new diagonal matrixC which represents the self-ownership only, such that matrix C = C −C codifies all cross-holdings. The equity value V i of institution i is defined by summing up its ownership of all assets and cross-holdings, V i = Σ k D ik p k + Σ j C ij V j . One thus obtains a matrix equation V = D p + C V , where equity value vector V ∈ R n is an n−dimensional vector. Accordingly, the market value is the equity value rescaled with its self-ownership, resulting in the n−dimensional market value vector v =C V . The solution to the linear matrix equation thus reads We introduce the nonlinear effect of panic in the model via a Heaviside-theta function Θ; if the market value v i drops below the critical value v i c , failure of institution i occurs and its equity value drops by β i ( p) which is governed by the price vector of assets. Once we define the failure vector b( v, p) = β( p) • (1 − Θ( v − v c )), where • denotes the Hadamard product, the market value vector with nonlinearity can be written as Mathematically, it is the nonlinearity of b( v, p) which makes financial networks so hard to be predicted. This drop may cause an institution's value to crash, a behavior which can infect other nodes in the network. Under our definition, a financial crash happens when, after a perturbation in the assets price, the market value of an institution considering the nonlinear term is lower than those pre-perturbation prices calculated with the linear model.

III. QUANTUM ANNEALING ALGORITHM
As proposed in Ref. [9] finding financial equilibrium can be represented as the minimization of an objective function, which is equivalent to finding the ground state of a classical spin Hamiltonian.
By squaring Eq. (2), we obtain an objective function that meets its minimum value when the market value state is set to be the equilibrium state Thus, our task is now to find the v that minimizes Obj( v) for a given financial network, which is an NP-Hard problem [28]. Next, we need to deal with the nonlinear terms (modeling failure) of the objective function. The reason is that once the objective function is transformed to a spin-1/2 Hamiltonian, it should ideally be made of polynomial terms only, due to the limitations of quantum annealers. Thus, one expands the failure terms with Heaviside-theta functions in terms of polynomials. This expansion is not unique, and here we choose the Legendre expansion [9], in the domain [−1, 1], with P l (x) to be the l-th Legendre polynomial. By setting Using this expansion as an example, we take the approximation The polynomial expansion removes the discontinuity while maintaining the strong nonlinearity of the network. We now encode the continuous variables v i with classical bits. This will allow rewriting the resulting objective function into digital form. The expansion is straightforward, and reads v i = ∞ α=−∞ x i,α 2 α . However, due to the limited resources in real-world devices, one must truncate this expansion, i.e., v i ≈ q α=−q x i,α 2 α , where x i,α are classical bits with binary values 0 or 1. In this way, the market value of institution i is encoded with 2q + 1 classical bits. The maximal market value v max i is given by 0≤α≤p 2 αmα x i,α , the resulting objective function is a polynomial of binary variables x i,α of degree 2r.
ij . To express it as a spin-1/2 Hamiltonian, we replace the binary variables x i,α by qubit operatorsx i,α with eigenvalues 0 and 1, i.e.,x i,α |0 = 0,x i,α |1 = |1 . The Pauli-z operator satisfiesx i,α = (1 +σ z i,α )/2, and therefore the Hamiltonian encodes the objective function but written with Pauli matrices, including all types of multispin interactions, up to 2r-body terms. The Hamiltonian obtained is appropriate for a quantum annealer that allows many-qubit interactions. However, state-of-the-art quantum annealers only accept inputs with at most twoqubit interactions. Finding the ground state of a spin-1/2 Hamiltonian is equivalent to solving a Quadratic Unconstrained Binary Optimization (QUBO) problem, which is the input of the quantum annealer. Thus, we should recast our quantum Hamiltonian into a modified, effective Hamiltonian with at most two-qubit interactions. Some protocols achieving exactly this are proposed in Refs. [29][30][31][32][33][34][35], in particular we base our protocol in Ref. [35], where k ancilla qubits are introduced to implement an effective k-qubit interaction. Suppose that there is a k-qubit interaction termĤ k = J k Π k i=1 σ z i with the same low-energy spectrum of another Hamiltonian termH k with at most two-qubit interactions. We can expressH k with k logical qubits and k extra ancilla qubits as as represented in Fig. 2. This two-qubit Hamiltonian has the same low-energy spectrum thanĤ k when J, J a , h and h a i are set to some appropriate values. As Ref. [35] suggested, this can be achieved once These conditions can be relaxed to |J k | < q 0 < J a and |J k | < J a −q 0 < J a if one aims at having the same ground state only, instead of the whole low-energy sector. We depict the low-energy spectrum of this twoqubit Hamiltonian for k = 3 logical qubits in Table I.

IV. IMPLEMENTATION IN A D-WAVE 2000Q QUANTUM ANNEALER
Once shown that it is possible to recast the problem of finding financial equilibrium into a language amenable to QUBO solvers, and in particular, quantum annealers, this section deals with its implementation using a state-of-the-art quantum annealer, namely, the D-Wave 2000Q. This quantum annealer consists of up to 2048 qubits connected according to the Chimera graph topology (see Fig. 3). It is designed to solve embedded Ising problems or QUBO problems.
Two simulations were produced: 1. A financial network without failure term, which is simple to solve on a classical computer in order to benchmark the performance of the quantum processor.

2.
A financial network with the inherently nonlinear risk of failure. We will perturb the asset price vector in this network to compute the new equilibrium configuration using the quantum annealing algorithm.
We initially generate a financial network with 10 institutions and 15 assets. To demonstrate the algorithm, we Chimera graph topology implemented by the D-Wave 2000Q quantum annealer. The 2048 qubits are partitioned into subgraphs of 8 qubits. The connection between subgraphs is sparse, in each of these subgraphs there are two sets of four qubits; each qubit connects to all qubits in the other set but to none in its own, forming a K4,4 bipartite graph.
randomize the ownership matrix D with a Dirichlet distribution that satisfies n i=1 D ij = 1, where D ij are random variables. The cross-holding matrix C is generated in a similar way but with a constraint that all diagonal elements should be larger than 0.5, ensuring that all institutions can make decisions according to their own wills. Thus, we randomizeC ii between 0.5 and 1 and randomize n i=1 C ij = 1−C jj with a rescaled Dirichlet distribution. The price vector p is also random, with p i ∈ [10,40]. The network configuration is shown in Figs. 4 a) and 4 b).
We can calculate the equilibrium state v q and the equity value vector V on a classical computer using which are linear equations that in fact can be implemented in a quantum annealer by using only 2-local terms, result of squaring the expression in a similar way to Eq. (3). The objective function shown in Eq. (3) was implemented, for benchmarking reasons, both in a quantum annealer and a classical simulator. Variables v i were encoded, v i = 6 α=0 2 α x i,α , on seven qubits. As such, this constrains the v i to be integers smaller than 127. A quantum implementation of this algorithm does not require ancilla qubits, as there are no many-qubit interactions.
The QUBO for this linear problem is a 70 × 70 matrix, with 210 couplers which cannot be solved directly due to the topology structure of the quantum annealer. D-Wave provides a software named qbsolv that allows to combine the classical computer with its quantum annealer by splitting the QUBO matrix into partition matrices that can be embedded in the quantum annealer. As a decomposing solver, it finds a minimum value of a large QUBO problem by splitting it into pieces solved either via a D-Wave system or a classical tabu solver (both approaches were considered here for comparison purposes). Since the D-Wave 2000Q processor is a quantum annealer, 20 results would be obtained from a qbsolv process with a default setting; these results should be handled by a correction process, e.g., majority voting, to help us identify the most plausible answer. The result of this QUBO problem is shown in Fig. 5, where the exact solution via solving a linear matrix equation, qbsolv solution with classical tabu solver, and qbsolv solution with D-Wave quantum annealer, are compared. It is straightforward to observe that comparing the individual equilibrium values, a quantum annealer provides a solution that present more accurate individual values of the assets than the prediction from the classical solver.
While the failure-free model only has linear and quadratic terms in v i , the nonlinear model has powers of v i up to order 2r. For large r, this can be extremely resource-consuming in terms of ancillary qubits due to the requested connectivity. An estimation of the number of qubits can be made by counting the number of interaction terms; Our HamiltonianĤ can have up to We randomize the ownership matrix D with a Dirichlet distribution that satisfies n i=1 Dij = 1. b) Cross-holding matrix C for the linear model that describes the cross-holdings and self-ownerships between institutions. Cross-holding matrix is generated in a similar way to ownership matrix but with a constraint that all diagonal elements should be larger than 0.5, ensuring that all institutions can make decisions according to their own wills. All of these data, as well of the asset prices, have been synthetically produced but following all constraint conditions proposed in the theoretical model [9]. 2r α=0 n(2q+1) α terms, where n(2q + 1) denotes the logical qubits that are required. In each term, 3-to-2r new ancilla qubits are needed, depending on the number of logical qubits in this term. Therefore, the number of necessary qubits grows rapidly with the degree of the polynomial expansion r. Note that the aforementioned QUBO problem is NP-hard for any n ≥ 2. In practice, this is an upper bound to the required resources, calculated assuming thatĤ has all possible terms up to order O(2r).
Here, we implement an enhanced model with failure terms on the basis of the linear model previously simulated. We perturb the vector of asset prices, leaving the ownership matrix D and cross-holding matrix C invariant, and recompute the equilibrium state. Specifically,

Market Values v_i
Equilibrium State we set the price of some random assets to zero (to simulate, e.g: the assets' destruction). In this study, we will use an expansion ofĤ to third order, which still characterizes the phenomenon of sudden drop near the critical value. Moreover, this approach provides strong nonlinearity while saving plenty of qubit resources. As a result, 70 logical qubits and 872,690 ancilla qubits are required, which leads to a QUBO matrix of 872,760x872,760 entries, although only a minority of them, 4,446,575 couplers, are non-zero. Storing this sparse matrix results in the requirement of about 6TB RAM memory, since each element has an accuracy of double float in qbsolv. Due to the limitations of state-of-the-art techniques, the network is reduced to three institutions and each market value v i is encoded by five qubits, bounding the maximum market value to be 31. New 3 × 7 ownership matrix D and 3 × 3 cross-holding matrix C are generated while the price vector p before perturbation is p = {8.43, 14.47, 6.75, 8.09, 19.11 , 11.32, 7.19} T . The network configuration is shown in Figs. 6 a) and 6 b). The equilibrium state before perturbation without nonlinearity is given as v q = {21.18 23.33, 30.83} T , and the critical value vector is still set to be 80% of the original equilibrium state, while the failure strength β is considered to be 30% of the original equity value. The corresponding perturbed price vector is given as p = {8.43, 14.47, 0, 8.09, 0 , 11.32, 7.19} T . Before calculating the new equilibrum state with nonlinearity and perturbation, some parameters, like J a and q 0 , must be set. For the minor embedding of a submatrix in the D-Wave quantum annealer, this is done by introducing a penalty function between qubits in the Chimera graph requiring J m ≥ J a , which means that the J a for mapping multi-qubit interactions to two-qubit interactions should be in the proper scale. Meanwhile, as we mentioned in the theory part, we need to sample out the thermal fluctuation by assuming that |Ĥ k | is much smaller than J a , or the protocol will break down because those ancilla qubits will not be in the corresponding ground state anymore. Thus, in the implementation we took J a = 20J k and q 0 = 10J k , such that this could ensure that either q 0 or J a − q 0 would be at least 10 times larger than J k .
For this problem, the QUBO matrix had the size of 8280 × 8280, with 15 logical qubits, 8265 ancilla qubits and 38,790 couplers. Remark that the available quantum annealer structure is not optimized for this problem and, also, that the translation is not efficient because of sparse connectivity of the quantum processor. Finally, we compare our results from the quantum annealer with the integer equilibrium solution calculated with a straightforward method by trying 32 3 times in Fig. 7, which shows a good agreement and the accuracy of the proposed method. Comparing the results after the perturbation with the pre-perturbation values, we can conclude that we have detected the financial crash.

V. RESULTS AND DISCUSSION
D-Wave is a quantum annealer designed to deal with Ising Model and QUBO problems. However, the problem faced in this paper, namely, financial crisis prediction with nonlinearity associated to panic, is not QUBO but HUBO instead, thus requiring multi-qubit interactions. In order to approximate this HUBO problem with twoqubit interactions, at the current stage of hardware and software we were limited to simulate a small financial network, made up of three institutions and cross-holdings.
An effective two-qubit quantum Hamiltonian could still not be read directly in D-Wave system which requires QUBO type input or Ising type input. Although some open-source software like pyqubo can generate it, the input size must be very small in order to avoid a stack overflow associated with recursion errors. A possible solution is to produce a Mathematica script that reads each term, write it as a string of coefficients and qubits in an input file for the D-Wave system. Once we generate the input for this problem, this is still too large to be embedded in the D-Wave 2000Q quantum annealer because of the graph structure. Thus, qbsolv is an inevitable option for us, which works by separating the large matrix to submatrices and solve them by a classical tabu solver or D-Wave solver. This kind of hybrid computation provides the possibility to solve the complicated problem but brings some new constraints, namely: (i) Local hardware. Once the QUBO matrix is provided, qbsolv allocates dynamic memory before separating it to submatrices with elements of double precision floats, by requiring a size of 8n × n bytes of memory. However, the bottleneck is not the memory size but the performance of CPU since a large QUBO matrix will consume exhaustive CPU time if one needs high accuracy of the optimized result; (ii) Algorithm. Instead of a real quantum annealing process for the whole matrix, qbsolv provides a tabu algorithm or D-Wave 2000Q quantum annealer for submatrices. The partition strategy for generating submatrices may get stuck in a local minimum instead of the global minimum that quantum annealing guarantees with high probability under ideal conditions, i.e. in absence of decoherence and in the adiabatic limit. Considering that the logical qubits only encode less than 1% in the QUBO matrix, the risk of getting stuck is still high, even if we sample over the thermal distribution or give a huge repeat limitation in the main loop to improve its accuracy. We would have to customize a random seed for the separation, and check the final result manually, to see whether the result is near from the equilibrium. Another option is that one may send the QUBO matrix to the solver many times and average the result to obtain the best solution. (iii) Quantum annealer. The submatrices will be sent to D-Wave 2000Q quantum annealing device for optimization after they are generated by Glover's algorithm [37]. In the quantum annealing process, magnetic fields are applied to the processors and the strength should be accurate because J k , J a in the QUBO matrix and J m for the embedding belong to different magnitudes. Any imprecision in the system preparation will cause significant deviation from the correct result.
In this implementation, the accuracy is not especially high, since we are not optimizing the objective function rigourously because the market values are integers v i ∈ [0, v max ] constrained for the qubits we take to encode them. The computation time is also long, considering that there is a straightforward but equivalent classical algorithm by testing the value of the objective function 32 3 times by brute force, corresponding to all the possible combinations. Although mapping it to a QUBO problem and optimizing it with a general quantum annealer is not efficient enough for current technology, we believe it is a valuable example of how one can solve an NP-hard problem via quantum computation. With quantum annealers designed for solving HUBO problems that allow the implementation of multi-qubit interactions, we would avoid the overhead of resources and we may obtain a speed up factor in forecasting the behavior of complex financial networks over the use of general purpose annealers. We expect this kind of quantum solver may be available in the near future. Meanwhile, D-Wave has recently announce its next generation of quantum annealers called the Advantage system [38]. It would consist of more than 5000 qubits connected with each other according to the Pegasus topology. In this manner, one could improve the number of qubits and the connectivity by a factor of 2.5.
Considering that a specialized quantum annealer for HUBO problems would not be available to the public anytime soon, we now analyze the possible ways to enhance the performance of D-Wave 2000Q quantum annealer on this problem. After compromising on the maximum two-qubit interactions in hardware, the subsequent The element D ik ≥ 0 corresponds to the percentage of asset k owned by institution i. We randomize the ownership matrix D with a Dirichlet distribution that satisfies n i=1 Dij = 1. b) Cross-holding matrix C for the implemented network with failure terms that describes the cross-holdings and self-ownerships between institutions. Cross-holding matrix is generated in a similar way to ownership matrix but with a constraint that all diagonal elements should be larger than 0.5, ensuring that all institutions can make decisions according to their own wills.
strategy will be reducing the number of ancilla qubits. With fewer ancilla qubits, the size and accuracy of a solvable network can be improved. As proposed in Ref. [35], is the equilibrium state without taking non-linearity terms (perturbations) into consideration, where the asset price is calculated by inverting the matrix of Eq. (1). The second group (center) is the equilibrium state after taking nonlinearity as the 'failure term', which is activated by a critical value vector 80% of the original equilibrium state calculated with a straightforward method by trying 32 3 times by brute force, corresponding to all the possible combinations. The third group (right) shows the outcome from qbsolv software in D-Wave 2000Q. The error bar characterizes a 95% confidence interval. The agreement between integer and annealer solutions confirms the feasibility and accuracy of the method. Additionally, by comparing both with the pre-perturbation values, we can conclude that we have detected the financial crash.
the multi-to-two mapping is a general method, but for three-to-two, for example, a more efficient mapping can be constructed with only one ancilla qubit. Suppose there is a sub-Hamiltonian of three-qubit interactionŝ A subgraph with full connectivity of three logical qubits and one ancilla qubit is shown in Fig. 8, where the equivalent Hamiltonian is given as (11) At variance with the previous protocol, J a = 2J > h and h a = 2h = 2J 3 . Also, for sampling out the thermal fluctuation, we take J a ≥ J 3 , to prevent the protocol to fail for the same reason. The ancilla qubits can be reduced to about 7000 with this method. Meanwhile, the partition method in qbsolv may cause the system to get stuck in local minima which requires a better algorithm in the main loop.

VI. CONCLUSION
We have implemented in a D-Wave quantum annealer the algorithm proposed in Ref. [9], to solve the equilibrium state of a complex financial network that predicts financial crashes. Although the size of the studied financial network is limited, this proof of principle is in agreement with the result of an exhaustive search. This result may be improved with the design of a customized 8. An efficient encoding of three qubits, making use of only one ancilla qubit. The multi-to-two interaction Hamiltonian mapping is a general method, but for three-to-two, a more efficient mapping can be constructed via a subgraph with full connectivity of three logical qubits and one ancilla.
"financial quantum annealer": a quantum processor with suitable connectivity for efficient embedding of this kind of problems. Such coherent quantum annealers can be built with current technology [39][40][41], providing convenient multi-qubit couplings.