Coverage for src/cvxcla/_leverage.py: 100%
127 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-09-29 05:24 +0000
1"""Leverage (gross-exposure) constraints ``||w||_1 <= c`` via a signed lift.
3The 1-norm is not linear in ``w``, but it is polyhedral: splitting each weight into
4a long and a short leg, ``w_i = u_i - v_i`` with ``u_i, v_i >= 0``, turns
5``||w||_1 <= c`` into the single linear row ``sum(u) + sum(v) <= c``. The lifted
6problem is a CLA over the legs with one extra inequality row, and every other
7constraint carries over by substituting ``w = P x``, where ``P`` maps each leg to
8its asset with sign ``+1`` (long) or ``-1`` (short).
10Only an asset whose box straddles zero (``lower < 0 < upper``) needs two legs. A
11long-only asset (``lower >= 0``) keeps one ``+1`` leg and a short-only asset
12(``upper <= 0``) one ``-1`` leg, so a long-only problem is not enlarged at all.
14The lifted covariance ``P.T Sigma P`` is singular (the direction that raises both
15legs of one asset is flat), but it only enters the solve through its free block,
16and that block is a signed principal submatrix of ``Sigma`` as long as no asset has
17both legs free. :class:`SignedLift` exposes it through the ``QuadraticForm``
18interface without forming ``P.T Sigma P``, so structured backends keep their
19advantage. Both legs are never free together on the traced path: while one leg is
20off its lower bound, the other leg's multiplier equals twice the leverage row's
21multiplier, which is non-negative, so the other leg can only become free at the
22same ``lambda`` as the row releases. :func:`mask_leg_events` drops that
23competing leg event so the row release wins the tie.
25Under ``Sigma = X^T X`` and ``mu = X^T y`` the capped program is the constrained
26LASSO read as a portfolio: its budget-indexed path is the LASSO path, and the tilt
27sweep traced here is that path rescaled (Schmelzer and Hastie, arXiv:2609.25704,
28Theorem 1 and Corollary 2). See :mod:`cvxcla.lasso`.
29"""
31from __future__ import annotations
33from dataclasses import dataclass
35import numpy as np
36from numpy.typing import NDArray
37from scipy.optimize import linprog # type: ignore[import-untyped]
39from .operators import QuadraticForm
40from .types import TurningPoint
43def _scale_rows(sign: NDArray[np.float64], v: NDArray[np.float64]) -> NDArray[np.float64]:
44 """Multiply the rows of a vector or a column-stacked matrix by ``sign``."""
45 result: NDArray[np.float64] = sign[:, None] * v if v.ndim == 2 else sign * v
46 return result
49class SignedLift(QuadraticForm):
50 """The lifted quadratic form ``P.T @ base @ P`` of a signed leg-to-asset map.
52 Leg ``k`` belongs to asset ``asset[k]`` with sign ``sign[k]``, so
53 ``(P x)_i = sum_{k: asset[k] = i} sign[k] x_k``. Products are routed through
54 the ``base`` operator on the distinct assets involved; a free-block solve is a
55 signed solve on the base's principal block and is only defined while no asset
56 has two legs in the free set.
58 Attributes:
59 base: The covariance over the assets.
60 asset: The asset each leg belongs to (length = number of legs).
61 sign: ``+1`` for a long leg, ``-1`` for a short leg.
62 """
64 def __init__(self, base: QuadraticForm, asset: NDArray[np.intp], sign: NDArray[np.float64]) -> None:
65 """Wrap ``base`` with the leg map ``(asset, sign)``."""
66 self.base = base
67 self.asset = asset
68 self.sign = sign
70 @property
71 def n(self) -> int:
72 """Number of legs."""
73 return int(self.asset.shape[0])
75 def matvec(self, x: NDArray[np.float64]) -> NDArray[np.float64]:
76 """Return ``P.T @ base @ P @ x``."""
77 legs = np.arange(self.n)
78 return self.block_matvec(legs, legs, x)
80 def block_matvec(self, rows: object, cols: object, v: NDArray[np.float64]) -> NDArray[np.float64]:
81 """Return the ``(rows, cols)`` block of the lifted form applied to ``v``.
83 The column legs are first summed onto their distinct assets (both legs of
84 one asset collapse to a single signed entry), the base block product is
85 taken over distinct assets, and the result is spread back onto the row legs.
86 """
87 rows = np.asarray(rows, dtype=np.intp)
88 cols = np.asarray(cols, dtype=np.intp)
89 col_assets, col_inverse = np.unique(self.asset[cols], return_inverse=True)
90 collapsed = np.zeros((col_assets.shape[0], *v.shape[1:]))
91 np.add.at(collapsed, col_inverse, _scale_rows(self.sign[cols], np.asarray(v, dtype=np.float64)))
92 row_assets, row_inverse = np.unique(self.asset[rows], return_inverse=True)
93 product = np.asarray(self.base.block_matvec(row_assets, col_assets, collapsed))
94 return _scale_rows(self.sign[rows], product[row_inverse])
96 def solve_free(self, free: object, rhs: NDArray[np.float64]) -> NDArray[np.float64]:
97 """Solve on the free block, a signed principal block of ``base``.
99 Raises:
100 numpy.linalg.LinAlgError: If an asset has both legs in ``free``, which
101 makes the block singular.
102 """
103 free = np.asarray(free, dtype=np.intp)
104 assets = self.asset[free]
105 if np.unique(assets).size != assets.size:
106 msg = "both legs of an asset are free, so the lifted free block is singular"
107 raise np.linalg.LinAlgError(msg)
108 sign = self.sign[free]
109 return _scale_rows(sign, np.asarray(self.base.solve_free(assets, _scale_rows(sign, rhs))))
111 def rcond_free(self, free: object) -> float:
112 """Reciprocal condition number of the free block (``0`` if an asset has both legs free).
114 Flipping signs leaves the spectrum unchanged, so this is the base's
115 conditioning of the distinct free assets.
116 """
117 assets = self.asset[np.asarray(free, dtype=np.intp)]
118 if np.unique(assets).size != assets.size:
119 return 0.0
120 return float(self.base.rcond_free(assets))
123@dataclass(frozen=True)
124class LeverageLift:
125 """The signed leg structure of a leverage-constrained problem.
127 Attributes:
128 asset: The asset each leg belongs to.
129 sign: ``+1`` for a long leg, ``-1`` for a short leg.
130 partner: For a leg of a two-legged asset, the index of its other leg;
131 ``-1`` for the single leg of a long-only or short-only asset.
132 lower: Lower bounds of the legs.
133 upper: Upper bounds of the legs.
134 """
136 asset: NDArray[np.intp]
137 sign: NDArray[np.float64]
138 partner: NDArray[np.intp]
139 lower: NDArray[np.float64]
140 upper: NDArray[np.float64]
142 @classmethod
143 def from_bounds(cls, lower: NDArray[np.float64], upper: NDArray[np.float64]) -> LeverageLift:
144 """Build the legs from the asset box ``lower <= w <= upper``.
146 An asset with ``lower >= 0`` gets one long leg on ``[lower, upper]``; one
147 with ``upper <= 0`` gets one short leg on ``[-upper, -lower]``; one with
148 ``lower < 0 < upper`` gets a long leg on ``[0, upper]`` and a short leg on
149 ``[0, -lower]``, stored next to each other.
150 """
151 asset: list[int] = []
152 sign: list[float] = []
153 leg_lower: list[float] = []
154 leg_upper: list[float] = []
155 partner: list[int] = []
156 for i, (lo, up) in enumerate(zip(lower.tolist(), upper.tolist(), strict=True)):
157 if lo >= 0.0:
158 asset.append(i)
159 sign.append(1.0)
160 leg_lower.append(lo)
161 leg_upper.append(up)
162 partner.append(-1)
163 elif up <= 0.0:
164 asset.append(i)
165 sign.append(-1.0)
166 leg_lower.append(-up)
167 leg_upper.append(-lo)
168 partner.append(-1)
169 else:
170 k = len(asset)
171 asset += [i, i]
172 sign += [1.0, -1.0]
173 leg_lower += [0.0, 0.0]
174 leg_upper += [up, -lo]
175 partner += [k + 1, k]
176 return cls(
177 asset=np.array(asset, dtype=np.intp),
178 sign=np.array(sign),
179 partner=np.array(partner, dtype=np.intp),
180 lower=np.array(leg_lower),
181 upper=np.array(leg_upper),
182 )
184 def columns(self, matrix: NDArray[np.float64]) -> NDArray[np.float64]:
185 """Return ``matrix @ P``: each asset column copied onto its legs with the leg's sign."""
186 result: NDArray[np.float64] = matrix[:, self.asset] * self.sign
187 return result
189 def to_assets(self, x: NDArray[np.float64], n: int) -> NDArray[np.float64]:
190 """Return the asset weights ``P @ x`` of the leg weights ``x``."""
191 weights = np.zeros(n)
192 np.add.at(weights, self.asset, self.sign * x)
193 return weights
195 def any_leg(self, mask: NDArray[np.bool_], n: int) -> NDArray[np.bool_]:
196 """Return, per asset, whether any of its legs is set in ``mask``."""
197 return np.bincount(self.asset, weights=mask.astype(np.float64), minlength=n) > 0
199 def with_cap(
200 self, g: NDArray[np.float64], h: NDArray[np.float64], cap: float | None
201 ) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
202 """Return the lifted ``G P``, ``h`` with the gross-exposure row ``sum(x) <= cap`` appended.
204 ``cap=None`` maps the rows without appending the cap (a cap already
205 resolved into tightened bounds).
206 """
207 lifted = self.columns(g)
208 if cap is None:
209 return lifted, h
210 return np.vstack([lifted, np.ones((1, self.asset.shape[0]))]), np.append(h, cap)
212 def to_turning_point(self, tp: TurningPoint, n: int, p: int) -> TurningPoint:
213 """Map a lifted turning point back to the ``n`` asset weights and the first ``p`` rows."""
214 return TurningPoint(
215 lamb=tp.lamb,
216 weights=self.to_assets(tp.weights, n),
217 free=self.any_leg(tp.free, n),
218 active_ineq=tp.active_ineq[:p],
219 )
221 def net(self, x: NDArray[np.float64]) -> NDArray[np.float64]:
222 """Net out overlapping legs: subtract ``min(u_i, v_i)`` from both legs of each asset.
224 The asset weights are unchanged and the gross exposure can only fall, so
225 every constraint that held for ``x`` still holds.
226 """
227 paired = self.partner >= 0
228 overlap = np.zeros_like(x)
229 overlap[paired] = np.minimum(x[paired], x[self.partner[paired]])
230 return x - overlap
233def mask_leg_events(
234 box: NDArray[np.float64], partner: NDArray[np.intp], at_lower: NDArray[np.bool_]
235) -> NDArray[np.float64]:
236 """Drop the leave-a-bound events of a leg whose partner is off its lower bound.
238 While one leg of an asset is free (or at its upper bound), the other leg's
239 multiplier is twice the leverage row's multiplier, so it reaches zero only
240 exactly when the row releases. Freeing that leg would put both legs in the free
241 set and make the lifted free block singular; the row release is the event that
242 must fire, so the leg's leave events (columns 2 and 3) are removed.
244 Args:
245 box: The ``(legs, 4)`` box-event matrix.
246 partner: The partner leg of every leg (``-1`` for none).
247 at_lower: Mask of the legs held at their lower bound.
249 Returns:
250 A copy of ``box`` with the suppressed events set to ``-inf``.
251 """
252 paired = partner >= 0
253 suppress = np.zeros(partner.shape[0], dtype=bool)
254 suppress[paired] = ~at_lower[partner[paired]]
255 box = box.copy()
256 box[suppress, 2:] = -np.inf
257 return box
260def _gross_lp(
261 lift: LeverageLift,
262 a: NDArray[np.float64],
263 b: NDArray[np.float64],
264 g: NDArray[np.float64],
265 h: NDArray[np.float64],
266 sense: float,
267) -> tuple[float, NDArray[np.float64]] | None:
268 """Optimise the gross exposure ``sum(x)`` over the lifted feasible set.
270 Minimises ``sense * sum(x)`` subject to ``A P x = b``, ``G P x <= h`` and the leg
271 box, via HiGHS.
273 Returns:
274 ``(sum(x), reduced costs of the leg lower bounds)`` at the optimum, or
275 ``None`` if the linear program has no solution (the caller's own first
276 vertex then reports the infeasibility).
277 """
278 n_legs = lift.asset.shape[0]
279 has_ineq = g.shape[0] > 0
280 result = linprog(
281 c=np.full(n_legs, sense),
282 A_eq=lift.columns(a),
283 b_eq=b,
284 A_ub=lift.columns(g) if has_ineq else None,
285 b_ub=h if has_ineq else None,
286 bounds=list(zip(lift.lower, lift.upper, strict=True)),
287 method="highs",
288 )
289 if not result.success:
290 return None
291 return float(np.sum(result.x)), np.asarray(result.lower.marginals, dtype=np.float64)
294def tighten_at_minimum_gross(
295 lower: NDArray[np.float64],
296 upper: NDArray[np.float64],
297 a: NDArray[np.float64],
298 b: NDArray[np.float64],
299 g: NDArray[np.float64],
300 h: NDArray[np.float64],
301 leverage: float,
302 tol: float,
303) -> tuple[NDArray[np.float64], NDArray[np.float64], bool]:
304 """Resolve a cap sitting at the smallest feasible gross exposure.
306 When ``leverage`` equals ``c_min = min ||w||_1`` over the feasible set, the cap
307 row is tight at every feasible point and implies a set of zero legs through the
308 other rows. With a fully-invested budget, ``leverage = 1`` forces every short leg
309 to zero. Carrying the cap row alongside those rows makes the maximum-return vertex
310 degenerate: the free set cannot span them all. This is also the one cap at which
311 Slater's condition fails.
313 The feasible set is then the optimal face of ``min sum(x)``. By complementary
314 slackness, a leg whose lower bound carries a positive reduced cost is zero on
315 that whole face. Pinning it there is therefore exact: the short leg of a split
316 asset pins ``lower = 0``, and its long leg pins ``upper = 0``. If the cap is
317 then implied by the tightened box, because ``max sum(x)`` over it is at most
318 ``leverage``, the row is dropped.
320 Args:
321 lower: Asset lower bounds.
322 upper: Asset upper bounds.
323 a: Equality-constraint matrix over the assets.
324 b: Equality-constraint right-hand side.
325 g: Inequality-constraint matrix over the assets.
326 h: Inequality-constraint right-hand side.
327 leverage: The gross-exposure cap ``c``.
328 tol: Tolerance for comparing the cap with ``c_min`` and for a positive
329 reduced cost.
331 Returns:
332 ``(lower, upper, keep_cap)``: the (possibly tightened) asset bounds and
333 whether the cap row is still needed. Unchanged bounds and ``True`` when the
334 cap exceeds ``c_min``.
335 """
336 lift = LeverageLift.from_bounds(lower, upper)
337 minimum = _gross_lp(lift, a, b, g, h, sense=1.0)
338 if minimum is None or leverage > minimum[0] + tol * max(1.0, minimum[0]):
339 return lower, upper, True
341 pinned = minimum[1] > tol
342 lower, upper = lower.copy(), upper.copy()
343 split = lift.partner >= 0
344 short = split & (lift.sign < 0) & pinned
345 long = split & (lift.sign > 0) & pinned
346 lower[lift.asset[short]] = 0.0
347 upper[lift.asset[long]] = 0.0
349 maximum = _gross_lp(LeverageLift.from_bounds(lower, upper), a, b, g, h, sense=-1.0)
350 keep_cap = maximum is None or maximum[0] > leverage + tol * max(1.0, leverage)
351 return lower, upper, keep_cap