From fac41ffdfcc2379d0ecc0dc8e87a23121a97132f Mon Sep 17 00:00:00 2001 From: herzogflorian Date: Wed, 4 Nov 2026 09:58:32 +0100 Subject: [PATCH] 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 --- README.md | 2 +- docs/ARCHITECTURE.md | 2 +- docs/METHODOLOGY.md | 29 ++ docs/USAGE.md | 6 + quantcore/api/schemas.py | 5 + quantcore/apicli.py | 5 +- quantcore/cli.py | 3 + quantcore/config.py | 15 +- quantcore/optimization/large_scale.py | 372 ++++++++++++++++++++++++++ quantcore/optimization/rowgen.py | 13 +- quantcore/optimization/service.py | 41 ++- quantcore/optimization/solvers.py | 4 + scripts/experiments/large_scale.py | 85 ++++++ tests/test_large_scale.py | 185 +++++++++++++ 14 files changed, 758 insertions(+), 9 deletions(-) create mode 100644 quantcore/optimization/large_scale.py create mode 100644 scripts/experiments/large_scale.py create mode 100644 tests/test_large_scale.py diff --git a/README.md b/README.md index aa8ec8e..9cc64d8 100644 --- a/README.md +++ b/README.md @@ -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 diff --git a/docs/ARCHITECTURE.md b/docs/ARCHITECTURE.md index 524bd16..1938a38 100644 --- a/docs/ARCHITECTURE.md +++ b/docs/ARCHITECTURE.md @@ -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 | diff --git a/docs/METHODOLOGY.md b/docs/METHODOLOGY.md index d323b44..abfb215 100644 --- a/docs/METHODOLOGY.md +++ b/docs/METHODOLOGY.md @@ -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 diff --git a/docs/USAGE.md b/docs/USAGE.md index d79fb8e..99b6978 100644 --- a/docs/USAGE.md +++ b/docs/USAGE.md @@ -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 diff --git a/quantcore/api/schemas.py b/quantcore/api/schemas.py index 8dd3a05..ef86ba1 100644 --- a/quantcore/api/schemas.py +++ b/quantcore/api/schemas.py @@ -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) diff --git a/quantcore/apicli.py b/quantcore/apicli.py index 9be3784..d883e21 100644 --- a/quantcore/apicli.py +++ b/quantcore/apicli.py @@ -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") diff --git a/quantcore/cli.py b/quantcore/cli.py index d1d0c7a..badfc67 100644 --- a/quantcore/cli.py +++ b/quantcore/cli.py @@ -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") diff --git a/quantcore/config.py b/quantcore/config.py index 71401f0..bc18fbc 100644 --- a/quantcore/config.py +++ b/quantcore/config.py @@ -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) --------------------------------------------------- diff --git a/quantcore/optimization/large_scale.py b/quantcore/optimization/large_scale.py new file mode 100644 index 0000000..31ff594 --- /dev/null +++ b/quantcore/optimization/large_scale.py @@ -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) diff --git a/quantcore/optimization/rowgen.py b/quantcore/optimization/rowgen.py index db0348c..8069e22 100644 --- a/quantcore/optimization/rowgen.py +++ b/quantcore/optimization/rowgen.py @@ -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): diff --git a/quantcore/optimization/service.py b/quantcore/optimization/service.py index 35395a1..1ac096b 100644 --- a/quantcore/optimization/service.py +++ b/quantcore/optimization/service.py @@ -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") diff --git a/quantcore/optimization/solvers.py b/quantcore/optimization/solvers.py index 9897814..cb9a009 100644 --- a/quantcore/optimization/solvers.py +++ b/quantcore/optimization/solvers.py @@ -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": { diff --git a/scripts/experiments/large_scale.py b/scripts/experiments/large_scale.py new file mode 100644 index 0000000..7a363fa --- /dev/null +++ b/scripts/experiments/large_scale.py @@ -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() diff --git a/tests/test_large_scale.py b/tests/test_large_scale.py new file mode 100644 index 0000000..3e93eea --- /dev/null +++ b/tests/test_large_scale.py @@ -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)