Large-scale strategy: relaxation, tail rows, screening, round-and-repair
The solver benchmark (GPU_Optimization/cuopt_benchmark, November 2026) showed
that the integer mean-CVaR MILP is an LP in disguise: HiGHS closes it at the
root node, only the tail scenario rows bind, and the cardinality constraint is
the only hard combinatorics. quantcore.optimization.large_scale exploits that:
LP relaxation by tail-row generation (evaluated at the continuous point, its
objective a valid bound), screening to the top names under a cardinality limit,
a boxed round-and-repair MILP, certification against the bound with escalation
to the integer row generation otherwise. Solver choice by size (HiGHS below
16 000 rows, cuOpt above).
Request option `strategy` ('auto' | 'milp' | 'large_scale'; auto from 200
candidates or 32 768 scenarios) in the service, the API schema and both CLIs.
Settings from the benchmark: cuOpt presolve off (its asynchronous GPU heuristic
otherwise holds every MILP for 5-25 s), dual simplex for pure LPs (concurrent
returns PDLP's inexact point on the L40S), 50 ms polling. rowgen takes an
initial working set.
Measured with HiGHS (M4 Max): N=300 S=4096 5.2 s -> 0.7 s, N=300 S=16384
45 s -> 2.9 s, N=750 S=16384 15.6 s -> 7.4 s, K=40 of 750 24.6 s -> 7.3 s,
objectives within the 1e-4 gap. Tests in tests/test_large_scale.py; benchmark
script scripts/experiments/large_scale.py; docs updated.
Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
This commit is contained in:
herzogflorian
parent
c057f8efb7
commit
fac41ffdfc
@ -17,7 +17,7 @@ never computes a number itself.
|
||||
| **PerformanceService** | `quantcore.measurement.performance` | Annualised return, volatility, Sharpe, Sortino, drawdown, Calmar, skew/kurtosis, hit ratio; rolling and EWM series; portfolios in shares. |
|
||||
| **Stock-price models** | `quantcore.models.simulation` | Geometric Brownian motion (multivariate normal) and the empirical-innovation model `eps_t = Sigma_{t-1}^{-1/2} r_t`, `r(j|T) = Sigma_T^{1/2} eps_j` with the latest rolling or EWM covariance. |
|
||||
| **Expected returns** | `quantcore.models.black_litterman` | Market-cap equilibrium returns (broad-market average, market weights), Black–Litterman on a 10-year weekly covariance with **Idzorek confidence** views (diagonal Ω), sentiment → view helper. |
|
||||
| **OptimizationService** | `quantcore.optimization` | Mean–CVaR MILP in integer share numbers: **new portfolio** (from cash) or **rebalancing** (existing holdings, transaction costs); constraints: max 10 % per stock, max per sector, no shorts, cost/turnover limits, optional cardinality; scenario sets with variance reduction (Sobol', moment matching, market-direction importance sampling) that stay identical when risk limits, constraints, costs or views change, tail-row generation for large sets, and an exact Gaussian (mean–variance-equivalent) mode; solved on cuOpt (GPU, over HTTP) with a HiGHS fallback (CPU, compiled into SciPy — runs in-process, nothing to deploy). |
|
||||
| **OptimizationService** | `quantcore.optimization` | Mean–CVaR MILP in integer share numbers: **new portfolio** (from cash) or **rebalancing** (existing holdings, transaction costs); constraints: max 10 % per stock, max per sector, no shorts, cost/turnover limits, optional cardinality; scenario sets with variance reduction (Sobol', moment matching, market-direction importance sampling) that stay identical when risk limits, constraints, costs or views change, tail-row generation for large sets, a relax-screen-round driver for large universes (`strategy`), and an exact Gaussian (mean–variance-equivalent) mode; solved on cuOpt (GPU, over HTTP) with a HiGHS fallback (CPU, compiled into SciPy — runs in-process, nothing to deploy). |
|
||||
| **HTTP API** | `quantcore.api` | FastAPI app (`/docs`) exposing all of the above; setup/update as background jobs. Runs as a container on the DGX Spark next to cuOpt, or as a Kubernetes Deployment. |
|
||||
|
||||
Documentation: **[docs/USAGE.md](docs/USAGE.md) (start here if you are building on this package: the API, the
|
||||
|
||||
@ -39,7 +39,7 @@ and reproducible through pinned data snapshots. The AI layer (news ingestion, `R
|
||||
| Market data | `quantcore.data` | Yahoo → `adj_close.parquet` (date × symbol, listing currency), `fx_chf.parquet` (date × currency); `MarketDataService.prices/returns(symbols, start, end, currency, freq)` → aligned frames |
|
||||
| Measurement | `quantcore.measurement` | returns frame → risk tables, series, matrices, `PortfolioRisk`; performance tables |
|
||||
| Models | `quantcore.models` | returns frame → covariance, `ScenarioSet` (S × N horizon log-returns, optional probabilities, fingerprint; `sampling`: Sobol'/moment matching/importance sampling), `BlackLittermanResult` |
|
||||
| Optimisation | `quantcore.optimization` | `OptimizationRequest` → `OptimizationResult` (shares, trades, cash, costs, CVaR, exposures, provenance of every model input, scenario hash, reproducibility level); `problem` (MILP, weighted scenarios or Gaussian cuts), `rowgen` (tail-row generation), `gaussian_exact` (outer approximation), `solvers` |
|
||||
| Optimisation | `quantcore.optimization` | `OptimizationRequest` → `OptimizationResult` (shares, trades, cash, costs, CVaR, exposures, provenance of every model input, scenario hash, reproducibility level); `problem` (MILP, weighted scenarios or Gaussian cuts), `rowgen` (tail-row generation), `large_scale` (relaxation, screening, round-and-repair for large universes), `gaussian_exact` (outer approximation), `solvers` |
|
||||
| Interface | `quantcore.api`, `quantcore.cli` | JSON/CLI surfaces of the above; no computation of their own |
|
||||
| Test bench | `quantcore/api/ui` | a static single-page client of the JSON API (`/ui/`): stock data, covariance and scenarios, views, constraints, optimisation with before/after. Static HTML + ES modules, no build step and no CDN, so it also runs in the offline container. It formats and charts what the services return and computes nothing itself |
|
||||
|
||||
|
||||
@ -271,6 +271,35 @@ checked, violated ones (`L_s > α*`) are added, and the loop stops when none is
|
||||
restricted problem is a relaxation, so the final solution is optimal for the full set (within the
|
||||
MIP gap). Measured: working sets ≈ 11 % of `S`, 2–3 iterations, 6–7× faster at 32 k–100 k scenarios.
|
||||
|
||||
**Large-scale strategy** (`quantcore.optimization.large_scale`, `strategy="auto"` from 200 candidates
|
||||
or 32 768 scenarios on, `"large_scale"` to force it): the solver benchmark of November 2026
|
||||
(`GPU_Optimization/cuopt_benchmark`) showed that HiGHS closes every mean–CVaR MILP at the root node —
|
||||
share counts in the hundreds round with a loss far below the MIP gap, so the branch-and-bound is paid
|
||||
for on every call and almost never used — and that the only hard combinatorics is the cardinality
|
||||
constraint. The driver therefore (1) solves the **LP relaxation** by tail-row generation, evaluated
|
||||
at the continuous solution; its objective is a valid bound for the MILP; (2) with a cardinality
|
||||
limit `K`, **screens** the universe to the `max(60, 4K)` names with the largest relaxed exposure
|
||||
(held names always stay); (3) solves the restricted model with every share count **boxed between the
|
||||
floor and the ceiling of its relaxed value** — a MILP with a tight relaxation that HiGHS solves at the
|
||||
root in a fraction of a second — and repeats the row generation on it until no scenario is violated;
|
||||
(4) **certifies** the result against the bound: within `mip_rel_gap` it is returned as `Optimal`
|
||||
with that gap, otherwise the ordinary integer row generation runs from the working set already
|
||||
found. Models below 16 000 rows go to HiGHS in-process, larger ones to cuOpt; with row generation
|
||||
the working sets stay small, so a 750-stock problem never leaves the CPU. The result's
|
||||
`solver_info["large_scale"]` carries the bound, the gap, the working set, the screened size and
|
||||
every solver call. Measured on a synthetic factor universe with Student-t scenarios
|
||||
(`scripts/experiments/large_scale.py`, HiGHS on an M4 Max, `mip_rel_gap = 1e-4`):
|
||||
|
||||
| case | integer MILP (row generation) | large-scale driver | gap to the bound |
|
||||
|---|---|---|---|
|
||||
| N = 300, S = 4 096 | 5.2 s | 0.7 s (3 solves, 485 rows) | 6.8e-5 |
|
||||
| N = 300, S = 16 384 | 45.4 s | 2.9 s (4 solves, 1 892 rows) | 9.8e-5 |
|
||||
| N = 750, S = 16 384 | 15.6 s | 7.4 s (4 solves, 1 842 rows) | 2.0e-5 |
|
||||
| N = 750, S = 16 384, K = 40 positions | 24.6 s | 7.3 s (5 solves, screened to 160 names) | 2.0e-5 |
|
||||
|
||||
The objectives agree within the gap in every case; the driver's time is the three or four
|
||||
working-set LPs of the relaxation, the boxed repair MILP takes 0.1–0.5 s.
|
||||
|
||||
**Exact Gaussian mode** (`scenario_model="gaussian_exact"`, `quantcore.optimization.gaussian_exact`):
|
||||
with normal horizon returns `R ~ N(m, C)` the CVaR is `TC − m'v + k_β √(v'Cv)` (`v = P∘n`,
|
||||
`k_β = φ(z_β)/(1−β)`); the mean–CVaR problem then coincides with the mean–variance problem (the
|
||||
|
||||
@ -302,6 +302,12 @@ that must be identical next week, freeze the data with a snapshot (`POST /snapsh
|
||||
for hours. Set it. 300 s is a sane default.
|
||||
* **The GPU solver is optional.** `solver="auto"` uses cuOpt if it answers and HiGHS otherwise; the
|
||||
result says which, and adds a warning. A cuOpt outage degrades speed, never correctness.
|
||||
* **Large universes do not need the integer MILP.** `strategy="auto"` switches to the large-scale
|
||||
driver (LP relaxation on the tail rows, screening for a cardinality limit, round-and-repair,
|
||||
certified against the relaxation bound) from 200 candidates or 32 768 scenarios on;
|
||||
`strategy="large_scale"` forces it, `"milp"` keeps the integer MILP. The result's
|
||||
`solver_info["large_scale"]["gap"]` is the proven gap; a `FeasibleFound` status there means the gap
|
||||
is above `mip_rel_gap`, not that a time limit was hit.
|
||||
* **`risk_drift="model"` (the default) measures the CVaR on scenarios that include the expected
|
||||
return.** At long horizons the drift grows with *h* while the spread grows with *√h*, so a bullish
|
||||
forecast can push the whole loss tail above zero and the risk limit stops binding — the optimiser
|
||||
|
||||
@ -205,6 +205,11 @@ class OptimizeRequest(BaseModel):
|
||||
"cannot relax the risk limit (the usual market-risk convention)")
|
||||
row_generation: str = Field("auto", pattern="^(auto|on|off)$",
|
||||
description="solve on the tail rows only (large scenario sets)")
|
||||
strategy: str = Field("auto", pattern="^(auto|milp|large_scale)$",
|
||||
description="'large_scale': relaxation, tail rows, screening for cardinality, "
|
||||
"round-and-repair certified against the relaxation bound (large universes "
|
||||
"or scenario sets); 'milp': the integer MILP; 'auto': large_scale from "
|
||||
"200 candidates or 32 768 scenarios on")
|
||||
covariance: dict = Field(default_factory=lambda: {"method": "ewm", "lam": 0.99, "shrinkage": "auto"})
|
||||
expected_returns: str = Field("black_litterman", pattern="^(black_litterman|historical|zero)$")
|
||||
views: list[ViewSpec] = Field(default_factory=list)
|
||||
|
||||
@ -261,7 +261,8 @@ def cmd_optimize(args) -> int:
|
||||
"scenario_model": args.scenario_model, "n_scenarios": args.n_scenarios, "seed": args.seed,
|
||||
"solver": args.solver, "time_limit": args.time_limit, "risk_drift": args.risk_drift,
|
||||
"expected_returns": args.expected_returns, "importance_sampling": args.importance_sampling,
|
||||
"row_generation": args.row_generation, "freq": args.freq, "symbols": args.symbols}
|
||||
"row_generation": args.row_generation, "strategy": args.strategy, "freq": args.freq,
|
||||
"symbols": args.symbols}
|
||||
for k, v in flags.items(): # the request file wins over a flag default
|
||||
if v is not None and k not in body:
|
||||
body[k] = v
|
||||
@ -621,6 +622,8 @@ def build_parser() -> argparse.ArgumentParser:
|
||||
"cannot relax the risk limit")
|
||||
op.add_argument("--importance-sampling", default="auto", choices=["auto", "market", "adaptive", "off"])
|
||||
op.add_argument("--row-generation", default="auto", choices=["auto", "on", "off"])
|
||||
op.add_argument("--strategy", default=None, choices=["auto", "milp", "large_scale"],
|
||||
help="large_scale: relaxation, tail rows, screening, round-and-repair (large universes)")
|
||||
op.add_argument("--max-weight", type=float, default=None, help="cap per stock (default 0.10)")
|
||||
op.add_argument("--max-sector-weight", type=float, default=None, help="cap per GICS sector (default 0.25)")
|
||||
op.add_argument("--view", action="append", metavar="SYM:SENTIMENT[:CONF[:STRENGTH]]", help="repeatable")
|
||||
|
||||
@ -167,6 +167,7 @@ def cmd_optimize(args) -> int:
|
||||
kwargs.setdefault("sampling", args.sampling)
|
||||
kwargs.setdefault("importance_sampling", args.importance_sampling)
|
||||
kwargs.setdefault("row_generation", args.row_generation)
|
||||
kwargs.setdefault("strategy", args.strategy)
|
||||
kwargs.setdefault("solver", args.solver)
|
||||
if args.holdings:
|
||||
kwargs["holdings"] = json.loads(args.holdings)
|
||||
@ -261,6 +262,8 @@ def build_parser() -> argparse.ArgumentParser:
|
||||
help="'market': tilt towards the tail of the cap-weighted portfolio; "
|
||||
"'adaptive': towards a first-pass solution")
|
||||
op.add_argument("--row-generation", default="auto", choices=["auto", "on", "off"])
|
||||
op.add_argument("--strategy", default="auto", choices=["auto", "milp", "large_scale"],
|
||||
help="large_scale: relaxation, tail rows, screening, round-and-repair (large universes)")
|
||||
op.add_argument("--solver", default="auto", choices=["auto", "cuopt", "highs"])
|
||||
op.add_argument("--snapshot", default=None)
|
||||
op.add_argument("--json", action="store_true", help="print the full result as JSON")
|
||||
|
||||
@ -57,7 +57,10 @@ class Settings(BaseSettings):
|
||||
description="cuOpt REST endpoint. The default is a cuOpt container on the same machine (the "
|
||||
"standalone DGX Spark setup, deploy/apptainer/); point it elsewhere for a remote one.")
|
||||
cuopt_time_limit: float = Field(default=120.0, description="MILP time limit in seconds.")
|
||||
cuopt_poll_interval: float = 1.0
|
||||
cuopt_poll_interval: float = Field(
|
||||
default=0.05, description="Seconds between polls of a submitted cuOpt job. The server has no blocking "
|
||||
"call, so this is a latency floor per solve (measured: 0.05 s per call at "
|
||||
"0.02 s, a full second at the old default of 1 s).")
|
||||
cuopt_timeout: float = Field(default=900.0, description="Max seconds to wait for a solution.")
|
||||
cuopt_request_timeout: float = Field(default=60.0, description="Seconds per HTTP call to cuOpt.")
|
||||
cuopt_upload_seconds_per_mb: float = Field(
|
||||
@ -68,6 +71,16 @@ class Settings(BaseSettings):
|
||||
default=True, description="Fall back to the local HiGHS MILP solver (scipy) when cuOpt is unreachable."
|
||||
)
|
||||
cuopt_threads: int | None = Field(default=None, description="CPU threads of the cuOpt MIP solver (None: server default).")
|
||||
cuopt_presolve: bool = Field(
|
||||
default=False, description="cuOpt MILP presolve. Off by default: with it on, cuOpt 26.08 starts an "
|
||||
"asynchronous GPU heuristic and waits for it, so a sub-second MILP takes "
|
||||
"5-25 s (measured 2026-11); off, the same models return in 0.5-2 s. Off, "
|
||||
"cuOpt aborted on 3 of 16 CVaR MILPs on one host (an error reply over REST; "
|
||||
"solve_with_fallback then solves on HiGHS) -- set True if that happens often.")
|
||||
cuopt_lp_method: int | None = Field(
|
||||
default=2, description="cuOpt method for pure LPs (no integer variables): 2 = dual simplex, the fastest "
|
||||
"exact method on the CVaR relaxations. None = cuOpt's concurrent default, which "
|
||||
"returns PDLP's inexact point (violated scenario rows) when PDLP finishes first.")
|
||||
cuopt_msgpack: bool = Field(default=True, description="Send cuOpt requests as msgpack when the package is available.")
|
||||
|
||||
# --- MILP tolerances (reproducibility) ---------------------------------------------------
|
||||
|
||||
372
quantcore/optimization/large_scale.py
Normal file
372
quantcore/optimization/large_scale.py
Normal file
@ -0,0 +1,372 @@
|
||||
"""Large integer mean-CVaR problems: relax, generate tail rows, round and repair, screen for cardinality.
|
||||
|
||||
Measured on the solver benchmark (GPU_Optimization/cuopt_benchmark, November 2026) the integer
|
||||
mean-CVaR MILP has three properties a general MILP solver cannot exploit on its own:
|
||||
|
||||
* only the tail scenarios bind (:mod:`quantcore.optimization.rowgen`): a working set of a few
|
||||
hundred rows replaces tens of thousands;
|
||||
* the integrality is weak -- share counts in the hundreds round with a loss far below the MIP gap,
|
||||
and HiGHS closes every CVaR MILP at the root node: the branch-and-bound machinery is paid for
|
||||
on every call and almost never used;
|
||||
* the cardinality constraint is the only hard combinatorics, and it is far cheaper on a screened
|
||||
subset of candidates than on the whole universe.
|
||||
|
||||
So this driver spends its time on a few small LPs instead of one big MILP:
|
||||
|
||||
1. **relaxation by row generation** -- the LP relaxation (continuous shares, binaries relaxed) is
|
||||
solved on a working set of scenario rows, the losses of all scenarios are evaluated at the
|
||||
continuous solution, violated rows are added, until none is. The converged relaxation is a
|
||||
relaxation of the full MILP, so its objective is a valid bound;
|
||||
2. **screening** (cardinality only) -- the names are ranked by relaxed exposure and the integer
|
||||
stage runs on the top ``max(screen_min, screen_multiple * K)`` of them (held names always
|
||||
stay in); the gap is still measured against the unscreened bound, so the heuristic is bounded;
|
||||
3. **round and repair** -- the integer stage is the restricted model with every share count boxed
|
||||
between the floor and the ceiling of its relaxed value: a MILP whose relaxation is tight and
|
||||
that HiGHS solves at the root in a fraction of a second; its solution is checked against all
|
||||
scenarios and the loop of step 1 repeats on it until no row is violated;
|
||||
4. **certify or escalate** -- the gap between the integer objective and the bound is known without
|
||||
any branching. Within ``mip_gap`` the result is returned as ``Optimal``; otherwise (or when the
|
||||
boxed MILP is infeasible) the ordinary integer row generation runs, starting from the working
|
||||
set already found.
|
||||
|
||||
Solver selection is by size: ``solve_small`` (HiGHS in the service) takes every model below
|
||||
``large_rows`` rows, ``solve_large`` (cuOpt, dual simplex for the LPs) the rest. With row
|
||||
generation the working sets stay small, so the large solver is reached only for very large tails.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import logging
|
||||
import time
|
||||
from collections.abc import Callable
|
||||
from dataclasses import dataclass, field, replace
|
||||
|
||||
import numpy as np
|
||||
|
||||
from quantcore.optimization.problem import MeanCVaRProblem, MILPModel, ProblemSpec
|
||||
from quantcore.optimization.rowgen import (
|
||||
reference_exposure,
|
||||
restricted_spec,
|
||||
scenario_losses,
|
||||
solve_by_row_generation,
|
||||
)
|
||||
from quantcore.optimization.solvers import SolverResult
|
||||
|
||||
log = logging.getLogger(__name__)
|
||||
|
||||
SolveFn = Callable[[MILPModel], SolverResult]
|
||||
|
||||
|
||||
@dataclass
|
||||
class LargeScaleInfo:
|
||||
"""What the driver did, for the result's ``solver_info``."""
|
||||
|
||||
calls: list[dict] = field(default_factory=list) # one entry per solver call
|
||||
bound: float | None = None # objective of the converged relaxation: a valid bound for the MILP
|
||||
relaxation_converged: bool = False
|
||||
relaxation_iterations: int = 0
|
||||
working_set: int = 0 # scenario rows in the final working set
|
||||
screened_to: int | None = None # names kept for the integer stage (cardinality only)
|
||||
integer_stage: str = "" # 'repair' | 'milp' | 'none'
|
||||
repair_iterations: int = 0
|
||||
objective: float | None = None
|
||||
gap: float | None = None # relative to the bound
|
||||
certified: bool = False # gap <= tolerance against the relaxation bound
|
||||
escalated: bool = False # the ordinary integer row generation was run as well
|
||||
total_time: float = 0.0
|
||||
|
||||
def as_dict(self) -> dict:
|
||||
by_stage: dict[str, dict] = {}
|
||||
for c in self.calls:
|
||||
s = by_stage.setdefault(c["stage"], {"calls": 0, "time_s": 0.0, "max_rows": 0})
|
||||
s["calls"] += 1
|
||||
s["time_s"] = round(s["time_s"] + c["time"], 3)
|
||||
s["max_rows"] = max(s["max_rows"], c["rows"])
|
||||
return {"bound": self.bound, "relaxation_converged": self.relaxation_converged,
|
||||
"relaxation_iterations": self.relaxation_iterations, "working_set": self.working_set,
|
||||
"screened_to": self.screened_to, "integer_stage": self.integer_stage,
|
||||
"repair_iterations": self.repair_iterations, "objective": self.objective, "gap": self.gap,
|
||||
"certified": self.certified, "escalated": self.escalated, "stages": by_stage,
|
||||
"solver_calls": len(self.calls), "total_time_s": round(self.total_time, 3)}
|
||||
|
||||
|
||||
# -- helpers ------------------------------------------------------------------------------------
|
||||
|
||||
def relaxed(model: MILPModel) -> MILPModel:
|
||||
"""The LP relaxation: every variable continuous."""
|
||||
return replace(model, var_types=["C"] * model.n_vars)
|
||||
|
||||
|
||||
def exposure_and_cost(problem: MeanCVaRProblem, x: np.ndarray) -> tuple[np.ndarray, float]:
|
||||
"""Money exposure per stock and total transaction cost of a solution vector, fractional or not
|
||||
(:meth:`MeanCVaRProblem.decode` rounds, which is wrong for a relaxed solution)."""
|
||||
value = x[problem.i_n] * problem.spec.prices
|
||||
cost = 0.0
|
||||
for idx, coef in problem._cost_terms():
|
||||
cost += float(coef @ x[idx])
|
||||
return value, cost
|
||||
|
||||
|
||||
def rounded(problem: MeanCVaRProblem, x: np.ndarray) -> np.ndarray:
|
||||
"""Integer variables snapped to the nearest integer (a MILP solver's integers carry 1e-6 noise)."""
|
||||
x = np.array(x, dtype=float)
|
||||
for idx in (problem.i_n, problem.i_b, problem.i_s, problem.i_y, problem.i_t):
|
||||
if idx is not None:
|
||||
x[idx] = np.rint(x[idx])
|
||||
return x
|
||||
|
||||
|
||||
def objective_of(problem: MeanCVaRProblem, x: np.ndarray) -> float:
|
||||
"""The model's objective value of ``x`` (``c'x``), without building the constraint matrix."""
|
||||
sp_ = problem.spec
|
||||
value, cost = exposure_and_cost(problem, x)
|
||||
pnl = float(sp_.mean_returns @ value) - cost
|
||||
if sp_.objective == "max_return":
|
||||
return pnl
|
||||
cvar_idx, cvar_coef = problem._cvar_expr()
|
||||
cvar = float(cvar_coef @ x[cvar_idx])
|
||||
return cvar if sp_.objective == "min_cvar" else pnl - sp_.risk_aversion * cvar
|
||||
|
||||
|
||||
def relative_gap(bound: float, objective: float, maximize: bool) -> float:
|
||||
g = (bound - objective) if maximize else (objective - bound)
|
||||
return max(float(g) / max(1.0, abs(bound)), 0.0)
|
||||
|
||||
|
||||
def _dispatcher(solve_small: SolveFn, solve_large: SolveFn | None, large_rows: int, info: LargeScaleInfo):
|
||||
def solve(model: MILPModel, stage: str) -> SolverResult:
|
||||
large = solve_large is not None and model.n_rows >= large_rows
|
||||
t0 = time.monotonic()
|
||||
res = (solve_large if large else solve_small)(model)
|
||||
info.calls.append({"stage": stage, "rows": model.n_rows, "vars": model.n_vars, "solver": res.solver,
|
||||
"status": res.status, "time": round(time.monotonic() - t0, 3)})
|
||||
return res
|
||||
return solve
|
||||
|
||||
|
||||
def _violated(spec: ProblemSpec, work: np.ndarray, value: np.ndarray, cost: float, alpha: float, tol: float
|
||||
) -> np.ndarray:
|
||||
scale = max(abs(spec.cvar_limit or 0.0), spec.wealth * 1e-3, 1.0)
|
||||
losses = scenario_losses(spec, value, cost)
|
||||
return np.setdiff1d(np.where(losses > alpha + tol * scale)[0], work)
|
||||
|
||||
|
||||
def _initial_work(spec: ProblemSpec, initial_multiple: float) -> np.ndarray:
|
||||
S = spec.n_scenarios
|
||||
losses0 = scenario_losses(spec, reference_exposure(spec))
|
||||
k0 = int(min(S, max(np.ceil(initial_multiple * (1.0 - spec.confidence) * S), 2)))
|
||||
return np.sort(np.argsort(losses0)[-k0:])
|
||||
|
||||
|
||||
# -- stage 1: the relaxation by row generation ---------------------------------------------------
|
||||
|
||||
def relax_by_row_generation(spec: ProblemSpec, solve, work: np.ndarray, max_iter: int = 30, tol: float = 1e-7,
|
||||
deadline: float | None = None
|
||||
) -> tuple[SolverResult, MeanCVaRProblem, np.ndarray, bool, int]:
|
||||
"""Tail-row generation on the LP relaxation, evaluated at the *continuous* solution. Returns
|
||||
the last result (restricted layout), its problem, the working set, convergence and the pass count."""
|
||||
result, problem_r = None, None
|
||||
for it in range(1, max_iter + 1):
|
||||
problem_r = MeanCVaRProblem(restricted_spec(spec, work))
|
||||
result = solve(relaxed(problem_r.build()), "relaxation")
|
||||
if not result.ok:
|
||||
return result, problem_r, work, False, it
|
||||
value, cost = exposure_and_cost(problem_r, result.x)
|
||||
viol = _violated(spec, work, value, cost, float(result.x[problem_r.i_alpha]), tol)
|
||||
log.info("large scale: relaxation pass %d, working set %d of %d, %d violated", it, len(work),
|
||||
spec.n_scenarios, len(viol))
|
||||
if len(viol) == 0:
|
||||
return result, problem_r, work, True, it
|
||||
if deadline is not None and time.monotonic() >= deadline:
|
||||
return result, problem_r, work, False, it
|
||||
work = np.union1d(work, viol)
|
||||
return result, problem_r, work, False, max_iter
|
||||
|
||||
|
||||
# -- stage 2: screening for the cardinality constraint ---------------------------------------------
|
||||
|
||||
def screen(spec: ProblemSpec, value: np.ndarray, keep: int) -> tuple[ProblemSpec, np.ndarray]:
|
||||
"""The spec on the ``keep`` names with the largest relaxed exposure (held names always included)."""
|
||||
n = len(spec.symbols)
|
||||
order = np.argsort(-value, kind="stable")
|
||||
idx = set(np.where(spec.holdings > 0)[0].tolist())
|
||||
for i in order:
|
||||
if len(idx) >= keep:
|
||||
break
|
||||
idx.add(int(i))
|
||||
idx = np.array(sorted(idx), dtype=int)
|
||||
if len(idx) >= n:
|
||||
return spec, np.arange(n)
|
||||
sub = replace(spec, symbols=[spec.symbols[i] for i in idx], prices=spec.prices[idx],
|
||||
scenarios=spec.scenarios[:, idx], sectors=[spec.sectors[i] for i in idx], holdings=spec.holdings[idx],
|
||||
expected_returns=spec.mean_returns[idx], probabilities=spec.probabilities)
|
||||
return sub, idx
|
||||
|
||||
|
||||
# -- stage 3: round and repair ---------------------------------------------------------------------
|
||||
|
||||
def round_and_repair(spec: ProblemSpec, x_rel: np.ndarray, problem_rel: MeanCVaRProblem, work: np.ndarray, solve,
|
||||
max_iter: int = 30, tol: float = 1e-7, deadline: float | None = None
|
||||
) -> tuple[SolverResult | None, MeanCVaRProblem | None, np.ndarray, int]:
|
||||
"""The restricted MILP with every share count boxed between floor and ceiling of its relaxed
|
||||
value, re-solved with violated scenario rows added until none is. ``None`` when the boxed model
|
||||
is infeasible or the solver failed (the caller then escalates)."""
|
||||
n_star = x_rel[problem_rel.i_n]
|
||||
lo = np.floor(n_star + 1e-9)
|
||||
hi = np.ceil(n_star - 1e-9)
|
||||
result, problem_r = None, None
|
||||
for it in range(1, max_iter + 1):
|
||||
problem_r = MeanCVaRProblem(restricted_spec(spec, work))
|
||||
model = problem_r.build()
|
||||
var_lb, var_ub = model.var_lb.copy(), model.var_ub.copy()
|
||||
var_lb[problem_r.i_n] = np.maximum(var_lb[problem_r.i_n], lo)
|
||||
var_ub[problem_r.i_n] = np.minimum(var_ub[problem_r.i_n], hi)
|
||||
var_lb[problem_r.i_n] = np.minimum(var_lb[problem_r.i_n], var_ub[problem_r.i_n])
|
||||
result = solve(replace(model, var_lb=var_lb, var_ub=var_ub), "repair")
|
||||
if not result.ok:
|
||||
return None, problem_r, work, it
|
||||
x = rounded(problem_r, result.x)
|
||||
result = replace(result, x=x, objective=float(model.c @ x + model.offset))
|
||||
value, cost = exposure_and_cost(problem_r, x)
|
||||
viol = _violated(spec, work, value, cost, float(x[problem_r.i_alpha]), tol)
|
||||
log.info("large scale: repair pass %d, working set %d, %d violated", it, len(work), len(viol))
|
||||
if len(viol) == 0:
|
||||
return result, problem_r, work, it
|
||||
if deadline is not None and time.monotonic() >= deadline:
|
||||
return None, problem_r, work, it
|
||||
work = np.union1d(work, viol)
|
||||
return None, problem_r, work, max_iter
|
||||
|
||||
|
||||
# -- embedding into the full layout ---------------------------------------------------------------
|
||||
|
||||
def embed(problem_full: MeanCVaRProblem, problem_sub: MeanCVaRProblem, x_sub: np.ndarray, idx: np.ndarray
|
||||
) -> np.ndarray:
|
||||
"""A restricted (and possibly screened) solution in the full problem's layout: ``u_s = max(L_s -
|
||||
alpha, 0)`` for every scenario, zeros for the screened-out names."""
|
||||
x = np.zeros(problem_full.n_vars)
|
||||
x[problem_full.i_n[idx]] = x_sub[problem_sub.i_n]
|
||||
if problem_full.rebalancing:
|
||||
x[problem_full.i_b[idx]] = x_sub[problem_sub.i_b]
|
||||
x[problem_full.i_s[idx]] = x_sub[problem_sub.i_s]
|
||||
alpha = float(x_sub[problem_sub.i_alpha])
|
||||
x[problem_full.i_alpha] = alpha
|
||||
value_sub, cost = exposure_and_cost(problem_sub, x_sub)
|
||||
value = np.zeros(problem_full.n)
|
||||
value[idx] = value_sub
|
||||
x[problem_full.i_u] = np.maximum(scenario_losses(problem_full.spec, value, cost) - alpha, 0.0)
|
||||
if problem_full.use_pos_bin:
|
||||
x[problem_full.i_y[idx]] = x_sub[problem_sub.i_y]
|
||||
if problem_full.use_trade_bin:
|
||||
x[problem_full.i_t[idx]] = x_sub[problem_sub.i_t]
|
||||
return x
|
||||
|
||||
|
||||
# -- the driver -----------------------------------------------------------------------------------
|
||||
|
||||
def solve_large_scale(spec: ProblemSpec, solve_small: SolveFn, solve_large: SolveFn | None = None, *,
|
||||
large_rows: int = 16_000, mip_gap: float = 1e-4, escalate: bool = True,
|
||||
screen_multiple: int = 4, screen_min: int = 60, initial_multiple: float = 2.0,
|
||||
max_iter: int = 30, tol: float = 1e-7, deadline: float | None = None
|
||||
) -> tuple[SolverResult, MeanCVaRProblem, LargeScaleInfo]:
|
||||
"""Solve the scenario MILP of ``spec`` by relaxation, tail-row generation, screening and rounding.
|
||||
|
||||
``solve_small`` / ``solve_large`` map a :class:`MILPModel` to a :class:`SolverResult`; models with
|
||||
at least ``large_rows`` rows go to ``solve_large`` when it is given. Returns the result in the
|
||||
*full* problem's layout (``MeanCVaRProblem(spec).decode`` works on it), the full problem and the
|
||||
driver's log. The result's ``mip_gap`` is the gap against the relaxation bound; its status is
|
||||
``Optimal`` when that gap is within ``mip_gap``.
|
||||
|
||||
``deadline`` is a :func:`time.monotonic` timestamp for the whole call; the solve functions are
|
||||
expected to carry their own time limits.
|
||||
"""
|
||||
if spec.is_gaussian:
|
||||
raise ValueError("the large-scale driver applies to scenario problems (use gaussian_exact)")
|
||||
t0 = time.monotonic()
|
||||
info = LargeScaleInfo()
|
||||
solve = _dispatcher(solve_small, solve_large, large_rows, info)
|
||||
full = MeanCVaRProblem(spec)
|
||||
maximize = spec.objective != "min_cvar"
|
||||
|
||||
def finish(result: SolverResult, problem_sub: MeanCVaRProblem, idx: np.ndarray, status: str | None = None,
|
||||
message: str = "") -> tuple[SolverResult, MeanCVaRProblem, LargeScaleInfo]:
|
||||
x_full = embed(full, problem_sub, result.x, idx)
|
||||
info.total_time = time.monotonic() - t0
|
||||
if info.bound is not None and result.objective is not None:
|
||||
info.gap = relative_gap(info.bound, result.objective, maximize)
|
||||
info.certified = info.gap <= mip_gap
|
||||
info.objective = result.objective
|
||||
if status is None:
|
||||
status = "Optimal" if info.certified else ("FeasibleFound" if result.status == "Optimal" else result.status)
|
||||
msg = (f"large scale: {message}; bound {info.bound:.6g}, gap {info.gap:.2e}" if info.gap is not None
|
||||
else f"large scale: {message}")
|
||||
return (replace(result, x=x_full, status=status, mip_gap=info.gap, message=msg,
|
||||
raw={**result.raw, "large_scale": info.as_dict()}),
|
||||
full, info)
|
||||
|
||||
# 1. relaxation on the full universe
|
||||
work = _initial_work(spec, initial_multiple)
|
||||
rel, problem_rel, work, converged, passes = relax_by_row_generation(spec, solve, work, max_iter, tol, deadline)
|
||||
info.relaxation_iterations = passes
|
||||
info.relaxation_converged = converged
|
||||
info.working_set = len(work)
|
||||
if not rel.ok:
|
||||
info.total_time = time.monotonic() - t0
|
||||
return replace(rel, message=f"large scale: relaxation failed ({rel.status} {rel.message})".strip()), full, info
|
||||
if converged:
|
||||
info.bound = float(rel.objective)
|
||||
value_rel, _ = exposure_and_cost(problem_rel, rel.x)
|
||||
|
||||
# 2. screening for a cardinality constraint
|
||||
sub_spec, idx = spec, np.arange(len(spec.symbols))
|
||||
k = spec.constraints.max_positions
|
||||
if k is not None:
|
||||
keep = int(min(len(spec.symbols), max(screen_min, screen_multiple * int(k))))
|
||||
sub_spec, idx = screen(spec, value_rel, keep)
|
||||
if len(idx) < len(spec.symbols):
|
||||
info.screened_to = len(idx)
|
||||
rel, problem_rel, work, _, passes = relax_by_row_generation(sub_spec, solve, work, max_iter, tol,
|
||||
deadline)
|
||||
info.relaxation_iterations += passes
|
||||
info.working_set = len(work)
|
||||
if not rel.ok:
|
||||
info.total_time = time.monotonic() - t0
|
||||
return replace(rel, message="large scale: relaxation on the screened universe failed"), full, info
|
||||
|
||||
# 3. round and repair
|
||||
rep, problem_rep, work, it = round_and_repair(sub_spec, rel.x, problem_rel, work, solve, max_iter, tol, deadline)
|
||||
info.repair_iterations = it
|
||||
info.working_set = len(work)
|
||||
if rep is not None:
|
||||
info.integer_stage = "repair"
|
||||
gap = relative_gap(info.bound, rep.objective, maximize) if info.bound is not None else None
|
||||
if gap is not None and gap <= mip_gap:
|
||||
return finish(rep, problem_rep, idx, message=f"rounded within the gap after {it} repair passes")
|
||||
if not escalate or (deadline is not None and time.monotonic() >= deadline):
|
||||
return finish(rep, problem_rep, idx,
|
||||
message=f"rounded; gap above tolerance, not escalated ({it} repair passes)")
|
||||
|
||||
# 4. the ordinary integer row generation, from the working set already found
|
||||
info.escalated = True
|
||||
info.integer_stage = "milp"
|
||||
res, _, rg = solve_by_row_generation(sub_spec, lambda m: solve(m, "milp"), initial_multiple, max_iter, tol,
|
||||
deadline, initial_work=work)
|
||||
info.working_set = max(info.working_set, rg.working_set_sizes[-1] if rg.working_set_sizes else 0)
|
||||
if res.ok:
|
||||
# the MILP returns in sub_spec's full layout: reduce to a restricted problem for embedding
|
||||
problem_r = MeanCVaRProblem(restricted_spec(sub_spec, np.arange(sub_spec.n_scenarios)))
|
||||
x = rounded(problem_r, res.x)
|
||||
res = replace(res, x=x, objective=objective_of(problem_r, x))
|
||||
better = rep is None or (res.objective >= rep.objective if maximize else res.objective <= rep.objective)
|
||||
if better:
|
||||
status = res.status if (info.bound is None or relative_gap(info.bound, res.objective, maximize) > mip_gap) \
|
||||
else "Optimal"
|
||||
return finish(res, problem_r, idx, status=status,
|
||||
message=f"integer row generation ({rg.iterations} passes) after rounding")
|
||||
if rep is not None:
|
||||
return finish(rep, problem_rep, idx, message="rounded; the escalated MILP did not improve on it")
|
||||
info.total_time = time.monotonic() - t0
|
||||
return (replace(res, x=(embed(full, MeanCVaRProblem(restricted_spec(sub_spec, np.arange(sub_spec.n_scenarios))),
|
||||
res.x, idx) if res.x is not None else None),
|
||||
message=f"large scale: no integer solution ({res.status} {res.message})".strip(),
|
||||
raw={**res.raw, "large_scale": info.as_dict()}), full, info)
|
||||
@ -88,7 +88,7 @@ def expand_solution(problem_r: MeanCVaRProblem, x_r: np.ndarray, problem_full: M
|
||||
|
||||
|
||||
def solve_by_row_generation(spec: ProblemSpec, solve: SolveFn, initial_multiple: float = 2.0, max_iter: int = 30,
|
||||
tol: float = 1e-7, deadline: float | None = None
|
||||
tol: float = 1e-7, deadline: float | None = None, initial_work: np.ndarray | None = None
|
||||
) -> tuple[SolverResult, MeanCVaRProblem, RowGenerationInfo]:
|
||||
"""Solve the scenario MILP of ``spec`` by tail-row generation.
|
||||
|
||||
@ -99,6 +99,8 @@ def solve_by_row_generation(spec: ProblemSpec, solve: SolveFn, initial_multiple:
|
||||
``deadline`` is a :func:`time.monotonic` timestamp: the *request's* time budget. Each pass is a
|
||||
MILP of its own, so without it a limit of ``T`` seconds would allow ``T`` per pass instead of
|
||||
``T`` in total. Once it passes, the current (feasible, not proven optimal) solution is returned.
|
||||
``initial_work`` replaces the reference portfolio's tail as the first working set (the
|
||||
large-scale driver passes the rows its relaxation found).
|
||||
"""
|
||||
if spec.is_gaussian:
|
||||
raise ValueError("row generation applies to scenario problems")
|
||||
@ -106,9 +108,12 @@ def solve_by_row_generation(spec: ProblemSpec, solve: SolveFn, initial_multiple:
|
||||
info = RowGenerationInfo()
|
||||
S = spec.n_scenarios
|
||||
full = MeanCVaRProblem(spec)
|
||||
losses0 = scenario_losses(spec, reference_exposure(spec))
|
||||
k0 = int(min(S, max(np.ceil(initial_multiple * (1.0 - spec.confidence) * S), 2)))
|
||||
work = np.sort(np.argsort(losses0)[-k0:])
|
||||
if initial_work is not None and len(initial_work) > 0:
|
||||
work = np.unique(np.asarray(initial_work, dtype=int))
|
||||
else:
|
||||
losses0 = scenario_losses(spec, reference_exposure(spec))
|
||||
k0 = int(min(S, max(np.ceil(initial_multiple * (1.0 - spec.confidence) * S), 2)))
|
||||
work = np.sort(np.argsort(losses0)[-k0:])
|
||||
scale = max(abs(spec.cvar_limit or 0.0), spec.wealth * 1e-3, 1.0)
|
||||
result: SolverResult | None = None
|
||||
for it in range(1, max_iter + 1):
|
||||
|
||||
@ -13,7 +13,9 @@ Pipeline (deterministic given the seed; the result carries the scenario hash and
|
||||
solution (``importance_sampling``), or -- ``scenario_model='gaussian_exact'`` -- no scenarios
|
||||
at all but the exact Gaussian (mean-variance-equivalent) problem by outer approximation;
|
||||
4. the mean-CVaR MILP (:mod:`quantcore.optimization.problem`) in "new" or "rebalancing" mode,
|
||||
for large scenario sets by tail-row generation (:mod:`quantcore.optimization.rowgen`);
|
||||
for large scenario sets by tail-row generation (:mod:`quantcore.optimization.rowgen`), and for
|
||||
large universes or scenario sets by relaxation, screening and rounding
|
||||
(:mod:`quantcore.optimization.large_scale`, ``strategy``);
|
||||
5. solved on cuOpt (HiGHS fallback), decoded into shares, trades, cash, costs, risk figures.
|
||||
"""
|
||||
|
||||
@ -43,6 +45,7 @@ from quantcore.models.covariance import CovarianceModel
|
||||
from quantcore.models.sampling import next_power_of_two, scenario_hash
|
||||
from quantcore.models.simulation import ScenarioSet, TiltSpec, make_simulator
|
||||
from quantcore.optimization.gaussian_exact import solve_gaussian_exact
|
||||
from quantcore.optimization.large_scale import solve_large_scale
|
||||
from quantcore.optimization.problem import (
|
||||
MeanCVaRProblem,
|
||||
OptimizationConstraints,
|
||||
@ -58,6 +61,12 @@ log = logging.getLogger(__name__)
|
||||
SCENARIO_MODELS = ("gbm", "empirical", "gaussian_exact")
|
||||
IMPORTANCE_SAMPLING_MODES = ("auto", "market", "adaptive", "off")
|
||||
ROW_GENERATION_MIN_SCENARIOS = 8192 # 'auto': row generation from this size on
|
||||
# strategy 'auto': the relax / round-and-repair driver (quantcore.optimization.large_scale) from this
|
||||
# many candidates or scenarios on; below, the integer MILP as before. Measured (cuopt_benchmark, 2026-11):
|
||||
# HiGHS closes every CVaR MILP at the root, so the driver loses nothing and saves the tree search.
|
||||
LARGE_SCALE_MIN_ASSETS = 200
|
||||
LARGE_SCALE_MIN_SCENARIOS = 32768
|
||||
LARGE_SCALE_CUOPT_ROWS = 16_000 # the driver sends models with at least this many rows to cuOpt
|
||||
IMPORTANCE_SAMPLING_MIN_TAIL = 128 # 'auto': tilt when fewer than this many scenarios would be in the tail
|
||||
STAGE1_SCENARIOS = 2048 # 'adaptive': first-pass solve that provides the tilt direction
|
||||
TILT_ALIGNMENT_WARN = 0.3 # warn below this correlation between the solution's loss and the tilt's
|
||||
@ -119,6 +128,11 @@ class OptimizationRequest:
|
||||
# objective keeps the expected return, which is the usual market-risk convention.
|
||||
risk_drift: str = "model" # 'model' | 'zero'
|
||||
row_generation: str = "auto" # 'auto' | 'on' | 'off': solve on the tail rows only
|
||||
# 'large_scale': relaxation, tail rows, screening for cardinality, round and repair, certified
|
||||
# against the relaxation bound (quantcore.optimization.large_scale); 'milp': the integer MILP
|
||||
# (with row generation per `row_generation`); 'auto': 'large_scale' from LARGE_SCALE_MIN_ASSETS
|
||||
# candidates or LARGE_SCALE_MIN_SCENARIOS scenarios on
|
||||
strategy: str = "auto" # 'auto' | 'milp' | 'large_scale'
|
||||
covariance: dict = field(default_factory=lambda: {"method": "ewm", "lam": 0.99, "shrinkage": "auto"})
|
||||
expected_returns: str = "black_litterman" # 'black_litterman' | 'historical' | 'zero'
|
||||
views: list[dict] = field(default_factory=list)
|
||||
@ -143,6 +157,8 @@ class OptimizationRequest:
|
||||
raise ValueError(f"importance_sampling must be one of {IMPORTANCE_SAMPLING_MODES}")
|
||||
if self.row_generation not in ("auto", "on", "off"):
|
||||
raise ValueError("row_generation must be 'auto', 'on' or 'off'")
|
||||
if self.strategy not in ("auto", "milp", "large_scale"):
|
||||
raise ValueError("strategy must be 'auto', 'milp' or 'large_scale'")
|
||||
if self.risk_drift not in ("model", "zero"):
|
||||
raise ValueError("risk_drift must be 'model' or 'zero'")
|
||||
if isinstance(self.constraints, dict):
|
||||
@ -169,6 +185,12 @@ class OptimizationRequest:
|
||||
return True
|
||||
return self.row_generation == "auto" and n_scenarios >= ROW_GENERATION_MIN_SCENARIOS
|
||||
|
||||
def uses_large_scale(self, n_assets: int, n_scenarios: int) -> bool:
|
||||
if self.strategy == "large_scale":
|
||||
return True
|
||||
return self.strategy == "auto" and (n_assets >= LARGE_SCALE_MIN_ASSETS
|
||||
or n_scenarios >= LARGE_SCALE_MIN_SCENARIOS)
|
||||
|
||||
def tilt_mode(self, n_scenarios: int) -> str:
|
||||
"""``'market'``, ``'adaptive'`` or ``'off'`` for this request.
|
||||
|
||||
@ -487,6 +509,14 @@ class OptimizationService:
|
||||
if weighted_cvar_objective and req.solver == "auto":
|
||||
info["solver_routing"] = "highs: CVaR objective with probability-weighted scenarios"
|
||||
solve = self._solve_fn(req, prefer_local=weighted_cvar_objective, deadline=deadline)
|
||||
if req.uses_large_scale(len(spec.symbols), spec.n_scenarios):
|
||||
# small working-set models on HiGHS locally (unless cuOpt was asked for explicitly), the
|
||||
# large ones on the configured solver; the CVaR-objective case stays off cuOpt entirely
|
||||
small = self._solve_fn(req, prefer_local=True, deadline=deadline)
|
||||
large = None if weighted_cvar_objective else solve
|
||||
sol, problem, ls = solve_large_scale(spec, small, large, large_rows=LARGE_SCALE_CUOPT_ROWS,
|
||||
mip_gap=self.settings.mip_rel_gap, deadline=deadline)
|
||||
return sol, problem, {**info, "large_scale": ls.as_dict()}
|
||||
if req.uses_row_generation(spec.n_scenarios):
|
||||
sol, problem, rg = solve_by_row_generation(spec, solve, deadline=deadline)
|
||||
return sol, problem, {**info, "row_generation": rg.as_dict()}
|
||||
@ -562,10 +592,14 @@ class OptimizationService:
|
||||
# FeasibleFound can come from the solver's own limit *or* from an iteration cap in the
|
||||
# outer approximation / row generation, which is not a time limit at all
|
||||
oa, rg = extra_info.get("outer_approximation") or {}, extra_info.get("row_generation") or {}
|
||||
ls = extra_info.get("large_scale") or {}
|
||||
if oa and not oa.get("converged"):
|
||||
why = f"the Gaussian outer approximation stopped after {oa.get('iterations')} cuts"
|
||||
elif rg and not rg.get("converged"):
|
||||
why = f"the tail-row generation stopped after {rg.get('iterations')} passes"
|
||||
elif ls and not ls.get("certified"):
|
||||
why = (f"the large-scale driver's gap against the relaxation bound is {ls.get('gap'):.2e}"
|
||||
if ls.get("gap") is not None else "the large-scale driver's relaxation did not converge")
|
||||
else:
|
||||
why = "the solver stopped at its limit"
|
||||
warnings.append(f"{why}; solution feasible but not proven optimal (gap {sol.mip_gap})")
|
||||
@ -617,6 +651,7 @@ class OptimizationService:
|
||||
if spec.cvar_limit is not None and d["cvar"] > spec.cvar_limit * (1 + 1e-6):
|
||||
rg = extra_info.get("row_generation") or {}
|
||||
oa = extra_info.get("outer_approximation") or {}
|
||||
ls = extra_info.get("large_scale") or {}
|
||||
over = d["cvar"] / spec.cvar_limit - 1
|
||||
# enough precision for a small overshoot: "by 0.0%" reads as no breach at all
|
||||
amount = f"{d['cvar'] - spec.cvar_limit:,.0f}" if over >= 5e-4 else f"{d['cvar'] - spec.cvar_limit:,.2f}"
|
||||
@ -629,6 +664,10 @@ class OptimizationService:
|
||||
cause = (f": the Gaussian outer approximation stopped after {oa.get('iterations')} cuts "
|
||||
f"without converging, and it bounds the CVaR from below until it does - raise "
|
||||
f"time_limit, or use scenario_model='gbm' for an exactly enforced limit")
|
||||
elif ls and not ls.get("relaxation_converged"):
|
||||
cause = (f": the large-scale driver's row generation stopped before every scenario was covered "
|
||||
f"({ls.get('working_set')} of {spec.n_scenarios} rows) - raise time_limit, or set "
|
||||
f"strategy='milp'")
|
||||
else:
|
||||
cause = (": the solver stopped within its MIP gap rather than at the true optimum - raise "
|
||||
"time_limit, or lower n_scenarios / max_assets")
|
||||
|
||||
@ -73,6 +73,8 @@ def cuopt_solver_config(settings: Settings, time_limit: float, extra_config: dic
|
||||
"tolerances": {"mip_relative_gap": settings.mip_rel_gap, "mip_absolute_gap": settings.mip_abs_gap}}
|
||||
if settings.cuopt_threads:
|
||||
cfg["num_cpu_threads"] = int(settings.cuopt_threads)
|
||||
if not settings.cuopt_presolve:
|
||||
cfg["presolve"] = 0
|
||||
if extra_config:
|
||||
tol = {**cfg["tolerances"], **(extra_config.get("tolerances") or {})}
|
||||
cfg = {**cfg, **extra_config, "tolerances": tol}
|
||||
@ -82,6 +84,8 @@ def cuopt_solver_config(settings: Settings, time_limit: float, extra_config: dic
|
||||
def to_cuopt_payload(model: MILPModel, time_limit: float, extra_config: dict | None = None,
|
||||
settings: Settings | None = None) -> dict:
|
||||
"""The cuOpt LP/MILP request body for a :class:`MILPModel`."""
|
||||
if settings is not None and settings.cuopt_lp_method is not None and "I" not in model.var_types:
|
||||
extra_config = {"method": int(settings.cuopt_lp_method), **(extra_config or {})} # an LP: dual simplex
|
||||
A = model.A.tocsr()
|
||||
payload = {
|
||||
"csr_constraint_matrix": {
|
||||
|
||||
85
scripts/experiments/large_scale.py
Normal file
85
scripts/experiments/large_scale.py
Normal file
@ -0,0 +1,85 @@
|
||||
"""The large-scale driver against the integer MILP on a synthetic universe (no market data needed).
|
||||
|
||||
python scripts/experiments/large_scale.py 300 16384 # N stocks, S scenarios
|
||||
python scripts/experiments/large_scale.py 750 16384,65536 --cardinality 40 --cuopt http://spark-6d4a:5000
|
||||
|
||||
Prints, per size: the full-MILP row generation (quantcore's existing path) and the driver, with wall
|
||||
time, objective, the driver's bound and gap, working sets and solver calls. ``--cuopt URL`` adds the
|
||||
cuOpt server as the large-model solver of the driver (and as the MILP solver of the reference run).
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import argparse
|
||||
import time
|
||||
|
||||
import numpy as np
|
||||
|
||||
from quantcore.config import Settings
|
||||
from quantcore.optimization.large_scale import solve_large_scale
|
||||
from quantcore.optimization.problem import OptimizationConstraints, ProblemSpec
|
||||
from quantcore.optimization.rowgen import solve_by_row_generation
|
||||
from quantcore.optimization.solvers import HighsSolver, solve_with_fallback
|
||||
|
||||
|
||||
def market(N: int, S: int, seed: int = 0, K: int = 6, horizon_days: int = 5):
|
||||
"""A K-factor universe with Student-t(4) scale-mixture scenarios (the benchmark's generator)."""
|
||||
rng = np.random.default_rng(seed)
|
||||
prices = np.exp(rng.uniform(np.log(10.0), np.log(600.0), N))
|
||||
sectors = [f"sector_{s}" for s in rng.integers(0, 10, N)]
|
||||
f_vol = np.concatenate([[0.18], rng.uniform(0.08, 0.20, K - 1)])
|
||||
F = np.diag((f_vol ** 2) * horizon_days / 252.0)
|
||||
B = np.column_stack([rng.normal(1.0, 0.3, N)] + [rng.normal(0.0, 0.5, N) for _ in range(K - 1)])
|
||||
D = (rng.uniform(0.15, 0.40, N) ** 2) * horizon_days / 252.0
|
||||
mu = rng.normal(0.06, 0.08, N) * horizon_days / 252.0
|
||||
rng2 = np.random.default_rng(seed + 7919)
|
||||
z = rng2.standard_normal((S, K)) @ np.linalg.cholesky(F).T
|
||||
e = rng2.standard_normal((S, N)) * np.sqrt(D)
|
||||
r = mu + z @ B.T + e
|
||||
w = np.sqrt(4.0 / rng2.chisquare(4.0, size=(S, 1)))
|
||||
r = mu + (r - mu) * w * np.sqrt(0.5)
|
||||
return prices, sectors, mu, r
|
||||
|
||||
|
||||
def main():
|
||||
ap = argparse.ArgumentParser()
|
||||
ap.add_argument("N", type=int)
|
||||
ap.add_argument("S", help="comma-separated scenario counts")
|
||||
ap.add_argument("--cardinality", type=int, default=None, help="max_positions (adds binaries)")
|
||||
ap.add_argument("--cuopt", default=None, help="cuOpt URL for the large models (default: HiGHS only)")
|
||||
ap.add_argument("--no-reference", action="store_true", help="skip the full-MILP reference run")
|
||||
ap.add_argument("--time-limit", type=float, default=600.0)
|
||||
ap.add_argument("--seed", type=int, default=0)
|
||||
a = ap.parse_args()
|
||||
settings = Settings(cuopt_url=a.cuopt or "http://127.0.0.1:9", mip_rel_gap=1e-4, solver_fallback=True)
|
||||
highs = lambda m: HighsSolver(settings).solve(m, time_limit=a.time_limit)
|
||||
cuopt = (lambda m: solve_with_fallback(m, settings, solver="cuopt", time_limit=a.time_limit)) if a.cuopt else None
|
||||
W = 1_000_000.0
|
||||
print(f"{'N':>4} {'S':>6} {'K':>4} | {'MILP s':>8} {'objective':>12} | {'driver s':>8} {'objective':>12} "
|
||||
f"{'bound':>12} {'gap':>8} {'stage':>7} {'rows':>6} {'calls':>5} {'escal':>5}")
|
||||
for S in (int(x) for x in a.S.split(",")):
|
||||
prices, sectors, mu, scen = market(a.N, S, a.seed)
|
||||
cons = OptimizationConstraints(max_weight=0.10, max_sector_weight=0.25, max_positions=a.cardinality)
|
||||
spec = ProblemSpec(symbols=[f"S{i}" for i in range(a.N)], prices=prices, scenarios=scen, sectors=sectors,
|
||||
cash=W, cvar_limit=0.05 * W, expected_returns=mu, constraints=cons)
|
||||
ref = "-"
|
||||
ref_obj = float("nan")
|
||||
if not a.no_reference:
|
||||
t0 = time.perf_counter()
|
||||
res, problem, _ = solve_by_row_generation(spec, cuopt or highs)
|
||||
ref = f"{time.perf_counter() - t0:8.1f}"
|
||||
ref_obj = problem.decode(res.x)["expected_pnl"] if res.ok else float("nan")
|
||||
t0 = time.perf_counter()
|
||||
res, problem, info = solve_large_scale(spec, highs, cuopt, mip_gap=1e-4)
|
||||
t_ls = time.perf_counter() - t0
|
||||
obj = problem.decode(res.x)["expected_pnl"] if res.ok else float("nan")
|
||||
print(f"{a.N:>4} {S:>6} {a.cardinality or 0:>4} | {ref:>8} {ref_obj:12.2f} | {t_ls:8.1f} {obj:12.2f} "
|
||||
f"{(info.bound or float('nan')):12.2f} {(info.gap if info.gap is not None else float('nan')):8.1e} "
|
||||
f"{info.integer_stage:>7} {info.working_set:>6} {len(info.calls):>5} {info.escalated!s:>5}", flush=True)
|
||||
for c in info.calls:
|
||||
print(f" {c['stage']:>10} rows={c['rows']:>6} vars={c['vars']:>7} {c['solver']:>6} {c['status']:>13} "
|
||||
f"{c['time']:7.2f}s")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
185
tests/test_large_scale.py
Normal file
185
tests/test_large_scale.py
Normal file
@ -0,0 +1,185 @@
|
||||
"""The large-scale driver (relaxation, tail rows, screening, round and repair) against the full MILP."""
|
||||
|
||||
import numpy as np
|
||||
import pytest
|
||||
|
||||
from quantcore.config import Settings
|
||||
from quantcore.models import sampling as smp
|
||||
from quantcore.optimization.large_scale import (
|
||||
exposure_and_cost,
|
||||
relative_gap,
|
||||
relaxed,
|
||||
solve_large_scale,
|
||||
)
|
||||
from quantcore.optimization.problem import MeanCVaRProblem, OptimizationConstraints, ProblemSpec
|
||||
from quantcore.optimization.rowgen import solve_by_row_generation
|
||||
from quantcore.optimization.service import OptimizationRequest, OptimizationService
|
||||
from quantcore.optimization.solvers import HighsSolver
|
||||
|
||||
SETTINGS = Settings(cuopt_url="http://127.0.0.1:9", mip_rel_gap=1e-6)
|
||||
|
||||
|
||||
def highs(model):
|
||||
return HighsSolver(SETTINGS).solve(model, time_limit=120)
|
||||
|
||||
|
||||
def inputs(seed=0, n=12, S=2048):
|
||||
"""A synthetic universe with moment-matched Sobol' scenarios and a fat-tailed scale mixture."""
|
||||
rng = np.random.default_rng(seed)
|
||||
symbols = [f"S{i}" for i in range(n)]
|
||||
prices = rng.uniform(20, 400, n)
|
||||
vol = rng.uniform(0.02, 0.06, n)
|
||||
corr = 0.4 * np.ones((n, n)) + 0.6 * np.eye(n)
|
||||
cov = corr * np.outer(vol, vol)
|
||||
mean = rng.uniform(0.0, 0.006, n)
|
||||
z = smp.moment_match(smp.standard_normal_sample(S, n, "sobol", seed))
|
||||
w = np.sqrt(4.0 / rng.chisquare(4.0, size=(S, 1))) * np.sqrt(0.5) # Student-t(4) scale mixture
|
||||
scen = mean + (z @ smp.pca_root(cov).T) * w
|
||||
sectors = ["Tech", "Fin", "Health"] * (n // 3) + ["Tech"] * (n % 3)
|
||||
return symbols, prices, scen, sectors
|
||||
|
||||
|
||||
def check_feasible(spec, problem, x, tol=1e-6):
|
||||
"""The returned vector satisfies the full model, and the decoded portfolio the business limits."""
|
||||
m = problem.build()
|
||||
Ax = m.A @ x
|
||||
assert (Ax >= m.row_lb - tol).all() and (Ax <= m.row_ub + tol).all()
|
||||
assert (x >= m.var_lb - tol).all() and (x <= m.var_ub + tol).all()
|
||||
ints = [i for i, t in enumerate(m.var_types) if t == "I"]
|
||||
assert np.abs(x[ints] - np.rint(x[ints])).max() < 1e-9
|
||||
d = problem.decode(x)
|
||||
assert d["cash_end"] >= -tol
|
||||
if spec.cvar_limit is not None:
|
||||
assert d["cvar"] <= spec.cvar_limit * (1 + 1e-7)
|
||||
return d
|
||||
|
||||
|
||||
def test_matches_full_milp_and_certifies_the_gap():
|
||||
# a portfolio of realistic size: share counts in the thousands, so rounding costs far less than the gap
|
||||
symbols, prices, scen, sectors = inputs(seed=2)
|
||||
spec = ProblemSpec(symbols=symbols, prices=prices, scenarios=scen, sectors=sectors, cash=5_000_000.0,
|
||||
cvar_limit=150_000.0, constraints=OptimizationConstraints(max_weight=0.3, max_sector_weight=0.6))
|
||||
full = MeanCVaRProblem(spec)
|
||||
res_full = highs(full.build())
|
||||
res, problem, info = solve_large_scale(spec, highs, mip_gap=1e-4)
|
||||
assert res.ok and res.status == "Optimal" and info.certified and not info.escalated
|
||||
assert info.relaxation_converged and info.working_set < 0.5 * spec.n_scenarios
|
||||
assert info.integer_stage == "repair"
|
||||
d = check_feasible(spec, problem, res.x)
|
||||
d_full = full.decode(res_full.x)
|
||||
# the relaxation bounds the MILP optimum, and the rounded portfolio is within the gap of it
|
||||
assert info.bound >= d_full["expected_pnl"] - 1e-6
|
||||
assert d["expected_pnl"] >= d_full["expected_pnl"] - 1e-4 * abs(info.bound) - 1e-6
|
||||
assert d["expected_pnl"] <= d_full["expected_pnl"] * (1 + 1e-6) + 1e-6
|
||||
assert res.mip_gap == info.gap and info.gap <= 1e-4
|
||||
assert res.raw["large_scale"]["solver_calls"] == len(info.calls)
|
||||
|
||||
|
||||
def test_relaxation_is_a_valid_bound_at_the_continuous_point():
|
||||
symbols, prices, scen, sectors = inputs(seed=5, n=9, S=1024)
|
||||
spec = ProblemSpec(symbols=symbols, prices=prices, scenarios=scen, sectors=sectors, cash=50_000.0,
|
||||
cvar_limit=2_000.0, constraints=OptimizationConstraints(max_weight=0.3, max_sector_weight=0.6))
|
||||
full = MeanCVaRProblem(spec)
|
||||
lp = highs(relaxed(full.build())) # the LP relaxation on all rows
|
||||
res, _, info = solve_large_scale(spec, highs, mip_gap=1e-4)
|
||||
assert res.ok and info.relaxation_converged
|
||||
assert abs(info.bound - lp.objective) <= 1e-6 * max(1.0, abs(lp.objective))
|
||||
value, cost = exposure_and_cost(full, lp.x)
|
||||
assert value.shape == (9,) and cost > 0
|
||||
|
||||
|
||||
def test_cardinality_screening_is_bounded():
|
||||
symbols, prices, scen, sectors = inputs(seed=7, n=30, S=1024)
|
||||
cons = OptimizationConstraints(max_weight=0.2, max_sector_weight=0.6, max_positions=5)
|
||||
spec = ProblemSpec(symbols=symbols, prices=prices, scenarios=scen, sectors=sectors, cash=200_000.0,
|
||||
cvar_limit=6_000.0, constraints=cons)
|
||||
full = MeanCVaRProblem(spec)
|
||||
res_full = highs(full.build())
|
||||
d_full = full.decode(res_full.x)
|
||||
res, problem, info = solve_large_scale(spec, highs, mip_gap=1e-4, screen_min=10, screen_multiple=2)
|
||||
assert res.ok and info.screened_to == 10
|
||||
d = check_feasible(spec, problem, res.x)
|
||||
assert int((d["shares"] > 0).sum()) <= 5
|
||||
# the heuristic never beats the true optimum, and the reported gap covers the real one
|
||||
assert d["expected_pnl"] <= d_full["expected_pnl"] * (1 + 1e-6) + 1e-6
|
||||
true_gap = relative_gap(info.bound, d["expected_pnl"], maximize=True)
|
||||
assert info.gap >= true_gap - 1e-9
|
||||
assert info.gap >= relative_gap(d_full["expected_pnl"], d["expected_pnl"], maximize=True) - 1e-9
|
||||
|
||||
|
||||
def test_min_cvar_objective_and_rebalancing():
|
||||
symbols, prices, scen, sectors = inputs(seed=3, n=9, S=1024)
|
||||
holdings = np.array([40, 0, 10, 0, 5, 0, 0, 8, 0], dtype=float)
|
||||
target = 0.001 * (20_000.0 + holdings @ prices) # a third of what a fully invested portfolio expects
|
||||
spec = ProblemSpec(symbols=symbols, prices=prices, scenarios=scen, sectors=sectors, cash=20_000.0,
|
||||
holdings=holdings, objective="min_cvar", target_pnl=target,
|
||||
constraints=OptimizationConstraints(max_weight=0.3, max_sector_weight=0.6))
|
||||
full = MeanCVaRProblem(spec)
|
||||
res_full = highs(full.build())
|
||||
d_full = full.decode(res_full.x)
|
||||
res, problem, info = solve_large_scale(spec, highs, mip_gap=1e-4)
|
||||
assert res.ok and info.relaxation_converged
|
||||
d = check_feasible(spec, problem, res.x)
|
||||
np.testing.assert_array_equal(d["shares"], holdings + d["buys"] - d["sells"])
|
||||
assert d["expected_pnl"] >= target - 1e-6
|
||||
assert info.bound <= d_full["cvar"] + 1e-6 # a lower bound on the minimal CVaR
|
||||
assert d["cvar"] <= d_full["cvar"] + 1e-4 * max(1.0, abs(info.bound)) + 1e-6
|
||||
|
||||
|
||||
def test_escalates_to_the_integer_milp_when_not_certified():
|
||||
symbols, prices, scen, sectors = inputs(seed=11, n=12, S=1024)
|
||||
spec = ProblemSpec(symbols=symbols, prices=prices, scenarios=scen, sectors=sectors, cash=100_000.0,
|
||||
cvar_limit=3_000.0, constraints=OptimizationConstraints(max_weight=0.3, max_sector_weight=0.6))
|
||||
res_rg, problem_rg, _ = solve_by_row_generation(spec, highs)
|
||||
d_rg = problem_rg.decode(res_rg.x)
|
||||
res, problem, info = solve_large_scale(spec, highs, mip_gap=1e-12) # nothing rounds this well
|
||||
assert res.ok and info.escalated and info.integer_stage == "milp"
|
||||
d = check_feasible(spec, problem, res.x)
|
||||
assert abs(d["expected_pnl"] - d_rg["expected_pnl"]) <= 1e-6 * abs(d_rg["expected_pnl"]) + 1e-6
|
||||
res_no, _, info_no = solve_large_scale(spec, highs, mip_gap=1e-12, escalate=False)
|
||||
assert res_no.ok and not info_no.escalated and res_no.status == "FeasibleFound" and not info_no.certified
|
||||
|
||||
|
||||
def test_solver_selection_by_size():
|
||||
symbols, prices, scen, sectors = inputs(seed=4, n=9, S=512)
|
||||
spec = ProblemSpec(symbols=symbols, prices=prices, scenarios=scen, sectors=sectors, cash=50_000.0,
|
||||
cvar_limit=2_000.0, constraints=OptimizationConstraints(max_weight=0.3, max_sector_weight=0.6))
|
||||
seen = []
|
||||
|
||||
def small(model):
|
||||
seen.append(("small", model.n_rows))
|
||||
return highs(model)
|
||||
|
||||
def large(model):
|
||||
seen.append(("large", model.n_rows))
|
||||
return highs(model)
|
||||
|
||||
res, _, info = solve_large_scale(spec, small, large, large_rows=100)
|
||||
assert res.ok
|
||||
assert all(w == ("large" if rows >= 100 else "small") for w, rows in seen)
|
||||
assert {c["stage"] for c in info.calls} <= {"relaxation", "repair", "milp"}
|
||||
with pytest.raises(ValueError):
|
||||
gspec = ProblemSpec(symbols=symbols[:3], prices=prices[:3], scenarios=None, sectors=sectors[:3], cash=1e4,
|
||||
cvar_limit=500.0, gaussian_mean=np.zeros(3), gaussian_cov=np.eye(3) * 1e-4)
|
||||
solve_large_scale(gspec, highs)
|
||||
|
||||
|
||||
BASE = {"cash": 500_000.0, "cvar_limit": 12_000.0, "solver": "highs", "horizon": "week", "scenario_model": "gbm",
|
||||
"n_scenarios": 2048, "importance_sampling": "off", "costs": {"proportional_buy": 0.0, "proportional_sell": 0.0},
|
||||
"constraints": {"max_weight": 0.15, "max_sector_weight": 0.4}}
|
||||
|
||||
|
||||
def test_service_routes_the_strategy(settings, md):
|
||||
from tests.conftest import SYMBOLS
|
||||
|
||||
svc = OptimizationService(settings, md)
|
||||
milp = svc.optimize(OptimizationRequest(symbols=SYMBOLS, strategy="milp", **BASE))
|
||||
ls = svc.optimize(OptimizationRequest(symbols=SYMBOLS, strategy="large_scale", **BASE))
|
||||
assert "large_scale" not in milp.solver_info and ls.solver_info["large_scale"]["certified"]
|
||||
assert ls.status == "Optimal" and ls.reproducible
|
||||
assert abs(ls.expected_pnl - milp.expected_pnl) <= 1e-4 * abs(milp.expected_pnl) + 1e-6
|
||||
assert ls.cvar <= 12_000.0 * (1 + 1e-7)
|
||||
req = OptimizationRequest(symbols=SYMBOLS, **BASE)
|
||||
assert not req.uses_large_scale(12, 2048) and req.uses_large_scale(200, 2048) and req.uses_large_scale(12, 32768)
|
||||
with pytest.raises(ValueError):
|
||||
OptimizationRequest(symbols=SYMBOLS, strategy="bogus", **BASE)
|
||||
Loading…
x
Reference in New Issue
Block a user