Portfolio Optimization

Portfolio Optimization#

QAOA and VQE for Markowitz portfolio optimization. Mapping portfolio selection to QUBO and solving on quantum hardware.

The Markowitz Problem#

Given \(N\) assets with expected returns \(\boldsymbol{\mu} \in \mathbb{R}^N\) and covariance matrix \(\boldsymbol{\Sigma} \in \mathbb{R}^{N \times N}\), find the portfolio weights \(\mathbf{w}\) that maximize return for a given risk tolerance \(\lambda\):

\[\max_{\mathbf{w}} \quad \boldsymbol{\mu}^T \mathbf{w} - \lambda \mathbf{w}^T \boldsymbol{\Sigma} \mathbf{w}\]

subject to \(\sum_i w_i = 1\) and \(w_i \geq 0\).

For discrete asset selection (buy/don’t buy), this becomes a binary optimization problem, which is NP-hard for large portfolios.

QUBO Formulation#

Binary portfolio selection maps naturally to QUBO. Let \(x_i \in \{0, 1\}\) indicate whether asset \(i\) is included:

\[\min_{\mathbf{x}} \quad \lambda \mathbf{x}^T \boldsymbol{\Sigma} \mathbf{x} - \boldsymbol{\mu}^T \mathbf{x} + A\left(\sum_i x_i - K\right)^2\]

where \(K\) is the target number of assets and \(A\) is a penalty weight enforcing the cardinality constraint.

Quantum Approach: QAOA#

The Quantum Approximate Optimization Algorithm (QAOA) solves QUBO by encoding the cost function into a quantum Hamiltonian and alternating between cost and mixer unitaries:

\[|\boldsymbol{\gamma}, \boldsymbol{\beta}\rangle = \prod_{p=1}^{P} e^{-i\beta_p H_M} e^{-i\gamma_p H_C} |+\rangle^{\otimes n}\]

where \(H_C\) encodes the QUBO cost and \(H_M = \sum_i X_i\) is the mixer Hamiltonian.

from qiskit_optimization import QuadraticProgram
from qiskit_optimization.algorithms import MinimumEigenOptimizer
from qiskit_algorithms import QAOA
from qiskit_algorithms.optimizers import COBYLA
from qiskit.primitives import Sampler

# Define portfolio problem
qp = QuadraticProgram("portfolio")
for i in range(4):
    qp.binary_var(f"asset_{i}")

# Objective: minimize risk - return
qp.minimize(
    linear=-expected_returns,
    quadratic=risk_factor * covariance_matrix
)

# Solve with QAOA
qaoa = QAOA(sampler=Sampler(), optimizer=COBYLA(maxiter=200), reps=2)
optimizer = MinimumEigenOptimizer(qaoa)
result = optimizer.solve(qp)
print(f"Selected assets: {result.x}")

import cirq
import numpy as np

# Encode QUBO as Ising Hamiltonian
qubits = cirq.LineQubit.range(4)

def qaoa_circuit(gamma, beta, depth=2):
    circuit = cirq.Circuit()
    # Initial superposition
    circuit.append(cirq.H.on_each(*qubits))
    for p in range(depth):
        # Cost unitary (ZZ interactions from covariance)
        for i in range(4):
            for j in range(i+1, 4):
                circuit.append(cirq.ZZPowGate(exponent=gamma[p] * cov[i,j] / np.pi).on(qubits[i], qubits[j]))
        # Mixer unitary
        circuit.append(cirq.rx(2 * beta[p]).on_each(*qubits))
    circuit.append(cirq.measure(*qubits, key='result'))
    return circuit

import pennylane as qml
from pennylane import qaoa
import numpy as np

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

# Define cost Hamiltonian from covariance matrix
cost_h = qml.Hamiltonian(
    coefficients=[cov[i,j] for i in range(n_assets) for j in range(i+1, n_assets)],
    observables=[qml.PauliZ(i) @ qml.PauliZ(j) for i in range(n_assets) for j in range(i+1, n_assets)]
)
mixer_h = qaoa.x_mixer(range(n_assets))

@qml.qnode(dev)
def cost_function(gamma, beta):
    for i in range(n_assets):
        qml.Hadamard(wires=i)
    for p in range(depth):
        qaoa.cost_layer(gamma[p], cost_h)
        qaoa.mixer_layer(beta[p], mixer_h)
    return qml.expval(cost_h)

from pyquil import Program
from pyquil.gates import H, CNOT, RZ, RX
from pyquil.quilbase import Declare

def qaoa_portfolio(gamma, beta, n_assets=4):
    p = Program()
    ro = p.declare('ro', 'BIT', n_assets)

    # Initial superposition
    for i in range(n_assets):
        p += H(i)

    # Cost layer (ZZ from covariance)
    for i in range(n_assets):
        for j in range(i+1, n_assets):
            p += CNOT(i, j)
            p += RZ(2 * gamma * covariance[i,j], j)
            p += CNOT(i, j)

    # Mixer layer
    for i in range(n_assets):
        p += RX(2 * beta, i)

    for i in range(n_assets):
        p += MEASURE(i, ro[i])
    return p

Current Limitations#

Factor

Status

Problem size

NISQ: ~20-50 assets; fault-tolerant: 1000+

Circuit depth

QAOA depth \(p \geq 3\) needed for quality solutions

Noise impact

Significant on current hardware; error mitigation helps

Classical competition

Gurobi/CPLEX solve up to ~1000 assets optimally

The quantum advantage for portfolio optimization is expected when portfolios exceed ~100 assets with complex constraints that make classical solvers struggle.

Further Reading#