Educational disclaimer. This is a hands-on starter tutorial in the beautiful field of quantum computing applications for finance, meant for learning. The data is synthetic and the models deliberately simplified; nothing here is complete or constitutes financial advice.
The qubit wall, and how to get past it
The free simulator tops out around 20 qubits, about 20 assets in one direct call. To handle 100 assets, or the paper's 250, you decompose. Meaning you split the universe into hardware-sized clusters, solve each on its own, and stitch the pieces back together. Following the pipeline of "Large-scale portfolio optimization on a trapped-ion quantum computer" (the Kipu/IonQ paper, summarized in The research behind it below) is part of the skillset that can make quantum optimization useful on near-term hardware at all.
Now, let's run a genuine version of that pipeline on a 100-asset universe (six sectors, hold exactly 50), a problem with 100-choose-50 ≈ 10²⁹ feasible portfolios, far beyond enumeration and far beyond one simulator call. The qubit budget is a hard limit: the free simulator sits at 20, and backends on Kipu's Quantum Hub and Qiskit Functions run larger (IBM Nighthawk: 120 qubits, Heron: 156). This influences the amount and size of clusters we need.
A note on the model. Unlike the paper, which is strictly quadratic (covariance is pairwise), we include the three-body terms from Part 1's encoding. A HUBO can express complex multi-variable interactions, and BF-DCQO solves them natively.
1. Decompose
Split a 100-asset universe into hardware-sized clusters, solve each on Miray, then recombine and repair to a feasible book.
2. Read & verify
Read the Miray results, brute-force check each cluster, and benchmark the portfolio against strong classical baselines.
3. Tune
Sweep the risk-aversion dial λ and watch the portfolio trace the efficient frontier.
Get up to speed
Part 2 is self-contained, but it assumes Part 1's two essentials: a Miray client and a HUBO formulation. Drop these cells in first. (New here? Part 1 builds both from scratch on an 8-asset book.)
1. The Miray client (same as Part 1, Step 2 Option B, Kipu Hub only). Hub docs: Quickstart, Service SDK, access tokens, using a service.
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,
)
2. The formulation, the advanced way. Part 1 hand-builds the HUBO dictionary to teach the encoding. Since you are here to scale, take the declarative path instead: state the model with qiskit-addon-opt-mapper (variables, objective, constraints) and let it expand the cardinality penalty for you. One small adapter folds the mapper's output into the plain dict Iskay and Miray expect, because binary x² = x, so each diagonal (i, i) term collapses into the linear one. This is the detailed version of the mapper that Part 1 only points to:
from qiskit_addon_opt_mapper import OptimizationProblem # pip install qiskit-addon-opt-mapper
from qiskit_addon_opt_mapper.converters import OptimizationProblemToHubo
# A tiny 6-asset stand-in for Part 1's book, formulated declaratively.
tickers = ["AAPL", "MSFT", "JPM", "XOM", "JNJ", "NVDA"]
returns = [18.0, 16.0, 12.0, 11.0, 8.0, 22.0]
correlation = {(0, 1): 0.70, (0, 5): 0.65, (1, 5): 0.68, (3, 4): 0.40}
triple = (0, 1, 5, 20.0) # AAPL-MSFT-NVDA joint-crash penalty
k, n = 3, 6 # hold exactly 3 of 6
risk_penalty, card_penalty = 8.0, 50.0 # Part 2's convention (Part 1's toy used 10.0)
# 1. MODEL declaratively: one binary variable per asset, in asset order.
op = OptimizationProblem("portfolio")
for t in tickers:
op.binary_var(name=t) # creation order == asset index
# 2. OBJECTIVE: reward return (linear), penalize correlated pairs (quadratic),
# add the tech-triple joint crash (higher_order). No cardinality term here.
op.minimize(
linear={tickers[i]: -returns[i] for i in range(n)},
quadratic={(tickers[i], tickers[j]): risk_penalty * rho
for (i, j), rho in correlation.items()},
higher_order={3: {(tickers[triple[0]], tickers[triple[1]], tickers[triple[2]]): triple[3]}},
)
# 3. CONSTRAINT: "hold exactly k" as a real equality. The mapper expands it into
# penalty terms, so you never hand-derive (Sum x - k)^2.
op.linear_constraint(linear={t: 1 for t in tickers}, sense="EQ", rhs=k, name="cardinality")
# 4. CONVERT to a HUBO with an explicit penalty (penalty=None auto-derives a
# safe-but-large value that drowns the risk signal).
obj = OptimizationProblemToHubo(penalty=card_penalty).convert(op).objective
# 5. ADAPT to the plain Iskay/Miray dict: binary x^2 = x, so fold each diagonal
# (i, i) quadratic term into the linear one.
def to_kipu_hubo(obj, n):
lin = {i: float(v) for i, v in obj.linear.to_dict().items()}
h = {"()": float(obj.constant)}
for (i, j), v in obj.quadratic.to_dict().items():
if i == j:
lin[i] = lin.get(i, 0.0) + float(v) # x_i^2 -> x_i
else:
h[f"({min(i, j)}, {max(i, j)})"] = float(v)
for i in range(n):
h[f"({i},)"] = lin.get(i, 0.0)
for order, expr in obj.higher_order.items(): # three-body and beyond
for idx, c in expr.to_dict().items():
key = tuple(sorted(idx))
h["(" + ", ".join(map(str, key)) + ")"] = float(c)
return h
hubo = to_kipu_hubo(obj, n)
print(len(hubo), "terms;", hubo["(0, 1, 5)"], "= the three-body penalty BF-DCQO takes natively")
# -> 23 terms; 20.0 = the three-body penalty BF-DCQO takes natively
3. Confirm the setup end to end on the free simulator (a 6-asset warm-up; the real 100-asset pipeline is Step 1):
execution = client.run(request={
"problem": hubo, "problem_type": "binary",
"shots": 2000, "num_iterations": 6, "num_greedy_passes": 1,
})
execution.wait_for_final_state()
print(execution.result().result["cost"], execution.result().result["mapped_solution"])
That to_kipu_hubo adapter is the bridge you reuse below: Step 1's per-cluster builder emits the very same dictionary shape, just specialized to one hardware-sized cluster at a time.
Step 1: Decompose and run the pipeline
Run this below the client you just set up. The universe is generated deterministically so the snippet is self-contained; in practice you load real returns and a denoised correlation matrix instead. For example from yfinance.
import itertools, random, math
# --- Illustrative synthetic universe: 100 names across 6 sectors. ---
def build_universe(seed=11):
random.seed(seed)
spec = [("TEC", 28, 20, 5, 0.72), ("FIN", 22, 13, 4, 0.55),
("ENE", 16, 11, 4, 0.68), ("HEA", 14, 9, 3, 0.50),
("IND", 12, 12, 4, 0.60), ("UTI", 8, 7, 2, 0.45)]
names, sector, returns = [], [], []
for tag, n_s, mu, sd, _ in spec:
for i in range(n_s):
names.append(f"{tag}{i+1}"); sector.append(tag)
returns.append(round(mu + random.uniform(-sd, sd), 1))
intra = {tag: rho for tag, _, _, _, rho in spec}
correlation = {}
for i in range(len(names)):
for j in range(i + 1, len(names)):
if sector[i] == sector[j]:
rho = round(min(0.92, max(0.2, intra[sector[i]] + random.uniform(-0.12, 0.12))), 3)
else:
rho = round(max(0.0, random.uniform(0.0, 0.18)), 3)
if rho:
correlation[(i, j)] = rho
return names, sector, returns, correlation
names, sector, returns, correlation = build_universe()
N, K, QMAX = len(names), 50, 20 # 100 assets, hold 50, 20-qubit budget
risk_penalty = 8.0 # risk appetite (lambda); Part 1's toy used 10.0, swept in Step 3
# Three-body terms: joint co-movement a QUBO cannot express. Decomposition keeps
# a triple only if all three members land in one cluster.
triples = [(2, 13, 9, 18.0), (2, 14, 18, 14.4), (36, 37, 48, 12.6)]
def abscorr(i, j):
return 0.0 if i == j else correlation.get((min(i, j), max(i, j)), 0.0)
def full_cost(bits):
e = -sum(returns[i] for i in range(N) if bits[i])
for (i, j), rho in correlation.items():
if bits[i] and bits[j]:
e += risk_penalty * rho
for (i, j, k, p) in triples:
if bits[i] and bits[j] and bits[k]:
e += p
return e
# 1. COMMUNITY DETECTION: threshold the correlation graph, take connected
# components (a stand-in for the paper's RMT denoise + Louvain).
def connected_components(threshold=0.30):
adj = {i: set() for i in range(N)}
for (i, j), rho in correlation.items():
if rho >= threshold:
adj[i].add(j); adj[j].add(i)
seen, comps = set(), []
for s in range(N):
if s in seen:
continue
stack, comp = [s], []
while stack:
u = stack.pop()
if u in seen:
continue
seen.add(u); comp.append(u); stack.extend(adj[u] - seen)
comps.append(sorted(comp))
return comps
# 2. GREEDY SPLIT (paper Algorithm 1): break any community larger than the qubit
# budget by seeding on its highest-degree node and grabbing its strongest
# neighbours, until every cluster fits QMAX.
def greedy_split(community, qmax=QMAX):
R, out = list(community), []
while R:
if len(R) <= qmax:
out.append(sorted(R)); break
deg = {i: sum(abscorr(i, j) for j in R) for i in R}
seed = max(R, key=lambda i: deg[i])
chosen = {seed}
for j in sorted(R, key=lambda j: -abscorr(seed, j)):
if len(chosen) >= qmax:
break
chosen.add(j)
out.append(sorted(chosen)); R = [i for i in R if i not in chosen]
return out
# 3. PER-CLUSTER HUBO: returns + intra-cluster risk + intra-cluster 3-body terms
# + a cardinality penalty (paper's adaptive Lambda = 2 * max influence).
def build_cluster_hubo(cluster, budget):
m = len(cluster)
glob2loc = {g: l for l, g in enumerate(cluster)}
lin = {l: -returns[cluster[l]] for l in range(m)}
quad = {(a, c): risk_penalty * abscorr(cluster[a], cluster[c])
for a in range(m) for c in range(a + 1, m) if abscorr(cluster[a], cluster[c])}
influence = [abs(lin[l]) + sum(abs(v) for (a, c), v in quad.items() if l in (a, c)) for l in range(m)]
Lam = 2 * max(influence)
hubo = {"()": Lam * budget ** 2}
for l in range(m):
hubo[f"({l},)"] = lin[l] + Lam * (1 - 2 * budget)
for a in range(m):
for c in range(a + 1, m):
hubo[f"({a}, {c})"] = quad.get((a, c), 0.0) + 2 * Lam
for (i, j, k, p) in triples:
if i in glob2loc and j in glob2loc and k in glob2loc:
key = tuple(sorted((glob2loc[i], glob2loc[j], glob2loc[k])))
hubo[f"({key[0]}, {key[1]}, {key[2]})"] = p
return hubo, glob2loc
# 4. SOLVE the clusters on Miray. The clusters are independent subproblems, so
# SUBMIT them all first (non-blocking), then COLLECT once they finish. Total
# wait is the slowest single cluster, not the sum of all eight.
def submit_to_miray(cluster, budget, client):
hubo, glob2loc = build_cluster_hubo(cluster, budget)
execution = client.run(request={
"problem": hubo, "problem_type": "binary",
"shots": 4000, "num_iterations": 10, "num_greedy_passes": 2,
})
return execution, glob2loc # queued, not awaited
def collect(execution, glob2loc):
execution.wait_for_final_state()
sol = execution.result().result["mapped_solution"] # local index -> 0/1
loc2glob = {l: g for g, l in glob2loc.items()}
return [loc2glob[int(l)] for l in sol if sol[l] == 1]
# 5. RECOMBINE the cluster picks, then REPAIR to exactly K with a
# cardinality-preserving swap local search on the FULL objective. This is
# what recovers the cross-cluster couplings the decomposition dropped.
def recombine(cluster_solutions):
bits = [0] * N
for picks in cluster_solutions:
for i in picks:
bits[i] = 1
return bits
def repair_swap(bits, K):
bits = bits[:]
flip = lambda b, i: b[:i] + [1 - b[i]] + b[i + 1:]
while sum(bits) < K:
outs = [i for i in range(N) if not bits[i]]
bits[min(outs, key=lambda i: full_cost(flip(bits, i)))] = 1
while sum(bits) > K:
ins = [i for i in range(N) if bits[i]]
bits[max(ins, key=lambda i: full_cost(flip(bits, i)))] = 0
improved = True
while improved:
improved = False
base = full_cost(bits)
ins = [i for i in range(N) if bits[i]]; outs = [i for i in range(N) if not bits[i]]
best_gain, best_swap = 0.0, None
for a in ins:
for b in outs:
bits[a], bits[b] = 0, 1
gain = base - full_cost(bits)
bits[a], bits[b] = 1, 0
if gain > best_gain + 1e-9:
best_gain, best_swap = gain, (a, b)
if best_swap:
a, b = best_swap; bits[a], bits[b] = 0, 1; improved = True
return bits
# --- run it ---
communities = connected_components(0.30)
clusters = []
for c in communities:
clusters += greedy_split(c) if len(c) > QMAX else [c]
clusters.sort(key=lambda c: -len(c))
budgets = [round(len(c) / 2) for c in clusters] # K_sub = size / 2, sums to 50
jobs = [submit_to_miray(c, b, client) for c, b in zip(clusters, budgets)] # fan out
solutions = [collect(ex, g2l) for ex, g2l in jobs] # gather
portfolio = repair_swap(recombine(solutions), K)
print(round(full_cost(portfolio), 2), "holding", sum(portfolio), "names")
The 100-asset universe decomposes into 8 clusters, none larger than 20 qubits.
cluster sector size budget carries 3-body terms
0 Tech 20 10 yes (two tech triples)
1 Financials 20 10 yes (one financials triple)
2 Energy 16 8
3 Health 14 7
4 Industrials 12 6
5 Tech (split) 8 4
6 Utilities 8 4
7 Financials (sp) 2 1
The two oversized communities get split by the greedy rule: Tech (28 names) into 20 + 8, Financials (22) into 20 + 2. Let's visualize the results.
import matplotlib.pyplot as plt
from collections import Counter
# Left: how the 100 assets split into clusters, each under the qubit budget,
# with the share of names each cluster contributed to the final portfolio.
labels = [f"{Counter(sector[i] for i in c).most_common(1)[0][0]} ({len(c)})" for c in clusters]
held_per_cluster = [sum(portfolio[i] for i in c) for c in clusters]
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
y = range(len(clusters))
ax1.barh(y, [len(c) for c in clusters], color="lightsteelblue", label="cluster size (qubits)")
ax1.barh(y, held_per_cluster, color="steelblue", label="names held")
ax1.axvline(QMAX, color="crimson", ls="--", lw=1.5, label=f"qubit budget = {QMAX}")
ax1.set_yticks(list(y)); ax1.set_yticklabels(labels); ax1.invert_yaxis()
ax1.set_xlabel("assets / qubits")
ax1.set_title(f"100 assets decomposed into {len(clusters)} clusters")
ax1.legend(fontsize=8)
# Right: the 50-name portfolio stays diversified (Tech has 28 names available
# but the optimizer refuses to over-concentrate in it).
sectors = sorted(set(sector))
avail = [sum(1 for i in range(N) if sector[i] == s) for s in sectors]
held = [sum(portfolio[i] for i in range(N) if sector[i] == s) for s in sectors]
x = range(len(sectors))
ax2.bar([i - 0.2 for i in x], avail, 0.4, color="lightgray", label="available")
ax2.bar([i + 0.2 for i in x], held, 0.4, color="seagreen", label="held")
ax2.set_xticks(list(x)); ax2.set_xticklabels(sectors); ax2.set_ylabel("names")
ax2.set_title(f"{K}-name portfolio: balanced across sectors")
ax2.legend(fontsize=8)
plt.tight_layout(); plt.show()
Step 2: Read & verify the results from Miray
The script above submits all independent clusters to Miray in parallel. In practice you often reserve quantum solvers for the clusters that need it or combine them with classical ones. For example, running the two largest (the 20-qubit Tech and Financials subproblems with higher-order terms) on the Miray simulator at 4000 shots, 10 iterations, 2 greedy passes, and solving the remaining six exactly. This describes a quantum-classical hybrid workflow, which is precisely what brings quantum closer to industrial usefulness.
Cluster 0 (Tech, 20 qubits, budget 10) Miray cost 50.548 = brute-force optimum (50.55)
held -> TEC2 TEC5 TEC10 TEC12 TEC14 TEC15 TEC19 TEC20 TEC21 TEC28
Cluster 1 (Financials, 20 qubits, budget 10) Miray cost 39.472 = brute-force optimum (39.47)
held -> FIN3 FIN4 FIN6 FIN9 FIN10 FIN11 FIN17 FIN18 FIN19 FIN22
BF-DCQO solved both 20-qubit subproblems, higher-order terms included, with no QUBO reduction and no auxiliary qubits, the higher-order capability demonstrated on IBM hardware in arXiv:2409.04477. Read them with response.result["mapped_solution"] (local index, not the transpiled bitstring). Then recombine all eight cluster solutions and repair to a feasible 50-name portfolio. Full comparison (lower is better):
Greedy (top-50 by return) +2313.04 over-concentrated
Quantum decomposition, recombined only +1066.87 eight cluster optima without repair
Random (best of 30,000 draws) +994.36
Equal-per-sector (top return) +978.61 naive diversification: ~8 per sector
Simulated annealing (60 restarts) +882.94 strong classical baseline
Quantum decomposition + repair +882.50 balanced across all six sectors
Reading the table top to bottom:
- Greedy (top-50 by return) ignores risk entirely: Part 1's concentration trap, at scale.
- Quantum decomposition, recombined only glues the eight cluster optima together with no global fix-up. Each piece is solved well, but every cross-cluster coupling was dropped, hence the gap.
- Random and equal-per-sector are naive diversification baselines.
- Simulated annealing (60 restarts) is the strong classical reference.
- Quantum decomposition + repair is the full pipeline, and the number to judge: solve each cluster on Miray, recombine the picks, then run a cardinality-preserving swap local search over the full objective to recover the cross-cluster couplings the split dropped.
To see this rather than rank it: we cannot enumerate all 10²⁹ portfolios, so we sample a cloud of random feasible 50-name picks and drop the two real portfolios onto it. The objective splits cleanly into the two plotted axes: cost = risk − return.
import matplotlib.pyplot as plt
import random as rnd
def port_return(bits):
return sum(returns[i] for i in range(N) if bits[i])
def port_risk(bits): # pairwise + 3-body penalties
r = sum(risk_penalty * rho for (i, j), rho in correlation.items() if bits[i] and bits[j])
return r + sum(p for (i, j, k, p) in triples if bits[i] and bits[j] and bits[k])
def topk_by_return(K): # the naive baseline
bits = [0] * N
for i in sorted(range(N), key=lambda i: -returns[i])[:K]:
bits[i] = 1
return bits
greedy = topk_by_return(K)
sampler = rnd.Random(7)
cloud = []
for _ in range(3000):
bits = [0] * N
for i in sampler.sample(range(N), K):
bits[i] = 1
cloud.append(bits)
plt.figure(figsize=(7, 5))
plt.scatter([port_risk(b) for b in cloud], [port_return(b) for b in cloud],
c="lightgray", s=10, label="random feasible (sampled)")
plt.scatter(port_risk(greedy), port_return(greedy), c="red", s=160, marker="X",
label="naive (top-50 return)", zorder=3)
plt.scatter(port_risk(portfolio), port_return(portfolio), c="green", s=160,
label="quantum decomposition + repair", zorder=3)
plt.xlabel("Correlation + 3-body risk score (lower is safer)")
plt.ylabel("Expected return")
plt.title("50-of-100 portfolios: risk vs return")
plt.legend(fontsize=8)
plt.tight_layout(); plt.show()
The naive portfolio sits alone in the high-risk corner: it wins on raw return but pays for it in risk, the concentration trap from Part 1. The quantum portfolio sits at the low-risk edge of the cloud, undominated: no sampled portfolio beats it on both axes at once.
Three honest takeaways:
- The decomposition ties a strong classical baseline and demolishes the naive ones. The repaired portfolio (+882.50) is level with 60-restart simulated annealing (+882.94) and less than half the cost of greedy (+2313), which takes all high-return assets and gets penalties from both correlations and three-body penalties. The model also beats careful diversification. The hand-built rule that holds roughly equal names per sector (+978.61) lands almost the same sector counts but still costs ~10% more, because the optimizer's real work is choosing which names within each sector to avoid the most mutually (pairs) and many-body-correlated ones (triplets).
- The repair step is not optional, and the pre-repair number is not a quantum failure. Recombining the eight cluster optima scores +1066.87, a 21% gap above best-known. That gap is not a solver-quality problem (each cluster was solved well); it is the cost of dropping cross-cluster couplings and of the fixed per-cluster budget split. The simple global swap search (post-processing) closes it from +1066.87 to +882.50. Decomposition needs both halves: a good per-cluster solver and a global repair to stitch the clusters back together.
- The three-body terms changed the final portfolio. In the Tech cluster the higher-order penalty pushed the optimizer to drop TEC3 (which sits in both tech triples) in favor of TEC12, avoiding exactly the joint co-movement a pairwise model does not include. BF-DCQO solved both higher-order subproblems natively, no QUBO reduction and no auxiliary qubits.
Verification. Each cluster is checked against its exact (brute-force) optimum, which is feasible because every cluster is ≤ 20 assets. The full 100-asset portfolio is not, at least on my Laptop in a reasonable amount of time (10²⁹ combinations). So we compare it against a high-quality feasible result from a strong classical heuristic without a proven global optimum. In practice, commercial and/or open-source solvers like CPLEX and Gurobi provide the reference the quantum-centric solution is measured against. Careful: not all commercial solver allow to be benchmarked, check your agreements.
Step 3: Tune risk aversion with one variable
The portfolio above looks cautious because it uses risk_penalty = 8.0, the risk side is about 70% of the objective, so the optimizer spreads across sectors. risk_penalty is the risk-aversion λ: a single coefficient that says how many units of return you will give up to remove one unit of correlation risk. If you write the objective with λ pulled out as a variable, you can sweep it. λ scales the correlation penalty, wheras the higher-order three-body term stays a fixed structural penalty (a joint crash is a "no no" at any risk appetite, you might want to control that assumption independently):
import matplotlib.pyplot as plt
def expected_return(bits):
return sum(returns[i] for i in range(N) if bits[i])
def correlation_risk(bits): # pairwise risk the dial trades against
return sum(rho for (i, j), rho in correlation.items() if bits[i] and bits[j])
def higher_order_risk(bits): # structural, NOT scaled by the dial
return sum(s for (i, j, k, s) in triples if bits[i] and bits[j] and bits[k])
def objective(bits, risk_aversion): # risk_aversion = 8.0 reproduces full_cost
return -expected_return(bits) + risk_aversion * correlation_risk(bits) + higher_order_risk(bits)
# Trace the frontier classically (fast): from the all-return greedy seed, run a
# swap local search at each risk_aversion. The quantum pipeline above solves at
# whichever single value you pick; here we just want to see the whole dial.
def seed():
bits = [0] * N
for i in sorted(range(N), key=lambda i: -returns[i])[:K]:
bits[i] = 1
return bits
def solve_at(risk_aversion):
bits = seed()
improved = True
while improved:
improved = False
base = objective(bits, risk_aversion)
held = [i for i in range(N) if bits[i]]
free = [i for i in range(N) if not bits[i]]
best_gain, best_swap = 1e-9, None
for a in held:
for b in free:
bits[a], bits[b] = 0, 1
gain = base - objective(bits, risk_aversion)
bits[a], bits[b] = 1, 0
if gain > best_gain:
best_gain, best_swap = gain, (a, b)
if best_swap:
a, b = best_swap; bits[a], bits[b] = 0, 1; improved = True
return bits
levels = [0, 1, 2, 4, 8, 16, 32]
frontier = [solve_at(ra) for ra in levels]
# The strategies from the comparison table, placed on the same axes so you can
# see which land ON the frontier and which are dominated (above or below it).
def equal_per_sector(): # naive diversification
sectors = sorted(set(sector))
bysec = {s: [i for i in range(N) if sector[i] == s] for s in sectors}
per, extra, held = K // len(sectors), K % len(sectors), []
for n, s in enumerate(sorted(sectors, key=lambda s: -len(bysec[s]))):
take = min(per + (1 if n < extra else 0), len(bysec[s]))
held += sorted(bysec[s], key=lambda i: -returns[i])[:take]
bits = [0] * N
for i in held[:K]: bits[i] = 1
return bits
strategies = {
"quantum decomp + repair": (portfolio, "green", "*", 240), # from the run above
"equal-per-sector": (equal_per_sector(), "orange", "s", 110),
"greedy (top-50 return)": (seed(), "red", "X", 130),
}
plt.figure(figsize=(7.5, 5.5))
xs = [correlation_risk(b) for b in frontier]
ys = [expected_return(b) for b in frontier]
plt.plot(xs, ys, "-o", color="steelblue", label="classical local search, swept risk_aversion")
for ra, x, y in zip(levels, xs, ys):
plt.annotate(f"λ={ra}", (x, y), textcoords="offset points", xytext=(6, 5), fontsize=8)
for name, (bits, color, marker, size) in strategies.items():
plt.scatter(correlation_risk(bits), expected_return(bits), c=color, marker=marker,
s=size, zorder=4, edgecolors="k", linewidths=0.4, label=name)
plt.xlabel("Correlation risk (lower is safer)")
plt.ylabel("Expected return")
plt.title("Strategies vs the efficient frontier")
plt.legend(fontsize=8)
plt.tight_layout(); plt.show()
λ = 0 return 862 risk 379 return-only: rebuilds the over-concentrated tech-heavy book
λ = 1 return 803 risk 259
λ = 2 return 747 risk 220
λ = 4 return 682 risk 198
λ = 8 return 633 risk 189 <- the default used above
λ = 16 return 608 risk 188
λ = 32 return 588 risk 187 diminishing safety, real return given up
The curve is the efficient frontier. At λ=0 the search maximizes raw return and rebuilds the over-concentrated greedy portfolio (the red marker in the high-risk corner). Crank λ up and it walks down the frontier, trading return for diversification, with sharply diminishing risk reduction past the default. Many naive strategies sit below the frontier: equal-per-sector offers similar return at higher risk, and the quantum decomposition + repair result lands right on the line near λ=8, which is the point of the whole exercise. I invite you to experiment with different values of λ.
| Quantum-Centric Optimization Pipeline |
|---|
| Define assets with correlation matrix |
| Random-matrix-theory denoising + community/sector detection |
| Decompose to hardware budget (qubit count or error budget) |
| Solve clusters, orchestrate quantum-classical workflow intelligently |
| Recombine low-energy candidates, repair with swap-local-search |
| Compare with a strong classical baseline |
This pipeline is algorithmic, not hand-tuned. In the paper every row is automated: random-matrix-theory denoising strips the noise eigenvalues (those inside the Marchenko-Pastur band) from the correlation matrix, then a correlation-guided greedy split keeps strongly-coupled names together so the couplings dropped across cluster cuts are the weak ones. That single design principle, strong correlations inside clusters and only weak coupling across the cuts, is what lets the global repair step recover what decomposition throws away. Step 1 above is a small, runnable version of exactly this.
On real, densely cross-coupled S&P data, growing the qubit budget from 36 to 64 means fewer clusters, fewer dropped cross-cluster couplings, and a final portfolio that sits closer to the global optimum. Our universe is deliberately easier. Its sectors split cleanly, so the couplings across the cuts are weak and decomposition plus repair is simpler. The results start to get interesting when the correlation structure does not split neatly, and subject to further research. In the end, better hardware is better for harder, denser problems. And hardware gets better every year.
The research behind it
This lab is a runnable miniature of Kipu Quantum and IonQ's February 2026 result, Large-scale portfolio optimization on a trapped-ion quantum computer (IonQ write-up):
- 250 S&P 500 names, selecting 125, decomposed onto a 64-qubit barium trapped-ion system, with BF-DCQO solving inside each cluster.
- Larger executable subproblems reduced decomposition error and improved the final risk-return trade-off, relative to randomized baselines under identical post-processing. The IonQ write-up additionally frames the work against the Gurobi classical solver.
You build the kernel; the paper built the machine that tiles that kernel across an index. The same engine drives Kipu's real-world finance work, including the Basel III / FRTB Capital Saving case, which cut regulatory capital consumption by more than 10% on real bank portfolios.
Your move
- Run the pipeline. It runs on the free Miray simulator: solve the two large clusters on the quantum optimizer, the rest exactly, then recombine and repair.
- Grow the qubit budget. Trade the 20-qubit simulator cap for a larger backend (IBM Nighthawk 120, Heron 156): fewer clusters, fewer dropped couplings, a tighter final portfolio.
- New here? Start with Part 1. Part 1 builds and runs the optimizer on 8 assets before you scale it.
- Bring your problem. If you run a finance or energy-trading book and want to assess practical fit on current hardware, the Kipu Quantum Academy is built for you.
Documentation: Kipu Quantum Hub | Hub docs | Iskay on IBM Quantum | Kipu Quantum Academy Research: Large-scale portfolio optimization (arXiv:2602.23976) | BF-DCQO algorithm (arXiv:2405.13898) | BF-DCQO for higher-order HUBO (arXiv:2409.04477) | DCQO portfolio on IonQ, 20 assets (arXiv:2308.15475) | Feature selection HUBO on IonQ (arXiv:2604.26834)
Tested live on the Kipu Quantum Hub Miray simulator: the two 20-qubit clusters each returned their brute-force-verified optimum (costs 50.548 and 39.472), and the recombined-and-repaired 100-asset portfolio (objective +882.50) ties a 60-restart simulated-annealing baseline (+882.94), while a naive top-50 pick scores +2313.04.
Happy Optimizing 🚀