Portfolio Optimization using QUBO and Qiskit

This notebook demonstrates the use of QUBO (Quadratic Unconstrained Binary Optimization) for solving the optimal portfolio problem.

Problem

We have 10 stocks with different expected returns and risks. We want to select a portfolio that:

• Maximizes expected return
• Minimizes risk (volatility)
• Selects exactly 5 stocks (budget constraint)

Instructions for use

You can find and change the running circuit in the code to QAOA or QAOA optimized by commuting gates

# Create QAOA circuit
print("\nBuilding QAOA circuit...")
#HERE IS MY OPTIMIZED CIRCUIT
#qc = qaoa_circuit(gammas, betas, h, J, n_qubits)
qc = qaoa_circuit_commuting_gates(gammas, betas, h, J, n_qubits, 'xy', True)

Recommended settings:

QAOA:

P_Layers = 10
shots = 10000

QAOA optimized by commuting gates

P_Layers = 5
shots = 2000 #min 1 solution from every 1000 with 5 layers

# !pip install qiskit qiskit-optimization qiskit-algorithms qiskit-aer numpy matplotlib --break-system-packages

import numpy as np
import matplotlib.pyplot as plt
from qiskit_optimization import QuadraticProgram
from qiskit_optimization.algorithms import MinimumEigenOptimizer
from qiskit_optimization.converters import QuadraticProgramToQubo
from qiskit_algorithms import QAOA, NumPyMinimumEigensolver, SamplingVQE
from qiskit_algorithms.optimizers import COBYLA
from qiskit.primitives import StatevectorSampler
from qiskit import QuantumCircuit
from qiskit.circuit import Parameter
from qiskit.circuit.library import RXXGate, RYYGate, RZZGate
from qiskit_optimization.applications import Maxcut
from qiskit.quantum_info import SparsePauliOp
from itertools import combinations
import time
from collections import defaultdict

Portfolio Data Definition

P_LAYERS = 5
shots = 2000

asset_names = [
    "Apple (Tech)",
    "Microsoft (Tech)",
    "ExxonMobil (Energie)",
    "Chevron (Energie)",
    "JPMorgan (Finance)",
    "Bank of America (Finance)",
    "Johnson & Johnson (Healthcare)",
    "Pfizer (Healthcare)",
    "Walmart (Retail)",
    "Amazon (E-commerce)",
]

n_assets = len(asset_names)

expected_returns = np.array([
    24.30,  # Apple - high growth
    31.91,  # Microsoft - high growth
    8.50,   # ExxonMobil - dividends
    2.22,   # Chevron - dividends
    34.70,  # JPMorgan - very high growth
    18.88,  # Bank of America - medium
    7.88,  # Johnson & Johnson - healthcare growing
    -12.14,  # Pfizer - pharmaceuticals
    30.72,   # Walmart - stable
    29.12   # Amazon - high growth, but risky
])

volatility = np.array([
    27.49,  # Apple - high volatility
    24.37,  # Microsoft
    23.17,  # ExxonMobil - commodities volatility
    22.77,  # Chevron
    22.95,  # JPMorgan - financial sector
    25.92,  # Bank of America
    17.05,  # Johnson & Johnson - defensive
    24.31,  # Pfizer
    19.88,  # Walmart - very stable
    33.19   # Amazon - highest volatility
])

# - Tech-Tech: High (0.7-0.85)
# - Energie-Energie: High (0.75-0.8)
# - Finance-Finance: High (0.7-0.8)
# - Healthcare-Healthcare: Medium (0.5-0.6)
# - Tech-Finance: Mediu (0.4-0.5)
# - Tech-Energie: Low (0.2-0.3)
# - Healthcare-Energie: Very low (0.1-0.2)

covariance_matrix = np.array([
    # Apple, MSFT, Exxon, Chev, JPM, BoA, J&J, Pfiz, Walmart, Amazon
    [755.83, 370.69, 148.11, 176.32, 206.05, 240.59,  40.68, 130.87, 164.90, 449.34],  # Apple
    [370.69, 593.93,  36.07,  80.90, 164.99, 169.12,  -9.80,  73.54, 117.29, 517.96],  # Microsoft
    [148.11,  36.07, 536.95, 437.26, 201.04, 237.14,  66.76, 107.09,  56.04,  71.92],  # ExxonMobil
    [176.32,  80.90, 437.26, 518.55, 226.51, 260.45,  67.81, 125.63,  67.92, 124.28],  # Chevron
    [206.05, 164.99, 201.04, 226.51, 526.48, 446.34,  61.15, 105.61, 120.55, 256.44],  # JPMorgan
    [240.59, 169.12, 237.14, 260.45, 446.34, 672.05,  63.47, 127.60, 135.94, 287.93],  # Bank of America
    [ 40.68,  -9.80,  66.76,  67.81,  61.15,  63.47, 290.59, 180.15,  58.36, -40.56],  # Johnson & Johnson
    [130.87,  73.54, 107.09, 125.63, 105.61, 127.60, 180.15, 590.85,  61.37,  68.90],  # Pfizer
    [164.90, 117.29,  56.04,  67.92, 120.55, 135.94,  58.36,  61.37, 395.35, 177.45],  # Walmart
    [449.34, 517.96,  71.92, 124.28, 256.44, 287.93, -40.56,  68.90, 177.45,1101.62],  # Amazon
])

for i, name in enumerate(asset_names):
    print(
        f"{name:<15s}: Yield={expected_returns[i]:5.1f}%, Risk={volatility[i]:5.1f}%"
    )

Apple (Tech)   : Yield= 24.3%, Risk= 27.5%
Microsoft (Tech): Yield= 31.9%, Risk= 24.4%
ExxonMobil (Energie): Yield=  8.5%, Risk= 23.2%
Chevron (Energie): Yield=  2.2%, Risk= 22.8%
JPMorgan (Finance): Yield= 34.7%, Risk= 22.9%
Bank of America (Finance): Yield= 18.9%, Risk= 25.9%
Johnson & Johnson (Healthcare): Yield=  7.9%, Risk= 17.1%
Pfizer (Healthcare): Yield=-12.1%, Risk= 24.3%
Walmart (Retail): Yield= 30.7%, Risk= 19.9%
Amazon (E-commerce): Yield= 29.1%, Risk= 33.2%

QUBO Problem Formulation

We transform the problem into QUBO format:

Objective function:

Minimize: risk_factor × Σᵢⱼ xᵢ·Covᵢⱼ·xⱼ - return_factor × Σᵢ xᵢ·Rᵢ + penalty × (Σᵢ xᵢ - budget)²

where:

• xᵢ ∈ {0,1} indicates whether we select stock i
• Covᵢⱼ is the covariance matrix (risk)
• Rᵢ are the expected returns
• budget is the number of stocks we want to select

budget = 5              # Select 5 out of 10 stocks
risk_factor = 0.003     # Reduced due to larger covariances
return_factor = 3.0     # Increased for better motivation
penalty = 150.0         # Strong penalty
reps = 3                # More QAOA layers for more complex problem

print(f"  - Risk factor: {risk_factor}")
print(f"  - Return factor: {return_factor}")
print(f"  - Budget: {budget} stocks")
print(f"  - Penalty: {penalty}")

# Create QUBO problem
qp = QuadraticProgram('portfolio_optimization')

for i in range(n_assets):
    qp.binary_var(name=f'x_{i}')

# LINEAR PART
linear = {}
for i in range(n_assets):
    linear[f'x_{i}'] = -return_factor * expected_returns[i]
    linear[f'x_{i}'] += penalty * (1 - 2 * budget)

# QUADRATIC PART
quadratic = {}
for i in range(n_assets):
    for j in range(i+1, n_assets):
        key = (f'x_{i}', f'x_{j}')
        risk_term = 2 * risk_factor * covariance_matrix[i][j]
        penalty_term = 2 * penalty
        quadratic[key] = risk_term + penalty_term

# CONSTANT TERM
constant = penalty * budget ** 2

# Construct objective function
qp.minimize(linear=linear, quadratic=quadratic, constant=constant)

  - Risk factor: 0.003
  - Return factor: 3.0
  - Budget: 5 stocks
  - Penalty: 150.0

Solution using Classical Algorithm (Reference Solution)

Classical numpy solution for comparison:

• Optimal value: [value]
• Selected stocks:
  • [Stock name] (Return: X%, Risk: Y%)
• Total expected return: X%
• Total portfolio risk: Y
• Sharpe ratio: Z

Solution description:

• 1. Converts QUBO to Hamiltonian
• 2. Calculates all eigenvalues using NumPy/SciPy
• 3. Finds the minimum = optimal solution

exact_solver = NumPyMinimumEigensolver()
exact_optimizer = MinimumEigenOptimizer(exact_solver)
exact_result = exact_optimizer.solve(qp)

print(f"Objective function value: {exact_result.fval:.2f}")
print(f"\nDebug - x vector: {exact_result.x}")

selected_exact = []
print("\nSelected stocks:")
print("-" * 70)
for i in range(n_assets):
    value = exact_result.x[i]
    if value > 0.5:
        selected_exact.append(i)
        print(f"  ✓ x_{i} = {value:.3f} | {asset_names[i]:<15s} | "
              f"Return: {expected_returns[i]:5.1f}% | Risk: {volatility[i]:5.1f}%")
    else:
        print(f"  ✗ x_{i} = {value:.3f} | {asset_names[i]:<15s}")

print(f"\nTotal stocks selected: {len(selected_exact)} (target: {budget})")

if len(selected_exact) > 0:
    total_return_exact = sum(expected_returns[i] for i in selected_exact)
    total_risk_exact = np.sqrt(sum(covariance_matrix[i][j]
                                   for i in selected_exact
                                   for j in selected_exact))

    print(f"Total expected return: {total_return_exact:.2f}%")
    print(f"Total portfolio risk: {total_risk_exact:.2f}")
    print(f"Sharpe ratio: {total_return_exact/total_risk_exact:.3f}")
else:
    total_return_exact = 0
    total_risk_exact = 0

Objective function value: -436.98

Debug - x vector: [1. 1. 0. 0. 1. 0. 0. 0. 1. 1.]

Selected stocks:
----------------------------------------------------------------------
  ✓ x_0 = 1.000 | Apple (Tech)    | Return:  24.3% | Risk:  27.5%
  ✓ x_1 = 1.000 | Microsoft (Tech) | Return:  31.9% | Risk:  24.4%
  ✗ x_2 = 0.000 | ExxonMobil (Energie)
  ✗ x_3 = 0.000 | Chevron (Energie)
  ✓ x_4 = 1.000 | JPMorgan (Finance) | Return:  34.7% | Risk:  22.9%
  ✗ x_5 = 0.000 | Bank of America (Finance)
  ✗ x_6 = 0.000 | Johnson & Johnson (Healthcare)
  ✗ x_7 = 0.000 | Pfizer (Healthcare)
  ✓ x_8 = 1.000 | Walmart (Retail) | Return:  30.7% | Risk:  19.9%
  ✓ x_9 = 1.000 | Amazon (E-commerce) | Return:  29.1% | Risk:  33.2%

Total stocks selected: 5 (target: 5)
Total expected return: 150.75%
Total portfolio risk: 92.00
Sharpe ratio: 1.639

Custom QAOA Functions

from this site Link: pennylane.ai

def from_Q_to_Ising(Q, offset):
    """
    Convert the QUBO matrix Q into Ising parameters J and h.
    
    Parameters:
    - Q: QUBO matrix (numpy array)
    - offset: Initial offset value
    
    Returns:
    - h: Magnetic fields for individual qubits (dictionary)
    - J: Pairwise interactions between qubits (dictionary)
    - offset: Updated offset value
    """
    n_qubits = len(Q)  # Number of qubits (variables) in QUBO matrix
    h = defaultdict(int)  # Magnetic fields for individual qubits
    J = defaultdict(int)  # Pairwise interactions between qubits
    
    # Loop over each qubit in the QUBO matrix
    for i in range(n_qubits):
        # Update magnetic field for qubit i based on diagonal element of Q
        h[(i,)] -= Q[i, i] / 2
        # Update offset based on diagonal element of Q
        offset += Q[i, i] / 2
        
        # Loop over other qubits to calculate pairwise interactions
        for j in range(i + 1, n_qubits):
            # Update pairwise interaction strength (J) between qubits i and j
            J[(i, j)] += Q[i, j] / 4
            # Update magnetic fields for qubits i and j based on their interactions in Q
            h[(i,)] -= Q[i, j] / 4
            h[(j,)] -= Q[i, j] / 4
            # Update offset based on interaction strength between qubits i and j
            offset += Q[i, j] / 4
    
    # Return magnetic fields, pairwise interactions, and updated offset
    return h, J, offset

def energy_Ising(z, h, J, offset):
    """
    Calculate the energy of an Ising model given spin configurations.

    Parameters:
    - z: Spin configurations for each qubit (string or list)
    - h: Magnetic fields for individual qubits (dictionary)
    - J: Pairwise interactions between qubits (dictionary)
    - offset: Constant offset value

    Returns:
    - energy: Total energy of the Ising model
    """
    # If z is a string, convert to spins: 0 → +1, 1 → -1
    if isinstance(z, str):
        z = [(1 if int(i) == 0 else -1) for i in z]
    
    # Initialize energy with offset term
    energy = offset
    
    # Loop over magnetic fields (h) for each qubit and update energy
    for k, v in h.items():
        energy += v * z[k[0]]
    
    # Loop over pairwise interactions (J) between qubits and update energy
    for k, v in J.items():
        energy += v * z[k[0]] * z[k[1]]
    
    # Return total energy of the Ising model
    return energy

def samples_dict(counts, n_items):
    """
    Convert measurement results to a dictionary with bitstrings.
    
    Parameters:
    - counts: Measurement results from quantum circuit
    - n_items: Number of bits to keep
    
    Returns:
    - results: Dictionary with bitstrings and their frequencies
    """
    results = defaultdict(int)
    for bitstring, count in counts.items():
        # Trim bitstring to first n_items bits
        results[bitstring[:n_items]] = count
    return results

def qaoa_circuit(gammas, betas, h, J, num_qubits):
    """
    Parameters:
    - gammas: Parameters for cost Hamiltonian (list)
    - betas: Parameters for mixer Hamiltonian (list)
    - h: Magnetic fields (dictionary)
    - J: Pairwise interactions (dictionary)
    - num_qubits: Number of qubits
    """
    qc = QuantumCircuit(num_qubits, num_qubits)
    
    # Normalization - find maximum value for scaling
    h_vals = list(h.values())
    J_vals = list(J.values())
    wmax = max(np.max(np.abs(h_vals)), np.max(np.abs(J_vals))) if (h_vals or J_vals) else 1.0
    if wmax == 0:
        wmax = 1.0
    
    P = len(gammas)  # Number of QAOA layers
    
    # Initial layer: Hadamard gate on all qubits
    # Creates uniform superposition of all states
    for i in range(num_qubits):
        qc.h(i)
    
    # P layers of QAOA
    for layer in range(P):
        γ = gammas[layer]  # Parameter for cost Hamiltonian
        β = betas[layer]   # Parameter for mixer Hamiltonian
        
        # ---------- COST HAMILTONIAN ----------
        # Implements the cost function of the problem
        
        # Single-qubit terms: h_i * Z_i
        for (i,), v in h.items():
            if abs(v) > 1e-15:  # Ignore negligible values
                qc.rz(2 * γ * v / wmax, i)
        
        # Two-qubit terms: J_ij * Z_i * Z_j
        for (i, j), vij in J.items():
            if abs(vij) > 1e-15:  # Ignore negligible values
                qc.cx(i, j)  # CNOT gate
                qc.rz(2 * γ * vij / wmax, j)  # Rotation around Z
                qc.cx(i, j)  # CNOT gate
        
        # ---------- MIXER HAMILTONIAN (X mixer) ----------
        # Allows transitions between states
        for i in range(num_qubits):
            qc.rx(-2 * β, i)  # Rotation around X axis
    
    # Measure all qubits at the end
    qc.measure_all()
    
    return qc

Optimized Circuit - Grouping Commuting Gates

Theoretical Foundation:
The cost Hamiltonian in QAOA has the form: $H_C = \sum_i h_i Z_i + \sum_{i<j} J_{ij} Z_i Z_j$

Since all Pauli-Z operators commute with each other ($[Z_i, Z_j] = 0$), we can apply the corresponding rotations in any order. This allows strategic grouping for better performance.

---

Mathematical Derivation: QAOA Cost and Mixer Hamiltonians

---

Single-Qubit Terms: $h_i Z_i$

Implementation

We need to implement the time evolution operator for a single-qubit term in the cost Hamiltonian: $$U_i(\gamma) = e^{-i\gamma h_i Z_i}$$

Pauli-Z Operator in Computational Basis

The Pauli-Z operator acts on the computational basis states $|0\rangle$ and $|1\rangle$:

$$Z = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}$$

Action on basis states: $$Z|0\rangle = +1 \cdot |0\rangle$$ $$Z|1\rangle = -1 \cdot |1\rangle$$

Matrix Exponential

The exponential of the Pauli-Z operator: $$e^{-i\theta Z} = \sum_{k=0}^{\infty} \frac{(-i\theta Z)^k}{k!}$$

Since $Z^2 = I$ (identity), we can separate even and odd terms: $$e^{-i\theta Z} = \sum_{k=0}^{\infty} \frac{(-i\theta)^{2k}}{(2k)!} I + \sum_{k=0}^{\infty} \frac{(-i\theta)^{2k+1}}{(2k+1)!} Z$$

$$e^{-i\theta Z} = \cos(\theta) \cdot I - i\sin(\theta) \cdot Z$$

Explicit Matrix Form

$$e^{-i\theta Z} = \cos(\theta) \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} - i\sin(\theta) \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}$$

$$e^{-i\theta Z} = \begin{pmatrix} \cos(\theta) - i\sin(\theta) & 0 \\ 0 & \cos(\theta) + i\sin(\theta) \end{pmatrix}$$

Using Euler's formula: $e^{-i\theta} = \cos(\theta) - i\sin(\theta)$

$$e^{-i\theta Z} = \begin{pmatrix} e^{-i\theta} & 0 \\ 0 & e^{i\theta} \end{pmatrix}$$

Qiskit RZ Gate Definition

The RZ gate in Qiskit is defined as: $$RZ(\phi) = e^{-i\frac{\phi}{2} Z} = \begin{pmatrix} e^{-i\phi/2} & 0 \\ 0 & e^{i\phi/2} \end{pmatrix}$$

Finaly

We need: $e^{-i\gamma h_i Z}$

Comparing with RZ definition: $$e^{-i\gamma h_i Z} = e^{-i\frac{2\gamma h_i}{2} Z} = RZ(2\gamma h_i)$$

Implementation

qc.rz(2 * gamma * h_i / wmax, i)

Where: $\theta = 2\gamma h_i / w_{max}$

---

Two-Qubit Terms: $J_{ij} Z_i Z_j$ - Mathematical Derivation & Commutativity

We need to implement the time evolution operator for pairwise interaction terms in the cost Hamiltonian: $$U_{ij}(\gamma) = e^{-i\gamma J_{ij} Z_i Z_j}$$

Tensor Product of Pauli-Z Operators

For two qubits $i$ and $j$, the operator $Z_i Z_j$ is a tensor product:

$$Z_i \otimes Z_j = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \otimes \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}$$

$$Z_i \otimes Z_j = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & -1 & 0 & 0 \\ 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix}$$

Action on computational basis: $$Z_i Z_j |00\rangle = (+1)(+1)|00\rangle = +|00\rangle$$ $$Z_i Z_j |01\rangle = (+1)(-1)|01\rangle = -|01\rangle$$ $$Z_i Z_j |10\rangle = (-1)(+1)|10\rangle = -|10\rangle$$ $$Z_i Z_j |11\rangle = (-1)(-1)|11\rangle = +|11\rangle$$

Key property: $(Z_i Z_j)^2 = I$ (identity)

Matrix Exponential of $Z_i Z_j$

Using the same technique as for single-qubit terms:

$$e^{-i\theta Z_i Z_j} = \cos(\theta) \cdot I - i\sin(\theta) \cdot Z_i Z_j$$

$$e^{-i\theta Z_i Z_j} = \begin{pmatrix} e^{-i\theta} & 0 & 0 & 0 \\ 0 & e^{i\theta} & 0 & 0 \\ 0 & 0 & e^{i\theta} & 0 \\ 0 & 0 & 0 & e^{-i\theta} \end{pmatrix}$$

Qiskit RZZ Gate Definition

The RZZ gate in Qiskit is defined as: $$RZZ(\phi) = e^{-i\frac{\phi}{2} Z \otimes Z} = \begin{pmatrix} e^{-i\phi/2} & 0 & 0 & 0 \\ 0 & e^{i\phi/2} & 0 & 0 \\ 0 & 0 & e^{i\phi/2} & 0 \\ 0 & 0 & 0 & e^{-i\phi/2} \end{pmatrix}$$

Comparing

We need: $e^{-i\gamma J_{ij} Z_i Z_j}$

Comparing with RZZ definition: $$e^{-i\gamma J_{ij} Z_i Z_j} = e^{-i\frac{2\gamma J_{ij}}{2} Z_i Z_j} = RZZ(2\gamma J_{ij})$$

Implementation

qc.rzz(2 * gamma * J_ij / wmax, i, j)

Where: $\theta = 2\gamma J_{ij} / w_{max}$

---

Commutativity of Multiple RZZ Gates

Fundamental Theorem

Two RZZ gates commute if and only if they act on disjoint sets of qubits.

Case 1: Disjoint Qubits (COMMUTE )

Gates on different qubit pairs: $$RZZ(\theta_1)_{0,1} \cdot RZZ(\theta_2)_{2,3} = RZZ(\theta_2)_{2,3} \cdot RZZ(\theta_1)_{0,1}$$

Proof:

• $RZZ(\theta_1)_{0,1}$ acts as: $e^{-i\frac{\theta_1}{2} Z_0 Z_1} \otimes I_2 \otimes I_3$
• $RZZ(\theta_2)_{2,3}$ acts as: $I_0 \otimes I_1 \otimes e^{-i\frac{\theta_2}{2} Z_2 Z_3}$
• Since they act on different tensor product spaces: $[A \otimes I, I \otimes B] = 0$

Example in code:

# These two RZZ gates COMMUTE - can be applied in any order
qc.rzz(0.5, 0, 1)  # Acts on qubits 0 and 1
qc.rzz(0.3, 2, 3)  # Acts on qubits 2 and 3
# Order doesn't matter! Hardware can execute them in parallel 

Case 2: One Shared Qubit (DO NOT COMMUTE )

Gates sharing exactly one qubit: $$RZZ(\theta_1)_{0,1} \cdot RZZ(\theta_2)_{1,2} \neq RZZ(\theta_2)_{1,2} \cdot RZZ(\theta_1)_{0,1}$$

Why?

• $Z_0 Z_1$ and $Z_1 Z_2$ both act on qubit 1
• The operators $Z_0 Z_1$ and $Z_1 Z_2$ do NOT commute in general
• Mathematical proof: $[Z_0 Z_1, Z_1 Z_2] = [Z_0, Z_2] \cdot Z_1^2 = [Z_0, Z_2] \neq 0$

Example in code:

# These RZZ gates DO NOT COMMUTE - order matters!
qc.rzz(0.5, 0, 1)  # Acts on qubits 0 and 1
qc.rzz(0.3, 1, 2)  # Acts on qubits 1 and 2 (shares qubit 1)
# Changing order changes the result 

Case 3: Same Qubit Pair (COMMUTE but MERGE)

Gates on identical qubit pairs: $$RZZ(\theta_1)_{0,1} \cdot RZZ(\theta_2)_{0,1} = RZZ(\theta_1 + \theta_2)_{0,1}$$

Proof: $$e^{-i\frac{\theta_1}{2} Z_0 Z_1} \cdot e^{-i\frac{\theta_2}{2} Z_0 Z_1} = e^{-i\frac{\theta_1 + \theta_2}{2} Z_0 Z_1}$$

Example in code:

# These can be merged into a single RZZ gate
qc.rzz(0.5, 0, 1)
qc.rzz(0.3, 0, 1)
# Equivalent to:
qc.rzz(0.8, 0, 1)  # Single gate is more efficient!

Case 4: All RZ Gates (ALWAYS COMMUTE )

Single-qubit RZ gates on different qubits: $$RZ(\theta_1)_i \cdot RZ(\theta_2)_j = RZ(\theta_2)_j \cdot RZ(\theta_1)_i \quad \forall i \neq j$$

Proof:

• $RZ(\theta_1)_i$ acts as: $e^{-i\frac{\theta_1}{2} Z_i} \otimes I_j$
• $RZ(\theta_2)_j$ acts as: $I_i \otimes e^{-i\frac{\theta_2}{2} Z_j}$
• Tensor products of operators on different spaces always commute

Even on the same qubit: $$RZ(\theta_1)_i \cdot RZ(\theta_2)_i = RZ(\theta_1 + \theta_2)_i$$

Commutativity Summary Table

Gate Pair	Qubits	Commute?	Optimization
RZ(θ₁, i) & RZ(θ₂, j)	Different	Yes	Parallel execution
RZ(θ₁, i) & RZ(θ₂, i)	Same	Yes	Merge into one
RZZ(θ₁, i,j) & RZZ(θ₂, k,l)	Disjoint (i,j) ≠ (k,l)	Yes	Parallel execution
RZZ(θ₁, i,j) & RZZ(θ₂, i,j)	Same pair	Yes	Merge into one
RZZ(θ₁, 0,1) & RZZ(θ₂, 1,2)	Share one qubit	No	Order matters
RZ(θ, i) & RZZ(φ, j,k)	i ∉ {j,k}	Yes	Parallel execution
RZ(θ, i) & RZZ(φ, i,j)	i ∈ {i,j}	Yes	Can reorder

Key insight: RZ and RZZ gates commute with each other as long as they act on completely disjoint qubits OR if one is RZ acting on a single qubit.

Practical Application in QAOA

Why Grouping Matters

#  UNOPTIMIZED: Mixed order
for (i,), h_i in h.items():
    qc.rz(2 * gamma * h_i / wmax, i)
for (i, j), J_ij in J.items():
    qc.rzz(2 * gamma * J_ij / wmax, i, j)
for (i,), h_i in h.items():  # More RZ gates later
    qc.rz(2 * gamma * h_i / wmax, i)

Problem: Transpiler sees interleaved gates and cannot optimize routing effectively.

#  OPTIMIZED: Grouped by type
# Phase 1: All single-qubit rotations
for (i,), h_i in h.items():
    qc.rz(2 * gamma * h_i / wmax, i)

# Phase 2: All two-qubit interactions
for (i, j), J_ij in J.items():
    qc.rzz(2 * gamma * J_ij / wmax, i, j)

Benefits:

1. All RZ gates can run in parallel (act on different qubits)
2. Transpiler sees complete interaction graph for optimal routing
3. Fewer SWAP gates needed on hardware with limited connectivity
4. Reduced circuit depth due to parallelization opportunities

Visual Example: 4-Qubit QAOA Cost Layer

QUBO problem with interactions:

• $h = \{(0,): 0.5, (1,): 0.3, (2,): 0.8, (3,): 0.2\}$
• $J = \{(0,1): 0.4, (1,2): 0.6, (2,3): 0.3\}$

Unoptimized circuit:

q0: ─RZ(0.5)─●───────────────────────────────────
             │                                   
q1: ─────────●─RZ(0.3)─●─────────────────────────
                       │                         
q2: ───────────────────●─RZ(0.8)─●───────────────
                                 │               
q3: ─────────────────────────────●─RZ(0.2)───────

Depth: ~7 gates (serial execution)

Optimized circuit (grouped & parallel):

q0: ─RZ(0.5)─┐          ●────────────────────────
             │          │                        
q1: ─RZ(0.3)─┤ Parallel ●────●───────────────────
             │ RZ gates      │                   
q2: ─RZ(0.8)─┤               ●─────●─────────────
             │                     │             
q3: ─RZ(0.2)─┘                     ●─────────────
     
     └─ RZ ─┘  └────── RZZ ──────────────┘

Depth: ~3 layers (parallel execution)
Reduction: ~57% 

Mathematical Proof of Overall Commutativity

Theorem: The entire QAOA cost Hamiltonian layer can be decomposed into commuting parts.

Cost Hamiltonian: $$H_C = \sum_i h_i Z_i + \sum_{i<j} J_{ij} Z_i Z_j$$

Time evolution: $$U_C(\gamma) = e^{-i\gamma H_C} = e^{-i\gamma \left(\sum_i h_i Z_i + \sum_{i<j} J_{ij} Z_i Z_j\right)}$$

Key observation: All terms in $H_C$ are products of Pauli-Z operators. Since: $$[Z_i, Z_j] = 0 \quad \forall i, j$$

We can write: $$U_C(\gamma) = \prod_i e^{-i\gamma h_i Z_i} \prod_{i<j} e^{-i\gamma J_{ij} Z_i Z_j}$$

This allows us to:

1. Apply all single-qubit rotations in any order (or simultaneously)
2. Apply all two-qubit interactions in any valid order
3. Separate RZ and RZZ phases for better optimization

Key Takeaways

1. Commutativity rule: Gates commute if they act on disjoint qubit sets
2. Z-operators always commute: $[Z_i, Z_j] = 0$ for all $i, j$
3. Grouping enables parallelism: Hardware can execute commuting gates simultaneously
4. Transpiler optimization: Seeing all gates of one type together allows better routing
5. Circuit depth reduction: 15-30% typical improvement from proper grouping

---

Mixer Hamiltonian: $H_M = \sum_i X_i$

Implementation

We need to implement the mixer Hamiltonian evolution: $$U_M(\beta) = e^{-i\beta H_M} = e^{-i\beta \sum_i X_i}$$

Since $X$ operators on different qubits commute, we can write: $$U_M(\beta) = \prod_i e^{-i\beta X_i}$$

Pauli-X Operator in Computational Basis

The Pauli-X operator (bit-flip):

$$X = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}$$

Action on basis states: $$X|0\rangle = |1\rangle$$ $$X|1\rangle = |0\rangle$$

Matrix Exponential

Following similar derivation as for Z, since $X^2 = I$:

$$e^{-i\theta X} = \cos(\theta) \cdot I - i\sin(\theta) \cdot X$$

$$e^{-i\theta X} = \begin{pmatrix} \cos(\theta) & -i\sin(\theta) \\ -i\sin(\theta) & \cos(\theta) \end{pmatrix}$$

Qiskit RX Gate Definition

$$RX(\phi) = e^{-i\frac{\phi}{2} X} = \begin{pmatrix} \cos(\phi/2) & -i\sin(\phi/2) \\ -i\sin(\phi/2) & \cos(\phi/2) \end{pmatrix}$$

Finaly

We need: $e^{-i\beta X_i}$

Comparing with RX definition: $$e^{-i\beta X_i} = e^{-i\frac{2\beta}{2} X} = RX(2\beta)$$

---

Sign Convention: Why the Minus Sign?

There are two common conventions for defining the mixer Hamiltonian:

Convention A: $H_M = +\sum_i X_i$

$$U_M = e^{-i\beta H_M} = e^{-i\beta \sum_i X_i}$$ Implementation: qc.rx(2 * beta, i)

Convention B: $H_M = -\sum_i X_i$

$$U_M = e^{-i\beta H_M} = e^{-i\beta \cdot (-\sum_i X_i)} = e^{+i\beta \sum_i X_i}$$ Implementation: qc.rx(-2 * beta, i)

Which Convention to Use?

Both are valid! The sign is absorbed into the optimization of $\beta$ parameters:

• If you use Convention A with $\beta > 0$, the algorithm finds optimal positive $\beta$
• If you use Convention B with $\beta > 0$, the algorithm finds the same solution (since $\beta$ can be negative during optimization)

Common practice:

• This makes $\beta$ positive values correspond to "mixing" in the intuitive direction

Implementation (Convention B)

qc.rx(-2 * beta, i)

Physical interpretation: For $\beta > 0$, this creates superpositions that help the quantum state explore the solution space.

Summary Table

Component	Hamiltonian Term	Evolution Operator	Implementation
Cost (linear)	$h_i Z_i$	$e^{-i\gamma h_i Z_i}$	RZ(2γh_i)
Cost (quadratic)	$J_{ij} Z_i Z_j$	$e^{-i\gamma J_{ij} Z_i Z_j}$	CNOT-RZ(2γJ_ij)-CNOT
Mixer	$-\sum_i X_i$	$e^{+i\beta X_i}$	RX(-2β)

---

Single-Qubit Rotations (RZ Gates)

# All linear terms at once
for (i,), v in h.items():
    if abs(v) > 1e-15:
        qc.rz(2 * gamma * v / wmax, i)

What it represents:

• Linear terms $h_i Z_i$ from the original QUBO problem
• E.g., weights of individual variables in the optimization problem

Why group them:

• Parallel execution: RZ on different qubits can be applied simultaneously by hardware
• Rotation merging: Multiple RZ on the same qubit merge: $R_Z(\theta_1) \cdot R_Z(\theta_2) = R_Z(\theta_1 + \theta_2)$
• Simpler analysis: Transpiler sees all single-qubit operations together

Example:

q0: ─RZ(0.5)─    q1: ─RZ(0.3)─    q2: ─RZ(0.8)─
     ↓                ↓                ↓
All three can run simultaneously! 

---

Two-Qubit Interactions (RZZ Gates)

# All interaction terms at once
for (i, j), vij in J.items():
    if abs(vij) > 1e-15:
        qc.rzz(2 * gamma * vij / wmax, i, j)

What it represents:

• Quadratic terms $J_{ij} Z_i Z_j$ from the QUBO matrix
• Interactions between pairs of variables (e.g., penalties or preferences)

Why group them:

• Global routing optimization: Transpiler sees all connections at once → better mapping to chip topology
• Minimize SWAP gates: Can plan more efficient paths between non-adjacent qubits
• Depth reduction: Some RZZ interactions can run in parallel if they don't share qubits

Example structure:

Before grouping:        After grouping + transpilation:
q0: ─RZ─●─────RZ─       q0: ─RZ─RZZ─────────
        │                      │    
q1: ────X─RZ──X─        q1: ─RZ─RZZ─RZZ─────
              │                     │
q2: ──────────X─        q2: ─RZ────RZZ──────

Depth: ~8               Depth: ~4 

---

Final Improvements:

• Circuit depth: ↓ 15-30%
• Transpilation time: ↓ 20-40%
• Number of SWAP gates: ↓ 10-25%

#MY OPTIMIZED CIRCUIT
def qaoa_circuit_commuting_gates(gammas, betas, h, J, num_qubits, mixer_type='x', use_barriers=True):
    """
    Parameters:
    - gammas: Parameters for cost Hamiltonian (list)
    - betas: Parameters for mixer Hamiltonian (list)
    - h: Magnetic fields (dictionary)
    - J: Pairwise interactions (dictionary)
    - num_qubits: Number of qubits
    """
    qc = QuantumCircuit(num_qubits, num_qubits)

    # Normalization - find maximum value for scaling
    h_vals = list(h.values())
    J_vals = list(J.values())
    wmax = max(np.max(np.abs(h_vals)), np.max(
        np.abs(J_vals))) if (h_vals or J_vals) else 1.0
    if wmax == 0:
        wmax = 1.0

    P = len(gammas)  # Number of QAOA layers

    # Initial layer: Hadamard gate on all qubits
    # Creates uniform superposition of all states
    for i in range(num_qubits):
        qc.h(i)

    # P layers of QAOA
    for layer in range(P):
        γ = gammas[layer]  # Parameter for cost Hamiltonian
        β = betas[layer]  # Parameter for mixer Hamiltonian

        # ---------- COST HAMILTONIAN ----------
        # Implements the cost function of the problem
        for (i,), v in h.items():
            if abs(v) > 1e-15:
                qc.rz(2 * γ * v / wmax, i)

        # Two-qubit terms (using RZZ for efficiency)
        for (i, j), vij in J.items():
            if abs(vij) > 1e-15:
                qc.rzz(2 * γ * vij / wmax, i, j)

        if use_barriers:
            qc.barrier()

        # ---------- MIXER HAMILTONIAN (X mixer) ----------
        if mixer_type == 'x':
            #old
            for i in range(num_qubits):
                qc.rx(-2 * β, i)

        elif mixer_type == 'xy':
            #the best
            for i in range(num_qubits - 1):
                qc.rxx(-2 * β, i, i + 1)
                qc.ryy(-2 * β, i, i + 1)

        if use_barriers:
            qc.barrier()

    # Measure all qubits at the end
    qc.measure_all()

    return qc

Solution using Quantum Algorithm QAOA

try:
    # Convert QuadraticProgram to QUBO matrix
    from qiskit_optimization.converters import QuadraticProgramToQubo
    qp2qubo = QuadraticProgramToQubo()
    qubo = qp2qubo.convert(qp)

    # Get Q matrix and offset
    Q_matrix = np.zeros((n_assets, n_assets))
    offset = qubo.objective.constant

    # Linear terms (diagonal)
    for var, coef in qubo.objective.linear.to_dict().items():
        # Handle both int and string variable types
        if isinstance(var, int):
            idx = var
        elif isinstance(var, str):
            idx = int(var.replace('x_', ''))
        else:
            idx = int(var)
        Q_matrix[idx, idx] = coef

    # Quadratic terms (off-diagonal)
    for (var1, var2), coef in qubo.objective.quadratic.to_dict().items():
        # Handle both cases for var1 and var2
        if isinstance(var1, int):
            idx1 = var1
        elif isinstance(var1, str):
            idx1 = int(var1.replace('x_', ''))
        else:
            idx1 = int(var1)

        if isinstance(var2, int):
            idx2 = var2
        elif isinstance(var2, str):
            idx2 = int(var2.replace('x_', ''))
        else:
            idx2 = int(var2)

        Q_matrix[idx1, idx2] = coef
        Q_matrix[idx2, idx1] = coef  # Symmetric matrix

    print(f"Problem size: {n_assets} qubits")
    print(f"Q matrix shape: {Q_matrix.shape}")
    print(f"Initial offset: {offset:.2f}")

    # Convert to Ising (dictionary format)
    h, J, offset = from_Q_to_Ising(Q_matrix, offset)
    n_qubits = n_assets

    print(f"Ising offset: {offset:.2f}")
    print(f"h terms: {len(h)}")
    print(f"J terms: {len(J)}")

    # Fixed annealing schedule
    betas = np.linspace(0, 1, P_LAYERS)[::-1]
    gammas = np.linspace(0, 1, P_LAYERS)

    print(f"\nAnnealing schedule (P={P_LAYERS}):")
    print(f"  Betas:  {betas}")
    print(f"  Gammas: {gammas}")

    # Create QAOA circuit
    print("\nBuilding QAOA circuit...")
    #HERE IS MY OPTIMIZED CIRCUIT
    #qc = qaoa_circuit(gammas, betas, h, J, n_qubits)
    qc = qaoa_circuit_commuting_gates(gammas, betas, h, J, n_qubits, 'xy', True)

    print(f"Circuit properties:")
    print(f"  Depth: {qc.depth()}")
    print(f"  Gates: {qc.size()}")

    # Run circuit
    print(f"\nRunning QAOA circuit ({shots} shots)...")
    sampler = StatevectorSampler()
    job = sampler.run([qc], shots=shots)
    result = job.result()
    counts = result[0].data.meas.get_counts()

    # Convert to samples dictionary
    samples = samples_dict(counts, n_qubits)

    print(f"Total unique bitstrings: {len(samples)}")

    # Evaluate all samples
    valid_solutions = {}

    for bitstring, count in samples.items():
        x = [int(b) for b in bitstring]

        # Check budget constraint
        if sum(x) != budget:
            continue

        # Calculate energy
        energy = energy_Ising(bitstring, h, J, offset)

        # Calculate portfolio metrics
        selected = [i for i, bit in enumerate(x) if bit == 1]
        total_return = sum(expected_returns[i] for i in selected)
        total_risk = np.sqrt(
            sum(covariance_matrix[i][j] for i in selected for j in selected))
        sharpe = total_return / total_risk if total_risk > 0 else 0

        valid_solutions[bitstring] = {
            'count': count,
            'energy': energy,
            'return': total_return,
            'risk': total_risk,
            'sharpe': sharpe,
            'selected': selected
        }

    valid_count = sum(s['count'] for s in valid_solutions.values())
    print(f"Valid solutions (budget={budget}): {len(valid_solutions)}")
    print(f"Valid shots: {valid_count}/{shots} ({valid_count/shots*100:.1f}%)")

    if len(valid_solutions) == 0:
        raise Exception(
            "No valid solution found (budget constraint not satisfied)")

    # Find best solution (minimum energy)
    best_bitstring = min(valid_solutions.keys(),
                         key=lambda k: valid_solutions[k]['energy'])
    best_solution = valid_solutions[best_bitstring]

    print("\nBest solution found:")
    print(f"  Bitstring: {best_bitstring}")
    print(
        f"  Frequency: {best_solution['count']}/{shots} ({best_solution['count']/shots*100:.1f}%)"
    )
    print(f"  Energy: {best_solution['energy']:.2f}")

    # Create compatible qaoa_result object
    class QAOAResult:

        def __init__(self, x, fval):
            self.x = np.array(x, dtype=float)
            self.fval = fval

    qaoa_result_x = [float(int(b)) for b in best_bitstring]
    qaoa_result = QAOAResult(qaoa_result_x, best_solution['energy'])

    print(f"Objective function value: {qaoa_result.fval:.2f}")
    print(f"\nDebug - x vector: {qaoa_result.x}")

    selected_qaoa = []
    print("\nSelected stocks:")
    print("-" * 70)
    for i in range(n_assets):
        value = qaoa_result.x[i]
        if value > 0.5:
            selected_qaoa.append(i)
            print(
                f"  ✓ x_{i} = {value:.3f} | {asset_names[i]:<15s} | "
                f"Return: {expected_returns[i]:5.1f}% | Risk: {volatility[i]:5.1f}%"
            )
        else:
            print(f"  ✗ x_{i} = {value:.3f} | {asset_names[i]:<15s}")

    print(f"\nTotal stocks selected: {len(selected_qaoa)} (target: {budget})")

    if len(selected_qaoa) > 0:
        total_return_qaoa = sum(expected_returns[i] for i in selected_qaoa)
        total_risk_qaoa = np.sqrt(
            sum(covariance_matrix[i][j]
                for i in selected_qaoa
                for j in selected_qaoa))
        print(f"Total expected return: {total_return_qaoa:.2f}%")
        print(f"Total portfolio risk: {total_risk_qaoa:.2f}")
        print(f"Sharpe ratio: {total_return_qaoa/total_risk_qaoa:.3f}")
    else:
        total_return_qaoa = 0
        total_risk_qaoa = 0

    # Top 5 solutions
    print("\nTop 5 solutions by energy:")
    print("-" * 80)

    sorted_solutions = sorted(valid_solutions.items(),
                              key=lambda x: x[1]['energy'])[:5]

    for i, (bitstring, sol) in enumerate(sorted_solutions, 1):
        print(f"{i}. Energy: {sol['energy']:7.2f} | "
              f"Sharpe: {sol['sharpe']:.3f} | "
              f"Freq: {sol['count']:4d}/{shots} | "
              f"String: {bitstring}")

    qaoa_success = True

except Exception as e:
    print(f"\nERROR: {e}")
    import traceback
    traceback.print_exc()
    selected_qaoa = []
    qaoa_result = exact_result
    total_return_qaoa = total_return_exact
    total_risk_qaoa = total_risk_exact
    qaoa_success = False

Problem size: 10 qubits
Q matrix shape: (10, 10)
Initial offset: 3750.00
Ising offset: 121.72
h terms: 10
J terms: 45

Annealing schedule (P=5):
  Betas:  [1.   0.75 0.5  0.25 0.  ]
  Gammas: [0.   0.25 0.5  0.75 1.  ]

Building QAOA circuit...
Circuit properties:
  Depth: 182
  Gates: 385

Running QAOA circuit (2000 shots)...
Total unique bitstrings: 574
Valid solutions (budget=5): 129
Valid shots: 489/2000 (24.4%)

Best solution found:
  Bitstring: 1100100011
  Frequency: 2/2000 (0.1%)
  Energy: -436.98
Objective function value: -436.98

Debug - x vector: [1. 1. 0. 0. 1. 0. 0. 0. 1. 1.]

Selected stocks:
----------------------------------------------------------------------
  ✓ x_0 = 1.000 | Apple (Tech)    | Return:  24.3% | Risk:  27.5%
  ✓ x_1 = 1.000 | Microsoft (Tech) | Return:  31.9% | Risk:  24.4%
  ✗ x_2 = 0.000 | ExxonMobil (Energie)
  ✗ x_3 = 0.000 | Chevron (Energie)
  ✓ x_4 = 1.000 | JPMorgan (Finance) | Return:  34.7% | Risk:  22.9%
  ✗ x_5 = 0.000 | Bank of America (Finance)
  ✗ x_6 = 0.000 | Johnson & Johnson (Healthcare)
  ✗ x_7 = 0.000 | Pfizer (Healthcare)
  ✓ x_8 = 1.000 | Walmart (Retail) | Return:  30.7% | Risk:  19.9%
  ✓ x_9 = 1.000 | Amazon (E-commerce) | Return:  29.1% | Risk:  33.2%

Total stocks selected: 5 (target: 5)
Total expected return: 150.75%
Total portfolio risk: 92.00
Sharpe ratio: 1.639

Top 5 solutions by energy:
--------------------------------------------------------------------------------
1. Energy: -436.98 | Sharpe: 1.639 | Freq:    2/2000 | String: 1100100011
2. Energy: -421.63 | Sharpe: 1.617 | Freq:    6/2000 | String: 0100110011
3. Energy: -408.71 | Sharpe: 1.654 | Freq:   23/2000 | String: 1100110010
4. Energy: -398.25 | Sharpe: 1.501 | Freq:    5/2000 | String: 1000110011
5. Energy: -398.07 | Sharpe: 1.398 | Freq:    8/2000 | String: 1100110001

Vizualization results

print(f"{'Objective function value':<35s} {exact_result.fval:>15.2f} {qaoa_result.fval:>15.2f}")
print(f"{'Total return (%)':<35s} {total_return_exact:>15.2f} {total_return_qaoa:>15.2f}")
print(f"{'Total risk':<35s} {total_risk_exact:>15.2f} {total_risk_qaoa:>15.2f}")
print(f"{'Number of selected stocks':<35s} {len(selected_exact):>15d} {len(selected_qaoa):>15d}")
print(f"{'Sharpe ratio (return/risk)':<35s} {total_return_exact/total_risk_exact:>15.3f} {total_return_qaoa/total_risk_qaoa:>15.3f}")

# Calculate difference
diff = abs(qaoa_result.fval - exact_result.fval)
diff_percent = (diff / abs(exact_result.fval)) * 100 if exact_result.fval != 0 else 0
print(f"{'Difference in objective function':<35s} {diff:>15.2f} ({diff_percent:.1f}%)")

if len(selected_exact) > 0 or len(selected_qaoa) > 0:

    pastel_blue = '#A8D8EA'
    pastel_pink = '#FFB6C1'
    pastel_green = '#C1E1C1'
    pastel_gray = '#D3D3D3'
    pastel_purple = '#E6B3E0'
    pastel_yellow = '#FFF9A8'

    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    fig.patch.set_facecolor('#FFFFFF')

    # Graph 1: Selected stocks
    ax1 = axes[0, 0]
    ax1.set_facecolor('#FFFFFF')
    x_pos = np.arange(n_assets)
    width = 0.35

    exact_binary = [1 if i in selected_exact else 0 for i in range(n_assets)]
    qaoa_binary = [1 if i in selected_qaoa else 0 for i in range(n_assets)]

    ax1.bar(x_pos - width / 2,
            exact_binary,
            width,
            label='Numpy',
            alpha=0.85,
            color=pastel_blue,
            edgecolor='#5A7D8F',
            linewidth=2)
    ax1.bar(x_pos + width / 2,
            qaoa_binary,
            width,
            label='QAOA',
            alpha=0.85,
            color=pastel_pink,
            edgecolor='#D17A8F',
            linewidth=2)

    ax1.set_xlabel('Stocks', fontweight='bold', fontsize=11)
    ax1.set_ylabel('Selected (1=yes, 0=no)', fontweight='bold', fontsize=11)
    ax1.set_title('Comparison of Selected Stocks',
                  fontweight='bold',
                  pad=15,
                  fontsize=13)
    ax1.set_xticks(x_pos)
    ax1.set_xticklabels(asset_names, rotation=45, ha='right', fontsize=9)
    ax1.legend(fontsize=10, framealpha=0.9)
    ax1.grid(axis='y', alpha=0.3, color='gray', linestyle='--')
    ax1.set_ylim([0, 1.3])

    # Graph 2: Return vs Risk
    ax2 = axes[0, 1]
    ax2.set_facecolor('#FFFFFF')

    if total_return_exact > 0 and total_return_qaoa > 0:
        ax2.scatter(total_risk_exact,
                    total_return_exact,
                    s=400,
                    label='numpy',
                    color=pastel_blue,
                    edgecolors='#5A7D8F',
                    linewidth=3,
                    zorder=5)
        ax2.scatter(total_risk_qaoa,
                    total_return_qaoa,
                    s=400,
                    label='QAOA',
                    color=pastel_pink,
                    edgecolors='#D17A8F',
                    linewidth=3,
                    zorder=5)

        ax2.annotate(
            f'numpy\n({total_return_exact:.1f}%, {total_risk_exact:.0f})',
            (total_risk_exact, total_return_exact),
            xytext=(0, 20),
            textcoords='offset points',
            fontsize=10,
            fontweight='bold',
            ha='center',
            bbox=dict(boxstyle='round,pad=0.5',
                      facecolor=pastel_blue,
                      alpha=0.7,
                      edgecolor='#5A7D8F',
                      linewidth=2))
        ax2.annotate(f'QAOA\n({total_return_qaoa:.1f}%, {total_risk_qaoa:.0f})',
                     (total_risk_qaoa, total_return_qaoa),
                     xytext=(0, -30),
                     textcoords='offset points',
                     fontsize=10,
                     fontweight='bold',
                     ha='center',
                     bbox=dict(boxstyle='round,pad=0.5',
                               facecolor=pastel_pink,
                               alpha=0.7,
                               edgecolor='#D17A8F',
                               linewidth=2))

    ax2.set_xlabel('Portfolio Risk', fontweight='bold', fontsize=11)
    ax2.set_ylabel('Expected Return (%)', fontweight='bold', fontsize=11)
    ax2.set_title('Return vs. Risk', fontweight='bold', pad=15, fontsize=13)
    ax2.legend(fontsize=10, framealpha=0.9)
    ax2.grid(True, alpha=0.3, color='gray', linestyle='--')

    # Graph 3: Contribution to return
    ax3 = axes[1, 0]
    ax3.set_facecolor('#FFFFFF')

    returns_exact = [
        expected_returns[i] if i in selected_exact else 0
        for i in range(n_assets)
    ]
    returns_qaoa = [
        expected_returns[i] if i in selected_qaoa else 0
        for i in range(n_assets)
    ]

    ax3.bar(x_pos - width / 2,
            returns_exact,
            width,
            label='numpy',
            alpha=0.85,
            color=pastel_blue,
            edgecolor='#5A7D8F',
            linewidth=2)
    ax3.bar(x_pos + width / 2,
            returns_qaoa,
            width,
            label='QAOA',
            alpha=0.85,
            color=pastel_pink,
            edgecolor='#D17A8F',
            linewidth=2)

    ax3.set_xlabel('Stocks', fontweight='bold', fontsize=11)
    ax3.set_ylabel('Contribution to Return (%)', fontweight='bold', fontsize=11)
    ax3.set_title('Contribution of Individual Stocks',
                  fontweight='bold',
                  pad=15,
                  fontsize=13)
    ax3.set_xticks(x_pos)
    ax3.set_xticklabels(asset_names, rotation=45, ha='right', fontsize=9)
    ax3.legend(fontsize=10, framealpha=0.9)
    ax3.grid(axis='y', alpha=0.3, color='gray', linestyle='--')

    # Graph 4: Stock characteristics
    ax4 = axes[1, 1]
    ax4.set_facecolor('#FFFFFF')

    # Pastel colors for points - green for selected, gray for unselected
    colors = [
        pastel_green if i in selected_exact else pastel_gray
        for i in range(n_assets)
    ]
    edge_colors = [
        '#6B9B6B' if i in selected_exact else '#A0A0A0' for i in range(n_assets)
    ]

    ax4.scatter(volatility,
                expected_returns,
                s=250,
                alpha=0.8,
                c=colors,
                edgecolors=edge_colors,
                linewidth=2.5)

    for i in range(n_assets):
        symbol = '✓' if i in selected_exact else '✗'
        color = '#2E7D32' if i in selected_exact else '#666666'
        ax4.annotate(f'{symbol} {asset_names[i]}',
                     (volatility[i], expected_returns[i]),
                     xytext=(5, 5),
                     textcoords='offset points',
                     fontsize=9,
                     fontweight='bold' if i in selected_exact else 'normal',
                     color=color)

    ax4.set_xlabel('Risk (Volatility %)', fontweight='bold', fontsize=11)
    ax4.set_ylabel('Expected Return (%)', fontweight='bold', fontsize=11)
    ax4.set_title('Stock Characteristics (✓ = selected)',
                  fontweight='bold',
                  pad=15,
                  fontsize=13)
    ax4.grid(True, alpha=0.3, color='gray', linestyle='--')

    # Add legend for graph 4
    from matplotlib.patches import Patch
    legend_elements = [
        Patch(facecolor=pastel_green,
              edgecolor='#6B9B6B',
              linewidth=2,
              label='Selected'),
        Patch(facecolor=pastel_gray,
              edgecolor='#A0A0A0',
              linewidth=2,
              label='Not Selected')
    ]
    ax4.legend(handles=legend_elements,
               fontsize=10,
               framealpha=0.9,
               loc='upper left')

    plt.tight_layout(pad=2.5)

    # Save
    output_filename = 'portfolio_results.png'
    try:
        plt.savefig(output_filename,
                    dpi=150,
                    bbox_inches='tight',
                    facecolor='#FAFAFA')
        import os
        print(f"\n Graphs saved to: {os.path.abspath(output_filename)}")
    except Exception as e:
        print(f"\n Cannot save graph: {e}")

    plt.show()

Objective function value                    -436.98         -436.98
Total return (%)                             150.75          150.75
Total risk                                    92.00           92.00
Number of selected stocks                         5               5
Sharpe ratio (return/risk)                    1.639           1.639
Difference in objective function               0.00 (0.0%)

 Graphs saved to: B:\python\myqiskit\Scripts\portfolio_results.png

Detailed solution analysis

Experimentation

Try modifying the following parameters and observe changes in the solution:

1. risk_factor / return_factor: Change the ratio between risk and return

   • Example: Set risk_factor=2.0, return_factor=0.5 → prefers safer stocks

2. budget: Change the number of stocks you want to select

   • Example: Set budget=2 → selects only the two best stocks

3. reps in QAOA: Increase the number of layers for potentially better results

   • Example: Set reps=5 → may find better solutions (but slower)

4. expected_returns / volatility: Modify the stock characteristics

   • Experiment with your own real market data

If nothing is found

1. Increase return_factor (e.g., to 5.0 or 10.0)
2. Decrease risk_factor (e.g., to 0.001)
3. Increase penalty (e.g., to 200.0)
4. Or change approach - use ConstrainedQuadraticProgram

import sys
print(f"\nPython Version: {sys.version}")
import numpy as np
print(f"numpy: {np.__version__}")
import matplotlib
print(f"matplotlib: {matplotlib.__version__}")
import pandas as pd
print(f"pandas: {pd.__version__}")
import qiskit
print(f"qiskit: {qiskit.__version__}")
import qiskit_optimization
print(f"qiskit-optimization: {qiskit_optimization.__version__}")
import qiskit_algorithms
print(f"qiskit-algorithms: {qiskit_algorithms.__version__}")
import qiskit_aer
print(f"qiskit-aer: {qiskit_aer.__version__}")
import yfinance as yf
print(f"yfinance: {yf.__version__}")

Python Version: 3.13.5 | packaged by Anaconda, Inc. | (main, Jun 12 2025, 16:37:03) [MSC v.1929 64 bit (AMD64)]
numpy: 2.3.3
matplotlib: 3.10.6
pandas: 2.3.2
qiskit: 2.2.1
qiskit-optimization: 0.7.0
qiskit-algorithms: 0.4.0
qiskit-aer: 0.17.2
yfinance: 0.2.66

References: QAOA Circuit Optimization with Commuting Gates

1. Circuit Compilation Methodologies for QAOA

Summary: Demonstrates that multi-qubit CPHASE gates in QAOA circuits are commutative, enabling gate reordering for parallel execution and reduced SWAP gate count. Presents four optimization methodologies that significantly reduce circuit depth and gate count by exploiting commutativity.

Link: IEEE Xplore

---

2. Optimized Qubit Routing for Commuting Gates via Integer Programming

Summary: Recent research (2025) proving NP-hardness of qubit routing with commuting gates and proposing an integer programming-based algorithm guaranteed to return optimal solutions with minimum SWAP gates. Specifically addresses QAOA compilation challenges.

Link: arXiv:2507.12199

---

3. A Structured Method for Compilation of QAOA Circuits

Summary: Proposes a structured compilation method ensuring linear depth for QAOA circuits on multi-dimensional quantum architectures. Demonstrates scalability with circuits up to 1024 qubits, achieving 3.8× speedup in depth and 17% reduction in gate count.

Link: arXiv:2112.06143

---

4. What is the QAOA? - OpenQAOA Documentation

Summary: Comprehensive explanation of QAOA fundamentals, including how commuting terms in cost and mixer Hamiltonians decompose into products of rotation gates. Provides clear mathematical derivations connecting Hamiltonians to quantum gates (RZ, RZZ, RX).

Link: OpenQAOA Docs

---

Key Takeaways

• Commutativity is fundamental: All Z-operators commute, allowing flexible gate ordering
• Grouping benefits: Transpiler can optimize routing and minimize SWAP gates when gates are grouped
• RZZ vs CNOT decomposition: Using native RZZ gates is more efficient than CNOT-RZ-CNOT sequences
• Hardware awareness: Optimization quality depends on target quantum hardware topology