Picking optimal portfolios using Quantum Computing

In this article, I aim to explain how to solve NP-hard problems like Portfolio Optimization using Quantum Computers. I use the prevalent algorithms VQE (that I will explain briefly) for the same.

Quantum Computers might be a thing of the future, but we can already do so much from them. IBMQ, Rigetti, and Google are providing us free cloud service of these multi-million dollar machines. There has never been a better time to be alive! Let us glimpse into how can we solve some challenging problems through the near term quantum devices.

What is Portfolio Optimization?

Portfolio Optimization is the art of selecting the best portfolios of stocks from all the given portfolios that maximizes the expected returns and minimizes the financial risk.

The expected return is the anticipated profit or loss of investment. It is determined by (new_price/old_price)-1

The financial risk is basically the volatility of the returns and is measured through statistical instruments like variance, standard deviation, and covariance.

Now we will try to operate the Variational Quantum Eigensolver (VQE) Algorithm to optimize the portfolio using Pennylane and Qiskit.

Understanding VQE

The Variational Quantum Eigensolver (VQE) algorithm is one of the essential algorithms to learn in the NISQ era. It makes the idea of the famous Quantum Phase Estimation Algorithm usable with the small number of qubits Quantum Computers we currently possess. VQE presents us with a hybrid quantum-classical approach to tackle the current hardware limits.

VQE exploits the fact that the ground state energy is always less than or equal to the expectation value of the Hamiltonian of the state. Hence, by minimizing the expectation value of the Hamiltonian, we can find the upper bound of the ground state energy.

Thus by varying the state, we can find the wavefunction on which the value of the expectation value is the lowest.

VQE, in a straightforward language, consists of three steps:

1. Converting the Hamiltonian into Pauli basis

We need to convert our Hamiltonian of the state into Pauli’s basis, i.e., in multiples of X, XY, XXX, etc. matrices.

For one qubit system:

We need to do so because the Pauli operators are commutable. We need commutable operators to measure the states without actually spoiling the states in both the basis.

There are multiple ways we can achieve this, WeightedPauliOperator in the Qiskit Library is one such technique.

2. Creation of a variational form

The algorithm goes over every state and checks which one has the minimum expectation value to determine our lowest eigenvalue of the state. For checking every state we have to create a parametrized circuit (ansatz) that covers each and every possible state. Picking a ‘good’ ansatz from all possible is an art.

3. Parameter Optimization

Now the easy classical work where you have to optimize your parameters for giving you the minimum expectation value. You can use any prebuilt classical optimizer (like Gradient Descent) for it or even create your own.

My implementation of Portfolio Optimization

My main inspiration for implementing Portfolio Optimization is Qiskit Tutorials. I have taken the same assumptions as given there to draft my algorithm. I will try to pick two stocks out of a total of 3 assets in this code. I will be basically implementing VQE from scratch using Pennylane and Qiskit.

Let’s get started,

Import the required libraries

import numpy as np
import matplotlib.pyplot as plt
import datetime
import pandas as pd
from nsepy import get_history
from docplex.mp.model import Model
from qiskit.optimization import QuadraticProgram
from qiskit.optimization.converters import LinearEqualityToPenalty
import pennylane as qml
from qiskit.optimization.applications.ising.common import sample_most_likely

I picked three random stocks, and I extract their closed prices for one month to find the expected returns and the covariance.

num_assets = 3
# Generate expected return and covariance matrix from (random) time-series
stocks = [("TICKER%s" % i) for i in range(num_assets)]
data = RandomDataProvider(tickers=stocks,
                 start=datetime.datetime(2016,1,1),
                 end=datetime.datetime(2016,1,30))
data.run()

Our objective is to maximize the following equation:

where:

We assume the following simplifications: — all assets have the same price (normalized to 1) — the full budget B has to be spent, i.e., one has to select exactly B assets.

Creating the Quadratic program for the following equation

Creating a Quadratic Program will help us to develop our Hamiltonian, which is the first step for the VQE.

#Calculate mean returns
mu = data.get_period_return_mean_vector() 
 
#Calculate mean covariance
sigma = data.get_period_return_covariance_matrix()
q = 0.5                   # set risk factor
budget = 2          # set budget (B)
penalty = num_assets      # set parameter to scale the budget penalty term
mdl = Model('docplex model')
x = mdl.binary_var_list(num_assets)
# set objective function: 
#
#     maximize { mu^T * x - q * x^T * sigma * x }
#
objective = mdl.sum([mu[i] * x[i] for i in range(num_assets)])  
# mu^T * x
objective -= q * mdl.sum([sigma[i][j]*x[i]*x[j] for i in range(num_assets) for j in range(num_assets)])
mdl.maximize(objective)
# add budget constraint: 
#
#     1^T * x == budget
#
cost = mdl.sum([x[i] for i in range(num_assets)])
mdl.add_constraint(cost == budget, ctname='budget')
# converting to Quadratic Program
mod = QuadraticProgram()        
mod.from_docplex(mdl)
#removing the constraint to create the QUBO
lineq2penalty = LinearEqualityToPenalty(penalty) 
qubo = lineq2penalty.convert(mod)
#converting QUBO to an Ising Hamiltonian
H, offset = qubo.to_ising()
H = H.to_legacy_op()  #coverting it to a legacy operator
H.print_details()

Here we have our Hamiltonian as the sum of Pauli Operators just like we wanted. But it doesn’t look obvious, so we will change it into a dictionary.

Change the ‘Hamiltonian’ into a ‘dictionary’ of Pauli operators just for the ease.

def pauli_operator_to_dict(pauli_operator):
    
    d = pauli_operator.to_dict()
    paulis = d['paulis']
    paulis_dict = {}
for x in paulis:
        label = x['label']
        coeff = x['coeff']['real']
        paulis_dict[label] = coeff
return paulis_dict
pauli_dict = pauli_operator_to_dict(H)
print(pauli_dict)

Now comes the hard part, i.e., creating an ansatz. I use the following two-layer ansatz with nine parameters.

def ansatz(theta):
    
    qml.RX(theta[0], wires=0)
    qml.RX(theta[1], wires=1)
    qml.RX(theta[2], wires=2)
    qml.CNOT(wires=[0,1])
    qml.CNOT(wires=[0,2])
    qml.CNOT(wires=[1,2])
    qml.RX(theta[3], wires=0)
    qml.RX(theta[4], wires=1)
    qml.RX(theta[5], wires=2)
    qml.CNOT(wires=[0,1])
    qml.CNOT(wires=[0,2])
    qml.CNOT(wires=[1,2])
    qml.RX(theta[6], wires=0)
    qml.RX(theta[7], wires=1)
    qml.RX(theta[8], wires=2)
    
    return ansatz

Functions for returning the ‘expectation value’ in a ‘particular basis’:

@qml.qnode(dev1)
def circuit_IIZ(params):
    ansatz(params)
    return qml.expval(qml.Identity(0)@qml.Identity(1)@qml.PauliZ(2))
@qml.qnode(dev1)
def circuit_IIZ(params):
    ansatz(params)
    return qml.expval(qml.Identity(0)@qml.Identity(1)@qml.PauliZ(2))
@qml.qnode(dev1)
def circuit_IZI(params):
    ansatz(params)
    return qml.expval(qml.Identity(0)@qml.PauliZ(1)@qml.Identity(2))
@qml.qnode(dev1)
def circuit_ZII(params):
    ansatz(params)
    return qml.expval(qml.PauliZ(0)@qml.Identity(1)@qml.Identity(2))
@qml.qnode(dev1)
def circuit_IZZ(params):
    ansatz(params)
    return qml.expval(qml.Identity(0)@qml.PauliZ(1)@qml.PauliZ(2))
@qml.qnode(dev1)
def circuit_ZIZ(params):
    ansatz(params)
    return qml.expval(qml.PauliZ(0)@qml.Identity(1)@qml.PauliZ(2))
@qml.qnode(dev1)
def circuit_ZZI(params):
    ansatz(params)
    return qml.expval(qml.PauliZ(0)@qml.PauliZ(1)@qml.Identity(2))
@qml.qnode(dev1)
def circuit_IZI(params):
    ansatz(params)
    return qml.expval(qml.Identity(0)@qml.PauliZ(1)@qml.Identity(2))
@qml.qnode(dev1)
def circuit_ZII(params):
    ansatz(params)
    return qml.expval(qml.PauliZ(0)@qml.Identity(1)@qml.Identity(2))
@qml.qnode(dev1)
def circuit_IZZ(params):
    ansatz(params)
    return qml.expval(qml.Identity(0)@qml.PauliZ(1)@qml.PauliZ(2))
@qml.qnode(dev1)
def circuit_ZIZ(params):
    ansatz(params)
    return qml.expval(qml.PauliZ(0)@qml.Identity(1)@qml.PauliZ(2))
@qml.qnode(dev1)
def circuit_ZZI(params):
    ansatz(params)
    return qml.expval(qml.PauliZ(0)@qml.PauliZ(1)@qml.Identity(2))

Now, we write our final function VQE (objective function) that we want to minimize for measuring the minimum expectation value.

def vqe(parameters):
    
    expval_ZII = pauli_dict['ZII'] * circuit_ZII(parameters)  
    expval_IIZ = pauli_dict['IIZ'] * circuit_IIZ(parameters)
    expval_IZI = pauli_dict['IZI'] * circuit_IZI(parameters)
    expval_ZZI = pauli_dict['ZZI'] * circuit_ZZI(parameters)
    expval_ZIZ = pauli_dict['ZIZ'] * circuit_ZIZ(parameters)
    expval_IZZ = pauli_dict['IZZ'] * circuit_IZZ(parameters)
    
    # summing the measurement results
    total_expval= expval_ZII + expval_IIZ + expval_ZIZ +  expval_ZZI + expval_ZIZ + expval_IZZ
    
    return total_expval

Finally, we optimize our objective function using a classical optimizer Gradient Descent present in the Pennylane library.

opt = qml.GradientDescentOptimizer(stepsize=0.5)
value = []
    # optimize parameters in objective
params = np.random.rand(9)
steps = 500
for i in range(steps):
    params = opt.step(vqe, params)
    value.append(vqe(params))
plt.plot(value)

Create the final circuit to find the most probable eigenstate:

@qml.qnode(dev1)
def final_circ(params):
    ansatz(params)
    return qml.probs(wires=[0,1,2])

Plot the eigenstates:

plt.bar(['000', '001', '010', '011', '100', '101', '110', '111'], final_circ(params[0]))

It clearly tells us that the [1, 0, 1] is the most optimal

We can also check our answer by a classical method:

exact_eigensolver = NumPyMinimumEigensolver(H)
result = exact_eigensolver.run()
sample_most_likely(result.eigenstate)

Voila! It matches! Our VQE works perfectly.

This is all for now; hope you gained some insights about the algorithms VQE and how they are used to solve challenging problems.

Resources if you want to dive deep

1. Medium article on classical portfolio optimization

2. Blog on VQE

References

1. Peruzzo, Alberto, et al. “A variational eigenvalue solver on a photonic quantum processor.” Nature communications 5 (2014): 4213
2. Farhi, Edward, Jeffrey Goldstone, and Sam Gutmann. “A quantum approximate optimization algorithm.” arXiv preprint arXiv:1411.4028 (2014)

I am a Physics undergrad student at Indian Institute of Technology, Roorkee. If you’d like to know more about me and my work, check out my website, https://rochisha0.github.io/