key: cord-0161552-9h5zu7wu
authors: Vesel'y, Martin
title: Application of Quantum Computers in Foreign Exchange Reserves Management
date: 2022-03-29
journal: nan
DOI: nan
sha: 2a990d93e30ec9ee3a08dfb0bb2721c5af0b4686
doc_id: 161552
cord_uid: 9h5zu7wu

The main purpose of this article is to evaluate possible applications of quantum computers in foreign exchange reserves management. The capabilities of quantum computers are demonstrated by means of risk measurement using the quantum Monte Carlo method and portfolio optimization using a linear equations system solver (the Harrow-Hassidim-Lloyd algorithm) and quadratic unconstrained binary optimization (the quantum approximate optimization algorithm). All demonstrations are carried out on the cloud-based IBM Quantum(TM) platform. Despite the fact that real-world applications are impossible under the current state of development of quantum computers, it is proven that in principle it will be possible to apply such computers in FX reserves management in the future. In addition, the article serves as an introduction to quantum computing for the staff of central banks and financial market supervisory authorities.

Computational techniques and instruments have always been an integral part of finance and trading. Even Fibonacci's famous Liber Abaci, besides introducing Arabian (Indian) numerals to medieval Europe, served mainly as a mathematical textbook for merchants and bankers. At the dawn of the modern computer era in the 1950s, corporations quickly recognized the opportunities these machines offered and started to use them not only to speed up their daily operations, but also to build sophisticated models of markets, solve complex optimization problems and so on. This led to the development of completely new branches of science and sectors of the economy. In the last decade, quantum computers have become available to the mass market, and we are now on the cusp of a revolution similar to that initiated by classical computers in the middle of the 20 th century. For this reason, we investigate possible applications of quantum computers in central banks, in particular their use for foreign exchange reserves management. This is the main aim of this paper. A secondary goal of the paper is to serve as an introduction to quantum computing for the staff of central banks and other financial market regulatory authorities.

Historically, there have been two main incentives to build a quantum computer. First, there are physical limits on the miniaturization of electronic parts preventing any further increase in the speed of classical computers at some point in time. Second, and more importantly, some tasks are difficult to solve on classical computers by definition (e.g., the traveling salesperson problem and factorization of large integers). In such cases, quantum algorithms offer significant speed-up (in other words, have lower complexity) in comparison with their classical counterparts.

The idea of a quantum computer was first proposed in the seminal paper [1] as a tool for simulating quantum systems, because this task is difficult to carry out on classical computers -in Feynman's own words: Nature isn't classical, dammit, and if you want to make a simulation of nature, you'd better make it quantum mechanical, and by golly it's a wonderful problem, because it doesn't look so easy. The first algorithm outside theoretical physics was presented in [2] . It is able to recognize if a binary function is either constant or balanced 1 with only one evaluation of the function, regardless of the number of possible inputs. On a classical computer, by contrast, the function values for a half of its inputs have to be evaluated. Although this algorithm is a pure academic exercise, it shows that there are other tasks besides quantum system simulations where quantum computers can outperform classical ones.

The first practical quantum algorithm allowing integers to be factorized in polynomial time was introduced in [3] , while the best performing classical algorithms for factorization have exponential complexity. 2 Another important practical algorithm was suggested in [6] to achieve quadratic speed-up of searching in unordered databases. Grover's algorithm was later modified in [7] for finding the extremes of a function. Besides optimization problems, systems of linear equations often appear in scientific and business tasks. A quantum solver for systems of linear equations (the "HHL algorithm") was proposed in [8] . For certain types of linear equations, it offers exponential speed-up in comparison with classical solvers. This list of algorithms is far from exhaustive and serves mainly as a historical insight into the development of quantum computing and the possibilities it offers. The interested reader can consult an extensive overview of quantum algorithms in [9] . We will also discuss other algorithms in the second section of this paper.

The above-mentioned algorithms and modifications thereof can be used in several tasks in finance, for example, portfolio optimization [10] , risk management [11] and [12] , interest rate derivatives pricing [13] , and machine learning [14] . Summaries of quantum computer applications in finance are provided in [15] , [16] , and [17] . On top of that, Sveriges Riksbank has published an introduction to quantum computing for economists containing a long list of algorithms useful in finance and business generally [18] . Importantly, that paper also discusses the concept of tamper-proof quantum money.

In this paper, we will focus on the application of quantum computers in foreign exchange reserves management. In particular, we will discuss portfolio optimization and risk measurement. In the case of portfolio optimization, we will focus first on finding the optimal balance between bonds and equities and second on selecting equities so as to maximize portfolio returns while minimizing risk. Concerning risk measures, we will present an evaluation of VaR and CVaR and some other statistical properties. All demonstrations will be carried out on the IBM Quantum TM computer. IBM offers a cloud-based platform with access to a real quantum computer for free, provided that it is used for educational and research purposes. IBM's quantum computer is programmed with a Python-based language called Qiskit. An introduction to the language is provided in [19] and details can be found in the official Qiskit documentation.

In this research, we found that quantum computers are able to solve some tasks connected with FX reserves management, in particular risk measurement using the quantum Monte Carlo method and Markowitz-like portfolio optimization using a quantum linear equations solver (the HHL algorithm) and a quantum quadratic unconstrained binary optimization method (the quantum approximate optimization algorithm -QAOA). Unfortunately, the current early stage of quantum hardware development restricted us to toy models only (for example, we used the QAOA for a task with only five binary variables). The technical imperfections of the quantum computers currently available also meant that only a few of our demonstrations were successful. On the one hand, we successfully carried out portfolio optimization with the QAOA employing five binary variables for all the cases tested. On the other hand, the application of the quantum Monte Carlo method to risk measurement generated only partially correct results, and the use of the HHL algorithm in portfolio optimization failed completely. However, running the algorithms discussed above on a quantum computer simulator (i.e., software simulating a quantum computer on a classical one) proved that our implementations are correct, hence we came to the conclusion that in principle nothing prevents us from deploying quantum computers in FX reserves management in the future. We only have to wait for technical progress in quantum hardware capabilities.

The rest of the paper is organized as follows. In the second section, we discuss the theoretical background of quantum computing, ranging from the very basics to the algorithms we will employ in this paper. The third section contains details of the FX reserves management tasks described above and the results of solving those tasks on IBM's quantum computer. The fourth section concludes. The appendixes contain a short description of the IBM Quantum TM environment and the Qiskit language, the Qiskit source codes for all the demonstrations presented in this paper, and other expanding materials.

In this section, we discuss the basics of quantum computing theory. First, we show what quantum bits (qubits) are and how we can use them to perform useful calculations. We also introduce the terms superposition and quantum entanglement. Not only are these terms crucial for understanding quantum computing, they are also the reason why quantum computers outperform classical ones in some tasks. After that, we introduce the algorithms we will employ in our demonstration of the application of quantum computers in FX reserves management. The first is a quantum method allowing us to calculate risk measures such as VaR and CVaR based on the historical empirical distribution of portfolio returns. The second is a quantum solver of systems of linear equations (the HHL algorithm), which will be employed in Markowitz-like portfolio optimization. The last algorithm to be discussed is the quantum approximate optimization algorithm (QAOA), which enables us to solve quadratic unconstrained binary optimization problems (e.g., portfolio optimization and the famous traveling salesperson problem).

This section covers the necessary minimum for understanding the theory behind quantum computing. We explain what quantum bits are, how quantum operations (or quantum gates) work, and how quantum gates are composed to form quantum circuits representing quantum algorithms. For readers interested in gaining a deeper understanding, we recommend the classical monograph on quantum computing [20] and the introduction to computer science (including quantum computing) [21] . A non-technical introduction to quantum mechanics can be found in [22] . The mathematical background of quantum mechanics is again discussed in [20] .

The basic building block of a quantum computer is the quantum bit (qubit). Similarly to a bit on a classical computer (i.e., a non-quantum computer such as we use in daily life), which is considered to be the smallest unit of information, the qubit is the smallest piece of information in the quantum sense.

A qubit is a quantum system which can be in two distinguishable states. Such states are 0 and 1, as on a classical computer. However, a qubit can be in both states at once, written mathematically as

where α and β are complex numbers. Symbols |0 and |1 represent states 0 and 1, and |q is the whole qubit. The symbol |q comes from the Dirac bra-ket notation. In fact, |q ("ket") is a column vector, and it holds that

The Dirac notation is mainly used as an economical way of writing down vectors. To have a complete picture of the Dirac notation, we also introduce the symbol q| ("bra"), which is the transpose conjugate of |q , i.e., a row vector composed of complex conjugates of the elements of |q . With (2), we can write the qubit defined in (1) as a vector:

Vector |q is clearly a linear combination of vectors |0 and |1 -we say that qubit |q is a superposition of states |0 and |1 . 3 States |0 and |1 are called basis states, because they form the basis for the vector space C 2 . Parameters α and β are called probability amplitudes. It also holds that the square of the absolute value of the probability amplitude 4 is the probability that the qubit is in a particular state after measurement (we will discuss measurement in detail later). It therefore holds that P (|0 ) = |α| 2 P (|1 ) = |β| 2 .

Since |α| 2 and |β| 2 are the probabilities of mutually exclusive outcomes, it holds that |α| 2 + |β| 2 = 1. This also means that the length of the vector describing qubit |q is unity. Therefore, |q is a proper quantum state if and only if it has unit Euclidean length. As α and β are complex numbers satisfying |α| 2 + |β| 2 = 1, we can use the exponential form of a complex number (i.e., α = |α|e iϕα ) and the trigonometric identity sin 2 θ 2 + cos 2 θ 2 = 1 to rewrite (1) as

The term ϕ α is called the global phase and can be ignored, as two qubits differing in the global phase only are physically indistinguishable and are considered to represent the same state of the qubit. The term ϕ β − ϕ α is the relative phase and is usually denoted by ϕ. The phase is a purely quantum feature of the qubit and we will see how it is used later. For the time being, consider it simply to be a parameter of the qubit. Neglecting the global phase, we arrive at the most general form of the qubit:

As can be seen from (6), the qubit is described by two angles θ and ϕ. Therefore, any qubit can be visualized on a unit sphere called the Bloch sphere. As an example, visualizations of qubits in states 1 √ 2 (|0 + |1 ), i.e., having θ = π 2 and ϕ = 0, and 1 √ 2 (|0 + i|1 ) (θ = π 2 and ϕ = π 2 ) on the Bloch sphere are shown in Figure 1 .

As can be seen in Figure 1 , both states are superpositions of |0 and |1 , because the vectors representing them lie outside the z axis, where states |0 and |1 are positioned. The states differ in phase ϕ. While state 1 √ 2 (|0 + |1 ) lies on the x axis, the other state is located on the y axis, i.e., it is rotated by π/2 radians. Hence, this example shows a graphical representation of the quantum phase as a rotation of a vector on the Bloch sphere. As we will see later, rotation is one of the most important operations allowing calculations to be performed on a quantum computer. (|0 + i|1 ) is on the right. Source: Author's own creation on IBM Quantum TM Now that we have introduced the single-qubit states, we can move to multi-qubit states. As two bits can form four different combinations (00, 01, 10, 11), the general two-qubit state is

where |α| 2 + |β| 2 + |γ| 2 + |δ| 2 = 1. In vector form, the basis states can be written as

Obviously, there has to be a connection between a two-qubit state and its partial qubits, i.e., the single-qubit states. This connection is established through an operation called the tensor product.

We show the definition of the tensor product on an example. Consider the state |01 , which is composed of two single qubits, the first being |0 and the second |1 . We can therefore write |01 = |0 ⊗ |1 , where symbol ⊗ is the tensor product defined as follows:

It should be emphasized that the tensor product is not a commutative operation. For example, |1 ⊗|0 = |10 is a clearly different state than |0 ⊗ |1 = |01 . Similarly, it is possible to define states composed of three or more qubits. 5 Not all multi-qubit states are tensor products of other states. For example, the Bell state

cannot be written (separated) as a tensor product of two single qubits. Such inseparable states are called entangled states. These states occupy a special place, because, together with superposition, they are responsible for the speed-up provided by quantum computers. Another special feature of entangled states is that it is possible to infer the whole state from only a few of the qubits the state is composed of. For example, in the Bell state (10), if the first qubit is |1 , we know that the whole state is |11 . 6 Examples of other entangled, generally n-qubit states include

• GHZ states: 1 √ 2 (|00 . . . 0 + |11 . . . 1 ) -for n = 2 we get the Bell state (10)

• W states: 1 √ n (| 00 . . . 01 n qubits + |00 . . . 10 + · · · + |01 . . . 00 + |10 . . . 00 )

Note that GHZ and W states can be prepared using the algorithm presented in [23] . So far, we have considered that the qubit is in superposition. However, according to the rules of quantum mechanics, it remains so only until a measurement is carried out. Measurement is the act of getting information (or knowledge) about the qubit. We have said that the qubit is in both states at once but that when we measure it we get a particular value, either 0 or 1. This process is also called wave function collapse, as (1) is also referred to as a wave function and we see that the qubit has "collapsed" from the superposition to one particular state. 7 The result of measurement is a probabilistic event -for the single qubit described by (1) the value 0 is measured with probability |α| 2 and 1 with probability |β| 2 . Similarly, for the two-qubit state (7) we get 00 with probability |α| 2 , 01 with probability |β| 2 , 10 with probability |γ| 2 , and 11 with probability |δ| 2 . This raises the question of how we can perform any useful computation when the results are probabilistic. First of all, we have to get enough samples to reconstruct the probability distribution governing the quantum state saved in the qubits. This means that the state is prepared and measured several times (such repetitions are called shots in quantum computing). Naturally, this approach partially cancels out the increase in performance gained by using a quantum computer. However, it is worth noting that suppressing the error rate in quantum computers (we will discuss this issue in detail later) reduces the number of shots needed to reconstruct the probability distribution. Moreover, there are some quantum calculations which in theory need only one evaluation, hence the high performance of quantum computers is preserved. In Section 2.2, we will show how the probability distribution of the qubits can be manipulated and translated to useful figures. 5 It is apparent that an n-qubit state is described by 2 n complex numbers, which means that the amount of memory needed to simulate it on a classical computer increases exponentially with the number of qubits involved. This has been the main incentive for the construction of quantum computers, as the amount of quantum memory necessary increases only linearly in the number of qubits simulated. 6 Remember Schrödinger's cat. The bottle is another quantum system in a superposition of states |broken and |unbroken . Clearly, the cat's state depends on the state of the bottle -we only have the combinations |broken |dead and |unbroken |alive , so the systems are entangled. We can see a clear resemblance to the Bell state |00 + |11 . However, this analogy is not perfect, as the bottle influences the cat but not vice versa. 7 Again recall Schrödinger's cat. Opening the box is like making a measurement and leads to the collapse of the wave function describing the cat, i.e., we find out whether the cat is dead or alive.

In the previous section, we considered qubits to be static objects. However, to carry out useful calculations we need a tool for changing the states of qubits. This tool is called a quantum gate, and a collection of such gates representing a quantum algorithm is known as a quantum circuit. 8 In the rest of this paper, we assume -in accordance with the criteria for the physical implementation of a quantum computer introduced in [24] -that any quantum computer is able to prepare a qubit (or qubits) in state |0 as the initial state for subsequent calculations.

In general, the time evolution of a quantum state |ψ is described by Schrödinger's equation

where is the reduced Planck constant and H is a Hermitian matrix 9 called the Hamiltonian. The Hamiltonian, also called the energy operator, describes how a quantum system behaves. For example, the Hamiltonian's eigenvectors represent the allowed energy levels and its eigenvalues represent the values of the energy on these levels. We will discuss Hamiltonians in detail in Section 2.2.

In the case of quantum computers, we can assume that the time evolution goes in discrete time steps. As a consequence, the change from state |ψ(t 1 ) to state |ψ(t 2 ) is described by the equation

where U is a unitary matrix. 10 The matrix U in expression (12) is the mathematical representation of a quantum gate. The simplest gate is an identity operator leaving a qubit unchanged. 11 The identity operator is represented by a unit matrix

Later, we will see why this operator is important. The first gate actually changing the state of a qubit is a quantum NOT gate, which converts state |0 to state |1 and vice versa. The NOT gate is described by matrix

An effect of gate X applied to states |0 and |1 can be easily verified by the direct calculation

In the previous part, we mentioned that superposition plays an important role in the high computational power of quantum computers. To prepare a qubit in superposition, we can employ, for example, the Hadamard gate, represented by the matrix

Applying H to |0 and |1 , we obtain the following states:

There are other models of quantum computing, for example, the single-purpose quantum annealers used in binary optimization, which do not work with gates. However, in this paper, we describe a gate-based model of a universal quantum computer as adopted by IBM and other companies such as Honeywell and Microsoft. 9 A matrix A is Hermitian if A † = A, where A † is the conjugate transpose matrix of A, i.e., a transposed matrix with all the elements switched for their complex conjugated values. 10 A matrix A is unitary if AA † = A † A = I, where I is a unit matrix. Clearly, a unitary matrix always has an inverse matrix. This means that any change on a quantum computer can be reversed, the exceptions being measurement and reset (i.e., setting a qubit to state |0 ). 11 An identity operator is similar to an idle (or empty) instruction in classical programming languages.

States |+ and |− are clearly in a superposition of |0 and |1 . The probability of obtaining 0 and 1 after measurement is 50% in both cases, since (1/ √ 2) 2 = 0.5. The states |+ and |− differ only in the relative phase, which is π in the case of |− (because −1 = e iπ ) and zero for |+ (1 = e i0 ). It is possible to convert state |+ to state |− and vice versa with the Z gate given by the matrix

It can be verified by direct calculation that Z|+ = |− and Z|− = |+ . In the previous section, we introduced the most general form of the single qubit (6) . To prepare this state on IBM's quantum computer, a gate denoted U3 is used

Clearly, operation U3(θ, ϕ, 0)|0 prepares the qubit in state (6) . There are also two specialized gates based on U3, namely, U1(λ) = U3(0, 0, λ) and U2(ϕ, λ) = U3(π/2, ϕ, λ). Note that these gates are specific to IBM. However, in this paper we will only use the IBM Quantum TM environment. At this point, we are able to change the state of a single qubit. To work with more qubits, we apply the single-qubit gates introduced above to more qubits. The matrix representation of such gates is given by the tensor product of the gates on the individual qubits. We demonstrate this on the example of two qubits. The Hadamard gate is applied to each of them

If we apply gate H ⊗ H to two-qubit state |00 , we get a new state 1 2 (|00 + |01 + |10 + |11 ). The probability of measuring each combination of 0s and 1s is 25%, since (1/2) 2 = 0.25. The same goes for the other two-qubit inputs to gate H ⊗ H. The states differ only in their relative phases.

Following this pattern, we are able to prepare any separable multi-qubit states, i.e., those that can be written as a tensor product of individual qubits. However, to exploit the high computational performance of quantum computers, we need a tool for preparing entangled states. This can be achieved with the aid of controlled gates. An example of such a gate is the two-qubit controlled NOT gate, or CNOT. When the first qubit (control qubit) of the gate is in state |0 , the second qubit (target qubit) is left unchanged. But when the control qubit is in state |1 , the NOT gate is applied to the target qubit. Note that the control qubit is unchanged. Consequently, the map of the two-qubit basis input states to the outputs has the form: |00 → |00 , |01 → |01 , |10 → |11 , |11 → |10 . 12 The CNOT gate can also be described by the matrix

We can see that the matrix X for a single-qubit NOT operation is hidden in (20) , because CNOT =

where O 2,2 is a two-by-two zero matrix. This is a general pattern for any two-qubit controlled gate. Similarly, we can write matrices for controlled Z, controlled U3, and other controlled gates.

A quantum gate can have more than one control qubit. An example of such a gate is the three-qubit Toffoli gate. In fact, the Toffoli gate is a double-controlled X gate, hence it is often denoted CCNOT. The third qubit is negated when both control qubits are in state |1 . 13 The matrix describing the Toffoli 12 Note that the output value of the target (second) qubit is a classical XOR function applied to both the control and target qubits, while the output value of the control (first) qubit is a copy of its input. 13 When the third qubit is in state |0 , after application of the Toffoli gate it contains a value equal to the classical AND function applied to both control qubits. If the third qubit is in state |1 , after the Toffoli gate is applied it contains the result of the NAND function applied to the control qubits. gate has following form: 

In the following text, we will introduce other gates. An overview of all the gates used is provided in Appendix B.

There are a number of relations between quantum gates. Complex multi-qubit and multi-controlled gates can be decomposed into simpler ones. An overview of those relations and decomposition techniques is provided in [25] . It is also worth noting that any quantum gate can be composed of members of a basic gates set, i.e., a small group of gates which are universal and enable any calculation to be performed on a quantum computer. 14 The algorithm for constructing a gate from the basic set is called the Solovay-Kitaev algorithm and is described in detail in [26] .

The composition of the quantum gates is called the quantum circuit, and the circuit represents a quantum algorithm. Here, we provide an example of a quantum circuit for the preparation of the Bell state (10) . 15 The circuit is composed of a Hadamard gate located on the first qubit, followed by a CNOT gate. While there is a Hadamard gate on the first qubit, there is nothing on the second one. This nothing is in fact the identity operator I (here we see the importance of the identity operator). The matrix form of the circuit is therefore

Clearly, for input |00 , the circuit returns the Bell state (10) . For inputs |01 , |10 , and |11 , we get the other Bell states |β 01 = 1 √ 2 (|01 +|10 ), |β 10 = 1 √ 2 (|00 −|11 ), and |β 11 = 1 √ 2 (|01 −|10 ), respectively. Quantum circuits are often visualized with diagrams. Examples of such diagrams describing circuit (22) and the output states for different inputs are provided in Figure 2. 14 Similarly, for classical computing, the logic NAND function constitutes a single-member basic functions set allowing a classical computer to perform any calculation. As the Toffoli gate can behave like a NAND function, we are able to implement any classical algorithm on a quantum computer, too. Hence, a gate-based quantum computer is universal from both the classical and quantum perspectives. 15 This circuit is sometimes considered to be the quantum analog of a Hello, World! program.

Note: The red rectangle with an H symbolizes a Hadamard gate, the blue circle is an X gate, and the blue circle with a line above it is an CNOT gate. Only H and CNOT are employed directly in circuit (22) . The X gate is used for setting the input state; for example, the input in the top right circuit is |01 , hence the X gate is located on the second qubit. The other inputs are |00 in the top left circuit, |10 in the bottom left circuit, and |11 in the bottom right circuit. The colors of the bars in the graphs represent the relative phase -blue for 0 and red-brown for π. Note that the order of the matrices in (22) is opposite to the order of the gates in the diagrams. Also note that the results were obtained on a quantum computer simulator (i.e., software simulating the behavior of a quantum computer on a classical one). Source: Author's own calculations on IBM Quantum TM simulator

A swap test is an algorithm allowing us to calculate the absolute value of the inner product of two quantum states |ψ and |φ , i.e., the value of | ψ|φ |. This algorithm is useful for post-processing of the results obtained when solving linear systems in general and when optimizing portfolios in particular. The algorithm is based on a three-qubit controlled swap gate (or Fredkin gate). If the control qubit of the gate is in state |1 , the states of the two target qubits are swapped. In other words, the gate's action is similar to switching two cables between two sockets. The matrix describing the Fredkin gate is the following (note that the red part of the matrix represents the uncontrolled version of the swap gate): 

A schematic of the swap test is shown in Figure 3 . For general n-qubit states |ψ and |φ , the "SWAP" block contains n controlled swap gates with targets on the corresponding qubits of states |ψ and |φ . All these swap gates are controlled with the uppermost qubit, which is initially set to state |0 . Source: Adapted from [27] The probability that the uppermost qubit is in state |0 is 16

hence for the inner product we have the formula | ψ|φ | = 2P (|0 ) − 1.

The main disadvantage of the swap test is its inability to determine the sign of the inner product for states with real probability amplitudes. Moreover, if the states have complex probability amplitudes, the information about the real and imaginary parts is lost. However, for the purposes of this paper, the method is sufficient, because the asset weights in portfolio optimization are always real, and if short positions are forbidden, the weights are non-negative.

In Figure 4 , we provide a practical implementation of the swap test on IBM Quantum TM for Bell states |β 00 and |β 01 . In this case ψ|φ = 0, because the states are orthogonal, hence the probability of measuring |0 in the control qubit is 50%.

In this part, we discuss the quantum Fourier transform (QFT), a very important part of many quantum algorithms. For example, in a quantum linear equations solver, the QFT is used for extracting the eigenvalues of a matrix, which are then used in subsequent calculations. Note that from this point on, for the sake of simplicity, we sometimes omit the symbol ⊗ in tensor products. For example, the expression |0 |1 has the same meaning as |0 ⊗|1 . For an n-qubit state having each qubit in either basis state |0 or basis state |1 , we introduce symbol |x , where x ∈ {0, 1, . . . 2 n − 1}, denoting the decimal representation of the binary number stored in the n-qubit state. For example, for n = 3, the expression |5 is equivalent to |101 = |1 |0 |1 . 17 We also introduce the notation 0.x 1 x 2 . . . x n , where x i ∈ {0; 1}, representing a binary fraction, i.e., 0.x 1 x 2 . . . x n is equivalent to x 1 /2 + x 2 /4 + · · · + x n /2 n in the decimal digits system.

Equipped with these new symbols, we define the quantum Fourier transform for the n-qubit basis state |x as an operation:

As this exact definition is not easy to read and understand, an equivalent formula is used in practice:

Input state |x can be rewritten as |x 1 |x 2 . . . |x n , where |x i ∈ {|0 ; |1 }, hence we see that the QFT encodes the input binary string and its sub-strings into quantum phases of the output qubits.

As the QFT is a quantum operation, it has an inverse -the inverse quantum Fourier transform, abbreviated as QFT −1 or QFT † . From (25) it follows that the inverse operation extracts the quantum phases of the input qubits and translates them into an output bit string. As mentioned in Section 2.1.1, some quantum algorithms encode the result of the calculation into the quantum phase. To get this result, we employ QFT −1 to translate the phase into a bit string.

A general circuit implementing the QFT on a quantum computer is shown in Figure 5 . The quantum gates R k in Figure 5 are defined as

The circuit for inverse quantum Fourier transformation is the circuit in Figure 5 taken from right to left, with all the rotation angles replaced by their opposite values, i.e., 2π/2 k in gate R k is substituted with −2π/2 k . The Hadamard gate is its own inverse, hence the H gates are left unchanged. Source: Adapted from [20] Note that R k is equivalent to gate U1(λ), where λ = 2π 2 k . Examples of three-qubit QFT and inverse QFT implementation on IBM Quantum TM are shown in Figure 6 . We also provide the Qiskit source code for QFT and QFT † in Appendix F. 

In this part, we will describe quantum algorithms that are useful for FX reserves management. First, we will discuss an algorithm for calculating risk metrics such as VaR and CVaR. Second, we will introduce algorithms for solving systems of linear equations and for finding optimum solutions to quadratic unconstrained binary optimization problems. We will then employ both algorithms to construct a portfolio.

Calculating risk indicators is an integral part of portfolio management. For this reason, we introduce the quantum algorithm for calculating VaR and CVaR proposed in [11] and [12] . An important benefit of this algorithm is that it provides quadratic speed-up in comparison with classical calculation techniques based on Monte Carlo simulation. In plain language, this means that we need, for example, only 100 operations instead of the 10,000 necessary on a classical computer. In more mathematical terms, while the classical Monte Carlo method converges as O(M −1/2 ), the quantum version does so as O(M −1 ), where M is the number of samples used.

In the rest of this section, we assume that we are dealing with discrete probability distributions with a finite number of outcomes. Without loss of generality, we also assume that the number of probability distribution outcomes is 2 n , where n ∈ N. 18 Where a continuous distribution is involved, discretization is carried out first.

Assume that X is a random variable with a discrete probability distribution having a probability function p satisfying

If we apply such operation on state |ψ |0 , we arrive at a new state

Clearly, the probability of measuring |1 in the last qubit of |ψ is the sum of the probabilities of all outcomes having shape |i |1 , i.e.,

which is the expected value E[f (X)] of the random variable f (X). Let us denote N = 2 n and f (x) = x/(N − 1). In this setting, (28) changes to P 1 = 1

. This means that we are able to translate the calculation of the expected value of a random variable to the measurement of a quantum state.

If

) 2 ) and the standard deviation, i.e., σ(X) = VAR(X), of the random variable X. A similar pattern in setting f (x) can be followed to get higher moments (for example, the skewness and kurtosis) of the random variable X.

A special setting of f (x) also enables us to evaluate a percentile of the random variable X. Since Value-at-Risk (VaR) is defined as VaR α (X) = min{x|P [X ≤ x] ≥ 1 − α}, it is in fact the 1 − α percentile. If, for a certain l ∈ N, we set f (i) = 1 for i ≤ l and f (i) = 0 otherwise, equation (28) transforms to P 1 = i≤l p i . Setting the initial value l = N − 1, we start an interval bisection and iterate until the probability of measuring state |1 in the last qubit of |ψ is higher than 1 − α. Once the iterative process stops, we save the value l, which is the desired 1 − α percentile or VaR α (X). Note that when evaluating VaR, we combine quantum and classical (i.e., interval bisection) calculations. This approach is called a hybrid algorithm.

Conditional VaR (CVaR) is defined as the expected value of outcomes below VaR, i.e.,

This means that we can employ an approach similar to the calculation of E[X]. We set f (i) = i/VaR α (X) for i ≤ VaR α (X) and f (i) = 0 otherwise. In this setting, (28) becomes

Since i≤VaRα(X) p i < 1, i.e., the probability distribution of values below VaR is not normalized, the probabilities p i have to be replaced by p i /P [X ≤ VaR α (X)]. Normalization, together with equation (30), leads to the final formula for calculating CVaR on a quantum computer:

18 If the number of outcomes is not equal to a power of two, the distribution can be extended by including additional outcomes with zero probability so that the total number of outcomes is equal to a power of two.

Note that VaR α (X) is known from the previous calculations and P [X ≤ VaR α (X)] = i≤VaRα(X) p i .

To carry out the calculations described above, we have to prepare state |ψ . We start by preparing the part |i ( 1 − f (i)|0 + f (i)|1 ). To do so, we design a gate with n controlling qubits changing the state of the n + 1 th qubit according to an n-qubit state |i and function f (x). Clearly, the uncontrolled version of this gate (i.e., with n = 0) is a single-qubit gate preparing state

Such gate is called y-rotation 19 and is described by matrix

Evidently, Ry(α)|0 = cos(α/2)|0 + sin(α/2)|1 . Setting α = 2 arccos(

If we calculate the rotation angles α for all f (i), we get a sequence α i = 2 arcsin[ f (i)] for N Ry(α i ) gates. However, we have to add controls to these gates to set the proper angle α i for state |i on the control qubits. To do so, we employ the method described in [28] . An example of a circuit implementing the preparation of |i

i.e., eight possible inputs |i , is shown in Figure 7 .

We easily see how the Ry and CNOT gates are arranged. Following this pattern, we can expand the circuit for any number of control qubits. It is important to emphasize that the rotational angles of the Ry gates are not α i based on f (i) as calculated before, but other angles θ i obtained by means of a transform. The reason for changing the angles from α i to θ i is that there are no native n-qubits controlled the Ry gate, but such qubits can be implemented with single-qubit Ry gates and two-qubit CNOTs, as shown in Figure 7 . Therefore, the rotational angles have to be adapted to this substituting configuration. The transformation from angles α i to angles θ i is given by the expression

where M ∈ R N,N with elements

where b(j) is the binary representation of column index j, g(i) is the Gray code representation of row index i (a method for converting a binary number into Gray code is described in Appendix C), and operation • is a bit-wise product modulo 2 defined as

where b k (j) is the k th bit of binary number b(j) and similarly for g(i).

To give an example of how the method described above works, we prepare a circuit for measuring the expected value of a distribution with eight outcomes (i.e., three qubits are involved). In this case, f (i) = i/7 and transformation matrix M 8 is

Function values f (i), rotational angles α, and transformed angles θ are shown in Table 1 . To prepare the whole state |ψ , we need state |ψ describing the probability distribution of random variable X. This state can again be obtained using the method introduced in [28] based on n-qubit controlled rotations discussed above. The interested reader can consult that paper to get an insight into how the method works. We also provide the Qiskit code implementing the method in Appendix F. Alternatively, it is possible to use the implementation of this method in the initialize function obtained from the Qiskit libraries. 20

Source: Adapted from [28] We now summarize the practical approach to calculating risk metrics on a quantum computer. The first step is to prepare an n-qubit state |ψ based on the probability distribution of random variable X. Then, a new qubit in state |0 is added to state |ψ to get an n + 1-qubit state |ψ = |ψ |0 . This new state is inputted into a circuit similar to that shown in Figure 7 and adapted according to the number of outcomes of random variable X. Finally, the n + 1 th qubit of state |ψ is measured and, based on the probability of measuring state |1 , the actual value of the risk metric is calculated classically. 

Systems of linear equations are often used in modeling and optimization. Harrow et al. designed a quantum algorithm for solving linear systems, commonly known as the HHL algorithm [8] . 21 . Later an application of the HHL algorithm in Markowitz-like portfolio optimization, which we will employ in the practical part of this paper, was proposed in [10] . The main advantage of the HHL algorithm is that it provides exponential speed-up in comparison with classical solvers. However, in some cases there is no speed-up at all, as emphasized in [30] and [31] (we will return to this discussion at the end of this part).

Note that we will only present an outline of the algorithm, as very technical details can hinder understanding of how it works. Moreover, in the practical part, we will use a debugged and optimized implementation of the algorithm taken from the Qiskit libraries. The interested reader can find technical 20 Note that a more sophisticated approach to preparing the initial state |ψ based on a probability distribution inferred from a huge dataset with quantum machine learning was proposed in [29] . 21 The abbreviation comes from the names of the creators of the algorithm, Harrow, Hassidim, and Lloyd. details in the above-mentioned papers. Examples of how to build a quantum circuit implementing the HHL algorithm for a particular linear system are provided in [32] , [33] (suitable for complete beginners), and [9] . Suppose we want to solve a system of linear equations A|x = |b , where A ∈ C N,N is a full-rank matrix, |b ∈ C N is a non-zero vector, and |x ∈ C N . Note that we use Dirac notation to make the formulas succinct. Without loss of generality, we assume that matrix A is Hermitian. 22 Such matrices have real eigenvalues, and their eigenvectors form the orthonormal basis of the vector space C N . 23 Moreover, Hermitian matrices always have a spectral decomposition, defined as 24

where λ i is the i th eigenvalue of matrix A and |u i is the respective eigenvector. With the spectral decomposition, we can introduce a matrix function f

Setting f (x) = x −1 , we can write the inverse matrix of A (we assume that the matrix is full-rank, hence it is invertible) as

Since the eigenvectors of matrix A form the basis of space C N , the right-side |b can be written as a linear combination of those eigenvectors |b = 25 With the help of (39), the solution of the linear system is

Because the eigenvectors are orthonormal, u i |u j = 0 for i = j and u i |u j = 1 for i = j, hence the solution of the linear system is

To solve A|x = |b on a quantum computer, we need to prepare solution (41) as a quantum state. First, we have to find eigenvalues λ i . To do so, we employ a phase estimation algorithm, which finds the phases of the unitary matrix eigenvalues. Since A is Hermitian, matrix U = e iA is unitary. 26 According to (38) , the eigenvalues of U have the form e iλj . This means that the eigenvalues of A are encoded into phases of the eigenvalues of matrix U. The first step of the phase estimation algorithm is shown in Figure 8 . 22 If A is non-Hermitian, we can switch from system

. Evidently, the matrix of the new system is Hermitian and part |x of the new solution is the same as the original one. 23 To be precise, the eigenvectors associated with distinct eigenvalues are orthogonal. However, it is always possible to get orthonormal vectors by normalizing such vectors by their Euclidian norm. Moreover, it holds that even in the case of degenerated eigenvalues, it is possible to find the orthogonal basis of C N composed of the matrix eigenvectors. 24 Any normal matrix A, i.e., one satisfying AA † = A † A, has a spectral decomposition. 25 Note that the coordinates of |b in the basis composed of the eigenvectors of A can be found by solving a linear system U|β = |b , where U is a matrix, the columns of which are created from the eigenvectors of A. Since these eigenvectors are orthonormal, U is unitary, hence |β = U † |b . The i th element of |β is the inner product of |b and the i th row of U † , which is a complex conjugated vector to the i th column of U (i.e., u i |). This gives us the expression β i = u i |b . 26 It holds that

Note: This part of the HHL algorithm is intended for finding the phases of the unitary matrix eigenvalues. Given a unitary matrix U and an eigenvector |u , the algorithm encodes phase θ of eigenvalue e iθ into quantum phases of the output qubits, while eigenvector |u remains unchanged. The number of qubits t depends on the required precision of the binary representation of phase θ. Since U is e iA in this case, phase θ is in fact eigenvalue λ of matrix A. Note that symbols U k represent the k th power of matrix U. Source: Adapted from [20] As can be seen in Figure 8 , after application of the first step of the phase estimation, eigenvalue λ corresponding to eigenvector |u of matrix A is encoded into phases of the output qubits. To obtain the phase (i.e., eigenvalue λ), we apply inverse QFT. If we put the right-side |b instead of eigenvector |u on the phase estimation input, we get a superposition of all the eigenvalues of A after inverse QFT, because |b is a linear combination (superposition) of all the eigenvectors of A. 27 The phase estimation does not change the input eigenvectors, so |b is not changed either. This all means that we have quantum state N i=1 β i |λ i |u i after the phase estimation. 28 To get the inverse of the eigenvalues, we add an auxiliary qubit (ancilla qubit) and apply controlled Ry gates with the target on the ancilla qubit and controlled with qubits containing the eigenvalues. The rotational angles of the Ry gates are set to θ i = 2 arcsin(C/λ i ), where constant C ensures that the resulting quantum state is normalized. 29 After application of the Ry gates, we have a new quantum state

The next step of the algorithm is to apply inverse operations to all the operations carried out so far, the Ry gates and the preparation of state |b being exceptions. As a result of these actions, all the qubits containing eigenvalues before are now in state |0 and we get quantum state

The final step of the algorithm is to measure the ancilla qubit. If the result of the measurement is |1 , neglecting qubits in state |0 . . . 0 , we come to state

which is equivalent to (41) up to the normalization constant C. This means that we do not get |x itself, but its normalized version. If the result of the ancilla qubit measurement is |0 , the algorithm has to be run again until |1 is measured, because only in this case is the solution of the linear system returned.

For better illustration, a schematic of the whole HHL algorithm is shown in Figure 9 .

At the beginning of this section, we noted that in some cases the HHL algorithm generates no speed-up over classical linear system solvers. In particular, it cannot outperform classical solvers if we want to know all the members of solution |x . Obtaining the solution involves measuring all n = log 2 N qubits representing the solution. Since n qubits can be in 2 n basis states, complete measurement is exponentially complex and entirely cancels out the speed-up provided by the HHL algorithm. 30 However, quantum state |x can be used as an input to another calculation. For example, if we are interested in the inner product of |x and another quantum state, we employ a swap test converting the inner product calculation to the measurement of only one qubit. As a result, the exponential speed-up of the HHL algorithm is preserved.

However, this is not the only pitfall of the HHL algorithm. The complexity of the algorithm is O[log(N )sκ/ ], where s is the sparsity of matrix A (telling us how many elements in a matrix row are non-zero on average), κ = |λ max |/|λ min | is the condition number, and is the desired accuracy of the solution. This means that the HHL algorithm only works well for sparse matrices, i.e., those where the majority of their elements equal zero. At the same, the matrices should be well conditioned, i.e., the difference between the highest and lowest eigenvalues should be low. Ideally, the condition number should be close to 1. 31 These requirements are not always satisfied, and this hinders the performance of the HHL algorithm. On top of that, if we require a high-accuracy solution, we again lose part of the performance of the HHL algorithm.

However, we now turn our attention to the application of the HHL algorithm for portfolio optimization (we will discuss how to deal with the described drawbacks later). Assume that we want to minimize the risk in a portfolio, achieve a return G, and spend the whole budget B. Let C ∈ R n,n be a covariance matrix of asset returns, |R ∈ R n a vector of expected asset returns, |P ∈ R n a vector of asset prices, and |w ∈ R n a vector of asset weights. The task can be formulated as min |w w|C|w , subject to R|w = G and P |w = B. Such optimization leads to minimization of the Lagrange function

where λ and µ are Lagrange multipliers. Since function L is quadratic, its derivatives are linear functions, and the optimization can be translated to the solution of a linear system 32

30 Note that the method used to reconstruct the whole quantum state is called quantum state tomography. The interested reader can find details in chapter 8.4.2 of [20] . 31 If a matrix is unitary, all its eigenvalues are located on a unit circle, hence |λ| = 1 for any eigenvalue and therefore κ = 1. 32 Note that ∂L ∂λ = R|w − G and ∂L ∂µ = P |w − B. Since covariance matrix C is symmetric (i.e., c ij = c ji ), the derivative with respect to w i is ∂L

For all weights, the derivative can be written as ∂L ∂|w = 2C|w + µ|P + λ|R . Setting all derivatives equal to zero leads to the linear system (46).

Note: In the first step, Hadamard gates are applied to t qubits (t depends on the desired precision of the representation of the eigenvalues λ i ). This is symbolized by gate Ht. After that, the controlled version of the powers of U = e iA is applied (symbolized by e iA ), followed by inverse QFT (denoted by QFT † ). These steps constitute the phase estimation extracting the eigenvalues of matrix A. Based on the eigenvalues λ i , rotational angles θ i are calculated (operation U λ ) and controlled Ry(θ) gates are applied to the ancilla qubit (the uppermost qubit). The next step (operation U † ) involves inverting the operations U λ , QFT † , e iA , and H k . Finally, the ancilla qubit is measured, and if the measurement result is |1 the right-side |b has been changed to solution |x . Source: Adapted from [32] System (46) can be solved with the HHL algorithm. However, as mentioned before, extracting the whole solution cancels out the speed-up. To exploit the high performance of the HHL algorithm, we have to use the portfolio optimization results (denoted |x opt ) in other ways, for example as suggested in [10] :

• To get the optimal exposure to a particular asset or sector: Assume that we are interested in the weight of the i th element of |x opt . To obtain it, we calculate the inner product of |x opt and a vector having its i th member equal to 1 and the rest as zeros. The latter vector is in fact a quantum state |i − 1 which can be easily prepared with X gates only. 33 Note that the first two members of state |x opt are in fact Lagrange multipliers λ and µ, therefore the first weight of an asset is the third member of |x opt .

• Compare the optimal solution with the current portfolio composition: To do so, we need to calculate the inner product S = x opt |x current , where |x current = N i=1

x current i |i − 1 is a quantum state representing the current composition of the portfolio. If S = 1, then |x opt = |x current , which means that our portfolio is optimal. In contrast, if S = 0, the current and optimal portfolios are orthogonal, i.e., totally dissimilar. Note that the Lagrange multipliers λ and µ should also be members of state |x current . Since we do not know the multipliers in advance, we first have to obtain them using the method described in the previous bullet point (i.e., we need to extract the first two members of x opt ). This may slightly attenuate the performance of the algorithm. However, the number of asset weights is much higher than two, hence we still get quantum speed-up.

Unfortunately, using the solution described above does not recover the whole speed-up. A quick check of system (46) reveals that its main component is a covariance matrix which is hardly sparse. As a result, s = N and the complexity of the HHL algorithm for portfolio optimization increases to O[N log(N )κ/ ]. This is still better than the best known classical algorithm introduced in [35] , which has complexity O(N 2.376 ), and the more often used algorithm designed in [36] , which has complexity O(N 2.8074 ). 34 33 Assume, for example, that our optimization involves 16 variables (i.e., four qubits are employed to represent the optimization results) and we are interested in the seventh one. This means that we need to calculate the inner product xopt|6 (the index runs from zero, hence the seventh position is represented by the number 6). In binary notation, 6 is represented by the string 0110 (the leading zero comes from the fact that we are dealing with four qubits), hence |6 = |0110 . Such state is prepared by applying X gates to the second and third qubits of the initial state |0000 . 34 The most popular algorithm for solving linear systems, the Gaussian elimination algorithm, has complexity O(N 3 ).

This part is dedicated to the use of quantum computers for solving a quadratic unconstrained binary optimization (QUBO) task. QUBO has several business applications, such as the traveling salesperson problem (TSP) and the knapsack problem. 35 Some QUBO problems can be adapted for application in finance. For example, the TSP is used for finding currency arbitrage opportunities. Importantly for FX reserves management, QUBO is also useful in portfolio optimization. In this part, we introduce the QUBO task quantum solver called the quantum approximate optimization algorithm (QAOA), designed in [38] . In the practical part, we employ the Qiskit implementation of the QAOA, so we only provide an outline of the algorithm. After introducing the QAOA, we will discuss binary portfolio optimization based on the work [39] , which we will use in the practical part of this paper. Note that besides the QAOA, there are other QUBO solvers, for example, the variational quantum eigensolver discussed in [40] and the solver based on Grover database search introduced in [41] . On top of that, D-Wave has designed a single-purpose quantum computer (a quantum annealer) which implements a quantum QUBO solver on the hardware level. 36 It is worth emphasizing that in contrast to the algorithms discussed previously, there is no proof that quantum QUBO solvers are superior to classical solvers. We will return to this issue at the end of this section.

In general, QUBO solvers aim to minimize the function

where x i ∈ {0; 1} is a binary variable and a ij , b i , and c are real parameters. With the help of the Dirac notation, function (47) is expressed in matrix form f (x 1 , . . . , x n ) = x|A|x + b|x + c. Note that QUBO is unconstrained optimization. Where constraints are involved, they are incorporated in the form of a penalization function (we will see this later). On a quantum computer, the QUBO task can be translated to finding the ground state (i.e., the minimum energy level) of a special class of Hamiltonians (Ising Hamiltonians) having the form 37

Term Z i denotes a quantum gate composed of a Z gate applied to the i th qubit and identity operators applied to the other qubits, and term Z i ⊗ Z j means a quantum gate composed of two Z gates applied to the i th and j th qubits while identity operators are applied to the other qubits. Q ij and c i are real numbers. There is a clear resemblance between the QUBO task formulation (47) and the Ising Hamiltonian (48) and they can be converted into each other. The energy level associated with a quantum state |x in a quantum system described by a Hamiltonian H is given by the expression x|H|x . The ground state |x ground minimizes this expression. Since our QUBO task has been converted into an Ising Hamiltonian, we are looking for its ground state, which in the end can be converted back into the optimal solution of the QUBO task. To find the ground state of the Ising Hamiltonian, we employ the QAOA. The algorithm is based on the adiabatic theorem, which states that if a quantum system is described by an initial Hamiltonian H 0 , the system is in its ground state, and when the system's Hamiltonian is slowly changed, the system remains in its ground state but this new ground state is associated with a new Hamiltonian. This means that having a Hamiltonian with a known ground state, we can find the ground state of H Ising . Such initial Hamiltonian is, for example,

where X i denotes a quantum gate composed of an X gate applied to the i th qubit and identity operators applied to other qubits. The ground state of Hamiltonian (49) is a uniformly distributed superposition

|i , which can be prepared by applying Hadamard gates to all n qubits formerly in state |0 35 Other interesting QUBO problems are discussed in [37] . 36 More information about this device can be found in the user manual [42] . 37 Ising Hamiltonians are connected with a special class of magnetic materials called spin glasses. Such materials lie between paramagnetic and ferromagnetic materials. While ferromagnetics (e.g., iron) preserve their magnetization even after the external magnetic field is switched off, paramagnetics (e.g. aluminum) lose their magnetization rapidly (exponentially). Spin glasses are similar to paramagnetic materials but lose their magnetization more slowly (almost linearly).

(such action is often denoted as H ⊗n |0 ⊗n ). A simulation of a slow change from H 0 to H Ising can be carried out using a quantum gate defined as U (β, γ) = p j=1 e −iβj H0 e −iγj HIsing (50) applied to state H ⊗n |0 ⊗n . 38 Parameter p is called the algorithm depth and is set by the user. Changing p can lead to better optimization results. Parameters β j and γ j have to be varied so that x|H Ising |x is minimized, which is done using a classical optimizer. We clearly see that the QAOA is a hybrid algorithm. The quantum part of the QAOA involves simulating the Hamiltonian encoded in gate (50) . The classical part deals with the calculation of the energy levels of H Ising , which can be carried out efficiently thanks to the special structure of Ising Hamiltonians. 39 After some iterations, the QAOA should find the optimal solution of the QUBO task. However, optimality of the solution is not guaranteed, because the algorithm can get stuck in a local minimum. An example of a QAOA circuit for a task with two binary variables and algorithm depth p = 1 is shown in Figure 10 . The gates denoted as Rx(α) and Rz(α) are rotations around the x and z axes, respectively, of the Bloch sphere and are described by the matrices

The Rz gates represent exponentiation of the Ising Hamiltonian composed of Z gates, and the Rx gates are connected with the initial Hamiltonian H 0 composed of X gates. 40 Parameters a ij , b i , and c in QUBO task (47) and β j and γ j in (50) are included in the rotational angles of the Rx and Rz gates.

As mentioned at the beginning of this part, there is no rigorous proof of the superiority of quantum QUBO solvers (including the QAOA). According to [43] , there is little hope that a proof will ever be found, but for some tasks, quantum QUBO solvers may be faster than classical ones. For example, the studies [44] concerning the prediction of economic downturns and [45] relating to maritime routing problems, provide empirical evidence that quantum QUBO solvers perform better. Additionally, [46] introduced the non-adiabatic quantum QUBO solver, which should be faster for tasks with no specific internal structure. While the QAOA operates with Hamiltonian H = (1 − t)H 0 + tH Ising , where t is slowly changed from 0 to 1 (the "adiabatic regime"), the non-adiabatic approach uses a Hamiltonian of the form H = H Ising − g t α H init , where g and α are user-defined parameters, H Init is the initial Hamiltonian 41 and t is slowly increased from zero theoretically to infinity. This work is so far at a theoretical stage and no practical implementation on a real quantum computer has been presented yet. On the other hand, [47] revealed that in many instances QUBO problems are specially tailored (in this case to the topology of D-Wave quantum processors) to be hard for classical solvers and easier for quantum ones. On top of that, the authors of the study criticized the fact that quantum QUBO solvers are often only compared with, in principle, related classical simulated annealing and/or derived algorithms, 42 while several other classical heuristics are neglected. They also showed that in an assignment problem, a classical exact minimum-weight perfect matching algorithm 43 outperforms the D-Wave quantum processor for instances with more than 20 variables. Overall, we should use quantum QUBO solvers as a complement to, not a substitute for, classical algorithms. This means we should try to employ both classical and quantum heuristics and, in the end, pick the best solution found regardless of which approach produced it.

The QAOA will be practically demonstrated on Markowitz-like portfolio optimization. Let us denote as x i a binary variable indicating whether the i th asset is a member of the portfolio (x i = 1) or not (x i = 0). In contrast to the optimization described in the previous part, we do not look for particular weights of assets, but only for information on whether to invest in some asset or not. Assume that the return on the i th asset is r i , the covariance of the i th and j th asset returns is c ij , M i is the maximum amount of money we are willing to invest in the i th asset, and B is the total budget. Our aim is simultaneously to maximize the portfolio return n i=1 r i x i , to minimize risk n i=1 n j=1 c ij x i x j , and to 38 Note that gate (50) is decomposed into basic gates using the method presented in [40] . 39 The method is outlined in [40] . 40 It holds that Ra(α) = e −i α 2 A , where A is either X or Z. 41 H Init is generally different from H 0 used in the QAOA and its form depends on the task being solved. 42 Simulated annealing was introduced in [48] and was originally employed in the unfortunate race to develop nuclear weapons (one of its creators was Edward Teller, father of the H-bomb). 43 See the detailed description of the algorithm in [49] .

invest the whole budget, i.e., n i=1 M i x i = B. These three aims have importance expressed by coefficients λ 1 , λ 2 , and λ 3 , respectively. The equivalent QUBO task is to minimize the function

The minus sign before the first term (return) indicates that this part of the function has to be maximized. The third term is a penalization function expressing the budget constraint. The square in the third term ensures that the constraint is fulfilled in both directions. As discussed in [50] , even binary optimization allows one to work with integer asset weights instead of only binary flags indicating whether to include an asset in a portfolio or not. The integer weights represent, for example, the amount of money invested in a particular asset or percentage weights. In order for integer variables to be incorporated into the QUBO task, they have to be expressed with binary variables as follows:

where x ij is a binary variable representing the j th bit in the binary representation of the i th weight.

It is important to emphasize that the number of bits used in representation (53) constitutes a natural ceiling on weight w i , because the maximum integer value represented with m bits is 2 m − 1. Hence, an upper limit can be introduced on a weight without the use of a penalization function. However, note that integer optimization consumes a good deal more qubits than the binary case.

In this part, we focus on applying the algorithms introduced in the previous section to FX reserves management problems. First, we will measure the VaR and CVaR of a chosen portfolio and compare the results with those obtained using classical methods. Second, we will apply the HHL algorithm to find the optimal ratio between equities and fixed-income instruments. Third, we will use the QAOA to find the best performing assets in an index while minimizing risk. In both cases of portfolio optimization, we will again compare the results obtained on a quantum computer with classical methods. In our demonstration, we will employ both real quantum computers and their simulators. We will also discuss the state of development of the underlying quantum technologies and their current limitations for practical deployment.

Before we proceed to actual applications of quantum computers to FX reserves management tasks, we have to discuss the limitations and hurdles that quantum computers suffer from in their current state of development. We will discuss the consequences of noise, the limited connectivity of qubits, and the unavailability of quantum memories. It is important to emphasize that in contrast to the inherent limitations of some of the algorithms discussed earlier (such as the HHL algorithm), the said hardware issues are merely questions of further research and development and will be overcome in the end. Therefore, we will also present a metric used to benchmark quantum processors, allowing us to compare processors with each other and track their development over time.

Noise and quantum state decoherence: A quantum system remains in its state until it is disturbed by an external factor, which can be intentional, that is, measuring or resetting a qubit to state |0 , or unwanted, such as thermal noise or cosmic radiation. Since natural noises are all around us, a quantum computer has to be shielded. In the case of IBM Quantum TM , the quantum processors are cooled to almost 15 mK to avoid thermal noise. 44 There is also shielding against electromagnetic fields and so on. Despite these measures, some noise is inevitable in a quantum processor, because its temperature is always above 0 K and the radioactive background cannot be shielded out perfectly. What is more, quantum computers are probabilistic machines and errors can occur without any external force because of ubiquitous quantum fluctuations. As a result, a calculation carried out by a quantum processor is prone to errors. There are two basic types of error -bit-flip and phase-flip (referred to together as quantum state decoherence). The former refers to state |0 switching to |1 or vice versa, while the latter is connected with a change of the phase of a state, for example, |+ switching spontaneously to |− or vice versa.

The resilience of a quantum processor to flips is expressed by the relaxation time T 1 , i.e., the expected time a qubit remains in state |1 before switching (relaxing) to |0 , and the dephasing time T 2 , i.e., the expected time for which a qubit retains its phase. 45 Obviously, increasing times T 1 and T 2 leads to higher quality of quantum processors and allows more complex calculations to be performed. Currently, the relaxation and dephasing times are aproaching 1 millisecond, which is still relatively low for solving real-world tasks. However, at the beginning of the millennium, the decoherence times of IBM-like processors were in the tens of nanoseconds. 46 To demonstrate decoherence in practice, we can prepare two simple circuits. The first consists of an X gate applied to a qubit followed by several identity operators I. These operators do not change the quantum state but introduce a time delay between the preparation and measurement of the state. In an ideal world, the qubit would remain in state |1 indefinitely, but in the real world, an increasing number of I gates heightens the probability of spontaneous relaxation, as demonstrated in Figure 11 (left). The second circuit consists of an H gate changing state |0 to |+ , again several I gates, and finally an H gate (switching |+ back to |0 ). In an ideal world, we should measure state |0 with 100% probability regardless of the number of I gates. However, the higher the number of I gates, the higher the probability that |+ changes to |− and in the end we measure |1 instead of the expected |0 . This is illustrated in Figure 11 (right).

Note that bit flips and phase flips can be eliminated with the error correction schemes discussed in detail, for example, in [20] . According to the threshold theorem, with a sufficient number of auxiliary qubits, it is possible to achieve any accuracy of a quantum computation. This means that once we have quantum processors with thousands of qubits, we can implement corrections to eliminate decoherence effects. 47 The number of qubits is currently in the higher tens/lower hundreds (IBM Washington has 127 qubits and Google Bristlecone has 72). 44 Note that the temperature of interstellar space is 2.7 K (i.e., 2 700 mK). To reach milli-Kelvin temperatures, a special device called a dilution refrigerator is used. The refrigerator is based on spontaneous separation of the helium isotopes 3 He and 4 He into two layers at temperatures below 1 K (above this temperature, "classical" cryogenic techniques are used). The 4 He layer contains a small amount of 3 He. Much higher partial pressure of 3 He allows to pump 3 He out of the 4 He layer and forces 3 He from 3 He layer to diffuse into the 4 He layer. This diffusion is an endothermic process, hence it drains energy from the surrounding environment and makes it colder. Further information on dilution refrigerators and other cryogenic techniques can be found in [51] . 45 Note that spontaneous relaxation and dephasing are governed by the exponential law. The probability of an event occurring in a time interval of length t is P (t) = 1 − e −t/T , where T is either T 1 or T 2 . 46 See the conclusion of [52] for a short discussion of the development of decoherence times over the last two decades. 47 [53] showed in practice on a Google quantum processor that an increasing number of qubits employed in an error correction scheme leads to exponential attenuation of the error rate. However, protecting a quantum processor against errors requires a good deal more qubits than the protected calculation. Limited qubit connectivity: To be able to implement controlled gates, qubits have to communicate with each other. However, quantum processors, especially those based on superconducting qubits, 48 are planar devices, which prevents connections being established between all qubits. An example of this limitation in the case of the Yorktown IBM Quantum TM processor is shown in Figure 12 . If we want to establish a line of communication between two qubits that are not directly connected, we have to use intermediate qubits and swap gates. For example, in the case of the Yorktown processor, if we want to control the gate acting on qubit 4 with qubit 0, we have to first swap the states of qubits 0 and 2, then put the controlled gate on qubits 2 (control) and 4 (target), and again swap the states of qubits 2 and 0. Adding new swap gates to the circuit makes it more complex and hence more prone to noise. Together with decoherence, limited connectivity is hindering the application of quantum computers to real-world problems for the time being. However, a design with all-to-all connectivity has been tested and examined, for example, in [56] , and a similar design for a quantum processor with all-to-all connectivity is being developed, for example, by Honeywell (note, however, that these processors have worse decoherence times than those based on superconducting technology).

Quantum Random Access Memory (qRAM): Like any other computer, a quantum computer also needs a memory for storing data and programs. However, in contrast to classical computers, data and programs are stored directly in a quantum processor, because there is currently no functional quantum equivalent to an operational memory (quantum RAM or qRAM). This means that the quantum processor has to be reprogrammed when the input data are changed. Besides its obvious ability to store data, qRAM is important for a quantum computer for another three reasons. First, qRAM is able to prepare quantum states based on stored data. It translates a superposition of addresses n i=1 a i |i to a superposition of useful data n i=1 a i |i |D i , where |D i is a quantum state stored at address i. While the method for preparing a general quantum state proposed in [28] shows an exponential increase in the number of gates with an increasing number of qubits, qRAM is able to prepare the general state efficiently. Without qRAM, we can only prepare several types of quantum states with polynomially increasing resources in the number of qubits involved, for example:

• any basis state, such as |1 or |101 (we initialize all qubits to state |0 and put X gates on those we want to have in state |1 )

• W states and GHZ states, using the method proposed in [23] • states representing a discretized Gaussian distribution, using the algorithm designed in [57] Second, some algorithms, for example the HHL one, assume that the input data are stored in qRAM, otherwise the expected speed-up is not necessarily achieved. As mentioned in [31] , in the HHL algorithm, qRAM would help with the preparation of the right-side |b and even with the controlled gate e iA . Without qRAM, the HHL algorithm is unusable for real-world problems. 49 Third, in some cases, algorithms equipped with qRAM are faster than their original versions -as shown in [58] , a Grover database searching algorithm employing qRAM can achieve exponential speed-up instead of only quadratic. Currently, one of the biggest problems preventing qRAM from working is decoherence of stored quantum states. New designs have been devised in [59] and [60] , but these remain rather experimental. Fortunately, Qunnect has introduced photonic qRAM for quantum communication networks 50 and the company plans to commercialize this product [61] . Quantum volume: We have said that it is mainly the decoherence of quantum states which is responsible for the limited capabilities of current quantum processors. However, the relaxation time T 1 and the dephasing time T 2 are not the only important factors for assessing the quality of quantum processes. The overall quality also depends on the implementation of gates, on their speed and their ability to perform the expected operations correctly (gate fidelity), on the connectivity between qubits, and on other factors. To take all these factors into account, a metric called quantum volume has been introduced in [62] . In simple terms, the quantum volume is the maximum size n of a quantum circuit consisting of n qubits with n gates on each of them which can be implemented without adverse effects on the calculation carried out by the circuit. Note that the gates are native ones (i.e., implemented on the hardware level of the quantum processor). The quantum volume is actually expressed as the number of basis states the circuit can operate with. So, if the quantum volume is, say, 128, then n is 7, because 7 qubits can be in 128 basis states (128 = 2 7 ). For example, the quantum volume of the free-to-use IBM Quantum TM processor is between 8 and 32, which means that circuits consisting of 3 to 5 qubits with a maximum of 3 to 5 gates on each can be run without major errors. Note that in demonstrating decoherence, we used circuits with 200 identity gates. However, an identity gate is a simple wait instruction, and no other action which could further contribute to decoherence is performed.

Note, however, that quantum volume is a somewhat artificial benchmark, because it is based on the generation of random quantum circuits. A quantum processor can behave differently in practical applications such as HHL or QAOA algorithms. Moreover, quantum volume is based on square circuits (i.e., the number of qubits is the same as the number of gates applied to each of them), a condition hardly fulfilled in practice (most circuits are rectangular). As proposed in [63] , to fully assess the capabilities of a quantum processor, its performance in different practical and theoretical tasks has to be evaluated and benchmarking should not be reduced to only one number.

A closing remark on the state of development of quantum hardware: It is worth emphasizing that it is possible to overcome some of the above-mentioned obstacles by making changes to the design of the algorithms. For example, [64] introduced a hybrid version of the HHL algorithm exploiting the iterative phase estimation designed in [34] . This HHL algorithm enhancement significantly reduces the number of qubits used, hence the noise and the inter-qubit connectivity requirements are reduced as well. Although this approach is not a panacea, as it does not remove other issues connected with the design of the HHL algorithm (for example, the requirement for matrix sparsity), it may help demonstrate the capabilities of quantum computers, enable them to be used to tackle some real-world problems sooner, and facilitate the path to mature quantum computers.

In this part, we present the results of the calculation of risk metrics using the algorithm introduced in Section 2.2.1. In particular, we calculate the expected value, the standard deviation, historical Value-at-Risk (VaR), and historical CVaR (both VaR and CVaR are calculated for the 95% and 99% significance levels). The underlying data are the daily profits and losses (P/L) of a portfolio expressed in basis points.

For our demonstration, we chose the USD fixed investment portfolio (excluding mortgage-backed securities). 51 The time series begins on January 2 nd , 2018 52 and ends on October 4 th , 2021 and comprises 927 observations in total. A histogram of the daily P/L observations is depicted in Figure 13 . Note that P/L values below -65 bp and above 65 bp are considered to be outliers, as they occur only four and eight times, respectively. Preparing a histogram is the first step in calculating risk measures on a quantum computer, because the input distribution for the calculation has to be discrete (the histogram is the discretized version of the continuous distribution). In the next step, an integer value (starting with 0) is assigned to each histogram bin, the integers are then expressed as binary numbers, and the numbers of occurrences are replaced by relative counts. In other words, the histogram is converted into an empirical probability function of a discrete probability distribution. We carried out the discretization for histograms with 8 and 16 bins, as these numbers are powers of two. In Table 2 , we show the expected values of the risk measures as calculated classically (in MS Excel), first with the assumption that the underlying distribution is continuous and second after discretization. Once the distribution is discretized, we get the value of each 51 The Czech National Bank reserves consist of two tranches -a liquidity one and an investment one. The former is used for monetary policy implementation and is composed of highly liquid short-term instruments denominated in EUR and USD only, while the latter is considered to be a gain generator. The investment tranche has two main parts -fixed-income (AUD, CAD, CNY, EUR, GBP, SEK, and USD) and equity (AUD, CAD, EUR, GBP, JPY, and USD). The investment tranche also comprises a small portion of USD-denominated mortgage-backed securities (MBS) and physical gold. The structure of the reserves is depicted in a diagram in Appendix A. 52 The tranching was introduced in August 2017, hence 2018 is the first full year it was active.

risk metric expressed as the number of the histogram bin. 53 Based on this number, and knowing the lower and upper bounds of the bins, we get an estimate of the value of the metric. As expected, the numbers in Table 2 show that finer discretization returns more accurate estimates of the values of the metrics. Note that in the case of the standard deviation, the bin number means "the size of the deviation expressed in bins", that is, to get the actual standard deviation, the bin number is multiplied by the bin width. The approach described above enables us to encode the discretized distribution into a quantum state |ψ = N −1 i=0 √ p i |i , where p i is the relative count in the i th bin and |i is a basis state representing the i th bin. This state is an input to the quantum calculation of the risk metrics. Note that |ψ is a three-qubit state for eight bins (8 = 2 3 ) and a four-qubit state for 16 bins (16 = 2 4 ). One additional qubit is necessary for the calculation, hence four and five qubits are used, respectively. First, we ran the calculation on the IBM Quantum TM simulator to check the correctness of the algorithm design. Note that the simulator behaves as an error-free quantum computer with all-to-all qubit connectivity. As can be seen in Table 3 , the results are in accordance with the classical calculation presented in Table 2 . Small deviations are caused by rounding and sampling errors in the simulation.

After checking on the simulator, we switched to actual quantum hardware. We employed four IBM Quantum TM processors available for free, namely, Yorktown, Lima, Belem, and Bogota, all of them five-qubit devices. 54 While Yorktown is a representative of an older generation of quantum processors (family code name Canary) and was retired shortly after our experiments were carried out, others are more modern IBM processors (family code name Falcon). 55 The quantum volume is 8 in the case of Yorktown and Lima, 16 for Belem, and 32 for Bogota. The decoherence times are close to 100 µs for the Falcon family processors and around 50 µs for the Yorktown processor. The processors also differ in qubit connectivity. While Bogota is a single row of qubits, Lima and Belem have their qubits arranged in a T-shaped lattice, and a schematic of the Yorktown connections is shown in Figure 12 .

The results of the calculations on the above-mentioned devices are shown in Table 3 . As can be seen, the expected value was calculated almost correctly for both the three-and four-qubit states, with the exception of the Lima processor. Unfortunately, the evaluation of the standard deviation and the 95% VaR is totally incorrect. As a result, CVaR is also evaluated wrongly. Although the 99% VaR for the three-qubit case seems to be correct, we are skeptical about this result, as it looks like the VaR is nested in the value 1. The most probable cause of this issue is the minimum possible change in the rotation angle of the Ry gates. As we approach lower percentiles, the rotational angles get smaller and the difference in the angles between consecutive iterations of the algorithm is below the lowest distinguishable level. This makes it impossible to distinguish between low percentiles of the distribution. Interestingly, if the 99% VaR was really equal to 1, the 99% CVaR calculation is not far from the correct value for the three-qubit case on the Lima and Belem processors. 53 For example, the expected value is calculated as N −1 i=0 ip i , where i is the number of the i th bin and p i is the probability that the P/L value belongs to the i th bin. 54 Note that other processors -Santiago, Manila, and Quito -are also available for free, but their characteristics are similar to the named ones so they were not used in this experiment. IBM also provides the Armonk processor, but it is a one-qubit device and is therefore unsuitable for our experiments employing more qubits. 55 The state-of-the-art IBM processor (not available for free) is the 127-qubit Washington (family code name Eagle). Based on these results, it seems meaningless to employ quantum computers to evaluate risk metrics, but, as discussed in Section 3.1, quantum computers are relatively new devices and some further development is necessary. However, we can extract several useful conclusions from the results:

• We provided a proof-of-concept that the quantum algorithm for calculating risk measures works.

We showed this on the simulator for all metrics and on real quantum processors for the expected value. This is the most important conclusion, because it shows that quantum computers have the potential to be employed in FX reserves management (or portfolio management in general).

• Qubit connectivity plays an important role in error mitigation. While the Falcon family processors have better decoherence times, the results obtained are often worse in comparison with the older Yorktown processor (for example, the expected value for the four-qubit state). As can be seen in Figure 12 , the number of connections is six in the case of Yorktown and only four in the case of the T-shaped (Belem and Lima) and row (Bogota) lattices. A lower number of connections leads to an increase in the number of swap gates, which makes the circuits more complex and hence prone to decoherence.

• The quantum volume of a processor has to be used carefully. Although the quantum volumes of the processors we employed differ significantly, the results are similar. One would expect Bogota, which has the highest quantum volume, to be the best processor, but this is not the case. The reason is that quantum volume assumes square circuits, whereas our circuit widths (the maximum number of gates on one qubit 56 ) are far higher than our circuit heights (the number of qubits). For example, the circuit for the expected value calculation (the three-qubit case) on the Lima processor has a width of 69 but a height of only 4.

We provide the complete source code for the risk metric calculation in Appendix F. Note that we implemented the preparation of the quantum state |ψ based on the method described in [28] . We omitted the part for setting the quantum phase, as this is not necessary for the risk metric calculation. One could also employ the initialize function available in the Qiskit libraries. However, for educational purposes, we wrote the quantum state preparation routine ourselves.

In this section, we present portfolio optimization based on solving a system of linear equations with the HHL algorithm. In particular, we demonstrate the balancing of the ratio between fixed-income assets and equities for a required expected return. For the demonstration, we used our USD fixed-income investment portfolio (excluding MBS), which consists of investment-grade bonds (issued by governments, agencies, and supranational bodies) and plain-vanilla interest rate swaps, and our USD equity portfolio, which replicates the S&P 500 index. First, we calculated the one-year moving returns for each day of the period from January 2 nd , 2019 to August 4 th , 2021 for both portfolios. For illustration, we present the returns on both portfolios in Figure 14 . The graph indicates very well the negative correlation between bond and equity returns. The market risk-off regime introduced at the beginning of the Covid-19 pandemic (March 2020) and the subsequent economic recovery in spring 2021 are clearly visible as well.

Based on the moving returns, we calculated the average yearly returns for fixed-income securities (5.86%) and equities (16.78%). These figures serve as the expected returns forming vector |R . The covariance matrix, again based on the moving returns, is

Note that 0.15% is the variance of the returns on the fixed-income portfolio and 2.46% is that of the returns on the equity portfolio. We set the price vector P | = (1, 1) and the budget B = 1 to impose the constraint w fix. inc. + w equity = 1. The required return G is a free parameter of the optimization task and is set by the user. Substituting all these figures into (46), we get a linear system describing the optimization task: 57     0.0000 0.0000 0.0586 0.1678 0.0000 0.0000 1.0000 1.0000 0.0586 1.0000 0.0015 −0.0043 0.1678 1.0000 −0.0043 0.0246

As can be seen in equation (54), the linear system breaks up into two separate systems in the case of only two assets. The first is composed of the equations 0.0586w fix. inc + 0.1678w equity = G and w fix. inc + w equity = 1. This allows us to evaluate the weights for a given expected return G regardless of the rest of (54), i.e., the risk part. For example, for G = 7% = 0.07, we get w fix. inc = 89.53% and w equity = 10.47%. Although this task is trivial and can be solved quickly by hand, its purpose is to demonstrate the capabilities and identify the bottlenecks of quantum computers and the HHL algorithm in portfolio optimization tasks.

As in the risk measurement case present in the previous section, we first run the algorithm on a quantum computer simulator. Note that despite the loss of the speed-up delivered by a quantum computer, we decided to extract the whole solution of (54) to assess the correctness of the implementation of the algorithm. Unfortunately, when we ran the HHL solver for G = 7%, we obtained totally wrong results -w fix. inc = 10.29% and w equity = 18.41%. To debug the code, we ran the HHL solver on a simple system with matrix A = diag (1, 2, 3, 4) and right-side b| = (1, 1, 1, 1) having the solution x| = (1, 1/2, 1/3, 1/4). In this case, the solver returned the correct result. Interestingly, for A = diag(−1, 2, 3, 4) and b| = (−1, 1, 1, 1) with the same solution x|, the HHL solver failed again, as the returned x| = (0.1429, 1/2, 1/3, 1/4). After several attempts, we realized that the correctness of the solution is linked with the matrix eigenvalues. In the case of a diagonal matrix, the eigenvalues are the diagonal elements. While in the former case they are all positive, in the latter there is one negative eigenvalue. To examine this conjecture further, we prepared matrix (H⊗H)diag(1, 2, 3, 4)(H⊗H), where H is a Hadamard gate. The resulting matrix is derived from the diagonal one by unitary transformation H ⊗ H, hence its eigenvalues are the same as those of matrix diag (1, 2, 3, 4) . Setting b| = (0, 1, 1, 0), we get the correct solution x| = (3/8, 5/8, 5/8, 3/8). However, changing the first element of the diagonal matrix to -1 led again to incorrect results. After some additional investigation, we realized that the current Qiskit implementation of the HHL algorithm is only able to work with matrices where all the eigenvalues are positive. 58 We calculated the eigenvalues of the matrix in system (54) and found that they are both positive and negative, specifically: -16.90792562, -0.36203117, 1.03367322, 18.84628357. As a result, we were unable to carry out the portfolio optimization task with the HHL algorithm successfully for the time being.

For illustration purposes, we tried to solve a simple system with the diagonal matrix diag(1, 2, 3, 4) on a real quantum processor. However, the results suffered from decoherence and were meaningless. Note that the circuit width was 159 gates. In the case of the portfolio optimization task (54), the width even reached 13,762 gates. Hence, regardless of the issue with the eigenvalues, the task would still be unsolvable on the current quantum hardware, because the circuit is too deep and complex.

To assess the current capabilities of free-to-use quantum processors in solving linear systems, we tried to run a circuit for the linear system 1.5 0.5 0.5 1.5

discussed in [32] and [14] . The implementation of the circuit is shown in Figure 15 . Note that the right side of (55) is the normalized quantum state and can be easily prepared with an Ry gate applied to state |0 . We run the circuit for several angles θ: π/3, π/4, π/6, and π/7. Regardless of the IBM quantum processor used (Lima, Quito, and Bogota were employed), the results were fully obscured because of decoherence. The circuit width after adaptation to the connectivity of the real processor and its basic gate set reached 102 gates for all the processors tested. Based on the results presented above, we conclude that the quantum hardware we tested is not able to perform the HHL algorithm successfully. The first issue is decoherence of the quantum states, caused by deep circuits and short decoherence times. The second issue is the missing implementation of the HHL algorithm for matrices with negative eigenvalues. The latter issue will be overcome soon, once IBM developers expand the capabilities of the Qiskit implementation of the HHL algorithm. However, the former issue could take more time to resolve, as it is connected with further research in the field of quantum processor construction. It is worth noting that we only employed quantum processors that are available to use free of charge. However, the decoherence times of all the IBM devices are close to 100 µs. 59 For educational purposes and hopefully future use, we provide the source code for portfolio optimization using the HHL algorithm and the circuit for solving system (55) in Appendix F.

This part is devoted to picking a subset of assets from a larger set, for example, a benchmark index, subject to some constraints and objectives. To do so, we employ quadratic unconstrained binary optimization and the QAOA algorithm. In our case, we want to choose three shares from a set containing five shares in total, and at the same time to maximize the return and minimize the market risk. The number of shares in the larger set is limited by the number of qubits on the quantum processors available to us (i.e., five). The shares we deal with in the demonstration are five semiconductor industry companies contained in the S&P 500 index, namely, AMD, Intel, Qualcomm, Analog Devices, and Texas Instruments. 60 Figure 15 : HHL Algorithm Circuit for a 2x2 System of Linear Equations Note: The circuit implements the HHL algorithm for system (55) with θ = π/3. The first part of the circuit implements the preparation of the right-side |b with the Ry(2θ) gate. The second part represents controlled rotations of the phase estimation algorithm. The algorithm continues with inverse QFT in the third part, and the fourth part implements the inverse of the eigenvalues. The last part provides the uncomputation of the working qubits q 0 and q 1 . Note that q 3 is the ancilla qubit and the solution is stored in qubit q 2 . Source: Author's own creation in Qiskit based on [32] and [14] To perform the task described above, we modified objective function (52) by setting the budget equal to the number of shares we want to invest in, i.e., B = m, and the maximum amount of money invested in the i th share was set to 1. Moreover, we need to rewrite (52) to the general form of a quadratic binary optimization objective function, i.e., f (x) = x|A|x + b|x + c. After algebraic manipulation, we end up with the objective function

where r| is a vector of share returns, C is a covariance matrix, m is the number of shares we want to invest in, 1 n,n is an n-by-n matrix having all elements equal to 1, and 1 n | is a row vector with all n members equal to 1. The returns vector and the covariance matrix are based on one-year moving returns calculated similarly to those in Section 3.3. The results of this calculation are shown in Table 4 . Parameters λ describing the importance and scale of each part of the objective function were set as follows: λ 1 = λ 3 = 1 and λ 2 = 4 (the risk term is scaled by λ 2 = 4 to be similar in magnitude to the return term).

To compare the results obtained on the quantum computer, we first calculated the minimum value of function (56) classically using the brute force method, i.e., all the solutions were listed and the objective function was evaluated. Since we are dealing with five equities, there are 2 5 = 32 possible solutions, each denoted with a bit string representing the equities included in the portfolio (the i th bit is equal to one 60 There is no particular reason for selecting these companies. In fact, the Czech National Bank's equity investments follow a passive index replication strategy and we do not choose particular shares to invest in. Therefore, this example serves only to demonstrate the capabilities of the QAOA on real market data and is not connected with the actual investment decision-making at the Czech National Bank.

if the i th share is included). However, only solutions having exactly three bits equal to one are feasible solutions, as we require m = 3. The optimal solution is x opt | = (1, 1, 0, 1, 0) , i.e., we should invest in AMD, Intel, and Analog Devices. All the possible solutions are presented in Table E1 in the Appendix. We carried out the calculation on all the available five-qubit real quantum processors with success. The only difference among the processors was the number of iterations needed to find the solution, as can be seen in Table 5 . The number of iterations depends on the quantum volume of the processor and is lowest for processors with quantum volume 32. It is worth noting that the number of iterations is higher than the number of possible solutions of the task (i.e., 32). This means that the QAOA seems to be less effective than the brute force method. However, we have to take into account the fact that the current quantum computers still suffer from noise and, as a result, more iterations are needed. Despite this issue, we found that the QAOA can be performed successfully on the current real quantum hardware. Interestingly, a quantum circuit for a five-qubit QAOA (see Appendix E) is comparable in terms of the number of gates to circuits performing the other algorithms discussed in this paper. However, in contrast to those, we recorded 100% success with the QAOA. This shows the importance of application-oriented benchmarking of quantum computers. In other words, we should assess the quality of the current quantum processors based on a particular application and not only on a single number like the quantum volume, even though it is also an important indicator.

In this article, we investigated possible applications of quantum computers in foreign exchange reserves management. We discussed risk measurement and portfolio construction. For the former task, we employed the quantum Monte Carlo method, while for the latter we used the HHL algorithm and quadratic unconstrained binary optimization with the QAOA. We also discussed technical issues connected with the currently available quantum computers. All tests and demonstrations were carried out on IBM Quantum TM processors. Our calculations were based on real data taken from the Czech National Bank's daily internal reporting on its foreign exchange reserves. Besides analyzing the capabilities of the current quantum computers, we also provided a detailed introduction to the mathematical background and basic building blocks of quantum computing, explained how several quantum algorithms work, and supplemented our discussion with examples. Moreover, we shared all the Qiskit source codes we used to help the reader with studying the basics of quantum computing on her/his own.

We found that, in principle, all the algorithms discussed are suitable for the purposes for which we employed them. First, we ran the algorithms on a quantum computer simulator and all the results obtained were in accordance with those obtained using classical methods.

After the simulations, we performed the tasks on several five-qubit real quantum processors. In this case, the results were mixed. Quadratic binary optimization for five variables using the QAOA returned the correct results on all the available processors. This means that even on the current noisy quantum computers, the QAOA ran successfully. Based on this result and several other studies mentioned in the theoretical part, it seems that binary optimization will be among the first real-world problems that quantum computers are used for. However, it is worth noting that the number of iterations needed to find the solution depended on the quantum volume of the quantum processor. In other words, noise still plays some role in algorithm efficiency.

On the other hand, portfolio optimization leveraging the HHL algorithm failed completely. This is because the current implementation of the HHL algorithm in Qiskit cannot work with matrices that have negative eigenvalues. Therefore, we tried to solve a simple linear system with a 2x2 matrix with all the eigenvalues positive. However, despite the fact only four qubits were used, wrong results were returned for all the right sides we used. The most probable reason for this is the complexity of the quantum circuit representing the HHL algorithm, especially the phase estimation part, which consists of many controlled quantum gates.

The quantum risk measurement algorithm returned the correct results for the expected value but failed for the other risk parameters evaluated. Note that we tried to calculate risk measures with two different levels of input probability distribution quantization -8 bins (three qubits) and 16 bins (four qubits). The performance of the algorithm was better in the case of three qubits because of a less complex circuit, which meant it was less impacted by noise. There is also an issue with limited accuracy of the rotational angles of the Ry gates, which prevented us from making a finer distinction between the percentiles in the Value-at-Risk calculation.

Interestingly, despite the fact we used "free" quantum processors, we had access to state-of-the-art technology, because even IBM processors with higher numbers of qubits show similar decoherence times and gate error rates to smaller five-qubit devices. For example, the quantum volume of the 127-qubit Washington processor is 32, i.e., the same as for the five-qubit Manila processor we used in our simulations. If we had access to a bigger processor, we would be able to implement, for example, binary optimization with more variables. However, noise would still negatively impact the results. Overall, this means that currently the biggest limiting factor of quantum processors is noise, not the number of qubits. Although the presented results and discussion could suggest that quantum computers are of low practical value, we should keep in mind that they represent a very new technology, and intensive research work continues to enhance the parameters of quantum processors. Recent advances in the quality of quantum processors allowed us to perform at least some calculations correctly. Had we carried out similar simulations a few years ago, we would have failed in all cases.

Let us briefly discuss future applications of quantum computers in the finance industry and their impact on central banks and other financial market regulatory authorities. The possible fields of application of quantum computers range from complex derivatives pricing to credit risk analysis (including loan underwriting). In these cases, the quantum Monte Carlo method can be deployed. Because of their promised high computational performance, quantum computers would be interesting to high-frequency traders. As such investors bother about every microsecond, they would again employ the quantum Monte Carlo or quick portfolio optimization methods (i.e., the HHL algorithm or the QAOA). Clearly, banks and other financial companies have incentives to exploit quantum computing. All the topics discussed in this article are therefore important not only for the FX reserves management arms of central banks, but also for their supervision and regulation departments and other supervisory authorities, because the latter have to examine financial firms' pricing, credit, and market models closely to prevent misconduct, both intentional and unintentional, and ensure the stability of financial systems.

It is true that many technological steps to advance quantum computers have to take place before the deployments discussed above can materialize. For example, decoherence times and the number of qubits have to be increased, gate fidelity has to be improved, and, above all, a functional qRAM has to be built.

Despite the rapid development, there is still a long way to a fully mature quantum computer. However, central banks and financial market regulators in general should be prepared for the massive deployment of quantum technologies, as it is not a question of if but when this happens. 

In this appendix, we present a list of quantum gates. This list is not exhaustive, but it contains the usually used gates and their widely accepted symbols. Note that in the case of controlled gates, the part of the matrix corresponding to the non-controlled version of the gate is highlighted in red.

• Pauli gates

Note that for any Pauli gate A, it holds that A † = A, hence A 2 = I.

• Square root of Pauli X gate √ X = 1 2

Note that the square root of X can be obtained with a spectral decomposition. Matrix X has eigenvalues 1 and -1 and respective eigenvectors |+ and |− .

• H, S, and T gates

Note that H † = H, hence H 2 = I.

• Rotations around axes of Bloch sphere R x (α) = cos(α/2) −i sin(α/2) −i sin(α/2) cos(α/2) R y (α) = cos(α/2) − sin(α/2) sin(α/2) cos(α/2)

Note that:

-R † a (α) = R a (−α) for any rotation around axis a ∈ {x, y, z} -R a (α)R a (β) = R a (α + β) for any rotation around axis a ∈ {x, y, z} Here we provide examples of universal gate sets, i.e., groups of gates allowing any unitary operation to be implemented on a quantum computer:

• H, T, and CNOT (sometimes enriched by gates X and S for easier implementation)

• R x (θ), R y (θ), R z (θ), and CNOT 

In this appendix, we describe the binary number encoding called Gray code. In comparison with standard binary number representation, the code has the special property that the immediate successor of a binary number differs from that number only in one bit. This property for four-bit numbers can be easily seen in Table C1 . Before we describe how to convert a binary number to Gray code and vice versa, we define an XOR function 61 denoted by symbol ⊕. For two input bits a and b, the XOR returns 1 if a = b and 0 if a = b, i.e., 0 ⊕ 0 = 0, 0 ⊕ 1 = 1, 1 ⊕ 0 = 1, and 1 ⊕ 1 = 0.

Let b be an n-bit binary number and b i its i th bit. Assume that b n is the leftmost bit. Then the representation of b in Gray code is given by the equations

In practice, the conversion can be carried out using the following procedure based on the equations above: 62

• take number b and the same number shifted by one bit to the right • remove the rightmost bit of the shifted number

• apply the XOR function to the equivalent bit positions in both numbers (bit-wise XOR)

We illustrate the procedure on the conversion of 1011 (i.e., the decimal number 11) to Gray code. Shifting the number to the right, we get 0101(1). The rightmost bit of the shifted number (i.e., 1) is forgotten, hence we are left with 0101. The next step is to apply the bit-wise XOR function to numbers 1011 and 0101, which delivers 1110, because the bits on the first (leftmost), second, and third places differ from each other, while the last bits (the rightmost ones) do not.

The conversion from Gray code back to standard binary representation is done as follows:

. . . b 2 = g n ⊕ g n−1 ⊕ · · · ⊕ g 2 b 1 = g n ⊕ g n−1 ⊕ · · · ⊕ g 2 ⊕ g 1 61 Sometimes denoted exclusive-OR or non-equivalence. 62 We employ this procedure in our risk metrics calculation. A Python function converting a binary number into Gray code is available in Appendix F.

IBM Quantum TM is a cloud-based quantum computing platform. It allows users to use real quantum processors and quantum computer simulators (i.e., software running on a classical computer and simulating a quantum computer). Some of the processors are available for free for education and research, provided that the results obtained are not used for commercial purposes. Currently, the fee-free processors comprise one single-qubit machine and six five-qubit processors. A 127-qubit processor is available on a partnership basis, and IBM plans to build a 1,000-qubit processor by 2023.

There are several ways users can program the quantum processors and simulators. The first is by using a quantum composer. This feature allows the user to build a quantum circuit using the drag-and-drop method. The composer also supports the QASM assembly language, which is a text description of the gates used in the circuits. The composer mode also provides visualization of the resulting quantum state and the probability distribution of the circuit's output. The composer is useful for the beginners in quantum computing and for circuit debugging.

A more sophisticated tool for the development of quantum algorithms is the Python-based language Qiskit. IBM provides a Quantum Lab, which is web-based environment using interactive Jupyter notebooks. Users can write Python code here, debug it, and run it on the quantum processors. All the classical Python libraries are also available to be combined with the Qiskit ones. IBM provides several modules for implementing the most important quantum algorithms, along with specialized libraries for finance, machine learning, optimization, and quantum chemistry. Qiskit can be installed on a local computer and included in Python development tools, too.

In this appendix, we provide a list of all the possible portfolios which can occur in the demonstration of binary optimization using the QAOA presented in Section 3.4. We also show a quantum circuit implementing one iteration of the QAOA for five qubits. 

# -------import libraries -----from qiskit import ClassicalRegister , QuantumRegister , Q uantumCi rcuit from qiskit import Aer , execute , IBMQ from qiskit . compiler import transpile , assemble import numpy as np import math as m 

Quantum mechanical computers

Rapid solutions of problems by quantum computation

Algorithms for quantum computation: Discrete logarithms and factoring

Experimental quantum cryptography

An integrated space-to-ground quantum communication network over 4,600 kilometers

A fast quantum mechanical algorithm for database search

A quantum algorithm for finding the minimum

Quantum algorithm for linear systems of equations

Quantum algorithm implementations for beginners

Quantum computational finance: Quantum algorithm for portfolio optimization

Quantum risk analysis

Credit risk analysis using quantum computers

Towards pricing financial derivatives with an ibm quantum computer

A survey on hhl algorithm: From theory to application in quantum machine learning

Quantum computing for finance: Overview and prospects

Quantum computing for finance: State of the art and future prospects

Prospects and challenges of quantum finance

Quantum technology for economists

Introduction to coding quantum algorithms: A tutorial series using qiskit

Quantum Computation and Quantum Information

Feynman Lectures on Computation

Particles, Fields and Forces: A Conceptual Guide to Quantum Field Theory and the Standard Model

Efficient quantum algorithms for ghz and w states, and implementation on the ibm quantum computer

The physical implementation of quantum computation

Elementary gates for quantum computation

The solovay-kitaev algorithm

Quantum fingerprinting

Transformation of quantum states using uniformly controlled rotations

Quantum versus classical generative modeling in finance

Equation solving by simulation

Quantum machine learning algorithms: Read the fine print

Quantum circuit design for solving linear systems of equations

Step-by-step hhl algorithm walkthrough to enhance the understanding of critical quantum computing concepts

Arbitrary accuracy iterative quantum phase estimation algorithm using a single ancillary qubit: A two-qubit benchmark

Matrix multiplication via arithmetic progressions

Gaussian elimination is not optimal

Ising formulations of many np problems

A quantum approximate optimization algorithm

Financial portfolio management using d-wave quantum optimizer: The case of abu dhabi securities exchange

Improving variational quantum optimization using cvar

Grover adaptive search for constrained polynomial binary optimization

Getting started with the d-wave system

Quantum computing in the nisq era and beyond

Towards prediction of financial crashes with a d-wave quantum computer

Formulating and solving routing problems on quantum computers

Analytical solution for low energy state estimation by quantum annealing to arbitrary ising spin hamiltonian

The pitfalls of planar spin-glass benchmarks: Raising the bar for quantum annealers (again)

Equation of state calculation by fast computing machines

Optimization Algorithms in Physics. Wiley-VCH

Quantum portfolio optimization with investment bands and target volatility

Basic operation of cryocoolers and related thermal machines

Moving beyond the transmon: Noise-protected superconducting quantum circuits

Superconducting qubits: Current state of play

Realization of efficient quantum gates with a superconducting qubit-qutrit circuit

Benchmarking an 11-qubit quantum computer

The efficient preparation of normal distributions in quantum registers

The effective solving of the tasks from np by a quantum computer

Quantum random access memory

Optically addressable nuclear spins in a solid with a six-hour coherence time

Qunnect announces sale of first commercial quantum memory

Validating quantum computers using randomized model circuits

A volumetric framework for quantum computer benchmarks

Advanced quantum supremacy using a hybrid algorithm for linear systems of equations

Qiskit documentation

F.4 Solving 2x2 Linear System with HHL from qiskit import ClassicalRegister , QuantumRegister , Q uantumCi rcuit from qiskit import Aer , execute , IBMQ from qiskit . compiler import transpile , assemble import numpy as np from numpy import linalg as la import math as m # m e a s u r e m e n t n = len ( returns ) A = lambda_risk * c o v a r i a n c e _ m a t r i x + l a m b d a _ a s s e t _ n u m b e r * np . ones (( n , n ) ) b = -( lambd a_return s * returns + l a m b d a _ a s s e t _ n u m b e r * 2 * assets_number * np . ones (( 1 , n ) ) ) # after c a l c u l a t i o n before a matrix 1 -by -n is returned , # however , the solver expects a