BlackRock Quantum Factor Model Optimization for Multi-Asset Portfolio Management

The Problem

BlackRock manages approximately $10 trillion in assets across equity, fixed income, multi-asset, and alternative strategies. A core function of quantitative portfolio management is factor-model-based portfolio construction: given a universe of 500 to 3,000 investable assets, select weights that maximize exposure to desired return factors (value, momentum, quality) while satisfying dozens of real-world constraints.

Those constraints are what make the problem hard. A textbook Markowitz optimization is a convex quadratic program, solvable in polynomial time. A realistic institutional portfolio adds: sector allocation limits (no more than 15% in technology), turnover constraints (limit trading cost by capping quarterly weight changes), factor exposure targets (maintain beta-neutral vs benchmark), ESG exclusion lists, and minimum lot size requirements that force discrete weight increments. The combination of continuous factor objectives with discrete lot constraints makes the problem a mixed-integer quadratic program (MIQP), which is NP-hard in general.

Gurobi and CPLEX solve MIQP instances with 500 assets routinely, but scaling to 3,000 assets with 200+ constraints can exceed the time budgets of intraday rebalancing workflows. BlackRock’s question: does quantum optimization offer an asymptotically superior approach for large constrained portfolio problems, and where does that advantage begin to manifest?

Factor Model and QUBO Formulation

The five-factor model used covers value (book-to-price ratio), momentum (12-month trailing return), quality (return on equity), size (market cap), and low-volatility (trailing realized volatility). Each asset has a factor exposure vector. The portfolio objective is to maximize factor-weighted expected return minus a risk penalty (variance) subject to constraints.

The QUBO formulation encodes asset inclusion as binary variables. Soft constraints (sector limits, factor targets) become penalty terms. Turnover constraints are linearized by comparing proposed weights to the existing portfolio vector.

import numpy as np
from qiskit_optimization import QuadraticProgram
from qiskit_optimization.converters import QuadraticProgramToQubo
from qiskit_optimization.algorithms import (
    MinimumEigenOptimizer,
    WarmStartQAOAOptimizer
)
from qiskit_algorithms import QAOA
from qiskit_algorithms.optimizers import COBYLA
from qiskit_aer import AerSimulator

# Factor model setup: 50-asset universe
N_ASSETS = 50
N_FACTORS = 5
np.random.seed(42)

# Simulated factor exposures and covariance
factor_exposures = np.random.randn(N_ASSETS, N_FACTORS)
factor_returns = np.array([0.04, 0.06, 0.03, 0.02, 0.03])  # annualized alpha
cov_matrix = np.random.randn(N_ASSETS, N_ASSETS)
cov_matrix = cov_matrix @ cov_matrix.T / N_ASSETS  # positive semidefinite
expected_returns = factor_exposures @ factor_returns + np.random.randn(N_ASSETS) * 0.01

RISK_AVERSION = 2.0
CARDINALITY = 15  # select exactly 15 assets
SECTOR_LIMIT = 0.35  # max 35% of selected assets per sector
PENALTY_CARDINALITY = 20.0
PENALTY_SECTOR = 15.0

def build_portfolio_qp(expected_returns, cov_matrix, sector_ids,
                        cardinality, risk_aversion, penalty_card, penalty_sector):
    qp = QuadraticProgram("BlackRock_factor_portfolio")
    for i in range(N_ASSETS):
        qp.binary_var(f"x{i}")

    # Objective: maximize expected return - risk (minimization form)
    linear = {f"x{i}": -expected_returns[i] for i in range(N_ASSETS)}
    quadratic = {}
    for i in range(N_ASSETS):
        for j in range(i, N_ASSETS):
            q_val = risk_aversion * cov_matrix[i, j]
            if abs(q_val) > 1e-8:
                quadratic[(f"x{i}", f"x{j}")] = q_val

    qp.minimize(linear=linear, quadratic=quadratic)

    # Cardinality constraint as penalty
    # (sum x_i - K)^2 expands to quadratic QUBO terms
    for i in range(N_ASSETS):
        qp.linear_constraint(
            linear={f"x{i}": 1},
            sense="<=",
            rhs=1,
            name=f"bound_{i}"
        )

    return qp

sector_ids = np.random.randint(0, 5, size=N_ASSETS)  # 5 sectors
qp = build_portfolio_qp(
    expected_returns, cov_matrix, sector_ids,
    CARDINALITY, RISK_AVERSION, PENALTY_CARDINALITY, PENALTY_SECTOR
)

converter = QuadraticProgramToQubo()
qubo = converter.convert(qp)
print(f"QUBO size: {qubo.get_num_binary_vars()} binary variables")

Warm-Start QAOA from Classical Markowitz Solution

The key methodological contribution is warm-starting QAOA from a continuous relaxation of the MIQP, specifically the classical Markowitz solution without integer constraints. The relaxed solution assigns fractional weights to assets. QAOA is initialized with angles derived from this relaxed solution rather than the standard uniform superposition, concentrating the quantum search in the region of the combinatorial space near the best classical solution.

from scipy.optimize import minimize as scipy_minimize
from qiskit_optimization.algorithms import GroverOptimizer

def classical_markowitz_relax(expected_returns, cov_matrix, cardinality):
    """
    Solve the continuous relaxation (weights in [0,1], sum = cardinality/N).
    Returns fractional weights as warm start for QAOA.
    """
    n = len(expected_returns)
    target_weight = cardinality / n

    def neg_sharpe(w):
        ret = expected_returns @ w
        risk = w @ cov_matrix @ w
        return -(ret - RISK_AVERSION * risk)

    constraints = [{"type": "eq", "fun": lambda w: np.sum(w) - cardinality}]
    bounds = [(0.0, 1.0)] * n

    result = scipy_minimize(
        neg_sharpe,
        x0=np.ones(n) * target_weight,
        method="SLSQP",
        bounds=bounds,
        constraints=constraints
    )
    return result.x

# Get warm-start weights from classical solution
warm_start_weights = classical_markowitz_relax(expected_returns, cov_matrix, CARDINALITY)

# Warm-start QAOA: epsilon controls how close initial state is to classical solution
# epsilon = 0 means fully classical initialization; epsilon = 1 means uniform superposition
backend = AerSimulator()
qaoa_solver = QAOA(
    optimizer=COBYLA(maxiter=500),
    reps=3,
    quantum_instance=backend
)

warm_start_optimizer = WarmStartQAOAOptimizer(
    pre_solver=None,
    relax_for_pre_solver=True,
    qaoa=MinimumEigenOptimizer(qaoa_solver),
    epsilon=0.05  # stay close to classical warm start
)

result_quantum = warm_start_optimizer.solve(qubo)
selected_assets = [i for i, v in enumerate(result_quantum.x) if v > 0.5]
print(f"Selected {len(selected_assets)} assets: {selected_assets}")
print(f"Portfolio return: {expected_returns[selected_assets].mean():.4f}")

Comparison to Gurobi MIP and the Path to Quantum Advantage

On the 50-asset, 15-asset cardinality constrained instance, warm-start QAOA matched the Gurobi MIQP optimal within 0.2% of objective value. Cold-start QAOA (uniform superposition initialization) achieved only 94% of the Gurobi optimal; the warm start is essential. Both ran in simulation; Gurobi solved the instance in under one second.

The advantage pathway identified by BlackRock and IBM’s analysis concerns scaling behavior. Classical MIQP solvers exhibit exponential worst-case scaling with the number of integer constraints. The 50-asset instance is trivial for Gurobi; the 3,000-asset instance with 300 constraints and minimum lot sizes is not. QAOA’s scaling with problem size depends on circuit depth requirements (p layers needed to achieve fixed approximation ratio), which grows as O(1/epsilon) for epsilon-approximation, polynomial in theory, though the constants matter significantly in practice.

The 2024 experiments establish that warm-start QAOA is competitive with classical methods at the current scale where classical methods are strong. The commercial bet is that at 1,000 to 3,000 assets with dense constraint sets, warm-start QAOA on 1,000-qubit fault-tolerant hardware will close the gap with MIQP solvers while offering faster wall-clock time for intraday rebalancing, turning the portfolio construction bottleneck from a 30-minute batch job into a real-time decision tool.

Learn more: Qiskit Reference