Experimental implementation of quantum-walk-based portfolio optimization

Quantum computation promises advantages for certain problems, such as prime factoring, exhaustive search and optimization [1–9]. Although only small quantum processors are currently available, there are tremendous expectations for quantum technology to assist with the solving of difficult real-life problems in the near future [10–12]. Modern financial trading relies on extensive computational resources for the analysis of historical data, high-frequency trading, portfolio optimization and risk analysis [13–16]. In the case of portfolio optimization, the computational complexity stems from the analysis of the balance between return and risk of portfolios involving many assets, motivating quantum computing to solve financial problems [17–23].

Recent works have combined financial services and quantum software to implement and test a financial portfolio rebalancing use case using hybrid quantum–classical algorithms [24]. Based on the mean-variance Markowitz model [25] with discrete asset constraints, Hodson et al explored the problem of portfolio optimization and rebalancing with the help of the quantum approximate optimization algorithm and quantum alternating operator ansatz, collectively known as QAOA [26–36].

Such heuristic, or approximation, algorithms are designed to identify high-quality solutions to optimization problems [37–40]. The problem solution space is mapped to a superposition of quantum states that is acted on by an interleaved sequence of solution-quality-dependent phase shits and mixing operations parameterised by classically tunable variational parameters. These are tuned using a classical algorithm that seeks to optimise the mean of the measured solution qualities. In doing so, constructive interference is produced at states corresponding to high-quality solutions.

Extending the QAOA schema to constrained optimization problems requires strategies for restricting or guiding probability amplitude to a subset of the problem solution space. The Quantum Walk-based Optimization Algorithm (QWOA) achieves this using an efficient indexing algorithm for the valid solutions and generalizing the QAOA mixing operation to a Continuous-Time Quantum Walk (CTQW) [41]. Slate et al proposed an application of the QWOA to the portfolio optimization problem, with numerical results indicating a significant advantage over previously explored methods [42]. Key to this was the ability of the QWOA to restrict the quantum search and, consequently, variational parameter optimization to a non-degenerate subspace of valid solutions.

This paper investigates the efficacy of the QWOA mixing operator as applied to the portfolio optimization problem. Using linear optical elements and single photons, we realize the first experimental implementation of the QWOA mixing unitary and demonstrate its reliable convergence to high-quality solutions over a wide range of quantum circuit depths. Our work provides strong evidence for the potential of quantum-walk-based algorithms to solve complex optimization problems of practical significance.

Combinatorial optimization problems seek solutions of the form

$$\begin{equation} \bar{s} = \left\{s \; | \; C\left(s\right) \in \text{min} \left\{C\left(s\right) \; | \; s \in \boldsymbol{\mathcal{S}}^\prime{}\right\}\right\}, \end{equation} \tag{ 1 }$$

where the problem cost-function C(s) is defined to map a solution s from an ordered set of problem solutions $\boldsymbol{\mathcal{S}}^{\prime} = \{s_i\}$ to $\mathbb{R}$, s is a k-permutation of discrete elements selected from a finite set and,

$$\begin{equation} \boldsymbol{\mathcal{S}}^{\prime} = \left\{s \; | \; s \in \boldsymbol{\chi}\right\} \end{equation} \tag{ 2 }$$

is the valid problem solution space where

$$\begin{equation} \boldsymbol{\chi} = \bigcup_i \left\{s \; | \; \chi_i\left(s\right) = a_i\right\} \end{equation} \tag{ 3 }$$

represents constraints imposed on $\bar{s}$, while $\boldsymbol{a} = \{a_i\}$ defines the constraints.

The QWOA [41, 43] begins by establishing an injective map from the complete solution space $\boldsymbol{\mathcal{S}}$ to a Hilbert space consisting of n qubits $\mathcal{H}$ where $2^n \unicode{x2A7E} |\boldsymbol{\mathcal{S}}|$ in order to define the quality operator

$$\begin{equation} Q = \sum_{i = 0}^{\left| \boldsymbol{\mathcal{S}} \right| - 1}q_i \vert i\rangle\langle i\vert, \end{equation} \tag{ 4 }$$

where $q_i = C(s_i)$.

We next require that there exists an efficient indexing algorithm from the complete solution space $\boldsymbol{\mathcal{S}}$ to the subset of valid solutions $\boldsymbol{\mathcal{S}}^{\prime}$, which makes possible an indexing unitary

$$\begin{equation} U_{\#}^{\dagger}|i\rangle = \left\{\begin{array}{l} \left|\mathrm{id}_{\boldsymbol{\chi}}\left(i\right)\right\rangle,|i\rangle \in\left|\boldsymbol{S}^{^{\prime}}\right\rangle, \\ |i\rangle, \text {otherwise, } \end{array}\right. \end{equation} \tag{ 5 }$$

where $U^{\dagger}_\#$ maps the states corresponding to valid solutions $\vert s^{\prime}\rangle$ to the indexed states $\vert \text{id}_{\boldsymbol{\chi}}(i)\rangle$. This enables the preparation of an initial state $\vert \psi_0\rangle$ that is equally superposed over qubits corresponding only to the valid solution space $\boldsymbol{\mathcal{S}}^{\prime}$,

$$\begin{equation} \vert \psi_0\rangle = \frac{1}{\sqrt{\left|\boldsymbol{\mathcal{S}}^{\prime}\right|}}\sum\limits_{i \in \boldsymbol{\mathcal{S}}^{\prime}}\vert i\rangle. \end{equation} \tag{ 6 }$$

The quality operator Q and the indexing unitary $U_\#^\dagger$ are used as part of a variational quantum circuit that works to amplify probability density associated with high-quality (low-cost) solutions. The variational quantum circuit implements a sequence of alternating operators. The first of these is the 'phase-shift' unitary

$$\begin{equation} U_Q\left(\gamma\right) = \text{exp}\left(- \text{i} \gamma Q\right), \end{equation} \tag{ 7 }$$

where $\gamma \in \mathbb{R}$, applies a quality-dependent phase shift to all quantum states $\vert i\rangle$.

Next is the indexed-walk unitary,

$$\begin{equation} U_w\left(t\right) = U_\# \text{exp}\left(-\text{i} t W\right) U_\#^\dagger, \end{equation} \tag{ 8 }$$

where $t \unicode{x2A7E} 0$ and $W = \sum_{i, j \in \boldsymbol{\mathcal{S}}^{\prime}} \vert \text{id}_{\boldsymbol{\chi}}(j)\rangle\langle \text{id}_{\boldsymbol{\chi}}(i)\vert$ for i ≠ j is equivalent to the adjacency matrix of a complete graph. The indexed-walk unitary U_w indexes $\vert \psi_0\rangle$ and induces a maximal unbiased coupling over the valid solution subspace which drives the mixing of probability amplitude between states in the subspace. Equivalent to a CTQW for time t, during this stage phase differences encoded by U_Q contribute to constructive and destructive interference. After competition of the CTQW, the indexing unitary un-indexes $U_\#^\dagger \vert i\rangle$ such that measurement is conducted with respect to the unindexed states $\vert i\rangle$.

An application of U_Q and U_w is carried out for p repetitions, with the final state of the quantum system being

$$\begin{equation} \vert \boldsymbol{\gamma}, \boldsymbol{t}\rangle = \prod_{i = 1}^{p} U_w\left(t_i\right)U_Q\left(\gamma_i\right)\vert \psi_0\rangle. \end{equation} \tag{ 9 }$$

With each of the p iterations, the potential for probably amplitude concentration at $\vert i\rangle$ corresponding to high-quality solutions increases—given suitable values for the variational parameters γ and t [44, 45]. Optimal variational parameters γ and t are obtained using a classical optimiser that minimises with respect to the objective function

$$\begin{equation} f\left(\boldsymbol{\gamma}, \boldsymbol{t}\right) = \langle \boldsymbol{\gamma}, \boldsymbol{t}\vert Q\vert \boldsymbol{\gamma},\boldsymbol{t}\rangle, \end{equation} \tag{ 10 }$$

such that a lowering of $f(\boldsymbol{\gamma}, \boldsymbol{t})$ entails an increased probability of measuring a quantum state that corresponds to a high-quality solution. In this sense, the QWOA is a hybrid algorithm whereby a small set of classical values parametrise the evolution of an exponentially larger quantum Hilbert space.

Consider an investor with a portfolio of assets $\{m_i\}$ of size M. For each asset, they must take one of three positions:

• (i)  
  Short position: the buying and short-term selling of an asset under the expectation that it will drop in value.
• (ii)  
  Long position: the buying and long-term holding of an asset with the expectation that it will increase in value.
• (iii)  
  No position: not investing in the asset (taking neither the long nor short position).

A quantum encoding of the possible solutions s_i to qubit states $\vert i\rangle$ uses two qubits per asset:

• (i)  
  $\vert 10\rangle \rightarrow \text{short position}$.
• (ii)  
  $\vert 01\rangle \rightarrow \text{long position}$.
• (iii)  
  $\vert 00\rangle$ or $\vert 11\rangle \rightarrow \text{no position}$.

Such that for a portfolio of size M we require $n = 2M$ qubits for complete representation of the solution space $\boldsymbol{\mathcal{S}}$.

The discrete mean-variance Markowitz model implements a cost function for a given combination of solutions that takes into account the historical behavior of the assets. It may be expressed as a minimization problem with respect to the above quantum mapping as a

$$\begin{equation} C\left(s\right) = \omega \sum_{i,j = 1}^{M} \sigma_{ij}z_iz_j - \left(1 - \omega\right) \sum_{i = 1}^{M} r_iz_i, \end{equation} \tag{ 11 }$$

where σ_ij is the covariance between assets i and j, r_i is the average return, $0 \unicode{x2A7D} \omega \unicode{x2A7D} 1$ and $z_i \in \{1,-1,0\}$ represents a choice of long, short and no position respectively. As $\omega \rightarrow 0$, the optimal portfolio is one providing maximum returns, while as $\omega \rightarrow 1$, the optimal portfolio is the one that minimizes asset volatility. Equation (11) is subject to the constraint

$$\begin{equation} a = \sum_{i = 1}^{M} z_i, \end{equation} \tag{ 12 }$$

which maintains the relative net position for a pre-existing portfolio.

Altogether, the key steps of the QWOA algorithm as applied to the portfolio optimization problem are:

• (i)  
  Preparation of $\vert \psi_0\rangle$ over $\boldsymbol{\mathcal{S}}^{\prime}$ (solutions satisfying constraint a).
• (ii)  
  Computation of $\vert \boldsymbol{\gamma}, \boldsymbol{t}\rangle$ by p applications of $U_Q(\gamma)$ with phase-shift quality operator $Q = \sum_{\boldsymbol{\mathcal{S}}}C(s)\vert s\rangle\langle s\vert$ and indexed-walk unitary $U_w(t)$.
• (iii)  
  Minimization of $f(\boldsymbol{\gamma}, \boldsymbol{t}) = \sum_{s \in \boldsymbol{\mathcal{S}}^{\prime}}C(s)\left|\langle s\,|\,\boldsymbol{\gamma}, \boldsymbol{t}\rangle\right|^2$ via classical optimization of variational parameters $\boldsymbol{\gamma} = (\gamma_0, \dots, \gamma_{p-1})$ and $\boldsymbol{t} = (t_0, \dots, t_{p-1})$.

The QWOA algorithm, originating from a quantum walk perspective [46], can be implemented using a sequence of quantum gates, as detailed in [41]. Extensive simulations [42] highlight QWOA's advantages over other QAOA-based algorithms, notably due to its significantly reduced search space, achieving high-quality portfolios with fewer iterations. Furthermore, the inherent global symmetry in the QWOA mixing operator contributes to an unbiased convergence to optimal solutions with sufficient quantum circuit depth. QWOA has also shown great promise in handling diverse optimization challenges with intricate constraints and solution domains. The portfolio optimization problem can be efficiently solved using the QWOA with a gate complexity of $\mathcal{O}(m^2)$ [42]. These theoretical results depend on the efficient encoding of Q into the phase of the $\vert i\rangle$ states, implementation of the indexed-mixing operator $U_w(t)$, and identification of (at least) locally optimal variational parameters γ and t.

In this work, we consider a partial physical implementation of the QWOA as applied to the portfolio optimization problem to experimentally verify the effectiveness of the alternating phase-shift and quantum-walk operations. We consider a restricted phase-shift operator $U_{Q^{\prime}}$ which encodes the cost of only the valid solutions over a Hilbert space of size $|\boldsymbol{\mathcal{S}^{\prime}}|$ and a corresponding quantum walk unitary $U_{w^{\prime}}$. Within this framework, solutions in $\boldsymbol{\mathcal{S}^{\prime}}$ are computed and indexed classically.

A portfolio consisting of three stocks (Google, IBM and Microsoft) is considered under the constraint a = 0, for which there are seven valid solutions ($|\boldsymbol{\mathcal{S}^{\prime}}| = 7$). We show the seven valid solutions in the supplemental material. The quality values q_i have been calculated classically from the adjusted closing prices between 1/1/2019 and 12/31/2020 [47]. These are used to simulate the QWOA using the software package QuOp_MPI to obtain optimal values of t and γ for p = 1 to 8 [48]. Optimization of γ and t is carried out using the BFGS algorithm [49] with initial values generated following a uniform distribution between 0 and 2π for 20 repeats at each p.

The restricted QWOA operators $U_{Q^{\prime}}$ and $U_{w^{\prime}}$ are implemented with single photons following the experimental setup as illustrated in figure 1(a). The seven computational basis states correspond to the QWOA indexes representing valid portfolios. In a 7-dimensional qudit, the basis states are encoded as $\vert 0\rangle = \vert V1\rangle$, $\vert 1\rangle = \vert H2\rangle$, $\vert 2\rangle = \vert V2\rangle$, ⋯, $\vert 5\rangle = \vert H4\rangle$, $\vert 6\rangle = \vert V4\rangle$. Here, i ($i = 1, \cdots, 4$) denotes the spatial modes of the single photons, and $H~(V)$ indicates the horizontal (vertical) polarization of the photons.

Zoom In Zoom Out Reset image size

In this experiment, the initial state of the qudit is prepared in an equal superposition state represented as $\vert s\rangle = \frac{1}{\sqrt{7}}\sum_{i = 0}^{6}\vert i\rangle$. Initially, the photons pass through a polarizing beam splitter (PBS) and a half-wave plate (HWP) set at $-20.4^\circ$, which divides the transmitted photons with different polarizations into two parallel paths using a beam displacer (BD). Next, in the upper and lower modes, two additional HWPs at $-27.4^\circ$ and $-22.5^\circ$ are inserted, respectively. Finally, the initial state is prepared as $\vert s\rangle$ by passing through the second BD, followed by four HWPs set at 90^∘, 22.5^∘, $-22.5^\circ$, and 22.5^∘.

We simulate the probability amplification process in the QWOA using the mixing operator $U_{w^{\prime}}(t)$ and the phase shift operator $U_{Q^{\prime}}(\gamma)$. The mixing operator is the CTQW $U_{w^{\prime}}(t) = e^{-iWt}$ on the complete graph K₇ with adjacency matrix W shown in figure 1(b), which connects all valid solutions. The phase shift operator is a diagonal unitary that applies a phase shift proportional to the variational parameter γ and the value of the mean-variance Markowitz model, $U_{Q^{\prime}}(\gamma) = e^{-i\gamma Q}$, where the quality operator $Q = \textrm{diag}(q_0, q_1, \cdots, q_6)$ is a diagonal operator defined by solution qualities. Thus, we apply the unitary operation $U_p = \prod_{i = 1}^{p}U_{w^{\prime}}(t_i)U_{Q^{\prime}}(\gamma_i)$ ($p = 1, \cdots, 6$) on the state $\vert s\rangle$. Details of the variational parameters t_i and γ_i can be found in supplemental material.

With the method introduced in [1, 50, 51], an arbitrary n × n unitary matrix can be decomposed into a product of $\frac{1}{2}n^2-\frac{1}{2}n$ two-level unitary matrices, which act non-trivially only on two-dimensional subspaces of the n-dimensional Hilbert space. In this experiment, each $7\times 7$ unitary operator U_p ($p = 1, \cdots, 6$) can be decomposed as

$$\begin{equation} U_p = U_{7,6} \cdots U_{7,1} U_{6,5}\cdots U_{6,1}\cdots U_{2,1}, \end{equation} \tag{ 13 }$$

where $U_{i,j}$ are two-level unitary matrices. A two-level unitary matrix takes the form with only four elements $E_{i,i}$, $E_{i,j}$, $E_{j,i}$ and $E_{j,j}$, which are neither 0 or 1. For the remaining elements, all the diagonal elements are set to 1, while all the off-diagonal elements are set to 0.

According to the requirement of realizing a two-level matrix $U_{i,j}$ [52], the horizontal and vertical polarization states of the photons in different paths are recombined into polarization states of the photons in the same path via BDs and 45^∘ HWPs. Then we apply a $2\times2$ unitary transformation $\left(\begin{array}{ll}E_{i, i} & E_{i, j} \\ E_{j, i} & E_{j, j}\end{array}\right)$ to the polarization state of this path, which can be realized via a set of wave plates (WPs), thus producing an exact correspondence with equations (4) and (8).

The complexity of our experimental setup is independent of the number of the iterations and only depends on the dimension of U_p which in our experiment is always 7. In this case, we employ a total of 10 beam displacers (BDs) to implement each U_p . The variational parameters, which can be considered as parameters of $U_{i,j}$, are adjusted by configuring the setting angles of WPs. A detailed resource analysis of our experimental implementation of QWOA mixing unitary with the bulk optics is included in Supplemental Material, together with the corresponding system-agnostic quantum circuits.

Once the transformation U_p is implemented, we determine the square of the overlap between the final state $\vert \psi\rangle$ and the basis states $\vert i\rangle$ through a projection measurement. To perform the projective measurement, a PBS is employed to map the basis states of the qudit state onto separate spatial modes. By calculating the proportion of photon counts in each spatial mode to the total photon counts, we estimate the probability of the photons being measured in each basis state.

Theoretical predictions and experimental results of the probability distribution for the final states are presented in figure 2, covering iterations p = 1 to 6. The probabilities of the final state projected onto basis state $\vert 2\rangle$ after the 1st, 2nd, ⋯, and 6th iterations are $0.324\pm 3$, $0.502\pm 3$, $0.763\pm 3$, $0.812\pm 2$, $0.816\pm 2$ and $0.797\pm 2$, respectively. These results align closely with the theoretical predictions of 0.363, 0.567, 0.876, 0.970, 0.998 and 0.998, respectively. Thus, the probability of the optimal portfolio occupying the valid solution space increases with the number of iterations and finally tends to be stable. The experimental results show that we will find the highest-quality portfolio with a probability close to 1.

Zoom In Zoom Out Reset image size

In addition, we measure the expectation value of the quality operator that encodes the objective function, given by $\langle \psi\vert Q\vert \psi\rangle = \sum_{i = 0}^{6}q_i|\langle i|\psi\rangle|^2$, as shown in figure 3. It can be seen that as the number of iterations increases, the value of the objective function is in a gradual downward trend, which indicates that the probability of the high-quality portfolio is gradually increasing. In our experiment, the imperfections of the results include the imperfections of the interferometers and inaccuracies of WPs. The former leads to the dephasing whose impact on the objective function can be estimated by assuming the dephasing rate is $\thicksim0.97$.

Zoom In Zoom Out Reset image size

The implementation of quantum optimization algorithms on current quantum processors is still constrained by the compounding effects of system noise [31–33]. In contrast, we present the realization of a variational search up to p = 6 with strong agreement with predictions for a noise-free system. At all considered p, minimization of the objective function amplifies probability at the basis state corresponding to the highest-quality solution.

We consider a 3 asset portfolio where, as two qubits are used to encode the position per asset, the complete Hilbert space is of size $2^{2\times 3} = 64$. However, the QWOA reduces the quantum search to a globally-symmetric sub-space of seven valid solutions. For the portfolio optimization problem, compared to the QAOA, this significantly reduces the search space and eliminates bias resulting from mixing asymmetry.

In this paper, we present compelling evidence that demonstrates the applicability of quantum algorithms in solving portfolio optimization problems. Our experimental approach is direct, flexible, and holds the potential for scalability. The exploration of quantum algorithms in practical applications is gaining momentum [53–55], even though they are currently in a preliminary stage. With the dedicated efforts of scientific researchers, we anticipate that quantum technology will soon be leveraged to tackle challenging real-life problems.

This work is supported by the National Key R\&D Program of China (Grant No. 2023YFA1406701) and the National Natural Science Foundation of China (Grant Nos. 92265209, 12025401, 12104009 and 12305008) and an Australian Government Research Training Program Scholarship at the University of Western Australia. D K Q acknowledge support from the China Postdoctoral Science Foundation (Grant No. 2023M730198) and the fellowship of China National Postdoctoral Program for Innovative Talents (Grant No. BX20230036). J B W wishes to express gratitude to Yuying Li and Song Wang for valuable and insightful discussions, particularly within the domain of computational finance.

No new data were created or analysed in this study.

The authors declare no conflict of interest.

The mean-variance Markowitz model expresses the problem of choosing the optimal set of positions through minimization of the objective function,

$$\begin{equation} C\left(s\right) = \omega \sum_{i,j = 1}^{M} \sigma_{ij}z_iz_j - \left(1 - \omega\right) \sum_{i = 1}^{M} r_iz_i, \end{equation} \tag{ A.1 }$$

which is subject to,

$$\begin{equation} a = \sum_{i = 1}^{M} z_i, \end{equation} \tag{ A.2 }$$

for the case of discrete asset constraints. Variable r_i represents the expected returns for a given asset, σ_ij represents the covariance between assets i and j, $z_i \in \{1,-1,0\}$ represents long, short and no position respectively, and ω is between 0 and 1 [42]. Equation (A.1) formalizes the idea of asset diversification. For a pair of assets in a given portfolio, the first term is minimized by taking the same position on both assets if they are anti-correlated or by taking a short position on one asset and a long position on the other if they are correlated. The second term is minimized by taking a long position on appreciating assets and a short position on depreciating assets. The ω parameter balances the influence of risk and potential return with the value chosen according to the preference of the portfolio manager. Equation (A.2) arises in the context of portfolio rebalancing, with constraint value a preserving the net position with respect to a pre-existing portfolio.

With M assets with three possible positions per asset (long, short or none), the size of the unconstrained solution space is 3^M . Introducing the constraint defined in equation (A.2), the number of valid solutions with net position a is

$$\begin{equation} \mathcal{M}\left(M, a\right) = \sum_{j = 0}^M\left(\begin{array}{c} M \\ j \end{array}\right)\left(\begin{array}{c} M-j \\ \frac{1}{2}\left(M+a-j\right) \end{array}\right), \end{equation} \tag{ A.3 }$$

where, if the bottom parameter of the rightmost binomial coefficient is not an integer, it is set to 0 [42]. While $\mathcal{M}(M, a)$ is always less than 3^M , for any possible a, the number of valid solutions nevertheless grows exponentially with M.

Consider a portfolio of three stocks: Google, IBM and Microsoft, under the constraint that we have the same number of long and short positions (a = 0). This results in a valid solution space spanning 7 computational basis states.

The daily closing and opening prices of the three stocks between the $1/1/2019$ and $12/31/2020$ were obtained from Yahoo Finance, and the expected returns $(r_i)$ and covariance matrix (σ_ij ) calculated see table A1. The 7 valid solution qualities, as given by the mean-variance Markowitz model with ω = 0.5 are shown in table A2. These define the diagonal, $\boldsymbol{q} = q_0, \ldots, q_6$, of the Q operator. Note that, for a problem of this size, the q_i are of the order 10⁻⁴. To assist in the optimization process all q_i have been scaled by a constant factor of 250.

Table A2. Quality vector, q , describes the 7 valid solutions under the given the constraint of a = 0. The position on each stock (Google, IBM and Microsoft) can be read from the 6 qubits, where $\vert 01\rangle$ denotes a long position, $\vert 10\rangle$ denotes a short position and $\vert 00\rangle$ is no position.

QWOA Index $\vert \text{id}_{\boldsymbol{\chi}}(i)\rangle$	quantum encoding $\vert i\rangle$	qi
$\vert 0\rangle$	$\vert 000000\rangle$	0.0
$\vert 1\rangle$	$\vert 000110\rangle$	0.19537791107281657
$\vert 2\rangle$	$\vert 001001\rangle$	−0.1119084583542568
$\vert 3\rangle$	$\vert 010010\rangle$	0.09592313540638248
$\vert 4\rangle$	$\vert 011000\rangle$	−0.038186379903069514
$\vert 5\rangle$	$\vert 100001\rangle$	−0.05124555365045395
$\vert 6\rangle$	$\vert 100100\rangle$	0.12193130046716749

The optimal variational parameters of t_i and γ_i are obtained by the software package QuOp_MPI, for which detailed instructions and simulation code for the portfolio optimization problem are provided in [48, 56]. We experimentally realize the unitary operator of p-iteration QWOA ($p = 1, \cdots, 6$) as follows

$$\begin{equation} \begin{array}{ll} U_1 = U_{w^{^{\prime}}}\left(0.688\right) U_{Q^{^{\prime}}}\left(6.915\right), \\ U_2 = U_{w^{^{\prime}}}\left(0.756\right) U_{Q^{^{\prime}}}\left(12.139\right) U_{w^{^{\prime}}}\left(0.606\right) U_{Q^{^{\prime}}}\left(6.304\right), \\ U_3 = U_{w^{^{\prime}}}\left(0.716\right) U_{Q^{^{\prime}}}\left(13.972\right) U_{w^{^{\prime}}}\left(0.669\right) U_{Q^{^{\prime}}}\left(13.806\right) U_{w^{^{\prime}}}\left(0.596\right) U_{Q^{^{\prime}}}\left(7.482\right), \\ U_4 = U_{w^{^{\prime}}}\left(0.760\right) U_{Q^{^{\prime}}}\left(15.512\right) U_{w^{^{\prime}}}\left(0.676\right) U_{Q^{^{\prime}}}\left(14.986\right) U_{w^{^{\prime}}}\left(0.635\right) U_{Q^{^{\prime}}}\left(13.501\right) \\ \quad \quad \,\, U_{w^{^{\prime}}}\left(0.536\right) U_{Q^{^{\prime}}}\left(6.957\right), \\ U_5 = U_{w^{^{\prime}}}\left(0.801\right) U_{Q^{^{\prime}}}\left(15.134\right) U_{w^{^{\prime}}}\left(0.724\right) U_{Q^{^{\prime}}}\left(15.861\right) U_{w^{^{\prime}}}(0.656) U_{Q^{^{\prime}}}(14.611) \\ \quad \quad \,\, U_{w^{^{\prime}}}(0.591) U_{Q^{^{\prime}}}(12.985) U_{w^{^{\prime}}}(0.528) U_{Q^{^{\prime}}}(6.078), \\ U_6 = U_{w^{^{\prime}}}(0.801) U_{Q^{^{\prime}}}(15.134) U_{w^{^{\prime}}}(0.724) U_{Q^{^{\prime}}}(15.860) U_{w^{^{\prime}}}(0.656) U_{Q^{^{\prime}}}(14.610) \\ \quad \quad \,\, U_{w^{^{\prime}}}(0.591) U_{Q^{^{\prime}}}(12.978) U_{w^{^{\prime}}}(1.383) U_{Q^{^{\prime}}}(0.005) U_{w^{^{\prime}}}(0.042) U_{Q^{^{\prime}}}(6.070) . \end{array} \end{equation} \tag{ B.1 }$$

The key to running QWOA on hardware is to implement the unitary matrix $U = \prod _{i = 1}^pU_w(t_i)U_Q(\gamma_i)$. For larger databases, the bulk optical setup can realize QWOA with the n × n arbitrary unitary matrix and the number of the beam displacers $N_\mathrm{BD}$ increases with n linearly, i.e. $N_\textrm{BD}\sim 2n-4$, as shown in figure C1. An arbitrary n × n unitary matrix can be decomposed into a product of $\frac{1}{2}n^2-\frac{1}{2}n$ two-level unitary matrices, which act non-trivially only on two-dimensional subspaces of the n-dimensional Hilbert space:

$$\begin{equation} U = U_{n,n-1} \cdots U_{n,1} U_{n-1,n-2}\cdots U_{n-1,1}\cdots U_{2,1}. \end{equation} \tag{ C.1 }$$

Therefore, such an implementation does make a significant contribution to the general scheme of QWOA.

Zoom In Zoom Out Reset image size

In this section, we present the quantum circuit of QWOA to demonstrate that the algorithm can also be implemented in other systems. The design of quantum circuits is based on the cosine-sine decomposition (CSD) [57, 58]. As an example, we describe the factorization of an arbitrary $8\times 8$ unitary matrix using the CSD. For the $8 \times 8$ unitary matrix U₈, it can be decomposed as $U_8 = \mathbb{L}_8(S_4\oplus \mathbb{I}_4)\mathbb{R}_8$, where $\mathbb{L}_8$ and $\mathbb{R}_8$ are block diagonal,

$$\begin{equation} \mathbb{L}_8 = \left[\begin{array}{c|c} L_4 & 0 \\ \hline 0 & L^{^{\prime}}_4 \end{array}\right], \quad \mathbb{R}_8 = \left[\begin{array}{c|c} R_4 & 0 \\ \hline 0 & R^{^{\prime}}_4 \end{array}\right], \end{equation} \tag{ D.1 }$$

and L₄, $L^{^{\prime}}_4$, R₄, and $R^{^{\prime}}_4$ are $4\times 4$ unitary matrix. The matrix S₄ is an orthogonal cosine-sine matrix,

$$\begin{equation} S_4 = \left[\begin{array}{cc|cc} \text{cos}\theta_1 & 0 &\text{sin}\theta_1 & 0 \\ 0 & \text{cos}\theta_2 & 0 & \text{sin}\theta_2 \\ \hline -\text{sin}\theta_1 & 0 & \text{cos}\theta_1 & 0 \\ 0 & -\text{sin}\theta_2 & 0 & \text{cos}\theta_2, \end{array}\right]. \end{equation} \tag{ D.2 }$$

For convenience, the unitary transformations $\mathbb{L}_8$ and $\mathbb{R}_8$ can be written as

$$\begin{equation} \begin{array}{ll} & \mathbb{L}_8 = |0\rangle\langle 0|\otimes L_4+|1\rangle\langle 1| \otimes L_4^{^{\prime}} \\ & \mathbb{R}_8 = |0\rangle\langle 0|\otimes R_4+|1\rangle\langle 1| \otimes R_4^{^{\prime}} . \end{array} \end{equation} \tag{ D.3 }$$

Then the $8\times 8$ unitary transformations can be implemented by three controlled three-qubit transformations in figure D1(a). In addition, the $4\times 4$ unitary matrix U₄ (L₄, $L^{^{\prime}}_4$, R₄, and $R^{^{\prime}}_4$) can be decomposed as $U_4 = \mathbb{L}S_4\mathbb{R}$, where

$$\begin{equation} \mathbb{L} = \left[\begin{array}{c|c} L & 0 \\ \hline 0 & L^{^{\prime}} \end{array}\right], \quad \mathbb{R} = \left[\begin{array}{c|c} R & 0 \\ \hline 0 & R^{^{\prime}} \end{array}\right], \end{equation} \tag{ D.4 }$$

are block-diagonal, and L, Lʹ, R, and Rʹ are $2\times 2$ unitary transformations. The $4\times 4$ unitary transformations $\mathbb{L}$, S₄, and $\mathbb{R}$ can be rewritten as

$$\begin{equation} \begin{array}{cc} &\mathbb{L} = |0\rangle\langle 0|\otimes L+|1\rangle\langle 1|\otimes L^{^{\prime}} \\ \nonumber &S_4 = S\otimes |0\rangle\langle 0|+S^{^{\prime}}\otimes |1\rangle\langle 1| \\ \nonumber &\mathbb{R} = |0\rangle\langle 0|\otimes R+|1\rangle\langle 1|\otimes R^{^{\prime}}, \end{array} \end{equation}$$

where $S = \left[\begin{array}{cc} \text{cos}\theta_1 & \text{sin}\theta_1 \\ -\text{sin}\theta_1 & \text{cos}\theta_1 \end{array}\right]$ and $S^{^{\prime}} = \left[\begin{array}{cc} \text{cos}\theta_2 & \text{sin}\theta_2 \\ -\text{sin}\theta_2 & \text{cos}\theta_2 \end{array}\right]$. Then these $4\times 4$ unitary transformations can be implemented by three controlled two-qubit transformations in figures D1(b) and (c). The decomposition method can be used to decomposed higher-dimensional unitary operators. Therefore, our idea can in principle design quantum circuit of any dimension, so as to realize the unitary matrix required by QWOA for generalization in other systems.

Zoom In Zoom Out Reset image size