Portfolio optimization with variational quantum eigensolver (VQE)-(2)

We will implement solving a real world data portfolio optimization with VQE in this post. Here we use the quantum simulation package provided by pennylane.

The first step is to download the stock price data from yahoo finance, and calculate the return vector and covariance matrix.

import pandas as pd
import yfinance as yf

assets = "AAPL MSFT AMZN TSLA GOOG BRK-B"
StockStartDate = '2018-01-01'
StockEndDate = '2018-12-31'
interval = '1d'

df = yf.download(assets, start=StockStartDate,\
                 end=StockEndDate, interval=interval)['Adj Close']

# daily return distribution
ret = df.pct_change().dropna()

# mean return vector & covariance matrix (annualized)
R = ret.mean()*252
Sigma = ret.cov()*252

Recall in the last post that we are minimizing the equation

and we have to convert x∈{0,1} to z∈{1,-1}, which end up with

you can find my calculation process at the end of this post.

import pennylane as qml
from pennylane import numpy as np

# define parameters in H
N = 6      # number of total assets
gamma = 1  # risk aversion coefficient
B = 3      # budget
P = 1.0    # penalty

ZZ = [qml.PauliZ(i)@qml.PauliZ(j) for i in range(N) for j in range(i+1,N)]
ZZ_coeff = [0.5*(gamma*Sigma.values[i][j] + P) for i in range(N) for j in range(i+1,N)]
Z = [qml.PauliZ(i) for i in range(N)]
Z_coeff = [-0.5*gamma*(sum(Sigma.values[i][:])) + 0.5*R[i] - 0.5*P*(N-2*B) for i in range(N)]
C = 0.25*gamma*(sum(sum(Sigma.values)) + np.trace(Sigma))- 0.5*sum(R) + 0.25*P*(N + (N-2*B)**2)

# Construct the problen Hamiltonian
obs = ZZ + Z
coeffs = ZZ_coeff + Z_coeff
H = qml.Hamiltonian(coeffs, obs, grouping_type="qwc")

Now we can start to build our VQE ansatz. Here we demonstrate with a heuristic ansatz, which is often called a hardware efficient ansatz, but we replace the U3 gate with just a RY gate.

# Design the ansatz
p = 2     # circuit repetitions
def ansatz(params, qubits, depth=p):
    for q in range(qubits):
            qml.RY(params[q], wires=q)
    for d in range(1,depth+1):
        for q in range(qubits-1):
            qml.CNOT(wires=[q,q+1])
        for q in range(qubits):
            qml.RY(params[d*qubits+q], wires=q)

dev = qml.device("default.qubit", wires=N)

# Set the cost function on dev
@qml.qnode(dev, diff_method="parameter-shift")
def cost(x):
    ansatz(x, qubits=N)
    return qml.expval(H)

# For analyze the optimized circuit
@qml.qnode(dev)
def probability_circuit(params):
    ansatz(params, qubits=N)
    return qml.probs(wires=range(N))

We then use classical optimizers to optimize the parameters in the ansatz, until we cannot minimize the cost function any more. Here we use the QNG optimizer, you can also try other optimizers provided by pennylane.

opt = qml.QNGOptimizer(stepsize=0.02)
max_steps=600
init_params = np.random.rand((p+1)*N)

params = init_params
old_cost = 9999.999999
for i in range(max_steps):
    params = opt.step(cost, params)
    obj_value = cost(params)
    if (i + 1) % 5 == 0:
        print("Cost after step {:5d}: {: .7f}".format(i + 1, obj_value))
        if np.round(old_cost, 7) == np.round(obj_value, 7):
            break
        else:
            old_cost = obj_value
    
print("Optimized parameters: {}".format(params))
print("Optimized objective function value: {}".format(obj_value))

The optimized parameter is then plug into the ansatz, and we can find our final solution by sampling from the optimized ansatz. If the optimization is successful, you should find the ground state having highest probability.

import matplotlib.pyplot as plt

probs = probability_circuit(params)
print("Final solution : {:06b}, with prob={:.5f}".format(np.argmax(probs), max(probs)))
plt.style.use("seaborn")
plt.bar(range(2 ** len(range(N))), probs)
plt.show()

Finally, let’s check the exact ground state solution via exhaustive search.

import itertools
all_combinations = itertools.product([0, 1], repeat=N)

E_g = 99999.99999
for x in all_combinations:
    E = gamma*np.dot(x, np.dot(Sigma, x)) - np.dot(R,x) + P*(sum(x)-B)**2
    if E < E_g:
        E_g = E
        sol = x
print("Exact solution:{}".format(sol))

We successfully find the ground state solution with VQE!

APPENDIX

The calculation process for converting x∈{0,1} to z∈{1,-1}, note that 1 and -1 is the eigenvalue for eigenstate |0> and |1> respectively, so we are converting 0 to 1 and 1 to -1, rather than 0 to -1 and 1 to 1.