Quantum Optimization Hands-On: Solve a Portfolio HUBO on Real Infrastructure

Educational disclaimer. Nothing here is complete or constitutes financial advice. Service names, endpoints and printed values are snapshots from the day this was written and will drift. Treat the live pages and the API references as the single source of truth.

This lab runs in four phases: build the model (Steps 1 and 2), submit it (Step 3), check and read the answer (Steps 4 and 5), then evaluate wheter a QPU chip could solve it (Step 6). Every task below carries collapsed hints: Hint 1 tells you where you could look, Hint 2 gives additional pieces of the answer. Before you open any hint or solution: what do you think is happening here? Pause a minute on your own, or discuss with a partner. The solutions reveal the complete information that solves the problem.

0Step 0: Set up the environment

uv init kipu-optimization-lab
cd kipu-optimization-lab
uv venv
uv add qhub-service

uv add matplotlib requests

npm install -g @quantum-hub/qhubctl

uv run python task1_hubo.py

1Step 1: State the problem in numbers

Formulation is the most important step in solving an optimization problem on a quantum computer. Its size (number of variables), density (non-zero terms), and order (linear, quadratic, cubic, and so on) drive the circuit depth: the amount of sequential operations, or "length", of a quantum program. Current Quantum Processing Units are constrained in size (number of qubits) and error budget (hardware layout, qubit connectivity and 2-qubit-gate errors). Circuit depth is a decisive factor in whether a program runs inside a given QPU's error budget.

In this tutorial we will slip into the role of a financial analyst. While trading expected return against the risk of holding correlated pairs you pick 4 assets out of 8 fictional S&P 500 stocks. This problem is small, but structured enough to show the behaviour and to teach you how to scale it.

Create a file with the given data:

tickers = ["AAPL", "MSFT", "JPM", "XOM", "JNJ", "PG", "NVDA", "CVX"]
sector  = ["Tech", "Tech", "Financials", "Energy", "Health", "Staples", "Tech", "Energy"]

# Expected annual return, in percent (illustrative)
returns = [18.0, 16.0, 12.0, 11.0, 8.0, 7.0, 22.0, 10.0]

# High pairwise correlations: holding both concentrates risk
correlation = {
    (0, 1): 0.70,  # AAPL - MSFT   (mega-cap tech)
    (0, 6): 0.65,  # AAPL - NVDA
    (1, 6): 0.68,  # MSFT - NVDA
    (3, 7): 0.88,  # XOM  - CVX     (energy)
    (4, 5): 0.55,  # JNJ  - PG      (defensives)
}

k = 4   # hold exactly four names

Kipu's quantum optimizers minimize an objective written as a flat dictionary or array of variable interactions: a HUBO, the higher-order generalization of a QUBO. Each binary variable x_i is 1 if asset i is held, 0 if not. Three relationships go into the objective in this step:

1. Reward return. A negative term -return_i for each asset, so the minimizer prefers holding high-return assets.
2. Penalize correlated pairs (risk). A positive quadratic (or higher) term on correlated pairs.
3. Enforce the count. A cardinality penalty (Σ x_i − k)² that is zero when exactly k assets are held.

Later we will assign weights for both risk and constraint penalties to balance their influence. Sometimes those values need to be rebalance.

# task1_hubo.py
risk_penalty = 10.0    # weight on correlation risk
card_penalty = 50.0    # must dominate return + risk, see the encoding rule below
n = len(returns)

hubo = {}
hubo["()"] = ...                              # constant
for i in range(n):
    hubo[f"({i},)"] = ...                     # linear
for i in range(n):
    for j in range(i + 1, n):
        rho = correlation.get((i, j), 0.0)    # 0.0 for pairs not listed
        hubo[f"({i}, {j})"] = ...             # pairwise

The loops run over every pair, because the count penalty applies to all.

hubo["()"]           = card_penalty * k**2
hubo[f"({i},)"]      = -returns[i] + card_penalty * (1 - 2 * k)
hubo[f"({i}, {j})"]  = risk_penalty * rho + card_penalty * 2

# task1_hubo.py
risk_penalty = 10.0   # weight on correlation risk
card_penalty = 50.0   # must dominate return + risk so "exactly k" always wins
n = len(returns)

hubo = {"()": card_penalty * k**2}                      # constant term

for i, r in enumerate(returns):
    hubo[f"({i},)"] = -r + card_penalty * (1 - 2 * k)   # reward return, count term

for i in range(n):
    for j in range(i + 1, n):
        rho = correlation.get((i, j), 0.0)
        hubo[f"({i}, {j})"] = risk_penalty * rho + card_penalty * 2

Encoding rule. Keep the cardinality penalty roughly 10× larger than the return-plus-risk range. Too small and the optimizer "cheats" by holding three or five names to dodge a correlation penalty; too large and the correlation-risk signal is drowned out. The constraint and the objective have to stay in balance.

2Step 2: Add the higher-order term

Everything in Step 1 is expressible as a QUBO: a pure means-and-pairwise-covariances model. During the 2008 financial crisis these pairwise models as we now know failed horribly. Means and pairwise covariances could not encode higher-order co-movement of multiple assets at all. Those models are computationally very expensive, and is the reason classical models are mostly pariwise-only. With quantum computers, we see the opportunity to encode the structure and eventually solve these problems buch better. Beware: the model is only as good as you formulate it!

BF-DCQO, the algorithm behind Kipu's optimizer, handles higher-order terms natively, with no quadratization and no auxiliary qubits, which is the subject of its own paper (arXiv:2409.04477). In dictionary terms it costs you exactly one more key.

hubo["(0, 1, 6)"] = 20.0

# Higher-order term: holding ALL of AAPL, MSFT, NVDA together is worse
# than the three pairwise penalties suggest (joint-crash co-movement).
# A QUBO cannot represent this; a HUBO takes it as one more dictionary key.
hubo["(0, 1, 6)"] = 20.0

Disclaimer: The pairwise-only model is what failed in 2008. While a cubic term lets you express joint co-movement,it does not make the model right. You chose the triple, the 20.0 and the correlations yourself, and the optimizer will faithfully return the best portfolio for the model you wrote. The model stays the binding limit.

from qiskit_addon_opt_mapper import OptimizationProblem
from qiskit_addon_opt_mapper.converters import OptimizationProblemToHubo

qp = OptimizationProblem(name="portfolio8")
x = [qp.binary_var(name=f"x{i}") for i in range(n)]

qp.minimize(
    linear={i: -returns[i] for i in range(n)},
    quadratic={(i, j): risk_penalty * rho for (i, j), rho in correlation.items()},
    higher_order={3: {(0, 1, 6): 20.0}},        # the cubic term, keyed by order
)
qp.linear_constraint(linear={i: 1 for i in range(n)}, sense="==", rhs=k, name="budget")

hubo_qp = OptimizationProblemToHubo(penalty=card_penalty).convert(qp)

3Step 3: Run it on Miray

Miray is Kipu's Advanced Quantum Optimizer on the Kipu Quantum Hub, running the BF-DCQO engine (arXiv:2405.13898), with a free simulator tier and a one-line switch to hardware versions later.

What you need before submitting

• SERVICE_ENDPOINT for Miray Advanced Quantum Optimizer - Simulator. You obtain this by subscribing to the service. See using a service.
• ACCESS_KEY_ID and SECRET_ACCESS_KEY, the credential pair a subscription issues. See Qhub Service SDK.
• Treat the keys like a password. An environment variable read with os.getenv is the least-bad option; keep them out of scripts you want to upload or share with others.
• You can delete and create new Application keys at any time in the Hub UI.

Reference pages for this step: Quickstart (account plus first run) and Service SDK (the client you are about to use).

client = HubServiceClient(service_endpoint=..., access_key_id=..., secret_access_key=...)
execution = client.run(request={"problem": hubo, "problem_type": ..., "shots": ..., ...})
execution.wait_for_final_state()
print(execution.result())

from qhub.service.client import HubServiceClient  # pip install qhub-service

client = HubServiceClient(
    service_endpoint=SERVICE_ENDPOINT,      # Miray Advanced Quantum Optimizer - Simulator
    access_key_id=ACCESS_KEY_ID,
    secret_access_key=SECRET_ACCESS_KEY,
)

execution = client.run(request={
    "problem": hubo,
    "problem_type": "binary",
    "shots": 2000,
    "num_iterations": 6,
    "num_greedy_passes": 1,
})
execution.wait_for_final_state()

response = execution.result()
print(response.result["cost"], response.result["mapped_solution"])
# Tip: increase shots and num_iterations for better results, or reduce for faster runtime.

cost = -56.5
mapped_solution = {0:1, 1:0, 2:1, 3:1, 4:0, 5:0, 6:1, 7:0}

from qiskit_ibm_catalog import QiskitFunctionsCatalog  # pip install qiskit-ibm-catalog

catalog = QiskitFunctionsCatalog(
    token=IBM_TOKEN, channel="ibm_quantum_platform", instance=INSTANCE_CRN,
)
iskay = catalog.load("kipu-quantum/iskay-quantum-optimizer")

job = iskay.run(
    problem=hubo,
    problem_type="binary",
    backend_name="ibm_fez",
    options={"shots": 2000, "num_iterations": 4},
)

result = job.result()
print(result["solution_info"]["cost"], result["solution"])
# Tip: increase shots and num_iterations for better results, or reduce for faster runtime.

4Step 4: Decode and verify

4a. Comparison strategies

At this complexity, the chepest comparison strategy is greedy: sort by expected return, take the top four. NVDA is the highest earner at 22%, AAPL next at 18%, MSFT at 16%, JPM at 12%. Done.

4b. Check the result

For a problem at this size we can confirm the global optimum by brute force. But beware: complexity grows quickly; brute force stops being an option well before a realistic portfolio size. Here you can compare results with other classical (heuristic) solvers.

solution = response.result["mapped_solution"]      # keys are string asset indices
selected = [tickers[int(i)] for i in solution if solution[i] == 1]
print(selected)
# -> ['AAPL', 'JPM', 'XOM', 'NVDA']

greedy = sorted(range(n), key=lambda i: -returns[i])[:k]
print([tickers[i] for i in sorted(greedy)])
# -> ['AAPL', 'MSFT', 'JPM', 'NVDA']

Verification, given

energy(x) recomputes your own objective in plain Python, and the brute force enumerates all 256 assignments.

from itertools import product

def energy(x):
    e = card_penalty * k**2
    for i in range(n):
        e += (-returns[i] + card_penalty * (1 - 2 * k)) * x[i]
    for i in range(n):
        for j in range(i + 1, n):
            rho = correlation.get((i, j), 0.0)
            e += (risk_penalty * rho + card_penalty * 2) * x[i] * x[j]
    e += 20.0 * x[0] * x[1] * x[6]   # the three-body tech-triple term
    return e

x = tuple(solution[str(i)] for i in range(n))

assert sum(x) == k                                              # exactly k names held
assert abs(energy(x) - response.result["cost"]) < 1e-3          # cost is reproducible
best = min(product([0, 1], repeat=n), key=energy)               # brute force, 256 assignments
assert energy(best) == energy(x)                                # it is the global optimum
print("Verified: global optimum, exactly k held, energy consistent.")

Brute force arrives at the same energy −56.5 matching the quantum solver. It is only available to you because the problem is tiny. At 50 assets it would need to evaluate 10^15 assignments, and grows infeasibly large. On the other hand, the greedy method reads only the return, while the pairwise and cubic terms are where you modelled the concentration risk and priced that in.

5Step 5: Plot the frontier

A single number ("cost −56.5") does not show why the pick is good. Plot every feasible portfolio in risk/return space and the answer is visual. There are only 70 ways to choose 4 of 8, so you can draw all of them.

This step needs matplotlib from Step 0.

import itertools
import matplotlib.pyplot as plt

GOLD, AMARANTH, BG, FG = "#b38f12", "#a52469", "#EDECEC", "#0f1319"

# Complete these two functions (Task 5) before running the cell:
# with the `...` stubs in place the plot below raises a TypeError.
def port_return(bits):
    ...        # Task 5

def port_risk(bits):
    ...        # Task 5

pts = [(port_risk(b), port_return(b)) for b in
       ([1 if i in c else 0 for i in range(n)]
        for c in itertools.combinations(range(n), k))]

quantum = [1, 0, 1, 1, 0, 0, 1, 0]                          # your verified solution
gcombo  = sorted(range(n), key=lambda i: -returns[i])[:k]
greedy  = [1 if i in gcombo else 0 for i in range(n)]

fig, ax = plt.subplots(figsize=(7, 5))
fig.patch.set_facecolor(BG); ax.set_facecolor(BG)
ax.scatter([p[0] for p in pts], [p[1] for p in pts], c="#cfcdcd", s=30, label="all portfolios")
ax.scatter(port_risk(quantum), port_return(quantum), c=GOLD, s=190, edgecolors="white", lw=1.5, zorder=3, label="Kipu quantum optimum")
ax.scatter(port_risk(greedy), port_return(greedy), facecolors="none", edgecolors=GOLD, s=190, lw=2.4, zorder=2, label="naive (top-4 return)")
ax.set_xlabel("Correlation risk score (lower is safer)", color=FG)
ax.set_ylabel("Expected return (%)", color=FG)
for sp in ("top", "right"): ax.spines[sp].set_visible(False)
ax.legend(frameon=False)
plt.show()

def port_return(bits): return sum(returns[i] for i in range(n) if bits[i])

def port_risk(bits):
    r = sum(rho for (i, j), rho in correlation.items() if bits[i] and bits[j])
    if bits[0] and bits[1] and bits[6]:
        r += 20.0 / risk_penalty            # three-body term, normalized risk scale

import itertools
import matplotlib.pyplot as plt

# Kipu brand colors (tokens.css)
GOLD, AMARANTH, BG, FG = "#b38f12", "#a52469", "#EDECEC", "#0f1319"

def port_return(bits):
    return sum(returns[i] for i in range(n) if bits[i])

def port_risk(bits):
    r = sum(rho for (i, j), rho in correlation.items() if bits[i] and bits[j])
    if bits[0] and bits[1] and bits[6]:
        r += 20.0 / risk_penalty
    return r

pts = []
for combo in itertools.combinations(range(n), k):
    bits = [1 if i in combo else 0 for i in range(n)]
    pts.append((port_risk(bits), port_return(bits)))

quantum = [1, 0, 1, 1, 0, 0, 1, 0]                          # verified solution
gcombo  = sorted(range(n), key=lambda i: -returns[i])[:k]   # naive pick
greedy  = [1 if i in gcombo else 0 for i in range(n)]

fig, ax = plt.subplots(figsize=(7, 5))
fig.patch.set_facecolor(BG); ax.set_facecolor(BG)
ax.scatter([p[0] for p in pts], [p[1] for p in pts], c="#cfcdcd", s=30, label="all portfolios")
ax.scatter(port_risk(quantum), port_return(quantum), c=GOLD, s=190, edgecolors="white", lw=1.5, zorder=3, label="Kipu quantum optimum")
ax.scatter(port_risk(greedy),  port_return(greedy),  facecolors="none", edgecolors=GOLD, s=190, lw=2.4, zorder=2, label="naive (top-4 return)")
ax.set_xlabel("Correlation risk score (lower is safer)", color=FG)
ax.set_ylabel("Expected return (%)", color=FG)
ax.set_title("4-of-8 portfolios: risk vs return", color=FG, fontweight="bold", loc="left")
for s in ("top", "right"): ax.spines[s].set_visible(False)
ax.legend(frameon=False)
plt.show()

An honest reading of that chart is more nuanced. The naive pick is not irrational: it sits on the efficiency frontier, at the maximum-return corner. It takes on large risk for the last bit of return. Going from the quantum optimum to the naive pick buys +5 points of return for +3.4 points of risk, while the step into the optimum bought +10 return for +0.65 risk. The risk appetite you encoded (risk_penalty = 10) controls where the optimizer lands and can be tuned. It is a modeling choice.

6Step 6: Will it run on hardware?

Everything so far ran on a simulator, and a simulator has no chip: no error rates, no wiring, no clock. Point the same dictionary at hardware and the question becomes whether your problem and the machine line up.

Three properties of your problem meet three properties of the hardware. Call it the 3 + 3.

The rows are not independent, which is the whole subtlety. Density and order both inflate the gate count, and the layout decides how much extra you pay to deliver those gates, so the bottom two rows meet in one number. Only the top row is a clean ceiling: 8 variables need 8 qubits, and every chip in the Hub list has hundreds.

That is why qubit count is the least interesting number on a spec sheet at these sizes: a sparse model with many variables can run where a dense model with far fewer variables fails, because interactions, not variables, drive the gate count.

Three quantities turn the bottom two rows into arithmetic.

1. What your model demands. Every interaction in your dictionary becomes two-qubit gates. A d-body term costs 2(d − 1) of them, so each pair costs 2 and your tech triple costs 4. Sum over the dictionary and you have the gate count your model asks for, before any chip has a say.

2. What the chip charges to deliver it. Two variables can only interact if their qubits are physically wired together. When they are not, the compiler shuttles them toward each other, and every move is more two-qubit gates. The wiring pattern therefore multiplies your gate count:

Two chips with identical gate quality can differ by a large factor purely from how they are wired. Whatever that factor is for your problem, you learn it by transpiling, not by looking it up: routing cost depends on your model's density as much as on the lattice, and a sparse model on heavy-hex can beat a dense one on a grid.

3. What the chip gets wrong. Every two-qubit gate has a failure probability, published per chip and refreshed daily in its calibration data.

Put the three together and you have one number:

λ  =  N_2q  ×  r  ×  p_2q

N_2q   two-qubit gates your model needs, summed as 2(d − 1) per d-body term
r      routing factor of the wiring, measured by transpiling, not looked up
p_2q   error per two-qubit gate, from the chip's live calibration

Task 6 computes this with the routing factor left at one. That gives you a floor, not a forecast: a real chip can only be worse, and how much worse is a transpiler question rather than a lookup.

λ is the expected number of errors in a single run, so the fraction of runs that come back clean is

P(clean)  =  e^(−λ)              λ = 1 → 37%      λ = 3 → 5%      λ = 5 → under 1%

and only clean runs teach the optimizer anything. Small λ and the histogram concentrates on one answer, which is what convergence looks like from the outside. Large λ and you are sampling noise: you may still stumble onto the optimum, but you will have no way to tell that you did.

Both quantities live on one picture, and it is the picture worth carrying out of this lab:

Read it in three moves. Your model's gate count sets your height, so every modelling choice (fewer variables, sparser couplings, lower degree) moves you down toward safety. The chip's quality sets your horizontal position, so a better chip moves you right and can handle a taller model. Whichever side of the band you sit on tells you what to do next: run it, lean on post-processing, or reshape the model.

The same split is how to read a hardware roadmap. Improvement is k_total = k_error × k_architecture: better gates and better wiring are separate levers, and the wiring lever is measured for your problem rather than quoted from a spec sheet.

This λ is the number you will defend when you ask for real QPU budget on a problem, so learn to get it before you spend.

cal = requests.get(f"{API}/{chip}/calibration", headers={"X-Auth-Token": TOKEN}).json()
p_2q = 1 - cal["median_two_qubit_fidelity"]["value"]     # also: median_t1, median_t2

import math, os, requests

TOKEN = os.environ["KQH_PERSONAL_ACCESS_TOKEN"]
API = "https://api.hub.kipu-quantum.com/quantum/backends"

def error_per_gate(backend_id):
    cal = requests.get(f"{API}/{backend_id}/calibration",
                       headers={"X-Auth-Token": TOKEN}, timeout=30).json()
    return 1 - cal["median_two_qubit_fidelity"]["value"], cal["calibrated_at"]

def degree(key):
    if key == "()":
        return 0
    return key.count(",") if key.endswith(",)") else key.count(",") + 1

def gates_needed(hubo):
    # a d-body term costs 2(d-1) two-qubit gates; constants and linear terms cost none
    return sum(2 * (degree(key) - 1) for key in hubo if degree(key) >= 2)

n_gates = gates_needed(hubo)
for chip in ("ibm.qpu.boston", "ibm.qpu.marrakesh"):
    p, when = error_per_gate(chip)
    lam = n_gates * p
    print(f"{chip:20} p_2q={p:.1e}  lambda={lam:.3f}  clean={math.exp(-lam):.0%}  (cal {when[:10]})")

To estimate gate count: for each term in the dictionary add 2(d − 1) for every key of degree d=2 or higher. Then λ = N_2q * p_2q, ignoring routing for now.

The porfolio hubo sits at minimum 60 two-qubit gates: 28 pairs at 2 each, plus the triple at 4. Against calibration from 24 August 2026, and before any routing overhead and asusming a perfect gate set:

Both sit comfortably below lambda < 2, which is expected for an 8-variable problem: it fits anywhere, and the QPU chip barely matters. When the model grows, and more terms are added, lambda increses rapidly. A fully connected quadratic model has n(n−1)/2 pairs:

Density, not variable count, is the important factor driving the gate count. To get around that limitation Part 2 introduces a decomposition technique for large problems. It cuts them up into hardware-sized pieces instead.

Caution: the calculation above assumes perfect connectivity (all-to-all) and a perfect gate sets. Actual circuits might be 4+ times larger once you account for routing and chip specifics. Transpilation_only runs can help estimate complexity and return the circuit information without submitting a job to the backend.

Documentation: Kipu Quantum Hub | Hub docs | Quickstart | Service SDK | Using a service | Access tokens | Iskay on IBM Quantum | Kipu Quantum Academy

Research: BF-DCQO, the optimizer's algorithm (arXiv:2405.13898) | BF-DCQO for higher-order (HUBO) optimization (arXiv:2409.04477) | DCQO portfolio optimization on IonQ, 20 assets (arXiv:2308.15475) | Large-scale portfolio optimization on a trapped-ion quantum computer (arXiv:2602.23976)